Numberical Analysis --- Numerical Differentiation and Integration

Numerical Differentiation

Given $x_0$, approximate $f'(x_0)$.
$$f'(x_0)=\lim\limits_{h\to 0}\frac{f(x_0+h)-f(x_0)}{h}$$
$$\Rightarrow f'(x_0)\approx \frac{f(x_0+h)-f(x_0)}{h}(\text{forward})$$
$$\Rightarrow f'(x_0)\approx \frac{f(x_0)-f(x_0-h)}{h}(\text{backward})$$

Approximate $f(x)$ by its Lagrange polynomial with interpolating points $x_0$ and $x_0+h$:
$$\begin{aligned} f(x)&=\frac{f(x_0)(x-x_0-h)}{x_0-x_0-h}+\frac{f(x_0+h)(x-x_0)}{x_0+h-x_0}\\ &+\frac{(x-x_0)(x-x_0-h)}{2}f''(\xi_x)\\ \end{aligned}$$

$$\begin{aligned} f'(x)&=\frac{f(x_0+h)-f(x_0)}{h}+\frac{2(x-x_0)-h}{2}f''(\xi_x)\\ &+\frac{(x-x_0)(x-x_0-h)}{2}\cdot \frac{\,d}{\,dx}[f''(\xi_x)]\\ \end{aligned}$$
$$\Rightarrow f'(x_0)=\frac{f(x_0+h)-f(x_0)}{h}-\frac{h}{2}f''(\xi)$$

Approximate $f(x)$ by its Lagrange polynomial with interpolating points $\{x_0,x_1,...,x_n\}$.
$$f(x)=\sum\limits_{k=0}^nf(x_k)L_k(x)+\frac{(x-x_0)...(x-x_n)}{(n+1)!}f^{(n+1)}(\xi_x)$$
$$f'(x_j)=\sum\limits_{k=0}^nf(x_k)L_k'(x_j)+\frac{f^{(n+1)}(\xi_j)}{(n+1)!}\prod\limits_{0\leq k\leq n,k\neq j}(x_j-x_k)$$

Note:

In general, more evaluation points produce greater accuracy.

On the other hand, the number of functional evaluations grows and the roudoff error increases. Hence the numerical differentiation is unstable!

Example: Given three points $x_0,x_0+h$ and $x_0+2h$, please derive the three-point formulae for each of the points. Then give the best three-point formulae for $f'(x)$.
$$f'(x_0)=\frac{1}{h}\left[-\frac{3}{2}f(x_0)+2f(x_0+h)-\frac{1}{2}f(x_0+2h)\right]+\frac{h^2}{3}f^{(3)}(\xi_0)$$
$$f'(x_0+h)=\frac{1}{h}\left[-\frac{1}{2}f(x_0)+\frac{1}{2}f(x_0+2h)\right]-\frac{h^2}{6}f^{(3)}(\xi_1)$$
$$f'(x_0)=\frac{1}{2h}\left[f(x_0+h)-f(x_0-h)\right]-\frac{h^2}{6}f^{(3)}(\xi_1)$$

Example: Find a way to approximate $f''(x_0)$.

Consider Taylor expansion of $f(x_0+h)$ and $f(x_0-h)$ at $x_0$:
$$f(x_0+h)=f(x_0)+f'(x_0)h+\frac{1}{2}f''(x_0)h^2+\frac{1}{6}f'''(x_0)h^3+\frac{1}{24}f^{(4)}(\xi_1)h^4$$
$$f(x_0-h)=f(x_0)-f'(x_0)h+\frac{1}{2}f''(x_0)h^2-\frac{1}{6}f'''(x_0)h^3+\frac{1}{24}f^{(4)}(\xi_{-1})h^4$$

$$f''(x_0)=\frac{1}{h^2}[f(x_0-h)-2f(x_0)+f(x_0+h)]-\frac{h^2}{12}f^{(4)}(\xi)$$

Elements of Numerical Integration

Integrate the Lagrange interpolating polynomial of $f(x)$ instead.

Select a set of distinct nodes $a\leq x_0<x_1<...<x_n\leq b$ from $[a,b]$. The Lagrange polynomial is $P_n(x)=\sum\limits_{k=0}^nf(x_k)L_k(x)$
$$\int_a^bf(x)\,dx\approx \sum\limits_{k=0}^nf(x_k)\int_a^bL_k(x)\,dx$$
$$A_k=\int_a^bL_k(x)\,dx=\int_a^b\prod\limits_{j\neq k}\frac{(x-x_j)}{(x_k-x_j)}\,dx$$
$$\begin{aligned} \text{Error} R[f]&=\int_a^bf(x)\,dx-\sum\limits_{k=0}^nA_kf(x_k)\\ &=\int_a^b[f(x)-P_n(x)]\,dx=\int_a^bR_n(x)\,dx\\ &=\int_a^b\frac{f^{(n+1)}(\xi_x)}{(n+1)!}\prod\limits_{k=0}^n(x-x_k)\,dx\\ \end{aligned}$$

Definition: The degree of accuracy, or precision, of a quadrature formulae is the largest positive integer $n$ such that the formulae is exact for $x^k$ for each $k=0,1,...,n$.

Example: Consider the linear interpolation on $[a,b]$, we have
$$P_1(x)=\frac{x-b}{a-b}f(a)+\frac{x-a}{b-a}f(b)$$
$$\Rightarrow A_1=A_2=\frac{b-a}{2}\Rightarrow \int_a^bf(x)\,dx\approx\frac{b-a}{2}[f(a)+f(b)]$$

Solution: Consider $x^k$ for each $k=0,1,...$
$$x^0=1:\int_a^b1\,dx=b-a=\frac{b-a}{2}[1+1]$$
$$x:\int_a^bx\,dx=\frac{b^2-a^2}{2}=\frac{b-a}{2}[a+b]$$
$$x^2:\int_a^bx^2\,dx=\frac{b^3-a^3}{3}\neq \frac{b-a}{2}[a^2+b^2]$$

For equally spaced nodes: $x_i=a+ih,h=\frac{b-a}{n},i=0,1,...,n$
$$\begin{aligned} A_i&=\int_{x_0}^{x_n}\prod\limits_{j\neq i}\frac{(x-x_j)}{(x_i-x_j)}\,dx\\ &=\int_0^n\prod\limits_{i\neq j}\frac{(t-j)h}{(i-j)h}\times h\,dt\\ &=\frac{(b-a)(-1)^{n-1}}{ni!(n-i)!}\int_0^n\prod\limits_{i\neq j}(t-j)\,dt\\ \end{aligned}$$

Note: Cotes coefficients does not depend on either $f(x)$ or $[a,b]$, and can be determined by $n$ and $i$ only. Hence we can find these coefficients from a table. The formulae are called Newton-Cotes formulae.

   $n=1$: $C_0^{(1)}=\frac{1}{2},C_1^{(1)}=\frac{1}{2}$ Precision=1
$$\text{Trapezoidal Rule: }\int_a^bf(x)\,dx\approx\frac{b-a}{2}[f(a)+f(b)]$$
$$\begin{aligned} R[f]&=\int_a^b\frac{f''(\xi_x)}{2!}(x-a)(x-b)\,dx\\ &=-\frac{1}{12}h^3f''(\xi)\\ \end{aligned}\quad \xi\in[a,b],h=b-a$$

   $n=2$: $C_0^{(2)}=\frac{1}{6},C_1^{(2)}=\frac{2}{3},C_2^{(2)}=\frac{1}{6}$ Precision=3.
$$\text{Simpson's Rule: }\int_a^bf(x)\,dx\approx \frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]$$
$$R[f]=-\frac{1}{90}h^5f^{(4)}(\xi),\quad \xi\in(a,b),h=\frac{b-a}{2}$$

   $n=3$: Simpson’s 3/8-Rule. Precision=3.
$$R[f]=-\frac{3}{80}h^5f^{(5)}(\xi)$$

   $n=4$: Cotes Rule. Precision=5.
$$R[f]=-\frac{8}{945}h^7f^{(6)}(\xi)$$

Theorem: For the $(n+1)$-point closed Newton-Cotes formulae, there exists $\xi\in(a,b)$ for which
$$\int_a^bf(x)\,dx=\sum\limits_{k=0}^nA_kf(x_k)+\frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!}\int_0^nt^2(t-1)...(t-n)\,dt$$
if $n$ is even and $f\in C^{n+2}[a,b]$, and
$$\int_a^bf(x)\,dx=\sum\limits_{k=0}^nA_kf(x_k)+\frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!}\int_0^nt(t-1)...(t-n)\,dt$$
if $n$ is odd and $f\in C^{n+1}[a,b]$

Composite Numerical Integration

Due to the osillatory nature of high-degree polynomials, piecewise interpolation is applied to approximate $f(x)\Rightarrow$ a piecewise approach that uses the low-order Newton-Cotes formulae.

Composite Trapezoidal Rule:$h=\frac{b-a}{n},x_k=a+kh\quad(k=0,...,n)$

Apply Trapezoidal Rule on each $[x_{k-1},x_k]$:
$$\int_{x_{k-1}}^{x_k}f(x)\,dx\approx \frac{x_k-x_{k-1}}{2}[f(x_{k-1})+f(x_k)],k=1,...,n$$
$$\int_a^bf(x)\,dx\approx\sum\limits_{k=1}^n\frac{h}{2}[f(x_{k-1})+f(x_k)]=\frac{h}{2}[f(a)+2\sum\limits_{k=1}^{n-1}f(x_k)+f(b)]=T_n$$
$$\begin{aligned} R[f]&=\sum\limits_{k=1}^n\left[-\frac{h^3}{12}f''(\xi_x)\right]=-\frac{h^2}{12}(b-a)\frac{\sum\limits_{k=1}^nf''(\xi_x)}{n}\\ &=-\frac{h^2}{12}(b-a)f''(\xi),\xi\in(a,b)\\ \end{aligned}$$

Composite Simpson;s Rule:$h=\frac{b-a}{n},x_k=a+kh\quad(k=0,...,n)$
$$\int_{x_k}^{x_{k+1}}f(x)\,dx\approx \frac{h}{6}[f(x_k)+4f(x_{k+\frac{1}{2}})+f(x_{k+1})]$$
$$\int_a^bf(x)\,dx=\frac{h}{6}[f(a)+4\sum\limits_{k=0}^{n-1}f(x_{k+\frac{1}{2}})+2\sum\limits_{k=0}^{n-2}f(x_{k+1})+f(b)]=S_n$$
$$R[f]=-\frac{b-a}{180}\left(\frac{h}{2}\right)^4f^{(4)}(\xi)$$

Note: To simplify the notation, we may let $n'=2n$. Then $h'= \frac{b-a}{n'}=\frac{h}{2},x_k=a+kh'$ and
$$S_n=\frac{h'}{3}[f(a)+4\sum\limits_{\text{odd } k}f(x_k)+2\sum\limits_{\text{even } k}f(x_k)+f(b)]$$

Composite integration techniques are all stable.

Example: Consider the Simpson’s Rule with $n$ subintervals on $[a,b]$. Assume that $f(x_i)$ is approximated by $f^*(x_i)$ such that $f(x_i)=f^*(x_i)+\epsilon_i$ for each $i=0,1,...,n$. Then the accumulated error $e(h)$ is
$$e(h)=\left|\frac{h}{3}[\epsilon_0+4\sum\limits_{\text{odd }k}\epsilon_k+2\sum\limits_{\text{even }k}\epsilon_k+\epsilon_n]\right|$$
If $|\epsilon_i|<\epsilon$ for all $i=0,1,...,n$, then
$$e(h)<\frac{h}{3}[\epsilon+4(n/2)\epsilon+2(n/2-1)\epsilon+\epsilon]=nh\epsilon=(b-a)\epsilon$$

When we refine the partition to ensure accuracy, the increased computation will NOT increase the roundoff error.

Example: Use Trapezoidal rule and Simpson’s rule with $n=8$ to approximate $\pi=\int_0^1\frac{4}{1+x^2}\,dx$.

Solution: Set $x_k=\frac{k}{8}$
$$T_8=\frac{1}{16}[f(0)+2\sum\limits_{k=1}^7f(x_k)+f(1)]=3.138988494$$
$$S_4=\frac{1}{24}[f(0)+4\sum\limits_{\text{odd}}f(x_k)+2\sum\limits_{\text{even}}f(x_k)+f(1)]=3.141592502$$

When programming we usually keep dividing the subintervals into 2 equally spaced smaller subintervals. That is, take $n=2^k$ for $k=0,1,...$. When $k=9,T_{512}=3.14159202$
$$\frac{4}{3}T_8-\frac{1}{3}T_4=3.141592502=S_4$$

Romberg Integration

Since the error of Trapezoidal rule is $R_n[f]=-\frac{h^2}{12}(b-a)f''(\xi)$, when we reduce the length of each subinterval into a half,
$$R_{2n}[f]=-\left(\frac{h}{2}\right)^2\frac{1}{12}f''(\xi')\approx \frac{1}{4}R_n[f]$$
$$\Rightarrow \frac{I-T_{2n}}{I-T_n}\approx \frac{1}{4}$$
Solve for $I$:
$$I\approx \frac{4T_{2n}-T_n}{4-1}=\frac{4}{3}T_{2n}-\frac{1}{3}T_n=S_n$$
In general:
$$\frac{4T_{2n}-T_n}{4-1}=S_n\quad \frac{4^2S_{2n}-S_n}{4^2-1}=C_n\quad \frac{4^3C_{2n}-C_n}{4^3-1}=R_n$$

Richardson’s Extrapolation

Generate high-accuracy results while using lower-order formulae. Suppose that for some $h\neq 0$, we have a formula $T_0(x)$ that approximates an unknown $I$, and that the truncation error has the form:
$$T_0(h)-I=\alpha_1h+\alpha_2h^2+...$$
Replace $h$ by half its value, we have
$$T_0\left(\frac{h}{2}\right)-I=\alpha_1(h/2)+\alpha_2(h/2)^2+...$$

How to improve the accuracy from $O(h)$ to $O(h^2)$?
$$\frac{2T_0(\frac{h}{2})-T_0(h)}{2-1}-I=-\frac{1}{2}\alpha_2h^2-\frac{3}{4}\alpha_3h^3-...$$

$$T_1(h)=\frac{2T_0(\frac{h}{2})-T_0(h)}{2-1}=I+\beta_1h^2+\beta_2h^3+...$$
$$T_2(h)=\frac{2^2T_1(\frac{h}{2})-T_1(h)}{2-1}=I+\gamma_1h^3+\gamma_4h^4+...$$
$$\Rightarrow T_m(h)=\frac{2^mT_{m-1}(\frac{h}{2})-T_{m-1}(h)}{2^m-1}=I+\delta_1h^{m+1}+\delta_2h^{m+2}+...$$

Adaptive Quadrature Methods

Predict the amount of functional variation and adapt the step size to the varying requirements.

Example: Consider the adaptive method based on Composite Simpson’s rule.
$$\int_a^bf(x)\,dx=S(a,b)-\frac{h^5}{90}f^{(4)}(\xi)$$
where
$$S(a,b)=\frac{h}{3}[f(a)+4f(a+h)+f(b)],\quad h=\frac{b-a}{2}$$

Let
$$S\left(a,\frac{a+b}{2}\right)=\frac{h}{6}[f(a)+4f\left(a+\frac{h}{2}\right)+f(a+h)]$$
$$S\left(\frac{a+b}{2},b\right)=\frac{h}{6}[f(a+h)+4f\left(a+\frac{3h}{2}\right)+f(b)]$$
Then
$$\int_a^bf(x)\,dx=S\left(a,\frac{a+b}{2}\right)+S\left(\frac{a+b}{2},b\right)-\frac{1}{16}\times \frac{h^5}{90}f^{(4)}(\eta)$$
$$\left|\int_a^bf(x)\,dx-S\left(a,\frac{a+b}{2}\right)-S\left(\frac{a+b}{2},b\right)\right|\approx \frac{1}{15}\left|S(a,b)-S\left(a,\frac{a+b}{2}\right)-S\left(\frac{a+b}{2},b\right)\right|<\epsilon$$

Gaussian Quadrature

Construct a formula $\int_a^bw(x)f(x)\,dx\approx\sum\limits_{k=0}^nA_kf(x_k)$ of precision degree $2n+1$ with $n+1$ points.

Determine $2n+2$ unknowns $x_0,...,x_n$ and $A_0,..,A_n$ such that the formula is accurate for $f(x)=1,x,x^2,...,x^{2n+1}$. The points $x_0,...,x_n$ are called Gaussian points, and the method is called Gaussian quadrature.

Example: Approximate $\int_0^1\sqrt{x}f(x)\,dx$ using Gaussian quadrature with $n=1$.

Solution: Assume $\int_0^1\sqrt{x}f(x)\,dx=A_0f(x_0)+A_1f(x_1)$
The formula must be accurate for $f(x)=1,x,x^2,x^3$
$$\begin{cases} \frac{2}{3}=A_0+A_1\\ \frac{2}{5}=A_0x_0+A_1x_1\\ \frac{2}{7}=A_0x_0^2+A_1x_1^2\\ \frac{2}{9}=A_0x_0^3+A_1x_1^3\\ \end{cases}\Rightarrow \begin{cases} x_0\approx 0.8212\\ x_1\approx 0.2899\\ A_0\approx 0.3891\\ A_1\approx 0.2776\\ \end{cases}$$

Theorem: $x_0,...,x_n$ are Gaussian points iff $W(x)=\prod\limits_{k=0}^n(x-x_k)$ is orthogonal to all the polynomials of degree no greater than $n$.

Proof:

1. If $x_0,...,x_n$ are Gaussian points, then the degree of precision of the formula $\int_a^bw(x)f(x)\,dx\approx\sum\limits_{k=0}^nA_kf(x_k)$ is at least $2n+1$. Then for any polynomial $P_m(x)$ with $m\leq n$, the degree of $P_m(x)W(x)$ must be no greater than $2n+1$. Hence the formula is accurate for $P_m(x)W(x)$. That is
   $$\int_a^bw(x)P_m(x)W(x)\,dx=\sum\limits_{k=0}^nA_kP_m(x_k)W(x_k)=0$$

2. To prove that $x_0,...,x_n$ are Gaussian points, we just have to prove that the formula is accurate for any polynomial $P_m(x)$ with $m\leq 2n+1$. Let $P_m(x)=W(x)q(x)+r(x)$. Then
   $$\begin{aligned} \int_a^bw(x)P_m(x)\,dx&=\int_a^bw(x)W(x)q(x)\,dx+\int_a^bw(x)r(x)\,dx\\ &=\sum\limits_{k=0}^nA_kr(x_k)=\sum\limits_{k=0}^nA_kP_m(x_k)\\ \end{aligned}$$

The set of orthogonal polynomials $\{\varphi_0,\varphi_1,...,\varphi_n,...\}$ is linearly independent and $\varphi_{n+1}$ is orthogonal to any polynomial $P_m(x)$ with $m\leq n$.$\Rightarrow$ If we take $\varphi_{n+1}$ to be $W(x)$, then the roots of $\varphi_{n+1}$ are the Gaussian points.

Example: Approximate $\int_0^1\sqrt{x}f(x)\,dx$ using Gaussian quadrature with $n=1$.

  Step 1:

Construct the orthogonal polynomial $\varphi_2$. Let $\varphi_0(x)=1,\varphi_1(x)=x+a,\varphi_2(x)=x^2+bx+c$
$$(\varphi_0,\varphi_1)=0\Rightarrow \int_0^1\sqrt{x}(x+a)\,dx=0\Rightarrow a=-\frac{3}{5}$$
$$(\varphi_0,\varphi_2)=0\Rightarrow \int_0^1\sqrt{x}(x^2+bx+c)\,dx=0$$
$$(\varphi_1,\varphi_2)=0\Rightarrow \int_0^1\sqrt{x}(x-\frac{3}{5})(x^2+bx+c)\,dx=0$$
$$\Rightarrow \begin{cases} b=-\frac{10}{9}\\ c=\frac{5}{21}\\ \end{cases}\Rightarrow \varphi_2(x)=x^2-\frac{10}{9}x+\frac{5}{21}$$

   Step 2:

Find the two roots of $\varphi_2$ which are the Gaussian points $x_0$ and $x_1$:
$$x_{0;1}=\frac{10/9\pm \sqrt{(10/9)^2-20/21}}{2}$$

   Step 3:

Since the formula must be accurate for $f(x)=1,x$, we can easily solve a linear system of equations for $A_0$ and $A_1$. The results are the same as we have obtained:
$$x_0\approx0.8212,x_1\approx0.2899,A_0\approx0.3891,A_1\approx0.2776$$

Use this formula to approximate $\int_0^1\sqrt{x}e^x\,dx\approx1.2555$

  Some special families of orthogonal functions

1.Legendre Polynomials: Defined on $[-1,1]$ with $w(x)\equiv 1$.
$$P_k(x)=\frac{1}{2^kk!}\frac{\,d^k}{\,dx^k}(x^2-1)^k\quad (P_k,P_l)=\begin{cases} 0\quad k\neq l\\ \frac{2}{2k+1}\quad k=l\\ \end{cases}$$

From $P_0=1$ and $P_1=x$ we have $(k+1)P_{k+1}=(2k+1)xP_k-kP_{k-1}$

The formula that uses the roots of $P_{n+1}$ is called the Gauss-Legendre quadrature formula.

2.Chebyshev Polynomials:

Defined on $[-1,1]$ with $w(x)=\frac{1}{\sqrt{1-x^2}}$

The roots of $T_{n+1}$ are $x_k=\cos \left(\frac{2k+1}{2n+2}\pi\right)$ for $k=0,...,n$.

The formula $$\int_{-1}^1\frac{1}{\sqrt{1-x^2}}f(x)\,dx=\sum\limits_{k=0}^nA_kf(x_k)$$
is called the Gauss-Chebyshev quadrature formula.