Differentiation

There are two situations where we can imagine needing to compute a derivative:

1. We have an analytic function, \(f(x)\), and we want to create a numerical approximation to its derivative

2. We have a function \(f(x)\) defined only at a finite set of (possibly regularly spaced) points, and we want to use that discrete data to estimate the derivative

For the first case, it is usually best to take the analytic derivative. In the previous notebook however, we did look at the effect of roundoff on computing a derivative.

We’ll focus on the second case here, but not that this can be applied to the case where the function is known analytically provided we evaluate it at a fixed set of points.

First order approximations

Consider a set of points labeled with an index \(i\), with the physical spacing between them denoted \(\Delta x\).

We’ll label the function value at \(x_i\) as \(f_i\), e.g., \(f_i = f(x_i)\).

We can use the result of the Taylor expansion we previously derived to write the derivative as:

\[\left . \frac{d f}{dx} \right |_i = \frac{f_{i+1} - f_i}{\Delta x} + \mathcal{O}(\Delta x)\]

where \(f_{i+1} = f(x_{i+1})\) is the data we have at the point \(i+1\).

As \(\Delta x \rightarrow 0\), this approaches the definition of the derivative from calculus. However, we are not free to choose \(\Delta x\)—it is a property of the discrete set of points we are given.

Note: we could alternately have used the point to the right of \(i\):

\[\left . \frac{d f}{dx} \right |_i = \frac{f_{i} - f_{i-1}}{\Delta x} + \mathcal{O}(\Delta x)\]

Second order approximation

Looking at the Taylor expansion of \(f_{i+1} = f(x_{i+1}) = f(x_i + \Delta x)\), we see

\[f_{i+1} = f_i + \Delta x \left .\frac{df}{dx} \right |_i + \frac{1}{2} \Delta x^2 \left . \frac{d^2f}{dx^2} \right |_i + \mathcal{O}(\Delta x^3)\]

likewise:

\[f_{i-1} = f_i - \Delta x \left .\frac{df}{dx} \right |_i + \frac{1}{2} \Delta x^2 \left . \frac{d^2f}{dx^2} \right |_i + \mathcal{O}(\Delta x^3)\]

subtracting these two expressions give:

\[f_{i+1} - f_{i-1} = 2 \Delta x \left . \frac{df}{dx} \right |_i + \mathcal{O}(\Delta x^3)\]

or

\[\left . \frac{df}{dx} \right |_i = \frac{f_{i+1} - f_{i-1}}{2 \Delta x} +\mathcal{O}(\Delta x^2)\]

This is called the centered-difference approximation to the first derivative. It is second order accurate.

Graphically, these different approximations appear as:

Example

Consider the function \(f(x) = \sin(x)\). The code below defines 10 equally spaced points and defines \(f(x)\) at each point.

Use this discrete data to estimate the derivative at x[3] and compute the error with respect to the analytic value.

import numpy as np

x = np.linspace(0, np.pi, 10, endpoint=False)
f = np.sin(x)

# first we'll write functions to evaluate each of the
# derivative approximations at a given index idx

def left_sided_deriv(x, f, idx):
    """return the left-sided derivative at x[idx]"""
    return (f[idx] - f[idx-1]) / (x[idx] - x[idx-1])

def right_sided_deriv(x, f, idx):
    """return the right-sided derivative at x[idx]"""
    return (f[idx+1] - f[idx]) / (x[idx+1] - x[idx])

def centered_deriv(x, f, idx):
    """return the left-sided derivative at x[idx]"""
    return (f[idx+1] - f[idx-1]) / (x[idx+1] - x[idx-1])

# always use x[ival] for the location of the derivative
ival = 3

dl = left_sided_deriv(x, f, ival)
dr = right_sided_deriv(x, f, ival)
dc = centered_deriv(x, f, ival)

print(f"left = {dl}, right = {dr}, centered = {dc}")
print(f"analytic value = {np.cos(x[ival])}")

left = 0.7042025064251414, right = 0.4521258405602084, centered = 0.5781641734926749
analytic value = 0.5877852522924731

Convergence

Now we can see if we get the convergence that truncation error suggests. We’ll look at the centered difference. As we cut \(h\) in half, the error should go down by \(1/4\).

def derivative(x0, f, h):
    """compute the derivative of the function f(x) at x0 using stepsize h"""
    return (f(x0 + h) - f(x0 - h)) / (2 * h)

def func(x):
    return np.sin(x)

def fprime(x):
    return np.cos(x)

h = 0.1
x0 = 1.0
while h > 1.e-5:
    err = np.abs(derivative(x0, func, h) - fprime(x0))
    print(f"{h:10.6g} : {err:20.10g}")
    
    h /= 2

       0.1 :      0.0009000536984
      0.05 :      0.0002250978217
     0.025 :      5.627973142e-05
    0.0125 :      1.407026262e-05
   0.00625 :      3.517586261e-06
  0.003125 :      8.793978465e-07
 0.0015625 :       2.19849555e-07
0.00078125 :      5.496242883e-08
0.000390625 :      1.374043412e-08
0.000195313 :      3.435006501e-09
9.76563e-05 :      8.588627587e-10
4.88281e-05 :       2.14826823e-10
2.44141e-05 :      5.339151343e-11
1.2207e-05 :       1.01905151e-11

We see the expected 2nd order convergence