# 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:

http://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... –  sebhofer Jul 19 '13 at 13:12
It's funny, because the documentation says that MatrixPower was last modified in version 6. –  sebhofer Jul 19 '13 at 13:13
In version 9, have you tried MatrixFunction ? –  b.gatessucks Jul 19 '13 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. –  nardol5 Jul 19 '13 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 point 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. –  nardol5 Jul 20 '13 at 23:14
Note that the r[x_] := R[[i]] doesn't make much sense... –  sebhofer Jul 21 '13 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]. –  sebhofer Jul 21 '13 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. –  Michael E2 Jul 21 '13 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: –  nardol5 Jul 19 '13 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 –  nardol5 Jul 19 '13 at 18:52
@nardol5 The error message says it all... Try MatrixFunction[Function[x, R[[1]]], T] –  sebhofer Jul 20 '13 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. –  nardol5 Jul 20 '13 at 3:12
That's the symbolic form. Put N[ ] around it to see a numerical version. –  bill s Jul 20 '13 at 3:22