smallest positive integer that is not exact in single pre- cision

Let $N$ be the smallest positive integer that is not exact. Now, I claim that

$$N-2=2^{23} \times 1.11 \ldots 1 \nonumber$$

and

$$N-1=2^{24} \times 1.00 \ldots 0 \nonumber$$

The integer $N$ would then require a one-bit in the $2^{-24}$ position, which is not available. Therefore, the smallest positive integer that is not exact is $2^{24}+1=16777217$. In MATLAB, single $\left(2^{24}\right)$ has the same value as single $\left(2^{24}+1\right) .$ Since single $\left(2^{24}+1\right)$ is exactly halfway between the two consecutive machine numbers $2^{24}$ and $2^{24}+2$, MATLAB rounds to the number with a final zero-bit in $f$, which is $2^{24} .$

Estimate $\sqrt{2}=1.41421356$ using Newton’s Method

The $\sqrt{2}$ is the zero of the function $f(x)=x^{2}-2$. To implement Newton’s Method, we use $f^{\prime}(x)=2 x$. Therefore, Newton’s Method is the iteration

$$x_{n+1}=x_{n}-\frac{x_{n}^{2}-2}{2 x_{n}} \nonumber$$

We take as our initial guess $x_{0}=1$. Then

$$\begin{aligned} &x_{1}=1-\frac{-1}{2}=\frac{3}{2}=1.5 \\ &x_{2}=\frac{3}{2}-\frac{\frac{9}{4}-2}{3}=\frac{17}{12}=1.416667, \\ &x_{3}=\frac{17}{12}-\frac{\frac{17^{2}}{12^{2}}-2}{\frac{17}{6}}=\frac{577}{408}=1.41426 . \end{aligned} \nonumber$$

Example of fractals using Newton’s Method

Consider the complex roots of the equation $f(z)=0$, where

$$f(z)=z^{3}-1 \nonumber$$

These roots are the three cubic roots of unity. With

$$e^{i 2 \pi n}=1, \quad n=0,1,2, \ldots \nonumber$$

the three unique cubic roots of unity are given by

$$1, \quad e^{i 2 \pi / 3}, \quad e^{i 4 \pi / 3} \nonumber$$

With

$$e^{i \theta}=\cos \theta+i \sin \theta, \nonumber$$

and $\cos (2 \pi / 3)=-1 / 2, \sin (2 \pi / 3)=\sqrt{3} / 2$, the three cubic roots of unity are

$$r_{1}=1, \quad r_{2}=-\frac{1}{2}+\frac{\sqrt{3}}{2} i, \quad r_{3}=-\frac{1}{2}-\frac{\sqrt{3}}{2} i \nonumber$$

The interesting idea here is to determine which initial values of $z_{0}$ in the complex plane converge to which of the three cubic roots of unity.

Newton’s method is

$$z_{n+1}=z_{n}-\frac{z_{n}^{3}-1}{3 z_{n}^{2}} \nonumber$$

If the iteration converges to $r_{1}$, we color $z_{0}$ red; $r_{2}$, blue; $r_{3}$, green. The result will be shown in lecture.

Newton’s Method

We start with Newton’s Method

$$x_{n+1}=x_{n}-\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)} \nonumber$$

Subtracting both sides from $r$, we have

$$r-x_{n+1}=r-x_{n}+\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)} \nonumber$$

Or

$$\epsilon_{n+1}=\epsilon_{n}+\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)} \nonumber$$

We use Taylor series to expand the functions $f\left(x_{n}\right)$ and $f^{\prime}\left(x_{n}\right)$ about the root $r$, using $f(r)=0$. We have

$$\begin{aligned} f\left(x_{n}\right) &=f(r)+\left(x_{n}-r\right) f^{\prime}(r)+\frac{1}{2}\left(x_{n}-r\right)^{2} f^{\prime \prime}(r)+\ldots, \\ &=-\epsilon_{n} f^{\prime}(r)+\frac{1}{2} \epsilon_{n}^{2} f^{\prime \prime}(r)+\ldots ; \\ f^{\prime}\left(x_{n}\right) &=f^{\prime}(r)+\left(x_{n}-r\right) f^{\prime \prime}(r)+\frac{1}{2}\left(x_{n}-r\right)^{2} f^{\prime \prime \prime}(r)+\ldots, \\ &=f^{\prime}(r)-\epsilon_{n} f^{\prime \prime}(r)+\frac{1}{2} \epsilon_{n}^{2} f^{\prime \prime \prime}(r)+\ldots \end{aligned} \nonumber$$

To make further progress, we will make use of the following standard Taylor series:

$$\frac{1}{1-\epsilon}=1+\epsilon+\epsilon^{2}+\ldots, \nonumber$$

which converges for $|\epsilon|<1 .$ Substituting $(2.2)$ into $(2.1)$, and using $(2.3)$ yields

$$\begin{aligned} \epsilon_{n+1} &=\epsilon_{n}+\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)} \\ &=\epsilon_{n}+\frac{-\epsilon_{n} f^{\prime}(r)+\frac{1}{2} \epsilon_{n}^{2} f^{\prime \prime}(r)+\ldots}{f^{\prime}(r)-\epsilon_{n} f^{\prime \prime}(r)+\frac{1}{2} \epsilon_{n}^{2} f^{\prime \prime \prime}(r)+\ldots} \\ &=\epsilon_{n}+\frac{-\epsilon_{n}+\frac{1}{2} \epsilon_{n}^{2} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}+\ldots}{1-\epsilon_{n} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}+\ldots} \\ &=\epsilon_{n}+\left(-\epsilon_{n}+\frac{1}{2} \epsilon_{n}^{2} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}+\ldots\right)\left(1+\epsilon_{n} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}+\ldots\right) \\ &=\epsilon_{n}+\left(-\epsilon_{n}+\epsilon_{n}^{2}\left(\frac{1}{2} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}-\frac{f^{\prime \prime}(r)}{f^{\prime}(r)}\right)+\ldots\right) \\ &=-\frac{1}{2} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)} \epsilon_{n}^{2}+\ldots \end{aligned} \nonumber$$

Therefore, we have shown that

$$\left|\epsilon_{n+1}\right|=k\left|\epsilon_{n}\right|^{2} \nonumber$$

as $n \rightarrow \infty$, with

$$k=\frac{1}{2}\left|\frac{f^{\prime \prime}(r)}{f^{\prime}(r)}\right| \nonumber$$

provided $f^{\prime}(r) \neq 0 .$ Newton’s method is thus of order 2 at simple roots.

Secant Method

Determining the order of the Secant Method proceeds in a similar fashion. We start with

$$x_{n+1}=x_{n}-\frac{\left(x_{n}-x_{n-1}\right) f\left(x_{n}\right)}{f\left(x_{n}\right)-f\left(x_{n-1}\right)} \nonumber$$

We subtract both sides from $r$ and make use of

$$\begin{aligned} x_{n}-x_{n-1} &=\left(r-x_{n-1}\right)-\left(r-x_{n}\right) \\ &=\epsilon_{n-1}-\epsilon_{n} \end{aligned} \nonumber$$

and the Taylor series

$$\begin{aligned} f\left(x_{n}\right) &=-\epsilon_{n} f^{\prime}(r)+\frac{1}{2} \epsilon_{n}^{2} f^{\prime \prime}(r)+\ldots, \\ f\left(x_{n-1}\right) &=-\epsilon_{n-1} f^{\prime}(r)+\frac{1}{2} \epsilon_{n-1}^{2} f^{\prime \prime}(r)+\ldots, \end{aligned} \nonumber$$

so that

$$\begin{aligned} f\left(x_{n}\right)-f\left(x_{n-1}\right) &=\left(\epsilon_{n-1}-\epsilon_{n}\right) f^{\prime}(r)+\frac{1}{2}\left(\epsilon_{n}^{2}-\epsilon_{n-1}^{2}\right) f^{\prime \prime}(r)+\ldots \\ &=\left(\epsilon_{n-1}-\epsilon_{n}\right)\left(f^{\prime}(r)-\frac{1}{2}\left(\epsilon_{n-1}+\epsilon_{n}\right) f^{\prime \prime}(r)+\ldots\right) \end{aligned} \nonumber$$

We therefore have

$$\begin{aligned} \epsilon_{n+1} &=\epsilon_{n}+\frac{-\epsilon_{n} f^{\prime}(r)+\frac{1}{2} \epsilon_{n}^{2} f^{\prime \prime}(r)+\ldots}{f^{\prime}(r)-\frac{1}{2}\left(\epsilon_{n-1}+\epsilon_{n}\right) f^{\prime \prime}(r)+\ldots} \\ &=\epsilon_{n}-\epsilon_{n} \frac{1-\frac{1}{2} \epsilon_{n} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}+\ldots}{1-\frac{1}{2}\left(\epsilon_{n-1}+\epsilon_{n}\right) \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}+\ldots} \\ &=\epsilon_{n}-\epsilon_{n}\left(1-\frac{1}{2} \epsilon_{n} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}+\ldots\right)\left(1+\frac{1}{2}\left(\epsilon_{n-1}+\epsilon_{n}\right) \frac{f^{\prime \prime}(r)}{f^{\prime}(r)}+\ldots\right) \\ &=-\frac{1}{2} \frac{f^{\prime \prime}(r)}{f^{\prime}(r)} \epsilon_{n-1} \epsilon_{n}+\ldots, \end{aligned} \nonumber$$

or to leading order

$$\left|\epsilon_{n+1}\right|=\frac{1}{2}\left|\frac{f^{\prime \prime}(r)}{f^{\prime}(r)}\right|\left|\epsilon_{n-1}\right|\left|\epsilon_{n}\right| \nonumber$$

The order of convergence is not yet obvious from this equation, and to determine the scaling law we look for a solution of the form

$$\left|\epsilon_{n+1}\right|=k\left|\epsilon_{n}\right|^{p} . \nonumber$$

From this ansatz, we also have

$$\left|\epsilon_{n}\right|=k\left|\epsilon_{n-1}\right|^{p} \nonumber$$

and therefore

$$\left|\epsilon_{n+1}\right|=k^{p+1}\left|\epsilon_{n-1}\right|^{p^{2}} \nonumber$$

Substitution into $(2.4)$ results in

$$k^{p+1}\left|\epsilon_{n-1}\right|^{p^{2}}=\frac{k}{2}\left|\frac{f^{\prime \prime}(r)}{f^{\prime}(r)}\right|\left|\epsilon_{n-1}\right|^{p+1} \nonumber$$

Equating the coefficient and the power of $\epsilon_{n-1}$ results in

$$k^{p}=\frac{1}{2}\left|\frac{f^{\prime \prime}(r)}{f^{\prime}(r)}\right| \nonumber$$

and

$$p^{2}=p+1 \nonumber$$

The order of convergence of the Secant Method, given by $p$, therefore is determined to be the positive root of the quadratic equation $p^{2}-p-1=0$, or

$$p=\frac{1+\sqrt{5}}{2} \approx 1.618 \nonumber$$

which coincidentally is a famous irrational number that is called The Golden Ratio, and goes by the symbol $\Phi$. We see that the Secant Method has an order of convergence lying between the Bisection Method and Newton’s Method.

Vandermonde polynomial

This Vandermonde polynomial is the most straightforward construction of the interpolating polynomial $P_{n}(x)$. We simply write

$$P_{n}(x)=c_{0} x^{n}+c_{1} x^{n-1}+\cdots+c_{n} . \nonumber$$

Then we can immediately form $n+1$ linear equations for the $n+1$ unknown coefficients $c_{0}, c_{1}, \ldots, c_{n}$ using the $n+1$ known points:

$$\begin{gathered} y_{0}=c_{0} x_{0}^{n}+c_{1} x_{0}^{n-1}+\cdots+c_{n-1} x_{0}+c_{n} \\ y_{2}=c_{0} x_{1}^{n}+c_{1} x_{1}^{n-1}+\cdots+c_{n-1} x_{1}+c_{n} \\ \vdots \\ y_{n}=c_{0} x_{n}^{n}+c_{1} x_{n}^{n-1}+\cdots+c_{n-1} x_{n}+c_{n} \end{gathered} \nonumber$$

The system of equations in matrix form is

$$\left(\begin{array}{ccccc} x_{0}^{n} & x_{0}^{n-1} & \cdots & x_{0} & 1 \\ x_{1}^{n} & x_{1}^{n-1} & \cdots & x_{1} & 1 \\ \vdots & \vdots & \ddots & \vdots & \\ x_{n}^{n} & x_{n}^{n-1} & \cdots & x_{n} & 1 \end{array}\right)\left(\begin{array}{c} c_{0} \\ c_{1} \\ \vdots \\ c_{n} \end{array}\right)=\left(\begin{array}{c} y_{0} \\ y_{1} \\ \vdots \\ y_{n} \end{array}\right) \nonumber$$

The matrix is called the Vandermonde matrix, and can be constructed using the MATLAB function vander.m. The system of linear equations can be solved in MATLAB using the \operator, and the MATLAB function polyval.m can used to interpolate using the $c$ coefficients. I will illustrate this in class and place the code on the website.

Lagrange polynomial

The Lagrange polynomial is the most clever construction of the interpolating polynomial $P_{n}(x)$, and leads directly to an analytical formula. The Lagrange polynomial is the sum of $n+1$ terms and each term is itself a polynomial of degree $n$. The full polynomial is therefore of degree $n$. Counting from 0 , the $i$ th term of the Lagrange polynomial is constructed by requiring it to be zero at $x_{j}$ with $j \neq i$, and equal to $y_{i}$ when $j=i$. The polynomial can be written as

$$\begin{array}{r} P_{n}(x)=\frac{\left(x-x_{1}\right)\left(x-x_{2}\right) \cdots\left(x-x_{n}\right) y_{0}}{\left(x_{0}-x_{1}\right)\left(x_{0}-x_{2}\right) \cdots\left(x_{0}-x_{n}\right)}+\frac{\left(x-x_{0}\right)\left(x-x_{2}\right) \cdots\left(x-x_{n}\right) y_{1}}{\left(x_{1}-x_{0}\right)\left(x_{1}-x_{2}\right) \cdots\left(x_{1}-x_{n}\right)} \\ +\cdots+\frac{\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n-1}\right) y_{n}}{\left(x_{n}-x_{0}\right)\left(x_{n}-x_{1}\right) \cdots\left(x_{n}-x_{n-1}\right)} . \end{array} \nonumber$$

It can be clearly seen that the first term is equal to zero when $x=x_{1}, x_{2}, \ldots, x_{n}$ and equal to $y_{0}$ when $x=x_{0}$; the second term is equal to zero when $x=x_{0}, x_{2}, \ldots x_{n}$ and equal to $y_{1}$ when $x=x_{1}$; and the last term is equal to zero when $x=x_{0}, x_{1}, \ldots x_{n-1}$ and equal to $y_{n}$ when $x=x_{n}$. The uniqueness of the interpolating polynomial implies that the Lagrange polynomial must be the interpolating polynomial.

Newton polynomial

The Newton polynomial is somewhat more clever than the Vandermonde polynomial because it results in a system of linear equations that is lower triangular, and therefore can be solved by forward substitution. The interpolating polynomial is written in the form

$$P_{n}(x)=c_{0}+c_{1}\left(x-x_{0}\right)+c_{2}\left(x-x_{0}\right)\left(x-x_{1}\right)+\cdots+c_{n}\left(x-x_{0}\right) \cdots\left(x-x_{n-1}\right) \nonumber$$

which is clearly a polynomial of degree $n$. The $n+1$ unknown coefficients given by the $c^{\prime}$ s can be found by substituting the points $\left(x_{i}, y_{i}\right)$ for $i=0, \ldots, n$ :

$$\begin{aligned} y_{0} &=c_{0} \\ y_{1} &=c_{0}+c_{1}\left(x_{1}-x_{0}\right) \\ y_{2} &=c_{0}+c_{1}\left(x_{2}-x_{0}\right)+c_{2}\left(x_{2}-x_{0}\right)\left(x_{2}-x_{1}\right) \\ \vdots & \vdots \\ y_{n} &=c_{0}+c_{1}\left(x_{n}-x_{0}\right)+c_{2}\left(x_{n}-x_{0}\right)\left(x_{n}-x_{1}\right)+\cdots+c_{n}\left(x_{n}-x_{0}\right) \cdots\left(x_{n}-x_{n-1}\right) \end{aligned} \nonumber$$

This system of linear equations is lower triangular as can be seen from the matrix form

$$\left(\begin{array}{ccccc} 1 & 0 & 0 & \cdots & 0 \\ 1 & \left(x_{1}-x_{0}\right) & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \left(x_{n}-x_{0}\right) & \left(x_{n}-x_{0}\right)\left(x_{n}-x_{1}\right) & \cdots & \left(x_{n}-x_{0}\right) \cdots\left(x_{n}-x_{n-1}\right) \end{array}\right)\left(\begin{array}{c} c_{0} \\ c_{1} \\ \vdots \\ c_{n} \end{array}\right) \nonumber$$

and so theoretically can be solved faster than the Vandermonde polynomial. In practice, however, there is little difference because polynomial interpolation is only useful when the number of points to be interpolated is small.

Midpoint rule

The midpoint rule makes use of only the first term in the Taylor series expansion. Here, we will determine the error in this approximation. Integrating,

$$\begin{aligned} I_{h}=h f(h / 2)+\int_{0}^{h}(&(x-h / 2) f^{\prime}(h / 2)+\frac{(x-h / 2)^{2}}{2} f^{\prime \prime}(h / 2) \\ &\left.+\frac{(x-h / 2)^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{(x-h / 2)^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d x . \end{aligned} \nonumber$$

Changing variables by letting $y=x-h / 2$ and $d y=d x$, and simplifying the integral depending on whether the integrand is even or odd, we have

$$\begin{aligned} I_{h}=& h f(h / 2) \\ &+\int_{-h / 2}^{h / 2}\left(y f^{\prime}(h / 2)+\frac{y^{2}}{2} f^{\prime \prime}(h / 2)+\frac{y^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{y^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d y \\ =& h f(h / 2)+\int_{0}^{h / 2}\left(y^{2} f^{\prime \prime}(h / 2)+\frac{y^{4}}{12} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d y \end{aligned} \nonumber$$

The integrals that we need here are

$$\int_{0}^{\frac{h}{2}} y^{2} d y=\frac{h^{3}}{24}, \quad \int_{0}^{\frac{h}{2}} y^{4} d y=\frac{h^{5}}{160} \nonumber$$

Therefore,

$$I_{h}=h f(h / 2)+\frac{h^{3}}{24} f^{\prime \prime}(h / 2)+\frac{h^{5}}{1920} f^{\prime \prime \prime \prime}(h / 2)+\ldots \nonumber$$

Trapezoidal rule

From the Taylor series expansion of $f(x)$ about $x=h / 2$, we have

$$f(0)=f(h / 2)-\frac{h}{2} f^{\prime}(h / 2)+\frac{h^{2}}{8} f^{\prime \prime}(h / 2)-\frac{h^{3}}{48} f^{\prime \prime \prime}(h / 2)+\frac{h^{4}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots, \nonumber$$

and

$$f(h)=f(h / 2)+\frac{h}{2} f^{\prime}(h / 2)+\frac{h^{2}}{8} f^{\prime \prime}(h / 2)+\frac{h^{3}}{48} f^{\prime \prime \prime}(h / 2)+\frac{h^{4}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots \ldots \nonumber$$

Adding and multiplying by $h / 2$ we obtain

$$\frac{h}{2}(f(0)+f(h))=h f(h / 2)+\frac{h^{3}}{8} f^{\prime \prime}(h / 2)+\frac{h^{5}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber$$

We now substitute for the first term on the right-hand-side using the midpoint rule formula:

$$\begin{array}{r} \frac{h}{2}(f(0)+f(h))=\left(I_{h}-\frac{h^{3}}{24} f^{\prime \prime}(h / 2)-\frac{h^{5}}{1920} f^{\prime \prime \prime \prime}(h / 2)\right) \\ +\frac{h^{3}}{8} f^{\prime \prime}(h / 2)+\frac{h^{5}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots, \end{array} \nonumber$$

and solving for $I_{h}$, we find

$$I_{h}=\frac{h}{2}(f(0)+f(h))-\frac{h^{3}}{12} f^{\prime \prime}(h / 2)-\frac{h^{5}}{480} f^{\prime \prime \prime \prime}(h / 2)+\ldots \nonumber$$

Simpson’s rule

To obtain Simpson’s rule, we combine the midpoint and trapezoidal rule to eliminate the error term proportional to $h^{3}$. Multiplying (6.3) by two and adding to (6.4), we obtain

$$3 I_{h}=h\left(2 f(h / 2)+\frac{1}{2}(f(0)+f(h))\right)+h^{5}\left(\frac{2}{1920}-\frac{1}{480}\right) f^{\prime \prime \prime \prime}(h / 2)+\ldots, \nonumber$$

or

$$I_{h}=\frac{h}{6}(f(0)+4 f(h / 2)+f(h))-\frac{h^{5}}{2880} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber$$

Usually, Simpson’s rule is written by considering the three consecutive points $0, h$ and $2 h$. Substituting $h \rightarrow 2 h$, we obtain the standard result

$$I_{2 h}=\frac{h}{3}(f(0)+4 f(h)+f(2 h))-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(h)+\ldots \nonumber$$

Trapezoidal rule

We suppose that the function $f(x)$ is known at the $n+1$ points labeled as $x_{0}, x_{1}, \ldots, x_{n}$, with the endpoints given by $x_{0}=a$ and $x_{n}=b$. Define

$$f_{i}=f\left(x_{i}\right), \quad h_{i}=x_{i+1}-x_{i} \nonumber$$

Then the integral of (6.1) may be decomposed as

$$\begin{aligned} \int_{a}^{b} f(x) d x &=\sum_{i=0}^{n-1} \int_{x_{i}}^{x_{i+1}} f(x) d x \\ &=\sum_{i=0}^{n-1} \int_{0}^{h_{i}} f\left(x_{i}+s\right) d s \end{aligned} \nonumber$$

where the last equality arises from the change-of-variables $s=x-x_{i}$. Applying the trapezoidal rule to the integral, we have

$$\int_{a}^{b} f(x) d x=\sum_{i=0}^{n-1} \frac{h_{i}}{2}\left(f_{i}+f_{i+1}\right) \nonumber$$

If the points are not evenly spaced, say because the data are experimental values, then the $h_{i}$ may differ for each value of $i$ and (6.6) is to be used directly.

However, if the points are evenly spaced, say because $f(x)$ can be computed, we have $h_{i}=h$, independent of $i$. We can then define

$$x_{i}=a+i h, \quad i=0,1, \ldots, n \nonumber$$

and since the end point $b$ satisfies $b=a+n h$, we have

$$h=\frac{b-a}{n} \nonumber$$

The composite trapezoidal rule for evenly space points then becomes

$$\begin{aligned} \int_{a}^{b} f(x) d x &=\frac{h}{2} \sum_{i=0}^{n-1}\left(f_{i}+f_{i+1}\right) \\ &=\frac{h}{2}\left(f_{0}+2 f_{1}+\cdots+2 f_{n-1}+f_{n}\right) \end{aligned} \nonumber$$

The first and last terms have a multiple of one; all other terms have a multiple of two; and the entire sum is multiplied by $h / 2$.

Simpson’s rule

We here consider the composite Simpson’s rule for evenly space points. We apply Simpson’s rule over intervals of $2 h$, starting from $a$ and ending at $b$ :

$$\begin{aligned} \int_{a}^{b} f(x) d x=\frac{h}{3}\left(f_{0}+4 f_{1}+f_{2}\right)+\frac{h}{3}\left(f_{2}+4 f_{3}+f_{4}\right)+& \ldots \\ &+\frac{h}{3}\left(f_{n-2}+4 f_{n-1}+f_{n}\right) \end{aligned} \nonumber$$

Note that $n$ must be even for this scheme to work. Combining terms, we have

$$\int_{a}^{b} f(x) d x=\frac{h}{3}\left(f_{0}+4 f_{1}+2 f_{2}+4 f_{3}+2 f_{4}+\cdots+4 f_{n-1}+f_{n}\right) \nonumber$$

The first and last terms have a multiple of one; the even indexed terms have a multiple of $2 ;$ the odd indexed terms have a multiple of $4 ;$ and the entire sum is multiplied by $h / 3$.

Initial value problem

One classic initial value problem is the $R C$ circuit. With $R$ the resistor and $C$ the capacitor, the differential equation for the charge $q$ on the capacitor is given by

$$R \frac{d q}{d t}+\frac{q}{C}=0 . \nonumber$$

If we consider the physical problem of a charged capacitor connected in a closed circuit to a resistor, then the initial condition is $q(0)=q_{0}$, where $q_{0}$ is the initial charge on the capacitor.

The differential equation (7.1) is separable, and separating and integrating from time $t=0$ to $t$ yields

$$\int_{q_{0}}^{q} \frac{d q}{q}=-\frac{1}{R C} \int_{0}^{t} d t \nonumber$$

which can be integrated and solved for $q=q(t)$ :

$$q(t)=q_{0} e^{-t / R C} . \nonumber$$

The classic second-order initial value problem is the $R L C$ circuit, with differential equation

$$L \frac{d^{2} q}{d t^{2}}+R \frac{d q}{d t}+\frac{q}{C}=0 . \nonumber$$

Here, a charged capacitor is connected to a closed circuit, and the initial conditions satisfy

$$q(0)=q_{0}, \quad \frac{d q}{d t}(0)=0 \nonumber$$

The solution is obtained for the second-order equation by the ansatz

$$q(t)=e^{r t} \nonumber$$

which results in the following so-called characteristic equation for $r$ :

$$L r^{2}+R r+\frac{1}{C}=0 \nonumber$$

If the two solutions for $r$ are distinct and real, then the two found exponential solutions can be multiplied by constants and added to form a general solution. The constants can then be determined by requiring the general solution to satisfy the two initial conditions. If the roots of the characteristic equation are complex or degenerate, a general solution to the differential equation can also be found.

Boundary value problems

The dimensionless equation for the temperature $y=y(x)$ along a linear heatconducting rod of length unity, and with an applied external heat source $f(x)$, is given by the differential equation

$$-\frac{d^{2} y}{d x^{2}}=f(x) \nonumber$$

with $0 \leq x \leq 1$. Boundary conditions are usually prescribed at the end points of the rod, and here we assume that the temperature at both ends are maintained at zero so that

$$y(0)=0, \quad y(1)=0 \nonumber$$

The assignment of boundary conditions at two separate points is called a twopoint boundary value problem, in contrast to the initial value problem where the boundary conditions are prescribed at only a single point. Two-point boundary value problems typically require a more sophisticated algorithm for a numerical solution than initial value problems.

Here, the solution of (7.2) can proceed by integration once $f(x)$ is specified. We assume that

$$f(x)=x(1-x) \nonumber$$

so that the maximum of the heat source occurs in the center of the rod, and goes to zero at the ends.

The differential equation can then be written as

$$\frac{d^{2} y}{d x^{2}}=-x(1-x) \nonumber$$

The first integration results in

$$\begin{aligned} \frac{d y}{d x} &=\int\left(x^{2}-x\right) d x \\ &=\frac{x^{3}}{3}-\frac{x^{2}}{2}+c_{1} \end{aligned} \nonumber$$

where $c_{1}$ is the first integration constant. Integrating again,

$$\begin{aligned} y(x) &=\int\left(\frac{x^{3}}{3}-\frac{x^{2}}{2}+c_{1}\right) d x \\ &=\frac{x^{4}}{12}-\frac{x^{3}}{6}+c_{1} x+c_{2} \end{aligned} \nonumber$$

where $c_{2}$ is the second integration constant. The two integration constants are determined by the boundary conditions. At $x=0$, we have

$$0=c_{2} \nonumber$$

and at $x=1$, we have

$$0=\frac{1}{12}-\frac{1}{6}+c_{1} \nonumber$$

so that $c_{1}=1 / 12$. Our solution is therefore

$$\begin{aligned} y(x) &=\frac{x^{4}}{12}-\frac{x^{3}}{6}+\frac{x}{12} \\ &=\frac{1}{12} x(1-x)\left(1+x-x^{2}\right) \end{aligned} \nonumber$$

The temperature of the rod is maximum at $x=1 / 2$ and goes smoothly to zero at the ends.

Eigenvalue problem

The classic eigenvalue problem obtained by solving the wave equation by separation of variables is given by

$$\frac{d^{2} y}{d x^{2}}+\lambda^{2} y=0 \nonumber$$

with the two-point boundary conditions $y(0)=0$ and $y(1)=0$. Notice that $y(x)=0$ satisfies both the differential equation and the boundary conditions. Other nonzero solutions for $y=y(x)$ are possible only for certain discrete values of $\lambda$. These values of $\lambda$ are called the eigenvalues of the differential equation.

We proceed by first finding the general solution to the differential equation. It is easy to see that this solution is

$$y(x)=A \cos \lambda x+B \sin \lambda x \nonumber$$

Imposing the first boundary condition at $x=0$, we obtain

$$A=0 \nonumber$$

The second boundary condition at $x=1$ results in

$$B \sin \lambda=0 \nonumber$$

Since we are searching for a solution where $y=y(x)$ is not identically zero, we must have

$$\lambda=\pi, 2 \pi, 3 \pi, \ldots \nonumber$$

The corresponding negative values of $\lambda$ are also solutions, but their inclusion only changes the corresponding values of the unknown $B$ constant. A linear superposition of all the solutions results in the general solution

$$y(x)=\sum_{n=1}^{\infty} B_{n} \sin n \pi x . \nonumber$$

For each eigenvalue $n \pi$, we say there is a corresponding eigenfunction $\sin n \pi x$. When the differential equation can not be solved analytically, a numerical method should be able to solve for both the eigenvalues and eigenfunctions.

Euler method

The Euler method is the most straightforward method to integrate a differential equation. Consider the first-order differential equation

$$\dot{x}=f(t, x), \nonumber$$

with the initial condition $x(0)=x_{0}$. Define $t_{n}=n \Delta t$ and $x_{n}=x\left(t_{n}\right)$. A Taylor series expansion of $x_{n+1}$ results in

$$\begin{aligned} x_{n+1} &=x\left(t_{n}+\Delta t\right) \\ &=x\left(t_{n}\right)+\Delta t \dot{x}\left(t_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) \\ &=x\left(t_{n}\right)+\Delta t f\left(t_{n}, x_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) \end{aligned} \nonumber$$

The Euler Method is therefore written as

$$x_{n+1}=x\left(t_{n}\right)+\Delta t f\left(t_{n}, x_{n}\right) \nonumber$$

We say that the Euler method steps forward in time using a time-step $\Delta t$, starting from the initial value $x_{0}=x(0)$. The local error of the Euler Method is $\mathrm{O}\left(\Delta t^{2}\right)$. The global error, however, incurred when integrating to a time $T$, is a factor of $1 / \Delta t$ larger and is given by $\mathrm{O}(\Delta t)$. It is therefore customary to call the Euler Method a first-order method.

Modified Euler method

This method is of a type that is called a predictor-corrector method. It is also the first of what are Runge-Kutta methods. As before, we want to solve (7.3). The idea is to average the value of $\dot{x}$ at the beginning and end of the time step. That is, we would like to modify the Euler method and write

$$x_{n+1}=x_{n}+\frac{1}{2} \Delta t\left(f\left(t_{n}, x_{n}\right)+f\left(t_{n}+\Delta t, x_{n+1}\right)\right) \nonumber$$

The obvious problem with this formula is that the unknown value $x_{n+1}$ appears on the right-hand-side. We can, however, estimate this value, in what is called the predictor step. For the predictor step, we use the Euler method to find

$$x_{n+1}^{p}=x_{n}+\Delta t f\left(t_{n}, x_{n}\right) \nonumber$$

The corrector step then becomes

$$x_{n+1}=x_{n}+\frac{1}{2} \Delta t\left(f\left(t_{n}, x_{n}\right)+f\left(t_{n}+\Delta t, x_{n+1}^{p}\right)\right) \nonumber$$

The Modified Euler Method can be rewritten in the following form that we will later identify as a Runge-Kutta method:

$$\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}\right) \\ k_{2} &=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{1}\right) \\ x_{n+1} &=x_{n}+\frac{1}{2}\left(k_{1}+k_{2}\right) \end{aligned} \nonumber$$

Second-order Runge-Kutta methods

We now derive all second-order Runge-Kutta methods. Higher-order methods can be similarly derived, but require substantially more algebra.

We consider the differential equation given by (7.3). A general second-order Runge-Kutta method may be written in the form

$$\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}\right) \\ k_{2} &=\Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta k_{1}\right) \\ x_{n+1} &=x_{n}+a k_{1}+b k_{2} \end{aligned} \nonumber$$

with $\alpha, \beta, a$ and $b$ constants that define the particular second-order Runge-Kutta method. These constants are to be constrained by setting the local error of the second-order Runge-Kutta method to be $\mathrm{O}\left(\Delta t^{3}\right)$. Intuitively, we might guess that two of the constraints will be $a+b=1$ and $\alpha=\beta$.

We compute the Taylor series of $x_{n+1}$ directly, and from the Runge-Kutta method, and require them to be the same to order $\Delta t^{2}$. First, we compute the Taylor series of $x_{n+1}$. We have

$$\begin{aligned} x_{n+1} &=x\left(t_{n}+\Delta t\right) \\ &=x\left(t_{n}\right)+\Delta t \dot{x}\left(t_{n}\right)+\frac{1}{2}(\Delta t)^{2} \ddot{x}\left(t_{n}\right)+\mathrm{O}\left(\Delta t^{3}\right) \end{aligned} \nonumber$$

Now,

$$\dot{x}\left(t_{n}\right)=f\left(t_{n}, x_{n}\right) . \nonumber$$

The second derivative is more complicated and requires partial derivatives. We have

$$\begin{aligned} \ddot{x}\left(t_{n}\right) &\left.=\frac{d}{d t} f(t, x(t))\right]_{t=t_{n}} \\ &=f_{t}\left(t_{n}, x_{n}\right)+\dot{x}\left(t_{n}\right) f_{x}\left(t_{n}, x_{n}\right) \\ &=f_{t}\left(t_{n}, x_{n}\right)+f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right) \end{aligned} \nonumber$$

Therefore,

$$x_{n+1}=x_{n}+\Delta t f\left(t_{n}, x_{n}\right)+\frac{1}{2}(\Delta t)^{2}\left(f_{t}\left(t_{n}, x_{n}\right)+f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)\right) \nonumber$$

Second, we compute $x_{n+1}$ from the Runge-Kutta method given by (7.5). Substituting in $k_{1}$ and $k_{2}$, we have

$$x_{n+1}=x_{n}+a \Delta t f\left(t_{n}, x_{n}\right)+b \Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta \Delta t f\left(t_{n}, x_{n}\right)\right) . \nonumber$$

We Taylor series expand using

$$\begin{aligned} f\left(t_{n}+\alpha \Delta t, x_{n}+\beta \Delta t f\left(t_{n}, x_{n}\right)\right) & \\ =f\left(t_{n}, x_{n}\right)+\alpha \Delta t f_{t}\left(t_{n}, x_{n}\right)+\beta \Delta t f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) \end{aligned} \nonumber$$

The Runge-Kutta formula is therefore

$$\begin{aligned} x_{n+1}=x_{n}+(a+b) & \Delta t f\left(t_{n}, x_{n}\right) \\ &+(\Delta t)^{2}\left(\alpha b f_{t}\left(t_{n}, x_{n}\right)+\beta b f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) \end{aligned} \nonumber$$

Comparing (7.6) and (7.7), we find

$$\begin{aligned} a+b &=1 \\ \alpha b &=1 / 2 \\ \beta b &=1 / 2 \end{aligned} \nonumber$$

There are three equations for four parameters, and there exists a family of secondorder Runge-Kutta methods.

The Modified Euler Method given by (7.4) corresponds to $\alpha=\beta=1$ and $a=$ $b=1 / 2$. Another second-order Runge-Kutta method, called the Midpoint Method, corresponds to $\alpha=\beta=1 / 2, a=0$ and $b=1 .$ This method is written as

$$\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}\right) \\ k_{2} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}\right) \\ x_{n+1} &=x_{n}+k_{2} \end{aligned} \nonumber$$

Higher-order Runge-Kutta methods

The general second-order Runge-Kutta method was given by (7.5). The general form of the third-order method is given by

$$\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}\right) \\ k_{2} &=\Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta k_{1}\right) \\ k_{3} &=\Delta t f\left(t_{n}+\gamma \Delta t, x_{n}+\delta k_{1}+\epsilon k_{2}\right) \\ x_{n+1} &=x_{n}+a k_{1}+b k_{2}+c k_{3} \end{aligned} \nonumber$$

The following constraints on the constants can be guessed: $\alpha=\beta, \gamma=\delta+\epsilon$, and $a+b+c=1$. Remaining constraints need to be derived.

The fourth-order method has a $k_{1}, k_{2}, k_{3}$ and $k_{4}$. The fifth-order method requires up to $k_{6}$. The table below gives the order of the method and the number of stages required.

Because of the jump in the number of stages required between the fourth-order and fifth-order method, the fourth-order Runge-Kutta method has some appeal. The general fourth-order method starts with 13 constants, and one then finds 11 constraints. A particularly simple fourth-order method that has been widely used is given by

$$\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}\right) \\ k_{2} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}\right) \\ k_{3} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}\right) \\ k_{4} &=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{3}\right) \\ x_{n+1} &=x_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right) \end{aligned} \nonumber$$

Adaptive Runge-Kutta Methods

As in adaptive integration, it is useful to devise an ode integrator that automatically finds the appropriate $\Delta t$. The Dormand-Prince Method, which is implemented in MATLAB’s ode45.m, finds the appropriate step size by comparing the results of a fifth-order and fourth-order method. It requires six function evaluations per time step, and constructs both a fifth-order and a fourth-order method from these function evaluations.

Suppose the fifth-order method finds $x_{n+1}$ with local error $\mathrm{O}\left(\Delta t^{6}\right)$, and the fourth-order method finds $x_{n+1}^{\prime}$ with local error $\mathrm{O}\left(\Delta t^{5}\right)$. Let $\varepsilon$ be the desired error tolerance of the method, and let $e$ be the actual error. We can estimate $e$ from the difference between the fifth-and fourth-order methods; that is,

$$e=\left|x_{n+1}-x_{n+1}^{\prime}\right| \nonumber$$

Now $e$ is of $\mathrm{O}\left(\Delta t^{5}\right)$, where $\Delta t$ is the step size taken. Let $\Delta \tau$ be the estimated step size required to get the desired error $\varepsilon$. Then we have

$$e / \varepsilon=(\Delta t)^{5} /(\Delta \tau)^{5} \nonumber$$

or solving for $\Delta \tau$,

$$\Delta \tau=\Delta t\left(\frac{\varepsilon}{e}\right)^{1 / 5} \nonumber$$

On the one hand, if $e<\varepsilon$, then we accept $x_{n+1}$ and do the next time step using the larger value of $\Delta \tau$. On the other hand, if $e>\varepsilon$, then we reject the integration step and redo the time step using the smaller value of $\Delta \tau$. In practice, one usually increases the time step slightly less and decreases the time step slightly more to prevent the waste of too many failed time steps.

System of differential equations

Our numerical methods can be easily adapted to solve higher-order differential equations, or equivalently, a system of differential equations. First, we show how a second-order differential equation can be reduced to two first-order equations. Consider

$$\ddot{x}=f(t, x, \dot{x}) . \nonumber$$

This second-order equation can be rewritten as two first-order equations by defining $u=\dot{x} .$ We then have the system

$$\begin{aligned} &\dot{x}=u, \\ &\dot{u}=f(t, x, u) . \end{aligned} \nonumber$$

This trick also works for higher-order equation. For another example, the thirdorder equation

$$\dddot x=f(t, x, \dot{x}, \ddot{x}), \nonumber$$

can be written as

$$\begin{aligned} &\dot{x}=u \\ &\dot{u}=v \\ &\dot{v}=f(t, x, u, v) . \end{aligned} \nonumber$$

Now, we show how to generalize Runge-Kutta methods to a system of differential equations. As an example, consider the following system of two odes,

$$\begin{aligned} &\dot{x}=f(t, x, y), \\ &\dot{y}=g(t, x, y), \end{aligned} \nonumber$$

with the initial conditions $x(0)=x_{0}$ and $y(0)=y_{0}$. The generalization of the commonly used fourth-order Runge-Kutta method would be

$$\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}, y_{n}\right) \\ l_{1} &=\Delta t g\left(t_{n}, x_{n}, y_{n}\right) \\ k_{2} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}, y_{n}+\frac{1}{2} l_{1}\right) \\ l_{2} &=\Delta t g\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}, y_{n}+\frac{1}{2} l_{1}\right) \\ y_{3} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}, y_{n}+\frac{1}{2} l_{2}\right) \\ l_{3} &=\Delta t g\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}, y_{n}+\frac{1}{2} l_{2}\right) \\ y_{n+1} &=y_{n}+\frac{1}{6}\left(l_{1}+2 l_{2}+2 l_{3}+l_{4}\right) \\ k_{4} &=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{3}, y_{n}+l_{3}\right) \\ l_{4} &=\Delta t g\left(t_{n}+\Delta t, x_{n}+k_{3}, y_{n}+l_{3}\right) \\ &=x_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right) \\ &=\Delta \end{aligned} \nonumber$$

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$$

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

Finite difference method

We first consider solving $(7.9)$ with the homogeneous boundary conditions $y(0)=$ $y(1)=0 .$ In this case, we have already shown that the eigenvalues of $(7.9)$ are given by $\lambda=\pi, 2 \pi, 3 \pi, \ldots$.

With $n$ interior points, we have $x_{i}=i h$ for $i=0, \ldots, n+1$, and $h=1 /(n+1)$. Using the centered-finite-difference approximation for the second derivative, (7.9) becomes

$$y_{i-1}-2 y_{i}+y_{i+1}=-h^{2} \lambda^{2} y_{i} \nonumber$$

Applying the boundary conditions $y_{0}=y_{n+1}=0$, the first equation with $i=1$, and the last equation with $i=n$, are given by

$$\begin{aligned} -2 y_{1}+y_{2} &=-h^{2} \lambda^{2} y_{1} \\ y_{n-1}-2 y_{n} &=-h^{2} \lambda^{2} y_{n} \end{aligned} \nonumber$$

The remaining $n-2$ equations are given by $(7.10)$ for $i=2, \ldots, n-1$ It is of interest to see how the solution develops with increasing $n$. The smallest possible value is $n=1$, corresponding to a single interior point, and since $h=1 / 2$ we have

$$-2 y_{1}=-\frac{1}{4} \lambda^{2} y_{1} \nonumber$$

so that $\lambda^{2}=8$, or $\lambda=2 \sqrt{2}=2.8284$. This is to be compared to the first eigenvalue which is $\lambda=\pi$. When $n=2$, we have $h=1 / 3$, and the resulting two equations written in matrix form are given by

$$\left(\begin{array}{rr} -2 & 1 \\ 1 & -2 \end{array}\right)\left(\begin{array}{l} y_{1} \\ y_{2} \end{array}\right)=-\frac{1}{9} \lambda^{2}\left(\begin{array}{l} y_{1} \\ y_{2} \end{array}\right) \nonumber$$

This is a matrix eigenvalue problem with the eigenvalue given by $\mu=-\lambda^{2} / 9$. The solution for $\mu$ is arrived at by solving

$$\operatorname{det}\left(\begin{array}{cc} -2-\mu & 1 \\ 1 & -2-\mu \end{array}\right)=0 \nonumber$$

with resulting quadratic equation

$$(2+\mu)^{2}-1=0 \nonumber$$

The solutions are $\mu=-1,-3$, and since $\lambda=3 \sqrt{-\mu}$, we have $\lambda=3,3 \sqrt{3}=5.1962$. These two eigenvalues serve as rough approximations to the first two eigenvalues $\pi$ and $2 \pi$.

With A an $n$-by-n matrix, the MATLAB variable mu=eig(A) is a vector containing the $n$ eigenvalues of the matrix A. The built-in function eig.m can therefore be used to find the eigenvalues. With $n$ grid points, the smaller eigenvalues will converge more rapidly than the larger ones.

We can also consider boundary conditions on the derivative, or mixed boundary conditions. For example, consider the mixed boundary conditions given by $y(0)=0$ and $y^{\prime}(1)=0$. The eigenvalues of (7.9) can then be determined analytically to be $\lambda_{i}=(i-1 / 2) \pi$, with $i$ a natural number.

The difficulty we now face is how to implement a boundary condition on the derivative. Our computation of $y^{\prime \prime}$ uses a second-order method, and we would like the computation of the first derivative to also be second order. The condition $y^{\prime}(1)=0$ occurs on the right-most boundary, and we can make use of the secondorder backward-difference approximation to the derivative that we have previously derived. This finite-difference approximation for $y^{\prime}(1)$ can be written as

$$y_{n+1}^{\prime}=\frac{3 y_{n+1}-4 y_{n}+y_{n-1}}{2 h} . \nonumber$$

Now, the $n$th finite-difference equation was given by

$$y_{n-1}-2 y_{n}+y_{n+1}=-h^{2} y_{n}, \nonumber$$

and we now replace the value $y_{n+1}$ using $(7.11)$; that is,

$$y_{n+1}=\frac{1}{3}\left(2 h y_{n+1}^{\prime}+4 y_{n}-y_{n-1}\right) . \nonumber$$

Implementing the boundary condition $y_{n+1}^{\prime}=0$, we have

$$y_{n+1}=\frac{4}{3} y_{n}-\frac{1}{3} y_{n-1} \nonumber$$

Therefore, the $n$th finite-difference equation becomes

$$\frac{2}{3} y_{n-1}-\frac{2}{3} y_{n}=-h^{2} \lambda^{2} y_{n} \nonumber$$

For example, when $n=2$, the finite difference equations become

$$\left(\begin{array}{rr} -2 & 1 \\ \frac{2}{3} & -\frac{2}{3} \end{array}\right)\left(\begin{array}{l} y_{1} \\ y_{2} \end{array}\right)=-\frac{1}{9} \lambda^{2}\left(\begin{array}{l} y_{1} \\ y_{2} \end{array}\right) \nonumber$$

The eigenvalues of the matrix are now the solution of

$$(\mu+2)\left(\mu+\frac{2}{3}\right)-\frac{2}{3}=0 \nonumber$$

Or

$$3 \mu^{2}+8 \mu+2=0 \nonumber$$

Therefore, $\mu=(-4 \pm \sqrt{10}) / 3$, and we find $\lambda=1.5853,4.6354$, which are approximations to $\pi / 2$ and $3 \pi / 2$.

Shooting method

We apply the shooting method to solve (7.9) with boundary conditions $y(0)=$ $y(1)=0$. The initial value problem to solve is

$$\begin{aligned} &y^{\prime}=z \\ &z^{\prime}=-\lambda^{2} y \end{aligned} \nonumber$$

with known boundary condition $y(0)=0$ and an unknown boundary condition on $y^{\prime}(0)$. In fact, any nonzero boundary condition on $y^{\prime}(0)$ can be chosen: the differential equation is linear and the boundary conditions are homogeneous, so that if $y(x)$ is an eigenfunction then so is $A y(x)$. What we need to find here is the value of $\lambda$ such that $y(1)=0 .$ In other words, choosing $y^{\prime}(0)=1$, we solve

$$F(\lambda)=0 \nonumber$$

where $F(\lambda)=y(1)$, obtained by solving the initial value problem. Again, an iteration for the roots of $F(\lambda)$ can be done using the Secant Method. For the eigenvalue problem, there are an infinite number of roots, and the choice of the two initial guesses for $\lambda$ will then determine to which root the iteration will converge.

For this simple problem, it is possible to write explicitly the equation $F(\lambda)=0$. The general solution to $(7.9)$ is given by

$$y(x)=A \cos (\lambda x)+B \sin (\lambda x) . \nonumber$$

The initial condition $y(0)=0$ yields $A=0 .$ The initial condition $y^{\prime}(0)=1$ yields

$$B=1 / \lambda \nonumber$$

Therefore, the solution to the initial value problem is

$$y(x)=\frac{\sin (\lambda x)}{\lambda} \nonumber$$

The function $F(\lambda)=y(1)$ is therefore given by

$$F(\lambda)=\frac{\sin \lambda}{\lambda} \nonumber$$

and the roots occur when $\lambda=\pi, 2 \pi, \ldots .$

If the boundary conditions were $y(0)=0$ and $y^{\prime}(1)=0$, for example, then we would simply redefine $F(\lambda)=y^{\prime}(1)$. We would then have

$$F(\lambda)=\frac{\cos \lambda}{\lambda} \nonumber$$

and the roots occur when $\lambda=\pi / 2,3 \pi / 2, \ldots$