35DEmthd.html

Partial Differential Equations PowerTool

by Dr. Jim Herod

Section 8.1: Numerical Solutions for Ordinary Differential Equations: Difference Methods

Maple Packages used in Section 8.1

> restart:

> with(plots):

Warning, the name changecoords has been redefined

> with(LinearAlgebra):

There are many methods for getting numerical solutions for ordinary differential equations. Through the years, with improving analysis, software and hardware, the methods have improved in accuracy. Euler's Method, Improved Euler's Method, and the Runge-Kutta Method are commonly taught in undergraduate classes as methods for obtaining numerical solutions for ordinary differential equations. Maple uses a Runge-Kutta-Fehlberg 4/5th order integration method, among others.

In this Section 8.1, we will follow closely the methods presented by Douglas Meade in Unit 18 of his notes for the Ordinary Differential Equations Powertool at the Maplesoft web site. See

http://www.mapleapps.com/powertools/des/des.shtml

In that material, Meade has a complete section, Part IV, on numerical methods for initial value value problems. Part V in his materials contains numerical methods for boundary value problems.

Of course, in partial differential equations, derivatives are made with respect to more than one variable, so the ability to simply "step across the interval" is missing. The methods shown in this 8th Chapter uses the techniques of finite difference methods . The methods could be called discrete replacement methods for, following the development in Unit 18 of Meade's work, we replace the derivatives by discrete approximations. With these methods, the job of getting an approximation for a solution of a differential equation is changed to a task in linear algebra.

The purpose of reviewing these ordinary differential equations methods is that similar techniques will be expanded for partial differential equations.

It should be no surprise that there are choices on what discrete replacements will be made for the derivatives. For example, there are forward difference replacements and backward difference replacements. We follow Meade as we make choices for the replacements. We use and apply the finite difference method to a two-point boundary value problem

Diff(y,`$`(x,2))+p(x)*Diff(y,x)+q(x)*y = f(x) .

There will be a boundary condition at each endpoint of an interval over which we solve the equation. Meade calls the endpoints a and b . The mesh points will be denoted [ , , ... , ].

We explain now how we will replace and diff(y,`$`(x,2)) .

Replacement for :

Recall Taylor's Series:

y( x + h) = y(x) + h y '(x) + 0 ( )

and

y( x - h) = y(x) - h y '(x) + 0 ( )

Thus,

y '(x) is replaced by the difference (y(x+h)-y(x-h))/h and this is good with order 0 ( ).

Replacement for diff(y,`$`(x,2)) :

Use Taylor's Series again:

y(x + h) = y(x) + h y '(x) + y ''(x) + y '''(x) + 0 ( ).
y(x - h) = y(x) - h y '(x) + y ''(x) - y '''(x) + 0 ( ).

Add these two equations and we we get

y( x + h) - y(x) + y(x - h) - y(x) = y '' (x) + 0 ( ).

Thus,

y ''(x) is replaced by the difference (y(x+h)-2*y(x)+y(x-h))/(h^2) and this is good with order 0 ( ).

We now work three problems. The first two were worked in Meade's Unit 18. We follow his syntax almost in the manner of "cut and paste." In this Section 8.1, before the problems are finished we will have restated each of them as a matrix problem.

Example 1

Consider the boundary value problem

> ode1:=diff(y(x),x$2)+2*diff(y(x),x)-10*y(x)=7*exp(-x)+4;
bc1:=y(0)=2,y(1)=-5;

with boundary conditions specified at the points

> a:=0;
b:=1;

To find the finite difference solution associated with this boundary value problem, Meade chooses 11 mesh points from a to b. He replaces the derivatives as we indicated above and shows graphically what a good approximation this makes to the exact solution. See his Unit 18.

Rather than repeat what he has done, we choose only 6 mesh points and emphasize the structure of the resulting equations to be solved. We will find there is a tridiagonal matrix associated with the analysis. We proceed. Choose N.

> N:=5;

This defines h and enables us to compute the X mesh.

> h:=(b-a)/N;
X:=k->a+k*h;

Next, we define the finite difference replacement terms.

> Yp:=k->(y[k+1]-y[k-1])/2/h;
Ypp:=k->(y[k+1]-2*y[k]+y[k-1])/h^2;

The boundary conditions provide two equations.

> eq[0]:=y[0]=rhs(bc1[1]);
eq[N]:=y[N]= rhs(bc1[2]);

Here are the remaining N-1 equations for the values at the interior mesh points. We make these equations by replacing

y '(x) and y ''(x).

>    for k from 1 to N-1 do
   eq[k]:=eval(ode1,{x=X(k),y(x)=y[k],
        diff(y(x),x)=Yp(k),diff(y(x),x$2)=Ypp(k)});
od;

>

This system of N+1 equations can be realized in a tridiagonal matrix equation. The matrix will need N+1 rows. The first one comes from the boundary condition at x = a.

> s[0]:=[1,seq(0,p=1..N)];

The next N-1 rows come from the interior mesh points.

> for n from 1 to N-1 do
s[n]:=[seq(0,p=1..n-1),
coeff(lhs(eq[n]),y[n-1]),
coeff(lhs(eq[n]),y[n]),
coeff(lhs(eq[n]),y[n+1]),
seq(0,p=1..N-1-n)];
od;

Finally, there is the last row that comes from the boundary condition at x = b.

> s[N]:=[seq(0,p=1..N),1];

Now, we write this as a matrix. We list out the rows together.

> rose:=[s[0],seq(s[n],n=1..N-1),s[N]];

We construct the matrix with these rows. Note that A is a tridiagonal matrix.

> A:=Matrix(rose);

Here are the unknown values , ... , .

> Y:=Vector([seq(y[p],p=0..N)]);

Here is the right hand side of the equations written as a column vector.

> V:=Vector([rhs(bc1[1]),seq(rhs(eq[p]),p=1..N-1),rhs(bc1[2])]);

Finally, here is the equation to be solve for y.

> A.Y=V;

There are many ways to solve matrix equations. We leave this up to Maple.

> Y:=LinearSolve(A,V);

In an attempt to make these numbers understandable by humans, we give a listing.

> result:=eval([seq([X(k-1),evalf(Y[k],5)],k=1..N+1)]);

It will be interesting to compare how good this approximation is with the exact solution. The approximation may be improved by rerunning this example and increasing the size of N.
We draw graphs. First, compute the exact solution.

> sol1:=combine(dsolve({ode1,bc1},y(x)));

Next, we superimpose the exact solution with this approximation.

>    plot([rhs(sol1),result],x=a..b,
  color=[BLACK,GREEN],thickness=[1,3],
  legend=[`exact`,`finite difference`]);

Now a bad fit, yes? I say again, you might go back and increase the size of N to improve the approximation.

Example 2

Here is the second example from Meade's Unit 18. We make modifications to see it as a tridiagonal matrix again.

The boundary value problem in this example does not have constant coefficients.

> ode2:=x^2*diff(y(x),x$2)+x*diff(y(x),x)+4*y(x)=-2*x+7;

> bc2:=y(1)=7,y(4)=-1;

The boundary conditions are given at

> a:=1;

> b:=4;

Here is where we choose the number of mesh points. The number of mesh points will be N+1.

> N:=5;

The interval over which we approximate this equation is [a, b] and for this interval, we have defined the step size.

> h:=(b-a)/N;

Here are the mesh points.

> X:=k->a+k*h;

The replacement operators are as they were explained in the discussion above.

> Yp:=k->(y[k+1]-y[k-1])/(2*h);

> Ypp:=k->(y[k+1]-2*y[k]+y[k-1])/h^2;

The boundary conditions provide two equations.

> eq[0]:=y[0]=rhs(bc2[1]);
eq[N]:=y[N]= rhs(bc2[2]);

Here are the remaining N-1 equations for the values at the interior mesh points. We make these equations by replacing

y '(x) and y ''(x).

>    for k from 1 to N-1 do
   eq[k]:=eval(ode2,{x=X(k),y(x)=y[k],
          diff(y(x),x)=Yp(k),diff(y(x),x$2)=Ypp(k)} );
od;

We realize this system of N+1 equations as tridiagonal matrix equation. The matrix will need N+1 rows. The first one comes from the boundary condition at x = a.

> s[0]:=[1,seq(0,p=1..N)];

The next N-1 rows come from the interior mesh points.

> for n from 1 to N-1 do
s[n]:=[seq(0,p=1..n-1),
coeff(lhs(eq[n]),y[n-1]),
coeff(lhs(eq[n]),y[n]),
coeff(lhs(eq[n]),y[n+1]),
seq(0,p=1..N-1-n)];
od;

Finally, there is the last row that comes from the boundary condition at x = b.

> s[N]:=[seq(0,p=1..N),1];

Now, we write this as a matrix. We list out the rows together.

> rose:=[s[0],seq(s[n],n=1..N-1),s[N]];

We construct the matrix with these rows.

> A:=Matrix(rose);

Here are the unknown values , ... , .

> Y:=Vector([seq(y[p],p=0..N)]);

Here is the right hand side of the equations written as a column vector.

> V:=Vector([rhs(bc2[1]),seq(rhs(eq[p]),p=1..N-1),rhs(bc2[2])]);

Finally, here is the equation to be solve for y.

> A.Y=V;

Maple will solve this system for us.

> Y:=LinearSolve(A,V);

We reprint for humans.

> result:=eval([seq([X(k-1),evalf(Y[k],5)],k=1..N+1)]);

Finally, we draw graphs. First, compute the exact solution.

> sol2:=combine(dsolve({ode2,bc2},y(x)));

Next, we superimpose the exact solution with this approximation.

>    plot([rhs(sol2),result],x=a..b,
   color=[BLACK,GREEN], thickness=[1,3],
   legend=[`exact`,`finite difference`]);

The fit with N = 5 is not as good as we would have liked. Going back and increasing N is irresistible.

Example 3

This example differs from the previous two in that the boundary condition at the right hand side involves the derivative of y, instead of the value of y. Also, we make the right hand side of the equation not continuous.

Here is the example: we seek a graph for a function y(x) that satisfies

- 3 y(x) = f(x)

with

boundary conditions y(-1) = 1 and y '(1) = -1. The function f is indicated below.

To the extent possible, we follow the pattern of the previous two examples.

We use a second order differential equation with forcing function f.

> ode3:=diff(y(x),x$2)-9*y(x)=f(x);
bc3:=y(-1)=1,D(y)(1)=-1;

> f:=x->(1-signum(x))/2;

The boundary conditions are given at

> a:=-1;

> b:=1;

Here is where we choose the number of mesh points. The number of mesh points will be N+1.

> N:=9;

The interval over which we approximate this equation is [a, b] and for this interval, we have defined the step size.

> h:=(b-a)/N;

Here are the mesh points.

> X:=k->a+k*h;

The replacement operators are as they were explained in the discussion above.

> Yp:=k->(y[k+1]-y[k-1])/(2*h);

> Ypp:=k->(y[k+1]-2*y[k]+y[k-1])/h^2;

Of course, the boundary condition at a provides an equation of the same character as in the previous two examples. The boundary condition at b needs some discussion.

> eq[0]:=y[0]=rhs(bc3[1]);

>

Here, we explain how to get equation N from the boundary condition at the right hand endpoint, b . Our replacement scheme for the derivative at b = 1 = should be

= -1.

However, the sequence of Y 's runs only from to . Thus we get in the following manner. From the previous line, we get that

= - 2 h .

If we followed the discrete replacement scheme at , we would have

- 9 = f( ).

Again, this has a value of Y , namely , that we do not use. It is replaced by the value of that we got from the boundary condition to yield

- 9 = f( ).
This defines .

> eq[N]:=simplify((y[N-1]-2*y[N]+y[N-1])/h^2-9*y[N]=
f(X(N))+2/h);

Here are the remaining N-1 equations for the values at the interior mesh points. We make these equations by replacing

y '(x) and y ''(x).

>    for k from 1 to N-1 do
   eq[k]:=eval(ode3,{x=X(k),y(x)=y[k],
          diff(y(x),x)=Yp(k),diff(y(x),x$2)=Ypp(k)} );
od;

We realize this system of N+1 equations as tridiagonal matrix equation. The matrix will need N+1 rows. The first one comes from the boundary condition at x = a.

> s[0]:=[1,seq(0,p=1..N)];

The next N-1 rows come from the interior mesh points.

> for n from 1 to N-1 do
s[n]:=[seq(0,p=1..n-1),
coeff(lhs(eq[n]),y[n-1]),
coeff(lhs(eq[n]),y[n]),
coeff(lhs(eq[n]),y[n+1]),
seq(0,p=1..N-1-n)];
od;

Finally, there is the last row that comes from the boundary condition at x = b.

> s[N]:=[seq(0,p=1..N-2),coeff(lhs(eq[N]),y[N-2]),
coeff(lhs(eq[N]),y[N-1]),coeff(lhs(eq[N]),y[N])];

Now, we write this as a matrix. We list out the rows together.

> rose:=[s[0],seq(s[n],n=1..N-1),s[N]];

We construct the matrix with these rows.

> A:=Matrix(rose);

Here are the unknown values , ... , .

> Y:=Vector([seq(y[p],p=0..N)]);

Here is the right hand side of the equations written as a column vector.

> V:=Vector([rhs(bc3[1]),seq(rhs(eq[p]),p=1..N-1),rhs(eq[N])]);

Finally, here is the equation to be solve for y.

> A.Y=V;

Maple will solve this system for us.

> Y:=LinearSolve(A,V);

We reprint for humans.

> result:=eval([seq([X(k-1),evalf(Y[k],5)],k=1..N+1)]);

Finally, we draw graphs. First, compute the exact solution.

> sol3:=combine(dsolve({ode3,bc3},y(x)));

Next, we superimpose the exact solution with this approximation.

>    plot([rhs(sol3),result],x=-1..1,
   color=[BLACK,GREEN], thickness=[1,3],
   legend=[`exact`,`finite difference`]);

>

If you don't think this fit is pretty good, go back and make N larger. Try N = 9.

This review of the finite difference methods, or discrete replacement methods, for boundary value problems in differential equations was made to suggest the techniques for partial differential equations. We examine these next.

EMAIL: herod@math.gatech.edu or jherod@tds.net

URL: http://www.math.gatech.edu/~herod