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`);**

`> `