2
$\begingroup$

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;
(* initial radius *)
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]
$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – xzczd Sep 24 at 6:22
  • $\begingroup$ 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 :) $\endgroup$ – charl1e Sep 24 at 8:06
3
$\begingroup$

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;
(*initial radius*)
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]

Figure 1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.