3
$\begingroup$

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

enter image description here

What should I do to plot a complete satellite trajectory?

The derivation process used is as follows: enter image description here enter image description here enter image description here enter image description here enter image description here

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.

$\endgroup$
11
$\begingroup$

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

Figure 1

$\endgroup$
2
  • $\begingroup$ Your answer is great.Thank you very much for your help. $\endgroup$ Feb 26 '20 at 4:07
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes You're welcome! $\endgroup$ Feb 26 '20 at 11:24
3
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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.

Mathematica graphics

$\endgroup$
2
  • 1
    $\begingroup$ 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}]). $\endgroup$ Feb 25 '20 at 6:44
  • 1
    $\begingroup$ @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 $\endgroup$
    – Nasser
    Feb 25 '20 at 6:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.