Mathematica Asked on June 11, 2021
I am trying to solve the following equations in Mathematica using NDSolve.
1.Given Data
(*Units*)
AU = 1.496*10^11;
Ms = 1.99*10^30;
pc = 3.086*10^16;
G = 6.67*10^(-11);
Subscript[K, B] = 1.39/10^23;
(*Given data*)
n0 = 10^19;
[Mu] = 2.33;
mp = 1.67/10^27;
(*Obtained Data using Given Data*)
[Rho]0 = n0*[Mu]*mp;
r0 = 5*AU;
2.Given functions
(*Given Functions*)
[Rho][r_] := [Rho]0/(1 + r/r0)^(2.2);
T[r_] := 1200*([Rho][r]/[Rho]0)^(0.1);
cs[r_] := ((5*Subscript[K, B]*T[r]/(3*[Mu]*mp)))^(0.5);
Menc[r_?NumericQ] :=
4*Pi*[Rho]0*r0^3 *NIntegrate[x^2/(1 + x)^(2.2), {x, 0, r/r0}];
3.Equations
eqns = {M'[t] ==
4*Pi*G^2*
M[t]^2*[Rho][
r[t]]/( cs[r[t]]^2 + ( r'[t]^2 + (r[t]*[Phi]'[t])^2))^(3/2),
r''[t] == -G*Menc[r[t]]/r[t]^2 + (r[t]*[Phi]'[t])^2/
r[t] - (4*Pi*G^2*M[t]*[Rho][r[t]])*
r'[t]/( cs[r[t]]^2 + ( r'[t]^2 + (r[t]*[Phi]'[t])^2 ))^(3/2),
[Phi]''[t] == -2*r'[t]*[Phi]'[t]/r[t] -
(4*Pi*G^2*M[t]*[Rho][r[t]])*[Phi]'[
t]/( cs[r[t]]^2 + ( r'[t]^2 + (r[t]*[Phi]'[t])^2 ))^(3/
2)};
4.Initial conditions
initcond = {M[0] == 0.38 Ms, r[0] == AU,
r'[0] == 7*10^3, [Phi]'[0] == (0.01*10^3)/r[0], [Phi][0] == 0};
5.NDSolve command
ClearAll[M, r, [Phi]];
s = NDSolve[{eqns, initcond}, {M, r, [Phi]},
{t, 0, 10^11}
]
The error generated is
NDSolve::ndsz: At t == 3.599524047826583`*^8, step size is effectively zero; singularity or stiff system suspected.
I have used Explicit Runge Kutta, StiffnessSwitching, Adams, and BDF method, still have not been able to overcome this error.
Most probably the issue is that in the first steps, the effective step size increases a lot as the motion is almost linear, but when the $r[t]$ starts decreasing again to a very small value, near $t=3.5*10^8$ maybe, the effective stepsize is to be decreased to a very small number again, which the AI is failing to do as going from stepsize~$10^9$(I checked the stepsize at this point used by mathematica) to say, $10^1$ or $10^2$ is maybe computationally difficult. I suspect something like this is going on.
How can I solve these equations here?
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP