[Tex/LaTex] serious scientific calculating with small numbers, using different methods

calculationspgfmathtikz-pgf

I got stuck in a problem with pgfmath module.
When trying to calculate the following functions values the result that gets printed out is always zero. I guess the problem is an underflow of pgfmaths co-domain, and the mathematical engine cuts off the small values.

This is an minimal example of what i was trying to explain in the introduction:

\documentclass[tikz]{standalone}
\usepackage[fleqn]{amsmath}
% physical constants:
\pgfmathdeclarefunction{m0}{0}{%
    \pgfmathparse{4*pi*1e-7}%
}

\pgfmathdeclarefunction{K}{0}{%
    \pgfmathparse{m0*pi/4}%
}
\pgfmathdeclarefunction{c1}{0}{% c1 = K/1.45
    \pgfmathparse{K/1.45}%
}
\pgfmathdeclarefunction{c2}{1}{% c2(gamma) = c1/gamma^2
    \pgfmathparse{c1/(#1)^2}%
}
\pgfmathdeclarefunction{c3}{2}{% c3(gamma,lambda)
    \pgfmathparse{K/((#1)*((#2)+0.45))}%
}
\pgfmathdeclarefunction{DL_rel}{2}{% DL_rel(lambda, gamma)
    \pgfmathparse{(2*sqrt((580*#1+261)*#2^3)+40*#1+18)/(29*#2^3-20*#1-9)}%
}

% the problematic pgfmath-function
\pgfmathdeclarefunction{N1}{5}{% N1(soluition, Lmin, gamma, lambda, d1)
    \pgfmathsetmacro\numa{(DL_rel(#4,#3)*#2)/2+#2}%
    \pgfmathsetmacro\numb{sqrt(\numa^2-c1*c3(#3,#4))/(4*c2(#3)^2)}%
    \pgfmathsetmacro\denom{2*c1*#5}%
    \pgfmathparse{(#1 -1) ?%
        (sqrt((\numa + \numb)/\denom))%
        :%
        (sqrt((\numa - \numb)/\denom))%
    }%
}

\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}
    \pgfmathparse{N1(1, 4e-4, 1.2, 0.466, 0.115)}\pgfmathresult\\
    \pgfmathparse{N1(2, 4e-4, 1.2, 0.466, 0.115)}\pgfmathresult
    \end{array}\right.
    \end{equation*}
\end{document}

It took me already way to long to get so far and it would be very frustrating for me to stop here, and get all the math done in a different application.

Maybe someone out there has a solution to my problem.

Anyway i want use this post to say thank you! this has been a wonderful place to learn Latex so far, good job everybody!

Best Answer

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 of 1e-101 ... the reason appears simple enough that \numb is a square root of a difference, and somehow this difference rather than being found to be zero turns out be 1e-101 due to a rounding error in the last, 92th digit of each term, which probably are of the order of 1e-10 !!! Yes, this should explain it for all levels N of float precision. I guess sometimes the difference is zero, sometimes it gives 1e-(9+N). For example with N=20 one can expect a difference of 1e-29, hence about 3e-15 on 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;

Best Answer

update to re-insist on numerical instability

update to comment an intriguing numerical curiosity/instability

original answer

Related Solutions

[Tex/LaTex] calculating polar coordinates with tikz

[Tex/LaTex] Calculating with different types of variables

Related Question