• Outline of Lesson 18
• Initialization
• 18.A Overview of the Finite Difference Method
• 18.B Example 1
• 18.C Example 2

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Lesson 18 -- Finite Difference Method

Prof. Douglas B. Meade

Industrial Mathematics Institute

Department of Mathematics

University of South Carolina

Columbia, SC 29208

URL:   http://www.math.sc.edu/~meade/

E-mail: meade@math.sc.edu

Copyright  2001  by Douglas B. Meade

All rights reserved

-------------------------------------------------------------------

Outline of Lesson 18

18.A Overview of the Finite Difference Method

18.B Example 1

18.C Example 2

>

Initialization

>

restart;

>	with( DEtools ):

>	with( plots ):

Warning, the name changecoords has been redefined

>	interface( rtablesize=25 );

>

18.A Overview of the Finite Difference Method

The finite difference method for the two-point boundary value problem

    

    

with Dirichlet boundary conditions seeks to obtain approximate values of the solution at a collection of nodes, , in the interval . Typically, the interval is uniformly partitioned into equal subintervals of length . In this case, for = 0, 1, ..., . Let , , denote the value of the solution at each node. Note that the boundary conditions imply and so there are actually only ( ) unknowns. To determine a system of ( ) equations satisfied by , , ..., , replace the first and second derivatives of the unknown function at the interior partion points with the central difference approximations

           

For linear problems, the resulting system of equations is linear. For nonlinear problems, an iterative method such as Gauss-Seidel is typically used to obtain an approximate solution.

For the Neumann boundary conditions

the discretization at each of the nodes yields equations in the unknowns , where the unknowns and are introduced by the discretizations at each of the endpoints.  To obtain two additional equations, apply the central difference approximation for the first derivative to the boundary conditions, resulting in the two additional equations

and

Clearly, the same strategy works for the Robin boundary conditions

as well as the general linear boundary conditions

or for any combination of one type condition at one end and a different type boundary condition at the other end.

>

18.B Example 1

Consider the BVP

>	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 to this BVP with

>

N := 10;

>

subdivisions, we have the nodal spacing

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

>

and the nodes

>	X   := k -> a+k*h: 'X'[k] = X(k);

>

To simplify the task of writing the nodal equations, define the finite-difference operators

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

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

>

>

so that the nodal equations at the interior nodes are

>	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)} );

>

end do;

>

The boundary conditions provide the two equations

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

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

>

Even though this system of equations is linear, the exact solution is so cumbersome it's best to obtain an approximate (numeric) solution via fsolve .  Thus, the nodal values are computed to be

>	fd_sol1 := fsolve( {seq( eq[k], k=0..N )}, {seq( y[k], k=0..N )} );

>

These results are reformatted as points (for plotting) via

>	fd_table1 := eval( seq([X(k),y[k]],k=0..N), fd_sol1 );

>

and as a table, via

>	Matrix([fd_table1]);

>

Because of the simple nature of this example, Maple is able to determine that the exact solution of this BVP is

>	infolevel[dsolve] := 3:

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

>	infolevel[dsolve] := 0:

`Methods for second order ODEs:`
`--- Trying classification methods ---`
`trying a quadrature`
`trying high order exact linear fully integrable`
`trying differential order: 2; linear nonhomogeneous with symmetry [0,1]\
`
`trying a double symmetry of the form [xi=0, eta=F(x)]`
`Try solving first the homogeneous part of the ODE`
`   -> Tackling the linear ODE "as given":`
`      checking if the LODE has constant coefficients`
`      <- constant coefficients successful`
`   <- successful solving of the linear ODE "as given"`
`   -> Determining now a particular solution to the non-homogeneous ODE`\
`      building a particular solution using variation of parameters`

`   <- solving first the homogeneous part of the ODE successful`

>

Figure 18.1 is used to compare the exact and finite-difference solutions.

>	plot( [rhs(exact_sol1),[fd_table1]], x=a..b,

>	color=[BLACK,GREEN], style=[line,point], symbolsize=16,

>	legend=[`exact`,`finite difference`], title="Figure 18.1" );

>

This also shows good agreement between the finite-difference and exact solutions.

>

18.C Example 2

Consider the BVP

>	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;

>

The finite difference method with

>

N := 10;

>

subdivisions means the nodal spacing is

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

>

and the nodes are

>	unassign('k'); X   := k -> a+k*h: 'X'[k] = X(k);

>

Again defining the finite-difference operators

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

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

>

the nodal equations at the interior nodes   are

>	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)} );

>

end do;

>

The boundary conditions provide the two equations

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

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

>

Even though this system of equations is linear, the exact solution is so cumbersome it's best to obtain an approximate (numeric) solution via fsolve .  Thus, the nodal values are computed to be

>	fd_sol2 := fsolve( {seq( eq[k], k=0..N )}, {seq( y[k], k=0..N )} );

>

These results are reformatted as points (for plotting) via

>	fd_table2 := eval( seq([X(k),y[k]],k=0..N), fd_sol2 );

>

and as a table, via

>	Matrix( [fd_table2] );

>

Because of the simple nature of this example, Maple is able to determine that the exact solution of this BVP is

>	infolevel[dsolve] := 3:

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

>	infolevel[dsolve] := 0:

`Methods for second order ODEs:`
`--- Trying classification methods ---`
`trying a quadrature`
`trying high order exact linear fully integrable`
`trying differential order: 2; linear nonhomogeneous with symmetry [0,1]\
`
`trying a double symmetry of the form [xi=0, eta=F(x)]`
`Try solving first the homogeneous part of the ODE`
`   -> Tackling the linear ODE "as given":`
`      checking if the LODE has constant coefficients`
`      checking if the LODE is of Euler type`
`      <- LODE of Euler type successful`
`   <- successful solving of the linear ODE "as given"`
`   -> Determining now a particular solution to the non-homogeneous ODE`\
`      trying a rational particular solution`
`      <- rational particular solution successful`
`   <- solving first the homogeneous part of the ODE successful`

>

Figure 18.2 is used to compare the exact and finite-difference solutions.

>	plot( [rhs(exact_sol2),[fd_table2]], x=a..b,

>	color=[BLACK,GREEN], style=[line,point], symbolsize=16,

>	legend=[`exact`,`finite difference`], title="Figure 18.2" );

>

This also shows good agreement between the finite-difference and exact solutions.

>

[Back to ODE Powertool Table of Contents]