# Plotting of an error estimate of NDSolve

So, in an earlier post, (see here), I asked a question regarding a surface of section for a system of coupled non-linear ODEs. What I am after now, is the numerical precision of the results I get back.

By solving the following system:

sol1 = With[{μ = 1/3},
ParametricNDSolve[{x''[t] == -Surd[x[t]^(Numerator@μ),
Denominator@μ] + Surd[(y[t] - x[t])^(Numerator@μ),Denominator@μ],
y''[t] == -Surd[(y[t] - x[t])^(Numerator@μ), Denominator@μ],
y[0] == 0, y'[0] == i , x[0] == 0,
x'[0] == Sqrt[2 - y'[0]^2 - ((Surd[(x[0])^(Numerator@(μ + 1)),
Denominator@(μ + 1)]) + (Surd[(y[0] - x[0])^(Numerator@(μ + 1)),
Denominator@(μ + 1)]))/(μ + 1)]}, {x, y}, {t, 100}, {i},
WorkingPrecision -> MachinePrecision]];


I get back x[t] and y[t] as ParametricFunction with parameters {i}. Next thing, I would like to plot the logarithm Log[10,Rer[t]] of the Relative Error function, defined as:

$$Rer(t)=\left|\frac{\left(\frac{x'(t)^2+y'(t)^2}{2}+\frac{x(t)^{\mu+1}+\left[ y(t)-x(t) \right]^{\mu+1}}{\mu+1}\right)-\left(\frac{x'(0)^2+y'(0)^2}{2}+\frac{x(0)^{\mu+1}+\left[ y(0)-x(0) \right]^{\mu+1}}{\mu+1}\right)}{\frac{x'(t)^2+y'(t)^2}{2}+\frac{x(t)^{\mu+1}+\left[ y(t)-x(t) \right]^{\mu+1}}{\mu+1}} \right|$$

versus the Log[10,t]. To that end, I asked another question, with regard of the LogLogPlot command (see here) and the values of the axes.

At the moment, my Mathematica skills are not quite sharpened, but I nevertheless, gave it a shot to try and get back the corresponding plot. My attempt is the following:

With[{μ = 1/3},
Show[LogLogPlot[
Evaluate[
Table[{Abs[(((x'[0.1*(i - 1)][t])^2 + (y'[0.1*(i - 1)][t])^2)/ 2 +
((x[0.1*(i - 1)][t])^(μ + 1) + (y[0.1*(i - 1)][t] -
x[0.1*(i - 1)][t])^(μ + 1))/(μ + 1) - (((x'[0.1*(i - 1)][0])^2 +
(y'[0.1*(i - 1)][0])^2)/2 + ((x[0.1*(i - 1)][0])^(μ + 1) +
(y[0.1*(i - 1)][0] - x[0.1*(i - 1)][0])^(μ + 1))/(μ + 1)))/
(((x'[0.1*(i - 1)][t])^2 + (y'[0.1*(i - 1)][t])^2)/2 +
((x[0.1*(i - 1)][t])^(μ + 1) + (y[0.1*(i - 1)][t] -
x[0.1*(i - 1)][t])^(μ + 1))/(μ + 1))] /. sol1}, {i, 2}]], {t, 1,100}],
FrameTicks -> {{ChartingScaledTicks[{#*Log[10] &, #/Log[10] &}],
ChartingScaledFrameTicks[{#*Log[10] &, #/Log[10] &}]},
{ChartingScaledTicks[{#*Log[10] &, #/Log[10] &}],
ChartingScaledFrameTicks[{#*Log[10] &, #/Log[10] &}]}},
ImageSize -> 600, GridLines -> Automatic, PlotRange -> All,
Frame -> True, Exclusions -> None]] // Timing


It actually took some time to learn how to use Evaluate with Table etc. The code runs, but the results look like this:

which does not correspond to whats really happening there, because a) the error estimate for the blue solution is huge and b) the second solution (for i=2) has zero error.

I am sure that I am doing something wrong, but I cannot figure out what. Any help would be much appreciated.

Update

After losing any hope into making it work, I thought of keeping it simple. Therefore, I defined the function:

h[t_] := Abs[(((x'[t])^2 + (y'[t])^2)/
2 + ((x[t])^(μ + 1) + (y[t] - x[t])^(μ + 1))/(μ + 1) - (((x'[0])^2 +
(y'[0])^2)/2 + ((x[0])^(μ + 1) +
(y[0] - x[0])^(μ + 1))/(μ + 1)))/(((x'[t])^2 + (y'[t])^2)/2 +
((x[t])^(μ + 1) + (y[t] - x[t])^(μ + 1))/(μ + 1))]


then solved the system for just one initial condition:

sol11 = With[{μ = 1/3},
ParametricNDSolve[{x''[t] == -Surd[x[t]^(Numerator@μ),
Denominator@μ] + Surd[(y[t] - x[t])^(Numerator@μ),Denominator@μ],
y''[t] == -Surd[(y[t] - x[t])^(Numerator@μ), Denominator@μ],
y[0] == 0, y'[0] == 0 , x[0] == 0,
x'[0] == Sqrt[2 - y'[0]^2 - ((Surd[(x[0])^(Numerator@(μ + 1)),
Denominator@(μ + 1)]) + (Surd[(y[0] - x[0])^(Numerator@(μ + 1)),
Denominator@(μ + 1)]))/(μ + 1)]}, {x, y, h}, {t, 100},
WorkingPrecision -> MachinePrecision]];


where I included h at {x,y,h}. Then I tried to see if I get logical results from h. It is estimating some kind of integration error and therefore it should be kelp at low levels. I tried t=10 and got a very logical result back:

With[{μ = 1/3}, Evaluate[h[10] /. sol11]]
{6.7405*10^-6}


And then I tried to plot it:

With[{μ = 1/3},
Plot[First[Evaluate[{h[t]} /. sol11]], {t, 0, 10}, PlotRange -> All]]


which of course returned nothing.. Not even to mention logarithmic scale axes.

Could anyone provide some lights here?

You have an evident problem with PlotRange. I don't know why do you use Show but putting the one more PlotRange->All inside the LogLogPlot changes the situation a lot:
• Thank you. Yes you are right, it does change the behavior, but still the numbers returned are way too huge and the second initial condition, the yellow one, is not supposed to behave like this at all. I think I am not using correctly the functions from NDSolve to plot the error of Rer(t). – Mitscaype Jun 5 '17 at 6:21