MathGroup Archive 2001

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Functional programming

  • To: mathgroup at
  • Subject: [mg27155] Re: [mg27131] Functional programming
  • From: Daniel Lichtblau <danl at>
  • Date: Thu, 8 Feb 2001 04:40:43 -0500 (EST)
  • References: <>
  • Sender: owner-wri-mathgroup at

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

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]],
    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

  • Prev by Date: Re: LogPlot/Plot-Identity
  • Next by Date: Re: Functional programming
  • Previous by thread: Re: Functional programming
  • Next by thread: Re: Functional programming