Matrices and polynomials: MMA 8 vs. 9?

3

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?

nardol5

Posted 2013-07-19T13:04:46.593

Reputation: 133

I think MatrixPower chokes on your matrix not being positive-definite... – sebhofer – 2013-07-19T13:12:20.720

It's funny, because the documentation says that MatrixPower was last modified in version 6. – sebhofer – 2013-07-19T13:13:55.557

In version 9, have you tried MatrixFunction ? – b.gates.you.know.what – 2013-07-19T13:25:52.183

@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 – 2013-07-19T14:52:57.850

Answers

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.

Michael E2

Posted 2013-07-19T13:04:46.593

Reputation: 190 928

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 – 2013-07-20T23:14:46.780

1Note that the r[x_] := R[[i]] doesn't make much sense... – sebhofer – 2013-07-21T00:16:07.180

1@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 – 2013-07-21T00:26:05.200

@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 – 2013-07-21T00:44:54.997

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.

bill s

Posted 2013-07-19T13:04:46.593

Reputation: 62 963

Unfortunately, when I replace applyPoly[r[x], x, T] with MatrixFunction[r[x],T] in the code above, it returns: – nardol5 – 2013-07-19T18:44:17.700

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 – 2013-07-19T18:52:52.437

@nardol5 The error message says it all... Try MatrixFunction[Function[x, R[[1]]], T] – sebhofer – 2013-07-20T00:21:22.567

@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 – 2013-07-20T03:12:21.987

That's the symbolic form. Put N[ ] around it to see a numerical version. – bill s – 2013-07-20T03:22:22.577

@bills Applying N[] does not change the expression. – nardol5 – 2013-07-20T03:28:21.760

@nardol5 I'm not quite sure what you are doing then because it gives me 13.1578 – sebhofer – 2013-07-20T09:18:09.027

@bills The output I'm getting includes several hundred Root expressions, with several polynomials in x thrown in for good measure. I tried placing N[ ] in my program inside the Print command, in several different ways. It did not lead to any understandable output either. – nardol5 – 2013-07-20T15:47:39.343

Why not build up your answer from the pieces you have, rather than trying to solve a single large problem all at once? – bill s – 2013-07-20T16:15:55.443

@bills I guess that's what I need to do. This is actually part of a larger program I was writing for a research project. I was just trying to make my code as compact as possible. – nardol5 – 2013-07-20T18:42:46.100

Also, I think my institution is upgrading their MMA license to 9 in a few days, so I won't be able to take advantage of the difference in 8 anymore. – nardol5 – 2013-07-20T18:45:17.643