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]
Sqrt
becomes negative betweent=16.46
andt=16.47
. If this should not happen, there's probably something wrong with the equation itself. $\endgroup$