\magnification=\magstep1
\overfullrule=0pt
\centerline{Mathematics 172 -- Numerical Analysis 2}
\centerline{Lab/Problem Set 2 -- ODE Methods; Maple Programming}
\centerline{February 10, 1999}
\bigskip
\noindent
{\bf Background}
\bigskip
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.
\bigskip
\noindent
{\bf Implementing the Modified Euler Method}
\bigskip
Recall that Modified Euler works as follows.  Choosing a
step size, $h$, we generate a sequence
of approximations $w_i \doteq y(t_i)$ to values of the exact solution
of $y' = f(t,y); y(t_0) =\alpha$ by the rule:
$$w_i = w_{i-1} 
+ \left( {f(t_{i-1},w_{i-1}) + f(t_{i-1}, w_{i-1} + f(t_{i-1},w_{i-1})h)\over 2}
\right)\cdot h\leqno(ME)$$
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:
$$\eqalign{w_i^{\rm Euler} &= w_{i-1} + f(t_{i-1},w_{i-1})h\cr
w_i &= w_{i-1} 
+ \left({f(t_{i-1},w_{i-1}) + f(t_i,w_i^{\rm Euler})\over 2}
\right)\cdot h}\leqno(ME')$$
To implement a Maple procedure for carrying out this method, we note 
that we need to supply the following inputs:
\bigskip
\item{1)}  the function $f(t,y)$ in the differential equation,
\item{2)}  the initial information $t_0, \alpha$, and
\item{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 = {t_n - t_0\over n}$.)
\bigskip
\noindent
The ``work'' done by the procedure will just be to carry out the 
steps in (ME') for each $i = 1,\ldots, 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
{\tt vals := [op(vals),[t,w]];}, which appends the new point 
{\tt [t,w]} to the end of the list {\tt vals}.

\vbox{\settabs 20 \columns
\+&\tt ModEuler := proc(f,t0,alpha,tf,n)\cr
\+\cr
\+&&\tt local i, h, vals, wE, t, w, lastt, lastw;\cr
\+\cr
\+&&\tt h := evalf((tf - t0)/n);\cr
\+&&\tt vals := [[t0,alpha]];\cr
\+&&\tt lastt := t0;\cr
\+&&\tt lastw := alpha;\cr
\+&&\tt for i to n do\cr
\+&&&\tt t := lastt + h;\cr
\+&&&\tt wE := lastw + f(lastt,lastw)*h;\cr
\+&&&\tt w := lastw + (f(lastt,lastw)+f(t,wE))*h/2;\cr
\+&&&\tt vals := [op(vals),[t,w]];\cr
\+&&&\tt lastw := w;\cr
\+&&&\tt lastt := t;\cr
\+&&\tt od;\cr
\+&&\tt RETURN(vals)\cr
\+&&\tt end:\cr}
\bigskip
\noindent
{\bf Lab Questions}
\bigskip
\item{A)}
The Maple code for the {\tt ModEuler} procedure
is available on the course homepage:  
\medskip
\centerline{\tt http://math.holycross.edu/~little/Num299.html}
\medskip
\item{}  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 {\tt ModEuler} procedure to find approximations
to the solution of 
$$\cases{y' = \sin(t)\cos(t+y^2)& \cr
         y(2) = 1& \cr}$$
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 {\tt plot} commands
with the {\tt display} command from the {\tt plots} package.  Different
colors are a good way to distinguish the different lists of approximate
solution points!)
\bigskip
\item{B)}  What would you need to change in {\tt ModEuler} to get
a procedure {\tt MidPoint}
that computes approximations to the solution of $y' = f(t,y)$
by the {\it midpoint method\/}?   Recall that the midpoint method
step consists of the following calculations:
$$\eqalign{dw^{\rm Euler} &= f(t_{i-1},w_{i-1})h\cr
w_i &= w_{i-1} +  f(t_{i-1} + h/2,w_{i-1} + dw^{\rm Euler}/2) h}\leqno(MP)$$
Note that the $dw^{\rm Euler}$ here is just the difference $w_i - w_{i-1}$
using the basic Euler method.  (Hint:  Only one or two lines 
in {\tt ModEuler} actually
need to be modified to get {\tt Midpoint}!)  Code this procedure also and repeat the calculations
from A using this new method.  (You can put the code for {\tt Midpoint}
in the same file as {\tt ModEuler} to start a small ``library'' of
numerical ODE methods.  You will add to the library below!)
\bigskip
\item{C)}  The Runge-Kutta 4th order method we discussed in class on
Tuesday can be described like this:  In each step we compute:
$$\eqalign{k_1 &= f(t_{i-1},w_{i-1})h\cr
           k_2 &= f\left(t_{i-1} + {h\over 2},w_{i-1} + {k_1\over 2}\right)h\cr
           k_3 &= f\left(t_{i-1} + {h\over 2},w_{i-1} + {k_2\over 2}\right)h\cr
           k_4 &= f(t_{i-t} + h, w_{i-1} + k_3)h\cr}$$
Then, 
$$w_i = w_{i-1} + {1\over 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.
\bigskip
\item{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.)
\bigskip
\item{E)}  Problem 15 from section 5.4 of Burden and Faires.
\bigskip
\item{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 ${1\over 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).
\bigskip
\noindent
{\bf Assignment}
\bigskip
Individual lab reports, due Thursday, February 18.
\bye