[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] &
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:
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.
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.
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.
-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 instead5.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=when the LHS contains no pattern and to use:=when the LHS contains a pattern. E.g.,slu :=...meanssluis recomputed (viaNSolve) every time you use/evaluate it. Wastes a lot of time. $\endgroup$ – Michael E2 Aug 6 '18 at 17:29