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:
\begin{equation} 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| \end{equation}
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 -> {{Charting`ScaledTicks[{#*Log[10] &, #/Log[10] &}],
Charting`ScaledFrameTicks[{#*Log[10] &, #/Log[10] &}]},
{Charting`ScaledTicks[{#*Log[10] &, #/Log[10] &}],
Charting`ScaledFrameTicks[{#*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?
