# Numeric solution of ODEs oscillates strongly and stops halfway？

I tried to reproduce the results of a paper.

The differential equations are given by

A1[t_] = 1/
3 (1/2 χ'[t]^2 + Ω^4 (1 + Cos[χ[t]/f]) +
3/2 ((ψ'[t] + α'[t] ψ[t])^2 + g^2 ψ[t]^4 ));
A2[t_] = α''[t] +
1/2 χ'[t]^2 + (ψ'[t] + α'[t] ψ[t])^2 +
g^2 ψ[t]^4 ;
C1[t_] = χ''[t] +
3 α'[t] χ'[t] - Ω^4/f Sin[χ[t]/f] +
3 g λ/f ψ[t]^2 (ψ'[t] + α'[t] ψ[t]);
P1[t_] = ψ''[t] +
3 α'[t] ψ'[
t] + (α''[t] + 2 α'[t]^2) ψ[t] +
2 g^2 ψ[t]^3 - g λ/f ψ[t]^2 χ'[t];


by using NDSolve with initial conditions

s1 = NDSolve[{A2[t] == 0, C1[t] == 0,
P1[t] == 0, α == -110, α' == Sqrt[
A1], ψ' == -1*10^-6*Sqrt[A1] , χ ==
5*10^-4, χ' == g*λ/f*10^-6 ((Ω^4 Sin[(5*10^-4)/f])/(
3 g λ Sqrt[A1]))^(2/3), ψ[
0] == ((Ω^4 Sin[(5*10^-4)/f])/(
3 g λ Sqrt[A1]))^(1/
3)}, {α, χ, ψ}, {t, 0, 8*10^10},
MaxSteps -> 20000000]


The parameters are

g = 2.0*10^-6; λ = 200; f = 0.01;
Ω = 3.16*10^-4;


The phase graph of

ParametricPlot[
Evaluate[{ψ[t], \!$$\*SubscriptBox[\(∂$$, $$t$$]$$(ψ[ t] Exp[α[t]])$$\) Exp[-α[t]] 10^3} /. s1], {t,
0, 10000000000}, PlotPoints -> 1000, PlotRange -> All]


shown in the paper is very nice: However, first I can not reproduce the graph. The program ceases at 10^9 because of the singularity. So I changed the parameter f from 0.01 to 0.03. It also ceases at 10^10 and the phase graph looks not so good because of its oscillation. It's confusing me. Why their graph looks so nice without oscillation? Why I had singularity even when I set the same conditions as them.

Seems that the default error estimation of NDSolve doesn't work well for your initial value problem (IVP), and this turns out to be a (relatively) rare case that AccuracyGoal option helps. With f = 0.01:

s1 = NDSolve[{A2[t] == 0, C1[t] == 0, P1[t] == 0,
α == -110,
α' == Sqrt[A1],
ψ == ((Ω^4 Sin[(5*10^-4)/f])/(3 g λ Sqrt[A1]))^(1/3),
ψ' == -1*10^-6*Sqrt[A1],
χ == 5*10^-4,
χ' == g*λ/f*10^-6 ((Ω^4 Sin[(5*10^-4)/f])/(3 g λ Sqrt[A1]))^(2/3)},
{α, χ, ψ}, {t, 0, 10^10},
MaxSteps -> Infinity, AccuracyGoal -> 16]; // AbsoluteTiming
(* {4.46549, Null} *)

ParametricPlot[
Evaluate[{ψ[t], D[ψ[t] Exp[α[t]], t]*Exp[-α[t]] 10^3} /. s1],
{t, 0, 10^10}, PlotRange -> All, AspectRatio -> 1, PlotPoints -> 400] Another choice is to set a proper MaxStepSize e.g. MaxStepSize -> 3 10^4. (The corresponding timing is 1.38926. )

• Thank you so much!! I tried to set the AccuracyGoal->8 before but still failed. I wonder how to detemine the accuracy? I mean, why '16', not '50' or other number ? Apr 9, 2021 at 4:22
• And what causes the oscillation? Apr 9, 2021 at 4:24
• @cc I myself have just determined the option value by trial and error, but Michael E2 is good at systematic analysis of this stuff, you may check his excellent posts here (don't miss the linked posts): mathematica.stackexchange.com/q/118249/1871 As to the oscillation, it's because the default adaptive step choosing of NDSolve doesn't work well on your IVP, and I just noticed MaxStepSize can also be used, see my update. Apr 9, 2021 at 4:36
• I'm trying my best to understant it. There's still lots of thing to learn. Anyway, thanks again! Apr 9, 2021 at 4:48