Using WhenEvent to calculate the period of an orbit

Question

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.

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. — bbgodfrey, Oct 05 '18 at 20:12

score 4 · Answer 1 · answered Oct 06 '18 at 09:36

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:

Rescale your quantities to "comparable" sizes (i.e. make a good choice of units) for stability
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

score 2 · Answer 2 · answered Oct 06 '18 at 06:14

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

Using WhenEvent to calculate the period of an orbit

2 Answers2