7.3: Numerical Methods - Boundary Value Problem

7.3.1. Finite difference method

We consider first the differential equation

\[-\frac{d^{2} y}{d x^{2}}=f(x), \quad 0 \leq x \leq 1 \nonumber \]

with two-point boundary conditions

\[y(0)=A, \quad y(1)=B \text {. } \nonumber \]

Equation (7.8) can be solved by quadrature, but here we will demonstrate a numerical solution using a finite difference method.

We begin by discussing how to numerically approximate derivatives. 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) \text {. } \nonumber \]

Sometimes 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) \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 \]

We now write a finite difference scheme to solve (7.8). We discretize \(x\) by defining \(x_{i}=i h, i=0,1, \ldots, n+1\). Since \(x_{n+1}=1\), we have \(h=1 /(n+1)\). The functions \(y(x)\) and \(f(x)\) are discretized as \(y_{i}=y\left(x_{i}\right)\) and \(f_{i}=f\left(x_{i}\right) .\) The second-order differential equation (7.8) then becomes for the interior points of the domain

\[-y_{i-1}+2 y_{i}-y_{i+1}=h^{2} f_{i}, \quad i=1,2, \ldots n, \nonumber \]

with the boundary conditions \(y_{0}=A\) and \(y_{n+1}=B\). We therefore have a linear system of equations to solve. The first and \(n\)th equation contain the boundary conditions and are given by

\[\begin{aligned} 2 y_{1}-y_{2} &=h^{2} f_{1}+A \\ -y_{n-1}+2 y_{n} &=h^{2} f_{n}+B \end{aligned} \nonumber \]

The second and third equations, etc., are

\[\begin{aligned} &-y_{1}+2 y_{2}-y_{3}=h^{2} f_{2} \\ &-y_{2}+2 y_{3}-y_{4}=h^{2} f_{3} \end{aligned} \nonumber \]

In matrix form, we have

\[\left(\begin{array}{rrrrrrrr} 2 & -1 & 0 & 0 & \cdots & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & -1 & 2 & -1 \\ 0 & 0 & 0 & 0 & \cdots & 0 & -1 & 2 \end{array}\right)\left(\begin{array}{c} y_{1} \\ y_{2} \\ y_{3} \\ \vdots \\ y_{n-1} \\ y_{n} \end{array}\right)=\left(\begin{array}{c} h^{2} f_{1}+A \\ h^{2} f_{2} \\ h^{2} f_{3} \\ \vdots \\ h^{2} f_{n-1} \\ h^{2} f_{n}+B \end{array}\right) \nonumber \]

The matrix is tridiagonal, and a numerical solution by Guassian elimination can be done quickly. The matrix itself is easily constructed using the MATLAB function diag.m and ones.m. As excerpted from the MATLAB help page, the function call ones \((\mathrm{m}, \mathrm{n})\) returns an m-by-n matrix of ones, and the function call \(\operatorname{diag}(\mathrm{v}, \mathrm{k})\), when \(\mathrm{v}\) is a vector with \(\mathrm{n}\) components, is a square matrix of order \(\mathrm{n}+\mathrm{abs}(\mathrm{k})\) with the elements of \(\mathrm{v}\) on the \(\mathrm{k}\)-th diagonal: \(k=0\) is the main diagonal, \(k>0\) is above the main diagonal and \(k<0\) is below the main diagonal. The \(n \times n\) matrix above can be constructed by the MATLAB code

\[M=\operatorname{diag}(-\text { ones }(n-1,1),-1)+\operatorname{diag}\left(2^{*} \text { ones }(n, 1), 0\right)+\operatorname{diag}(-\text { ones }(n-1,1), 1) ; \nonumber \]

The right-hand-side, provided \(\mathrm{f}\) is a given n-by-1 vector, can be constructed by the MATLAB code

\[\mathrm{b}=\mathrm{h}^{\wedge} 2^{*} \mathrm{f} ; \mathrm{b}(1)=\mathrm{b}(1)+\mathrm{A} ; \mathrm{b}(\mathrm{n})=\mathrm{b}(\mathrm{n})+\mathrm{B} \nonumber \]

and the solution for \(u\) is given by the MATLAB code

\[\mathrm{y}=\mathrm{M} \backslash \mathrm{b} \nonumber \]

7.3.2. Shooting method

The finite difference method can solve linear odes. For a general ode of the form

\[\frac{d^{2} y}{d x^{2}}=f(x, y, d y / d x) \nonumber \]

with \(y(0)=A\) and \(y(1)=B\), we use a shooting method. First, we formulate the ode as an initial value problem. We have

\[\begin{aligned} &\frac{d y}{d x}=z \\ &\frac{d z}{d x}=f(x, y, z) \end{aligned} \nonumber \]

The initial condition \(y(0)=A\) is known, but the second initial condition \(z(0)=b\) is unknown. Our goal is to determine \(b\) such that \(y(1)=B\).

In fact, this is a root-finding problem for an appropriately defined function. We define the function \(F=F(b)\) such that

\[F(b)=y(1)-B \nonumber \]

In other words, \(F(b)\) is the difference between the value of \(y(1)\) obtained from integrating the differential equations using the initial condition \(z(0)=b\), and \(B\). Our root-finding routine will want to solve \(F(b)=0\). (The method is called shooting because the slope of the solution curve for \(y=y(x)\) at \(x=0\) is given by \(b\), and the solution hits the value \(y(1)\) at \(x=1\). This looks like pointing a gun and trying to shoot the target, which is \(B\).)

To determine the value of \(b\) that solves \(F(b)=0\), we iterate using the Secant method, given by

\[b_{n+1}=b_{n}-F\left(b_{n}\right) \frac{b_{n}-b_{n-1}}{F\left(b_{n}\right)-F\left(b_{n-1}\right)} \nonumber \]

We need to start with two initial guesses for \(b\), solving the ode for the two corresponding values of \(y(1)\). Then the Secant Method will give us the next value of \(b\) to try, and we iterate until \(|y(1)-B|<\) tol, where tol is some specified tolerance for the error