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.