Matrices and polynomials: MMA 8 vs. 9?

Question

I want a function to take a polynomial from a list, plug a certain matrix T into that polynomial, and return the answer as another matrix. My program uses the applyPoly function suggested by a user in response to:

https://stackoverflow.com/questions/17394035/in-mathematica-how-calculate-pa-where-p-is-a-polynomial-and-a-is-a-square-mat

Here is my code:

applyPoly[poly_, var_, A_?MatrixQ] := 
  With[{c = CoefficientList[poly, var]}, 
    c.MapIndexed[MatrixPower[A, #2[[1]] - 1] &, c]]

Compute[a_, b_, c_] := (
  p[x_] := x^4 - c*x^3 - b*x^2 - a*x; 
  T := {{0, 0, 0, 0}, {1, 0, 0, a}, {0, 1, 0, b}, {0, 0, 1, c}}; 
  L := FactorList[p[x]]; s = Length[L]; R = {}; 
  Do[R = Append[R, Simplify[p[x]/L[[i + 1, 1]]]], {i, s - 1}]; 
  u = Length[R]; 
  Do[r[x_] := R[[i]]; Print[applyPoly[r[x], x, T]], {i, u}])

Compute[1,2,3]

In Mathematica 8, this returns the correct matrix output:

{{-1,0,0,0},{-2,0,0,0},{-3,0,0,0},{1,0,0,0}}

{{0,0,0,0},{1,0,0,1},{0,1,0,2},{0,0,1,3}}

But in Mathematica 9, I get an error:

Dot::rect: Nonrectangular tensor encountered. >>

{-1,-2,-3,1}.{MatrixPower[{{0,0,0,0}, {1,0,0,1},{0,1,0,2},{0,0,1,3}},0],{{0,0,0,0}, {1,0,0,1},{0,1,0,2},{0,0,1,3}}, {{0,0,0,0}, {0,0,1,3},{1,0,2,7}, {0,1,3,11}}, {{0,0,0,0}, {0,1,3,11}, {0,2,7,25}, {1,3,11,40}}}

Dot::rect: Nonrectangular tensor encountered. >>

{0,1}.{MatrixPower[{{0,0,0,0},{1,0,0,1},{0,1,0,2},{0,0,1,3}},0],{{0,0,0,0},{1,0,0,1},{0,1,0,2},{0,0,1,3}}}

I understand that 8 and 9 must handle matrix or list objects in different ways. What is the difference, and how can I alter my program so that it works in 9?

Answer 1

Don't use Print for output; change Do to Table. Then the output will be the matrices (with the Root objects), and N will work.

applyPoly[poly_, var_, A_?MatrixQ] :=  MatrixFunction[poly /. var -> # &, A]

Compute[a_, b_, c_] := (p[x_] := x^4 - c*x^3 - b*x^2 - a*x;
  T := {{0, 0, 0, 0}, {1, 0, 0, a}, {0, 1, 0, b}, {0, 0, 1, c}};
  L := FactorList[p[x]]; s = Length[L]; R = {};
  Do[R = Append[R, Simplify[p[x]/L[[i + 1, 1]]]], {i, s - 1}];
  u = Length[R];
  Table[(*r[x_] := R[[i]];*) applyPoly[(*r[x]*)R[[i]], x, T], {i, u}])

Edit: I commented out the (unnecessary) definition of r, thanks to sebhofer, and replaced it with R[[i]].

---

Example

Compute[1, 2, 3] // N

(* {{ {-1. + 0. I, 0., 0., 0.},
      {-2. + 0. I, 2.53644*10^-16 + 0. I, 3.20011*10^-16 + 0. I, 1.68815*10^-15 + 0. I},
      {-3. - 4.93038*10^-32 I, 7.28119*10^-16 + 0. I, 8.93665*10^-16 + 0. I, 3.69631*10^-15 + 0. I},
      {1. + 0. I, 3.20011*10^-16 + 0. I, 1.68815*10^-15 + 0. I, 5.95812*10^-15 + 0. I}},
    { {0., 0., 0., 0.},
      {1. + 0. I, 0. + 0. I, 0. + 1.38778*10^-17 I, 1. + 0. I},
      {-3.40006*10^-16 + 0. I, 1. + 5.55112*10^-17 I, 0. + 0. I, 2. + 2.77556*10^-17 I},
      {-4.996*10^-16 + 0. I, 2.77556*10^-17 + 0. I, 1. + 0. I, 3. + 0. I}}} *)

---

There are other ways (1, 2, 3) to turn a polynomial poly into a function, but the above was the easiest, given your setup.

---

The above was rather too long to put into a comment, but I think bill s deserves the credit for pointing out MatrixFunction.

Answer 2

MatrixFunction is the way to go... for example:

MatrixFunction[#^5 + 2 #^2 + 1 &, {{a, 1}, {0, b}}]

gives the polynomial function x^5+2x^2+1 when x is the symbolic matrix {{a, 1}, {0, b}}. It also works for more complicated functions such as trigonometric functions. And of course, it works on numerical matrices as well. It does require that the matrices be square, which makes sense because otherwise the x^n terms don't make any sense.