• Parameterizing elliptical motion
• Relative motions
• Finding a path determined by gravity

Planetary Motion

Worksheet by Mike May, S.J.- maymk@slu.edu

> restart;

Section 16.5 deals with planetary motion. That is actually a good example to study since one could make a good case for the claim that Newton invented calculus to solve problems of planetary motion.

This is a demonstration worksheet without exercises.

Parameterizing elliptical motion

The first approach to planetary motion in the class is to look at trying to parameterize the elliptical paths of the planets. We want to start with parameterizing circular motion and modifying it.

If we want to move around a circle of randius r and period omega the easiest parameterization is

x(t) = a*cos(2*Pi*t/omega), y(t) = a*sin(2*Pi*t/omega).

For the earth, r=93 in millions of miles, and omega = 1 in years.

> a := 93: omeg := 1:
x := t -> a*cos(2*Pi*t/omeg): y := t ->a*sin(2*Pi*t/omeg):
plot([x(t), y(t), t=0..1] );

Next we want to take into account that the orbit of the earth is flattened into an ellipse. The eccentricity of the earth's orbit is .02.

That means that the sum is shifted c = e*a = .02*93 million miles off center.

It also means that the short axis is b = sqrt(a^2-c^2) = a*sqrt(1-e^2) million miles.

> e := .02:
c := a*e;
b := a*sqrt(1-e^2);

This leads to a new parameterizaton of the earth's motion.

> a := 93: omeg := 1:
e := .02: c := a*e: b := a*sqrt(1-e^2):
x := t -> a*cos(2*Pi*t/omeg) + c:
y := t -> b*sin(2*Pi*t/omeg):
plot([x(t), y(t), t=0..1] );

This still looks like a circle to the unaided eye.

If we look at the orbit of Mercury, a is 36 million miles and e is .21.
We also need to compute the period, which will be (36/93)^(3/2). (The period is the size of the orbit, measured in AUs (The size of the earth's orbit) raised to the 3/2 power.)

> a := 36: omeg := evalf((a/93.0)^(3/2));
e := .21: c := a*e: b := a*sqrt(1-e^2):
x := t -> a*cos(2*Pi*t/omeg) + c:
y := t -> b*sin(2*Pi*t/omeg):
plot([x(t), y(t), t=0..omeg], scaling=constrained);

For comparison we plot a circle at the same center.

> xc := t -> a*cos(2*Pi*t/omeg) + c:
yc := t -> a*sin(2*Pi*t/omeg):
plot({[xc(t), yc(t), t=0..omeg], [x(t), y(t), t=0..omeg]}, scaling=constrained);

It is still very close to circular in orbit, even if noticably off center.

For Halley's commet the numbers are that a=1680 and e=.97.

> a := 1680: omeg := evalf((a/93.0)^(3/2));
e := .97: c := a*e: b := a*sqrt(1-e^2):
x := t -> a*cos(2*Pi*t/omeg) + c:
y := t -> b*sin(2*Pi*t/omeg):
plot([x(t), y(t), t=0..omeg], scaling=constrained);

This orbit is clearly an ellipse. It is worth noting that our model would have the period of Mercury be 88 days and Haley's comet be 77 years.

Relative motions

The next question of interest is to plot relative motion fo two planets. In parametric form, we simply take the difference of the two parameterizations. Not that we want to take our time period to be at least the longer of the two periods.

Consider Mercury and the earth together.

> am := 36: omegm := evalf((am/93.0)^(3/2)):
em := .21: cm := am*em: bm := am*sqrt(1-em^2):
xm := t -> am*cos(2*Pi*t/omegm) + cm:
ym := t -> bm*sin(2*Pi*t/omegm):
ae := 93: omege := evalf((ae/93.0)^(3/2)):
ee := .02: ce := ae*ee: be := ae*sqrt(1-ee^2):
xe := t -> ae*cos(2*Pi*t/omege) + ce:
ye := t -> be*sin(2*Pi*t/omege):
tmax := max(2*omegm, 2*omege);
plot([xm(t) - xe(t), ym(t) - ye(t), t=0..tmax], scaling=constrained);

Looking from one planet to the other it appears that the planet is stopping and starting with some backing up as well.

Consider what happens if we try the earth and halley's comet.

> am := 1680: omegm := evalf((am/93.0)^(3/2)):
em := .97: cm := am*em: bm := am*sqrt(1-em^2):
xm := t -> am*cos(2*Pi*t/omegm) + cm:
ym := t -> bm*sin(2*Pi*t/omegm):
ae := 93: omege := evalf((ae/93.0)^(3/2)):
ee := .02: ce := ae*ee: be := ae*sqrt(1-ee^2):
xe := t -> ae*cos(2*Pi*t/omege) + ce:
ye := t -> be*sin(2*Pi*t/omege):
tmax := max(2*omegm, 2*omege);
plot([xm(t) - xe(t), ym(t) - ye(t), t=0..tmax], scaling=constrained);

The motion is obviously wrong. We have neglected that fact that planets move at different speeds on their orbits. They go faster when they are near the sun.

>

Finding a path determined by gravity

If we are to look at planetary motion from calculus, what we know is that the only force that acts on the planets is gravity and that it is in the direction of the sun with a force inversely proportional to the distance.

In terms of differential equations this becomes:

x'' = -x/r^3, y'' = -y/r^3.

(Note that x^2 + y^2 = r^2 and that sqrt((x'')^2+(y'')^2)=1/r^2.)

We have only worked with vector fields that are first order differential equations and this system is second order. We handle that problem by defining u=x' and v=y' to get a system of 4 first order equations.

> x:='x': y:= 'y': u :='u': v:= 'v': t:= 't':

> ODEa := diff(x(t),t)=u(t):
ODEb := diff(y(t),t)=v(t):
ODEc := diff(u(t),t)=-x(t)/(x(t)^2+y(t)^2)^(3/2):
ODEd := diff(v(t),t)=-y(t)/(x(t)^2+y(t)^2)^(3/2):
eqns := [ODEa, ODEb, ODEc, ODEd];

We then specify initial conditions and use the command dsolve to produce a solution function.

> ICx := x(0) = 1:
ICy := y(0) = 0:
ICu := u(0) = 0:
ICv := v(0) = 1:
solcurve := dsolve({ODEa, ODEb, ODEc, ODEd, ICx, ICy, ICu, ICv},
[x(t), y(t), u(t), v(t)], numeric);

Then we can look at where the planet is at times Pi/2, Pi, 3Pi/2, and 2Pi.

> P1 := solcurve(Pi/2);
P2 := solcurve(Pi);
P3 := solcurve(3*Pi/2);
P4 := solcurve(2*Pi);

We are really interested in the values of x and y at those 4 times.

> P1a := [op(2,P1[2]),op(2,P1[3])];
P2a := [op(2,P2[2]),op(2,P2[3])];
P3a := [op(2,P3[2]),op(2,P3[3])];
P4a := [op(2,P4[2]),op(2,P4[3])];

This lets us see that with 6 decimal points of accuracy, the path is a circle with period 2Pi. We can graph the results to see this more clearly.

> with(plots):
pointplot([seq([op(2,solcurve(t*.05)[2]),op(2,solcurve(t*.05)[3])],
t=0..125)], connect=true);

Warning, the name changecoords has been redefined

We can also see what happens if we change the eccentricity.

> e := 0.8; a := 1; tmax := 2*Pi*a*sqrt(a):
ICx := x(0) = a*(1-e);
ICy := y(0) = 0:
ICu := u(0) = 0:
ICv := v(0) = evalf(sqrt((1+e)/(a*(1-e))));
solcurve := dsolve({ODEa, ODEb, ODEc, ODEd, ICx, ICy, ICu, ICv},
[x(t), y(t), u(t), v(t)], numeric);
pointplot([seq([op(2,solcurve(t*tmax/100)[2]),
op(2,solcurve(t*tmax/100)[3])],
t=0..100)], connect=true, scaling=constrained);

We see that we get an ellipse for the orbit.

We are now ready to return to the question of relative orbits.

Consider first the problem of the earth and Mercury. We will measure in astronomical units.

> ee := 0.02; ae := 1; tmaxe := 2*Pi*ae*sqrt(ae);
ICxe := x(0) = ae*(1-ee);
ICye := y(0) = 0:
ICue := u(0) = 0:
ICve := v(0) = evalf(sqrt((1+ee)/(ae*(1-ee))));
solcurvee := dsolve({ODEa, ODEb, ODEc, ODEd, ICxe, ICye, ICue, ICve},
[x(t), y(t), u(t), v(t)], numeric);
em := 0.21; am := 36/93; tmaxm := 2*Pi*am*sqrt(am);
ICxm := x(0) = am*(1-em);
ICym := y(0) = 0:
ICum := u(0) = 0:
ICvm := v(0) = evalf(sqrt((1+em)/(am*(1-em))));
solcurvem := dsolve({ODEa, ODEb, ODEc, ODEd, ICxm, ICym, ICum, ICvm},
[x(t), y(t), u(t), v(t)], numeric);

> tmax := 2*max(tmaxe, tmaxm);
pointplot([seq(
[op(2,solcurvee(t*tmax/200)[2])-op(2,solcurvem(t*tmax/200)[2]),
op(2,solcurvee(t*tmax/200)[3])-op(2,solcurvem(t*tmax/200)[3])],
t=0..400)], connect=true, scaling=constrained);

Notice the slight change in the relative path of the two planets. Some of the inner loops are clearly smaller than others.

Now we try the same trick with Halley's comet.

> ee := 0.02; ae := 1; tmaxe := 2*Pi*ae*sqrt(ae);
ICxe := x(0) = ae*(1-ee);
ICye := y(0) = 0:
ICue := u(0) = 0:
ICve := v(0) = evalf(sqrt((1+ee)/(ae*(1-ee))));
solcurvee := dsolve({ODEa, ODEb, ODEc, ODEd, ICxe, ICye, ICue, ICve},
[x(t), y(t), u(t), v(t)], numeric);
em := 0.97; am := 1680/93; tmaxm := evalf(2*Pi*am*sqrt(am));
ICxm := x(0) = am*(1-em);
ICym := y(0) = 0:
ICum := u(0) = 0:
ICvm := v(0) = evalf(sqrt((1+em)/(am*(1-em))));
solcurvem := dsolve({ODEa, ODEb, ODEc, ODEd, ICxm, ICym, ICum, ICvm},
[x(t), y(t), u(t), v(t)], numeric);
tmax := 2*max(tmaxe, tmaxm);
pointplot([seq(
[op(2,solcurvee(t*tmax/1000)[2])-op(2,solcurvem(t*tmax/1000)[2]),
op(2,solcurvee(t*tmax/1000)[3])-op(2,solcurvem(t*tmax/1000)[3])],
t=-100..100)], connect=true, scaling=constrained);

Some things to notice from the computations.

At its closest point, Halley's comet is inside the earth's orbit and moving almost twice as fast.

If a year is 2Pi, then Halley's commet orbits once every (482/(2Pi) = 76.7 years.

The part of the orbit close in is worth looking at in more detail.

> pointplot([seq(
[op(2,solcurvee(t*tmax/5000)[2])-op(2,solcurvem(t*tmax/5000)[2]),
op(2,solcurvee(t*tmax/5000)[3])-op(2,solcurvem(t*tmax/5000)[3])],
t=-100..100)], connect=true, scaling=constrained);

The reason for the odd shape of the path is that we have the commet and the Earth closest to the sun at the same time and on the same side. This has the comet passing the earth on the inside.

If we have the earth and Haley's comet on different sides of the sun at closest approach, we get a smoother curve.

> pointplot([seq(
[op(2,solcurvee(t*tmax/5000+Pi)[2])-op(2,solcurvem(t*tmax/5000)[2]),
op(2,solcurvee(t*tmax/5000+Pi)[3])-op(2,solcurvem(t*tmax/5000)[3])],
t=-100..100)], connect=true, scaling=constrained);

>

>