Topic 15.1: Finite-Difference Method (Theory)

Contents Previous Chapter Start of Chapter Previous Topic Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple Next Topic Next Chapter

Given the linear 2nd-order ordinary differential equation

c2(t) y(2)(t) + c1(t) y(1)(t) + c0(t) y(t) = g(t)

where the ci(t) and g(t) are known functions, a boundary-value problem is a problem where the solution y(t) is subject to two constraints, each given at a different value of t:

y(a) = ya
y(b) = yb

Thus, at each interior point on [a, b], the solution y(t) must satisfy the differential equation. We can break the interval [a, b] into n sub-intervals by defining h = (b − a)/n and letting ti = a + ih. Therefore, we see that:

t0 = a
tn = b

We would like to find approximations of y(ti), so let us represent the approximation of this by yi ≈ y(ti). From the initial conditions, we already have two such values:

y0 = ya
yn = yb

Recall from Topic 13 (Differentiation) that we may approximate the first and second derivatives by the centred divided-difference formulae:

If we substitute these two approximations into our ODE, we get its a discrete form, that is, a difference equation:

If we multiply through by 2h2 and collect with respect to the three values yi + 1, yi, and yi − 1, we get the form:

If we evaluate this for each value of i = 1, 2, ..., n − 1, we get n − 1 equations:

Notice that in the first and last equations, y0 and yn are known. Therefore, we have n − 1 unknowns, namely, y1, ..., yn − 1.

We may solve this system of linear equations using any of the tools we have available from linear algebra.

Special Case: Constant Coefficients

In the special case where the ODE has constant coefficients, this greatly simplifies the definition of the system of linear equations. If the matrix is defined as

c2 y(2)(t) + c1 y(1)(t) + c0 y(2)(t) = g

Writing out the equations, we note that all of the coefficients are identical:

Therefore, if we define the four variables l, d, u, and v (for lower, diagonal, upper, and vector, respectively):

l = 2c2hc1
d = 2h2c0 − 4c2
u = 2c2 + hc1
v = 2h2g

We can rewrite the equations as:

       dy1 + uy2 = vly0
ly1 + dy2 + uy3 = v
ly2 + dy3 + uy4 = v
              .
              .
              .
lyn − 2 + dyn − 1 = vuyn

We can write this as the system of linear equations:

which we can easily solve for y.

Once we have the vector y, we have approximations for each of the interior points y1, ..., yn − 1. If we needed a point between these, we could use interpolation.

Copyright ©2005 by Douglas Wilhelm Harder. All rights reserved.