Differential Equations: Fall 2020 Online Projects
4.3 Numerical Methods for Systems 259 PART TWO
Inordertohelpincreaseyourunderstandingofthetechniques and concepts taught in this course, and gain some knowledge how to write and implement a computer program, you are being assigned three projects.
1. Additionally,writeout1-2page(s)responsepaper onthe following:
The projects come in two parts; I have outlined the parts and their requirements below.
a. What you did to solve each problem/the solution itself; b. What each project taught you;
c. Any critical insights you have into each project (strengths,weaknesses, etc…);
PARTONE
1. Goto4.3ApplicationofNumericalMethodsforSystems(see
d. Any questions you still have(optional). 2. Theformatshouldbeasfollows:
below) Topic: Comets and Spacecraft 2. Project#1:SpacecraftLanding;
a. 1-2 pages (PDF);
b. Double-spaced;
c. 12 pt with font like Times New Roman; d. Default margins;
e. Name, ID#, project numbers, Course.
Project #2:K:epler’s Law of Planetary(Satellite) Motion ;
Project #3:Your Own Comet.
3. Write each D.E. as First-Order System , and solve it by using
3. Uploadeachprojectasone-file(PDF)onthe Blackboard—-> Assignments;
Numerical Methods for Systems from Maple (details about commands to use see Course Material (Blackboard) and Maple Help).
4. Advice:Startsoonaspossible,becausetaketimetobe familiarized with Maple.
4. WorkonthethreeprojectsandprintouttheresultsinPDF format.
4.3 Application
Comets and Spacecraft
5. DEADLINE: Dec 19(Saturday) )
Figure 4.3.13 lists TI-Nspire CX CAS and BASIC versions of the two-dimensional Runge–Kutta program RK2DIM. You should note that it closely parallels the one- dimensional Runge–Kutta program listed in Fig. 2.6.11, with a single line there replaced (where appropriate) with two lines here to calculate a pair of x- and y- values or slopes. Note also that the notation used is essentially that of Eqs. (13) and (14) in this section. The first several lines define the functions and initial data needed for Example 1. (The TI-Nspire program can be downloaded from www.pearsonhighered.com/mathstatsresources/.)
Figure 4.3.14 exhibits an n-dimensional M A T L A B implementation of the Runge–Kutta method. The M A T L A B function f defines the vector of right-hand sides of the differential equations in the system x0 D f.t; x/ to be solved. The rkn function then takes as input the initial t-value t, the column vector x of initial x- values, the final t -value t1, and the desired number n of subintervals. As output it produces the resulting column vector T of t -values and the matrix X whose rows give the corresponding x-values. For instance, with f as indicated in the figure, the MATLAB command
[T,X] = rkn(0, [0;1], 5, 50)
then generates the data shown in the table of Fig. 4.3.1 (which lists only every fifth value of each variable).
You can use Examples 1 through 3 in this section to test your own imple- mentation of the Runge–Kutta method. Then investigate the comet and space- craft problems described next. Additional application material at the web site www.pearsonhighered.com/mathstatsresources/ describes additional nu- merical ODE investigations ranging from batted baseballs to the Apollo orbits shown in Figs. 4.3.9 and 4.3.10.
Your Spacecraft Landing
Your spacecraft is traveling at constant velocity V , approaching a distant earthlike planet with mass M and radius R. When activated, your deceleration system pro- vides a constant thrust T until impact with the surface of the planet. During the period of deceleration, your distance x.t/ from the center of the planet satisfies the
260
Chapter 4 Introduction to Systems of Differential Equations
TI-Nspire CX CAS
Define rk2dim()=Prgm
f(t,x,y):=y
g(t,x,y):=-x n:=50
t:=0.0 x:=0.0 y:=1.0 t1:=5.0 h:=(t1-t)/n For i,1,n
t0:=t
x0:=x
y0:=y
f1:=f(t,x,y)
g1:=g(t,x,y)
t:=t0+h/2
x:=x0+(h*f1)/2
y:=y0+(h*g1)/2
f2:=f(t,x,y)
g2:=g(t,x,y)
x:=x0+(h*f2)/2
y:=y0+(h*g2)/2
f3:=f(t,x,y)
g3:=g(t,x,y)
t:=t0+h
x:=x0+h*f3
y:=y0+h*g3
f4:=f(t,x,y)
g4:=g(t,x,y)
fa:=(f1+2*f2+2*f3+f4)/6
ga:=(g1+2*g2+2*g3+g4)/6
x:=x0+h*fa
BASIC
Program RK2DIM
Comment
y:=y0+h*ga
Disp t,x,y
EndFor
EndPrgm
DEF FN
DEF FN N=50
T=0
X=0
Y=1
T1=5
H = (T1-T)/N FORI=1TON T0=T X0=X Y0=Y
=Y = -X
Program title
Define function f
Define function g
No. of steps
Initial t
Initial x
Initial y
Final t
Step size
Begin loop
Save previous t
Save previous x
Save previous y
First f-slope
First g-slope
Midpoint t
Midpt x-predictor
Midpt y-predictor
Second f-slope
Second g-slope
F(T,X,Y)
G(T,X,Y)
F1 = FNF(T,X,Y)
G1 = FNG(T,X,Y) T=T0+H/2 X=X0+H*F1/2 Y=Y0+H*G1/2
F2 = FNF(T,X,Y)
G2 = FNG(T,X,Y) X=X0+H*F2/2 Y=Y0+H*G2/2
F3 = FNF(T,X,Y)
G3 = FNG(T,X,Y) T=T0+H X=X0+H*F3 Y=Y0+H*G3
F4 = FNF(T,X,Y)
G4 = FNG(T,X,Y)
FA = (F1+2*F2+2*F3+F4)/6 GA = (G1+2*G2+2*G3+G4)/6 X = Y0 + H*FA
Y = Y0 + H*GA
PRINT T,X,Y
NEXT I
Midpt
Midpt
Third
Third
New t
Endpt
Endpt
Fourth f-slope
Fourth g-slope
Average f-slope
Average g-slope
x-corrector
y-corrector
Display results
End loop
x-predictor
y-predictor
f-slope
g-slope
x-predictor
y-predictor
FIGURE 4.3.13.
TI-Nspire CX CAS and BASIC two-dimensional Runge–Kutta programs.
Here, you may rewrite
the MATLAB program to a
Maple program.
FIGURE 4.3.14. MATLAB implementation of the Runge–Kutta method.
differential equation
d2x D T GM ; (1) dt2 x2
where G 6:6726 1011 N(m=kg)2 as in Example 3. Your question is this: At what altitude above the surface should your deceleration system be activated in order to achieve a soft touchdown? For a reasonable problem, you can take
M D 5:97 1024 (kg); R D 6:38 106 (m); V D p 104 (km=h); T D g C q (m=s2)
where g D GM=R2 is the surface gravitational acceleration of the planet. Choose p to be the smallest nonzero digit and q the next-to-smallest nonzero digit in your ID number. Find the “ignition altitude” accurate to the nearest meter and the resulting “descent time” accurate to the nearest tenth of a second.
Kepler’s Law of Planetary (or Satellite) Motion
Consider a satellite in elliptical orbit around a planet of mass M , and suppose that physical units are so chosen that GM D 1 (where G is the gravitational constant). If the planet is located at the origin in the xy-plane, then the equations of motion of the satellite are
d2x D x ; d2y D y : (2) dt2 .x2 C y2/3=2 dt2 .x2 C y2/3=2
h=(t1 T = t; X=x’; fori=
– t)/n;
% step size
% initial t
% initial x-vector % begin loop
% first k-vector
% second k-vector % third k-vector
% fourth k-vector % average k-vector %newt
%newx
% update t-column % update x-matrix % end loop
1:n
k1 = f(t,x);
k2 = f(t+h/2,x+h*k1/2); k3 = f(t+h/2,x+h*k2/2); k4 = f(t+h ,x+h*k3 ); k = (k1+2*k2+2*k3+k4)/6; t=t+h;
x =x+h*k; T = [T;t];
X = [X;x’]; end
4.3 Numerical Methods for Systems 261
function xp = f(t,x) xp = x;
xp(1) = x(2);
xp(2) = -x(1);
function [T,Y] = rkn(t,x,t1,n)
262
Chapter 4 Introduction to Systems of Differential Equations
z
Halley’s comet last reached perihelion (its point of closest approach to the sun at the origin) on February 9, 1986. Its position and velocity components at that time were
p0 D.0:325514;0:459460;0:166229/ and v0 D.9:096111;6:916686;1:305721/
(respectively), with position in AU (astronomical units, in which the unit of distance is the major semiaxis of the earth’s orbit) and time in years. In this system, the three- dimensional equations of motion of the comet are
d2xDx; d2yDy; d2 ́D ́ (4) dt2 r3 dt2 r3 dt2 r3
Let T denote the period of revolution of the satellite. Kepler’s third law says that the square of T is proportional to the cube of the major semiaxis a of its elliptical orbit. In particular, if GM D 1, then
T2 D42a3: (3)
(For details, see Section 11.6 of Edwards and Penney, Calculus: Early Transcen- dentals, 7th ed. (Upper Saddle River, NJ: Prentice Hall, 2008).) If the satellite’s x- andy-componentsofvelocity,x3Dx0Dx10 andx4Dy0Dx20,areintroduced,then the system in (2) translates into a system of four first-order differential equations having the form of those in Eq. (22) of this section.
(a) Solve this 4 4 system numerically with the initial conditions x.0/D1; y.0/D0; x0.0/D0; y0.0/D1
that correspond theoretically to a circular orbit of radius a D 1, so Eq. (3) gives T D 2. Is this what you get?
(b) Now solve the system numerically with the initial conditions
x.0/D1; y.0/D0; x0.0/D0; y0.0/D 1p6
that correspond theoretically to an elliptical orbit with major semiaxis a D 2, so
2
Eq. (3) gives T D 4 Halley’s Comet
p
2. Is this what you get?
5 10 15 20 25 30 y
–5 –10
where
Solve the equations in (4) numerically to verify the appearance of the y ́-projection
of the orbit of Halley’s comet shown in Fig. 4.3.15. Plot the xy- and x ́-projections as well.
Figure 4.3.16 shows the graph of the distance r.t/ of Halley’s comet from the sun. Inspection of this graph indicates that Halley’s comet reaches a maximum distance (at aphelion) of about 35 AU in a bit less than 40 years and returns to perihelion after about three-quarters of a century. The closer look in Fig. 4.3.17 indicates that the period of revolution of Halley’s comet is about 76 years. Use your numerical solution to refine these observations. What is your best estimate of the calendar date of the comet’s next perihelion passage?
p
D42 and rD x2Cy2C ́2:
FIGURE 4.3.15. y ́-projection of the orbit of Halley’s comet.
r
4.3 Numerical Methods for Systems
263
30 20 10
FIGURE 4.3.16.
cusp near t D 75?
15
10
5
74 75 76 77 78t FIGURE 4.3.17. A closer look at Halley’s perihelion
passage after about 76 years.
Your Own Comet
The night before your birthday in 2007 you set up your telescope on a nearby moun- taintop. It was a clear night, and you had a stroke of luck: At 12:30 A.M. you spotted a new comet. After repeating the observation on successive nights, you were able to calculate its solar system coordinates p0 D .x0; y0; ́0/ and its velocity vector v0 D .vx0; vy0; v ́0/ on that first night. Using this information, determine the following:
the comet’s perihelion (point nearest the sun) and aphelion (point farthest from the sun),
the comet’s velocity at perihelion and at aphelion,
the comet’s period of revolution around the sun, and
the comet’s next two dates of perihelion passage.
Using units of length in AU and time in earth years, the equations of motion of your comet are given in (4). For your personal comet, begin with random initial position and velocity vectors with the same order of magnitude as those of Halley’s comet. Repeat the random selection of initial position and velocity vectors, if nec- essary, until you get a plausible eccentric orbit that ranges well outside the earth’s orbit (as most real comets do).
200 t 200-year plot of the distance r.t/ of Halley’s comet from the sun. Is there a
25 50
75 100 125 150 175
r