# Bad result from evaluating a determinant

I have a problem in calculating the determinant of the matrix in the picture.

If one applies Laplace formula for computing determinants, multiplying all the elements of the first column with their minors, a result linear in the variables Ns and Nst must be there. However, when I calculate the determinant with Mathematica, it gives also terms proportional to Ns^2 (see lower part of the pictures) which should not be there at all. With smaller matrices the problem is not there, but with the one in the figure (which is 12 x 12) the problem arises and it stays also for larger matrices.

I have no clue why this happens, if anyone can help to clarify and/or suggest some other method to compute determinants it would be great!

attached there is also the code to generate the matrix:

n = 12;
(*Analitical definitions*)
ω = Table[ωω[i], {i, 1, n}];
k = Table[kk[i], {i, 1, n - 1}];
d = Table[dd[i], {i, 1, n}];
e = Table[ee[i], {i, 1, n}];

(*Numerical definitions*)
Do[ωω[i] = 0.1, {i, 1, n}];
Do[kk[i] = 1., {i, 1, n}];
γ = 1.;
Do[ee[i] = 1., {i, 1, n}];
Do[dd[i] = 0., {i, 1, n}];
(*Definition of matrix M and vector h*)
f = Table["f" <> ToString[i] // ToExpression, {i, 1, n}];
h = Table["f" <> ToString[i] // ToExpression, {i, 1, n}];
x = Table["x" <> ToString[i] // ToExpression, {i, 1, n}];
f[[1]] = s d[[1]] - I (k[[1]] e[[2]] + ω[[1]] e[[1]]);
Do[
f[[i]] =
s d[[i]] -
I (k[[i]] e[[i + 1]] + ω[[i]] e[[i]] +
k[[i - 1]] e[[i - 1]]);, {i, 2, n - 2}];
f[[n - 1]] =
s d[[n - 1]] -
I (k[[n - 1]] e[[n]] + ω[[n - 1]] e[[n - 1]] +
k[[n - 2]] e[[n - 2]]);
f[[n]] = s d[[n]] - I (ω[[n]] e[[n]] + k[[n - 1]] e[[n - 1]]);
Do[h[[i]] = f[[i]], {i, 1, n - 2}];
h[[n - 1]] = f[[n - 1]] - I k[[n - 1]] (Ns + Nst);
h[[n]] = f[[n]] + (s - I ω[[n]]) Ns - (s + I ω[[n]]) Nst;
M = Table[
"M" <> ToString[i] <> ToString[j] // ToExpression, {i, 1, n}, {j,
1, n}];
M[[1, 1]] = s^2 + ω[[1]]^2 + k[[1]]^2;
M[[1, 2]] = k[[1]] (ω[[1]] + ω[[2]]);
M[[1, 3]] = k[[1]] k[[2]];
Do[M[[1, j]] = 0, {j, 4, n}];
M[[2, 1]] = k[[1]] (ω[[1]] + ω[[2]]);
M[[2, 2]] = s^2 + ω[[2]]^2 + k[[2]]^2 + k[[1]]^2;
M[[2, 3]] = k[[2]] (ω[[2]] + ω[[3]]);
If[n > 3, M[[2, 4]] = k[[2]] k[[3]]];
Do[M[[2, j]] = 0, {j, 5, n}];
If[n > 3, M[[n - 1, n - 3]] = k[[n - 2]] k[[n - 3]]];
M[[n - 1, n - 2]] = k[[n - 2]] (ω[[n - 2]] + ω[[n - 1]]);
M[[n - 1, n - 1]] =
s^2 + ω[[n - 1]]^2 + k[[n - 1]]^2 + k[[n - 2]]^2;
M[[n - 1, n]] = k[[n - 1]] (ω[[n - 1]] + ω[[n]]);
Do[M[[n - 1, j]] = 0, {j, 1, n - 4}];
M[[n, n - 2]] = k[[n - 1]] k[[n - 2]];
M[[n, n - 1]] = k[[n - 1]] (ω[[n - 1]] + ω[[n]]);
M[[n, n]] = s (s + γ) + ω[[n]]^2 + k[[n - 1]]^2;
Do[M[[n, j]] = 0, {j, 1, n - 3}];
If[n > 3,
Do[M[[i, j]] =
KroneckerDelta[i, j + 2] k[[i - 1]] k[[i - 2]] +
KroneckerDelta[i,
j + 1] k[[i - 1]] (ω[[i - 1]] + ω[[i]]) +
KroneckerDelta[i,
j] (s^2 + ω[[i]]^2 + k[[i]]^2 + k[[i - 1]]^2) +
KroneckerDelta[i,
j - 1] k[[i]] (ω[[i]] + ω[[i + 1]]) +
KroneckerDelta[i, j - 2] k[[i]] k[[i + 1]] +
KroneckerDelta[i, n] KroneckerDelta[j,
n] ((s (s + γ) + ω[[n]]^2 +
k[[n - 1]]^2) - (s^2 + ω[[i]]^2 + k[[i]]^2 +
k[[i - 1]]^2));
, {i, 3, n - 2}, {j, 1, n}]]
MatrixForm[h];
MatrixForm[M];
(*Definition of Subscript[M, 1]*)
Subscript[M, 1] =
Table["M" <> ToString[i] <> ToString[j] // ToExpression, {i, 1,
n}, {j, 1, n}];
Subscript[M, 1] = M;
Do[Subscript[M, 1] =
ReplacePart[Subscript[M, 1], {i, 1} -> h[[i]]], {i, 1, n}]
MatrixForm[Subscript[M, 1]]
(*Determinant of Subscript[M, 1] without rationalize*)
Collect[Det[Subscript[M,
1]], Ns](*Determinant of Subscript[M, 1] with rationalize*)
Collect[Det[Rationalize[Subscript[M, 1]]], Ns]


UPDATE:

I have tried to use the command "Rationalize" before computing the determinants and that solved the issue! Thanks to everyone for the help.

• Can you provide the matrix in a copy/pasted form, or something? Nobody is going to type a 12x12 matrix in by hand. – evanb Jan 11 '18 at 14:59
• Difficult to guess what is going on. I can't even see that example, let alone replicate it for diagnosing. – Daniel Lichtblau Jan 11 '18 at 15:52
• I believe the problem is being caused by the inexactness of machine arithmetic. Try using Mathematica's arbitrary precision arithmetic instead. I predict you will get a better result. – m_goldberg Jan 11 '18 at 18:24
• What about Det[Rationalize[M2]]? – anderstood Jan 11 '18 at 18:29
• Chop[Det[Subscript[M, 1]]] will remove terms that amount to numerical fuzz. – Daniel Lichtblau Jan 12 '18 at 16:37

I think the reason is that a clever implementation of Det would not use the Laplace formula for its computational complexity is immense. For purely numerical matrices, the determinant can be obtained quite accurately from certain factorizations (e.g. LU-decomposition). Maybe that's why a rounding error might be introduced. Note that the coefficient of Ns^2 is quite small, so this seems to be indeed a rounding error (as such, it is rather large, the reason might be the mixture of symbolic and numeric entries which prevents suitable row permutations.)
You can try to apply Rationalize to your matrix and apply the determinant afterwards. This should give you a more precise answer.