# Matrices and polynomials: MMA 8 vs. 9?

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?

• I think MatrixPower chokes on your matrix not being positive-definite... Jul 19, 2013 at 13:12
• It's funny, because the documentation says that MatrixPower was last modified in version 6. Jul 19, 2013 at 13:13
• In version 9, have you tried MatrixFunction ? Jul 19, 2013 at 13:25
• @b.gatessucks I did try MatrixFunction but I seem to remember it had problems with using list entries. I'll get back to you later when I can look at those files. Jul 19, 2013 at 14:52

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.

• I agree, this seems like an adequate solution. I don't have enough rep to vote up bill s but I appreciate his help. Jul 20, 2013 at 23:14
• Note that the r[x_] := R[[i]] doesn't make much sense... Jul 21, 2013 at 0:16
• @nardol5 I just want to point out that MatrixFunction shows the same behaviour as MatrixPower if asked to compute T^0 (where T is not positive definite), which caused your code to fail in the first place! Try MatrixFunction[#^0 &,T]. Jul 21, 2013 at 0:26
• @sebhofer Quite right. I didn't look at the code closely enough. One can simply have applyPoly[R[[i]], x, T] in the Table. Jul 21, 2013 at 0:44

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.

• Unfortunately, when I replace applyPoly[r[x], x, T] with MatrixFunction[r[x],T] in the code above, it returns: Jul 19, 2013 at 18:44
• MatrixFunction::nunipf: -1-2 x-3 x^2+x^3 is not a univariate pure function. I suspect that MatrixFunction has issues with list elements, but I don't see how I can avoid using a list in my program. @bill s Jul 19, 2013 at 18:52
• @nardol5 The error message says it all... Try MatrixFunction[Function[x, R[[1]]], T] Jul 20, 2013 at 0:21
• @sebhofer Your suggestion works for individual cases, but in my program it produces a long list of expressions similar to: Root[-1-2#1-3#1^2+#1^3 &, 1]^2. I'm not really sure what's going on (being new to programming in MMA) but I might be able to make Function work for me somehow. Jul 20, 2013 at 3:12
• That's the symbolic form. Put N[ ] around it to see a numerical version. Jul 20, 2013 at 3:22