Latex

Last updated

May 31, 2022
Save as PDF
- Licensing
- 1: IEEE Arithmetic

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

$\newcommand{\vectorA}[1]{\vec{#1}} % arrow$

$\newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow$

$\newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$

$\newcommand{\vectorC}[1]{\textbf{#1}}$

$\newcommand{\vectorD}[1]{\overrightarrow{#1}}$

$\newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}}$

$\newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}}$

$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$

$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$

$\newcommand{\avec}{\mathbf a}$

$\newcommand{\bvec}{\mathbf b}$

$\newcommand{\cvec}{\mathbf c}$

$\newcommand{\dvec}{\mathbf d}$

$\newcommand{\dtil}{\widetilde{\mathbf d}}$

$\newcommand{\evec}{\mathbf e}$

$\newcommand{\fvec}{\mathbf f}$

$\newcommand{\nvec}{\mathbf n}$

$\newcommand{\pvec}{\mathbf p}$

$\newcommand{\qvec}{\mathbf q}$

$\newcommand{\svec}{\mathbf s}$

$\newcommand{\tvec}{\mathbf t}$

$\newcommand{\uvec}{\mathbf u}$

$\newcommand{\vvec}{\mathbf v}$

$\newcommand{\wvec}{\mathbf w}$

$\newcommand{\xvec}{\mathbf x}$

$\newcommand{\yvec}{\mathbf y}$

$\newcommand{\zvec}{\mathbf z}$

$\newcommand{\rvec}{\mathbf r}$

$\newcommand{\mvec}{\mathbf m}$

$\newcommand{\zerovec}{\mathbf 0}$

$\newcommand{\onevec}{\mathbf 1}$

$\newcommand{\real}{\mathbb R}$

$\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}$

$\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}$

$\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}$

$\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}$

$\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}$

$\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}$

$\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}$

$\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}$

$\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}$

$\newcommand{\laspan}[1]{\text{Span}\{#1\}}$

$\newcommand{\bcal}{\cal B}$

$\newcommand{\ccal}{\cal C}$

$\newcommand{\scal}{\cal S}$

$\newcommand{\wcal}{\cal W}$

$\newcommand{\ecal}{\cal E}$

$\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}$

$\newcommand{\gray}[1]{\color{gray}{#1}}$

$\newcommand{\lgray}[1]{\color{lightgray}{#1}}$

$\newcommand{\rank}{\operatorname{rank}}$

$\newcommand{\row}{\text{Row}}$

$\newcommand{\col}{\text{Col}}$

$\renewcommand{\row}{\text{Row}}$

$\newcommand{\nul}{\text{Nul}}$

$\newcommand{\var}{\text{Var}}$

$\newcommand{\corr}{\text{corr}}$

$\newcommand{\len}[1]{\left|#1\right|}$

$\newcommand{\bbar}{\overline{\bvec}}$

$\newcommand{\bhat}{\widehat{\bvec}}$

$\newcommand{\bperp}{\bvec^\perp}$

$\newcommand{\xhat}{\widehat{\xvec}}$

$\newcommand{\vhat}{\widehat{\vvec}}$

$\newcommand{\uhat}{\widehat{\uvec}}$

$\newcommand{\what}{\widehat{\wvec}}$

$\newcommand{\Sighat}{\widehat{\Sigma}}$

$\newcommand{\lt}{<}$

$\newcommand{\gt}{>}$

$\newcommand{\amp}{&}$

$\definecolor{fillinmathshade}{gray}{0.9}$

Chapter 1 Arithmetic

Definitions

$\begin{array}{lll} \text { Bit } & =0 \text { or } 1 & \\ \text { Byte } & =8 \text { bits } & \\ \text { Word } & =\text { Reals: } & 4 \text { bytes (single precision) } \\ & & 8 \text { bytes (double precision) } \\ & =\text { Integers: } & 1,2,4, \text { or } 8 \text { byte signed } \\ & & 1,2,4, \text { or } 8 \text { byte unsigned } \end{array} \nonumber$

Numbers with a decimal or binary point

$\begin{array}{lccccccccc} & \square & \square & \square & \square & \square & \square & \square & \square & \square \\ \text { Decimal: } & 10^{3} & 10^{2} & 10^{1} & 10^{0} & & 10^{-1} & 10^{-2} & 10^{-3} & 10^{-4} \\ \text { Binary: } & 2^{3} & 2^{2} & 2^{1} & 2^{0} & & 2^{-1} & 2^{-2} & 2^{-3} & 2^{-4} \end{array} \nonumber$

Examples of binary numbers

$\begin{array}{cc}\text { Decimal } & \text { Binary } \\ 1 & 1 \\ 2 & 10 \\ 3 & 11 \\ 4 & 100 \\ 0.5 & 0.1 \\ 1.5 & 1.1\end{array}$

numbers

$\{0,1,2,3, \ldots, 9,10,11,12,13,14,15\}=\{0,1,2,3 \ldots \ldots . .9, \mathrm{a}, \mathrm{b}, \mathrm{c}, \mathrm{d}, \mathrm{e}, \mathrm{f}\} \nonumber$

unsigned integers as hex numbers

$\begin{array}{ccc}\text { Decimal } & \text { Binary } & \text { Hex } \\ 1 & 0001 & 1 \\ 2 & 0010 & 2 \\ 3 & 0011 & 3 \\ \vdots & \vdots & \vdots \\ 10 & 1010 & \text { a } \\ \vdots & \vdots & \vdots \\ 15 & 1111 & \text { f }\end{array}$

single precision format:

where $s=\operatorname{sign}$

$e=\text { biased exponent } \nonumber$

$\mathrm{p}=\mathrm{e}-127=\text { exponent } \nonumber$

$\text { 1.f = significand (use binary point) } \nonumber$

Special numbers

Smallest exponent: Largest exponent:

reserved reserved

Smallest positive normal number

bin: 00000000000000000000000000000000

hex: 00000000

Examples of computer numbers

What is $1.0,2.0 \& 1 / 2$ in hex?

$1.0=(-1)^{0} \times 2^{(127-127)} \times 1.0$

Therefore, $s=0, e=01111111, f=00000000000000000000000$

bin: 00111111100000000000000000000000

hex: 3 f 800000

$2.0=(-1)^{0} \times 2^{(128-127)} \times 1.0$

Therefore, $s=0, e=10000000, f=00000000000000000000000$

bin: 010000000100000000000000000000000

hex: 40000000

$1 / 2=(-1)^{0} \times 2^{(126-127)} \times 1.0$

Therefore, $s=0, e=01111110, f=00000000000000000000000$

bin: 00111111000000000000000000000000

hex: 3 f00 0000

Inexact numbers

Example:

$\frac{1}{3}=(-1)^{0} \times \frac{1}{4} \times\left(1+\frac{1}{3}\right) \text {, } \nonumber$

so that $p=e-127=-2$ and $e=125=128-3$ , or in binary, $e=01111101$ . How is $f=1 / 3$ represented in binary? To compute binary number, multiply successively by 2 as follows:

$\begin{array}{lr}0.333 \ldots & 0 . \\ 0.666 \ldots & 0.0 \\ 1.333 \ldots & 0.01 \\ 0.666 \ldots & 0.010 \\ 1.333 \ldots & 0.0101 \\ \text { etc. } & \end{array}$

so that $1 / 3$ exactly in binary is $0.010101 \ldots$ With only 23 bits to represent $f$ , the number is inexact and we have

$f=01010101010101010101011 \nonumber$

where we have rounded to the nearest binary number (here, rounded up). The machine number $1 / 3$ is then represented as

$00111110101010101010101010101011 \nonumber$

or in hex

$\text { Зeaaaaab. } \nonumber$

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

Machine epsilon

Machine epsilon $\left(\epsilon_{\mathrm{mach}}\right)$ is the distance between 1 and the next largest number. If $0 \leq \delta<\epsilon_{\mathrm{mach}} / 2$ , then $1+\delta=1$ in computer math. Also since

$x+y=x(1+y / x) \nonumber$

if $0 \leq y / x<\epsilon_{\operatorname{mach}} / 2$ , then $x+y=x$ in computer math.

Find $\epsilon_{\text {mach }}$

The number 1 in the IEEE format is written as

$1=2^{0} \times 1.000 \ldots 0, \nonumber$

with $230^{\prime}$ s following the binary point. The number just larger than 1 has a 1 in the 23 rd position after the decimal point. Therefore,

$\epsilon_{\text {mach }}=2^{-23} \approx 1.192 \times 10^{-7} \nonumber$

What is the distance between 1 and the number just smaller than 1? Here, the number just smaller than one can be written as

$2^{-1} \times 1.111 \ldots 1=2^{-1}\left(1+\left(1-2^{-23}\right)\right)=1-2^{-24} \nonumber$

Therefore, this distance is $2^{-24}=\epsilon_{\operatorname{mach}} / 2$ .

The spacing between numbers is uniform between powers of 2, with logarithmic spacing of the powers of 2 . That is, the spacing of numbers between 1 and 2 is $2^{-23}$ , between 2 and 4 is $2^{-22}$ , between 4 and 8 is $2^{-21}$ , etc. This spacing changes for denormal numbers, where the spacing is uniform all the way down to zero.

Find the machine number just greater than 5

A rough estimate would be $5\left(1+\epsilon_{\text {mach }}\right)=5+5 \epsilon_{\text {mach }}$ , but this is not exact. The exact answer can be found by writing

$5=2^{2}\left(1+\frac{1}{4}\right) \nonumber$

so that the next largest number is

$2^{2}\left(1+\frac{1}{4}+2^{-23}\right)=5+2^{-21}=5+4 \epsilon_{\text {mach }} \nonumber$

double precision format

Most computations take place in double precision, where round-off error is reduced, and all of the above calculations in single precision can be repeated for double precision. The format is

where $\mathrm{s}=\operatorname{sign}$

$\mathrm{e}=$ biased exponent

$\mathrm{p}=\mathrm{e}-1023=$ exponent

1.f = significand (use binary point)

$1.12$ Roundoff error example

Consider solving the quadratic equation

$x^{2}+2 b x-1=0 \nonumber$

where $b$ is a parameter. The quadratic formula yields the two solutions

$x_{\pm}=-b \pm \sqrt{b^{2}+1} \nonumber$

Consider the solution with $b>0$ and $x>0$ (the $x_{+}$ solution) given by

$x=-b+\sqrt{b^{2}+1} \nonumber$

As $b \rightarrow \infty$

$\begin{aligned} x &=-b+\sqrt{b^{2}+1} \\ &=-b+b \sqrt{1+1 / b^{2}} \\ &=b\left(\sqrt{1+1 / b^{2}}-1\right) \\ & \approx b\left(1+\frac{1}{2 b^{2}}-1\right) \\ &=\frac{1}{2 b} . \end{aligned} \nonumber$

Now in double precision, realmin $\approx 2.2 \times 10^{-308}$ and we would like $x$ to be accurate to this value before it goes to 0 via denormal numbers. Therefore, $x$ should be computed accurately to $b \approx 1 /(2 \times$ realmin $) \approx 2 \times 10^{307}$ . What happens if we compute (1.1) directly? Then $x=0$ when $b^{2}+1=b^{2}$ , or $1+1 / b^{2}=1 .$ That is $1 / b^{2}=\epsilon_{\text {mach }} / 2$ , or $b=\sqrt{2} / \sqrt{\epsilon_{\text {mach }}} \approx 10^{8} .$ For a subroutine written to compute the solution of a quadratic for a general user, this is not good enough. The way for a software designer to solve this problem is to compute the solution for $x$ as

$x=\frac{1}{b\left(1+\sqrt{1+1 / b^{2}}\right)} \nonumber$

In this form, if $1+1 / b^{2}=1$ , then $x=1 / 2 b$ which is the correct asymptotic form.

Chapter 2 Finding

Solve $f(x)=0$ for $x$ , when an explicit analytical solution is impossible.

Bisection Method

The bisection method is the easiest to numerically implement and almost always works. The main disadvantage is that convergence is slow. If the bisection method results in a computer program that runs too slow, then other faster methods may be chosen; otherwise it is a good choice of method.

We want to construct a sequence $x_{0}, x_{1}, x_{2}, \ldots$ that converges to the root $x=r$ that solves $f(x)=0$ . We choose $x_{0}$ and $x_{1}$ such that $x_{0}<r<x_{1}$ . We say that $x_{0}$ and $x_{1}$ bracket the root. With $f(r)=0$ , we want $f\left(x_{0}\right)$ and $f\left(x_{1}\right)$ to be of opposite sign, so that $f\left(x_{0}\right) f\left(x_{1}\right)<0$ . We then assign $x_{2}$ to be the midpoint of $x_{0}$ and $x_{1}$ , that is $x_{2}=\left(x_{0}+x_{1}\right) / 2$ , or

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

The sign of $f\left(x_{2}\right)$ can then be determined. The value of $x_{3}$ is then chosen as either the midpoint of $x_{0}$ and $x_{2}$ or as the midpoint of $x_{2}$ and $x_{1}$ , depending on whether $x_{0}$ and $x_{2}$ bracket the root, or $x_{2}$ and $x_{1}$ bracket the root. The root, therefore, stays bracketed at all times. The algorithm proceeds in this fashion and is typically stopped when the increment to the left side of the bracket (above, given by $\left(x_{1}-\right.$ $\left.x_{0}\right) / 2$ ) is smaller than some required precision.

$2.2$ Newton’s Method

This is the fastest method, but requires analytical computation of the derivative of $f(x)$ . Also, the method may not always converge to the desired root.

We can derive Newton’s Method graphically, or by a Taylor series. We again want to construct a sequence $x_{0}, x_{1}, x_{2}, \ldots$ that converges to the root $x=r$ . Consider the $x_{n+1}$ member of this sequence, and Taylor series expand $f\left(x_{n+1}\right)$ about the point $x_{n}$ . We have

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

To determine $x_{n+1}$ , we drop the higher-order terms in the Taylor series, and assume $f\left(x_{n+1}\right)=0 .$ Solving for $x_{n+1}$ , we have

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

Starting Newton’s Method requires a guess for $x_{0}$ , hopefully close to the root $x=r .$

Secant Method

The Secant Method is second best to Newton’s Method, and is used when a faster convergence than Bisection is desired, but it is too difficult or impossible to take an analytical derivative of the function $f(x)$ . We write in place of $f^{\prime}\left(x_{n}\right)$ ,

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

Starting the Secant Method requires a guess for both $x_{0}$ and $x_{1}$ .

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.

Order of convergence

Let $r$ be the root and $x_{n}$ be the $n$ th approximation to the root. Define the error as

$\epsilon_{n}=r-x_{n} \nonumber$

If for large $n$ we have the approximate relationship

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

with $k$ a positive constant, then we say the root-finding numerical method is of order $p$ . Larger values of $p$ correspond to faster convergence to the root. The order of convergence of bisection is one: the error is reduced by approximately a factor of 2 with each iteration so that

$\left|\epsilon_{n+1}\right|=\frac{1}{2}\left|\epsilon_{n}\right| . \nonumber$

We now find the order of convergence for Newton’s Method and for the Secant Method.

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$

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

Chapter 3 Systems of equations

Consider the system of $n$ linear equations and $n$ unknowns, given by

$\begin{array}{cc} a_{11} x_{1}+a_{12} x_{2}+\cdots+a_{1 n} x_{n} & =b_{1} \\ a_{21} x_{1}+a_{22} x_{2}+\cdots+a_{2 n} x_{n} & =b_{2} \\ \vdots & \vdots \\ a_{n 1} x_{1}+a_{n 2} x_{2}+\cdots+a_{n n} x_{n} & =b_{n} \end{array} \nonumber$

We can write this system as the matrix equation

$\mathrm{A} \mathbf{x}=\mathbf{b} \nonumber$

with

$\mathrm{A}=\left(\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} \\ a_{21} & a_{22} & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n} \end{array}\right), \quad \mathbf{x}=\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{array}\right), \quad \mathbf{b}=\left(\begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{array}\right) \nonumber$

Gaussian Elimination

The standard numerical algorithm to solve a system of linear equations is called Gaussian Elimination. We can illustrate this algorithm by example.

Consider the system of equations

$\begin{aligned} &-3 x_{1}+2 x_{2}-x_{3}=-1 \\ &6 x_{1}-6 x_{2}+7 x_{3}=-7 \\ &3 x_{1}-4 x_{2}+4 x_{3}=-6 \end{aligned} \nonumber$

To perform Gaussian elimination, we form an Augmented Matrix by combining the matrix $A$ with the column vector $b$ :

$\left(\begin{array}{rrrr} -3 & 2 & -1 & -1 \\ 6 & -6 & 7 & -7 \\ 3 & -4 & 4 & -6 \end{array}\right) \nonumber$

Row reduction is then performed on this matrix. Allowed operations are (1) multiply any row by a constant, (2) add multiple of one row to another row, (3) interchange the order of any rows. The goal is to convert the original matrix into an upper-triangular matrix

We start with the first row of the matrix and work our way down as follows. First we multiply the first row by 2 and add it to the second row, and add the first row to the third row:

$\left(\begin{array}{rrrr} -3 & 2 & -1 & -1 \\ 0 & -2 & 5 & -9 \\ 0 & -2 & 3 & -7 \end{array}\right) \nonumber$

We then go to the second row. We multiply this row by $-1$ and add it to the third row:

$\left(\begin{array}{rrrr} -3 & 2 & -1 & -1 \\ 0 & -2 & 5 & -9 \\ 0 & 0 & -2 & 2 \end{array}\right) \nonumber$

The resulting equations can be determined from the matrix and are given by

$\begin{aligned} -3 x_{1}+2 x_{2}-x_{3} &=-1 \\ -2 x_{2}+5 x_{3} &=-9 \\ -2 x_{3} &=2 \end{aligned} \nonumber$

These equations can be solved by backward substitution, starting from the last equation and working backwards. We have

$\begin{aligned} &-2 x_{3}=2 \rightarrow x_{3}=-1 \\ &-2 x_{2}=-9-5 x_{3}=-4 \rightarrow x_{2}=2 \\ &-3 x_{1}=-1-2 x_{2}+x_{3}=-6 \rightarrow x_{1}=2 \end{aligned} \nonumber$

Therefore,

$\left(\begin{array}{l} x_{1} \\ x_{2} \\ x_{3} \end{array}\right)=\left(\begin{array}{r} 2 \\ 2 \\ -1 \end{array}\right) \nonumber$

U\) decomposition

The process of Gaussian Elimination also results in the factoring of the matrix A to

$\mathrm{A}=\mathrm{LU}, \nonumber$

where $L$ is a lower triangular matrix and U is an upper triangular matrix. Using the same matrix A as in the last section, we show how this factorization is realized. We have

$\left(\begin{array}{rrr} -3 & 2 & -1 \\ 6 & -6 & 7 \\ 3 & -4 & 4 \end{array}\right) \rightarrow\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & -2 & 3 \end{array}\right)=\mathrm{M}_{1} \mathrm{~A} \nonumber$

where

$\mathrm{M}_{1} \mathrm{~A}=\left(\begin{array}{rrr} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 1 & 0 & 1 \end{array}\right)\left(\begin{array}{rrr} -3 & 2 & -1 \\ 6 & -6 & 7 \\ 3 & -4 & 4 \end{array}\right)=\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & -2 & 3 \end{array}\right) \nonumber$

Note that the matrix $\mathrm{M}_{1}$ performs row elimination on the first column. Two times the first row is added to the second row and one times the first row is added to the third row. The entries of the column of $\mathrm{M}_{1}$ come from $2=-(6 /-3)$ and $1=-(3 /-3)$ as required for row elimination. The number $-3$ is called the pivot.

The next step is

$\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & -2 & 3 \end{array}\right) \rightarrow\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & 0 & -2 \end{array}\right)=\mathrm{M}_{2}\left(\mathrm{M}_{1} \mathrm{~A}\right) \nonumber$

where

$\mathrm{M}_{2}\left(\mathrm{M}_{1} \mathrm{~A}\right)=\left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -1 & 1 \end{array}\right)\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & -2 & 3 \end{array}\right)=\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & 0 & -2 \end{array}\right) \nonumber$

Here, $\mathrm{M}_{2}$ multiplies the second row by $-1=-(-2 /-2)$ and adds it to the third row. The pivot is $-2$ .

We now have

$\mathrm{M}_{2} \mathrm{M}_{1} \mathrm{~A}=\mathrm{U} \nonumber$

$\mathrm{A}=\mathrm{M}_{1}^{-1} \mathrm{M}_{2}^{-1} \mathrm{U} \nonumber$

The inverse matrices are easy to find. The matrix $\mathrm{M}_{1}$ multiples the first row by 2 and adds it to the second row, and multiplies the first row by 1 and adds it to the third row. To invert these operations, we need to multiply the first row by $-2$ and add it to the second row, and multiply the first row by $-1$ and add it to the third row. To check, with

$\mathrm{M}_{1} \mathrm{M}_{1}^{-1}=\mathrm{I} \nonumber$

we have

$\left(\begin{array}{lll} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 1 & 0 & 1 \end{array}\right)\left(\begin{array}{rrr} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1 & 0 & 1 \end{array}\right)=\left(\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right) \nonumber$

Similarly,

$\mathrm{M}_{2}^{-1}=\left(\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \end{array}\right) \nonumber$

Therefore,

$\mathrm{L}=\mathrm{M}_{1}^{-1} \mathrm{M}_{2}^{-1} \nonumber$

is given by

$L=\left(\begin{array}{rll} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1 & 0 & 1 \end{array}\right)\left(\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \end{array}\right)=\left(\begin{array}{rrr} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1 & 1 & 1 \end{array}\right) \nonumber$

which is lower triangular. The off-diagonal elements of $\mathrm{M}_{1}^{-1}$ and $\mathrm{M}_{2}^{-1}$ are simply combined to form L. Our LU decomposition is therefore

$\left(\begin{array}{rrr} -3 & 2 & -1 \\ 6 & -6 & 7 \\ 3 & -4 & 4 \end{array}\right)=\left(\begin{array}{rrr} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1 & 1 & 1 \end{array}\right)\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & 0 & -2 \end{array}\right) . \nonumber$

Another nice feature of the LU decomposition is that it can be done by overwriting A, therefore saving memory if the matrix A is very large.

The LU decomposition is useful when one needs to solve $A \mathbf{x}=\mathbf{b}$ for $\mathbf{x}$ when A is fixed and there are many different b’s. First one determines $\mathrm{L}$ and $\mathrm{U}$ using Gaussian elimination. Then one writes

$(\mathrm{LU}) \mathbf{x}=\mathrm{L}(\mathrm{U} \mathbf{x})=\mathbf{b} . \nonumber$

We let

$\mathbf{y}=\mathbf{U} \mathbf{x} \nonumber$

and first solve

$\mathrm{Ly}=\mathbf{b} \nonumber$

for $y$ by forward substitution. We then solve

$\mathrm{Ux}=\mathbf{y} \nonumber$

for $\mathbf{x}$ by backward substitution. When we count operations, we will see that solving $(\mathrm{LU}) \mathbf{x}=\mathbf{b}$ is significantly faster once $\mathrm{L}$ and $\mathrm{U}$ are in hand than solving $\mathrm{A} \mathbf{x}=\mathbf{b}$ directly by Gaussian elimination.

We now illustrate the solution of $L U \mathbf{x}=\mathbf{b}$ using our previous example, where

$\mathrm{L}=\left(\begin{array}{rrr} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1 & 1 & 1 \end{array}\right), \quad \mathrm{U}=\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & 0 & -2 \end{array}\right), \quad \mathbf{b}=\left(\begin{array}{l} -1 \\ -7 \\ -6 \end{array}\right) \nonumber$

With $\mathbf{y}=\mathbf{U x}$ , we first solve $\mathbf{L} \mathbf{y}=\mathbf{b}$ , that is

$\left(\begin{array}{rrr} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1 & 1 & 1 \end{array}\right)\left(\begin{array}{l} y_{1} \\ y_{2} \\ y_{3} \end{array}\right)=\left(\begin{array}{l} -1 \\ -7 \\ -6 \end{array}\right) \nonumber$

Using forward substitution

$\begin{aligned} &y_{1}=-1 \\ &y_{2}=-7+2 y_{1}=-9 \\ &y_{3}=-6+y_{1}-y_{2}=2 \end{aligned} \nonumber$

We now solve $U x=y$ , that is

$\left(\begin{array}{rrr} -3 & 2 & -1 \\ 0 & -2 & 5 \\ 0 & 0 & -2 \end{array}\right)\left(\begin{array}{l} x_{1} \\ x_{2} \\ x_{3} \end{array}\right)=\left(\begin{array}{r} -1 \\ -9 \\ 2 \end{array}\right) \nonumber$

Using backward substitution,

$\begin{aligned} &-2 x_{3}=2 \rightarrow x_{3}=-1 \\ &-2 x_{2}=-9-5 x_{3}=-4 \rightarrow x_{2}=2 \\ &-3 x_{1}=-1-2 x_{2}+x_{3}=-6 \rightarrow x_{1}=2 \end{aligned} \nonumber$

and we have once again determined

$\left(\begin{array}{l} x_{1} \\ x_{2} \\ x_{3} \end{array}\right)=\left(\begin{array}{r} 2 \\ 2 \\ -1 \end{array}\right) \nonumber$

Partial pivoting

When performing Gaussian elimination, the diagonal element that one uses during the elimination procedure is called the pivot. To obtain the correct multiple, one uses the pivot as the divisor to the elements below the pivot. Gaussian elimination in this form will fail if the pivot is zero. In this situation, a row interchange must be performed.

Even if the pivot is not identically zero, a small value can result in big roundoff errors. For very large matrices, one can easily lose all accuracy in the solution. To avoid these round-off errors arising from small pivots, row interchanges are made, and this technique is called partial pivoting (partial pivoting is in contrast to complete pivoting, where both rows and columns are interchanged). We will illustrate by example the LU decomposition using partial pivoting. Consider

$A=\left(\begin{array}{rrr} -2 & 2 & -1 \\ 6 & -6 & 7 \\ 3 & -8 & 4 \end{array}\right) \nonumber$

We interchange rows to place the largest element (in absolute value) in the pivot, or $a_{11}$ , position. That is,

$\mathrm{A} \rightarrow\left(\begin{array}{rrr} 6 & -6 & 7 \\ -2 & 2 & -1 \\ 3 & -8 & 4 \end{array}\right)=\mathrm{P}_{12} \mathrm{~A} \nonumber$

where

$P_{12}=\left(\begin{array}{lll} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right) \nonumber$

is a permutation matrix that when multiplied on the left interchanges the first and second rows of a matrix. Note that $P_{12}^{-1}=P_{12}$ . The elimination step is then

$\mathrm{P}_{12} \mathrm{~A} \rightarrow\left(\begin{array}{rrc} 6 & -6 & 7 \\ 0 & 0 & 4 / 3 \\ 0 & -5 & 1 / 2 \end{array}\right)=\mathrm{M}_{1} \mathrm{P}_{12} \mathrm{~A} \nonumber$

where

$M_{1}=\left(\begin{array}{ccc} 1 & 0 & 0 \\ 1 / 3 & 1 & 0 \\ -1 / 2 & 0 & 1 \end{array}\right) \nonumber$

The final step requires one more row interchange:

$\mathrm{M}_{1} \mathrm{P}_{12} \mathrm{~A} \rightarrow\left(\begin{array}{rrc} 6 & -6 & 7 \\ 0 & -5 & 1 / 2 \\ 0 & 0 & 4 / 3 \end{array}\right)=\mathrm{P}_{23} \mathrm{M}_{1} \mathrm{P}_{12} \mathrm{~A}=\mathrm{U} \nonumber$

Since the permutation matrices given by P are their own inverses, we can write our result as

$\left(P_{23} M_{1} P_{23}\right) P_{23} P_{12} A=U \nonumber$

Multiplication of $M$ on the left by $P$ interchanges rows while multiplication on the right by $P$ interchanges columns. That is,

$P_{23}\left(\begin{array}{cll} 1 & 0 & 0 \\ 1 / 3 & 1 & 0 \\ -1 / 2 & 0 & 1 \end{array}\right) P_{23}=\left(\begin{array}{ccc} 1 & 0 & 0 \\ -1 / 2 & 0 & 1 \\ 1 / 3 & 1 & 0 \end{array}\right) P_{23}=\left(\begin{array}{ccc} 1 & 0 & 0 \\ -1 / 2 & 1 & 0 \\ 1 / 3 & 0 & 1 \end{array}\right) \nonumber$

The net result on $M_{1}$ is an interchange of the nondiagonal elements $1 / 3$ and $-1 / 2$ .

We can then multiply by the inverse of $\left(\mathrm{P}_{23} \mathrm{M}_{1} \mathrm{P}_{23}\right)$ to obtain

$P_{23} P_{12} A=\left(P_{23} M_{1} P_{23}\right)^{-1} U \nonumber$

which we write as

$\mathrm{PA}=\mathrm{LU} \nonumber$

Instead of L, MATLAB will write this as

$\mathrm{A}=\left(\mathrm{P}^{-1} \mathrm{~L}\right) \mathrm{U} \nonumber$

For convenience, we will just denote $\left(\mathrm{P}^{-1} \mathrm{~L}\right)$ by $\mathrm{L}$ , but by $\mathrm{L}$ here we mean a permutated lower triangular matrix.

For example, in MATLAB, to solve $A \mathbf{x}=\mathbf{b}$ for $\mathbf{x}$ using Gaussian elimination, one types

$\mathrm{x}=\mathrm{A} \backslash \mathrm{b} \nonumber$

where \solves for $\mathbf{x}$ using the most efficient algorithm available, depending on the form of $\mathrm{A}$ . If $\mathrm{A}$ is a general $n \times n$ matrix, then first the LU decomposition of $\mathrm{A}$ is found using partial pivoting, and then $x$ is determined from permuted forward and backward substitution. If $A$ is upper or lower triangular, then forward or backward substitution (or their permuted version) is used directly.

If there are many different right-hand-sides, one would first directly find the LU decomposition of A using a function call, and then solve using . That is, one would iterate for different $\mathbf{b}^{\prime}$ s the following expressions:

$\begin{aligned} &{[\mathrm{LU}]=l u(\mathrm{~A})} \\ &\mathrm{y}=\mathrm{L} \backslash \mathrm{b} \\ &\mathrm{x}=\mathrm{U} \backslash y \end{aligned} \nonumber$

where the second and third lines can be shortened to

$\mathrm{x}=\mathrm{U} \backslash(\mathrm{L} \backslash \mathrm{b}) \nonumber$

where the parenthesis are required. In lecture, I will demonstrate these solutions in MATLAB using the matrix $A=[-2,2,-1 ; 6,-6,7 ; 3,-8,4] ;$ which is the example in the notes.

$3.4$ Operation counts

To estimate how much computational time is required for an algorithm, one can count the number of operations required (multiplications, divisions, additions and subtractions). Usually, what is of interest is how the algorithm scales with the size of the problem. For example, suppose one wants to multiply two full $n \times n$ matrices. The calculation of each element requires $n$ multiplications and $n-1$ additions, or say $2 n-1$ operations. There are $n^{2}$ elements to compute so that the total operation count is $n^{2}(2 n-1)$ . If $n$ is large, we might want to know what will happen to the computational time if $n$ is doubled. What matters most is the fastest-growing, leading-order term in the operation count. In this matrix multiplication example, the operation count is $n^{2}(2 n-1)=2 n^{3}-n^{2}$ , and the leading-order term is $2 n^{3}$ . The factor of 2 is unimportant for the scaling, and we say that the algorithm scales like $\mathrm{O}\left(n^{3}\right)$ , which is read as "big Oh of n cubed." When using the big-Oh notation, we will drop both lower-order terms and constant multipliers.

The big-Oh notation tells us how the computational time of an algorithm scales. For example, suppose that the multiplication of two large $n \times n$ matrices took a computational time of $T_{n}$ seconds. With the known operation count going like $\mathrm{O}\left(n^{3}\right)$ , we can write

$T_{n}=k n^{3} \nonumber$

for some unknown constant $k$ . To determine how much longer the multiplication of two $2 n \times 2 n$ matrices will take, we write

$\begin{aligned} T_{2 n} &=k(2 n)^{3} \\ &=8 k n^{3} \\ &=8 T_{n} \end{aligned} \nonumber$

so that doubling the size of the matrix is expected to increase the computational time by a factor of $2^{3}=8$ .

Running MATLAB on my computer, the multiplication of two $2048 \times 2048$ matrices took about $0.75 \mathrm{sec}$ . The multiplication of two $4096 \times 4096$ matrices took about $6 \mathrm{sec}$ , which is 8 times longer. Timing of code in MATLAB can be found using the built-in stopwatch functions tic and toc.

What is the operation count and therefore the scaling of Gaussian elimination? Consider an elimination step with the pivot in the $i$ th row and $i$ th column. There are both $n-i$ rows below the pivot and $n-i$ columns to the right of the pivot. To perform elimination of one row, each matrix element to the right of the pivot must be multiplied by a factor and added to the row underneath. This must be done for all the rows. There are therefore $(n-i)(n-i)$ multiplication-additions to be done for this pivot. Since we are interested in only the scaling of the algorithm, I will just count a multiplication-addition as one operation.

To find the total operation count, we need to perform elimination using $n-1$ pivots, so that

$\begin{aligned} \text { op. counts } &=\sum_{i=1}^{n-1}(n-i)^{2} \\ &=(n-1)^{2}+(n-2)^{2}+\ldots(1)^{2} \\ &=\sum_{i=1}^{n-1} i^{2} \end{aligned} \nonumber$

Three summation formulas will come in handy. They are

$\begin{aligned} &\sum_{i=1}^{n} 1=n \\ &\sum_{i=1}^{n} i=\frac{1}{2} n(n+1) \\ &\sum_{i=1}^{n} i^{2}=\frac{1}{6} n(2 n+1)(n+1) \end{aligned} \nonumber$

which can be proved by mathematical induction, or derived by some tricks.

The operation count for Gaussian elimination is therefore

$\begin{aligned} \text { op. counts } &=\sum_{i=1}^{n-1} i^{2} \\ &=\frac{1}{6}(n-1)(2 n-1)(n) \end{aligned} \nonumber$

The leading-order term is therefore $n^{3} / 3$ , and we say that Gaussian elimination scales like $\mathrm{O}\left(n^{3}\right)$ . Since LU decomposition with partial pivoting is essentially Gaussian elimination, that algorithm also scales like $\mathrm{O}\left(n^{3}\right)$ . However, once the LU decomposition of a matrix A is known, the solution of $\mathrm{A} \mathbf{x}=\mathbf{b}$ can proceed by a forward and backward substitution. How does a backward substitution, say, scale? For backward substitution, the matrix equation to be solved is of the form

$\left(\begin{array}{ccccc} a_{1,1} & a_{1,2} & \cdots & a_{1, n-1} & a_{1, n} \\ 0 & a_{2,2} & \cdots & a_{2, n-1} & a_{2, n} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & a_{n-1, n-1} & a_{n-1, n} \\ 0 & 0 & \cdots & 0 & a_{n, n} \end{array}\right)\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{n-1} \\ x_{n} \end{array}\right)=\left(\begin{array}{c} b_{1} \\ b_{2} \\ \vdots \\ b_{n-1} \\ b_{n} \end{array}\right) \nonumber$

The solution for $x_{i}$ is found after solving for $x_{j}$ with $j>i$ . The explicit solution for $x_{i}$ is given by

$x_{i}=\frac{1}{a_{i, i}}\left(b_{i}-\sum_{j=i+1}^{n} a_{i, j} x_{j}\right) \nonumber$

The solution for $x_{i}$ requires $n-i+1$ multiplication-additions, and since this must be done for $n$ such $i^{\prime}$ s, we have

$\begin{aligned} \text { op. counts } &=\sum_{i=1}^{n} n-i+1 \\ &=n+(n-1)+\cdots+1 \\ &=\sum_{i=1}^{n} i \\ &=\frac{1}{2} n(n+1) \end{aligned} \nonumber$

The leading-order term is $n^{2} / 2$ and the scaling of backward substitution is $\mathrm{O}\left(n^{2}\right)$ . After the LU decomposition of a matrix A is found, only a single forward and backward substitution is required to solve $A \mathbf{x}=\mathbf{b}$ , and the scaling of the algorithm to solve this matrix equation is therefore still $\mathrm{O}\left(n^{2}\right)$ . For large $n$ , one should expect that solving $\mathbf{A x}=\mathbf{b}$ by a forward and backward substitution should be substantially faster than a direct solution using Gaussian elimination.

$3.5$ System of nonlinear equations

A system of nonlinear equations can be solved using a version of Newton’s Method. We illustrate this method for a system of two equations and two unknowns. Suppose that we want to solve

$f(x, y)=0, \quad g(x, y)=0 \nonumber$

for the unknowns $x$ and $y .$ We want to construct two simultaneous sequences $x_{0}, x_{1}, x_{2}, \ldots$ and $y_{0}, y_{1}, y_{2}, \ldots$ that converge to the roots. To construct these sequences, we Taylor series expand $f\left(x_{n+1}, y_{n+1}\right)$ and $g\left(x_{n+1}, y_{n+1}\right)$ about the point $\left(x_{n}, y_{n}\right) .$ Using the partial derivatives $f_{x}=\partial f / \partial x, f_{y}=\partial f / \partial y$ , etc., the twodimensional Taylor series expansions, displaying only the linear terms, are given by

$\begin{aligned} f\left(x_{n+1}, y_{n+1}\right)=f\left(x_{n}, y_{n}\right)+\left(x_{n+1}-x_{n}\right) f_{x}\left(x_{n}, y_{n}\right) & \\ &+\left(y_{n+1}-y_{n}\right) f_{y}\left(x_{n}, y_{n}\right)+\ldots \end{aligned} \nonumber$

$g\left(x_{n+1}, y_{n+1}\right)=g\left(x_{n}, y_{n}\right)+\left(x_{n+1}-x_{n}\right) g_{x}\left(x_{n}, y_{n}\right)$

$+\left(y_{n+1}-y_{n}\right) g_{y}\left(x_{n}, y_{n}\right)+\ldots \nonumber$

To obtain Newton’s method, we take $f\left(x_{n+1}, y_{n+1}\right)=0, g\left(x_{n+1}, y_{n+1}\right)=0$ and drop higher-order terms above linear. Although one can then find a system of linear equations for $x_{n+1}$ and $y_{n+1}$ , it is more convenient to define the variables

$\Delta x_{n}=x_{n+1}-x_{n}, \quad \Delta y_{n}=y_{n+1}-y_{n} . \nonumber$

The iteration equations will then be given by

$x_{n+1}=x_{n}+\Delta x_{n}, \quad y_{n+1}=y_{n}+\Delta y_{n} \nonumber$

and the linear equations to be solved for $\Delta x_{n}$ and $\Delta y_{n}$ are given by

$\left(\begin{array}{ll} f_{x} & f_{y} \\ g_{x} & g_{y} \end{array}\right)\left(\begin{array}{l} \Delta x_{n} \\ \Delta y_{n} \end{array}\right)=\left(\begin{array}{l} -f \\ -g \end{array}\right) \nonumber$

where $f, g, f_{x}, f_{y}, g_{x}$ , and $g_{y}$ are all evaluated at the point $\left(x_{n}, y_{n}\right) .$ The twodimensional case is easily generalized to $n$ dimensions. The matrix of partial derivatives is called the Jacobian Matrix.

We illustrate Newton’s Method by finding the steady state solution of the Lorenz equations, given by

$\begin{array}{r} \sigma(y-x)=0 \\ r x-y-x z=0 \\ x y-b z=0 \end{array} \nonumber$

where $x, y$ , and $z$ are the unknown variables and $\sigma, r$ , and $b$ are the known parameters. Here, we have a three-dimensional homogeneous system $f=0, g=0$ , and $h=0$ , with

$\begin{aligned} &f(x, y, z)=\sigma(y-x) \\ &g(x, y, z)=r x-y-x z \\ &h(x, y, z)=x y-b z . \end{aligned} \nonumber$

The partial derivatives can be computed to be

$\begin{array}{lll} f_{x}=-\sigma, & f_{y}=\sigma, & f_{z}=0, \\ g_{x}=r-z, & g_{y}=-1, & g_{z}=-x, \\ h_{x}=y, & h_{y}=x, & h_{z}=-b . \end{array} \nonumber$

The iteration equation is therefore

$\left(\begin{array}{ccc} -\sigma & \sigma & 0 \\ r-z_{n} & -1 & -x_{n} \\ y_{n} & x_{n} & -b \end{array}\right)\left(\begin{array}{c} \Delta x_{n} \\ \Delta y_{n} \\ \Delta z_{n} \end{array}\right)=-\left(\begin{array}{c} \sigma\left(y_{n}-x_{n}\right) \\ r x_{n}-y_{n}-x_{n} z_{n} \\ x_{n} y_{n}-b z_{n} \end{array}\right) \nonumber$

with

$\begin{aligned} &x_{n+1}=x_{n}+\Delta x_{n} \\ &y_{n+1}=y_{n}+\Delta y_{n} \\ &z_{n+1}=z_{n}+\Delta z_{n} \end{aligned} \nonumber$

The MATLAB program that solves this system is contained in newton_system.m.

Chapter 4 Least-squares approximation

The method of least-squares is commonly used to fit a parameterized curve to experimental data. In general, the fitting curve is not expected to pass through the data points, making this problem substantially different from the one of interpolation.

We consider here only the simplest case of the same experimental error for all the data points. Let the data to be fitted be given by $\left(x_{i}, y_{i}\right)$ , with $i=1$ to $n$ .

Fitting a straight line

Suppose the fitting curve is a line. We write for the fitting curve

$y(x)=\alpha x+\beta . \nonumber$

The distance $r_{i}$ from the data point $\left(x_{i}, y_{i}\right)$ and the fitting curve is given by

$\begin{aligned} r_{i} &=y_{i}-y\left(x_{i}\right) \\ &=y_{i}-\left(\alpha x_{i}+\beta\right) . \end{aligned} \nonumber$

A least-squares fit minimizes the sum of the squares of the $r_{i}$ ’s. This minimum can be shown to result in the most probable values of $\alpha$ and $\beta$ .

We define

$\begin{aligned} \rho &=\sum_{i=1}^{n} r_{i}^{2} \\ &=\sum_{i=1}^{n}\left(y_{i}-\left(\alpha x_{i}+\beta\right)\right)^{2} \end{aligned} \nonumber$

To minimize $\rho$ with respect to $\alpha$ and $\beta$ , we solve

$\frac{\partial \rho}{\partial \alpha}=0, \quad \frac{\partial \rho}{\partial \beta}=0 \nonumber$

Taking the partial derivatives, we have

$\begin{aligned} &\frac{\partial \rho}{\partial \alpha}=\sum_{i=1}^{n} 2\left(-x_{i}\right)\left(y_{i}-\left(\alpha x_{i}+\beta\right)\right)=0 \\ &\frac{\partial \rho}{\partial \beta}=\sum_{i=1}^{n} 2(-1)\left(y_{i}-\left(\alpha x_{i}+\beta\right)\right)=0 \end{aligned} \nonumber$

These equations form a system of two linear equations in the two unknowns $\alpha$ and $\beta$ , which is evident when rewritten in the form

$\begin{aligned} \alpha \sum_{i=1}^{n} x_{i}^{2}+\beta \sum_{i=1}^{n} x_{i} &=\sum_{i=1}^{n} x_{i} y_{i} \\ \alpha \sum_{i=1}^{n} x_{i}+\beta n &=\sum_{i=1}^{n} y_{i} \end{aligned} \nonumber$

These equations can be solved either analytically, or numerically in MATLAB, where the matrix form is

$\left(\begin{array}{ccc} \sum_{i=1}^{n} & x_{i}^{2} & \sum_{i=1}^{n} \\ \sum_{i=1}^{n} & x_{i} & n \end{array}\right)\left(\begin{array}{c} \alpha \\ \beta \end{array}\right)=\left(\begin{array}{c} \sum_{i=1}^{n} x_{i} y_{i} \\ \sum_{i=1}^{n} y_{i} \end{array}\right) . \nonumber$

A proper statistical treatment of this problem should also consider an estimate of the errors in $\alpha$ and $\beta$ as well as an estimate of the goodness-of-fit of the data to the model. We leave these further considerations to a statistics class.

Fitting to a linear combination of functions

Consider the general fitting function

$y(x)=\sum_{j=1}^{m} c_{j} f_{j}(x) \nonumber$

where we assume $m$ functions $f_{j}(x)$ . For example, if we want to fit a cubic polynomial to the data, then we would have $m=4$ and take $f_{1}=1, f_{2}=x, f_{3}=x^{2}$ and $f_{4}=x^{3}$ . Typically, the number of functions $f_{j}$ is less than the number of data points; that is, $m<n$ , so that a direct attempt to solve for the $c_{j}^{\prime}$ s such that the fitting function exactly passes through the $n$ data points would result in $n$ equations and $m$ unknowns. This would be an over-determined linear system that in general has no solution.

We now define the vectors

$\mathbf{y}=\left(\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array}\right), \quad \mathbf{c}=\left(\begin{array}{c} c_{1} \\ c_{2} \\ \vdots \\ c_{m} \end{array}\right) \nonumber$

and the $n$ -by-m matrix

$\mathrm{A}=\left(\begin{array}{cccc} f_{1}\left(x_{1}\right) & f_{2}\left(x_{1}\right) & \cdots & f_{m}\left(x_{1}\right) \\ f_{1}\left(x_{2}\right) & f_{2}\left(x_{2}\right) & \cdots & f_{m}\left(x_{2}\right) \\ \vdots & \vdots & \ddots & \vdots \\ f_{1}\left(x_{n}\right) & f_{2}\left(x_{n}\right) & \cdots & f_{m}\left(x_{n}\right) \end{array}\right) \nonumber$

With $m<n$ , the equation $A \mathbf{c}=\mathbf{y}$ is over-determined. We let

$\mathbf{r}=\mathbf{y}-\mathrm{A} \mathbf{c} \nonumber$

be the residual vector, and let

$\rho=\sum_{i=1}^{n} r_{i}^{2} \nonumber$

The method of least squares minimizes $\rho$ with respect to the components of c. Now, using $T$ to signify the transpose of a matrix, we have

$\begin{aligned} \rho &=\mathbf{r}^{T} \mathbf{r} \\ &=(\mathbf{y}-\mathrm{A} \mathbf{c})^{T}(\mathbf{y}-\mathrm{A} \mathbf{c}) \\ &=\mathbf{y}^{T} \mathbf{y}-\mathbf{c}^{T} \mathrm{~A}^{T} \mathbf{y}-\mathbf{y}^{T} \mathrm{~A} \mathbf{c}+\mathbf{c}^{T} \mathrm{~A}^{T} \mathrm{~A} \mathbf{c} . \end{aligned} \nonumber$

Since $\rho$ is a scalar, each term in the above expression must be a scalar, and since the transpose of a scalar is equal to the scalar, we have

$\mathbf{c}^{T} \mathrm{~A}^{T} \mathbf{y}=\left(\mathbf{c}^{T} \mathrm{~A}^{T} \mathbf{y}\right)^{T}=\mathbf{y}^{T} \mathrm{~A} \mathbf{c} . \nonumber$

Therefore,

$\rho=\mathbf{y}^{T} \mathbf{y}-2 \mathbf{y}^{T} \mathrm{~A} \mathbf{c}+\mathbf{c}^{T} \mathrm{~A}^{T} \mathrm{~A} \mathbf{c} . \nonumber$

To find the minimum of $\rho$ , we will need to solve $\partial \rho / \partial c_{j}=0$ for $j=1, \ldots, m$ . To take the derivative of $\rho$ , we switch to a tensor notation, using the Einstein summation convention, where repeated indices are summed over their allowable range. We can write

$\rho=y_{i} y_{i}-2 y_{i} \mathrm{~A}_{i k} c_{k}+c_{i} \mathrm{~A}_{i k}^{T} \mathrm{~A}_{k l} c_{l} \nonumber$

Taking the partial derivative, we have

$\frac{\partial \rho}{\partial c_{j}}=-2 y_{i} \mathrm{~A}_{i k} \frac{\partial c_{k}}{\partial c_{j}}+\frac{\partial c_{i}}{\partial c_{j}} \mathrm{~A}_{i k}^{T} \mathrm{~A}_{k l} c_{l}+c_{i} \mathrm{~A}_{i k}^{T} \mathrm{~A}_{k l} \frac{\partial c_{l}}{\partial c_{j}} \nonumber$

Now,

$\frac{\partial c_{i}}{\partial c_{j}}= \begin{cases}1, & \text { if } i=j \\ 0, & \text { otherwise }\end{cases} \nonumber$

Therefore,

$\frac{\partial \rho}{\partial c_{j}}=-2 y_{i} \mathrm{~A}_{i j}+\mathrm{A}_{j k}^{T} \mathrm{~A}_{k l} c_{l}+c_{i} \mathrm{~A}_{i k}^{T} \mathrm{~A}_{k j} \nonumber$

Now,

$\begin{aligned} c_{i} \mathrm{~A}_{i k}^{T} \mathrm{~A}_{k j} &=c_{i} \mathrm{~A}_{k i} \mathrm{~A}_{k j} \\ &=\mathrm{A}_{k j} \mathrm{~A}_{k i} c_{i} \\ &=\mathrm{A}_{j k}^{T} \mathrm{~A}_{k i} c_{i} \\ &=\mathrm{A}_{j k}^{T} \mathrm{~A}_{k l} c_{l} \end{aligned} \nonumber$

Therefore,

$\frac{\partial \rho}{\partial c_{j}}=-2 y_{i} \mathrm{~A}_{i j}+2 \mathrm{~A}_{j k}^{T} \mathrm{~A}_{k l} c_{l} \nonumber$

With the partials set equal to zero, we have

$\mathrm{A}_{j k}^{T} \mathrm{~A}_{k l} c_{l}=y_{i} \mathrm{~A}_{i j} \nonumber$

$\mathrm{A}_{j k}^{T} \mathrm{~A}_{k l} c_{l}=\mathrm{A}_{j i}^{T} y_{i} \nonumber$

In vector notation, we have

$\mathrm{A}^{T} \mathrm{~A} \mathbf{c}=\mathrm{A}^{T} \mathbf{y} . \nonumber$

Equation (4.2) is the so-called normal equation, and can be solved for c by Gaussian elimination using the MATLAB backslash operator. After constructing the matrix A given by (4.1), and the vector y from the data, one can code in MATLAB

$c=\left(A^{\prime} A\right) \backslash\left(A^{\prime} y\right) \nonumber$

But in fact the MATLAB back slash operator will automatically solve the normal equations when the matrix $A$ is not square, so that the MATLAB code

$c=A \backslash y \nonumber$

yields the same result.

Chapter 5 Interpolation

Consider the following problem: Given the values of a known function $y=f(x)$ at a sequence of ordered points $x_{0}, x_{1}, \ldots, x_{n}$ , find $f(x)$ for arbitrary $x .$ When $x_{0} \leq$ $x \leq x_{n}$ , the problem is called interpolation. When $x<x_{0}$ or $x>x_{n}$ the problem is called extrapolation.

With $y_{i}=f\left(x_{i}\right)$ , the problem of interpolation is basically one of drawing a smooth curve through the known points $\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots,\left(x_{n}, y_{n}\right)$ . This is not the same problem as drawing a smooth curve that approximates a set of data points that have experimental error. This latter problem is called least-squares approximation.

Here, we will consider three interpolation algorithms: (1) polynomial interpolation; (2) piecewise linear interpolation, and; (3) cubic spline interpolation.

Polynomial interpolation

The $n+1$ points $\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots,\left(x_{n}, y_{n}\right)$ can be interpolated by a unique polynomial of degree $n .$ When $n=1$ , the polynomial is a linear function; when $n=2$ , the polynomial is a quadratic function. There are three standard algorithms that can be used to construct this unique interpolating polynomial, and we will present all three here, not so much because they are all useful, but because it is interesting to learn how these three algorithms are constructed.

When discussing each algorithm, we define $P_{n}(x)$ to be the unique $n$ th degree polynomial that passes through the given $n+1$ data points.

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.

Piecewise linear interpolation

Instead of constructing a single global polynomial that goes through all the points, one can construct local polynomials that are then connected together. In the the section following this one, we will discuss how this may be done using cubic polynomials. Here, we discuss the simpler case of linear polynomials. This is the default interpolation typically used when plotting data.

Suppose the interpolating function is $y=g(x)$ , and as previously, there are $n+1$ points to interpolate. We construct the function $g(x)$ out of $n$ local linear polynomials. We write

$g(x)=g_{i}(x), \quad \text { for } x_{i} \leq x \leq x_{i+1} \nonumber$

where

$g_{i}(x)=a_{i}\left(x-x_{i}\right)+b_{i} \nonumber$

and $i=0,1, \ldots, n-1$ .

We now require $y=g_{i}(x)$ to pass through the endpoints $\left(x_{i}, y_{i}\right)$ and $\left(x_{i+1}, y_{i+1}\right)$ . We have

$\begin{aligned} y_{i} &=b_{i} \\ y_{i+1} &=a_{i}\left(x_{i+1}-x_{i}\right)+b_{i} . \end{aligned} \nonumber$

Therefore, the coefficients of $g_{i}(x)$ are determined to be

$a_{i}=\frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}}, \quad b_{i}=y_{i} \nonumber$

Although piecewise linear interpolation is widely used, particularly in plotting routines, it suffers from a discontinuity in the derivative at each point. This results in a function which may not look smooth if the points are too widely spaced. We next consider a more challenging algorithm that uses cubic polynomials.

Cubic spline interpolation

The $n+1$ points to be interpolated are again

$\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots\left(x_{n}, y_{n}\right) \nonumber$

Here, we use $n$ piecewise cubic polynomials for interpolation,

$g_{i}(x)=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i}, \quad i=0,1, \ldots, n-1, \nonumber$

with the global interpolation function written as

$g(x)=g_{i}(x), \quad \text { for } x_{i} \leq x \leq x_{i+1} . \nonumber$

To achieve a smooth interpolation we impose that $g(x)$ and its first and second derivatives are continuous. The requirement that $g(x)$ is continuous (and goes through all $n+1$ points) results in the two constraints

$\begin{gathered} g_{i}\left(x_{i}\right)=y_{i}, \quad i=0 \text { to } n-1 \\ g_{i}\left(x_{i+1}\right)=y_{i+1}, \quad i=0 \text { to } n-1 \end{gathered} \nonumber$

The requirement that $g^{\prime}(x)$ is continuous results in

$g_{i}^{\prime}\left(x_{i+1}\right)=g_{i+1}^{\prime}\left(x_{i+1}\right), \quad i=0 \text { to } n-2 \nonumber$

And the requirement that $g^{\prime \prime}(x)$ is continuous results in

$g_{i}^{\prime \prime}\left(x_{i+1}\right)=g_{i+1}^{\prime \prime}\left(x_{i+1}\right), \quad i=0 \text { to } n-2 \nonumber$

There are $n$ cubic polynomials $g_{i}(x)$ and each cubic polynomial has four free coefficients; there are therefore a total of $4 n$ unknown coefficients. The number of constraining equations from (5.1)-(5.4) is $2 n+2(n-1)=4 n-2$ . With $4 n-2$ constraints and $4 n$ unknowns, two more conditions are required for a unique solution. These are usually chosen to be extra conditions on the first $g_{0}(x)$ and last $g_{n-1}(x)$ polynomials. We will discuss these extra conditions later.

We now proceed to determine equations for the unknown coefficients of the cubic polynomials. The polynomials and their first two derivatives are given by

$\begin{aligned} g_{i}(x) &=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i} \\ g_{i}^{\prime}(x) &=3 a_{i}\left(x-x_{i}\right)^{2}+2 b_{i}\left(x-x_{i}\right)+c_{i} \\ g_{i}^{\prime \prime}(x) &=6 a_{i}\left(x-x_{i}\right)+2 b_{i} \end{aligned} \nonumber$

We will consider the four conditions (5.1)-(5.4) in turn. From (5.1) and (5.5), we have

$d_{i}=y_{i}, \quad i=0 \text { to } n-1, \nonumber$

which directly solves for all of the $d$ -coefficients.

To satisfy $(5.2)$ , we first define

$h_{i}=x_{i+1}-x_{i} \nonumber$

and

$f_{i}=y_{i+1}-y_{i} \nonumber$

Now, from (5.2) and (5.5), using (5.8), we obtain the $n$ equations

$a_{i} h_{i}^{3}+b_{i} h_{i}^{2}+c_{i} h_{i}=f_{i}, \quad i=0 \text { to } n-1 \nonumber$

From (5.3) and (5.6) we obtain the $n-1$ equations

$3 a_{i} h_{i}^{2}+2 b_{i} h_{i}+c_{i}=c_{i+1}, \quad i=0 \text { to } n-2 \nonumber$

From (5.4) and (5.7) we obtain the $n-1$ equations

$3 a_{i} h_{i}+b_{i}=b_{i+1} \quad i=0 \text { to } n-2 \nonumber$

It is will be useful to include a definition of the coefficient $b_{n}$ , which is now missing. (The index of the cubic polynomial coefficients only go up to $n-1 .$ ) We simply extend (5.11) up to $i=n-1$ and so write

$3 a_{n-1} h_{n-1}+b_{n-1}=b_{n} \nonumber$

which can be viewed as a definition of $b_{n}$ .

We now proceed to eliminate the sets of $a$ - and $c$ -coefficients (with the $d$ -coefficients already eliminated in (5.8)) to find a system of linear equations for the $b$ -coefficients. From (5.11) and (5.12), we can solve for the $n a$ -coefficients to find

$a_{i}=\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right), \quad i=0 \text { to } n-1 \nonumber$

From (5.9), we can solve for the $n$ c-coefficients as follows:

$\begin{aligned} c_{i} &=\frac{1}{h_{i}}\left(f_{i}-a_{i} h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ &=\frac{1}{h_{i}}\left(f_{i}-\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right) h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ &=\frac{f_{i}}{h_{i}}-\frac{1}{3} h_{i}\left(b_{i+1}+2 b_{i}\right), \quad i=0 \text { to } n-1 . \end{aligned} \nonumber$

We can now find an equation for the $b$ -coefficients by substituting (5.8), (5.13) and $(5.14)$ into $(5.10)$ :

$\begin{gathered} 3\left(\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right)\right) h_{i}^{2}+2 b_{i} h_{i}+\left(\frac{f_{i}}{h_{i}}-\frac{1}{3} h_{i}\left(b_{i+1}+2 b_{i}\right)\right) \\ =\left(\frac{f_{i+1}}{h_{i+1}}-\frac{1}{3} h_{i+1}\left(b_{i+2}+2 b_{i+1}\right)\right) \end{gathered} \nonumber$

which simplifies to

$\frac{1}{3} h_{i} b_{i}+\frac{2}{3}\left(h_{i}+h_{i+1}\right) b_{i+1}+\frac{1}{3} h_{i+1} b_{i+2}=\frac{f_{i+1}}{h_{i+1}}-\frac{f_{i}}{h_{i}} \nonumber$

an equation that is valid for $i=0$ to $n-2$ . Therefore, (5.15) represent $n-1$ equations for the $n+1$ unknown $b$ -coefficients. Accordingly, we write the matrix equa- tion for the $b$ -coefficients, leaving the first and last row absent, as

$\left(\begin{array}{cccccccc} \cdots & \cdots & \cdots & \cdots & \text { missing } & \cdots & \cdots \\ \frac{1}{3} h_{0} & \frac{2}{3}\left(h_{0}+h_{1}\right) & \frac{1}{3} h_{1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & \frac{1}{3} h_{n-2} & \frac{2}{3}\left(h_{n-2}+h_{n-1}\right) & \frac{1}{3} h_{n-1} \\ \cdots & \cdots & \cdots & \cdots & \text { missing } & \cdots & \cdots \end{array}\right)\left(\begin{array}{c} b_{0} \\ b_{1} \\ \vdots \\ b_{n-1} \\ b_{n} \end{array}\right)\left(\begin{array}{c} \operatorname{missing} \\ \frac{f_{1}}{h_{1}}-\frac{f_{0}}{h_{0}} \\ \vdots \\ \frac{f_{n-1}}{h_{n-1}}-\frac{f_{n-2}}{h_{n-2}} \\ \operatorname{missing} \end{array}\right) \nonumber$

Once the missing first and last equations are specified, the matrix equation for the $b$ -coefficients can be solved by Gaussian elimination. And once the $b$ -coefficients are determined, the $a$ - and $c$ -coefficients can also be determined from (5.13) and (5.14), with the $d$ -coefficients already known from (5.8). The piecewise cubic polynomials, then, are known and $g(x)$ can be used for interpolation to any value $x$ satisfying $x_{0} \leq x \leq x_{n}$

The missing first and last equations can be specified in several ways, and here we show the two ways that are allowed by the MATLAB function spline.m. The first way should be used when the derivative $g^{\prime}(x)$ is known at the endpoints $x_{0}$ and $x_{n}$ ; that is, suppose we know the values of $\alpha$ and $\beta$ such that

$g_{0}^{\prime}\left(x_{0}\right)=\alpha, \quad g_{n-1}^{\prime}\left(x_{n}\right)=\beta . \nonumber$

From the known value of $\alpha$ , and using $(5.6)$ and $(5.14)$ , we have

$\begin{aligned} \alpha &=c_{0} \\ &=\frac{f_{0}}{h_{0}}-\frac{1}{3} h_{0}\left(b_{1}+2 b_{0}\right) \end{aligned} \nonumber$

Therefore, the missing first equation is determined to be

$\frac{2}{3} h_{0} b_{0}+\frac{1}{3} h_{0} b_{1}=\frac{f_{0}}{h_{0}}-\alpha . \nonumber$

From the known value of $\beta$ , and using $(5.6),(5.13)$ , and $(5.14)$ , we have

$\begin{aligned} \beta &=3 a_{n-1} h_{n-1}^{2}+2 b_{n-1} h_{n-1}+c_{n-1} \\ &=3\left(\frac{1}{3 h_{n-1}}\left(b_{n}-b_{n-1}\right)\right) h_{n-1}^{2}+2 b_{n-1} h_{n-1}+\left(\frac{f_{n-1}}{h_{n-1}}-\frac{1}{3} h_{n-1}\left(b_{n}+2 b_{n-1}\right)\right), \end{aligned} \nonumber$

which simplifies to

$\frac{1}{3} h_{n-1} b_{n-1}+\frac{2}{3} h_{n-1} b_{n}=\beta-\frac{f_{n-1}}{h_{n-1}} \nonumber$

to be used as the missing last equation.

The second way of specifying the missing first and last equations is called the not- $a$ -knot condition, which assumes that

$g_{0}(x)=g_{1}(x), \quad g_{n-2}(x)=g_{n-1}(x) \nonumber$

Considering the first of these equations, from (5.5) we have

$\begin{aligned} a_{0}\left(x-x_{0}\right)^{3}+b_{0}\left(x-x_{0}\right)^{2}+c_{0}(x&\left.-x_{0}\right)+d_{0} \\ &=a_{1}\left(x-x_{1}\right)^{3}+b_{1}\left(x-x_{1}\right)^{2}+c_{1}\left(x-x_{1}\right)+d_{1} \end{aligned} \nonumber$

Now two cubic polynomials can be proven to be identical if at some value of $x$ , the polynomials and their first three derivatives are identical. Our conditions of continuity at $x=x_{1}$ already require that at this value of $x$ these two polynomials and their first two derivatives are identical. The polynomials themselves will be identical, then, if their third derivatives are also identical at $x=x_{1}$ , or if

$a_{0}=a_{1} \nonumber$

From (5.13), we have

$\frac{1}{3 h_{0}}\left(b_{1}-b_{0}\right)=\frac{1}{3 h_{1}}\left(b_{2}-b_{1}\right) \nonumber$

or after simplification

$h_{1} b_{0}-\left(h_{0}+h_{1}\right) b_{1}+h_{0} b_{2}=0 \nonumber$

which provides us our missing first equation. A similar argument at $x=x_{n}-1$ also provides us with our last equation,

$h_{n-1} b_{n-2}-\left(h_{n-2}+h_{n-1}\right) b_{n-1}+h_{n-2} b_{n}=0 . \nonumber$

The MATLAB subroutines spline.m and ppval.m can be used for cubic spline interpolation (see also interp1.m). I will illustrate these routines in class and post sample code on the course web site.

Multidimensional interpolation

Suppose we are interpolating the value of a function of two variables,

$z=f(x, y) \nonumber$

The known values are given by

$z_{i j}=f\left(x_{i}, y_{j}\right) \nonumber$

with $i=0,1, \ldots, n$ and $j=0,1, \ldots, n$ . Note that the $(x, y)$ points at which $f(x, y)$ are known lie on a grid in the $x-y$ plane.

Let $z=g(x, y)$ be the interpolating function, satisfying $z_{i j}=g\left(x_{i}, y_{j}\right) .$ A twodimensional interpolation to find the value of $g$ at the point $(x, y)$ may be done by first performing $n+1$ one-dimensional interpolations in $y$ to find the value of $g$ at the $n+1$ points $\left(x_{0}, y\right),\left(x_{1}, y\right), \ldots,\left(x_{n}, y\right)$ , followed by a single one-dimensional interpolation in $x$ to find the value of $g$ at $(x, y)$ .

In other words, two-dimensional interpolation on a grid of dimension $(n+1) \times$ $(n+1)$ is done by first performing $n+1$ one-dimensional interpolations to the value $y$ followed by a single one-dimensional interpolation to the value $x$ . Twodimensional interpolation can be generalized to higher dimensions. The MATLAB functions to perform two-and three-dimensional interpolation are interp2.m and interp3.m. 5.4. MULTIDIMENSIONAL INTERPOLATION

Chapter 6 Integration

We want to construct numerical algorithms that can perform definite integrals of the form

$I=\int_{a}^{b} f(x) d x \nonumber$

Calculating these definite integrals numerically is called numerical integration, numerical quadrature, or more simply quadrature.

Elementary formulas

We first consider integration from 0 to $h$ , with $h$ small, to serve as the building blocks for integration over larger domains. We here define $I_{h}$ as the following integral:

$I_{h}=\int_{0}^{h} f(x) d x . \nonumber$

To perform this integral, we consider a Taylor series expansion of $f(x)$ about the value $x=h / 2$ :

$\begin{aligned} f(x)=f(h / 2)+(x-h / 2) f^{\prime} &(h / 2)+\frac{(x-h / 2)^{2}}{2} f^{\prime \prime}(h / 2) \\ &+\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 \end{aligned} \nonumber$

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$

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

Composite rules

We now use our elementary formulas obtained for (6.2) to perform the integral given by $(6.1)$

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

Local versus global error

Consider the elementary formula (6.4) for the trapezoidal rule, written in the form

$\int_{0}^{h} f(x) d x=\frac{h}{2}(f(0)+f(h))-\frac{h^{3}}{12} f^{\prime \prime}(\xi) \nonumber$

where $\xi$ is some value satisfying $0 \leq \xi \leq h$ , and we have used Taylor’s theorem with the mean-value form of the remainder. We can also represent the remainder as

$-\frac{h^{3}}{12} f^{\prime \prime}(\xi)=\mathrm{O}\left(h^{3}\right) \nonumber$

where $\mathrm{O}\left(h^{3}\right)$ signifies that when $h$ is small, halving of the grid spacing $h$ decreases the error in the elementary trapezoidal rule by a factor of eight. We call the terms represented by $\mathrm{O}\left(h^{3}\right)$ the Local Error.

More important is the Global Error which is obtained from the composite formula (6.7) for the trapezoidal rule. Putting in the remainder terms, we have

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

where $\xi_{i}$ are values satisfying $x_{i} \leq \xi_{i} \leq x_{i+1}$ . The remainder can be rewritten as

$-\frac{h^{3}}{12} \sum_{i=0}^{n-1} f^{\prime \prime}\left(\xi_{i}\right)=-\frac{n h^{3}}{12}\left\langle f^{\prime \prime}\left(\xi_{i}\right)\right\rangle \nonumber$

where $\left\langle f^{\prime \prime}\left(\xi_{i}\right)\right\rangle$ is the average value of all the $f^{\prime \prime}\left(\xi_{i}\right)^{\prime}$ s. Now,

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

so that the error term becomes

$\begin{aligned} -\frac{n h^{3}}{12}\left\langle f^{\prime \prime}\left(\xi_{i}\right)\right\rangle &=-\frac{(b-a) h^{2}}{12}\left\langle f^{\prime \prime}\left(\xi_{i}\right)\right\rangle \\ &=\mathrm{O}\left(h^{2}\right) \end{aligned} \nonumber$

Therefore, the global error is $\mathrm{O}\left(h^{2}\right)$ . That is, a halving of the grid spacing only decreases the global error by a factor of four.

Similarly, Simpson’s rule has a local error of $\mathrm{O}\left(h^{5}\right)$ and a global error of $\mathrm{O}\left(h^{4}\right)$ .

Figure 6.1: Adaptive Simpson quadrature: Level $1 .$

Adaptive integration

The useful MATLAB function quad.m performs numerical integration using adaptive Simpson quadrature. The idea is to let the computation itself decide on the grid size required to achieve a certain level of accuracy. Moreover, the grid size need not be the same over the entire region of integration.

We begin the adaptive integration at what is called Level $1 .$ The uniformly spaced points at which the function $f(x)$ is to be evaluated are shown in Fig. 6.1. The distance between the points $a$ and $b$ is taken to be $2 h$ , so that

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

Integration using Simpson’s rule (6.5) with grid size $h$ yields

$I=\frac{h}{3}(f(a)+4 f(c)+f(b))-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(\xi) \nonumber$

where $\xi$ is some value satisfying $a \leq \xi \leq b$ .

Integration using Simpson’s rule twice with grid size $h / 2$ yields

$I=\frac{h}{6}(f(a)+4 f(d)+2 f(c)+4 f(e)+f(b))-\frac{(h / 2)^{5}}{90} f^{\prime \prime \prime \prime}\left(\xi_{l}\right)-\frac{(h / 2)^{5}}{90} f^{\prime \prime \prime \prime}\left(\xi_{r}\right) \nonumber$

with $\xi_{l}$ and $\xi_{r}$ some values satisfying $a \leq \xi_{l} \leq c$ and $c \leq \xi_{r} \leq b$ .

We now define

$\begin{aligned} &S_{1}=\frac{h}{3}(f(a)+4 f(c)+f(b)) \\ &S_{2}=\frac{h}{6}(f(a)+4 f(d)+2 f(c)+4 f(e)+f(b)), \\ &E_{1}=-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(\xi), \\ &E_{2}=-\frac{h^{5}}{2^{5} \cdot 90}\left(f^{\prime \prime \prime \prime}\left(\xi_{l}\right)+f^{\prime \prime \prime \prime}\left(\xi_{r}\right)\right) . \end{aligned} \nonumber$

Now we ask whether $S_{2}$ is accurate enough, or must we further refine the calculation and go to Level 2? To answer this question, we make the simplifying approximation that all of the fourth-order derivatives of $f(x)$ in the error terms are equal; that is,

$f^{\prime \prime \prime \prime}(\xi)=f^{\prime \prime \prime \prime}\left(\xi_{l}\right)=f^{\prime \prime \prime \prime}\left(\xi_{r}\right)=C . \nonumber$

CHAPTER 6. INTEGRATION

Then

$\begin{aligned} &E_{1}=-\frac{h^{5}}{90} C \\ &E_{2}=-\frac{h^{5}}{2^{4} \cdot 90} C=\frac{1}{16} E_{1} \end{aligned} \nonumber$

Then since

$S_{1}+E_{1}=S_{2}+E_{2} \nonumber$

and

$E_{1}=16 E_{2} \nonumber$

we have for our estimate for the error term $E_{2}$ ,

$E_{2}=\frac{1}{15}\left(S_{2}-S_{1}\right) \nonumber$

Therefore, given some specific value of the tolerance tol, if

$\left|\frac{1}{15}\left(S_{2}-S_{1}\right)\right|<\text { tol } \nonumber$

then we can accept $S_{2}$ as $I$ . If the tolerance is not achieved for $I$ , then we proceed to Level $2 .$

The computation at Level 2 further divides the integration interval from $a$ to $b$ into the two integration intervals $a$ to $c$ and $c$ to $b$ , and proceeds with the above procedure independently on both halves. Integration can be stopped on either half provided the tolerance is less than tol/2 (since the sum of both errors must be less than tol). Otherwise, either half can proceed to Level 3, and so on.

As a side note, the two values of $I$ given above (for integration with step size $h$ and $h / 2$ ) can be combined to give a more accurate value for I given by

$I=\frac{16 S_{2}-S_{1}}{15}+\mathrm{O}\left(h^{7}\right) \nonumber$

where the error terms of $\mathrm{O}\left(h^{5}\right)$ approximately cancel. This free lunch, so to speak, is called Richardson’s extrapolation.

Chapter 7 Ordinary differential equations

We now discuss the numerical solution of ordinary differential equations. These include the initial value problem, the boundary value problem, and the eigenvalue problem. Before proceeding to the development of numerical methods, we review the analytical solution of some classic equations.

Examples of analytical solutions

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.

Numerical methods: initial value problem

We begin with the simple Euler method, then discuss the more sophisticated RungeKutta methods, and conclude with the Runge-Kutta-Fehlberg method, as implemented in the MATLAB function ode45.m. Our differential equations are for $x=$ $x(t)$ , where the time $t$ is the independent variable, and we will make use of the notation $\dot{x}=d x / d t$ . This notation is still widely used by physicists and descends directly from the notation originally used by Newton.

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.

order	2	3	4	5	6	7	8
minimum # stages	2	3	4	6	7	9	11

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$

Numerical methods: boundary value problem

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

Numerical methods: eigenvalue problem

For illustrative purposes, we develop our numerical methods for what is perhaps the simplest eigenvalue ode. With $y=y(x)$ and $0 \leq x \leq 1$ , this simple ode is given by

$y^{\prime \prime}+\lambda^{2} y=0 \nonumber$

To solve (7.9) numerically, we will develop both a finite difference method and a shooting method. Furthermore, we will show how to solve $(7.9)$ with homogeneous boundary conditions on either the function $y$ or its derivative $y^{\prime}$ .

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$

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

Search

Text Color

Text Size

Margin Size

Font Type

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

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

Example of fractals using Newton’s Method

Newton’s Method

Secant Method

Vandermonde polynomial

Lagrange polynomial

Newton polynomial

Midpoint rule

Trapezoidal rule

Simpson’s rule

Trapezoidal rule

Simpson’s rule

Initial value problem

Boundary value problems

Eigenvalue problem

Euler method

Modified Euler method

Second-order Runge-Kutta methods

Higher-order Runge-Kutta methods

Adaptive Runge-Kutta Methods

System of differential equations

Finite difference method

Shooting method

Finite difference method

Shooting method

Chapter 1

Arithmetic

Definitions

Numbers with a decimal or binary point

Examples of binary numbers

numbers

unsigned integers as hex numbers

single precision format:

Special numbers

Examples of computer numbers

Inexact numbers

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

Machine epsilon

double precision format

1.121.12 Roundoff error example

Chapter 2

Finding

Bisection Method

2.22.2 Newton’s Method

Secant Method

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

Example of fractals using Newton’s Method

Order of convergence

Newton’s Method

Secant Method

Chapter 3

Systems of equations

Gaussian Elimination

U\) decomposition

Partial pivoting

3.43.4 Operation counts

3.53.5 System of nonlinear equations

Chapter 4

Least-squares approximation

Fitting a straight line

Fitting to a linear combination of functions

Chapter 5

Interpolation

Polynomial interpolation

Vandermonde polynomial

Lagrange polynomial

Newton polynomial

Piecewise linear interpolation

Cubic spline interpolation

Multidimensional interpolation

Chapter 6

Integration

Elementary formulas

Midpoint rule

Trapezoidal rule

Simpson’s rule

Composite rules

Trapezoidal rule

Simpson’s rule

Local versus global error

Adaptive integration

CHAPTER 6. INTEGRATION

Chapter 7

Ordinary differential equations

Examples of analytical solutions

Initial value problem

Boundary value problems

Eigenvalue problem

Numerical methods: initial value problem

Euler method

Modified Euler method

Second-order Runge-Kutta methods

Higher-order Runge-Kutta methods

Adaptive Runge-Kutta Methods

System of differential equations

Numerical methods: boundary value problem

Finite difference method

Shooting method

Numerical methods: eigenvalue problem

Finite difference method

Shooting method

Support Center

How can we help?

$1.12$ Roundoff error example

$2.2$ Newton’s Method

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

$3.4$ Operation counts

$3.5$ System of nonlinear equations