Mathematics 172 – Numerical Analysis 2 Lab/Problem Set 2 – ODE Methods; Maple Programming February 10, 1999

Background

Over the last few days, we have introduced the Modified Euler, Midpoint, and 4th order Runge-Kutta methods for generating approximate solutions of ODE initial value problems. In today’s lab and problem set, we will implement the Modified Euler and Runge-Kutta Methods and compare them on some examples.

Implementing the Modified Euler Method

Recall that Modified Euler works as follows. Choosing a step size, h, we generate a sequence of approximations w_{i} ≐ y(t_{i}) to values of the exact solution of y′ = f(t, y); y(t_{0}) = α by the rule:

(ME)  w_{i} = w_{i−1} + (\frac{f(t_{i−1}, w_{i−1}) + f(t_{i−1}, w_{i−1} + f(t_{i−1}, w_{i−1})h)}{2}) · h

We can also break the computation into two steps by realizing that the w_{i−1}+f(t_{i−1}, w_{i−1})h is just one ordinary Euler’s method step. So (ME) is equivalent to the sequence of two operations:

(ME′)

w_{i}^{Euler} = w_{i−1} + f(t_{i−1}, w_{i−1})h
w_{i} = w_{i−1} + (\frac{f(t_{i−1}, w_{i−1}) + f(t_{i}, w_{i}^{Euler})}{2}) · h

To implement a Maple procedure for carrying out this method, we note that we need to supply the following inputs:

1) the function f(t, y) in the differential equation,

2) the initial information t_{0}, α, and

3) information that will determine the step size to be used. (One possible way to do this is to give a final t-value t = t_{n}, and the number of steps n. Then h = \frac{t_{n}−t_{0}}{n}.)

The “work” done by the procedure will just be to carry out the steps in (ME’) for each i = 1, …, n and report the results. We can output the list of pairs (t_{i}, w_{i}) as the results of the calculation – that contains all the information we are likely to need about what happened!

Here is a rudimentary implementation of the Modified Euler method, taking these choices into account. The list is built up one point [t_{i}, w_{i}] at a time in the statement vals := [op(vals),[t,w]];, which appends the new point [t,w] to the end of the list vals.

ModEuler := proc(f,t0,alpha,tf,n)
    local i, h, vals, wE, t, w, lastt, lastw;
    h := evalf((tf - t0)/n);
    vals := [[t0,alpha]];
    lastt := t0;
    lastw := alpha;
    for i to n do
       t := lastt + h;
       wE := lastw + f(lastt,lastw)*h;
       w := lastw + (f(lastt,lastw)+f(t,wE))*h/2;
       vals := [op(vals),[t,w]];
       lastw := w;
       lastt := t;
    od;
    RETURN(vals)
    end:

Lab Questions

A) The Maple code for the ModEuler procedure is available on the course homepage:

           http://math.holycross.edu/ little/Num299.html

Copy and paste from the Netscape browser window into a text file, save the procedure file, then read it into your Maple session. Use the ModEuler procedure to find approximations to the solution of

\{ \begin{matrix}y′ = sin(t) cos(t + y^{2}) \\ y(2) = 1 \end{matrix}

on the interval [2, 5]. Use n = 10, 20, 100 steps and compare the approximations to y(5). (Suggestion: one nice way to do this is to plot the lists of points returned by the procedures. Recall that you can do this in Maple by combining the output of several plot commands with the display command from the plots package. Different colors are a good way to distinguish the different lists of approximate solution points!)

B) What would you need to change in ModEuler to get a procedure MidPoint that computes approximations to the solution of y′ = f(t, y) by the midpoint method? Recall that the midpoint method step consists of the following calculations:

(MP) dw^{Euler} = f(t_{i−1}, w_{i−1})h w_{i} = w_{i−1} + f(t_{i−1} + h/2, w_{i−1} + dw^{Euler}/2)h

Note that the dw^{Euler} here is just the difference w_{i} − w_{i−1} using the basic Euler method. (Hint: Only one or two lines in ModEuler actually need to be modified to get Midpoint!) Code this procedure also and repeat the calculations from A using this new method. (You can put the code for Midpoint in the same file as ModEuler to start a small “library” of numerical ODE methods. You will add to the library below!)

C) The Runge-Kutta 4th order method we discussed in class on Tuesday can be described like this: In each step we compute:

\begin{aligned}k_{1} & = f(t_{i−1}, w_{i−1})h \\ k_{2} & = f (t_{i−1} + \frac{h}{2}, w_{i−1} + \frac{k_{1}}{2} ) h \\ k_{3} & = f (t_{i−1} + \frac{h}{2}, w_{i−1} + \frac{k_{2}}{2} ) h \\ k_{4} & = f(t_{i−t} + h, w_{i−1} + k_{3})h \end{aligned}

Then,

w_{i} = w_{i−1} + \frac{1}{6}(k_{1} + 2k_{2} + 2k_{3} + k_{4}).

Code a Runge-Kutta procedure and again test it on the initial value problem from questions A and B.

D) Problems 4d, 8d, 11d from section 5.4 of Burden and Faires. (Plot the graph of the exact solution together with the approximate values generated by each method.)

E) Problem 15 from section 5.4 of Burden and Faires.

F) (non-lab problem) Problem 16 from section 5.4 of Burden and Faires. Once again, recall that B+F’s “local truncation” error is what we would call \frac{1}{h} times the one-step error (i.e. in the notation of the problem the error in w_{i+1}, assuming w_{i} is exact).

Assignment

Individual lab reports, due Thursday, February 18.