5.3. The finite-difference method

The finite-difference method for solving a boundary value problem replaces the derivatives in the ODE with finite-difference approximations derived from the Taylor series. This results in linear system of algebraic equations that can be solved to give an approximation of the solution to the BVP.

The derivatives of \(y_i\) are approximated using values of neighbouring nodes \(y_{i+1}\), \(y_{i-1}\) etc. using expressions derived by truncating the Taylor series and rearranging to make the derivative term the subject. Some common finite-difference approximations are tabulated below.

Finite-difference approximation	Error	Name
\(y' = \dfrac{y_{i+1} - y_i}{h}\)	\(O(h)\)	Forward difference
\(y' = \dfrac{y_i - y_{i-1}}{h}\)	\(O(h)\)	Backward difference
\(y' = \dfrac{y_{i+1} - y_{i-1}}{2h}\)	\(O(h^2)\)	Central difference
\(y'' = \dfrac{y_{i-1} - 2 y_i + y_{i+1}}{h^2}\)	\(O(h^2)\)	Symmetric difference

The solution to a boundary value problem using then finite-difference method is determined by approximating the derivatives in the ODE using finite-differences. Consider the following boundary value problem

\[ \begin{align*} y'' = f(t,y), \qquad t \in [t_0, t_{\max}], \qquad y(t_0) = a, \qquad y(t_{\max}) = b. \end{align*} \]

Using the symmetric difference to approximate \(y''\) in the ODE we have

\[\begin{split} \begin{align*} \frac{y_{i-1} - 2y_i + y_{i+1}}{h^2} &= f(t_i ,y_i) \\ y_{i-1} - 2y_i + y_{i+1} &= h^2f(t_i, y_i). \end{align*} \end{split}\]

Since we know that \(y_0 = a\) and \(y_n = b\) then

\[\begin{split} \begin{align*} y_0 &= a, \\ y_0 - 2y_1 + y_2 &= h^2 f(t_1, y_1), \\ y_1 - 2y_2 + y_3 &= h^2 f(t_2, y_2), \\ & \vdots \\ y_{n-3} - 2y_{n-2} + y_{n-1} &= h^2 f(t_{n-2}, y_{n-2}), \\ y_{n-2} - 2y_{n-1} + y_{n} &= h^2 f(t_{n-1}, y_{n-1}), \\ y_{n} &= b, \end{align*} \end{split}\]

which can be written as the matrix equation

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 \\ 1 & -2 & 1 \\ & 1 & -2 & 1 \\ & & \ddots & \ddots & \ddots \\ & & & 1 & -2 & 1 \\ & & & & 1 & -2 & 1 \\ & & & & & 0 & 1 \end{pmatrix} \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-2} \\ y_{n-1} \\ y_n \end{pmatrix} = \begin{pmatrix} a \\ f(t_1, y_1) \\ f(t_2, y_2) \\ \vdots \\ f(t_{n-2}, y_{n-2} \\ f(t_{n-1}, y_{n-1}) \\ b \end{pmatrix}. \end{align*} \end{split}\]

---

5.3.1. The Thomas algorithm

The coefficient matrix for the linear system that results from the application of the second-order symmetric finite-difference approximation of \(y''\) is tri-diagonal matrix (the only non-zero elements are on the main, lower and upper diagonals). The usual methods for solving linear systems of equations such as Gaussian elimination and LU decomposition could be applied here, however there is a more efficient method that can solve tri-diagonal systems called the Thomas algorithm).

Definition 5.2 (The Thomas algorithm)

Given a tri-diagonal linear system of equations of the form \(A \mathbf{x} = \mathbf{b}\)

\[\begin{split} \begin{align*} \begin{pmatrix} b_0 & c_0 \\ a_1 & b_1 & c_1 \\ & \ddots & \ddots & \ddots \\ & & a_{n-1} & b_{n-1} & c_{n-1} \\ & & & a_n & b_n \end{pmatrix} \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1} \\ x_n \end{pmatrix} = \begin{pmatrix} d_0 \\ d_1 \\ \vdots \\ d_{n-1} \\ d_n \end{pmatrix}, \end{align*} \end{split}\]

the solution can be found by performing a forward sweep on the \(b_i\) and \(d_i\) coefficients

\[\begin{split} \begin{align*} b_i &= b_i - c_{i-1} \left( \frac{a_i}{b_{i-1}} \right), & i = 1, 2, \ldots, n, \\ d_i &= d_i - d_{i-1} \left( \frac{a_i}{b_{i-1}} \right), & i = 1, 2, \ldots, n, \end{align*} \end{split}\]

and then calculate the values of \(x_i\) using back substitution

\[\begin{split} \begin{align*} x_n &= \frac{d_n}{b_n}, \\ x_i &= \frac{d_i - c_i x_{i+1}}{b_i}, \qquad i = n-1, n-2, \ldots, 1. \end{align*} \end{split}\]

Code

def tridiagonal_solver(A, b):

    n = len(b)
    for i in range(1, n):
        A[i,i] -= A[i-1,i] * A[i,i-1] / A[i-1,i-1]
        b[i] -= b[i-1] * A[i,i-1] / A[i-1,i-1]
    
    x = np.zeros(n)
    x[-1] = b[-1] / A[-1,-1]
    for i in range(n-2, -1, -1):
        x[i] = (b[i] - A[i,i+1] * x[i+1]) / A[i,i]

    return x

Example 5.4

Use the finite-difference method to solve the boundary below using a step length of \(h = 0.2\)

\[ \begin{align*} y'' - y' - y = 0, \qquad y(0) = 0, \qquad y(1) = 2. \end{align*} \]

---

Solution

Using a forward difference to approximate \(y'\) and a symmetric difference to approximate \(y''\) we have

\[ \begin{align*} \frac{y_{i-1} -2y_i +y_{i+1} }{h^2 }-\frac{y_{i+1} -y_i }{h}-y_i = 0, \end{align*} \]

which can be simplified to

\[ \begin{align*} y_{i-1} +(h-h^2 -2)y_i +(1-h)y_{i+1} = 0. \end{align*} \]

This gives the following system of linear equations

\[\begin{split} \begin{align*} y_0 &= 0, \\ y_0 + (h - h^2 - 2) y_1 + (1 - h) y_2 &= 0, \\ y_1 + (h - h^2 - 2) y_2 + (1 - h) y_3 &= 0, \\ &\vdots \\ y_{n-3} + (h - h^2 - 2) y_{n-2} + (1 - h) y_n &= 0, \\ y_{n-2} + (h - h^2 - 2) y_{n-1} + (1 - h) y_n &= 0, \\ y_n &= 2. \end{align*} \end{split}\]

which we can write the linear system as the matrix equation

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 & & & & \\ 1 & h-h^2 -2 & 1-h & & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & h-h^2 -2 & 1-h \\ & & & 0 & 1 \end{pmatrix} \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \\ y_n \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 2 \end{pmatrix}. \end{align*} \end{split}\]

If we use a step length of \(h=0.2\) we have

\[ \begin{align*} \vec{t} = (0, 0.2, 0.4, 0.6, 0.8, 1), \end{align*} \]

and

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 \\ 1 & -1.84 & 0.8 \\ & 1 & -1.84 & 0.8 \\ & & 1 & -1.84 & 0.8 \\ & & & 1 & -1.84 & 0.8 \\ & & & & 0 & 1 \end{pmatrix} \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 2 \end{pmatrix}. \end{align*} \end{split}\]

Performing the forward sweep of the Thomas algorithm

\[\begin{split} \begin{align*} b_1 &= 1, \\ b_2 &= b_2 - c_1 \left( \frac{a_2}{b_1} \right) = -1.84 - 0 \left( \frac{1}{1} \right) = -1.84, \\ b_3 &= b_3 - c_2 \left( \frac{a_3}{b_2} \right) = -1.84 - 0.8 \left( \frac{1}{-2.64} \right) = -1.405217, \\ b_4 &= b_4 - c_3 \left( \frac{a_4}{b_3} \right) = -1.84 - 0.8 \left( \frac{1}{-1.405217} \right) = -1.270693, \\ b_5 &= b_5 - c_4 \left( \frac{a_5}{b_4} \right) = -1.84 - 0.8 \left( \frac{1}{-1.270693} \right) = -1.210422, \\ b_6 &= b_6 - c_5 \left( \frac{a_6}{b_5} \right) = 1 - 0.8 \left( \frac{0}{-1.210422} \right) = 1, \\ d_1 &= 0, \\ d_2 &= d_2 - d_1 \left( \frac{a_2}{b_1} \right) = 0 - 0 \left( \frac{1}{1} \right) = 0, \\ d_3 &= d_3 - d_2 \left( \frac{a_3}{b_2} \right) = 0 - 0 \left( \frac{1}{-2.64} \right) = 0, \\ d_4 &= d_4 - d_3 \left( \frac{a_4}{b_3} \right) = 0 - 0 \left( \frac{1}{-1.405217} \right) = 0, \\ d_5 &= d_5 - d_4 \left( \frac{a_5}{b_4} \right) = 0 - 0 \left( \frac{1}{-1.270693} \right) = 0, \\ d_6 &= d_6 - d_5 \left( \frac{a_6}{b_5} \right) = 2 - 0 \left( \frac{0}{-1.270693} \right) = 2, \\ \end{align*} \end{split}\]

Performing the back substitution

\[\begin{split} \begin{align*} y_6 &= \frac{d_6}{b_6} = \frac{2}{1} = 2, \\ y_5 &= \frac{d_5 - c_5 y_6}{b_5} = \frac{0 - 0.8(2)}{-1.210422} = 1.321853, \\ y_4 &= \frac{d_4 - c_4 y_5}{b_4} = \frac{0 - 0.8(1.321853)}{-1.270693} = 0.832209, \\ y_3 &= \frac{d_3 - c_3 y_4}{b_3} = \frac{0 - 0.8(0.832209)}{-1.405217} = 0.473782, \\ y_2 &= \frac{d_2 - c_2 y_3}{b_2} = \frac{0 - 0.8(0.473782)}{-2.64} = 0.205992, \\ y_1 &= \frac{d_1 - c_1 y_2}{b_1} = \frac{0 - 0(0.205992)}{1} = 0, \\ \end{align*} \end{split}\]

The solutions to this boundary value problem using the finite-difference method are tabulated below.

\(t_n\)	\(y_n\)	Exact	Error
0.0	0.000000	0.000000	0.00e+00
0.2	0.205992	0.221296	1.53e-02
0.4	0.473782	0.501419	2.76e-02
0.6	0.832209	0.865840	3.36e-02
0.8	1.321853	1.349412	2.76e-02
1.0	2.000000	2.000000	0.00e+00

---

5.3.2. Code

The code below calculates the solution to the boundary value problem in Example 5.4. The coefficient arrays a, b and c along with the constant vector d are defined for this problem and the tridiagonal_solver() function is used to solve the linear system. Note that the values in a, b, c and d are specific to the problem being solved and will need to be changed when solving other boundary value problems.

Python

import numpy as np
import matplotlib.pyplot as plt

def tridiagonal_solver(a, b, c, d):
    n = len(b)
    for i in range(1, n):
        w = a[i] / b[i-1]
        b[i] -= w * c[i-1]
        d[i] -= w * d[i-1]
    
    x = np.zeros(n)
    x[-1] = d[-1] / b[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (d[i] - c[i] * x[i+1]) / b[i]
    
    return x
    

def exact(t):
    return (2 * np.exp((1 - np.sqrt(5)) * (t - 1) / 2) * (np.exp(np.sqrt(5) * t) -1)) / (np.exp(np.sqrt(5)) - 1)
 

# Define BVP parameters
tspan = [0, 1]  # boundaries of the t domain
bvals = [0, 2]  # boundary values
h = 0.2         # step length

# Define linear system
n = int((tspan[1] - tspan[0]) / h) + 1
t = np.arange(n) * h
a, b, c, d = [np.zeros(n) for _ in range(4)]
b[::n-1] = 1
d[::n-1] = bvals
for i in range(1, n - 1):
    a[i] = 1
    b[i] = h - h ** 2 - 2
    c[i] = 1 - h
    d[i] = 0
        
# Solve linear system
y = tridiagonal_solver(a, b, c, d)

# Print table of solution values
print("|  t  |    y     |  Exact   |  Error   |")
print("|:---:|:--------:|:--------:|:--------:|")
for n in range(len(t)):
    print(f"| {t[n]:0.1f} | {y[n]:6.6f} | {exact(t[n]):6.6f} | {abs(exact(t[n]) - y[n]):6.2e} |")
    
# Plot solution
t_exact = np.linspace(0, tspan[1], 200)
y_exact = exact(t_exact)
fig, ax = plt.subplots()
plt.plot(t_exact, y_exact, "k-", label="Exact")
plt.plot(t, y, "bo-", label=f"Finite-difference method")
plt.xlabel("$t$", fontsize=14)
plt.ylabel("$y$", fontsize=14)
plt.legend(fontsize=12)
plt.show()

MATLAB

% Define ODE
f = @(t, y) [y(2), y(1) + y(2)];

% Define exact solution
exact = @(t, y) (2 .* exp((1 - sqrt(5)) .* (t - 1) / 2) .* (exp(sqrt(5) .* t) - 1)) / (exp(sqrt(5)) - 1);

% Solve BVP using finite-difference method
n = floor((tspan(2) - tspan(1)) / h) + 1;
t = (0 : n - 1) * h;
a = zeros(n, 1);
b = zeros(n, 1);
c = zeros(n, 1);
d = zeros(n, 1);
b([1, end]) = 1;
d([1, end]) = bvals;
for i = 2 : n - 1
    a(i) = 1;
    b(i) = h - h ^ 2 - 2;
    c(i) = 1 - h;
    d(i) = 0;
end
y = tridiagonal_solver(a, b, c, d);

% Calculate exact solution
t_exact = linspace(tspan(1), tspan(2), 200);
y_exact = exact(t_exact);

% Plot solutions
plot(t_exact, y_exact, "k", LineWidth=2)
hold on
plot(t, y, "b-o", LineWidth=2, MarkerFaceColor="b")
hold off
axis padded
xlabel("$t$", FontSize=14, Interpreter="latex")
ylabel("$y$", FontSize=14, Interpreter="latex")
legend("Exact", "Finite-difference method", Location="northwest")

%% ----------------------------------------------------------------------------
function x = tridiagonal_solver(a, b, c, d)

n = length(b);
for i = 2 : n
    w = a(i) / b(i-1);
    b(i) = b(i) - w * c(i-1);
    d(i) = d(i) - w * d(i-1);
end

x = zeros(n, 1);
x(n) = d(n) / b(n);
for i = n - 1 : -1 : 1
    x(i) = (d(i) - c(i) * x(i+1)) / b(i);
end

end

5.3.3. Order vs number of nodes

The solutions seen in Example 5.4 seem to show that the finite-difference method produces reasonably accurate results for this boundary value problem. One way to improve on the accuracy of our solution is to increase the number of nodes used.

Consider the solution of the following BVP using the finite-difference method

\[ y'' + 3ty' + 7y = \cos (2t), \qquad y(0) = 1, \qquad y(3) = 0. \]

Using forward and symmetric differences to approximate \(y'\) and \(y''\) respectively gives

\[\begin{split} \begin{align*} \frac{y_{i-1} - 2y_i + y_{i+1}}{h^2} + 3t_i \left( \frac{y_{i+1} - y_i}{h}\right) + 7y_i &= \cos(2t_i) \\ y_{i-1} + (7h^2 - 3ht_i - 2)y_i + (1 + 3ht_i)y_{i+1} &= h^2 \cos(2t_i). \end{align*} \end{split}\]

So the linear system is

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 \\ 1 & 7h^2 - 3ht_1 - 2 & 1 + 3ht_1 \\ & \ddots & \ddots & \ddots \\ && 1 & 7h^2 - 3ht_{N-1} - 2 & 1 + 3ht_{N-1} \\ &&& 0 & 1 \end{pmatrix} \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_{N-1} \\ y_N \end{pmatrix} \\ = \begin{pmatrix} 1 \\ h^2 \cos(2t_1) \\ \vdots \\ h^2 \cos(2t_{N-1}) \\ 0 \end{pmatrix}. \end{align*} \end{split}\]

The solution of this BVP is shown in bvp-fdm-accuracy-figure-1 for \(h=0.05\) and \(h = 0.005\). Since we used a first-order approximation for \(y'\) the error for this method is \(O(h)\) and we expect the solution using \(h=0.005\) to be more accurate.

To obtain a more accurate solution, instead of increasing the number of nodes we could use the central difference approximation (which is second-order accurate) to approximate \(y'\)

\[\begin{split} \begin{align*} \frac{y_{i-1} - 2y_i + y_{i+1}}{h^2} + 3t_i\left( \frac{y_{i+1} - y_{i-1}}{2h} \right) + 7y_i &= \cos(2t_i) \\ (2 - 3ht_i)y_{i-1} + (14h^2 - 4)y_i + (2 + 3ht_i)y_{i+1} &= 2h^2 \cos(2t_i). \end{align*} \end{split}\]

So the linear system for the second-order finite-difference method is

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 \\ 2 - 3ht_1 & 14h^2 - 4 & 2 + 3ht_1 \\ & \ddots & \ddots & \ddots \\ && 2 - 3ht_{N-1} & 14h^2 - 4 & 2 + 3ht_{N-1} \\ &&& 0 & 1 \end{pmatrix} \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_{N-1} \\ y_N \end{pmatrix} \\ = \begin{pmatrix} 1 \\ 2h^2 \cos(2t_1) \\ \vdots \\ 2h^2 \cos(2t_{N-1}) \\ 0 \end{pmatrix}. \end{align*} \end{split}\]

The solution using the second-order finite difference method with \(h=0.05\) has been plotted against the first-order solution using \(h=0.05\) and \(h=0.005\) in bvp-fdm-accuracy-figure-2. The second order solution gives good agreement with the first order solution using 10 times fewer nodes.