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.