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.04935751leads 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.
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'.
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.
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.
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!
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.
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 - 4000where 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 - 4000If 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.
Look at the function make_movie.m and finish it up. Run a few movies to be sure things are working right.
~cs371/submit lab10to send them to me. These programs will be due by 5:00pm on Monday, April 29.