TransWikia.com

Numeric solution of ODEs oscillates strongly and stops halfway?

Mathematica Asked by cchongb on May 19, 2021

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:

enter image description here

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.

enter image description here

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.

One Answer

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]

enter image description here

Another choice is to set a proper MaxStepSize e.g. MaxStepSize -> 3 10^4. (The corresponding timing is 1.38926. )

Answered by xzczd on May 19, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP