Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a 15x15 Matrix with all polynomial entries. I want to calculate the determinant of the matrix. To my understanding the determinant should be a (albeit high order) polynom, too. And the paper, I am following, expects this, too.

But Mathematica does not seem to agree with me here.

Anyone had this?

If I take only the upper left submatrix with rows/columns 1 to 11, the determinant is a polynom. But after that using upper left submatrix with rows/columns 1 to 12 and higher results in non-polynomial determinants.

Code Explanation: Until the line where the array T appears first, the code is just constructing polynomials of three variables (x,y,w) according to a paper ("Closed-Form Dynamic Equations of the General Stewart Platform through the Newton–Euler Approach"). The two For-loops afterwards extract the x polynomials such that the inner product of the i'th row of M15x15 with T corresponds to the i'th polynomial (pi). The Matrix M15x15 then only consists of univariate polynomials.

alpha = N[8/180*Pi];
beta = N[12.8/180*Pi];
ra = N[15.0];
rb = N[9.45];
xi  = ra*Cos[{0, 2*alpha, 2/3*Pi, 2/3*Pi + 2*alpha, 4/3*Pi, 4/3*Pi + 2*alpha}]
yi = ra*Sin[{0, 2*alpha, 2/3*Pi, 2/3*Pi + 2*alpha, 4/3*Pi, 4/3*Pi + 2*alpha}]
zi = {0, 0, 0, 0, 0, 0}
pi = rb*Cos[1/3*Pi + {0, 2*beta, 2/3*Pi, 2/3*Pi + 2*beta, 4/3*Pi, 4/3*Pi + 2*beta}]
qi = rb*Sin[1/3*Pi + {0, 2*beta, 2/3*Pi, 2/3*Pi + 2*beta, 4/3*Pi, 4/3*Pi + 2*beta}]
ri = {0, 0, 0, 0, 0, 0};
Li = {20.0, 20.0, 20.0, 20.0, 20.0, 20.0};
mi = 1/2*(Li^2 - xi^2 - yi^2 - pi^2 - qi^2);
M6x10 = Array[0, {6, 10}];
For[i = 1, i <= 6, i++,
 M6x10[[i, ;;]] = {pi[[i]]*xi[[i]], pi[[i]]*yi[[i]], -pi[[i]], 
   qi[[i]]*xi[[i]], qi[[i]]*yi[[i]], -qi[[i]], xi[[i]], yi[[i]], -0.5,
    mi[[i]]}
]
a0 = Det[ M6x10[[;; , 1 ;; 6]] ];
a6x4 = Array[0, {6, 4}];
M6x6 =  M6x10[[;; , 1 ;; 6]];
For[i = 1, i <= 6, i++,
 For[j = 1, j <= 4, j++,
  tmp = M6x6;
  tmp[[;; , i]] = M6x10[[;; , j + 6]];
  a6x4[[i, j]] = Det[tmp];
 ]
]
AA = a6x4[[1, 1]]*x + a6x4[[1, 2]]*y + a6x4[[1, 3]]*w + a6x4[[1, 4]];
BB = a6x4[[2, 1]]*x + a6x4[[2, 2]]*y + a6x4[[2, 3]]*w + a6x4[[2, 4]];
CC = a6x4[[3, 1]]*x + a6x4[[3, 2]]*y +                  a6x4[[3, 4]];
DD = a6x4[[4, 1]]*x + a6x4[[4, 2]]*y + a6x4[[4, 3]]*w + a6x4[[4, 4]];
FF = a6x4[[5, 1]]*x + a6x4[[5, 2]]*y + a6x4[[5, 3]]*w + a6x4[[5, 4]];
GG = a6x4[[6, 1]]*x + a6x4[[6, 2]]*y +                  a6x4[[6, 4]];
f1 = 1 - a0^(-2)*(AA^2 + BB^2);
f2 = 1 - a0^(-2)*(DD^2 + FF^2);
f3 = 1 - a0^(-2)*(AA*DD + BB*FF);
f4 = a0^(-1)*(-CC*AA*x + BB*y);
f5 = a0^(-1)*(-GG*DD*x + FF*y);
f6 = w - x^2 - y^2;
p1 = f1*f6 - f4^2;
p2 = f2*f6 - f5^2;
p3 = f3*f6 - f4*f5;
p4 = f1*f5 - f3*f4;
p5 = f2*f4 - f3*f5;
p6 = f1*f2 - f3^2; 
p7 = -DD*p1 + AA*p3 - x*p4;
p8 = FF*p1 - BB*p3 + y*p4;
p9 = DD*p3 - AA*p2 - x*p5;
p10 = -FF*p1 + BB*p3 + y*p5;
p11 = AA*p5 + DD*p4 + x*p6;
p12 = BB*p5 + FF*p4 + y*p6;
p13 = -FF*p7 - BB*p9 + y*p11;
p14 = -CC*p5 + DD*p7 + BB*p10 - a0^2*p2 - x*p11;
p15 = p14 - DD*p7 - FF*p8 - AA*p9 - BB*p10 + x*p11 + y*p12;
T = {w^4, w^3 y, w^2 y^2, w y^3, y^4, w^3, w^2 y, w y^2, y^3, w^2, w y, y^2, w, y, 1};
pAll = {p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14, 
   p15};
M15x15 = Array[0, {15, 15}];
For[i = 1, i <= 15, i++,
 For[j = 1, j <= 14, j++,
  M15x15[[i, j]] = Coefficient[ pAll[[i]], T[[j]] ] /. {y -> 0, w -> 0};
  ]
]
For[i = 1, i <= 15, i++,
 M15x15[[i, 15]] = ( pAll[[i]] /. {y -> 0, w -> 0});
]
determinante = Det[M15x15];
share|improve this question

closed as off-topic by rasher, belisarius, bobthechemist, Yves Klett, ubpdqn Apr 11 at 2:37

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – rasher, belisarius, bobthechemist, Yves Klett, ubpdqn
If this question can be reworded to fit the rules in the help center, please edit the question.

2  
Please add (copyable) code that reproduces your problem. Try to keep it as small as possible. –  Sjoerd C. de Vries Apr 8 at 20:11
    
ignore this one –  alexander Apr 8 at 20:21
    
What do you mean? Would you like to retract your question? –  Sjoerd C. de Vries Apr 8 at 20:24
    
sry just messed up, original post is more detailed now –  alexander Apr 8 at 20:36
    
With exact arithmetic the divisions done internally would exactly cancel. Not likely to happen when coefficients are approximate. On the other hand, with exact coefficients the process might well hang anyway. You might try using PolynomialDivision to get quotient and remainder pairs. If coeffs of remainder terms are suitably small, just drop them. –  Daniel Lichtblau Apr 8 at 21:20

1 Answer 1

Well, this gives a polynomial answer:

determinante = N@Expand@Det[SetPrecision[M15x15, Infinity]]
(*
  -8.85105*10^55 - 5.50608*10^69 x + 2.79901*10^80 x^2 - 
   5.66691*10^90 x^3 + 3.30477*10^100 x^4 + 4.27517*10^110 x^5 - 
   ...
   7.20261*10^254 x^51 + 4.72779*10^252 x^52 + 6.39829*10^250 x^53 - 
   4.74757*10^247 x^54 - 4.47661*10^245 x^55
*)

The size of the coefficients makes me worry about the numerics. If we set the precision of the parameters that go into making M15x15 higher, the first six digits of the result settle down when the precision gets to 40 digits or more.

determinante = N@Expand@Det[SetPrecision[M15x15, Infinity]]
(*
  -8.88605*10^55 - 5.5058*10^69 x + 2.7989*10^80 x^2 - 
   5.66674*10^90 x^3 + 3.30466*10^100 x^4 + 4.27514*10^110 x^5 - 
   ...
   3.8117*10^255 x^51 + 4.84531*10^252 x^52 + 6.37814*10^250 x^53 - 
   5.03708*10^247 x^54 - 2.77218*10^245 x^55
*)

But that's a rather naive approach to the numerics.

share|improve this answer
    
Excellent answer. –  Daniel Lichtblau Apr 9 at 16:06

Not the answer you're looking for? Browse other questions tagged or ask your own question.