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