# NDSolve Maxsteps Reached Error

I'm trying to solve a system of Equations. The code is commented below, I am mainly conserved with the last few lines, specifically the NDSolve line. No matter what I do to the Maxsteps option, the limit is reached at the same point t=16.486...

I'm not sure why NDSolve stops at this point. Any help would be greatly appreciated!

Edit: I've attempted to use AccuracyGoal and PrecisionGoal, they've stopped the error of maxsteps, but when I plot r(t) it only plots up until the time limit above.

(*all quantities are in their natural units, per unit kg*)
(* angular momentum of star in J (conserved i.e constant) *)
L = 3.086*10^25;
(* V is a constant in km/s *)
V = 200*10^3;
(* initial radial velocity (in polar coordinates) *)
vr = 50*10^3;
(* initial tangential velocity (in polar coordinates) *)
vt = 10^5;
r0 = 10*3.086*10^19;
(* seconds in a mega year (units of time I want to plot in *)
megayear = 3.154*10^7*10^6;
(* energy of the star system calculated from initial conditions \
(conserved i.e constant energy), the V^2 Log(r0) term is the \
spherical potential *)
energy = vr^2/2 + L^2/(2*r0^2) + V^2*Log[r0];
r =.
(* solve for the max and min radii of the orbit *)
Solve[L^2/(2 r^2) + V^2 Log[r] == energy, r]

(* Problem = If I increase the range over t (say, 100 instead of 10) \
I get a maxstep error *)
s = NDSolve[{r'[t] ==
Sqrt[2*(energy - V^2*Log[r[t]] - L^2/(2*r[t]^2))*
megayear^2], \[Theta]'[t] == L/r[t]^2*megayear, \[Theta][0] ==
vt/r0, r[0] == r0}, {r, \[Theta]}, {t, 0, 100}, MaxSteps -> 100]

• It's because the expression inside Sqrt becomes negative between t=16.46 and t=16.47. If this should not happen, there's probably something wrong with the equation itself. Sep 24 '19 at 6:22
• This might be the case. I'll see if I've gone wrong there and edit my post if that's the case, thank you :) Sep 24 '19 at 8:06

Sign r'[t] changes when reflected from points rm = r /. Solve[L^2/(2 r^2) + V^2 Log[r] == energy-vr^2/2, r]. Therefore, the model should be subject to the sign r'[t]. But it’s better to use a second order equation

(*all quantities are in their natural units,per unit kg*)(*angular \
momentum of star in J (conserved i.e constant)*)L = 3.086*10^25;
(*V is a constant in km/s*)
V = 200*10^3;
(*initial radial velocity (in polar coordinates)*)
vr = 50*10^3;
(*initial tangential velocity (in polar coordinates)*)
vt = 10^5;
r0 = 10*3.086*10^19;
(*seconds in a mega year (units of time I want to plot in*)
megayear = 3.154*10^7*10^6;
(*energy of the star system calculated from initial conditions \
(conserved i.e constant energy),the V^2 Log(r0) term is the spherical \
potential*)
energy = vr^2/2 + L^2/(2*r0^2) + V^2*Log[r0];
rm = r /.
Solve[L^2/(2 r^2) + V^2 Log[r] == energy - vr^2/2,
r];(*solve for the max and min radii of the orbit*)

s = NDSolve[{r''[t] == (-V^2/r[t] + L^2/r[t]^3)*megayear^2, \[Theta][
0] == vt/r0, r[0] == r0, r'[0] == vr}, {r, \[Theta]}, {t, 0,
400}]
Plot[{r[t] /. s, rm}, {t, 0, 400}, PlotRange -> All,
AxesOrigin -> {0, 0}, AxesLabel -> Automatic]


• Wow, that's brilliant. Thank you very much! Sep 24 '19 at 8:49
• @charl1e you're welcome! Sep 24 '19 at 8:51
• Why is the equation modified in this way? Sep 24 '19 at 9:15
• @xzczd We have r'[t]==Sqrt[...], then r'[t]^2==(...), and therefore  r'[t]r''[t]==D[(...),r]r'[t]. Sep 24 '19 at 9:21
• I see. The ODE is deduced based on energy conservation, right? Sep 24 '19 at 9:26