(b). set up the forward Euler step for the system in (a) with h = O.1.

(c). write a Matlab script that computes 1O Euler steps, to obtain the value ① 10 – ① (1).

9.4    Higher Order ODEs with various methods

①、、- t① = O）       ①(O) = 1）  ①、(O) = 1）

rewrite it as a system of irst order equations.

compute ① (O.1) and ① (O.2) with 2 time steps using h = O.1, using the following methods:

a) Euler,s method,

b) A 2nd order Runge-kutta method,

c) A 4th order Runge-kutta method,

d) The 2nd order Adams-Bashforth-Moulton method.  Note that this is a multi-step method.  For the  2nd initial value ① 1, you can use the solution ① 1  from b).  For this method, please compute ① (O.2) and ① (O.3).

NB! Do not write Matlab codes for these computations. You may use Matlab as a fancy calculator.

9.5 Non-uniqueness and Taylor series method

consider the initial value problem

①、(t) = ①1/3）     ①(O) = O.

verify that the following two functions are both solutions

2 3

① 1(t) = O）      ①2(t) = (3 t) 2 ） for t > O

Now, apply Taylor series method, and try with m = 1 and m = 2. what happens？

9.6

(a).

Ill-conditioned and sensitive ODEs

consider the ODE

①地 = 1O① + 11t - 5t2  - 1,

① (O) = e.

show that the exact solution is：

① (t) = ee10t  + t2 /2 - t.

Now, set e = O.  solve the ODE on the time interval t =  [O, 3] with your RK4 method with h = 2-8.  plot the numerical solution and the exact solution on O 三 t 三 3.  what do you observe？what do you think accounts for the discrepancy between the numerical and exact solution？

(b). Dahlquist and Bj曾ork are credited for this pathological example.   consider  the

ODE

①地 (t) = 1OO(sint - ①), ①(O) = O.

solve it with your RK4 method on O 三 t 三 3, with

h = O.O15,  O.O2O, O.O25, O.O3O.

plot all these solutions. Observe the numerical instability!

9.7    The Lorenz system； A study in chaos

consider the Lorenz system of irst order ODEs

=   T① - g - ①从                                                (2)

从(.)   =    ①g - b从                                                      (3)

with the parameters

σ = 1O,       b = 8/3,       T = 28.

(i) show that ① = O, g = O, 从 = O is an equilibrium solution.

(ii) solve the system numerically in Matlab, using the function ode45, over the time spam t = [O, 8O], and with initial condition ① = O, g = 1, 从 = O.

Make the following plots：

. plot the 3D curve; (can you ind the Matlab command？)

. plot the solution functions ①(t), g(t), 从(t);

.  plot ① against g;

.  plot ① against 从 ;

.  plot g against 从 .

You are welcome to explore with diferent initial conditions and parameters.

what to submit： The Matlab iles (for the function and the script), and all the plots (clearly labelled).

9.8    solving the Airy Equation

In Matlab, solve the Airy diferential equation

①、、= t①

with initial conditions

① (O) = O.355O28O53887817） ①、(O) = -O.2588194O37928O7）

on the interval [-4.5）4.5] using RK4 method.

(Hint：solve the intervals [-4.5）O] and [O）4.5] separately.)

plot the numerical solution ①(t)）①、(t) on the interval [-4.5）4.5].

A point to verify your answer： The value ① (4.5) = O.OOO33O25O34 is correct.

9.9    Matlab solvers:  A Case study in Mechanics

suppose we have two objects orbiting in space, with masses 1 - μ and μ, rotating around each other.   For  example,  think of the earth and the moon, where the moon moves around the earth at distance 1.  (of course, here both the masses and the distance are normalized.)  A third object, which is relatively much smaller and does not afect the motion of the irst two, is also orbiting in space.  Think of this as a space-ship, or a meteorite.

To simplify the analysis, we assume that all trajectories lie in the same plane. we choose the barycenter of mass as the origin of our coordinate system.  Moreover, we adopt a rotating frame of coordinates.  In this moving frame, the earth and the moon are always located at the points

(g1(E)）g2(E)) = (-μ) O)） (g1(M)）g2(M)) = (1 - μ) O)）

and their barycenter is at (O）O).  we denote by  (g1 (t)）g2 (t)) the position of the space ship at time t, w.r.t. this rotating frame of coordinates.

The equations of motion are:

g1(、、)   =   g1 + 2g2(、) - Ⅴ - μ

g2(、、)   =   g2 - 2g1(、) - Ⅴ - μ

dE     =   ((g1  + μ)2 + g2(2))1/2

dM     =   ((g1  - Ⅴ )2 + g2(2))1/2

μ   =   O.O12277471）

Ⅴ   =   1 - μ .

Notice that, in each of above ODEs, theirst two terms account for the coriolis force, due to the fact that the rotating frame is not an inertial frame. The last two terms account for the gravitational pull of the earth and the moon, respectively. The numbers dE）dM measure the distance of the space ship from the earth and from the moon, respectively. consider the initial conditions

g1 (O) = O.994）

g1(、)(O) = O）

g2 (O) = O）

g2(、)(O) = -2.OO15851O6379O825.

For this particular initial data the solution turns out to be periodic.  Indeed, at time t = 17.O65211656 the rotating object goes back to its initial coniguration, and the entire motion is repeated over again.

You need to do the following:

(a) Rewrite the system of 2 second order equations as a system of 4 irst order ODEs.

(b) use Matlab function‘ode45, to solve this system of ODEs over one period.  Draw the orbit.

Hand in your Matlab script and functions, and the outputs and plots.

(c) write a program that would solve the system with forward Euler,s method, using 24OOO uniform time steps.

Hand in your  Matlab script and function of Euler method, and the outputs and plots.

(d) write a program that would solve the system with the classical 4th order Runge- kutta method, but reduce the step number to 6OOO.

Hand in your Matlab script and function of RK4 method, and the outputs and plots.

(e) The total amount of computations in (c) and (d) is about the same.  compare these two solutions. Do you think any of them is acceptable？ Explain why.

(f) Find out how many steps‘ode45, used for computing the solution over one period. can you explain why ,ode45, uses much fewer steps but still gets a much better solution than the other two methods？

(g) Even though the exact solution is periodic, a small error (for example, some per- turbation in initial condition) could throw the object out of its orbit. This means that the numerical error, which in other cases might be invisible, could have big consequences in this case.  In order to see this efect, you could try to solve the system with‘ode45, over 3 periods.  what accuracy (error tolerance) you must require to get an acceptable result？  (In this case, “acceptable” means that the plot of the solution looks reasonable to you.)

Hand in your Matlab script and the outputs/plots. And your com- ments.

Hand in whatever else you think is relevant. Have fun!

9.10    Midpoint Method；Bonus

consider the initial value problem

①地 (t) = -①)       ①(O) = 1.

we see that ①(t) = e-t  is the exact solution.

we wish to compute the numerical solution of using the Midpoint method

地 ①n+1 - ①n-1

①n  =         2h ）

which gives

①n+1  = ①n-1 + 2h①n(地) = ①n-1  - 2h①n）

where h is the mesh size. we set h = O.2 and consider O < t < 5.

(a). write a Matlab code that carry out the computation with the starting values

①0 = 1）       ① 1   = 1 - h.

Note that ① 1  is chosen to be numerical solution of the forward Euler step. plot the numerical solution together with the exact solution.

(b). write a Matlab code that carry out the computation with the starting values

①0 = 1）       ① 1 = -h +“1 + h2）

plot the numerical solution together with the exact solution.