MA/CS 371 - Lab 10
Solving A Simplified Three Body Problem with a Predictor-Corrector Method


Section 1: Introduction

In this lab you will

Section 2: Topics for Lecture

I will briefly discuss MATLAB's movie function. We will be using it in this lab to show the motion of the objects through time.

Section 3: The System of ODEs

The following differential equations describe the motion of a body in orbit about two much heavier bodies. An example would be a space shuttle in orbit around the earth and moon.

The coordinate system is a little tricky if we try to consider the moon's orbit around the earth at the same time as we consider the shuttle moving around the whole system. Let's simplify things a bit by saying that the center of the earth-moon system is the origin and that the positive x-axis points towards the center of mass of the moon. (This means that the coordinate system is rotating as the moon revolves around the earth.) Then if u represents the relative mass of the moon to the total mass of the earth-moon system, the center of mass of the earth is at (-u,0) and the center of mass of the moon is at (1-u,0). Note that the origin is inside the earth but not at the center of the earth! We can then consider the position of the shuttle as a function of time, (x(t),y(t)). See the figure below.

Approximate values for some constants include u = 0.0121285, distance between earth and moon approximately 238,000 miles, radius of the earth approximately 4000 miles, and radius of the moon approximately 900 miles.

The equations of motion are derived from Newton's Laws and the inverse square law of gravitation. The first derivatives in the equation come from the rotating coordinate system and from a frictional term, which is assumed to be proportional to velocity with proportionality constant f:

	                (1-u)(x+u)     u(x+u-1)
	x'' = 2y' * x - ----------  -  --------  - fx'
			   r1^3          r2^3

			  (1-u)y      uy
	y'' = -2x' + y -  ------  -  ---- - fy'
                           r1^3      r2^3
where r1 = sqrt((x+u)^2 + y^2) and r2 = sqrt((x+u-1)^2+y^2).

Although a great deal is known about these equations, it is not possible to find closed-form solutions in general. One interesting class of problems involves the study of periodic solutions in the absence of friction. It is known that the initial conditions

	x(0) = 1.2,  x'(0) = 0,  y(0) = 0,  y'(0) = -1.04935751
leads to a solution which is periodic with period T = 6.19216933 when f = 0. This means that the space shuttle starts on the far side of the moon with an altitude of about 0.2 times the earth-moon distance and a certain initial velocity. The resulting orbit brings the space shuttle close to the earth, out in a big loop on the opposite side of the earth from the moon, close to the earth again, and finally back to its original position on the far side of the moon.

Your task in this lab will be to write a driver program which, when given a set of initial conditions to this problem, will use the classical fourth-order Runge-Kutta method and the Adams-Bashforth/ Adams-Moulton predictor-corrector method to solve the space shuttle's motion.

Section 4: Vectorizing the System and Data Dictionary

In order to express the system as a vector of state variables, we will adopt the data dictionary technique described in section 9.2 of the text. Define a vector X such that:
Old Variable   New Variable	Differential Equation
     t             X(1)         X(1)' = 1
     x             X(2)         X(2)' = X(4)
     y             X(3)         X(3)' = X(5)
     x'            X(4)         X(4)' = 2*X(5) + X(2) - 
                                 (1-u)*(X(2)+u) * ((X(2)+u)^2+X(3)^2)^(-2/3) -
				 u * (X(2)+u-1) * ((X(2)+u-1)^2+X(3)^2)^(-2/3)
				 - f * X(4)
     y'            X(5)         X(5)' = -2*X(4) + X(3) - 
	                         (1-u) * X(3) * ((X(2)+u)^2+X(3)^2)^(-2/3) -
			         u * X(3) * (X(2)-u+1)^2+X(3)^2)^(-2/3) 
				 - f * X(5)
Now we can simplify discussion of the system by referring to the vector X instead of the individual variables t, x, y, x', and y'.

Section 5: Vectorizing the Classical Fourth-Order Runge-Kutta Method

You will be using the technique in Section 9.3 of the text to generate the solution of this system of ODEs. In order to initialize things, you will need to run the classical fourth-order Runge-Kutta method for three steps to generate the initial conditions. As in Lab 9, I have given you the function func.m which you can call instead of having to hard code the differential equations into your solver. Your first task is to write a function rk4.m whose calling sequence is as follows:
	function Z = rk4(X,h,f)
where Z and X are five-by-one column vectors, and h and f are scalars. This function will apply the classical fourth-order Runge-Kutta method to X as shown on page 365 of the text, with step size h and friction constant f, to generate the answer Z.

Section 6: Implementing the Predictor-Corrector Method

After you get things started off with the Runge-Kutta method, you can switch over to the Adams-Bashforth/Adams-Moulton predictor-corrector method. You should implement this as one function:
	function Z = ab_am(XX,h,f)
where Z is a five-by-one column vector, h and f are scalars, and XX is a five-by-four matrix whose first column is X(t-3h), second column is X(t-2h), third column is X(t-h), and fourth column is X(t). This function should apply the predictor (equation (2) on page 375) to generate a temporary answer and then apply the corrector (equation (3) on page 375) to get the answer Z.

Section 7: Putting the Pieces Together

Once you have implemented vectorized versions of the Runge-Kutta and AB/AM methods, you can put them together into the actual ODE solver. The solver should look like this:
	function XX = ode_solver(Xinit,h,f,tfinal)
where Xinit is the five-by-one column vector of initial conditions, h is the size of one time step, f is the friction constant, tfinal is the time to end the simulation, and XX is a five-by-(tfinal/h plus one) matrix whose first column is Xinit, second column is X(h), third column is X(2h), etc.

IMPORTANT NOTE: Your code will run a lot faster if the first thing you do is to initialize XX to zeros(5,(tfinal/h)+1) so that MATLAB will not have to continually allocate memory while your program is running. Even then, don't panic when the ODE solver takes a couple minutes to run!

Section 8: Checking Your Program

To test your program, try out the model problem given above with a small step size:
	Xinit = [ 0 ; 1.2 ; 0 ; 0 ; -1.0493575 ];
	XX = ode_solver(Xinit,.001,0,6.2);
Note: This will take a few minutes to run.

To plot your results for now, try plotting the earth as a circle near the origin and the moon as an x, and plot the shuttle's path as a line:

hold off; plot(-u,0,'o'); hold on; plot(1-u,0,'x'); plot(XX(2,:),XX(3,:))
Recall that XX(2,:) is the second row of the XX matrix which is the x position of the shuttle through time and XX(3,:) is the y position of the shuttle through time. Your plot should look like this:

If your plot doesn't look anywhere near this, you have a bug somewhere.

Section 9: Did it Hit the Earth?

One thing we haven't considered yet is whether the shuttle crashed to the earth. You can use the output of the ODE solver and a little geometry to determine whether this happened.

Take another look at the picture that explains the coordinate system. We can't just say something like sqrt(x^2+y^2) < 4000, because the earth is not centered at the origin. The expression is a little more complicated:

	dist = sqrt ( (x+u)^2 + y^2 ) * 238000 - 4000
where u is defined as above (u = 0.0121285).

Recall that x is the second row of the output of the ODE solver and y is the third row, you can compute the minimum distance that the shuttle approaches the earth with the one-liner

	min_dist = min(sqrt((XX(2,:)+u).^2+XX(3,:).^2)) * 238000 - 4000
If this number is negative, the shuttle has crashed on the earth and the rest of the ODE solver's solution is invalid.

You will need the above formula to answer a couple of the questions for this lab.

Section 10: Making Movies

Let's make the graphical output a little more interesting by plotting the shuttle's orbit not as a static plot but rather as a motion picture. MATLAB's movie command will do that for us.

Look at the function make_movie.m and finish it up. Run a few movies to be sure things are working right.

Section 11: Turning in Your Work

When your rk4, ab_am, ode_solver, and make_movie functions are ready, run
	~cs371/submit lab10
to send them to me. These programs will be due by 5:00pm on Monday, April 29.

This lab's questions