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

Question

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*);
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[{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.

Answer 1

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}]

Answer 2

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.