Boundary Value Problems by Finite Differences

A boundary value problem fixes $y$ at both ends of an interval instead of giving $y$ and $y'$ at the start, so stepping forward from one end does not work directly. Instead, put a grid on the interval, replace $y''$ by the second difference, and solve for all the unknown grid values at once with a linear system, using Solving with PA = LU.

For $y'' = f(x)$ on $[a, b]$ with $y(a) = \alpha$, $y(b) = \beta$, take $x_i = a + ih$, $i = 0, \ldots, N$, with $h = \frac{b - a}{N}$. The unknowns are $y_1, \ldots, y_{N-1}$; $y_0 = \alpha$ and $y_N = \beta$ are known. At each inside point, $$y_{i-1} - 2y_i + y_{i+1} = h^{2}f(x_i).$$

Key idea

The equations form a tridiagonal system: $1, -2, 1$ around the diagonal. In the first and last equations the known boundary values move to the right-hand side. The error at the grid points is $O(h^{2})$, and the system is solved in $O(N)$ operations.

Why it works

The second difference has truncation error $\frac{h^{2}}{12}y^{(4)}$, so the discrete equations differ from the true ones by $O(h^{2})$; for this kind of problem that carries over to the computed values. Each equation involves only three neighbors, so the matrix is tridiagonal. It is also symmetric and diagonally dominant (weakly), so elimination needs no pivoting and costs about $8N$ operations instead of $\frac{2}{3}N^{3}$.

Worked example: A problem the method solves exactly

Solve $y'' = 6x$, $y(0) = 0$, $y(1) = 1$ with $h = \frac{1}{3}$.

Inside points $\frac{1}{3}, \frac{2}{3}$; $h^{2} = \frac{1}{9}$.
$i = 1$: $y_0 - 2y_1 + y_2 = \frac{1}{9} \cdot 2$, and $y_0 = 0$: $-2y_1 + y_2 = \frac{2}{9}$.
$i = 2$: $y_1 - 2y_2 + y_3 = \frac{1}{9} \cdot 4$, and $y_3 = 1$ moves right: $y_1 - 2y_2 = \frac{4}{9} - 1 = -\frac{5}{9}$.
Solve: $y_1 = \frac{1}{27}$, $y_2 = \frac{8}{27}$. The exact solution is $x^{3}$, and these are exact, because the second difference has no error for cubics.

Answer: $y_1 = \frac{1}{27}$, $y_2 = \frac{8}{27}$

Worked example: A term in y

Set up and solve $y'' - 4y = 0$, $y(0) = 0$, $y(1) = 1$ with $h = \frac{1}{4}$.

$\frac{y_{i-1} - 2y_i + y_{i+1}}{h^{2}} - 4y_i = 0$ becomes $y_{i-1} - 2.25y_i + y_{i+1} = 0$.
System: $\begin{bmatrix} -2.25 & 1 & 0 \\ 1 & -2.25 & 1 \\ 0 & 1 & -2.25 \end{bmatrix}\begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ -1 \end{bmatrix}$ (the $y_4 = 1$ moved to the right).
Solution: $(0.1451, 0.3265, 0.5896)$. Exact: $\frac{\sinh 2x}{\sinh 2}$ gives $(0.1437, 0.3240, 0.5871)$: errors about $0.002$, and halving $h$ would cut them by about 4.

Answer: $\approx (0.1451, 0.3265, 0.5896)$

The finite-difference solution of y″ − 4y = 0, y(0) = 0, y(1) = 1. Red squares are the known boundary values; blue dots are the unknowns that the tridiagonal system finds.

The same recipe handles $y'' + p(x)y' + q(x)y = f(x)$: replace $y'$ by the central difference $\frac{y_{i+1} - y_{i-1}}{2h}$ as well, and the matrix stays tridiagonal. A boundary condition on $y'$ instead of $y$ adds an equation (or a "ghost point") at that end. Very fine grids are cheap thanks to the tridiagonal structure; the limit is round-off, which grows like $1/h^{2}$ as in The Best Step Size.

Common mistake

Treating the known boundary values as unknowns, or leaving them on the left side. $y_0$ and $y_N$ are numbers, not unknowns; the system has $N - 1$ equations in $y_1, \ldots, y_{N-1}$, and $\alpha$ and $\beta$ appear only on the right side of the first and last equations. Check: plug your answer back in, including the boundary values, and every equation should hold.

Sources and further reading

Tobin A. Driscoll and Richard J. Braun, Fundamentals of Numerical Computation (University of Delaware, free online edition): chapter "Boundary-value problems", sections on finite differences and collocation for linear problems

Autar Kaw, Holistic Numerical Methods (University of South Florida): Ordinary Differential Equations, "Finite Difference Method"

Boundary Value Problems by Finite Differences

Practice

Builds on

Leads to

Sources and further reading