ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 18 -- Finite Difference Method

Outline of Unit 18

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

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

.

seeks to obtain approximate values of the solution at a collection of points, , in the interval . Typically, the interval is uniformly partitioned into equal subintervals of length . In this case, for = 0, 1, ..., . Let , = 0, 1, ..., , 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 centered difference quotient 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.

>

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:

>

we have

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

>

and

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

>

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;

>

The boundary conditions provide the two conditions

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

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

>

and the equations for the values 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;

>

Even though this system of equations is linear, it's simpler to use fsolve to find an approximate solution.

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

>

Reformatting these results for plotting

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

>

and for tabular presentation

> stackmatrix( fd_table1 );

>

Because of the simple nature of this example, Maple is able to determine the exact solution to this problem.

> infolevel[dsolve] := 3:

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

> infolevel[dsolve] := 0:

Methods for second order ODEs:

Trying to isolate the derivative d^2y/dx^2...

Successful isolation of d^2y/dx^2

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

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

>

To compare the solutions, plot the exact and finite difference solutions

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

> color=[BLACK,GREEN], thickness=[1,3],

> legend=[`exact`,`finite difference`] );

>

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

>

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;

>

means the stepsize is

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

>

and

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

>

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;

>

The boundary conditions provide the two conditions

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

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

>

and the equations for the values 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;

>

Even though this system of equations is linear, it's simpler to use fsolve to find an approximate solution.

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

>

Reformatting these results for plotting

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

>

and for tabular presentation

> stackmatrix( fd_table2 );

>

Again, because of the simple nature of this example, Maple is able to determine the exact solution to this problem.

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

>

To compare the solutions, plot the exact and finite difference solutions

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

> color=[BLACK,GREEN], thickness=[1,3],

> legend=[`exact`,`finite difference`] );

>

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

>