Numerical Differentiation and Integration

Numerical Differentiation

Basic Problems

  • Problem 1: Derive a formula that approximates the derivative of a function using a linear combination of function values. This is applicable when the function is known.
  • Problem 2: Approximate the value of a derivative of a function defined by discrete data.

Solution Approaches

  • Use Taylor Series Expansion.
  • Pass a polynomial through the given data and differentiate the interpolating polynomial.

Applications

  • Solving Ordinary and Partial Differential Equations.

First Derivative

  • For a function f : [a, b] \rightarrow R, the derivative f'(c) at a point c is defined as: f'(c) = \lim_{h \to 0} \frac{f(c + h) - f(c)}{h}, for all c \in (a, b)
  • Geometrically, f'(c) represents the slope of the tangent to the curve f(x) at x = c.

Taylor Series Derivative

  • Taylor's Theorem: If f has two continuous derivatives, then: f(x0 + h) = f(x0) + hf'(x0) + \frac{h^2}{2}f''(\theta), where \theta \in (x0, x_0 + h)
  • Forward Formula: f'(x0) \approx \frac{f(x0 + h) - f(x_0)}{h}

Absolute Error

  • Error is defined as: |\text{true value} - \text{approximate value}|
  • Error bound for the Forward Formula: |ED(f)| \leq \max{\theta \in [a, b]} \frac{h}{2} |f''(\theta)|
  • Example: Find the derivative of f(x) = x^2 at x = 1 with h = 0.1 using the Taylor series.

Other Formulae

  • Backward Formula: f'(x0) = \frac{f(x0) - f(x_0 - h)}{h}
  • Central Formula: f'(x0) = \frac{f(x0 + h) - f(x_0 - h)}{2h}
  • Second Derivative Approximation: f''(x0) \approx \frac{f(x0 + h) - 2f(x0) + f(x0 - h)}{h^2}

Derivative for Discrete Data

  • Given data points: (x0, f(x0)), (x1, f(x1)), (x2, f(x2)), …, (xn, f(xn))
  • Assumption: Data points are equispaced, i.e., xi - x{i-1} = h
  • Remark 1: Newton's divided difference (or Forward/Backward) formula can be used, depending on the point's location.
  • Remark 2: If data is not equispaced, Lagrange interpolating polynomial can be used.

Newton's Forward Difference Formula

  • Formula: f(x) \approx Pn(x) = f(x0) + \sum{k=1}^{n} \binom{s}{k} \Delta^k f(x0), where x = x_0 + sh
  • We use Pn(x) to calculate the derivatives of f, i.e., f'(x) \approx P'n(x) for all x \in [x0, xn]

Derivative Calculation Using Newton's Forward Difference

  • First derivative: f'(x) \approx P'(x) = \frac{dPn}{dx} = \frac{1}{h} \frac{dPn}{ds}
  • Second derivative: f''(x) \approx \frac{d^2Pn}{dx^2} = \frac{1}{h^2} \frac{d^2Pn}{ds^2}
  • In general: f^{(K)}(x) = \frac{1}{h^K} \frac{d^KP_n}{ds^K}
  • Example: Compute the first and second derivatives at x = 2 for the function defined by the following data using Taylor series expansion and Newton's forward divided difference.
    • x: 1, 2, 3, 4
    • f(x): 2, 5, 7, 10

Solution using Taylor Series

  • Here h = 1
  • f'(2) = \frac{f(2 + h) - f(2)}{h} = \frac{f(3) - f(2)}{1} = 2
  • f''(2) = \frac{f(2 + h) - 2f(2) + f(2 - h)}{h^2} = \frac{f(3) - 2f(2) + f(1)}{1} = 1

Solution using Newton's Forward Difference Formula

  • Difference table:

    • x: 1, 2, 3, 4
    • f(x): 2, 5, 7, 10
    • \Delta f: 3, 2, 3
    • \Delta^2 f: -1, 1
    • \Delta^3 f: 2

  • f(x) \approx Pn(x) = f(x0) + s\Delta f(x0) + \frac{s(s - 1)}{2!} \Delta^2 f(x0) + \frac{s(s - 1)(s - 2)}{3!} \Delta^3 f(x_0)

  • f'(x) \approx \frac{1}{h} \frac{dPn}{ds} = \frac{1}{h} \left[\Delta f(x0) + \frac{2s - 1}{2!} \Delta^2 f(x0) + \frac{3s^2 - 6s + 2}{3!} \Delta^3 f(x0)\right]

  • Here x = 2, x_0 = 1, s = 1, h = 1

  • f'(2) = 3 - \frac{1}{2} + \frac{2}{6} = \frac{13}{6}

  • f''(x) \approx \frac{1}{h^2} \frac{d^2Pn}{ds^2} = \frac{1}{h^2} \left[\Delta^2 f(x0) + (s - 1)\Delta^3 f(x_0)\right]

  • f''(2) = -1

  • Example 3: Calculate f^{(4)}(0.15)

    • Given data:
      • x: 0.1, 0.2, 0.3, 0.4, 0.5, 0.6
      • f(x): 0.425, 0.475, 0.400, 0.450, 0.525, 0.575

Numerical Integration

Basics

  • If f : [a, b] \rightarrow R is integrable, we obtain a new function F : [a, b] \rightarrow R defined by: F(x) = \int_{a}^{x} f(t)dt, for all x \in [a, b]
  • If f is a nonnegative function, then \int_{a}^{b} f(x)dx represents the area under the curve f(x).

Antiderivative

  • If f = F', then F is called an antiderivative of f.
  • Fundamental Theorem of Calculus: If f : [a, b] \rightarrow R is integrable and has an antiderivative F, then: \int_{a}^{b} f(x)dx = F(b) - F(a)

Basic Problems

  • Difficult to find an antiderivative (e.g., f(x) = e^{-x^2}).
  • Function is given in discrete form (data points).

Newton-Cotes Methods/Formulae

  • Based on polynomial interpolation.
  • Replace f by Pn(x) and evaluate \int{a}^{b} P_n(x)dx
  • \int{a}^{b} f(x)dx \approx \int{a}^{b} Pn(x)dx = \int{a}^{b} \sum{i=0}^{n} Li(x)f(xi)dx = \sum{i=0}^{n} f(xi) \int{a}^{b} Li(x)dx = \sum{i=0}^{n} Ai f(xi)
  • Where Ai = \int{a}^{b} L_i(x)dx is called the weight.

Types of Newton-Cotes Formulae

  • Quadrature formulae or rule
  • n=1: Trapezoidal Rule
  • n=2: Simpson’s rule
  • n=3: Simpson’s Three-Eighths rule
  • n=4:
  • Generalized formula

Quadrature Formulae

  • If x0, …, xn are distinct numbers in [a, b]
  • And f \in C^{n+1}[a, b]. Then, for all x \in [a, b], there exists \xi(x) \in (a, b) such that:
    f(x) = Pn(x) + \frac{f^{(n+1)}(\xi(x))}{(n + 1)!} (x - x0)(x - x1)…(x - xn)
  • Where Pn(x) is the interpolating polynomial defined as: Pn(x) = \sum{i=0}^{n} f(xi)L_i(x)
  • Quadrature Formula:
    f(x) \approx \sum{i=0}^{n} Ai f(x_i)
  • Error Term:
    E(f) = \frac{1}{(n + 1)!} \int{a}^{b} f^{(n+1)}(\xi(x)) \prod{i=0}^{n} (x - x_i) dx

Trapezoidal Rule

  • Objective: To approximate F = \int_{a}^{b} f(x)dx
  • First order Lagrange polynomial: n = 1, x0 = a, x1 = b, h = (b - a)
  • Pn(x) = \frac{x - x1}{x0 - x1} f(x0) + \frac{x - x0}{x1 - x0} f(x_1)
  • \int{a}^{b} f(x)dx = \int{x0}^{x1} \left[\frac{x - x1}{x0 - x1} f(x0) + \frac{x - x0}{x1 - x0} f(x1)\right]dx + \frac{1}{2} \int{x0}^{x1} f''(\xi(x))(x - x0)(x - x_1)dx

Trapezoidal Rule Error

  • The term (x - x0)(x - x1) does not change sign in [x0, x1] therefore applying weighted mean value theorem:
    \int{x0}^{x1} f''(\xi(x))(x - x0)(x - x1)dx= f''(\xi) \int{x0}^{x1} (x - x0)(x - x1)dx = f''(\xi) \Big[\frac{x^3}{3} - (x1 + x0) \frac{x^2}{2} + x0x1x \Big]{x0}^{x_1} = - \frac{h^3}{6} f''(\xi)

Trapezoidal Rule Formula

  • \int{a}^{b} f(x)dx=\frac{(x - x1)^2}{2(x0 - x1)}f(x0) + \frac{(x - x0)^2}{2(x1 - x0)}f(x1)\Big|{x0}^{x1} − \frac{h^3}{12}f''(\xi) = \frac{(x1 - x0)}{2} (f(x0) + f(x1)) - \frac{h^3}{12}f''(\xi)
  • With h = x1 - x0 the Trapezoidal rule is:
    \int{a}^{b} f(x)dx= \frac{h}{2} (f(x0) + f(x_1)) - \frac{h^3}{12}f''(\xi)
  • Exact for polynomials of degree one or less as the error term involves f''(\xi)

Simpson's 1/3 Rule

  • Objective: To approximate F = \int_{a}^{b} f(x)dx
  • Second order Lagrange polynomial: n = 2, x0 = a, x2 = b, x_1 = a + h, h = (b - a)/2
  • \int{a}^{b} f(x)dx = \int{x0}^{x2} Pn(x)dx + \int{x0}^{x2} (x - x0)(x - x1)(x - x2) \frac{f'''(\xi(x))}{6}dx
  • Third Taylor polynomial about x1: f(x) = f(x1) + f'(x1)(x - x1) + f''(x1)\frac{(x - x1)^2}{2} + f'''(x1)\frac{(x - x1)^3}{6} + \frac{f^{(4)}(\xi(x))}{24}(x - x_1)^4
  • Integration of the Taylor polynomial:
    \int{x0}^{x2} f(x)dx = 2hf(x1) + \frac{h^3}{3} f''(x1) + \frac{f^{(4)}(\xi1)}{60}h^5
  • Replacing f''(x1) by its approximation: \int{x0}^{x2} f(x)dx = \frac{h}{3} [f(x0) + 4f(x1) + f(x2)] - \frac{h^5}{12} (\frac{1}{3} f^{(4)}(\xi2) - \frac{1}{5} f^{(4)}(\xi_1) )

Simpson's 1/3 Rule Formula

  • \int{x0}^{x2} f(x)dx=\frac{h}{3} [f(x0) + 4f(x1) + f(x2)] - \frac{h^5}{90} f^{(4)}(\xi)

Simpson's 3/8 Rule

  • Objective: To approximate F = \int_{a}^{b} f(x)dx

  • Third order Lagrange polynomial: n = 3, x0 = a, x3 = a + 3h = b, x1 = a + h, x2 = a + 2h, h = (b - a)/3

  • \int{x0}^{x3} f(x)dx = \frac{3h}{8} [f(x0) + 3f(x1) + 3f(x2) + f(x_3)] - \frac{3h^5}{80} f^{(4)}(\xi), with \xi \in (a, b)

  • Example 4: Using Trapezoidal, Simpson 1/3 and 3/8 rules, find \int{0}^{2} x^4 dx and \int{0}^{2} sin(x) dx and find the upper bound for the error.

Composite Rules

  • If the interval [a, b] is large, the error using Trapezoidal or Simpson's rule will be large.
  • Idea: Divide [a, b] into equal subintervals and apply Trapezoidal or Simpson's rules in each subinterval to reduce error.

Composite Simpson's Rule

  • With h = (b - a)/n and x_j = a + jh, for j = 0, 1, …, n
  • \int{a}^{b} f(x)dx = \sum{j=1}^{n/2} \int{x{2j-2}}^{x{2j}} f(x)dx = \sum{j=1}^{n/2} \frac{h}{3} [f(x{2j-2}) + 4f(x{2j-1}) + f(x{2j})] - \frac{h^5}{90} f^{(4)}(\xij)
  • \int{a}^{b} f(x)dx= \frac{h}{3} f(x0) + 2 \sum{j=1}^{(n/2)-1} f(x{2j}) + 4\sum{j=1}^{n/2} f(x{2j-1}) + f(xn) − \frac{h^5}{90} \sum{j=1}^{n/2} f^{(4)}(\xi_j)

Composite Simpson's Rule Error

  • Error: E(f) = - \frac{h^5}{90} \sum{j=1}^{n/2} f^{(4)}(\xij) for \xij \in (x{2j-2}, x_{2j}) and j = 1, …, n/2
  • If f \in C^4[a, b], there exists \mu \in (a, b) such that:
    • E(f) = - \frac{(b - a)}{180} h^4 f^{(4)}(\mu)

Composite Simpson's Rule Theorem

  • Let f \in C^4[a, b] and n be even, h = (b - a)/n, xj = a + jh, for j = 0, 1, …, n. There exists a \mu \in (a, b) for which the Composite Simpson’s rule for n subintervals can be written with its error term as: \int{a}^{b} f(x)dx=\frac{h}{3} f(a) + 2 \sum{j=1}^{(n/2)-1} f(x{2j}) + 4\sum{j=1}^{n/2} f(x{2j-1}) + f(b) − \frac{(b - a)}{180} h^4 f^{(4)}(\mu)

Composite Trapezoidal Rule Theorem

  • Let f \in C^2[a, b] and h = (b - a)/n, xj = a + jh, for j = 0, 1, …, n. There exists a \mu \in (a, b) for which the Composite Trapezoidal rule for n subintervals can be written with its error term as: \int{a}^{b} f(x)dx=\frac{h}{2} f(a) + 2 \sum{j=1}^{(n-1)} f(x{j}) + f(b) − \frac{(b - a)}{12} h^2 f''(\mu)

Composite Midpoint Rule Theorem

  • Let f \in C^2[a, b] and n an even number, h = (b - a)/(n + 2), xj = a + (j + 1)h, for j = -1, 0, 1, …, n + 1. There exists a \mu \in (a, b) for which the Composite Midpoint rule for n + 2 subintervals can be written with its error term as: \int{a}^{b} f(x)dx= 2h\sum{j=0}^{n/2} f(x{2j}) + \frac{(b - a)}{6} h^2 f''(\mu)

  • Example 6: Determine values of n that will ensure an approximation error of less than 0.00002 when approximating \int_{0}^{\pi} sin(x)dx employing:

    • (a) Composite Trapezoidal rule
    • (b) Composite Simpson’s rule

  • For each case, approximate the integral with the smallest possible n.