0

Having studied the answers to this question on how to control the error of NDSolve, I then tried to apply it to the following system of second-order nonlinear ODEs:

lowsol = With[{μ = 1/3}, 
NDSolve[{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}, {t, 0, 10}, 
   Method -> "ExplicitRungeKutta", WorkingPrecision -> 61, 
   InterpolationOrder -> All, MaxSteps -> 2*^6, 
   StartingStepSize -> 1*^-8, MaxStepSize -> 1*^-4]] // Timing

which returns a pair of Interpolating Function. The system has a first integral of motion, which is set by its energy $E$, given by: \begin{equation} E=\frac{x'^2+y'^2}{2}+V(x,y) \end{equation}

I then tried to plot this energy and see if it remains a constant in the following way:

With[{μ = 1/3}, 
Plot[{(((x'[t])^2 + (y'[t])^2)/2 + (3/4)*((x[t])^(4/3) + 
(y[t] - x[t])^(4/3)))} /. lowsol // RealExponent // Evaluate, {t, 0, 10}, 
PlotRange -> All]]

but what was returned, is the following error:

ReplaceAll::rmix: Elements of {513.969,
{{x->InterpolatingFunction[{{0,10.0000000000000000000000000000000000000000000000
0000000000000}},{5,1,7,{327114},
{4},{InterpolatingFunction[<<5>>]},0,0,0,Automatic,{},{},False},
{{0,5.086860500065041295097264042754353058189185580664200997582414*10^-
20,<<47>>,5.890559344799673487983779799641879661003747176870946110552147*10^-
15,<<327064>>}},
{{0,1.414213562373095048801688724209698078569671875376948073176680},<<49>>,
<<327064>>},{{{<<327115>>},
{<<327115>>}}}],y->InterpolatingFunction[{{0,10.00000000000000000000000000000000
000000000000000000000000000}},<<3>>,{{{<<327115>>},{<<327115>>}}}]}}} are a 
mixture of lists and nonlists.

which states that there are elements which are a mixture of lists and nonlists. I would like to ask the following:

  • Why did this error occur and how can I avoid it?
  • How is it possible to reduce the computation time? In the question linked above, the computation time is significantly less than the one I need to acquire the results for my system.
  • Is it possible to do the same for a number of initial conditions? I tried ParametricNDSolve but the plotting issue would remain.

Update

After running the piece of code provided by @zhk, I found out that the numerical integration returns a huge error estimate for t=10^4 time units as it can be seen below: enter image description here

How can I reduce this error?

Bazinga
  • 737
  • 4
  • 15
  • You're setting lowsol equal to the result of Timing. You probably want to set it equal to the result of NDSolve. Try moving it inside the With: With[{..}, lowsol = NDSolve[..]] // Timing. – Michael E2 Jun 26 '17 at 23:04
  • @MichaelE2 I will try it now and let you know. – Bazinga Jun 27 '17 at 09:05
  • It seems like you should be using the Projection or SymplecticPartitionedRungeKutta methods for NDSolve, which explicitly check that the invariant is conserved during the integration process. But I'm not getting that to work either. – SPPearce Jun 27 '17 at 09:43
  • @KraZug I read about those methods on the documentation but I was not able to make it run, though I have not idea why. – Bazinga Jun 27 '17 at 09:46
  • I integrate it with SPRK, but the invariant is not conserved: \[Mu] = 1/3; ode1 = x''[t] == -x[t]^\[Mu] + (-x[t] + y[t])^\[Mu]; ode2 = y''[t] == -(-x[t] + y[t])^\[Mu]; ics = {y[0] == 0, y'[0] == 0, x[0] == 0, x'[0] == Sqrt[2]}; invariant = (x'[t]^2 + y'[t]^2)/ 2 + (3/4)*(x[t]^(4/3) + (y[t] - x[t])^(4/3)); projerksol = NDSolve[Flatten[{ode1, ode2, ics}], {x, y}, {t, 0, 1}, Method -> {"SymplecticPartitionedRungeKutta", "DifferenceOrder" -> 2, "PositionVariables" -> {x, y}}, StartingStepSize -> 0.001][[1]]; Plot[Im[invariant /. projerksol], {t, 0, 1}] – SPPearce Jun 27 '17 at 09:48
  • Can you explain how you get the invariant? – SPPearce Jun 27 '17 at 09:50
  • @KraZug The kinetic energy is always $ \sum (x_i^2+y_^2)$ while for the potential I just used the existence of a gradient: $F=-\nabla V$, where $V=V(x,y)$ is a scalar function. – Bazinga Jun 27 '17 at 10:02
  • @KraZug It does conserve for the first t=10^3 time steps with only 10^-5 relative error but then it goes out of hand. Long time numerical integration seems to pose a problem – Bazinga Jun 27 '17 at 10:13
  • Note that there is an imaginary part that grows (the plot command that I included). – SPPearce Jun 27 '17 at 10:14
  • @KraZug Why would that play a role? I am using Surd to acquire only the corresponding real roots. – Bazinga Jun 27 '17 at 10:18

1 Answers1

1

Something like this?

μ = 1/3;

ode1 = x''[t] == -x[t]^μ + (-x[t] + y[t])^μ;

ode2 = y''[t] == -(-x[t] + y[t])^μ;

lowsol = NDSolve[{ode1, ode2, y[0] == 0, y'[0] == 0, x[0] == 0, 
   x'[0] == Sqrt[2 - (x[0]^(1 + μ) + (-x[0] + y[0])^(1 + μ))/(1 + μ) - 
   y'[0]^2]}, {x, y}, {t, 0, 10}, Method -> {"FixedStep", Method -> "ExplicitEuler"}];

xy = Flatten /@Table[{t, Re[{(((x'[t])^2 + (y'[t])^2)/
      2 + (3/4)*((x[t])^(4/3) + (y[t] - x[t])^(4/3)))} /. First@lowsol]}, {t, 0, 10, 0.1}];

ListLinePlot[xy, PlotRange -> All]

enter image description here

zhk
  • 11,939
  • 1
  • 22
  • 38
  • Looks good. Perhaps Flatten was what I was looking for? But why "ExplicitEuler"? Also, from the plot I can see an error estimate of 10^-3 for the first 10 time steps. I would like to bound it under 10^-4 since I am running 10^4 time steps integration. How can I apply the options: WorkingPrecision -> 61, InterpolationOrder -> All, MaxSteps -> 2*^6, StartingStepSize -> 1*^-8, MaxStepSize -> 1*^-4 as above? – Bazinga Jun 26 '17 at 17:27