5
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ Correction: I get the output 0.118714 - 0.0526506 I (as the second output) and 0.106201 - 0.0979004 I (as the third output); sorry, used a different matrix. But the problem still stands. $\endgroup$
    – Nils
    Commented Mar 10, 2019 at 3:30
  • 2
    $\begingroup$ This looks like a polynomial eigenvalue problem. $\endgroup$ Commented Mar 10, 2019 at 8:24

2 Answers 2

8
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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]

plot of eigenvalue error

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]

another plot of eigenvalue error

$\endgroup$
5
  • 1
    $\begingroup$ (Hopefully, when I can get my computer fixed, I'll be able to devote more time on this and other problems!) $\endgroup$ Commented Mar 11, 2019 at 6:43
  • $\begingroup$ Thanks, that's very useful! Though, does this imply that I can't assume the first few solutions actually solve the problem with reasonable accuracy? $\endgroup$
    – Nils
    Commented Mar 12, 2019 at 3:00
  • $\begingroup$ @Nils, yes, for some reason (which I still need to do further research on), the large eigenvalues are not being computed accurately by both procedures (i.e. the normal one and the one with reversed coefficients). But, let me suggest an experiment for you to try while I'm not at a computer: evaluate mh = Map[HornerForm, mat, {2}], and then try both ListLinePlot[Log10[Abs[Det[mh /. h -> #] & /@ eigs]], Axes -> None, Frame -> True, PlotRange -> All] and ListLinePlot[Log10[Abs[Det[mh /. h -> #] & /@ eigs2]], Axes -> None, Frame -> True, PlotRange -> All]. $\endgroup$ Commented Mar 12, 2019 at 3:09
  • $\begingroup$ The problem seems to be numerically unstable, as the characteristic polynomial has reasonably large order. Hence, one subtracts two relatively large numerical values, which results in a loss of precision especially for large roots. I will try to find the roots by searching for local minima of the characteristic polynomial and see what happens. $\endgroup$
    – Nils
    Commented Mar 12, 2019 at 15:20
  • $\begingroup$ @Nils, James Wilkinson has long warned about ever trying to do numerical stuff involving the characteristic polynomial, so a direct expansion into coefficients is out. $\endgroup$ Commented Mar 12, 2019 at 22:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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