Numeric solution of ODEs oscillates strongly and stops halfway？

Question

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, α[0] == -110, α'[0] == Sqrt[
    A1[0]], ψ'[0] == -1*10^-6*Sqrt[A1[0]] , χ[0] == 
    5*10^-4, χ'[0] == g*λ/f*10^-6 ((Ω^4 Sin[(5*10^-4)/f])/(
      3 g λ Sqrt[A1[0]]))^(2/3), ψ[
     0] == ((Ω^4 Sin[(5*10^-4)/f])/(
     3 g λ Sqrt[A1[0]]))^(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.

Answer 1

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,
     α[0] == -110, 
     α'[0] == Sqrt[A1[0]],
     ψ[0] == ((Ω^4 Sin[(5*10^-4)/f])/(3 g λ Sqrt[A1[0]]))^(1/3),
     ψ'[0] == -1*10^-6*Sqrt[A1[0]], 
     χ[0] == 5*10^-4, 
     χ'[0] == g*λ/f*10^-6 ((Ω^4 Sin[(5*10^-4)/f])/(3 g λ Sqrt[A1[0]]))^(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. )