Approximation Methods

MA 507

October 25, 2000

Richard Hitt

Introduction

We look at methods of numerically approximating solutions to IVPs of the form , . For a concrete example we will use the exponential decay differential equation with initial condition .

Euler's Method

> restart;

> f:=(x,y)->y;

>

> y||0:=1;

>

If we want to approximate the solution of this IVP when t=2, we can partition the interval [0,1] into N equal subintervals and use the linear shooting method to move in steps out to a t-value of 2. Using N=10 we first calculate the t-values

> N:=10; a:=0; b:=2; h:=evalf((b-a)/N);
t||0:=a;
for i from 1 to N do
t||i:=h + t||(i-1);
od;

>

We can now successively approximate the value of that correspond to the value of using the tangent line approximation based at the point ( ).

> y||0:=1;
for i from 1 to N do
y||i:= y||(i-1)+h*f(t||(i-1),y||(i-1));
od;

>

To get a visualization for this approximation, we plot the points.

> points:=[seq([t||i,y||i],i=0..N)];

> plot(points,style=point,symbol=cross,view=[0..b,0..7.4]);

>

We can compare the approximations to the actual solutions in this case since we can solve the IVP explicitly.

> dsolve({diff(y(t),t)=f(t,y(t)), y(0)=y||0},y(t));

> exactsol:=unapply(rhs(%),t);

> J:=plot(points,style=POINT,symbol=CROSS,view=[0..b,0..exactsol(b)]):
K:=plot(exactsol(t),t=0..b,color=blue):
plots[display]({J,K});

>

Runge-Kutta Method

The Runge-Kutta method provides us with much more accuracy to the point where we can use a much smaller value of N and still be a better result than with the Euler method.

> N:=4; a:=0; b:=2; h:=evalf((b-a)/N); t||0:=a;
for i from 1 to N do
t||i:=h + t||(i-1);
od;

> y||0:=1:
for i from 1 to N do
k1 := f(t||(i-1),y||(i-1)):
k2 := f(t||(i-1)+h/2,y||(i-1)+h/2*k1):
k3 := f(t||(i-1)+h/2,y||(i-1)+h/2*k2):
k4 := f(t||i,y||(i-1)+h*k3):
y||i:= y||(i-1) + h/6*( k1 + 2*k2 + 2*k3 + k4);
od;

> rkpoints :=[seq([t||i,y||i],i=0..N)];

We can plot the exact solution together with the Euler and Runge-Kutta approximations together. Note how the Runge-Kutta does a much better job in this case even though we are using fewer points.

> L:=plot(rkpoints,style=POINT,symbol=DIAMOND,
view=[0..b,0..exactsol(b)],color=red):
plots[display]({J,K,L});

>

>

Maple has numerical methods built-in

Rather than doing the numerical methods by hand (well almost) as above, we can take advantage of Maple's built-in ability to handle this.

> restart:
with( DEtools ):
with( plots ):

Warning, the name changecoords has been redefined

Let's work on the IVP

,

We can use the Euler method, the improved Euler method (Heun method), and the Runge-Kutta method and compare the results.

> eul :=dsolve(
{diff(y(t),t) = t + y(t), y(0)=1}, y(t),
type=numeric, method=classical[foreuler], stepsize=.33
);

> heun :=dsolve(
{diff(y(t),t)=t+y(t), y(0)=1}, y(t),
type=numeric, method=classical[heunform], stepsize=.33
);

> rk :=dsolve(
{diff(y(t),t)=t+y(t), y(0)=1}, y(t),
type=numeric, method=classical[rk4], stepsize=.33
);

> plot1 := odeplot(eul,[t,y(t)],color=red):
plot2 := odeplot(heun,[t,y(t)],color=green):

> plot3 := odeplot(rk,[t,y(t)],color=blue):

> display(plot1,plot2,plot3,view=[0..3,0..30]);

You can experiment with the stepsize on the Euler method to see how small you have to make it to get an approximation to the solution that is as good as the Runge-Kutta method. Note that the Heun method gives much better results in this case that the regular Euler method.

The "method=classical" option used in the dsolve command instructs Maple to use a static (fixed) step size. Left to its defaults, Maple would employ an adaptive stepsize making is smaller as needed to gain accuracy.

When Things Go Wrong

An Example

Here is a simple example that illustrates the possible pathology in numerical solutions.

We solve the IVP , . We can get an exact solution using dsolve.

> sol:=dsolve({2*diff(y(t),t)*y(t)=-1, y(0)=1},y(t));

> plot(rhs(sol),t=0..2);

If we use DEplot to get a plot of a numerical solution, the default method Maple uses is classical[rk4]. This method is not sensitive to vertical tangents on the solution, so it blindly computes right across them producing potentially rediculous results.

> DEplot(2*diff(y(t),t)*y(t)=-1, y(t), t=0..2, [[0,1]]);

We can get slightly less rediculous results by reducing the stepsize.

> DEplot(2*diff(y(t),t)*y(t)=-1, y(t), t=0..2, [[0,1]], stepsize=0.01);

But for better results, we need to let Maple use an adaptive stepsize.

> DEplot(2*diff(y(t),t)*y(t)=-1, y(t), t=0..2, [[0,1]], method=rkf45);

Error, (in DEtools/DEplot/drawlines) Stopping integration due to rkf45 is unable to achieve requested accuracy

The numerical method recognizes that it cannot subdivide the stepsize enough to achieve its default accuracy and simply reports that rather than drawing an erroneous graph.

>

>

Another Example

> DEplot(diff(y(x),x) = x^2/(1-y(x)^2),y(x),x=-2..2,y=-2..2,[[0,2]]);

>

Another Example

> ode:=diff(y(t),t$2)+y(t)=0;

> plot1:=DEplot(ode,y(t),t=0..100,[[y(0)=1,D(y)(0)=0]],stepsize=.6,color=red):

> plot2:=plot(rhs(dsolve({ode,y(0)=1,D(y)(0)=0},y(t))),t=0..100,color=blue):

> display({plot1,plot2});

>