Re: Functional programming
- To: mathgroup at smc.vnet.net
- Subject: [mg27155] Re: [mg27131] Functional programming
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 8 Feb 2001 04:40:43 -0500 (EST)
- References: <200102070712.CAA29571@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Arturo Rania wrote:
>
> I'm kind of new to Mathematica and I'm still too procedural in my
> programming. I would like to get some advice in order to make my
> programming more efficient, and especially more Mathematica oriented. In
> order to do this I invite anyone who want to critize or improve a
> program I wrote that finds symbolic matrix powers when MatrixPower fails
> (when the matrix is non-diagonalizable). I think that it is a perfect
> example of unprofesional programming style, but it's improvement can be
> a great source of knowledge for me. Thanks for your time and good luck.
>
> NDMatrixPower[M_/;And[MatrixQ[M], Length[M]] == Length[M[[1]]], t_
> /; Or[NumericQ[t], Head[t] == Symbol]] :=
> Module[{v, w, r, EigenTable, Multiplicidad, A, B, Sol, k},
> Clear[v, w, r, EigenTable, Multiplicidad, Sol, k];
> v = Eigenvalues[M];
> w = Union[v, v];
> w[[1]];
> EigenTable = Table[{w[[i]], Count[v, w[[i]]]}, {i, 1,
> Length[w]}];
> Multiplicidad = Join@@Table[i, {j, 1, Length[EigenTable]},
> {i, 1, EigenTable[[j, 2]]}] - 1;
> A = Table[D[Power[r, j - 1], {r, Multiplicidad[[i]]}] /. r
> -> v[[i]], {i, 1, Length[M]}, {j, 1, Length[M]}];
> B = Table[D[Power[r, k], {r, Multiplicidad[[i]]}] /. r ->
> v[[i]], {i, 1, Length[M]}, {j, 1, 1}];
> Sol = Transpose[LinearSolve[A, B]][[1]];
> Plus@@Table[MatrixPower[M, i - 1]*Sol[[i]], {i, 1,
> Length[v]}] /. k -> t];
>
> Any comments are welcome,
>
> Arturo Rania
The code is basically fine. There are some small improvements you can
make though they are mostly for readability.
The Clear[] is not necessary, and the w[[1]] line is also not necessary.
Union needs only one argument for your purposes.
I think you should allow non-symbol powers e.g. t[4], hence only avoid
non-integer exact NumberQ powers.
The predicates might be separated from the main module.
Likewise the run-length encoding of eigenvalues with multiplicities.
You can use Range to get lists from 1 to multiplicity of eigenvalue, and
then flatten these lists into one master list. This is slightly cleaner
than what you have and, I think, a bit simpler to follow.
You'll have a bug if the eigenvalues are not initially sorted so that
repeats are non-adjacent: A can get messed up. So I'd explicitly sort
them.
The B computation does not need a double indexed Table, and the
LinearSolve can be made simpler. I am not certain I changed this
correctly although it does do the right thing in the simple example I
show below.
Plus @@ Table can be relpaced by Sum.
Putting this all together, here is what I get.
squareMatrixQ[m_] := MatrixQ[m] && Length[m]==Length[m[[1]]]
badPower[t_] := !NumberQ[t] || (!IntegerQ[t] && Precision[t]===Infinity)
runLength[vals_] := Map[{#[[1]],Length[#]}&, Split[Sort[vals]]]
NDMatrixPower[M_?squareMatrixQ, t_?badPower] :=
Module[{v, r, EigenTable, Multiplicidad, A, B, Sol, k},
v = Sort[Eigenvalues[M]];
EigenTable = runLength[v];
Multiplicidad = Flatten[Map[Range[#[[2]]]&, EigenTable]] - 1;
A = Table[D[r^(j-1),{r,Multiplicidad[[i]]}] /. r->v[[i]],
{i,Length[M]}, {j,Length[M]}];
B = Table[D[r^k,{r,Multiplicidad[[i]]}] /. r->v[[i]],
{i,Length[M]}];
Sol = LinearSolve[A, B];
Sum[MatrixPower[M,i-1]*Sol[[i]], {i,Length[v]}] /. k->t
]
In[40]:= NDMatrixPower[{{2, 2}, {0, 2}}, t]
t t t t t
Out[40]= {{-(2 (-1 + t)) + 2 t, 2 t}, {0, -(2 (-1 + t)) + 2 t}}
MatrixPower for such examples works in our development version. So we
can readily check this simple example.
In[41]:= MatrixPower[{{2, 2}, {0, 2}}, t]
t t t
Out[41]= {{2 , 2 t}, {0, 2 }}
In[42]:= Together[%-%%]
Out[42]= {{0, 0}, {0, 0}}
While I get a simpler form, I would not be terribly surprised if your
code is a bit faster on larger examples.
Daniel Lichtblau
Wolfram Research
- References:
- Functional programming
- From: "Arturo Rania" <abasida@hotmail.com>
- Functional programming