Jtdylktuy

I have a large matrix with numerical components and want to set the determinant to zero using the parameter h (see below). Naively, I would have expected that h sets the determinant to (approximately) zero, which isn't the case. On top of that, the order of applying the rule sol seems to affects the final outcome for a reason to don't see.

My output of the code below is:

{h -> -0.744736 + 4.42008 I}



0.0445865 - 0.0285418 I



0.0545654 - 0.114258 I

I am not familiar with how Mathematica handles floating point numbers so that's probably where my error lies. I have also tried to increase the precision with SetPrecision, but without success.

mat={{0.16 - (0.36 + 0.001 I) h - (1.35808 - 

  0.00120116 I) h^2 - (0.49603 - 0.00137214 I) h^3 - (0.11307 - 

  0.00105331 I) h^4 + (0.249794 - 0.000384238 I) h^5 - 

0.39204 h^6, -0.1711 h^2 ((-0.143205 + 

   0.000186623 I) - (0.36 + 0.001 I) h - (1.15528 + 

    0.00267142 I) h^2 - (0.637164 - 0.0009801 I) h^3 + 

 1. h^4), (0.0000353051 - 1.67323*10^-6 I) h^4, 

0}, {-0.1711 h^2 ((-0.143205 + 

   0.000186623 I) - (0.36 + 0.001 I) h + (19.6394 - 

    0.00267142 I) h^2 - (0.637164 - 0.0009801 I) h^3 + 1. h^4), 

0.16 - (0.36 + 0.001 I) h - (11.3534 - 

  0.00119507 I) h^2 - (0.484268 - 0.00140481 I) h^3 - (5.0714 - 

  0.00114074 I) h^4 + (0.27061 - 0.000416258 I) h^5 - 

0.42471 h^6, -0.223386 h^2 ((-0.143205 + 

   0.000186623 I) - (0.36 + 0.001 I) h + (4.95742 - 

    0.00267502 I) h^2 - (0.637164 - 0.0009801 I) h^3 + 

 1. h^4), (0.0000484431 - 2.29589*10^-6 I) h^4}, {(0.0000353051 - 

 1.67323*10^-6 I) h^4, -0.223386 h^2 ((-0.143205 + 

   0.000186623 I) - (0.36 + 0.001 I) h + (41.4016 - 

    0.00267502 I) h^2 - (0.637164 - 0.0009801 I) h^3 + 1. h^4), 

0.16 - (0.36 + 0.001 I) h - (29.348 - 

  0.00118803 I) h^2 - (0.470698 - 0.00144251 I) h^3 - (13.9095 - 

  0.00124106 I) h^4 + (0.294629 - 0.000453204 I) h^5 - 

0.462406 h^6, -0.234771 h^2 ((-0.143205 + 

   0.000186623 I) - (0.36 + 0.001 I) h + (19.0123 - 

    0.00267319 I) h^2 - (0.637164 - 0.0009801 I) h^3 + 

 1. h^4)}, {0, (0.0000484431 - 

 2.29589*10^-6 I) h^4, -0.234771 h^2 ((-0.143205 + 

   0.000186623 I) - (0.36 + 0.001 I) h + (71.32 - 

    0.00267319 I) h^2 - (0.637164 - 0.0009801 I) h^3 + 1. h^4), 

0.16 - (0.36 + 0.001 I) h - (55.3462 - 

  0.00118568 I) h^2 - (0.466163 - 0.0014551 I) h^3 - (26.6556 - 

  0.00127449 I) h^4 + (0.302655 - 0.000465551 I) h^5 - 

0.475003 h^6}};

sol = Part[NSolve[Det[%] == 0, h], 1]

Det[mat /. sol]

Det[mat] /. sol

As you suspected when you mentioned SetPrecision, you are encountering numerical errors, probably catastrophic loss of precision when calculating the determinant; your calculations do in fact need to be carried out at higher precision.

If possible, you would want to use exact numbers in your matrix, or take advantage of the arbitrary-precision capabilities of Mathematica. For instance, we can convert all machine-precision numbers to arbitrary-precision ones with a number of digits of precision equal to that of common machine-precision numbers on your machine using SetPrecision (see also $MachinePrecision in the documentation):

det = Det[SetPrecision[mat, $MachinePrecision]];

sol = NSolve[det == 0, h];

det /. sol // PossibleZeroQ



(* Out: 

 {True, True, True, True, True, True, True, True, True, True, True,

  True, True, True, True, True, True, True, True, True, True, True,

  True, True} 

*)

As you can see, all those values of $h$ do bring your determinant reasonably close to zero, within machine-precision approximations.

As already noted, this is a polynomial eigenproblem. First, let's use SeriesCoefficient to extract the coefficient matrices of your $h$-matrix:

coeffs = Table[SeriesCoefficient[mat, {h, 0, k}], {k, Max[Exponent[mat, h]], 0, -1}];

You can then use PolynomialEigenvalues to find the eigenvalues of your $h$-matrix:

eigs = PolynomialEigenvalues[coeffs];

To check the eigenvalues returned:

ListLinePlot[Log10[Abs[Det[mat /. h -> #] & /@ eigs]], 

             Axes -> None, Frame -> True, PlotRange -> All]

and we see that the small eigenvalues are computed accurately, but the large eigenvalues are not. (Why this is so, I still have to do some research on.)

So, as a sanity check, let us compute the eigenvalues a little differently:

eigs2 = Reverse[1/PolynomialEigenvalues[Reverse[coeffs]]];

ListLinePlot[Log10[Abs[Det[mat /. h -> #] & /@ eigs2]],

             Axes -> None, Frame -> True, PlotRange -> All]