Maple for Euler's Method

Prévia do material em texto

21/04/2017 Maple for Euler's Method
http://homepages.math.uic.edu/~hanson/MAPLE/euler.html 1/3
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;
21/04/2017 Maple for Euler's Method
http://homepages.math.uic.edu/~hanson/MAPLE/euler.html 2/3
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);
21/04/2017 Maple for Euler's Method
http://homepages.math.uic.edu/~hanson/MAPLE/euler.html 3/3
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`);
>