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:
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,
α[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
. )
Answered by xzczd on May 19, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP