# Coupled differential equation with Sqrt [duplicate]

I'm trying to solve the equation below, but I get:

e = -0.02;
L = 4.0;

U[r_] := -1/r;
r0 = NSolve[U[r] + L^2/(2 r^2) == e, r, Reals]
r0 = Max[r /. r0]

NDSolve[
{r'[t] == Sqrt[2 (e - U[r[t]]) - L^2/r[t]^2], ϕ'[t] == L/r[t]^2,
r[0] == r0, ϕ[0] == 0},
{r, ϕ}, {t, 0, 100}]


During evaluation of In[86]:= NDSolve::mxst: Maximum number of 98930 steps reached at the point t == 0.05406112246188643.

{{r -> InterpolatingFunction[{{0., 0.043953}}, <>],
ϕ -> InterpolatingFunction[{{0., 0.043953}}, <>]}}


How can I solve this?

• Try to find the problem: zed={r[t],φ[t]} /. NDSolve[{r'[t]==Sqrt[2 (-0.02 + 1/r[t]) - 16./r[t]^2], φ'[t]==16./r[t]^2, r[0]==40., φ[0]==0}, {r[t],φ[t]}, {t,0,100}][[1]] followed by Table[zed, {t, 0, .04, .005}] and it appears your r goes to infinity. Are you sure you do not have a sign error? Or another error in the equation? – Bill Apr 25 '17 at 0:03
• This is the equation of motion of a central force with potential given by U(r). I copied this from a Classical Mechanics book. I really don't know what is happening :/ – AVemado Apr 26 '17 at 18:05
• @AVemado Any link to that book> – zhk Apr 27 '17 at 3:15
• It's called Classical Dynamics of Particles and Systems - Marion, Thornton. – AVemado Apr 29 '17 at 17:21
• To be more specific, the problem is caused by Sqrt, the 2 methods suggested in the post linked above are all applicable to your problem i.e. you can either add e.g. WorkingPrecision -> 16 to NDSolve or add a Re@ before Sqrt to resolve the problem. – xzczd May 4 '17 at 16:43

Comment

There are two ways I know which can help to avoid NDSolve::mxst: But the results may not be accurate.

Firstly,

sol1 = NDSolve[{r'[t] == Sqrt[2 (e - U[r[t]]) - L^2/r[t]^2], ϕ'[t] == L/r[t]^2,
r[0] == r0, ϕ[0] == 0}, {r, ϕ}, {t, 0, 100},
Method -> {"FixedStep", Method -> "ExplicitEuler"}, StartingStepSize -> 0.01];


Secondly,

sol2 = NDSolve[{r'[t] == Sqrt[2 (e - U[r[t]]) - L^2/r[t]^2], ϕ'[t] == L/r[t]^2,
r[0] == r0, ϕ[0] == 0}, {r, ϕ}, {t, 0, 100},
Method -> "BDF", PrecisionGoal -> 0, AccuracyGoal -> 0, MaxStepSize -> 1. 10^-2]
`

You can try your luck with other methods.

• Thanks for the answer, but it didn't work. I found another way to calculate this problem. :) – AVemado Apr 26 '17 at 22:26
• @AVemado its works for me and fix the warring. – zhk Apr 27 '17 at 2:50