1
$\begingroup$

I'm having trouble trying to have Mathematica tell me when the period of the orbit of my system is. Namely, I'm trying to use Mathematica to find the orbital period of the Earth around the Sun. Here's my code:

eq = {r''[t] == b/(r[t].r[t]) r[t]/Norm[r[t]]};

inits = {r[0] == {Ro, 0}, r'[0] == {0, Vo}};

event = WhenEvent[r[t] == {Ro, -0.001}, Period = t];

SEperiod = NDSolve[{eq,inits,event},r,{t,0,9.4*10^7}]

It seemed to run fine, but Period is never given a value from t. I tried it again by putting it in NDSolve explicitly:

SEperiod = 
 NDSolve[
  {
   r''[t] == b/(r[t].r[t]) r[t]/Norm[r[t]], 
   r[0] == {Ro, 0}, 
   r'[0] == {0, Vo}, 
   WhenEvent[r[t] == {Ro, -0.001}, Period = t]
   }, 
  r,
  {t, 0, 9.4*10^7}
  ]

Which came with no error messages; but once again, Period wasn't saved. I chose the condition r[t]=={Ro,-0.001} as the triggering for the event because I expect it to be at the position approximately right before the Earth completes a full orbit. However changing around the values still sets no value for Period. How do I get it to be assigned a value like I want?

Edit: Just to make things clear, Ro 1 AU in meters, b is the product of the Gravitational Constant, the mass of the sun, and Vo is the initial velocity of the Earth.

$\endgroup$
3
  • 2
    $\begingroup$ Please provide values for b, Ro and Vo. $\endgroup$
    – Carl Woll
    Oct 5, 2018 at 19:15
  • 1
    $\begingroup$ I would guess that the problem is associated with r[t] == {Ro, -0.001}. Try using something like Norm[r[t]] < 10^-3. Some experimentation may be necessary to find the right number for comparison. $\endgroup$
    – bbgodfrey
    Oct 5, 2018 at 20:12
  • $\begingroup$ A related question. $\endgroup$ Oct 6, 2018 at 9:50

2 Answers 2

4
$\begingroup$

I've been wanting to post an updated version of the technique I used here, so this is a good opportunity. You are right about wanting to use WhenEvent[], but in addition:

  1. Rescale your quantities to "comparable" sizes (i.e. make a good choice of units) for stability
  2. Use a symplectic method to integrate an ODE with a Hamiltonian

With these considerations:

v0 = QuantityMagnitude[UnitConvert[PlanetData["Earth", "AverageOrbitVelocity"],
                       "AstronomicalUnit"/"Years"]]
   6.2783251843074416842`3.

Gm = QuantityMagnitude[UnitConvert[Quantity[1, "GravitationalConstant"]
                                   StarData["Sun", "Mass"], 
                       "AstronomicalUnit"^3/"Years"^2]]
   39.4221026742864858787`4.031532467441943

{xf, yf} = NDSolveValue[Join[
             Thread[{x''[t], y''[t]} == -Gm {x[t], y[t]}/Norm[{x[t], y[t]}]^3], 
             Thread[Through[{x, y}[0]] == {1, 0}], Thread[Through[{x', y'}[0]] == {0, v0}],
             {WhenEvent[{x'[t], y'[t]}.({x[t], y[t]} - {1, 0}) > 0, "StopIntegration"]}],
             {x, y}, {t, 0, ∞}, Method -> {"SymplecticPartitionedRungeKutta",
                                           "PositionVariables" -> {x[t], y[t]}}];

where we find that

Through[{xf, yf}["Domain"]]
   {{{0., 1.00053}}, {{0., 1.00053}}}

Plot a complete period:

ParametricPlot[{xf[t], yf[t]}, {t, 0, xf["Domain"][[1, 2]]}]

Earth orbit for one revolution

$\endgroup$
2
$\begingroup$

For an approximate solution, you can use this code.

r = {x[t], y[t]};
r'' = {x''[t], y''[t]};
eq = r'' + b/(r.r) r/Norm[r] == {0, 0};

inits = {{x[0], y[0]} - {Ro, 0} == {0, 
     0}, {x'[0], y'[0]} - {0, Vo} == {0, 0}};

event = WhenEvent[{y[t] > 0}, Print["t = ", t, " r[t] = ", r]];



G = 6.67408*10^(-11);
mSun = 1.988435`7.*^30;

b = G*mSun;
Ro = 149597870700;
Vo = 29785.1;
SEperiod = 
 NDSolve[{eq, inits, event}, {x, y}, {t, 0, 9.4*10^7}] // Quiet

t = 3.15597*10^7 r[t] = {1.49598*10^11,-0.0000391006}

t = 6.31194*10^7 r[t] = {1.49598*10^11,0.0000686646}

This can be compared to the orbital period of the Earth.

QuantityMagnitude[
 Entity["Planet", "Earth"]["OrbitPeriod"], "SI"]

Out[]= 3.1558149*10^7
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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