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 - cx^3 - bx^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?

Community
  • 1
nardol5
  • 133
  • 2

2 Answers2

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.

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
bill s
  • 68,936
  • 4
  • 101
  • 191
  • 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 00: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 03:12
  • That's the symbolic form. Put N[ ] around it to see a numerical version. – bill s Jul 20 '13 at 03:22
  • @bills Applying N[] does not change the expression. – nardol5 Jul 20 '13 at 03:28
  • @nardol5 I'm not quite sure what you are doing then because it gives me 13.1578 – sebhofer Jul 20 '13 at 09:18
  • @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 Jul 20 '13 at 15:47
  • 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 Jul 20 '13 at 16:15
  • @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 Jul 20 '13 at 18:42
  • 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 Jul 20 '13 at 18:45
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.

Community
  • 1
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • 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
  • 1
    Note that the r[x_] := R[[i]] doesn't make much sense... – sebhofer Jul 21 '13 at 00:16
  • 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 Jul 21 '13 at 00: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 00:44