Finite Difference Solving Higher Order, Ordinary Differential Equations

Nathan Collier, Autar Kaw
University of South Florida
United States of America
kaw@eng.usf.edu

Introduction

The following worksheet demonstrates how the finite difference method may be used to solve higher order differential equations.

Section 1: Input

This differential equation used here to illustrate this method is of the form 
The user may specify values for A1, A2, A3, and A4 which may be functions of x.

A1[x_] := 1;
A2[x_] := 2;
A3[x_] := 3;
A4[x_] := 4*x;

The boundary condition must be specified.

x0 = 0.0;
y0 = 10.0;

xf = 1;
yf = 12;

Section 2: Exact Solution

The following uses Maple's built-in tools to find the exact solution for the ODE.

soln = 
      DSolve[{A1[x]*y''[x] + A2[x]*y'[x] + A3[x]*y[x] + 
              A4[x] == 0, y[xf] == yf, y[x0] == y0}, y, 
        x]

plot1 = 
        Plot[Evaluate[y[x]  /.  soln], {x, x0, xf}, 
          PlotStyle -> {{Thickness[0.012], 
                RGBColor[1, 0, 0]}}];

Section 3: Approximate Solution

This section uses the finite difference method to approximate the solution of the differential equation. The solution will be found by dividing the interval into 5 pieces. Here we find the length of each piece.

h = (xf - x0)/5.0

We also need to setup the X locations going to used.

For[i = 0, 
      i <= 5, i++, 
X[i] = x0 + i*h]

Now we setup an equation for each division and use approximations for the derivatives found in the differental equation. The y values at the beginning and end nodes are known due to the boundary conditions. The approximations used for the derivatives are (y[i+1]-2*y[i]+y[i-1])/(h^2) for the second derivative term and (y[i+1]-y[i-1])/(2*h) for the first derivative term. The following are the system of equations set up which are written for each division we made in the solution.

Equation 1: boundary condition at node 0

eq1 := Y[0]

eq2 := 
      A1[X[1]]*(Y[2] - 2*Y[1] + Y[0])/h^2 + 
        A2[X[1]]*(Y[2] - Y[0])/2/h + A3[X[1]]*Y[1] + 
        A4[X[1]]

eq3 := 
      A1[X[2]]*(Y[3] - 2*Y[2] + Y[1])/h^2 + 
        A2[X[2]]*(Y[3] - Y[1])/2/h + A3[X[2]]*Y[2] + 
        A4[X[2]]

eq4 := 
      A1[X[3]]*(Y[4] - 2*Y[3] + Y[2])/h^2 + 
        A2[X[3]]*(Y[4] - Y[2])/2/h + A3[X[3]]*Y[3] + 
        A4[X[3]]

eq5 := 
      A1[X[4]]*(Y[5] - 2*Y[4] + Y[3])/h^2 + 
        A2[X[4]]*(Y[5] - Y[3])/2/h + A3[X[4]]*Y[4] + 
        A4[X[4]]

eq6 := Y[5]

soln := 
      NSolve[{eq1 == y0, eq2 == 0, eq3 == 0, 
          eq4 == 0, eq5 == 0, eq6 == yf}, {Y[0], Y[1], 
          Y[2], Y[3], Y[4], Y[5]}]

YY[0] = Y[0]  /.  soln[[1]]
YY[1] = Y[1]  /.  soln[[1]]
YY[2] = Y[2]  /.  soln[[1]]
YY[3] = Y[3]  /.  soln[[1]]
YY[4] = Y[4]  /.  soln[[1]]
YY[5] = Y[5]  /.  soln[[1]]

Graphical depiction of the solution

data = Table[{X[i], YY[i]}, {i, 0, 5}];
plot2 = 
        ListPlot[
          data, 
PlotJoined -> 
            True, 
PlotStyle -> {Thickness[0.012], 
              RGBColor[0, 1, 0]}, 
DisplayFunction -> 
            Identity];

Show[plot1, plot2]

References

Chapra & Canale, Numerical Methods for Engineering: With Programming Applications. 3rd ed. WCB/McGraw-Hill. 1998.

Conclusions

Finite differences is a powerful way of solving higher order differential equations. Although accurate, it requires a way of solving systems of equations which is a time consuming operation. The other methods (Runge-Kutta and Euler) do not involve solving a system of equations.