Numerical Differentiation

Forward and Backward Differences

Inspired by the definition of derivative: \[ \begin{aligned} f'(x_0) = \lim_{h\rightarrow 0} \frac{f(x_0+h)-f(x_0)}{h}, \end{aligned} \] choose a small \(h\) and approximate \[ \begin{aligned} f'(x_0) \approx \frac{f(x_0+h)-f(x_0)}{h} \end{aligned} \] The error term for the linear Lagrange polynomial gives: \[ \begin{aligned} f'(x_0) = \frac{f(x_0+h)-f(x_0)}{h} -\frac{h}{2}f''(\xi) \end{aligned} \] Also known as the forward-difference formula if \(h>0\) and the backward-difference formula if \(h<0\)

General Derivative Approximations

Differentiation of Lagrange Polynomials

Differentiate \[ \begin{aligned} f(x)=\sum_{k=0}^n f(x_k)L_k(x) + \frac{(x-x_0)\cdots (x-x_n)}{(n+1)!} f^{(n+1)}(\xi(x)) \end{aligned} \] to get \[ \begin{aligned} f'(x_j) = \sum_{k=0}^n f(x_k)L'_k(x_j)+ \frac{f^{(n+1)}(\xi(x_j))}{(n+1)!}\prod_{k\ne j} (x_j-x_k) \end{aligned} \] This is the \((n+1)\)-point formula for approximating \(f'(x_j)\).

Commonly Used Formulas

Using equally spaced points with \(h=x_{j+1}-x_j\), we have the three-point formulas \[ \begin{aligned} f'(x_0) &= \frac{1}{2h}[-3f(x_0)+4f(x_0+h)-f(x_0+2h)]+\frac{h^2}{3}f^{(3)}(\xi_0) \\ f'(x_0) &= \frac{1}{2h}[-f(x_0-h)+f(x_0+h)] - \frac{h^2}{6}f^{(3)}(\xi_1) \\ f'(x_0) &= \frac{1}{2h}[f(x_0-2h)-4f(x_0-h)+3f(x_0)] + \frac{h^2}{3}f^{(3)}(\xi_2) \\ 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) \end{aligned} \] and the five-point formula \[ \begin{aligned} f'(x_0) &= \frac{1}{12h}[f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)] \\ &+ \frac{h^4}{30} f^{(5)}(\xi) \end{aligned} \]

Optimal \(h\)

• Consider the three-point central difference formula: \[ \begin{aligned} f'(x_0) = \frac{1}{2h}[f(x_0+h)-f(x_0-h)] - \frac{h^2}{6}f^{(3)}(\xi_1) \end{aligned} \]
• Suppose that round-off errors \(\varepsilon\) are introduced when computing \(f\). Then the approximation error is \[ \begin{aligned} \left|f'(x_0)-\frac{\tilde{f}(x_0+h)-\tilde{f}(x_0-h)}{2h}\right| \le \frac{\varepsilon}{h}+\frac{h^2}{6}M = e(h) \end{aligned} \] where \(\tilde{f}\) is the computed function and \(|f^{(3)}(x)|\le M\)
• Sum of truncation error \(h^2M/6\) and round-off error \(\varepsilon/h\)
• Minimize \(e(h)\) to find the optimal \(h=\sqrt[3]{3\varepsilon/M}\)

Richardson’s Extrapolation

• Suppose \(N(h)\) approximates an unknown \(M\) with error \[ \begin{aligned} M-N(h) = K_1h + K_2h^2+K_3h^3+\cdots \end{aligned} \] then an \(O(h^j)\) approximation is given for \(j=2,3,\ldots\) by \[ \begin{aligned} N_j(h) = N_{j-1}\left(\frac{h}{2}\right) + \frac{N_{j-1}(h/2)-N_{j-1}(h)}{2^{j-1}-1} \end{aligned} \]
• The results can be written in a table:

\[ \begin{array}{cccc} \hline O(h) & O(h^2) & O(h^3) & O(h^4) \\ \hline \textbf{1:} N_1(h)\equiv N(h) \\ \textbf{2:} N_1(\frac{h}{2})\equiv N(\frac{h}{2}) & \textbf{3:} N_2(h) \\ \textbf{4:} N_1(\frac{h}{4})\equiv N(\frac{h}{4}) & \textbf{5:} N_2(\frac{h}{2}) & \textbf{6:} N_3(h) \\ \textbf{7:} N_1(\frac{h}{8})\equiv N(\frac{h}{8}) & \textbf{8:} N_2(\frac{h}{4}) & \textbf{9:} N_3(\frac{h}{2}) & \textbf{10:} N_4(h) \end{array} \]

Richardson’s Extrapolation

• If some error terms are zero, different and more efficient formulas can be derived
• Example: If \[ \begin{aligned} M-N(h) = K_2h^2+K_4h^4+\cdots \end{aligned} \] then an \(O(h^{2j})\) approximation is given for \(j=2,3,\ldots\) by \[ \begin{aligned} N_j(h) = N_{j-1}\left(\frac{h}{2}\right) + \frac{N_{j-1}(h/2)-N_{j-1}(h)}{4^{j-1}-1} \end{aligned} \]

Numerical Quadrature

Integration of Lagrange Interpolating Polynomials

Select \(\{x_0,\ldots,x_n\}\) in \([a,b]\) and integrate the Lagrange polynomial \(P_n(x)=\sum_{i=0}^n f(x_i)L_i(x)\) and its truncation error term over \([a,b]\) to obtain \[ \begin{aligned} \int_a^b f(x)\, dx = \sum_{i=0}^n a_i f(x_i) + E(f) \end{aligned} \] with \[ \begin{aligned} a_i = \int_a^b L_i(x)\, dx \end{aligned} \] and \[ \begin{aligned} E(f) = \frac{1}{(n+1)!} \int_a^b \prod_{i=0}^n (x-x_i) f^{(n+1)} (\xi(x))\,dx \end{aligned} \]

Trapezoidal and Simpson’s Rules

The Trapezoidal Rule

Linear Lagrange polynomial with \(x_0=a\), \(x_1=b\), \(h=b-a\), gives \[ \begin{aligned} \int_a^b f(x)\,dx = \frac{h}{2}[f(x_0)+f(x_1)]-\frac{h^3}{12} f''(\xi) \end{aligned} \]

Simpson’s Rule

Second Lagrange polynomial with \(x_0=a\), \(x_2=b\), \(x_1=a+h\), \(h=(b-a)/2\) gives \[ \begin{aligned} \int_{x_0}^{x_2} f(x)\, dx = \frac{h}{3}[f(x_0)+4f(x_1)+f(x_2)]- \frac{h^5}{90}f^{(4)}(\xi) \end{aligned} \]

Definition

The degree of accuracy, or precision, of a quadrature formula is the largest positive integer \(n\) such that the formula is exact for \(x^k\), for each \(k=0,1,\ldots, n\).

The Newton-Cotes Formulas

The Closed Newton-Cotes Formulas

Use nodes \(x_i=x_0+ih\), \(x_0=a\), \(x_n=b\), \(h=(b-a)/n\): \[ \begin{aligned} \int_a^b f(x)\,dx &\approx \sum_{i=0}^n a_i f(x_i) \\ a_i=\int_{x_0}^{x_n} L_i(x)\,dx &= \int_{x_0}^{x_n} \prod_{j\ne i} \frac{(x-x_j)}{(x_i-x_j)}\,dx \end{aligned} \] \(n=1\) gives the Trapezoidal rule, \(n=2\) gives Simpson’s rule.

The Open Newton-Cotes Formulas

Use nodes \(x_i=x_0+ih\), \(x_0=a+h\), \(x_n=b-h\), \(h=(b-a)/(n+2)\). Setting \(n=0\) gives the Midpoint rule: \[ \begin{aligned} \int_{x_{-1}}^{x_1} f(x)\,dx = 2hf(x_0)+\frac{h^3}{3} f''(\xi) \end{aligned} \]

Composite Rules

Theorem

Let \(f\in C^2[a,b]\), \(h=(b-a)/n\), \(x_j=a+jh\), \(\mu\in(a,b)\). The Composite Trapezoidal rule for \(n\) subintervals is \[ \begin{aligned} \int_a^b f(x)\,dx &= \frac{h}{2} \left[ f(a)+2\sum_{j=1}^{n-1} f(x_j) + f(b) \right] -\frac{b-a}{12} h^2 f''(\mu) \end{aligned} \]

Theorem

Let \(f\in C^4[a,b]\), \(n\) even, \(h=(b-a)/n\), \(x_j=a+jh\), \(\mu\in(a,b)\). The Composite Simpson’s rule for \(n\) subintervals is \[ \begin{aligned} \int_a^b f(x)\,dx &= \frac{h}{3} \left[ f(a)+2\sum_{j=1}^{(n/2)-1} f(x_{2j}) + 4\sum_{j=1}^{n/2} f(x_{2j-1}) + f(b) \right] \\ &-\frac{b-a}{180} h^4 f^{(4)}(\mu) \end{aligned} \]

Romberg Integration

• Compute a sequence of \(n\) integrals using the Composite Trapezoidal rule, where \(m_1=1, m_2=2, m_3=4, \ldots\) and \(m_n=2^{n-1}\).
• The step sizes are then \(h_k=(b-a)/m_k=(b-a)/2^{k-1}\)
• The Trapezoidal rule becomes \[ \begin{aligned} \int_a^b f(x)\,dx &= \frac{h_k}{2} \left[f(a) + f(b) + 2\left(\sum_{i=1}^{2^{k-1}-1} f(a+ih_k)\right)\right] \\ & -\frac{(b-a)}{12}h_k^2 f''(\mu_k) \end{aligned} \]

Romberg Integration

• Let \(R_{k,1}\) denote the trapezoidal approximation, then \[ \begin{aligned} R_{1,1} &= \frac{h_1}{2}[f(a)+f(b)]=\frac{(b-a)}{2}[f(a)+f(b)] \\ R_{2,1} &= \frac{1}{2}[R_{1,1}+h_1f(a+h_2)] \\ R_{3,1} &= \frac{1}{2}\{R_{2,1}+h_2[f(a+h_3)+f(a+3h_3)]\} \\ R_{k,1} &= \frac{1}{2}\left[ R_{k-1,1}+h_{k-1} \sum_{i=1}^{2^{k-2}} f(a+(2i-1)h_k)\right] \end{aligned} \]
• Apply Richardson extrapolation to these values: \[ \begin{aligned} R_{k,j}=R_{k,j-1}+\frac{R_{k,j-1}-R_{k-1,j-1}}{4^{j-1}-1} \end{aligned} \]

Romberg Integration

MATLAB Implementation

function R = romberg(f, a, b, n)
% Compute integral of f(x) from a to b using Romberg integration.

h = b-a;
R = zeros(n,n);
R(1,1) = h/2 * (f(a) + f(b));
for i = 2:n
    R(i,1) = 1/2 * (R(i-1,1) + h*sum(f(a + ((1:2^(i-2))-0.5)*h)));
    for j = 2:i
        R(i,j) = R(i,j-1) + (R(i,j-1)-R(i-1,j-1)) / (4^(j-1)-1);
    end
    h = h/2;
end

Error Estimation

• The error term in Simpson’s rule requires knowledge of \(f^{(4)}\): \[ \begin{aligned} \int_a^b f(x)\,dx = S(a,b) - \frac{h^5}{90}f^{(4)}(\mu) \end{aligned} \]
• Instead, apply it again with step size \(h/2\): \[ \begin{aligned} \int_a^b f(x)\,dx = S\left(a,\frac{a+b}{2}\right) +S\left(\frac{a+b}{2},b\right) - \frac{1}{16} \left(\frac{h^5}{90}\right) f^{(4)}(\tilde{\mu}) \end{aligned} \]
• The assumption \(f^{(4)}(\mu)\approx f^{(4)}(\tilde{\mu})\) gives the error estimate \[ \begin{aligned} \left|\int_a^b f(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| \end{aligned} \]

Adaptive Quadrature

• To compute \(\int_a^b f(x)\,dx\) within a tolerance \(\varepsilon>0\), first apply Simpson’s rule with \(h=(b-a)/2\) and with \(h/2\)
• If \[ \begin{aligned} \left|S(a,b) - S\left(a,\frac{a+b}{2}\right) - S\left(\frac{a+b}{2},b\right) \right| < 15\varepsilon \end{aligned} \] then the integral is sufficiently accurate
• If not, apply the technique to \([a,(a+b)/2]\) and \([(a+b)/2,b]\), and compute the integral within a tolerance of \(\varepsilon/2\)
• Repeat until each portion is within the required tolerance

Gaussian Quadrature

• Basic idea: Calculate both nodes \(x_1,\ldots,x_n\) and coefficients \(c_1,\ldots,c_n\) such that \[ \begin{aligned} \int_a^b f(x)\,dx \approx \sum_{i=1}^n c_i f(x_i) \end{aligned} \]
• Since there are \(2n\) parameters, we might expect a degree of precision of \(2n-1\)
• Example: \(n=2\) gives the rule \[ \begin{aligned} \int_{-1}^1 f(x)\,dx \approx f\left(\frac{-\sqrt{3}}{3}\right) +f\left(\frac{\sqrt{3}}{3}\right) \end{aligned} \] with degree of precision 3

Legendre Polynomials

• The Legendre polynomials \(P_n(x)\) have the properties
  1. For each \(n\), \(P_n(x)\) is a monic polynomial of degree \(n\) (leading coefficient 1)
  2. \(\int_{-1}^1 P(x)P_n(x)\,dx=0\) when \(P(x)\) is a polynomial of degree less than \(n\)
• The roots of \(P_n(x)\) are distinct, in the interval \((-1,1)\), and symmetric with respect to the origin.
• Examples: \[ \begin{aligned} P_0(x)&=1, & P_1(x)&=x \\ P_2(x)&=x^2-\frac{1}{3} & P_3(x)&=x^3-\frac{3}{5}x \\ P_4(x)&=x^4-\frac{6}{7}x^2+\frac{3}{35} \end{aligned} \]

Gaussian Quadrature

Theorem

Suppose \(x_1,\ldots,x_n\) are roots of \(P_n(x)\) and \[ \begin{aligned} c_i = \int_{-1}^1 \prod_{j\ne i}^n \frac{x-x_j}{x_i-x_j}\,dx \end{aligned} \] If \(P(x)\) is any polynomial of degree less than \(2n\), then \[ \begin{aligned} \int_{-1}^1 P(x)\,dx = \sum_{i=1}^n c_i P(x_i) \end{aligned} \]

Computing Gaussian Quadrature Coefficients

MATLAB Implementation

function [x, c] = gaussquad(n)
% Compute Gaussian quadrature points and coefficients.

P = zeros(n+1,n+1);
P([1,2],1) = 1;
for k = 1:n-1
    P(k+2,1:k+2) = ((2*k+1)*[P(k+1,1:k+1) 0] - ...
                    k*[0 0 P(k,1:k)]) / (k+1);
end
x = sort(roots(P(n+1,1:n+1)));

A = zeros(n,n);
for i = 1:n
    A(i,:) = polyval(P(i,1:i),x)';
end
c = A \ [2; zeros(n-1,1)];

Arbitrary Intervals

Transform integrals \(\int_a^b f(x)\,dx\) into integrals over \([-1,1]\) by a change of variables: \[ \begin{aligned} t=\frac{2x-a-b}{b-a} \Leftrightarrow x=\frac{1}{2}[(b-a)t+a+b] \end{aligned} \] Gaussian quadrature then gives \[ \begin{aligned} \int_a^b f(x)\,dx = \int_{-1}^1 f\left(\frac{(b-a)t+(b+a)}{2}\right) \frac{(b-a)}{2}\,dt \end{aligned} \]

Double Integrals

• Consider the double integral \[ \begin{aligned} \iint_R f(x,y)\,dA,\qquad R=\{(x,y)\,|\,a\le x\le b,c\le y \le d \} \end{aligned} \]
• Partition \([a,b]\) and \([c,d]\) into even number of subintervals \(n,m\)
• Step sizes \(h=(b-a)/n\) and \(k=(d-c)/m\)
• Write the integral as an iterated integral \[ \begin{aligned} \iint_R f(x,y)\,dA = \int_a^b \left(\int_c^d f(x,y)\,dy\right)\,dx \end{aligned} \] and use any quadrature rule in an iterated manner.

Composite Simpson’s Rule Double Integration

The Composite Simpson’s rule gives \[ \begin{aligned} \int_a^b \left(\int_c^d f(x,y)\,dy\right)\,dx = \frac{hk}{9} \sum_{i=0}^n \sum_{j=0}^m w_{i,j} f(x_i,y_j) + E \end{aligned} \] where \(x_i=a+ih\), \(y_j=c+jk\), \(w_{i,j}\) are the products of the nested Composite Simpson’s rule coefficients (see below), and the error is \[ \begin{aligned} E=-\frac{(d-c)(b-a)}{180}\left[h^4\frac{\partial^4 f}{\partial x^4} (\bar{\eta},\bar{\mu}) + k^4\frac{\partial^4 f}{\partial y^4} (\hat{\eta},\hat{\mu})\right] \end{aligned} \]

\[ \begin{array}{rlllll} d & \bullet 1\ & \bullet 4 & \bullet 2 & \bullet 4 & \bullet 1 \\ \ \\ & \bullet 4 & \bullet 16 & \bullet 8 & \bullet 16 & \bullet 4 \\ \ \\ c & \bullet 1\ & \bullet 4 & \bullet 2 & \bullet 4 & \bullet 1 \\ & a & & & & b \end{array} \]

Non-Rectangular Regions

The same technique can be applied to double integrals of the form \[ \begin{aligned} \int_a^b \int_{c(x)}^{d(x)} f(x,y)\,dy\,dx \end{aligned} \] The step size for \(x\) is still \(h=(b-a)/n\), but for \(y\) it varies with \(x\): \[ \begin{aligned} k(x) = \frac{d(x)-c(x)}{m} \end{aligned} \]

Gaussian Double Integration

• For Guassian integration, first transform the roots \(r_{n,j}\) from \([-1,1]\) to \([a,b]\) and \([c,d]\), respectively
• The integral is then \[ \begin{aligned} \int_a^b \int_c^d f(x,y)\,dy\,dx \approx \frac{(b-a)(d-c)}{4}\sum_{i=1}^n \sum_{j=1}^n c_{n,i} c_{n,j} f(x_i,y_j) \end{aligned} \]
• Similar techniques can be used for non-rectangular regions

Improper Integrals with a Singularity

The improper integral below, with a singularity at the left endpoint, converges if and only if \(0<p<1\) and then \[ \begin{aligned} \int_a^b \frac{1}{(x-a)^p}dx = \left.\frac{(x-a)^{1-p}}{1-p}\right|_a^b =\frac{(b-a)^{1-p}}{1-p} \end{aligned} \] More generally, if \[ \begin{aligned} f(x)=\frac{g(x)}{(x-a)^p},\quad 0<p<1,\quad g\text{ continuous on }[a,b], \end{aligned} \] construct the fourth Taylor polynomial \(P_4(x)\) for \(g\) about \(a\): \[ \begin{aligned} P_4(x)&=g(a)+g'(a)(x-a)+\frac{g''(a)}{2!}(x-a)^2 \\ &+\frac{g'''(a)}{3!}(x-a)^3 + \frac{g^{(4)}(a)}{4!}(x-a)^4 \end{aligned} \]

Improper Integrals with a Singularity

and write \[ \begin{aligned} \int_a^b f(x)\,dx = \int_a^b \frac{g(x)-P_4(x)}{(x-a)^p}dx +\int_a^b \frac{P_4(x)}{(x-a)^p}dx \end{aligned} \] The second integral can be computed exactly: \[ \begin{aligned} \int_a^b \frac{P_4(x)}{(x-a)^p}dx &= \sum_{k=0}^4 \int_a^b \frac{g^{(k)}(a)}{k!}(x-a)^{k-p}\,dx \\ &= \sum_{k=0}^4 \frac{g^{(k)}(a)}{k!(k+1-p)}(b-a)^{k+1-p} \end{aligned} \]

Improper Integrals with a Singularity

For the first integral, use the Composite Simpson’s rule to compute the integral of \(G\) on \([a,b]\), where \[ \begin{aligned} G(x)= \begin{cases} \frac{g(x)-P_4(x)}{(x-a)^p}, & \text{if }a<x\le b \\ 0, & \text{if }x=a \end{cases} \end{aligned} \] Note that \(0<p<1\) and \(P_4^{(k)}(a)\) agrees with \(g^{(k)}(a)\) for each \(k=0,1,2,3,4\), so \(G\in C^4[a,b]\) and Simpson’s rule can be applied.

Singularity at the Right Endpoint

• For an improper integral with a singularity at the right endpoint \(b\), make the substitution \(z=-x\), \(dz=-dx\) to obtain \[ \begin{aligned} \int_a^b f(x)\,dx = \int_{-b}^{-a} f(-z)\,dz \end{aligned} \] which has its singularity at the left endpoint
• For an improper integral with a singularity at \(c\), where \(a<c<b\), split into two improper integrals \[ \begin{aligned} \int_a^b f(x)\,dx = \int_a^c f(x)\,dx + \int_c^b f(x)\,dx \end{aligned} \]

Infinite Limits of Integration

An integral of the form \(\int_a^\infty \frac{1}{x^p}dx\), with \(p>1\), can be converted to an integral with left endpoint singularity at \(0\) by the substitution \[ \begin{aligned} t=x^{-1},\quad dt=-x^{-2}\,dx, \text{ so }dx=-x^2 dt=-t^{-2}dt \end{aligned} \] which gives \[ \begin{aligned} \int_a^\infty \frac{1}{x^p} dx = \int_{1/a}^0 -\frac{t^p}{t^2}dt=\int_0^{1/a}\frac{1}{t^{2-p}}dt \end{aligned} \] More generally, this variable change converts \(\int_a^\infty f(x)\,dx\) into \[ \begin{aligned} \int_a^\infty f(x)\,dx = \int_0^{1/a} t^{-2} f\left(\frac{1}{t}\right)\,dt \end{aligned} \]