Execises (3)

Assigned: February 7, 2013
Due: February 14, 2013

Setup your Solution Notebook (in-class)

In order to solve these execises, you need need to define the equations of motion and a particular solution. Copy the equations of motion for an orbiting planet and demonstrate that your equations are correct by recomputing and plotting Halley's orbit.

Solution

p = {x[t], y[t], vx[t], vy[t]};

dpdt = {D[x[t],t] == vx[t],
D[y[t],t] == vy[t],
D[vx[t],t] == -4 (π^2) x[t]/r^3,
D[vy[t],t] == -4 (π^2) y[t]/r^3} /. r -> Sqrt[x[t]^2 + y[t]^2]

p0 = {x[0] == x0, y[0] == y0,vx[0] == vx0, vy[0] == vy0}

eqs = dpdt ~ Join ~p0

halley = NDSolve[eqs /.x0 -> 0.59 /. y0 -> 0 /. vy0->11.46/. vx0 -> 0, p, {t, 0, 76},MaxSteps -> 2000] // Flatten

ParametricPlot[Evaluate[{x[t],y[t]}/.halley],{t,0,76},AspectRatio->Automatic]

Conservation of Energy (in-class)

Planetary motion conserves the total energy as well as the angular momentum. The total energy is the sum of the kinetic energy and the gravitational potential energy.
	Energy  = (1/2) (v^2) - 4(π^2)/ r
Show that the total energy of an orbiting body is constant.

Hint

For this problem, you need to define a function that returns the total energy for Halley's comet as a function of time.

Start by defining (and testing!) a function of the form:
	energyHalley[t_] := (...) /. halley

This function uses the solution to the differential equations (labeled by "halley") to substitute for your expressions within the parentheses.

Once your function has been defined, you can plot the error as a function of time by writing...
	Plot[100 (energyHalley[t]-energyHalley[0])/energyHalley[0], {t, 0, 76}, 
  		PlotLabel->"Percent Error of Total Energy"]

Solution

energyHalley[t_] = (vx[t]*vx[t] + vy[t]*vy[t])/2 - 4(π^2)/Sqrt[x[t]*x[t] + y[t]*y[t]] /. halley

Plot[((100 (energyHalley[t]-energyHalley[0]))/(energyHalley[0])),{t,0,76},PlotLabel->"Percent Error of Total Energy"]

The Planets (in-class)

The orbital parameters of the planets of our solar system are

pname = {"Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"};
pradius = {0.39, 0.72, 1.00, 1.52, 5.20, 9.54, 19.19, 30.06, 39.53};
peccentricity = {0.206, 0.007, 0.017, 0.093, 0.048, 0.056, 0.046, 0.010, 0.248};

TableForm[{pname,pradius, peccentricity},TableDirections->{Row,Column},TableHeadings->{{" ", "Radius","Eccentricity"}, Automatic}]

Find the harmonic constant, (T^2)/(a^3), for each of the nine planets.

In order to solve this problem, you need to know the initial conditions for each planet. As Kepler showed, the initial conditions determine the geometry of the orbit and vice-versa.

The function that determines the velocity of an object orbiting the sun at the closest distance to the sun is
	vc[a_, e_] := 2π (Sqrt[((1+e)/(a - a e))])

where "a" is the semimajor radius and "e" is the eccentricity.

Hint

For simplifity, take the planet's mean radius to be equal to the planet's semimajor radius, a. Then, find the planet's orbital period (including the planet's eccentricity). Remember, in our astronomical units,  (T^2)/(a^3) = 1.

You may solve this problem in many ways, but one that I like begins by finding the period to orbit from one side of the sun to the other. This is one-half of the period, T.

Let's use a Module[...] and define a local variable, sol, that is the solution to the equations of motion. Then, use FindRoot[...], to find the time when the orbit crosses the x-axis. For example,
	findT[a_,e_] :=Module[{sol},
    		sol =  NDSolve[...] // Flatten;
    		2t /. FindRoot[Evaluate[vx[t]/.sol],{t,  a^(3/2)/2}]]

With this function, the orbital periods of the planets in years is obtained with...
	pT =  Table[findT[pradius[[i]], peccentricity[[i]]],  {i,  1,  9}]

Further proof of Kepler's Third Law is found with...
	pHarmonic = Table[((pT[[i]])^2)/((pradius[[i]])^3),{i,1,9}]

Solution

First, we define the initial condition from the geometry of the orbit.

vc[a_, e_] := 2π Sqrt[(1+e)/(a - a e)]

Next, we'll take the planet's mean radius to be equal to the planet's semimajor radius, a. NOTICE: we must start the planet at it's closest distance to the sun, equal to a (1 - e). Then, we must find the planet's orbital period (including the planet's eccentricity). We note in our astronomical units,  (T^2)/(a^3) = 1. We use the FindRoot[ ] command to find when the planet crosses the x-axis at one-half of an orbital period.

findT[a_,e_] :=Module[{sol, tGuess},
tGuess = Sqrt[a^3];sol =NDSolve[eqs /.x0 ->a (1 - e) /. y0 -> 0 /.vy0->vc[a,e]/. vx0 -> 0, p, {t, 0, tGuess},MaxSteps -> 2000] // Flatten;2t /. FindRoot[Evaluate[vx[t]/.sol], {t, tGuess/2}]]

The orbital periods of the planets in years...

(pT= Table[findT[pradius[[i]],peccentricity[[i]]], {i, 1, 9}] )// TableForm

Further proof of Kepler's Third Law:

(pHarmonic = Table[((pT[[i]])^2)/((pradius[[i]])^3),{i,1,9}] )// TableForm

TableForm[{pname,pradius, peccentricity, pT, pHarmonic},TableDirections-> {Row,Column},TableHeadings-> {{" ", "Radius","Eccentricity", "Period", "T^2/a^3"}, Automatic}]

Using Euler-Cromer

Write a procedure to integrate the equations of motion for an orbiting object using the Euler-Cromer method. Use your procedure to integrate the orbits for the planet Earth and the comet Halley. Examine your solution for a few values of the time step, dt, and discuss the accuracy of this method (with a fixed and constant time step) to the adaptive time-step method used by NDSolve[…].

Hint

Define your solution on the (x, y) plane. This implies that you have four unknown functions of time: x(t), y(t), vx(t), and vy(t). Your procedure needs to map a list of four numbers, {x_, y_, vx_, vy_}, to four new numbers corresponding to the new orbital position of the of the object. Your procedure will look something like:

	orbitStep[{x_, y_, vx_, vy_}] := Module[{vxNext, vyNext, …},
		…;
		{list of four nunbers}]

Of course, you also need to remember the form of the acceleration due to gravity. You must define the radius between the star and the obit, r = (Sqrt[(x^2) + (y^2)]), and use the equations of motion discussed in class. 

Finally, remember the initial conditions for Earth are (x(0), y(0), vx(0), vy(0)) = (1.0, 0.0, 0.0, 2π), and the initial conditions for Halley are (x(0), y(0), vx(0), vy(0)) = (0.59, 0.0, 0.0, 11.46).

Good values of dt to try are dt = 0.001 and dt = 0.01, but be careful: when dt is small, it takes a very long time (about 1 minute) to compute the 100,000 steps needed to represent Halley's orbit!