Numerical Integration Notes

Numerical Integration

Goal and Approach

• Goal: Approximate the definite integral of a function f(x) from a to b numerically: \int_{a}^{b} f(x) dx, where f : \mathbb{R} \rightarrow \mathbb{R}.
• Quadrature Formula: Approximate the definite integral using a weighted sum of function values at specific points: \int{a}^{b} f(x) dx \approx \omega0 f(x0) + \cdots + \omegan f(x_n)
  • x0, \dots, xn are the quadrature nodes.
  • \omega0, \dots, \omegan are the quadrature weights.
• Definition: A quadrature rule is exact for polynomials of degree n if it gives the exact value when f(x) is a polynomial of degree less than or equal to n.

Change of Variables

• Transformation: Use the change of variable y = \frac{x-a}{b-a} to transform the integral over the interval [a, b] to an integral over [0, 1]:
  \int{a}^{b} f(x) dx = (b - a) \int{0}^{1} g(y) dy, where g(y) = f(a + (b - a)y).


• Interval Transformation: Quadrature formulas on [0, 1] can be transformed to any interval [a, b]:
  

• Idea: Replace f(x) with an interpolating polynomial Pn(x) using n+1 equispaced nodes in [a, b]: Pn(x) = \sum{k=0}^{n} f(xk) \ellk(x), where \ellk(x) = \prod{j \neq k} \frac{x - xj}{xk - xj}, and x_k = a + k \cdot \frac{b - a}{n}.
• Closed Quadrature Newton-Cotes Formulas: Approximate the integral as:
  \int{a}^{b} f(x) dx \approx \sum{k=0}^{n} f(xk) \alphak, where \alphak = \int{a}^{b} \ell_k(x) dx.

Weights of the NC Quadrature Rules

• Lemma: The coefficients \\alphak are unique such that \int{a}^{b} p(x) dx = \sum{k=0}^{n} p(xk) \alpha_k for any polynomial p(x) of degree less than or equal to n.
• Calculating Weights: To find the values of \\alphak, take pj(x) = x^j for j = 0, \dots, n. Since \int{0}^{1} x^j dx = \frac{1}{j+1} for j = 0, \dots, n, we solve the following linear system: \begin{aligned} \alpha0 + \alpha1 + \cdots + \alpha{n-1} + \alphan &= \frac{1}{1} \ \frac{1}{n}\alpha1 + \cdots + \frac{n-1}{n} \alpha{n-1} + \alphan &= \frac{1}{2} \
  \vdots \
  \left(\frac{1}{n}\right)^n \alpha1 + \cdots + \left(\frac{n-1}{n}\right)^n \alpha{n-1} + \alpha_n &= \frac{1}{n+1}
  \end{aligned}

Basic Closed NC Quadrature Rules for n = 1, 2

• Trapezoid's Rule (n = 1):
  \int_{a}^{b} f(x) dx \approx (b - a) \cdot \frac{f(a) + f(b)}{2}
• Simpson's Rule (n = 2):
  \int_{a}^{b} f(x) dx \approx \frac{b - a}{6} \cdot \left[ f(a) + 4 \cdot f\left(\frac{a+b}{2}\right) + f(b) \right]

Error in the Basic NC Formulas

• Recall: The error in polynomial interpolation:
  f(x) = Pn(x) + \frac{\prod{i=0}^{n} (x - x_i) f^{(n+1)}(\xi(x))}{(n+1)!}, where a \leq \xi(x) \leq b.
• Integrating: Integrating the error term gives:
  \int{a}^{b} f(x) dx = \int{a}^{b} \left[ Pn(x) + \frac{\prod{i=0}^{n} (x - xi) f^{(n+1)}(\xi(x))}{(n+1)!} \right] dx = \int{a}^{b} Pn(x) dx + \int{a}^{b} \frac{\prod{i=0}^{n} (x - xi) f^{(n+1)}(\xi(x))}{(n+1)!} dx
  = \sum{k=0}^{n} f(xk) \alphak + \int{a}^{b} \frac{\prod{i=0}^{n} (x - xi) f^{(n+1)}(\xi(x))}{(n+1)!} dx.
• Error Bound: The error is thus bounded by
  En(f) := \left| \int{a}^{b} f(x) dx - \sum{k=0}^{n} f(xk) \alphak \right| = \left| \int{a}^{b} \frac{\prod{i=0}^{n} (x - xi) f^{(n+1)}(\xi(x))}{(n + 1)!} dx \right|
  \leq \int{a}^{b} \frac{\prod{i=0}^{n} |x - xi| |f^{(n+1)}(\xi(x))|}{(n + 1)!} dx \leq \frac{(b - a)^{n+1} \max{a \leq x \leq b} |f^{(n+1)}(x)|}{(n + 1)!} \int_{a}^{b} 1 dx.
• Error Bound for Simple Closed Newton-Cotes Formula: En(f) \leq \frac{(b - a)^{n+2}}{(n+1)!} \cdot \max{a \leq x \leq b} |f^{(n+1)}(x)| .
• Even n Error: If n is even, then E{n+1}(f) \leq cn \cdot (b - a)^{n+3} \cdot \max{a \leq x \leq b} |f^{(n+2)}(x)|, where cn = \frac{1}{(n+2)!} \left| \int_{0}^{n} t^2 (t - 1) \cdots (t - n) dt \right|.

Composite NC Formulas: Basic Setting

• Idea: Split the interval [a, b] into subintervals [a, x1] \cup [x1, x2] \cup \cdots \cup [x{p-1}, b], integrate over each subinterval, and apply a Newton-Cotes formula with n+1 nodes to each.
  \int{a}^{b} f(x) dx = \int{a}^{x1} f(x) dx + \int{x1}^{x2} f(x) dx + \cdots + \int{x{p-1}}^{b} f(x) dx,
• Ingredients:
  • N = np (total number of intervals).
  • h_N = \frac{b - a}{N} (integration step).
  • xj = a + j \cdot hN for j = 0, \dots, N (integration nodes).

Composite NC Formulas for n = 1, 2

• Composite Trapezoid's Rule (n = 1 \implies N = p):
  \int{a}^{b} f(x) dx \approx \frac{hN}{2} \left( f(a) + 2f(x1) + \cdots + 2f(x{N-1}) + f(b) \right)
• Composite Simpson's Rule (n = 2 \implies N = 2p):
  \int{a}^{b} f(x) dx \approx \frac{hN}{3} \left( f(a) + 4f(x1) + 2f(x2) + 4f(x3) + \cdots + 2f(x{N-2}) + 4f(x{N-1}) + f(b) \right) \int{a}^{b} f(x) dx \approx \frac{h_N}{3} (E + 4I + 2P)

Error in Composite NC Formulas

• Error Bound: For n odd:
  E(f) \leq (b - a) \left( \frac{b - a}{N} \right)^{n+1} \cdot \frac{n^{n+1}}{(n+1)!} \cdot \max_{a \leq x \leq b} |f^{(n+1)}(x)|
• Error Bound: For n even:
  E(f) \leq cn \cdot (b - a) \left( \frac{b - a}{N} \right)^{n+2} \cdot n^{n+2} \cdot \max{a \leq x \leq b} |f^{(n+2)}(x)|
• Note: n is fixed. This method can be too expensive.

Richardson Extrapolation

• Concept: Given an approximation A(h) to some quantity A, dependent on a stepsize h:
  A(h) = A + c h^n + O(h^{n+1}) (where c may depend on n but not on h).
• Example: NC for n odd is an approximation to A = \int_{a}^{b} f(x) dx with order n+1 on h = \frac{b - a}{N}.
• Richardson Extrapolation Formula:
  R(h, t) := \frac{t^n A(h/t) - A(h)}{t^n - 1}
  has convergence order n+1, since
  R(h, t) = \frac{t^n \left( A + c \left( \frac{h}{t} \right)^n + O(h^{n+1}) \right) - (A + c h^n + O(h^{n+1}))}{t^n - 1} = A + O(h^{n+1}).

Extrapolated Simpson Rule

• Apply Richardson extrapolation to the simple Simpson rule, S(h), with t = 2 in [a, b], where h = b - a.
• Basic Simpson's Rule: A(h) = S(h) = \frac{h}{6} \left[ f(a) + 4f \left( \frac{a+b}{2} \right) + f(b) \right]
• Composite Simpson's Rule: A(\frac{h}{2}) = S_2(h) = \frac{h}{12} \left[ f(a) + 4f \left( \frac{3a+b}{4} \right) + 2f \left( \frac{a+b}{2} \right) + 4f \left( \frac{a+3b}{4} \right) + f(b) \right]
• Extrapolated Simpson Rule: Q(h) = \frac{16 S_2(h) - S(h)}{15} = \frac{h}{90} \left[ 7 f(a) + 32 f \left( \frac{3a + b}{4} \right) + 12 f \left( \frac{a + b}{2} \right) + 32 f \left( \frac{a + 3b}{4} \right) + 7 f(b) \right].
• Let h = b - a:
  Q(h) = \frac{h}{90} \left[ 7 f(a) + 32 f \left(a + \frac{h}{4} \right) + 12 f \left(a + \frac{h}{2} \right) + 32 f \left(a + \frac{3h}{4} \right) + 7 f(b) \right]
• This is the simple closed Newton-Cotes formula for 5 nodes and is exact for polynomials with degree less than or equal to 5.

Adaptive Integration

• Goal: Obtain a good approximation to \int_{a}^{b} f(x) dx while minimizing the number of function evaluations.
• Idea: Break [a, b] into pieces. Approximate the integral in each piece using two different methods (e.g., simple and composite formulas). If the difference is small enough, use Richardson extrapolation. Otherwise, subdivide the interval and repeat.
• Steps:
  1. Compute S and S_2 (simple and composite Simpson's rule with 2 subintervals).
  2. Evaluate E = |S - S_2|, which estimates the error.
  3. If E \leq tol, then \int{a}^{b} f(x) dx \approx \frac{16S2 - S}{15} = Q.
  4. If E > tol, repeat steps 1 and 2 in:
     • Step 3.1: the interval \left[a, \frac{a+b}{2} \right], and
     • Step 3.2: the interval \left[\frac{a+b}{2}, b\right].
  5. Add up the values of Q obtained in all steps at each iteration.

MATLAB Commands

• quad: MATLAB's professional method for approximating integrals, using an extrapolated Simpson's rule in an adaptive method.
• quadtx: version of quad in the NCM book/folder.

Example

• f(x) = e^{-4x} \sin(2x) in [0, 4]

Exercise

• I := 4 \int_{0}^{1} \frac{1}{1+x^2} dx.
  1. What is the value of I?
  2. Approximate I using the trapezoid rule.
  3. Approximate I using Simpson's rule.
  4. Approximate I using composite Simpson's rule with 8 intervals.
  5. How many subintervals are needed to approximate I to within 10^{-16} using the composite Simpson's formula?