Numerical Solution of ODE - Euler's Method and Improved Euler
ODE- IVP : y' = f(x,y), y(0) = 1
Goal: Generate a direction field graph and
compare the exact solution to 2 numerical approximations.
Turn in graphs for the following the IVP that illustrate the use of
two numerical methods and the exact solution along with the direction
field plot.
Computer Problem 1: f(x,y) = x- y^2
> restart;
Define the function f(x,y) in the RHS of ODE
> f := (x,y) -> y;
Generate a directionfield plot and sketch the solution if
x(0)=0 and y(0)=1
>
DEtools[DEplot](diff(y(x),x)=f(x,y(x)),y(x),x=0..8,
y=0..3,arrows=slim,title=`Your Name Here`);
Numerical Approximations
The sequence of approximations is (x[n],y[n]), n=0,1,2,...
Define initial conditions and step size
> x[0]:= 0; y[0] := 1.0; h:= 0.1;
Euler's method or algorithm:
>
for n from 1 to 20 do
x[n] := n*h;
y[n] := y[n-1] + h*f(x[n-1],y[n-1]);
od:
The next command generates sequence of of approximatons:
> data := [seq([x[n],y[n]],n=0..20)]:
Generate plot which is not displayed but instead stored under the name t1
> t1 := plot(data,style=point,color=red):
Improved Euler's method - use different names for approximation
>
xx[0]:=0: yy[0]:=1:
for n from 1 to 20 do
xx[n] := n*h;
ystar := y[n-1] + h*f(x[n-1],y[n-1]);
yy[n] := yy[n-1] +
h/2.0*(f(xx[n-1],yy[n-1])+f(xx[n],ystar)):
od:
Generate the sequence of approximations for the Improved Euler method
> data_improve := [seq([xx[n],yy[n]],n=0..20)]:
Generate plot which is not displayed but instead stored under the name t3
> t3 := plot(data_improve,style=point,color=blue):
Now have Maple construct the exact solution, if possible.
We write symbolic form of ODE using a new function u(x)
> eqn := diff(u(z),z) = f(z,u(z));
Try to construct exact solution of IVP using dsolve
>
### WARNING: `dsolve` has been extensively rewritten, many new result forms can occur and options are slightly different, see help page for details
dsolve({eqn,u(0)=1},u(z));
We need to define u as a function
> u := unapply(rhs(%),z);
Graph without displaying the exact solution
> t2 := plot(u(z),z=0..2):
Plot three graphs on same axes but insert your name as the title.
> plots[display]({t1,t2,t3},title=`your name`);
>