serious scientific calculating with small numbers, using different methods
You could try going the lualatex route and then have all the precision of the math in lua, so no pgf at all. The following may show the wrong answer, but that isn't due to the mathematical capabilities of lua, more likely just a rather hasty translation of the required formula:
\documentclass[preview,border=5]{standalone}
\usepackage[fleqn]{amsmath}
\usepackage{luacode}
\begin{luacode*}
pi = math.pi
sqrt = math.sqrt
m0 = 4 * pi * 1e-7
K = m0 * pi / 4
c1 = K / 1.45
c2 = function (x) return c1 / (x^2); end
c3 = function (x, y) return K / (x * (y + 0.45)); end
DLrel = function (x, y)
return (2 * sqrt((580 * x + 261) * y^3) + 40 * x + 18) /
(29 * y^3 - 20 * x - 9)
end
N1 = function (s, Lmin, g, l, d1)
dL = DLrel(l, g) * Lmin
nm = dL / 2 + Lmin
rt = (nm^2 - c1 * c3(g, l) * 0.25 * c2(g)^-2 * dL^2)
dn = 2 * c1 * d1
s = -(s % 2) * 2 + 1
return sqrt((nm + s * sqrt(rt)) / dn)
end
\end{luacode*}
\def\luaprint#1{\directlua{tex.print(#1)}}
\begin{document}
\begin{equation*}
N_{1_{1,2}} = \left\{
\begin{array}{l}
\luaprint{N1(1, 4e-4, 1.2, 0.44, 0.115)}
\\
\luaprint{N1(2, 4e-4, 1.2, 0.44, 0.115)}
\end{array}\right.
\end{equation*}
\end{document}
update to re-insist on numerical instability
Here is from the log output, using gamma=1.2, lambda=0.466 and Lmin=4e-6, which were the values from the original OP. Similar, but differing result are observed with gamma=1.2, lambda=0.44.
4: -1.00000e-13
5: 0
6: 0
7: 0
8: 0
9: 0
10: 2.00000e-19
11: 1.00000e-20
12: 2.00000e-21
13: 1.00000e-22
14: 1.00000e-23
15: 1.00000e-24
16: 0
17: 0
18: 0
19: 0
20: 1.00000e-29
21: 1.00000e-30
22: 0
23: 0
24: 1.00000e-33
25: 1.00000e-34
26: 0
27: 2.00000e-36
28: 0
29: 0
30: 0
31: -1.00000e-40
32: 2.00000e-41
33: 0
34: 2.00000e-43
35: 0
36: 1.00000e-45
37: 1.00000e-46
38: 0
39: -1.00000e-48
40: 0
41: 0
42: -1.00000e-51
43: 1.00000e-52
44: 1.00000e-53
45: 1.00000e-54
46: 0
47: 1.00000e-56
48: 1.00000e-57
49: 0
50: 1.00000e-59
51: 2.00000e-60
52: 0
53: 0
54: 2.00000e-63
55: 1.00000e-64
56: 2.00000e-65
57: 0
58: 0
59: 1.00000e-68
60: 1.00000e-69
61: 0
62: -1.00000e-71
63: 0
64: -1.00000e-73
65: 0
66: 0
67: 1.00000e-76
68: 0
69: 0
70: -1.00000e-79
71: 1.00000e-80
72: 0
73: 0
74: 0
75: 0
76: 0
77: 0
78: -1.00000e-87
79: 1.00000e-88
80: 1.00000e-89
81: 1.00000e-90
82: 1.00000e-91
83: 1.00000e-92
84: 1.00000e-93
85: 0
86: 0
87: 0
88: 1.00000e-97
89: 1.00000e-98
90: 1.00000e-99
91: 0
92: 1.00000e-101
Note: I have set constant K to value 1 to skip the computations with Pi. This modifies the result here as the rounding float operations are not exactly the same.
Source code to generate the above.
\documentclass{article}
\usepackage{xintexpr}% tested with 1.2e release
%\xintverbosetrue
\usepackage[fleqn]{amsmath}
\begin{document}
This is the equation im trying to solve:
\begin{equation} \label{eq:N1}
N_{1_{1,2}} =\sqrt{
\frac{\frac{\Delta L}{2} +L_{min} \pm \sqrt{\left(\frac{\Delta L}{2} + L_{min}\right)^2 - \frac{c_1 c_3}{4\cdot c_2^2}\cdot \Delta L^2}}
{2\cdot c_1 d_1}
}
\end{equation}
The $\Delta L/L_{min}$, $c_1$, $c_2$, $c_3$ are functions of some parameters
($L_{min}$ is only an overall scaling) and these functions are set-up in such
a way that actually the square root at the numerator of the fraction vanishes
exactly, independently of the values of the parameters. Numerically though it
may be found non zero and this will then induce a catastrophic drop in
precision in the final result. And it may even happen that it will be found
negative, thus potentially raising an error in the complete evaluation.
See the compilation log for this problematic $\left(\frac{\Delta L}{2} +
L_{min}\right)^2 - \frac{c_1 c_3}{4\cdot c_2^2}\cdot \Delta L^2$ evaluated
with the float precision set to various values. You will see it turns negative
at times.
% \xintdeffloatvar pi:=
% 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342;
\xintFor* #1 in {\xintSeq{4}{92}}\do
{
\xintDigits := #1;
% constants
% \xintdeffloatvar m0:= 4pi*1e-7;
% \xintdeffloatvar K := m0*pi/4;
\xintdeffloatvar K := 1;
\xintdeffloatvar c1:= K/1.45;
% functions
\xintdeffloatfunc c2(u) := c1/u^2;
\xintdeffloatfunc c3(u,v):= K/(u(v+0.45));
% attention arguments like in OP-update, permuted compared to the OP-original
\xintdeffloatfunc DL_rel(u,v):= (2sqrt((580v+261)u^3)+40v+18)/(29u^3-20v-9);
\xintdeffloatfunc DL(t,u,v) := DL_rel(u,v)*t;
% Notice that t=Lmin acts only as an overall scaling factor.
\xintdeffloatfunc numbsquared(t,u,v):=
subs((Z/2+t)^2-c1*c3(u,v)/(4*c2(u)^2)*Z^2, Z=DL(t,u,v));
\typeout{#1: \xintthefloatexpr [6] numbsquared (4e-6, 1.2, 0.466)\relax }
}
\end{document}
I have also tested with Maple. Again setting K to 1.
numbsquared := proc (N)
local m0, K, c1, c2, c3, DL_rel, DL, localnumbsquared;
Digits:=N;
# m0 := 4*Pi*1e-7;
# K := m0*Pi/4;
K := 1;
c1 := K/1.45;
c2 := u->c1/u^2;
c3 := (u,v)->K/(u*(v+0.45));
DL_rel := (u,v)->(2*sqrt((580*v+261)*u^3)+40*v+18)/(29*u^3-20*v-9);
DL := (t,u,v)->DL_rel(u,v)*t;
localnumbsquared := (t,u,v)->subs(Z=DL(t,u,v),(Z/2+t)^2-c1*c3(u,v)/(4*c2(u)^2)*Z^2);
return localnumbsquared(4e-6, 1.2, 0.466)
end proc:
for N from 4 to 92 do printf("%2d, %e\n", N, numbsquared(N)) end do;
The results are of the same type but with both coincidences and differences.
4, 0.000000e+00
5, 0.000000e+00
6, 0.000000e+00
7, 0.000000e+00
8, 0.000000e+00
9, 1.000000e-18
10, 2.000000e-19
11, 1.000000e-20
12, -1.000000e-21
13, 1.000000e-22
14, -1.000000e-23
15, 1.000000e-24
16, 0.000000e+00
17, 0.000000e+00
18, 0.000000e+00
19, 0.000000e+00
20, 1.000000e-29
21, 1.000000e-30
22, 0.000000e+00
23, 1.000000e-32
24, -1.000000e-33
25, 1.000000e-34
26, 0.000000e+00
27, 2.000000e-36
28, 1.000000e-37
29, 0.000000e+00
30, 0.000000e+00
31, -1.000000e-40
32, 2.000000e-41
33, 0.000000e+00
34, 2.000000e-43
35, -1.000000e-44
36, 1.000000e-45
37, 1.000000e-46
38, 0.000000e+00
39, -1.000000e-48
40, 0.000000e+00
41, 0.000000e+00
42, -1.000000e-51
43, 1.000000e-52
44, -1.000000e-53
45, 1.000000e-54
46, 0.000000e+00
47, -1.000000e-56
48, 1.000000e-57
49, 0.000000e+00
50, 1.000000e-59
51, 0.000000e+00
52, 0.000000e+00
53, 0.000000e+00
54, 2.000000e-63
55, 1.000000e-64
56, 0.000000e+00
57, 0.000000e+00
58, 0.000000e+00
59, 1.000000e-68
60, 1.000000e-69
61, 0.000000e+00
62, -1.000000e-71
63, 0.000000e+00
64, -1.000000e-73
65, 0.000000e+00
66, 0.000000e+00
67, 1.000000e-76
68, -2.000000e-77
69, -1.000000e-78
70, -1.000000e-79
71, 1.000000e-80
72, 0.000000e+00
73, 0.000000e+00
74, 0.000000e+00
75, 0.000000e+00
76, 0.000000e+00
77, 0.000000e+00
78, 0.000000e+00
79, 1.000000e-88
80, -1.000000e-89
81, -1.000000e-90
82, 1.000000e-91
83, 0.000000e+00
84, -1.000000e-93
85, 0.000000e+00
86, 0.000000e+00
87, 0.000000e+00
88, -1.000000e-97
89, -2.000000e-98
90, 1.000000e-99
91, -2.000000e-100
92, 1.000000e-101
update to comment an intriguing numerical curiosity/instability
First version of OP asked for N1(1, 4e-6, 1.2, 0.466, 0.115), which is treated below. Turns out apparently that the \numb in that case is zero. But depending on the float precision it may be found or not to be zero. The \numa being only about 1e-5 the precision of final result may be drastically reduced by small but non zero \numb.
I compared xint with maple and obtained similar results (differ only in last digit) for 16, 20, 24 digits of precision (takes time to use maple, more testing on the way). For Pi, I did my tests in xint starting with 94 decimals and multiplying first by 1.0 to reduce to stated \xintDigits. Thus I ran the code below with
\xintDigits := 16; % or 20, 24, 28, ...
\xintdeffloatvar pi:= 1.*3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342;
On the Maple side, Pi is treated as a symbol until the final evalf.
Here is how the results look like on the xint side with higher and higher float precision:
8.038962683509860
8.0389626847662074875
8.03896268352242165100730
8.038962683509858157707745289
8.0389626835098594140570752438847
8.03896268350985817027123858823346261
8.038962683509858157707745288681429203705
8.0389626835098581577090016380113844070447859
8.03896268350985815770775785217472875573691848736
(this is for the solution with \numa+\numb). Notice how the result with 16 digits is much better than the one with 20 digits!!!!!!! and it is even quite better than the one for 24 digits! (but not as good as 28 digits). This is due to the fact that at 16 digits, both xint and maple find \numb to be zero, but about 3e-15 at 20 digits which induces a big error in the sum, as \numa is about 1e-5.
With 92 digits of precision one finds \numb about 3e-51 numerically. It the exact value is zero, this means that it spoils the digits of the result after about 46 of them...
With 92 digits of precision, Maple finds \numb to be
> evalf(Q(4e-6, 1.2, 0.466, 0.115));
0.3162277660168379331998893544432718533719555139325216826857504852792594438\
-50
6392382213442481084 10
And xint obtains
3.1622776601683793319988935444327185337195551393252168268575048527925944386392382213442481084e-51
you can see it matches up to the end. ;-)
update what sort of fool I am! How many years before recognizing the square root of
10???? notice that the above is essentially the square root of1e-101... the reason appears simple enough that\numbis a square root of a difference, and somehow this difference rather than being found to be zero turns out be1e-101due to a rounding error in the last, 92th digit of each term, which probably are of the order of1e-10!!! Yes, this should explain it for all levelsNof float precision. I guess sometimes the difference is zero, sometimes it gives1e-(9+N). For example withN=20one can expect a difference of1e-29, hence about3e-15on square root which is exactly what is observed. Strangely the numerics never seem to give a difference of-1e-(9+N)which would raise on error on square root.
As the returned values decrease with increasing precision perhaps I can trust the exact one is zero (I have not done the algebra). If the exact value is really zero, adding or subtracting the above to something which is 0.000010117182975... will corrupt it after about 46 significant digits, destroying the 92 digits float evaluations to get it...
Very surprising ! (but do read quoted block above)
This should be taken into account when comparing with any other math engine: the formula is numerically instable due to possible catastrophic cancellation in a subtraction.
original answer
Here is an approach using another math engine. It knows only square root, but this is enough here. Note that in the given example \numb turns out to be exactly zero.
\documentclass[tikz]{standalone}
\usepackage{xintexpr}% tested with 1.2e release
\usepackage[fleqn]{amsmath}
% constants
\xintdeffloatvar pi:= 3.14159265358979323846;
\xintdeffloatvar m0:= 4pi*1e-7;
\xintdeffloatvar K := m0*pi/4;
\xintdeffloatvar c1:= K/1.45;
% functions
\xintdeffloatfunc c2(x) := c1/x^2;
\xintdeffloatfunc c3(x,y):= K/(x(y+0.45));
\xintdeffloatfunc DL_rel(u,v):=
(2*sqrt((580*u+261)*v^3)+40*u+18)/(29*v^3-20*u-9);
% This is allowed by xint parser also (tacit multiplications):
% \xintdeffloatfunc DL_rel(u,v):= (2sqrt((580u+261)v^3)+40u+18)/(29v^3-20u-9);
% Of course we could simplify here by defining more intermediate functions.
% We could define "numa" and "numb" functions, and set them up as functions
% of an already computed "DL_rel" which serves in both.
% It is possible to use the "subs(expression, x=...)" syntax.
% Limitation is that the dummy parameter must be a single letter.
% Also, the inner-most subs will have the last defined thing, and the
% outer-most subs the first defined thing.
\xintdeffloatfunc N1(a,t,u,v,w):=
subs(subs(subs(subs(
if(a=1, sqrt((P+Q)/D), sqrt((P-Q)/D)),
% debugging because something is strange with Q = \numb which is zero
% (P, sqrt(c1*c3(u,v))/c2(u)*X ),
% well after all it was CORRECT that Q was zero with these numerics
Q = sqrt(P^2-c1*c3(u,v)/(c2(u)^2)*X^2)% =\numb,
),
P = X+t % P=\numa, and I think t is Lmin
),
X = DL_rel(v,u)*t/2 % X= DeltaL/2
),
D = 2c1*w % D=\denom
)% must use single letters in subs
;%
\begin{document}
This is the equation im trying to solve:
\begin{equation} \label{eq:N1}
N_{1_{1,2}} =\sqrt{
\frac{\frac{\Delta L}{2} +L_{min} \pm \sqrt{\left(\frac{\Delta L}{2} + L_{min}\right)^2 - \frac{c_1 c_3}{4\cdot c_2^2}\cdot \Delta L^2}}
{2\cdot c_1 d_1}
}
\end{equation}
This should print out its solution:
\begin{equation*}
N_{1_{1,2}} =
\left\{ \begin{array}{l}
\xintthefloatexpr N1(1, 4e-6, 1.2, 0.466, 0.115)\relax\\
\xintthefloatexpr N1(2, 4e-6, 1.2, 0.466, 0.115)\relax
\end{array}\right.
\end{equation*}
\end{document}
In this example, the two solutions are the same, as the \numb vanishes ...
Notice that the solutions are computed expandably. That matters to some (fools...).
If you want more precision start with:
\xintDigits := 32;
\xintdeffloatvar pi:= 3.141592653589793238462643383279503;