TransWikia.com

step size effectively 0 error in ndsolve

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?

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