# NDSolve: why is the solution error reaching a limit larger than machine precision?

Let us take the following differential equation example: $$\dot u^2=\mu^2(1-u^2)(1-m u^2),\quad u(0)=0.$$

According to Wikipedia, the analytical solution can be established with Jacobi elliptic functions as: $$u(t) = \mathrm{sn}\left(\mu\,t,m\right)$$

For a fixed step size $$h$$, I wish to follow the $$\ell_\infty$$ absolute error evolution, between the NDSolve solution and the analytical one, as $$h$$ step size decreases.

My test is using a fifth-order Bogacki-Shampin Runge-Kutta variant which is readily available through NDSolve as described in Mathematica's manual.

The following function evaluates the error:

errors[stepSize_] :=
Module[{m = 0.99999,
\[Mu]sq = 25,
tf = 1.4},
h = stepSize*tf;
uEx[t_] := JacobiSN[Sqrt[\[Mu]sq]*t, m];

(*Bogacki-Shampine 5th order*)
BS54 = {"ExplicitRungeKutta",
"Coefficients" -> "EmbeddedExplicitRungeKuttaCoefficients",
"DifferenceOrder" -> 5, "StiffnessTest" -> False};
uNum[t_] =
NDSolve[{(u'[t])^2 == \[Mu]sq (1 - u[t]^2) (1 - m*u[t]^2),
u[0] == 0}, u[t], {t, 0, tf},
Method -> {"FixedStep", Method -> BS54},
StartingStepSize -> h][[2, 1, 2]];
uExact = uEx[Range[0, tf, h]];
uNumeric = Re@uNum[Range[0, tf, h]];
error = Max[Abs[uExact - uNumeric]];
Clear[uEx, BS54, uNum, uExact, uNumeric];
{error}
]


Then, when running the desired test, for the step size going from $$h=1e-1$$ down to $$h=1e-6$$, we see that the error seems to attain a limit value which is above machine precision. More specifically, when plotting these values we see that the method is converging very fast at the beginning but then its rate of convergence tends to zero in a portion of the curve:

Values for the above curve are obtained with these instructions

initPower = 1;
endPower = 6;
errVal =
Table[Join[{10^(-i)}, errors[10^(-i)]], {i, initPower, endPower,
0.25}]


and plotted with

ListLogLogPlot[
Transpose[{errVal[[;; , 1]], errVal[[;; , 2]]}],
PlotLabel -> "max|\!$$\*SubscriptBox[ StyleBox[\"u\",\nFontSlant->\"Italic\"], \ \(exact$$]\)-\!$$\*SubscriptBox[ StyleBox[\"u\",\nFontSlant->\"Italic\"], \(numeric$$]\)|",
AxesLabel -> {"step size",
"\[LeftDoubleBracketingBar]e\!$$\*SubscriptBox[\(\ \[RightDoubleBracketingBar]$$, $$\[Infinity]$$]\)"},
Joined -> True,
Mesh -> Automatic,
PlotMarkers -> Automatic,
ImageSize -> Large
]


The questions are:

• Why is the error not reaching machine precision?
• Is the behavior in the curve typical of the BS54 method?

• On simpler examples, like approximating $$\pi \cos u$$ which is the solution to $$\ddot u + u =0,\ u(0)=\pi,\ \dot u(0)=0$$, the error does seem to reach machine precision; at least it gets much closer to it. More importantly, there are no strange jumps in the error (log-log) curve.

• It converges to machine precision for $t$ up to about $0.71428\pm$. Then there's a sudden jump, it seems. Don't know why. Commented Jun 17 at 4:23
• @MichaelE2 Do you think it has to do with other internal functions being truncated below machine precision? Or is it the effect of accumulated rounding errors? Commented Jun 17 at 8:15
• My guess is slight stiffness. The errors oscillate in sign after 0.714 up until 1. or a bit later, when it settles down somewhat. Commented Jun 17 at 13:20
• This, too: Module[{m = 0.99999100, \[Mu]sq = 25}, Precision[ Sqrt[\[Mu]sq (1 - u[t]^2) (1 - m*u[t]^2)] /. u[t] -> 0.998100] - 100 ] shows more than a 2-digit loss of precision around the point where the error jumps. And if you raise WorkingPrecision, you get a consistent decrease in error as step size decreases. Still not sure how to explain the apparent suddenness of the jump in error at machine precision at a value of $t$ that seems not to depend on $h$. Commented Jun 17 at 14:23

The solution of the pendulum equation for $$x=\sin \phi(t)$$ is not in general the $$\text{sn}$$ function, that indeed has the differential

$$\text{sn}' = \sin(\text{am}(\omega t,m))' = \omega\ \text{cn}\ \text{am}' = \omega\ \text{cn}\ \text{dn} = \omega \ \sqrt{1-\text{sn}^2}\ \sqrt{(1- m \ \text{sn}^2)} ; \quad m=k^2$$

This is true only for $$0<=m <1.$$ It's the case of a wobbling rotation.

In the rotational mode the kinetic energy $$\frac{\omega^2}{2}$$ at $$am=0$$ must be larger than the potential energy $$1- cn=2$$ at the top position.

For $$m=1$$ the solution is an $$\tan^{-1}$$ and for m>1 its a more complicated elliptic function, where the velocity function $$\text {dn}(\omega t, \frac{1}{\omega^2})$$ is replaced by a $$\text {cn}( , \omega^2 )$$ that has a more complicated integral.

The parameters have to be specified for the different modes and the NIntegrate works only for a quarter period of monotony for $$x=\sin \phi$$, as it is the case always for a first order equation with square roots and zeros of the velocity expression.

$$dt = \frac{ \mu \ du}{\sqrt{1-u^2}\ \sqrt{1-m ^2 u^2}}$$

At $$\frac{du}{dt}=0$$, time stalls.

• I have been giving some thought to your answer. In the studied case, $m<1$ and the system is not the pendulum. Do you mean that the issue boils down to NIntegrate stalling at some point? (And so it affects sn which calls it?). Or is it that the issue comes from the sn function not being an appropriate exact solution? Commented Jul 2 at 13:05