O dear Ophelia, I am ill at these numbers. I have not
   	art to reckon my groans.
   	Hamlet 2.2 (William Shakespeare)

Examples of Numerical Methods for 
First Order Initial Value Problems

Steve Dunbar

Department of Mathematics and Statistics
University of Nebraska-Lincoln
Lincoln, NE 68588-0323
sdunbar@mathlab01.unl.edu

Math 221H, Fall 1993

Covers Nagle & Saff, Section 3.5

Before Starting

We  need to access the programs for the numerical methods discussed in class since these programs are not built into Mathematica as pre-loaded programs.  The following commands will do this necessary task.

<</LocalLibrary/MathLabNotes/M221HNotes/Euler.m

<</LocalLibrary/MathLabNotes/M221HNotes/ImprovEuler.m

<</LocalLibrary/MathLabNotes/M221HNotes/RKsimple.m

Explanation of Notation for the Interested!

The Mathematica notation <<.....    tells Mathematica to load a text file containing a program of further Mathemtica commands.

Mathematica searches a set  of standard places to find programs when asked to load a program.  The  set of locations is kept in the Path environment variable (just as in the UNIX operating system for those who know the lingo) which is

$Path

Since some of the programs are simple minded example programs that I wrote, I keep them in the MathLab Notes/NotesM221H directory , rather than the system locations.  Since that directory is not in the standard set of locations for Mathematica, I need to specify the full  location to retrieve them.

Analytic Solution for a simple problem

With these I can solve simple first order linear differential equations without trouble so that I have the exact solution available for comparison.  For instance solve,

y' = 4 x - 2y,     y(0) = 2

solnrule = DSolve[ {y'[x] == 4*x - 2*y[x], y[0] == 2},
 y[x], x]

yexact = First[ y[x] /. solnrule]

Numerical Solution and Comparison with Analytic Solution

N[ yexact /. x -> 5]

?Euler

eulersoln =
 Euler[{1, 4 x - 2 y}, {x,y}, {0,2}, {5, 0.1}]

improveulersoln = 
ImprovEuler[{1, 4 x - 2 y}, {x,y}, {0,2}, {5,0.1}]

rksimple = 
RKsimple[{1, 4 x - 2 y}, {x,y}, {0,2}, {5, 0.1}]

Look at the number of points in the Euler solution

Dimensions[eulersoln]

exactplot = Plot[yexact, {x,0,2},
	DisplayFunction -> Identity]

eulerplot = ListPlot[eulersoln, 
	PlotJoined -> True,
	PlotStyle -> {GrayLevel[0.1]},
	DisplayFunction -> Identity]

improveulerplot = ListPlot[improveulersoln,
	PlotJoined -> True,
	PlotStyle -> {GrayLevel[0.3]},
	DisplayFunction -> Identity]

rksimpleplot = ListPlot[rksimple,
	PlotJoined -> True,
	PlotStyle -> {GrayLevel[0.5]},
	DisplayFunction -> Identity]

Now display all of these plots similtaneously with the more accurate succcessively drawn darker.

Show[exactplot, eulerplot, improveulerplot, rksimpleplot,
	DisplayFunction -> $DisplayFunction,
	PlotRange -> {{0,2}, {1,2}}]

Here's a zoom to a region of the plot near the minimum.

Show[exactplot, eulerplot, improveulerplot, rksimpleplot,
	DisplayFunction -> $DisplayFunction,
	PlotRange -> {{0.4,0.6}, {0.8,1.2}}]

An example with no analytic solution

Now try the equation from the class notes which does not have a solution in terms of elementary functions

DSolve[ y'[x] == 1 - x - y[x] + x*y[x]^2, y[x] , x]

ImprovEuler[{1, 1 - x - y + x*y^2}, {x,y}, {0,0}, {5, 0.1}];

ListPlot[%, PlotJoined -> True]

Exercises

#4D. Experiment 2.6, page 61, Borelli, Coleman, Boyce.  
(a)  Find the explicit, exact solution to the (separable) IVP : y' = y sin(3x), y(0) = 1; on the interval 0 ≤ x ≤ 4.  Evaluate at the right endpoint of the interval x = 4.

(b)  Solve the  IVP approximately on the given interval user Euler's method with 8 steps, then 16 steps then 32 steps.  

(c)  Calculate the error in the approximate Euler solution from the exact value at x = 4.
By what factor is the error decreasing each time?  Is this conssistent with the prediction of the text?

#17, page 108, Nagle & Saff:  Use Euler's Method with a step size h = 0.1 to approximate the solution to the initial value problem
	y' = -20 y, y(0) = 1
on the interval 0 <= x < = 1.  Compare your answers to the true solution y(x) = exp(-20).
What went wrong?  Try drawing the direction field to get an understanding of what went wrong in the Euler method.  Next, try the step size h = 0.025 and also h = 0.2.  What conclusion can you draw concerning the step size in relation to the steepness of the direction field?

Extra Credit:  Example 3, page 58, Borelli, Coleman, Boyce.  Solve and plot the approximate solution to the IVP y' = 1 - x sin(y), y(0) = 1, on the interval 0 ≤ x ≤ 30 using h = 0.15.  Is this equation solvable by methods we have?  What do you suspect goes wrong?