Numerical Differentiation and Integration

5.1 Numerical Differentiation

1) Numerical Differentiation

Numerical Differentiation: The approximation of a derivative using function values at a finite set of points.

It is required when a function is given by measured data, tabulated values, or a simulation instead of an explicit formula.

2) Finite Difference

Finite Difference: A difference between function values used to approximate a derivative.

The derivative definition is

$$f'(x) = \lim_{h\to0} \frac{f(x+h)-f(x)}{h}.$$

Numerical differentiation replaces the limit with a small finite step size $h$.

3) Forward-Difference Formula

Forward Difference: A first-derivative approximation using $f(x)$ and $f(x+h)$.

$$f'(x) \approx \frac{f(x+h)-f(x)}{h}.$$

Its error formula is

$$f'(x) = \frac{f(x+h)-f(x)}{h} - \frac{h}{2}f''(c),$$

where $c\in(x,x+h)$.

Therefore, the truncation error is

$$O(h).$$

The forward difference is a first-order method.

4) Backward-Difference Formula

Backward Difference: A first-derivative approximation using $f(x)$ and $f(x-h)$.

$$f'(x) \approx \frac{f(x)-f(x-h)}{h}.$$

Its truncation error is

$$O(h).$$

Backward differences are useful near the right endpoint of an interval.

5) Centered-Difference Formula

Centered Difference: A first-derivative approximation using points on both sides of $x$.

$$f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}.$$

The error formula is

$$f'(x) = \frac{f(x+h)-f(x-h)}{2h} - \frac{h^2}{6}f'''(c),$$

where $c\in(x-h,x+h)$.

Therefore, the truncation error is

$$O(h^2).$$

The centered formula is more accurate than the forward or backward formula for the same step size.

6) Second-Derivative Formula

The three-point centered approximation to the second derivative is

$$f''(x) \approx \frac{ f(x-h)-2f(x)+f(x+h) }{h^2}.$$

Its error formula is

$$f''(x) = \frac{ f(x-h)-2f(x)+f(x+h) }{h^2} - \frac{h^2}{12}f^{(4)}(c).$$

Therefore, the truncation error is

$$O(h^2).$$

7) Order of Accuracy

Order of Accuracy: A numerical approximation has order $p$ if its truncation error satisfies

$$E(h)=O(h^p).$$

This means that, for sufficiently small $h$,

$$|E(h)|\leq Ch^p$$

for some constant $C$ independent of $h$.

If $h$ is halved, an order-$p$ error is approximately reduced by

$$2^p.$$

8) Truncation Error

Truncation Error: The error caused by replacing an exact mathematical process with a finite approximation.

For the centered-difference formula,

$$E_{\mathrm{trunc}} = O(h^2).$$

Truncation error generally decreases as $h$ becomes smaller.

9) Rounding Error in Differentiation

Finite-difference formulas subtract nearly equal function values.

For very small $h$, this causes cancellation and loss of significant digits.

The rounding-error contribution is approximately

$$E_{\mathrm{round}} = O\left( \frac{\epsilon_{\mathrm{mach}}}{h} \right).$$

Rounding error increases as $h$ becomes smaller.

10) Total Differentiation Error

For the centered-difference formula, the total error can be modeled by

$$E(h) \approx C_1h^2 + C_2\frac{\epsilon_{\mathrm{mach}}}{h}.$$

The first term is truncation error, and the second is rounding error.

Therefore, decreasing $h$ indefinitely does not indefinitely improve accuracy.

11) Optimal Step Size

For the centered-difference approximation, if

$$M\approx|f'''(x)|,$$

the estimated optimal step size is

$$h_{\mathrm{opt}} \approx \left( \frac{3\epsilon_{\mathrm{mach}}}{M} \right)^{1/3}.$$

Ignoring the scale factor, this is approximately

$$h_{\mathrm{opt}} \sim \epsilon_{\mathrm{mach}}^{1/3}.$$

For double precision, this is typically near

$$10^{-5}.$$

12) Richardson Extrapolation

Richardson Extrapolation: A method that combines approximations with different step sizes to cancel the leading error term.

Suppose

$$Q = F(h)+Ch^p+O(h^{p+1}).$$

Then

$$Q = F(h/2) + C\frac{h^p}{2^p} + O(h^{p+1}).$$

Eliminating the $Ch^p$ term gives

$$Q \approx \frac{ 2^pF(h/2)-F(h) }{ 2^p-1 }.$$

The extrapolated approximation has a higher order than the original formula.

13) Five-Point Centered Difference

Applying extrapolation to the second-order centered difference gives

$$f'(x) \approx \frac{ f(x-h) -8f(x-h/2) +8f(x+h/2) -f(x+h) }{ 6h }.$$

This formula has truncation error

$$O(h^4).$$

An equivalent form using points spaced by $h$ is

$$f'(x) \approx \frac{ f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h) }{ 12h }.$$

5.2 Newton-Cotes Formulas for Numerical Integration

1) Numerical Integration

Numerical Integration: The approximation of a definite integral using a weighted sum of function values.

A quadrature rule has the form

$$\int_a^b f(x)\,dx \approx \sum_{i=0}^n w_if(x_i),$$

where:

• $x_i$ are quadrature nodes.
• $w_i$ are quadrature weights.

2) Quadrature

Quadrature: Another name for numerical integration.

Quadrature rules approximate the area under a function using values at selected points.

3) Newton-Cotes Formula

Newton-Cotes Formula: A quadrature formula obtained by integrating an interpolating polynomial through equally spaced nodes.

A closed Newton-Cotes formula includes both endpoints.

An open Newton-Cotes formula excludes the endpoints.

4) Trapezoid Rule

Trapezoid Rule: A quadrature rule obtained by integrating the line through $(a,f(a))$ and $(b,f(b))$.

Let

$$h=b-a.$$

Then

$$\int_a^b f(x)\,dx \approx \frac{h}{2} \left[ f(a)+f(b) \right].$$

The error is

$$E_T = -\frac{h^3}{12}f''(c)$$

for some $c\in(a,b)$.

The Trapezoid Rule has degree of precision $1$.

5) Simpson’s Rule

Simpson’s Rule: A quadrature rule obtained by integrating the quadratic interpolating polynomial through the endpoints and midpoint.

Let

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

Then

$$\int_a^b f(x)\,dx \approx \frac{h}{3} \left[ f(a) + 4f(a+h) + f(b) \right].$$

Equivalently,

$$\int_a^b f(x)\,dx \approx \frac{b-a}{6} \left[ f(a) + 4f\left(\frac{a+b}{2}\right) + f(b) \right].$$

The error is

$$E_S = -\frac{h^5}{90}f^{(4)}(c)$$

for some $c\in(a,b)$.

Simpson’s Rule has degree of precision $3$.

6) Midpoint Rule

Midpoint Rule: An open Newton-Cotes formula that evaluates the function at the interval midpoint.

$$\int_a^b f(x)\,dx \approx (b-a) f\left( \frac{a+b}{2} \right).$$

Its error is

$$E_M = \frac{(b-a)^3}{24}f''(c)$$

for some $c\in(a,b)$.

The Midpoint Rule has degree of precision $1$.

7) Degree of Precision

Degree of Precision: The largest integer $k$ such that a quadrature method integrates every polynomial of degree at most $k$ exactly.

For example:

• Trapezoid Rule: degree $1$
• Midpoint Rule: degree $1$
• Simpson’s Rule: degree $3$

8) Composite Quadrature

Composite Quadrature: Applying a quadrature rule separately over several subintervals and adding the results.

Composite rules improve accuracy without requiring a single high-degree interpolating polynomial.

9) Composite Trapezoid Rule

Divide $[a,b]$ into $m$ equal subintervals:

$$h=\frac{b-a}{m},$$ $$x_i=a+ih.$$

Then

$$\int_a^b f(x)\,dx \approx h \left[ \frac{1}{2}f(x_0) + \sum_{i=1}^{m-1}f(x_i) + \frac{1}{2}f(x_m) \right].$$

The error is

$$E_{CT} = -\frac{(b-a)h^2}{12}f''(c).$$

Therefore, the Composite Trapezoid Rule has global error

$$O(h^2).$$

10) Composite Simpson’s Rule

Divide $[a,b]$ into $2m$ equal subintervals:

$$h=\frac{b-a}{2m}.$$

Then

$$\begin{aligned} \int_a^b f(x)\,dx \approx \frac{h}{3} \Bigg[ &f(x_0) + 4\sum_{i=1}^{m}f(x_{2i-1})\\ &+ 2\sum_{i=1}^{m-1}f(x_{2i}) + f(x_{2m}) \Bigg]. \end{aligned}$$

The error is

$$E_{CS} = -\frac{(b-a)h^4}{180}f^{(4)}(c).$$

Therefore, the Composite Simpson’s Rule has global error

$$O(h^4).$$

11) Composite Midpoint Rule

Divide $[a,b]$ into $m$ panels of width

$$h=\frac{b-a}{m}.$$

Then

$$\int_a^b f(x)\,dx \approx h \sum_{i=0}^{m-1} f\left( a+\left(i+\frac{1}{2}\right)h \right).$$

Its global error is

$$O(h^2).$$

5.3 Romberg Integration

1) Romberg Integration

Romberg Integration: A method that applies Richardson extrapolation repeatedly to Composite Trapezoid Rule approximations.

It produces a triangular table of increasingly accurate approximations.

2) First Romberg Column

Let

$$h_j=\frac{b-a}{2^{j-1}}.$$

The first-column value

$$R_{j,1}$$

is the Composite Trapezoid Rule using step size $h_j$.

The first value is

$$R_{1,1} = \frac{b-a}{2} \left[ f(a)+f(b) \right].$$

3) Incremental Trapezoid Update

The next trapezoid approximation can reuse the previous function values:

$$R_{j,1} = \frac{1}{2}R_{j-1,1} + h_j \sum_{i=1}^{2^{j-2}} f\left( a+(2i-1)h_j \right).$$

Only newly introduced midpoint values must be evaluated.

4) Romberg Extrapolation

Higher columns are computed by

$$R_{j,k} = R_{j,k-1} + \frac{ R_{j,k-1}-R_{j-1,k-1} }{ 4^{k-1}-1 },$$

where

$$2\leq k\leq j.$$

Equivalently,

$$R_{j,k} = \frac{ 4^{k-1}R_{j,k-1} - R_{j-1,k-1} }{ 4^{k-1}-1 }.$$

5) Romberg Error Order

The first column has error

$$O(h^2).$$

Successive columns cancel the leading even powers of $h$:

• First column: $O(h^2)$
• Second column: $O(h^4)$
• Third column: $O(h^6)$
• $k$th column: $O(h^{2k})$

6) Romberg Stopping Criterion

A common stopping criterion is

$$|R_{j,j}-R_{j-1,j-1}| < \mathrm{TOL}.$$

The difference between successive diagonal entries estimates convergence.

5.4 Adaptive Quadrature

1) Adaptive Quadrature

Adaptive Quadrature: A numerical integration method that automatically uses smaller subintervals where the integrand is difficult and larger subintervals where it is smooth.

It concentrates computation where the estimated error is largest.

2) Local Error Estimation

Let $S[a,b]$ be a quadrature approximation on $[a,b]$.

Split the interval at

$$c=\frac{a+b}{2}$$

and compute

$$S[a,c]+S[c,b].$$

The difference

$$S[a,b] - \left( S[a,c]+S[c,b] \right)$$

estimates the local error.

3) Adaptive Trapezoid Rule

The Trapezoid Rule has second-order composite error.

When an interval is divided into two halves, the error is approximately reduced by a factor of $4$.

The difference between the coarse and refined approximations is approximately three times the refined error.

Thus,

$$E_{\mathrm{refined}} \approx \frac{ S[a,c]+S[c,b]-S[a,b] }{3}.$$

Accept the refined result when its estimated error is below the local tolerance.

4) Adaptive Simpson’s Rule

Let

$$S=S[a,b]$$

be Simpson’s Rule over the full interval, and let

$$S_2=S[a,c]+S[c,b]$$

be the sum of Simpson’s Rule over the two halves.

Because Simpson’s Rule is fourth order,

$$E_{\mathrm{refined}} \approx \frac{S_2-S}{15}.$$

A common acceptance condition is

$$|S_2-S| < 15\,\mathrm{TOL}.$$

A more conservative implementation may use

$$|S_2-S| < 10\,\mathrm{TOL}.$$

5) Corrected Simpson Estimate

When the interval is accepted, Richardson correction gives

$$I \approx S_2 + \frac{S_2-S}{15}.$$

This is more accurate than using $S_2$ alone.

6) Recursive Subdivision

If the estimated error is too large:

1. Divide $[a,b]$ into $[a,c]$ and $[c,b]$.
2. Assign half the tolerance to each subinterval.
3. Apply the same test recursively.
4. Add the accepted subinterval results.

Adaptive quadrature should include a maximum recursion depth or minimum interval width.

7) Advantages of Adaptive Quadrature

Adaptive quadrature is effective for functions with:

• Localized peaks
• Rapid changes
• High curvature
• Uneven smoothness
• Expensive function evaluations

It avoids wasting evaluations in smooth regions.

5.5 Gaussian Quadrature

1) Orthogonal Functions

Orthogonal Functions: Nonzero functions $p_j$ and $p_k$ satisfying

$$\int_a^b p_j(x)p_k(x)\,dx = 0 \qquad\text{for }j\neq k.$$

A set of functions is orthogonal if every distinct pair is orthogonal.

2) Orthogonal Polynomials

Orthogonal Polynomials: Polynomials that form an orthogonal set under an integral inner product.

If

$$\deg(p_i)=i,$$

then

$$\{p_0,p_1,\ldots,p_n\}$$

forms a basis for the polynomial space of degree at most $n$.

3) Legendre Polynomials

Legendre Polynomials: A sequence of polynomials orthogonal on $[-1,1]$.

They can be defined by Rodrigues’ formula:

$$P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} \left[ (x^2-1)^n \right].$$

The first Legendre polynomials are

$$P_0(x)=1,$$ $$P_1(x)=x,$$ $$P_2(x)=\frac{1}{2}(3x^2-1),$$ $$P_3(x)=\frac{1}{2}(5x^3-3x).$$

4) Legendre Roots

The degree-$n$ Legendre polynomial has $n$ distinct roots in

$$(-1,1).$$

These roots are used as the nodes of Gauss-Legendre Quadrature.

5) Gauss-Legendre Quadrature

Gauss-Legendre Quadrature: A quadrature rule using the roots of a Legendre polynomial as optimally selected nodes.

On $[-1,1]$,

$$\int_{-1}^{1}f(x)\,dx \approx \sum_{i=1}^n w_if(x_i),$$

where:

• $x_i$ are the roots of $P_n(x)$.
• $w_i$ are the corresponding quadrature weights.

6) Gaussian Weights

If $L_i(x)$ is the Lagrange basis polynomial associated with node $x_i$, then

$$w_i = \int_{-1}^{1}L_i(x)\,dx.$$

The weights are positive and symmetric.

7) Two-Point Gauss-Legendre Rule

For $n=2$,

$$x_1=-\frac{1}{\sqrt{3}}, \qquad x_2=\frac{1}{\sqrt{3}},$$

and

$$w_1=w_2=1.$$

Therefore,

$$\int_{-1}^{1}f(x)\,dx \approx f\left(-\frac{1}{\sqrt{3}}\right) + f\left(\frac{1}{\sqrt{3}}\right).$$

8) Three-Point Gauss-Legendre Rule

For $n=3$,

$$x_1=-\sqrt{\frac{3}{5}}, \qquad x_2=0, \qquad x_3=\sqrt{\frac{3}{5}},$$

and

$$w_1=\frac{5}{9}, \qquad w_2=\frac{8}{9}, \qquad w_3=\frac{5}{9}.$$

Thus,

$$\begin{aligned} \int_{-1}^{1}f(x)\,dx \approx{}& \frac{5}{9} f\left(-\sqrt{\frac{3}{5}}\right)\\ &+ \frac{8}{9}f(0) + \frac{5}{9} f\left(\sqrt{\frac{3}{5}}\right). \end{aligned}$$

9) Degree of Precision of Gaussian Quadrature

An $n$-point Gauss-Legendre rule has degree of precision

$$2n-1.$$

It integrates every polynomial of degree at most $2n-1$ exactly.

This is the highest possible degree of precision for a quadrature rule using $n$ nodes.

10) Change of Interval

To transform an integral over $[a,b]$ to $[-1,1]$, use

$$x = \frac{a+b}{2} + \frac{b-a}{2}t.$$

Then

$$dx = \frac{b-a}{2}\,dt,$$

so

$$\int_a^b f(x)\,dx = \frac{b-a}{2} \int_{-1}^{1} f\left( \frac{a+b}{2} + \frac{b-a}{2}t \right)dt.$$

The Gaussian approximation becomes

$$\int_a^b f(x)\,dx \approx \frac{b-a}{2} \sum_{i=1}^n w_i f\left( \frac{a+b}{2} + \frac{b-a}{2}x_i \right).$$

Comparison of Numerical Integration Methods

Method	Nodes	Degree of Precision	Global Error
Composite Midpoint	Midpoints	1	$O(h^2)$
Composite Trapezoid	Equally spaced	1	$O(h^2)$
Composite Simpson	Equally spaced	3	$O(h^4)$
Romberg	Reused equally spaced nodes	Increasing	$O(h^{2k})$ in column $k$
Adaptive Simpson	Adaptive nodes	3 locally	Controlled by tolerance
$n$-point Gaussian	Legendre roots	$2n-1$	Depends on higher derivatives

Essential Concepts

1. Numerical Differentiation: Approximates derivatives using finite function differences.
2. Forward Difference: A first-order one-sided derivative formula.
3. Centered Difference: A second-order symmetric derivative formula.
4. Second-Derivative Difference: Uses three centered function values.
5. Order of Accuracy: Describes how quickly truncation error decreases with $h$.
6. Truncation Error: Results from replacing an exact process with a finite formula.
7. Rounding Error: Grows in differentiation when nearly equal values are subtracted.
8. Optimal Step Size: Balances truncation and rounding errors.
9. Richardson Extrapolation: Cancels the leading error term using two step sizes.
10. Quadrature: Numerical approximation of a definite integral.
11. Newton-Cotes Formula: Integrates a polynomial through equally spaced nodes.
12. Trapezoid Rule: Integrates a linear interpolant.
13. Simpson’s Rule: Integrates a quadratic interpolant.
14. Degree of Precision: The highest polynomial degree integrated exactly.
15. Composite Rule: Applies a quadrature rule across multiple panels.
16. Romberg Integration: Repeatedly extrapolates Composite Trapezoid approximations.
17. Adaptive Quadrature: Refines only intervals with large estimated errors.
18. Orthogonal Polynomial: Is orthogonal to other basis polynomials under an integral inner product.
19. Legendre Polynomial: Provides the nodes for Gauss-Legendre Quadrature.
20. Gaussian Quadrature: Achieves degree of precision $2n-1$ using $n$ nodes.