Math 3280 - Hampton - Lab 3- Due 6/19

Lab reports should conform to the guidelines explained on the handout.  Collaboration during the lab is strongly encouraged, although the write-up of the lab report should be done individually.

Numerical methods for solving ODEs and an introduction to matrices with Mathematica.

Use the functions below to solve the exercises and in numerical homework assignments.  

In the first three, xStart is the initial x-value, yStart is the initial y-value, and xFinish is the ending x-value. (In the notation used in class, xStart = (Subscript[x, 0]), yStart = (Subscript[y, 0]), and xFinish = (Subscript[x, end]) for the initial value problem y((Subscript[x, 0])) = (Subscript[y, 0]), y'=f(x,y).)  n is the number of steps, and f is the function on the right-hand side of the ODE in standard form.  Some examples are given below the definitions.

EulerMethod[xStart_,yStart_,xFinish_,n_,f_]:=
Module[{h,xCurrent,yCurrent,StepTable}, (* temporary variables *)
h=(xFinish-xStart)/n; (* Stepsize = h *)
xCurrent=xStart; 
yCurrent=yStart;
StepTable={{xStart,yStart}};
(* now the while-loop will crank through Euler's Method *) 
While[xCurrent<xFinish-h/2,
(* stepping forward in y: *)
yCurrent=yCurrent+h*(f/.{x->xCurrent,y->yCurrent}); 
(* stepping forward in x: *)
xCurrent=xCurrent+h;
StepTable = Join[StepTable,{{xCurrent,yCurrent}}];
];
StepTable (* returns the y-value at xFinish *)
];

ImprovedEulerMethod[xStart_,yStart_,xFinish_,n_,f_]:=
Module[{h,xCurrent,yCurrent,k1,k2,StepTable}, (* temporary variables *)
h=(xFinish-xStart)/n; (* Stepsize = h *)
xCurrent=xStart; 
yCurrent=yStart;
StepTable={{xStart,yStart}};
While[xCurrent<xFinish-h/2, (* begin the main loop *)
k1 = (f/.{x->xCurrent,y->yCurrent});  (* slope at current point *)
k2 = (f/.{x->xCurrent + h,y->yCurrent+h k1}); (* slope at Euler Method point *)
yCurrent=yCurrent+h*(k1 + k2)/2; 
(* stepping forward in x: *)
xCurrent=xCurrent+h;
StepTable = Join[StepTable,{{xCurrent,yCurrent}}];
];
StepTable (* returns the table of {x, y} values at each step *)
];

RungeKutta4[xStart_,yStart_,xFinish_,n_,f_]:=
Module[{h,xCurrent,yCurrent,k1,k2,k3,k4,StepTable}, (* temporary variables *)
h=(xFinish-xStart)/n; (* Stepsize = h *)
xCurrent=xStart; 
yCurrent=yStart;
StepTable={{xStart,yStart}};
While[xCurrent<xFinish-h/2, (* begin the main loop *)
k1 = (f/.{x->xCurrent,y->yCurrent});
k2 = (f/.{x->xCurrent + h/2,y->yCurrent+h k1/2});
k3 = (f/.{x->xCurrent + h/2,y->yCurrent+h k2/2});
k4 =  (f/.{x->xCurrent + h,y->yCurrent+h k3});
yCurrent=yCurrent+h*(k1 +2  k2+ 2 k3 + k4)/6; 
(* stepping forward in x: *)
xCurrent=xCurrent+h;
StepTable = Join[StepTable,{{xCurrent,yCurrent}}];
];
StepTable (* returns the table of {x, y} values at each step *)
];

FunctionTable[xStart_,xFinish_,n_,f_]:=Table[{xStart+i(xFinish-xStart)/n,f/.{x->xStart+i(xFinish-xStart)/n}},{i,0,n}];

Here is an example of these functions in use for the initial value problem y'[x]==x-y[x], y[0]==1.  First we'll get the exact solution using DSolve.

DSolve[{y'[x]==x-y[x],y[0]==1},y[x],x]

Now lets look at the some approximations from x=0 to x=1. Note that in plugging in for "f" in my functions you don't need to show the argument of y:

EulerMethod[0.,1,1,10,x-y]

RungeKutta4[0.,1,1,10,x-y]

FunctionTable[0.,1,10,(E^(-x)) (2-(E^x)+(E^x) x)]

Let's compare the error at the last value of Euler's method and the Runge-Kutta method.

(E^(-x)) (2-(E^x)+(E^x) x)/.{x->1} (* This is the exact endpoint value: *)

In the expression below, we use the function Last to get the last x- and y-value from the Euler approximations.  The [[2]] at the end picks out the y-value.

LastEulerValue=Last[EulerMethod[0.,1,1,10,x-y]][[2]]

N[(2/E)-LastEulerValue]

LastRungeKutta4=Last[RungeKutta4[0.,1,1,10,x-y]][[2]]

N[(2/E)-LastRungeKutta4]

As another example of manipulating these lists, lets look at the final cumulative error for Euler's method as a function of n, the number of steps.

FinalYFromEuler[n_]:=Last[EulerMethod[0.,1,1,n,x-y]][[2]]

Use the Table command to make a list of error values:

EulerErrorTable=Table[Abs[(2/E)-FinalYFromEuler[n]],{n,10,100}];

ListPlot[EulerErrorTable,PlotJoined->True]

A log-log plot will reveal the order m of the method - the error should start behaving as Error = C*(1/n)^m for some constant C.  After taking logs, this looks like Log(Error) = Log(C) - m Log(n), so the order m is -1 times the slope of a fitting line:

LogLogEulerErrorTable=Table[{Log[n],Log[Abs[(2/E)-FinalYFromEuler[n]]]},{n,10,100}];
ListPlot[LogLogEulerErrorTable,PlotJoined->True]

Exercise

1.  For the initial value problem y'[x] == y[x] Sin[x], y[0]==1, first find the exact solution.  Then make plots of the error versus n (the number of steps) at x=Pi for Euler's method and the fourth-order Runge-Kutta method.  

2.  Repeat the above exercise for the IVP y'[x] == y[x] + Sin[x], y[0] == -.48.