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]
Improve this question
$\endgroup$
2
  • 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$ Commented Sep 24, 2019 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$ Commented Sep 24, 2019 at 8:06

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
6
  • $\begingroup$ Wow, that's brilliant. Thank you very much! $\endgroup$ Commented Sep 24, 2019 at 8:49
  • $\begingroup$ @charl1e you're welcome! $\endgroup$ Commented Sep 24, 2019 at 8:51
  • $\begingroup$ Why is the equation modified in this way? $\endgroup$ Commented Sep 24, 2019 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$ Commented Sep 24, 2019 at 9:21
  • $\begingroup$ I see. The ODE is deduced based on energy conservation, right? $\endgroup$ Commented Sep 24, 2019 at 9:26

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.