1
$\begingroup$

t is a generated matrix with a parameter kz. Through solving Det[t]==0, the solutions, i.e. kz, are obtained. When I substitute any one solution in Det[t], the result is a large number, why? How to solve this problem? Many thanks. The codes are as following:

Clear["`*"]
s = 2;
m = 1;
th = Pi/4;
fi = Pi/6;
vh = 16;
mu = 11;
delta = 8;
HBAR = 1.05457266*10^(-34);
ME = 9.1093897*10^(-31);
ELEC = 1.60217733*10^(-19);
Kh = 2.11*10^(10);
vKh[1] := {0, 0, 0}
vKh[2] := {Kh, 0, 0};
vKh[3] := {-Kh, 0, 0};
vKh[4] := {0, Kh, 0};
vKh[5] := {0, -Kh, 0};
vKh[0] := {0, -Kh, 0};
vKh[i_] := vKh[Mod[i, 5]];
kc = Sqrt[2*ME*ELEC/HBAR^2];
ku := kc*Sqrt[mu + delta];
kd := kc*Sqrt[mu - delta];
a3 = {Pi/Kh, Pi/Kh, Sqrt[2]*Pi/Kh};
k := {-ku*Sin[th]*Cos[fi], -ku*Sin[th]*Sin[fi], kz};
vkz[i_] := 
  If[Mod[i, 5] != 0, {0, 0, (i - Mod[i, 5])/5*Kh*Sqrt[2]/(m + 1)}, {0,
     0, (i - Mod[i, 5] - 5)/5*Kh*Sqrt[2]/(m + 1)}];
f[i_, i_] := 
  Total[(k + vKh[i])^2] - ku^2 - kz^2 + (kz + Total[vkz[i]])^2;
f[i_, j_] := 
  If[i == j, f[i, i], 
   kc^2*vh*Total[
     Table[Exp[I*n*Total[(vKh[j] + vkz[j] - vKh[i] + vkz[i])*a3]], {n,
        0, m}]]];
t := Array[f, {5*s, 5*s}];
slu := Select[kz /. NSolve[Det[t] == 0, kz], 
   Re[#] >= 0 && Im[#] >= 0 &];
td[i_] := t /. kz -> slu[[i]];
Det[td[1]]
$\endgroup$
  • $\begingroup$ The number 10^209 is not from real life. Such type questions are scholasticism. $\endgroup$ – user64494 Aug 4 '18 at 4:58
  • $\begingroup$ Thank you very much! $\endgroup$ – user59546 Aug 4 '18 at 5:22
  • $\begingroup$ As a general guideline, especially when doing complex calculations like this, it's best if you can convert all your Machine-Precision nos. (those with decimals) to Exact nos. (look up "Machine-Precision Numbers" in the docs for more info.). For instance, your code outputs -1.18221*10^194 + 3.28171*10^192 I. But if you convert your four Machine-Precision numbers to Exact numbers (e.g., HBAR = Rationalize[1.05457266`50, 0]*10^(-34);), your output is instead 5.90304*10^192, indicating the imaginary component in the first result was likely an artifact of precision loss during calculation. $\endgroup$ – theorist Aug 4 '18 at 17:42
  • 1
    $\begingroup$ You have constants that differ by dozens of orders of magnitude; numerical results will be catastrophically bad. Please, work in natural units, where all parameters are of order one. $\endgroup$ – AccidentalFourierTransform Aug 6 '18 at 1:06
  • 2
    $\begingroup$ Tangential programming tip: See mathematica.stackexchange.com/questions/8829/… -- Basically it is normally better to use = when the LHS contains no pattern and to use := when the LHS contains a pattern. E.g., slu :=... means slu is recomputed (via NSolve) every time you use/evaluate it. Wastes a lot of time. $\endgroup$ – Michael E2 Aug 6 '18 at 17:29
7
$\begingroup$

[Edit notice: Rewritten, restructured and expanded, but basically still the same explanation.]

I think the talk of "not from real life," "ridiculously large," and "better to use exact" is a bit overblown, because there is nothing very wrong with the answer. In fact, it's quite good. We can show that the solutions have relative errors that range from 6*10^-15 to 4*10-13, which is typical of machine-precision solutions. Machine precision is fast, and if it is good enough, then you should use it.

The main issue seems to be understanding that the residuals of the polynomial Det[t], which are numbers from 10^194 to 10^197, are numbers that are quite close to zero:

Abs@*Det@*td /@ Range@7
(*
  {1.18266*10^194, 1.85298*10^194, 3.06391*10^194, 4.105*10^195, 
   2.68997*10^196, 4.17866*10^196, 2.3451*10^197}
*)

Compared to "real life" numbers they don't look close to zero, but in fact they are very close to zero. By "very," I mean in relation to the magnitude Det[t]. If you subtract a number from the floating-point number 1.5, the smallest difference you can make is 2.2 * 10^-16, or $MachineEpsilon. This is because of the discrete nature of floating-point numbers. Likewise, if you subtract a number from 2.2791*10^209, the smallest difference from zero you can make is 1.82498*10^193. That is the same as a single error in the last bit (cf. Note 1). The number 2.2791*10^209 is the largest term of the expanded polynomial Det[t] evaluated at kz -> slu[[1]]. On comparison, the value of Det[t] at slu[[1]] is about 5 * 10^-16 times the largest coefficient, which is pretty small. Now td[1] effectively evaluates Det[t /. kz -> slu[[1]]], which is not the same as Det[t] /. kz -> slu[[1]]. The magnitude 10^209 may not be the best choice to compare the residual with, but it suggests a context in which it might be appropriately considered close to zero. Let's investigate further.

On reason to be concerned is that such large coefficients tend to cause problems with evaluating polynomials (e.g. Funny behaviour when plotting a polynomial of high degree and large coefficients, NSolve for high degree univariate polynomials). Also, it has been known since Wilkinson that the form of the polynomial can perturb the roots significantly. Yet another is the size of the determinant: as a number, it is "not from real life," as @user64494 remarks in a comment.

None of these concerns are particularly troubling in the OP's problem. The evaluation at arbitrary inputs is not going to matter unless it interferes with finding the roots. The perturbation of the roots can be tested, and the roots slu can be compared with high-precision results. These will show that the results in slu are accurate to 12-14 digits, which is more than the physical constants on which the calculation is based. Finally, the largest numbers involved seem to be around 10^212 at most, although I cannot inspect all the internal computations. Still, for 16-digit computations in floating-point, that's well within their dynamic range of 10^-308 ~ 10^308.

One can rescale by choosing different units to make the numbers closer to 1 in magnitude, but it won't help one whit. It will scale the entries in the matrix t by the same factor and have no effect on the numerics since floating-point scales admirably within its limits (except for the effect of scaling a binary number system by a powers of ten, which effect is negligible on a stable problem). The solutions kz are approximately 1/L, where the length L is around 1 Å = 10 nm = 10^-10 m or less, which is about the lengths of the smallest molecular bonds. The dimensions of the entries of t are the same as 1/L^2, and therefore the dimensions of Det[t] are 1/L^20. So if we change the units from meters to Ångströms, we should divide the entries of t by 10^20. The results are scaled down, Det[t] is divided by 10^200, but all it does is shift the problem down the scale.

scaled = NSolve[
  Det[t/10.^20 /. kz -> 10^10 kz] == 0 && Re[kz] >= 0 && Im[kz] >= -10^-3, kz]
kz /. scaled
(*  {1.05945, 1.26717, 1.79276, 1.84599, 2.72977, 2.95544, 3.35895}  *)

(10^10 % - slu)/slu // Abs // Max
(*  3.24024*10^-13  *)

Det[t/10.^20 /. kz -> 10^10 kz] /. scaled // Abs
(*
  {2.42512*10^-7, 2.17557*10^-6, 5.72205*10^-6, 
   2.86104*10^-6, 0.00012207, 0.000610352, 0.00146484}
*)

Are some of these residuals small and some not small? Relative to what? In fact they are just as good as before relative to the size of the scaled determinant in an interval containing the roots (see Note 2). To get a deceptively small-looking residuals (< 10^-42), change the units to picometers. The point is to understand the significance of an error you have to consider it relative to the smallest possible change in the output and to the effect of a small change in the input. Notice also that the roots did change slightly, a relative shift up to 3.24024*10^-13. The relative error is consistent with the round-off error in computing the determinant, which will be shown later when we discuss the perturbation of the entries in the matrix t.

Let's explore the accuracy and stability of the OP's problem. First, let's compute the results as if t were exact, using Rationalize to convert the entries to exact fractions. We increase the working precision to 300 digits, which is probably about 250 more than we need, and compute the relative errors.

With[{t0 = Rationalize[t, 0]},
  obj[z_?NumericQ] := Det[t0 /. kz -> z]
  ];
s1 = FindRoot[obj[kz], {kz, #},
     PrecisionGoal -> 20, AccuracyGoal -> 0, WorkingPrecision -> 300] & /@ 
       SetPrecision[slu, 300];

(slu - (kz /. s1))/(kz /. s1) // Abs
(*
  {6.48381*10^-15, 1.30952*10^-14, 3.65987*10^-14, 4.09576*10^-13, 
   3.33989*10^-14, 2.34915*10^-14, 7.04122*10^-15}
*)

We see that the roots for kz in slu are estimated to have around 12 to 14 digits of accuracy. I suppose that should be enough, but let's examine the stability of the solution.

Below we show that when the entries of matrix t are perturbed by random amounts of decreasing magnitude, the solutions of Det[t] == 0 converge to the OP's slu, until the solutions get within 10^-13 of slu. This is a "round-off plateau," at which level the round-off errors in Det[t] perturb its roots. It should be noted that when n = -16, the matrix t is altered in only some of the entries by one unit-last-place (ulp).

SeedRandom[0];
Table[    (* max. rel. err. in perturbations 10^n, n = -8,...,-16 *)
  Table[  (* rel. err. in 20 sol'ns *)
    With[{perturbed = kz /. 
       NSolve[Det[t (1 + RandomReal[10^-n {-1, 1}, {10, 10}])] == 0 && 
        Re[kz] >= 0 && Im[kz] >= -10^-3, kz]},
     (perturbed - slu)/slu // Abs // Max
     ], {20}] // Max,
  {n, 8, 16}] //
 ListLogPlot[#, DataRange -> {8, 16}, Joined -> True, Mesh -> All] &

Mathematica graphics

Finally, it might be worth showing the difference between these:

Det[t] /. kz -> sol  (* symbolic det., then plug in *)
Det[t /. kz -> sol]  (* plug in, then numeric det.  *)

The first formula, which computes the symbolic determinant, a polynomial, has a couple of issues. First there is some round-off error in computing the coefficients. More importantly, the terms of the polynomial are of a similar size and theoretically should cancel out. Consequently the rounding error will be large compared to the theoretical result, which is zero. That's what always happens when you compute zero. It has nothing to do with the numbers being very large; it has to do with them being about the same size, having opposite signs, and adding up to a comparatively small sum. The second formula computes the determinant of a numeric matrix, which is better behaved numerically, at least in this case.

Let's examine the first and fourth solutions. The determinant is plotted at consecutive floating-point numbers near each solution. The determinant is calculate in each of the ways above (blue, gold resp.) as well as with high precision (green). The fourth root slu[[4]] has the worst error. The first root slu[[1]] is located where Det[t] is crossing zero, but the fourth is not.

With[{sol = slu[[1]], r = Range[-100, 100]},
 With[{e = 2^(Last@MantissaExponent[sol, 2] - 54) r}, (* consecutive floats *)
  ListLinePlot[{
    Det[t] /. Table[{kz -> sol + δ}, {δ, e}] // Re,
    Det /@ (t /. Table[Evaluate@{kz -> sol + δ}, {δ, e}]) // Re,
    Det@SetPrecision[t, 300] /. 
      Table[{kz -> SetPrecision[sol, 300] + δ}, {δ, e}] // Re},
   DataRange -> r[[{1, -1}]], PlotRange -> {-1.5*^194, 1.5*^194}]
  ]]

With[{sol = slu[[4]], r = Range[-5000, 5000]},
 With[{e = 2^(Last@MantissaExponent[sol, 2] - 54) r}, (* consecutive floats *)
  ListLinePlot[{
    Det[t] /. Table[{kz -> sol + δ}, {δ, e}] // Re,
    Det /@ (t /. Table[Evaluate@{kz -> sol + δ}, {δ, e}]) // Re,
    Det@SetPrecision[t, 500] /. 
      Table[{kz -> SetPrecision[sol, 500] + δ}, {δ, e}] // Re},
   DataRange -> r[[{1, -1}]], PlotRange -> {-0.5*^196, 0.5*^196}]
  ]]

The noise of the first formula (blue) is much worse than the second (gold). Further, the second tracks the high-precision determinant closely.

It's not clear whether it is worth improving the result slu of NSolve, but it can be done with machine precision using FindRoot and the second formula. The residual value of Det[t] is a few orders of magnitude closer to zero. One can also see that how many floating-point numbers lie between slu[[4]] and the more accurate root s2[[4]] corresponds to the graph (on the right above; see also Note 3). (If you don't like the warning messages and Quiet, add a reasonable AccuracyGoal, like AccuracyGoal -> -197.)

obj[z_?NumericQ] := Block[{kz}, Det[t /. kz -> z]];
s2 = Quiet@FindRoot[obj[kz], {kz, #}] & /@ slu;
Det[t /. s2[[4]]] // Abs
(*  2.03164*10^192  *)

(* how many FP numbers to shift slu[[4]] over to the left *)
Re[slu[[4]] - (kz /. s2[[4]])]/2^(Last@MantissaExponent[slu[[4]], 2] - 54)
(*  3964.  *)

Notes:

  1. Consider f[x_] := 0.5*^209 x^2 - 1*^209. Then the residuals of NSolve[f[x] == 0, x] will be ±1.82498*10^193, which seem quite large. Are the solutions returned by NSolve wrong? They are ±Sqrt[2.], which are as close to the true solutions ±Sqrt[2] as one can get using machine precision.

  2. Suppose you found the residual of a root to be 10^-7 but the function was 10^-7 Cos[t]. Is that a small residual? No, it's the maximum of the function.

  3. Throughout, one sometimes finds that the roots have small imaginary parts, roughly on the order of $MachineEpsilon relative to the magnitude of the root. These are caused by rounding errors in computing the determinant of the matrix t, which contains complex numbers. They should be chopped of with Re[], but for most part, I ignored them.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Let us compare $10^{200}$ with the the Avogadro constant $ 6.022 10^{23} $ and the Planck constant $ 6.626 10^{-34} $. and the atom size (a ten-billionth of a meter) . $\endgroup$ – user64494 Aug 5 '18 at 5:52
  • 1
    $\begingroup$ @user64494 Why? Specious. $\endgroup$ – Michael E2 Aug 5 '18 at 13:04
  • $\begingroup$ The numbers may not be ridiculously large, but the scale difference is certainly bad, relative to their precision. In this situation it will be essentially impossible to avoid truncation error. And Det is going to be more susceptible than most other functions (I think you and others noted this in an earlier version of this endlessly morphing thread). $\endgroup$ – Daniel Lichtblau Aug 5 '18 at 14:24
  • $\begingroup$ @DanielLichtblau It seems to me that in this particular problem, scale is not an issue. All terms of Det[t] would be scaled by the same factor with a change of units, and the relative error would be unaffected. I don't think you get an excessive loss of precision except near the zeros, which happens with polynomials. The roots stably converge under perturbation of t until we run into a round-off plateau when the relative perturbation is less than 10^-13. I don't think there's any numerical problem here except seeing that 10^194 == 10^209*10^-15 is a pretty good approximation to zero. $\endgroup$ – Michael E2 Aug 5 '18 at 17:56
  • $\begingroup$ @MichaelE2 Thank you very much! My codes are edited to solve a physical problem. Only when the parameters s and m are some cases (e.g. s=5,m=2), the NullSpace[t] gives a {} but Det[t]==0 (at least a very small number). Several experts give me some helpful methods (e.g. Rationalize, SetPrecision, Chop...), including you. These methods are valid when some case of s and m. However, when s=5 and m=2, these methods are invalid. This case has been described in my newest question. $\endgroup$ – user59546 Aug 6 '18 at 5:19
4
$\begingroup$

If I had to guess, it's because your numbers are ridiculously large

Look at this:

With[{d1 = Det[t] // Re // ComplexExpand},
 Plot[d1, {kz, -10^11, 10^11}, PlotRange -> {All, {Automatic, 10^209}}]
 ]

asd

Approximate numerical things don't work well all the well with huge numbers like that. We can try to scale it down and just look at the real part:

With[{d1 = Det[t] // Re // ComplexExpand},
 Plot[d1/10^207, {kz, 1.0594*10^10, 1.0595*10^10}]
 ]

and you can FindRoot on that:

FindRoot[Re@Det[t]/10^207, {kz, 1.059*10^10}]

{kz -> 1.05945*10^10}

And that worked almost to MachinePrecision:

Re[Det[t]]/10^207 /. {kz -> 1.0594478650070786`*^10}

1.82498*10^-14

But then you realize your numbers are ridiculous and so you scale that almost perfect result by 10^207:

1.82498*10^-14*10^207.

1.82498*10^193

And it looks like your result is terrible, when in fact it's just the scale of your system that is terrible.

Unless you're going to carry about 208 digits of precision even little round off errors will scale terribly.

I'd suggest working in atomic units, i.e. with:

HBAR = 1.;
ME = 1.;
ELEC = 1.;
Kh = 1.;

But I don't think this will save you, as your numbers are still on the order of 10^10

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you very much! $\endgroup$ – user59546 Aug 4 '18 at 5:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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