6: Finite Difference Approximation

Last updated
Save as PDF

Page ID: 93752

$ \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } $ $ \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} $$\newcommand{\id}{\mathrm{id}}$ $ \newcommand{\Span}{\mathrm{span}}$ $ \newcommand{\kernel}{\mathrm{null}\,}$ $ \newcommand{\range}{\mathrm{range}\,}$ $ \newcommand{\RealPart}{\mathrm{Re}}$ $ \newcommand{\ImaginaryPart}{\mathrm{Im}}$ $ \newcommand{\Argument}{\mathrm{Arg}}$ $ \newcommand{\norm}[1]{\| #1 \|}$ $ \newcommand{\inner}[2]{\langle #1, #2 \rangle}$ $ \newcommand{\Span}{\mathrm{span}}$ $\newcommand{\id}{\mathrm{id}}$ $ \newcommand{\Span}{\mathrm{span}}$ $ \newcommand{\kernel}{\mathrm{null}\,}$ $ \newcommand{\range}{\mathrm{range}\,}$ $ \newcommand{\RealPart}{\mathrm{Re}}$ $ \newcommand{\ImaginaryPart}{\mathrm{Im}}$ $ \newcommand{\Argument}{\mathrm{Arg}}$ $ \newcommand{\norm}[1]{\| #1 \|}$ $ \newcommand{\inner}[2]{\langle #1, #2 \rangle}$ $ \newcommand{\Span}{\mathrm{span}}$$\newcommand{\AA}{\unicode[.8,0]{x212B}}$

We introduce here numerical differentiation, also called finite difference approximation. This technique is commonly used to discretize and solve partial differential equations.

Finite difference formulas

Consider the Taylor series approximation for $y(x+h)$ and $y(x-h)$, given by

\[\begin{aligned} &y(x+h)=y(x)+h y^{\prime}(x)+\frac{1}{2} h^{2} y^{\prime \prime}(x)+\frac{1}{6} h^{3} y^{\prime \prime \prime}(x)+\frac{1}{24} h^{4} y^{\prime \prime \prime \prime}(x)+\ldots, \\ &y(x-h)=y(x)-h y^{\prime}(x)+\frac{1}{2} h^{2} y^{\prime \prime}(x)-\frac{1}{6} h^{3} y^{\prime \prime \prime}(x)+\frac{1}{24} h^{4} y^{\prime \prime \prime \prime}(x)+\ldots . \end{aligned} \nonumber \]

The standard definitions of the derivatives give the first-order approximations

\[\begin{aligned} &y^{\prime}(x)=\frac{y(x+h)-y(x)}{h}+\mathrm{O}(h), \\ &y^{\prime}(x)=\frac{y(x)-y(x-h)}{h}+\mathrm{O}(h) . \end{aligned} \nonumber \]

The more widely-used second-order approximation is called the central-difference approximation and is given by

\[y^{\prime}(x)=\frac{y(x+h)-y(x-h)}{2 h}+\mathrm{O}\left(h^{2}\right) . \nonumber \]

The finite difference approximation to the second derivative can be found from considering

\[y(x+h)+y(x-h)=2 y(x)+h^{2} y^{\prime \prime}(x)+\frac{1}{12} h^{4} y^{\prime \prime \prime \prime}(x)+\ldots, \nonumber \]

from which we find

\[y^{\prime \prime}(x)=\frac{y(x+h)-2 y(x)+y(x-h)}{h^{2}}+\mathrm{O}\left(h^{2}\right) . \nonumber \]

Often a second-order method is required for $x$ on the boundaries of the domain. For a boundary point on the left, a second-order forward difference method requires the additional Taylor series

\[y(x+2 h)=y(x)+2 h y^{\prime}(x)+2 h^{2} y^{\prime \prime}(x)+\frac{4}{3} h^{3} y^{\prime \prime \prime}(x)+\ldots \nonumber \]

We combine the Taylor series for $y(x+h)$ and $y(x+2 h)$ to eliminate the term proportional to $h^{2}$ :

\[y(x+2 h)-4 y(x+h)=-3 y(x)-2 h y^{\prime}(x)+\mathrm{O}\left(\mathrm{h}^{3}\right) \text {. } \nonumber \]

Therefore,

\[y^{\prime}(x)=\frac{-3 y(x)+4 y(x+h)-y(x+2 h)}{2 h}+\mathrm{O}\left(h^{2}\right) . \nonumber \]

For a boundary point on the right, we send $h \rightarrow-h$ to find

\[y^{\prime}(x)=\frac{3 y(x)-4 y(x-h)+y(x-2 h)}{2 h}+\mathrm{O}\left(h^{2}\right) . \nonumber \]

Example: the Laplace equation

As an example of the finite difference technique, let us consider how to discretize the two dimensional Laplace equation

\[\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right) \Phi=0 \nonumber \]

on the rectangular domain $[0,1] \times[0,1]$. We assume here that the values of $\Phi$ are known on the boundaries of the domain. We form a two-dimensional grid with $N$ intervals in each direction together with end points that lie on the boundaries by writing

\[\begin{array}{ll} x_{i}=i h, & i=0,1, \ldots, N, \\ y_{j}=j h, & j=0,1, \ldots, N, \end{array} \nonumber \]

where $h=1 / N$. We will also denote the value of $\Phi\left(x_{i}, y_{j}\right)$ by $\Phi_{i, j}$. The finite difference approximation for the second derivatives at the interior point $\left(x_{i}, y_{j}\right)$ then results in an equation that we write in the form

\[4 \Phi_{i, j}-\Phi_{i+1, j}-\Phi_{i-1, j}-\Phi_{i, j+1}-\Phi_{i, j-1}=0, \nonumber \]

valid for $i=1,2, \ldots, N-1$ and $j=1,2, \ldots, N-1$. The boundary values $\Phi_{0, j}, \Phi_{N, j}, \Phi_{i, 0}$, and $\Phi_{i, N}$ are assumed to be given.

One can observe that (6.10) represents a system of $(N-1)^{2}$ linear equations for $(N-1)^{2}$ unknowns. We can write this as a matrix equation if we decide how to order the unknowns $\Phi_{i, j}$ into a vector. The standard ordering is called natural ordering, and proceeds from the bottom of the grid along rows moving towards the top. This orders the unknowns as

\[\Phi=\left[\Phi_{1,1}, \Phi_{2,1}, \ldots, \Phi_{(N-1), 1}, \ldots, \Phi_{1,(N-1)}, \Phi_{2(N-1)}, \ldots, \Phi_{(N-1),(N-1)}\right]^{T} . \nonumber \]

To illustrate the construction of the matrix equation, we consider the case $N=3$, with two interior points in each direction. The four resulting linear equations with the boundary terms written on the right-hand-side are

\[\begin{aligned} &4 \Phi_{1,1}-\Phi_{2,1}-\Phi_{1,2}=\Phi_{0,1}+\Phi_{1,0} \\ &4 \Phi_{2,1}-\Phi_{1,1}-\Phi_{2,2}=\Phi_{3,1}+\Phi_{2,0} \\ &4 \Phi_{1,2}-\Phi_{2,2}-\Phi_{1,1}=\Phi_{0,2}+\Phi_{1,3} \\ &4 \Phi_{2,2}-\Phi_{1,2}-\Phi_{2,1}=\Phi_{3,2}+\Phi_{2,3} \end{aligned} \nonumber \]

and the corresponding matrix equation is

\[\left(\begin{array}{rrrr} 4 & -1 & -1 & 0 \\ -1 & 4 & 0 & -1 \\ -1 & 0 & 4 & -1 \\ 0 & -1 & -1 & 4 \end{array}\right)\left(\begin{array}{l} \Phi_{1,1} \\ \Phi_{2,1} \\ \Phi_{1,2} \\ \Phi_{2,2} \end{array}\right)=\left(\begin{array}{l} \Phi_{0,1}+\Phi_{1,0} \\ \Phi_{3,1}+\Phi_{2,0} \\ \Phi_{0,2}+\Phi_{1,3} \\ \Phi_{3,2}+\Phi_{2,3} \end{array}\right) . \nonumber \]

The pattern here may not be obvious, but the Laplacian matrix decomposes into 2 -by- 2 block matrices.

With some thought, a general result can be obtained for an $N_{x}$-by- $N_{y}$ grid of internal points (with $N_{x}+1$ and $N_{y}+1$ intervals in the $x$ - and y-directions, respectively). The Laplacian matrix then decomposes into $N_{x}$-by- $N_{x}$ block matrices. The diagonal contains $N_{y}$ of these $N_{x}$-by- $N_{x}$ block matrices, each of which are tridiagonal with a 4 on the diagonal and a $-1$ on the off-diagonals. Immediately above and below these block matrices on the diagonal are $N_{y}-1$ block matrices also of size $N_{x}$-by- $N_{x}$, each of which are diagonal with $-1$ along the diagonal.

MATLAB code for the Laplacian matrix can be found on the web in the function sp_laplace.m. This code was written by Bordner and Saied in 1995, and I have written a more modern and faster version of this code in sp_laplace_new.m.

An alternative solution method, which we will later make use of in $\$ 17.2 .4$, includes the boundary values in the solution vector. If one of the boundary conditions is at position $n$ in this vector, then row $n$ of the left-hand-side matrix will have just a one on the diagonal with all other elements equal to zero (i.e, a row of the identity matrix), and the corresponding element in the right-hand-side vector will have the value at the boundary. Here, with the total number of grid points (including the boundary points) in the $x$ - and $y$ directions given by $N_{x}$ and $N_{y}$, the left-hand-side matrix is then generated using $A=$ sp_laplace_new (N_X, N_Y), and all the rows corresponding to the boundary values are replaced with the corresponding rows of the identity matrix. This formulation may be slightly easier to code, and also easier to incorporate other more general boundary conditions.