Question 6.2: Show how the error in the finite difference approximation de...
Show how the error in the finite difference approximation depends on the number of nodes n for the reaction–diffusion equation
\frac{d^2c}{dx^2} = \phi ^2c (6.3.3)
with c(0) = 0 and c(1) = 1. In chemical engineering lingo, the quantity φ is the Thiele modulus for this problem.
This problem is solved analytically using the characteristic polynomial,
r^2 − \phi ^2 = 0 (6.3.4)
which has roots r = ±\phi . The solution is then
c(x) = a \sinh \phi x + b \cosh \phi x (6.3.5)
From the left boundary condition we get b = 0 since sinh 0 = 0. From the right boundary condition we get a^{−1} = \sinh \phi . Thus, the exact solution is
c(x) =\frac{ \sinh \phi x}{\sinh \phi } (6.3.6)
To solve this problem by finite differences, we approximate the second derivative with finite differences,
\frac{c_{i+1} − 2c_i + c_{i−1}}{Δx^2} − \phi ^2c_i = 0 (6.3.7)
For the boundary conditions, we have c_1 = 0 and c_n = 1.
This leads to a problem of the form
Ac = b (6.3.8)
or, more specifically,
\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & ··· & 0& 0 & 0 \\ 1 & α & 1 & 0 & 0 & ··· & 0 & 0 & 0\\ 0 & 1 & α & 1 & 0 & ··· & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 0 & ··· & 1 & α & 1 \\ 0 & 0 & 0 & 0 & 0 &···& 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3\\\vdots \\c_{n−1} \\c_n \end{bmatrix} =\begin{bmatrix} 0 \\ 0 \\ 0\\\vdots \\0 \\1 \end{bmatrix} (6.3.9)
where α = −2 − \phi ^2 Δx^2. We can then solve this system using any of the methods from Chapter 2. Since the finite difference method leads to a tridiagonal band structure (p = q = 1) for Dirichlet boundary conditions, we should take advantage of a banded solver.
Let’s take a look at the MATLAB program to solve this problem.
1 function c = bvp_linear(n,L)
2 dx = 1/(n-1);
3 alpha = -2 - (L*dx)^2;
4 A = zeros(n);
5 b = zeros(n,1);
6 A(1,1) = 1;
7 for i = 2:n-1
8 A(i,i-1) = 1;
9 A(i,i) = alpha;
10 A(i,i+1) = 1;
11 end
12 A(n,n) = 1;
13 b(n) = 1;
14 A = sparse(A);
15 c = A\b;
The program is just setting up the system of linear equations produced by the finite difference method and then solving them using a banded solver. Note that lines 4 and 5 declare the matrix A and the vector b to be all zeros. We thus only have to fill in the entries that are non-zero. For the left boundary condition, we have
A_{1,1} = 1, b_1 = 0 (6.3.10)
since the corresponding equation is c_1 = 0. Since we already set all the entries in b to be zero at the start, we do not need to reassign them to zero here. For the interior nodes, we have the repeating pattern [1, α, 1]. The for loop in lines 7–11 writes these parts of the matrix. The differential equation is homogeneous, so again we do not need to do anything about the vector b for the interior nodes. The only time we have to address the vector b is in the nth row, where we set it equal to 1 in line 13. Note that line 14 declares the matrix A to be sparse. This will force MATLAB to use the banded solver.
Figure 6.2 compares the exact solution in Eq. (6.3.6) to the numerical solution. Even with just a few nodes, we get a reasonably good prediction of the solution.
