# Why this differential equation can only find half of the satellite's trajectory

I calculate the orbit of the satellite according to the gravity, but the following code can only find half of the elliptical orbit and output the warning message:

G = 6.672*(10^-11)(*Gravitational constant \
N·m\.b2/kg\.b2*); M =
5.965*10^24(*The mass property of Earth kg*); m = 10(*The mass \
property of satellite kg*);
vθ =
1.3*7.9*10^3;(*First cosmic speed of near earth satellite 7.9km/s*)

vr = 0;

(*ε=1.2;
p=0.5;*)
L = m*r0*vθ;(*Initial angular momentum of satellite*)

\[DoubleStruckCapitalE] =
1/2 m (vr^2 + vθ^2) -
G*(M*m)/r0;(*Initial total energy of satellite*)
p = L^2/(G*M*m^2);
ε = Sqrt[
1 + (2 \[DoubleStruckCapitalE]*L^2)/(G^2*M^2*m^3)];
sol = NDSolve[{r'[θ] ==
r[θ]^2*
Sqrt[(ε/p)^2 - (1/r[θ] - 1/p)^2],
r[0] == r0}, r[θ], {θ, 0, 4 Pi}] // FullSimplify
PolarPlot[r[θ] /. First[sol], {θ, 0, 2 Pi}]


What should I do to plot a complete satellite trajectory?

The derivation process used is as follows:

According to the analysis of the actual process, it is likely that when the satellite turns, that is, when x changes to $$\pi$$, the radial velocity will switch to 0, which leads to the difficulty of solution.

There is a typical typo here when in equation $$r'=r^2\sqrt {…}$$ only one branch is used, while there are two branches $$r'=\pm r^2\sqrt{…}$$. To take both branches into account, we must square the equation and differentiate by $$\theta$$ , then we obtain

G = 6.672*(10^-11)(*Gravitational constant \
N·m\.b2/kg\.b2*); M =
5.965*10^24(*The mass property of Earth kg*); m = 10(*The mass \
property of satellite kg*);
r0 = 6.371*10^6;(*Earth radius 6371km*)vθ =
1.3*7.9*10^3;(*First cosmic speed of near earth satellite \
7.9km/s*)vr = 0;

(*ε=1.2;
p=0.5;*)
L = m*r0*vθ;(*Initial angular momentum of satellite*)\
\[DoubleStruckCapitalE] =
1/2 m (vr^2 + vθ^2) -
G*(M*m)/r0;(*Initial total energy of satellite*)p = L^2/(G*M*m^2);
ε =
Sqrt[1 + (2 \[DoubleStruckCapitalE]*L^2)/(G^2*M^2*m^3)];
sol = NDSolve[{2 r''[θ] ==
D[r[θ]^4 ((ε/p)^2 - (1/r[θ] -
1/p)^2), r[θ]], r[0] == r0,
r'[0] == r0^2*Sqrt[(ε/p)^2 - (1/r0 - 1/p)^2]},
r[θ], {θ, 0, 4 Pi}] // FullSimplify
PolarPlot[r[θ] /. First[sol], {θ, 0, 2 Pi}]


• @PleaseCorrectGrammarMistakes You're welcome! Feb 26 '20 at 11:24

Your ODE is stiff. There is a limit of 10,000 steps build in NDSolve and it reached that at Pi. You can see that from the output

Notice the upper limit is 3.16 and not 2 Pi which is 6.28319

You can improve this by using Method -> "StiffnessSwitching"

sol = NDSolve[{r'[θ] == r[θ]^2*Sqrt[(ε/p)^2 - (1/r[θ] - 1/p)^2], r[0] == r0}, r,
{θ, 0, 2 Pi},   Method -> "StiffnessSwitching"]


And now it goes to 5.47

  Plot[Abs@Evaluate[r[θ] /. sol], {θ, 0, 2 Pi}]


Noticed also the solution becomes complex near Pi. This tells me may be your ODE is not physically correct. You might want to double check the physics.

• Thank you very much, but I can't show the other half of the ellipse after you use polar coordinates(sol = NDSolve[{r'[θ] == r[θ]^2* Sqrt[(ε/p)^2 - (1/r[θ] - 1/p)^2], r[0] == r0}, r, {θ, 0, 2 Pi}, Method -> "StiffnessSwitching"]; PolarPlot[Evaluate[r[θ] /. sol], {θ, Pi, 2 Pi}]). Feb 25 '20 at 6:44
• @PleaseCorrectGrammarMistakes that is correct. It does not show it because near Pi the solution becomes complex as I mentioned above. So you need to do something like PolarPlot[ Re@First[Evaluate[r[\[Theta]] /. sol]], {\[Theta], Pi, 2 Pi}] or PolarPlot[ Abs@First[Evaluate[r[\[Theta]] /. sol]], {\[Theta], Pi, 2 Pi}] But I think the ODE itself is not valid or the physics/numbers used are not right. Because if the physics is correct, one should not get complex solution after Pi Feb 25 '20 at 6:49