Boundary Value Problems: The Finite Difference Method

Many techniques exist for the numerical solution of BVPs. A discussion of such methods is beyond the scope of our course. However, we would like to introduce, through a simple example, the finite difference (FD) method which is quite easy to implement. Moreover, it illustrates the key differences between the numerical solution techniques for the IVPs and the BVPs.

Let's consider the linear BVP describing the steady state concentration profile C(x) in the following reaction-diffusion problem in the domain $0\leq x\leq 1$. The BVP can be stated as

$$\frac{d^2C}{dx^2} - C = 0$$ C(x=0) = 1

$$\frac{dC}{dx}(x=1) = 0.$$ (30)

The analytical solution to the BVP above is simply given by $C^e(x) = \frac{\exp(2-x)+\exp(x)}{1+e^2}$.

We are interested in solving the above equation using the FD technique. The first step is to partition the domain [0,1] into a number of sub-domains or intervals of length h. So, if the number of intervals is equal to n, then nh = 1. We denote by x_i the interval end points or nodes, with x₁ =0 and x_n+1 = 1. In general, we have x_i = (i-1)h, $i=1,2,\cdots n+1$. Let us denote the concentration at the ith node by C_i. The second step is to express the differential operator d²C/dx² in a discrete form. This can be accomplished using finite difference approximations to the differential operators. In this problem, we will use the approximation

$${(\frac{d^2C}{dx^2})}_i = \frac{C_{i+1}-2C_{i}+ C_{i-1}}{h^2}.$$ (31)

For the discretization of the no flux boundary condition at x=1, we will use the discretization given by

$${(\frac{dC}{dx})}_i = \frac{C_{i+1}-C_{i-1}}{2h}.$$ (32)

The finite difference discretizations given above are referred to as the central difference approximations. The local truncation error (LTE) associated with either of the approximations given above is O(h²).

Let's now derive the discretized equations. First of all, we have two boundary conditions to be implemented. The boundary condition at x=0 gives

C₁ = 1 (33)

The implementation of the no flux condition at x=1 is somewhat tricky. Note that according to Eq. 33, we should write (with i=n+1)

$$\frac{C_{n+2}-C_{n}}{2h}=0 \:\:\: {\mbox{or}} \:\:\: C_{n+2}-C_{n} = 0.$$ (34)

However, C_n+2 is not in our domain thereby creating a difficulty. What do we do here? The solution is to extend the original domain from (0,nh) to (0, (n+1)h). This results in n+2 nodes, the n+2th one being a fictitious one. We solve for $C_1, C_2, \cdots , C_{n+1}$ and the additional variable introduced due to the fictitious node C_n+2 and discard C_n+2 from the final solution.

So far, we have supplied 2 equations for the n+2 unknowns, the remaining n equations are obtained by writing the discretized ODE for nodes $2,3,\cdots , n+1$. Application of Eq. 32 and the use of the boundary conditions lead to the following system of linear equations for C_i, $i=1,2,\cdots , n+2$.

C_n+2 - C_n = 0. (35)

We can express this system compactly using matrices. The $(n+2)\times (n+2)$coefficient matrix, say ${\bf {A}}$, corresponding to the system of equations given above is

$${\bf {A}} = \left( \begin{array}{cccccccc} 1 & 0 & 0 & 0 & \c... ...1 \\ 0 & 0 & 0 & 0 & \cdots & -1 & 0 & 1 \end{array}\right).$$ (36)

The right hand side vector ${\bf {b}}={(1,0,0,\cdots , 0)}^T$. This means that the solution to the system ${\bf {A}}.{\bf {c}} = {\bf {b}}$ gives the solution vector ${\bf {c}} = {(C_1, C_2, \cdots , C_{n+2})}^T$. The system of equations can be solved using Gaussian elimination or more typically using a special linear system solver designed to take advantage of the tridiagonal structure of the coefficient matrix.

In Figure 5, the FD solution with h=0.1 and h=0.05 are presented along with the exact solution to the BVP of Eq. 31. A very good agreement between the exact and the computed solutions can be seen from there.

How does the FD scheme above converge to the exact solution as h is decreased? The absolute error at the center of the domain (x=0.5) for three different values of h are plotted vs. h in Figure 6 on a log-log plot. It can be seen from there that the error decreases as O(h²). Hence, the FD approximation used here has quadratic convergence. This is because the discretization errors in the approximation of the first and second derivative operators (see Eqs. 32 and 33) are O(h²). Indeed, the convergence characteristics can be improved by using more accurate discretization of the differential operators.