Calculus: Differentiation

Numerical Methods

David Mayerich

Scalable Tissue Imaging and Modeling (STIM) Laboratory

Department of Electrical and Computer Engineering

Cullen College of Engineering

University of Houston

David Mayerich

STIM Laboratory, University of Houston

Finite Difference Methods

Discrete Functions

Taylor Series Approximations to Derivatives

Forward, Backward, and Central Differences

David Mayerich

STIM Laboratory, University of Houston

Specifying Functions

Formulaic or algorithmic:

David Mayerich

STIM Laboratory, University of Houston

f(x, y, n) = \cos(e^x) + \sin(y)j_n(x)

double f(double x, double y, int n) {
	return cos(exp(x)) + sin(y) * jn(n, x);
}

x, y, f \in \mathbb{R} \quad n \in \mathbb{N}

Look-up tables:

X \in \mathbb{R}^n = \{x_1, x_2, x_3, \cdots, x_n\}

F \in \mathbb{R}^n = \{f_1, f_2, f_3, \cdots, f_n\}

f(x_k)=F[k] = f_k

It is common to use a uniform spacing between samples: \(x_{n+1} = x_n + x_\Delta\)

double* cos_lut(int n) {
	double* lut = new double[n];    
    double dx = (2.0 * M_PI) / (double)n;
    for(int i = 0; i < n; i++) {
        double x = i * dx;
        lut[i] = cos(x);
    }
    return lut;
}

Specifying Functions

Formulaic or algorithmic:

David Mayerich

STIM Laboratory, University of Houston

f(x, y, n) = \cos(e^x) + \sin(y)j_n(x)

double f(double x, double y, int n) {
	return cos(exp(x)) + sin(y) * jn(n, x);
}

x, y, f \in \mathbb{R} \quad n \in \mathbb{N}

Look-up tables:

X \in \mathbb{R}^n = \{x_1, x_2, x_3, \cdots, x_n\}

F \in \mathbb{R}^n = \{f_1, f_2, f_3, \cdots, f_n\}

f(x_k)=F[k] = f_k

It is common to use a uniform spacing between samples: \(x_{n+1} = x_n + x_\Delta\)

double* cos_lut(int n) {
	double* lut = new double[n];    
    double dx = (2.0 * M_PI) / (double)n;
    for(int i = 0; i < n; i++) {
        double x = i * dx;
        lut[i] = cos(x);
    }
    return lut;
}

int N = 100;
double theta = 1.23;
double* cos = cos_lut(N);
int xi = theta / (2 * M_PI) * N;
double cos_x = cos[xi];

Numerical Differentiation

May not have access to the algorithm used to generate \(f(x)\)
Some algorithms and optimization (ex. Newton's method) rely on \(f'(x)\)
Can we approximate derivatives from samples or implementations of a function?

David Mayerich

STIM Laboratory, University of Houston

Approximations from Taylor Expansion

Start from a first-order Taylor series:

David Mayerich

STIM Laboratory, University of Houston

f(x)\approx f(c)+f'(c)(x-c) + O\left[(x-c)^2\right]

Given two known points \(c=x_i\) and \(x=x_{i+1}=x_i + x_\Delta\):

f(x_{i+1}) \approx f(x_i)+f'(x_i)(x_{i+1}-x_i) + O\left[(x_{i+1}-x_i)^2\right]

f(x_{i+1}) \approx f(x_i)+f'(x_i)(x_\Delta) + O\left(x_\Delta^2\right)

Solve for \(f'(x)\):

f'(x_i) \approx \frac{f(x_{i+1})-f(x_i) + O\left(x_\Delta^2 \right)}{x_\Delta}

Note that the order in the error term is reduced:

f'(x_i) \approx \frac{f(x_{i+1})-f(x_i)}{x_\Delta} + O\left(x_\Delta\right)

Forward Finite Difference

David Mayerich

STIM Laboratory, University of Houston

\frac{d}{dt}f(x_2)\approx \frac{f(x_3)-f(x_2)}{x_3 - x_2}

f'(x_i) \approx \frac{f(x_{i+1})-f(x_i)}{x_\Delta} + O\left(x_\Delta\right)

f(x_0)

f(x_1)

f(x_3)

f(x_4)

f(x_2)

Backwards Differences

Consider evaluating the function at a new position \(f(x_{i-1})\), where \(x_{i-1}=x_i - x_\Delta\):

David Mayerich

STIM Laboratory, University of Houston

f(x)\approx f(c)+f'(c)(x-c) + O\left[(x-c)^2\right]

f(x_{i-1}) \approx f(x_i)+f'(x_i)(x_{i-1}-x_i) + O\left[(x_{i-1}-x_i)^2\right]

f(x_{i+1}) \approx f(x_i)-f'(x_i)(x_\Delta) - O\left(x_\Delta^2\right)

f(x_{i+1}) \approx f(x_i)-f'(x_i)(x_\Delta) + O\left(x_\Delta^2\right)

Solve for \(f'(x)\):

f'(x_i) \approx \frac{f(x_{i})-f(x_{i-1}) + O\left(x_\Delta^2 \right)}{x_\Delta}

f'(x_i) \approx \frac{f(x_{i})-f(x_{i-1})}{x_\Delta} + O\left(x_\Delta\right)

Forward and Backward Differences

David Mayerich

STIM Laboratory, University of Houston

\frac{d}{dt}f(x_2)\approx \frac{f(x_3)-f(x_2)}{x_3 - x_2}

f'(x_i) \approx \frac{f(x_{i+1})-f(x_i)}{x_\Delta} + O\left(x_\Delta\right)

f'(x_i) \approx \frac{f(x_{i})-f(x_{i-1})}{x_\Delta} + O\left(x_\Delta\right)

\frac{d}{dt}f(x_3)\approx \frac{f(x_3)-f(x_2)}{x_3 - x_2}

f(x_0)

f(x_1)

f(x_3)

f(x_4)

f(x_2)

Second Order Approximations

Start from a 2nd-order Taylor expansion:

David Mayerich

STIM Laboratory, University of Houston

f(x)\approx f(c)+f'(c)(x-c) + \frac{f''(c)}{2}(x-c)^2 + O\left[(x-c)^3\right]

Substitute \(x=x_{i+1}\), \(x_i=c\), and \(x_\Delta = x_{i+1} - x_i\):

f(x_{i+1})\approx f(x_i)+f'(x_i)x_\Delta + \frac{f''(x_i)}{2}x_\Delta^2 + O\left[x_\Delta^3\right]

Substitute \(x=x_{i-1}\), \(x_i=c\), and \(x_\Delta = x_{i-1} - x_i\):

f(x_{i-1})\approx f(x_i)-f'(x_i)x_\Delta + \frac{f''(x_i)}{2}x_\Delta^2 + O\left[x_\Delta^3\right]

This formulation requires the higher-order derivative \(f''\)

Central Difference Approximation

Given the second-order approximation for \(f(x_{i-1})\) and \(f(x_{i+1})\) from the previous slide:

David Mayerich

STIM Laboratory, University of Houston

f(x_{i+1})\approx f(x_i)+f'(x_i)x_\Delta + \frac{f''(x_i)}{2}x_\Delta^2 + O\left[x_\Delta^3\right]

f(x_{i-1})\approx f(x_i)-f'(x_i)x_\Delta + \frac{f''(x_i)}{2}x_\Delta^2 + O\left[x_\Delta^3\right]

Remove \(f''\) by calculating \(f(x_{i+1})-f(x_{i-1})\):

\begin{split} f(x_{i+1}) - f(x_{i-1}) &\approx f(x_i)+f'(x_i)x_\Delta + \frac{f''(x_i)}{2}x_\Delta^2 + O\left[x_\Delta^3\right]\\ &- \left( f(x_i)-f'(x_i)x_\Delta + \frac{f''(x_i)}{2}x_\Delta^2 + O\left[x_\Delta^3\right] \right) \end{split}

f(x_{i+1}) - f(x_{i-1}) \approx 2f'(x_i)x_\Delta + O\left[x_\Delta^3\right]

f'(x_i) \approx \frac{f(x_{i+1}) - f(x_{i-1})}{2 x_\Delta} + O\left[x_\Delta^2\right]

solve for \(f'\)

Finite Differences

David Mayerich

STIM Laboratory, University of Houston

\frac{d}{dt}f(x_2)\approx \frac{f(x_3)-f(x_2)}{x_3 - x_2}

\frac{d}{dt}f(x_3)\approx \frac{f(x_3)-f(x_2)}{x_3 - x_2}

f'(x_i) \approx \frac{f(x_{i+1}) - f(x_{i-1})}{2 x_\Delta}

f(x_0)

f(x_1)

f(x_2)

f(x_3)

f(x_4)

Standard Finite Difference Methods

First forward finite difference approximation:

David Mayerich

STIM Laboratory, University of Houston

f'(x) \approx \frac{f(x + x_\Delta)-f(x)}{x_\Delta} + O\left(x_\Delta\right)

First backward finite difference approximation:

f'(x) \approx \frac{f(x)-f(x - x_\Delta)}{x_\Delta} + O\left(x_\Delta\right)

First central finite difference approximation:

f'(x) \approx \frac{f(x + x_\Delta) - f(x - x_\Delta)}{2 x_\Delta}

at \(O(x_\Delta)\) error

at \(O(x_\Delta^2)\) error

Mean Value Theorem

The mean value theorem can be derived from the central difference approximation:

David Mayerich

STIM Laboratory, University of Houston

f'(x) \approx \frac{f(x + x_\Delta) - f(x - x_\Delta)}{2 x_\Delta}

Assume we want to evaluate \(f(b)-f(a)\) where \(b > a\):

b = x_{i+1}

a = x_{i-1}

f'(x_i) \approx \frac{f(b) - f(a)}{b-a}

b - a = 2 x_\Delta

f(b) - f(a) \approx f'(\theta)(b-a)

where \(\theta = x_i = \frac{b - a}{2}\)

Examples

Compute the derivative of \(\log{x}\) using forward and central differences with 3 significant digits, \(x=100\), and \(x_\Delta=10\). How many bits are lost? What is the relative error?

David Mayerich

STIM Laboratory, University of Houston

\frac{\log{110} - \log{100}}{10} \approx \frac{2.04 - 2.00}{10}

=\frac{0.04}{10} = 0.004

\frac{- \log{\left(1 - \frac{\log{100}}{\log{110}}\right)}}{\log{2}} \approx 5.62

Forward

Central

\frac{\log{110} - \log{90}}{20} \approx \frac{2.04 - 1.95}{20}

=\frac{0.09}{20} = 0.0045

\frac{- \log{\left(1 - \frac{\log{90}}{\log{110}}\right)}}{\log{2}} \approx 4.55

\epsilon_r = \frac{|0.00434 - 0.004|}{|0.00434|} = 7.8 \%

\epsilon_r = \frac{|0.00434 - 0.0045|}{|0.00434|} = 3.7 \%

Examples

Compute the derivative using central differences and 3 significant figures

David Mayerich

STIM Laboratory, University of Houston

\(f(x)=e^x\) for \(x=1\) and \(x_\Delta=0.5\)

\frac{e^{1.5} - e^{0.5}}{2(0.5)} \approx \frac{4.48-1.65}{1.00} = 2.83

What is the relative error?

\epsilon_r = \frac{|2.83 - 2.73|}{|2.72|} \approx 4\%

f(x)=2.72

How many bits are lost if \(x_\Delta = 0.01\)

\frac{-\log{\left(1 - \frac{e^{0.99}}{e^{1.01}} \right)}}{\log{2}}

\approx \frac{-\log{\left(1 - \frac{2.69}{2.74} \right)}}{\log{2}}

\approx \frac{-\log{\left(1 - 0.98 \right)}}{\log{2}}

\approx 5.66

between \(5\) and \(6\) bits are lost

Examples

Calculate \(f'(x)\) where \(f(x)=|x|\) at \(x=0.5\) and \(x_\Delta = 0.75\)

David Mayerich

STIM Laboratory, University of Houston

f'(x)=1.00

f_f'(x)=\frac{1.25-0.50}{0.75} = 1.00

f_b'(x)=\frac{0.50-0.25}{0.75} = 0.333

f_c'(x)=\frac{1.25-0.25}{1.50} = 0.667

\epsilon_r = 0\%

\epsilon_r = \frac{|1.00-0.333|}{|1.00|} = 66.7\%

\epsilon_r = \frac{|1.00-0.667|}{|1.00|} = 33.3\%

x - x_\Delta

x + x_\Delta