Derivatives

Hide code cell content

 
import mmf_setup;mmf_setup.nbinit()
import logging;logging.getLogger('matplotlib').setLevel(logging.CRITICAL)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt

This cell adds /home/docs/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:

  • Choose "Trust Notebook" from the "File" menu.
  • Re-execute this cell.
  • Reload the notebook.

Derivatives#

Hide code cell source

 
from phys_581.testing import Functions
L = 10.0
x = np.linspace(0, L)
f = Functions(L=L, eta=0.5)
fig, ax = plt.subplots(figsize=(4,3))
ax.plot(x, f(x))
ax.plot(x, f(x, d=1))
ax.plot(x, f(x, d=2));
../_images/8ac2bf19db1934a9c85a6caa66e77c676ec1b46fef7b5bfd7eee9bc6b0ea4910.png

Finite Differences#

A straightforward way of computing derivatives is to use finite difference formula, derive from the Taylor series:

\[\begin{gather*} f(x \pm h) = f(x) \pm h f'(x) + \frac{h^2}{2} f''(x) \pm \frac{h^3}{3!}f'''(x) + \cdots. \end{gather*}\]

The simplest of these are:

\[\begin{alignat*}{3} f_{R}'(x) &= \frac{f(x+h) - f(x)}{h} & &+ \frac{h}{2}f''(x^{*}) & &+ \epsilon\frac{2\bar{f}}{h},\\ f_{C}'(x) &= \frac{f(x+h) - f(x-h)}{2h} & &+ \frac{h^2}{6}f'''(x^{*}) & &+ \epsilon\frac{2\bar{f}}{2h},\\ f_{C}''(x) &= \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} & &+ \frac{h^2}{12}f''''(x^{*}) & &+ \epsilon\frac{4\bar{f}}{h^2}. \end{alignat*}\]

5-point stencil for \(f''(x)\).

To improve precision, a 5-point stencil [-1, 16, -30, 16, -1]/12 for computing \(f''(x)\) is useful. To simplify notation, we write \(f_{n} = f(x+nh)\):

\[\begin{gather*} f''_{5}(x) = \frac{-(f_{2}+f_{-2}) + 16(f_1 + f_{-1}) - 30 f_{0}}{12h^2} + \frac{h^4}{90}f^{(6)}(x^*) + \epsilon\frac{4\bar{f}}{h^2},\\ h_{5} = \sqrt[6]{180\epsilon \frac{\bar{f}}{\bar{f}^{(6)}}}. \end{gather*}\]
import sympy
x, h = sympy.var('x, h')
f = sympy.Function('f')
f2 = (-(f(x+2*h) + f(x-2*h)) + 16 * (f(x+h) + f(x-h)) - 30*f(x))/12/h**2
sympy.series(f2, h, 0, 6).simplify()
\[\displaystyle \frac{d^{2}}{d x^{2}} f{\left(x \right)} - \frac{h^{4} \frac{d^{6}}{d x^{6}} f{\left(x \right)}}{90} + O\left(h^{6}\right)\]

The first error term is due to truncation, while the second estimates the round-off error. Minimizing the error gives an estimate of the optimal step size:

\[\begin{gather*} h_{1R} \approx \sqrt{4\epsilon\frac{\bar{f}}{\bar{f}''}},\qquad h_{1C} \approx \sqrt[3]{3\epsilon\frac{\bar{f}}{\bar{f}'''}},\qquad h_{2C} \approx \sqrt[4]{48\epsilon\frac{\bar{f}}{\bar{f}''''}}. \end{gather*}\]

Do It! Prove the Lagrange remainder theorem.

Prove that under suitable conditions on \(f(x)\) (state these):

\[\begin{gather*} f'_R(x) = \frac{f(x + h) - f(x)}{h} + \frac{h}{2}f''(x^*), \qquad x^* \in [x, x+h], \end{gather*}\]

and similarly for the other derivatives (assuming no round-off error). Hint: use the mean-value theorem.

Solution

A simple proof of Taylor’s theorem is due to [Cox, 1851] who considers the following quantity

\[\begin{align*} Q_{n}(x) = & f(a+x) - f(a) - xf'(a) - \frac{x^2}{2!}f''(a) - \dots \frac{x^{n-1}}{(n-1)!}f^{(n-1)}(a)\\ -\frac{x^n}{h^n}\Bigl( &f(a+h) - f(a) - hf'(a) - \frac{h^2}{2!}f''(a) - \dots \frac{h^{n-1}}{(n-1)!}f^{(n-1)}(a)\Bigr). \end{align*}\]

Cox argues that, since \(Q_n(0) = Q_n(h) = 0\), by the mean-value theorem, \(Q'_n(x_1) = 0\) for some intermediary point \(x_1 \in (0, h)\). Continuing, since \(Q'_{n}(0) = 0\), there is another intermediary point \(x_2 \in (0, x_1)\) where \(Q''_{n}(x_2) = 0\) vanishes, and so forth until finally \(Q^{(n)}_{n}(x_{n}) = 0\) for some \(x_n \in (0, h)\). This implies

\[\begin{multline*} 0 = Q_{n}^{(n)}(x_{n}) = f^{(n)}(a+x_n) \\ - \frac{n!}{h^{n}}\Bigl( f(a+h) - f(a) - hf'(a) - \frac{h^2}{2!}f''(a) - \dots \frac{h^{n-1}}{(n-1)!}f^{(n-1)}(a)\Bigr), \end{multline*}\]

establishing the Lagrange form of the remainder in Taylor’s theorem

\[\begin{gather*} f(a+h) = f(a) + hf'(a) + \frac{h^2}{2!}f''(a) + \dots \frac{h^{n-1}}{(n-1)!}f^{(n-1)}(a) + \frac{h^n}{n!}f^{(n)}(x^*) \end{gather*}\]

for some \(x^* = a+x_n \in (a, a+h)\). For this argument to hold, all \(n\) derivatives of \(f(x)\) must be continuous. It remains a simple exercise to verify that \(Q_{n}(0) = Q_{n}(h) = 0\) and that all derivatives \(Q_{n}^{(i)}(0) = 0\) for \(i < n\).

Hide code cell source

 
from phys_581.testing import Functions
L = 10.0
x = np.linspace(0, L)
f = Functions(L=L, eta=0.5)
fig, ax = plt.subplots()

hs = 10**np.linspace(-9, 1)
errs = np.array([(abs((f(x+h) - f(x))/h - f(x, d=1)).max(), 
         abs((f(x+h) - f(x-h))/2/h - f(x, d=1)).max(), 
         abs((f(x+h) + f(x-h) - 2*f(x)) / h**2 - f(x, d=2)).max())
         for h in hs]).T

eps = np.finfo(float).eps
f1 = abs(f(x)).mean()
f2 = abs(f(x, d=2)).max()
h = (3*eps)**(1/3)
f3 = abs(f(x+h, d=2) - f(x-h, d=2)).max() / 2 / h
h = (12*eps)**(1/4)
f4 = abs(f(x+h, d=2) + f(x-h, d=2) - 2*f(x, d=2)).max() / h**2

h_opts = [
    (4*eps * f1/f2)**(1/2),
    (3*eps * f1/f3)**(1/3),
    (48*eps * f1/f4)**(1/4)
]


truncation_errs = [
    hs/2*f2,
    hs**2/6*f3,
    hs**2/12*f4,
]

roundoff_errs = [
    eps * 2*f1/hs,
    eps * 2*f1/2/hs,
    eps * 4*f1/hs**2,
]

labels = [
    "Forward difference",
    "Centered difference",
    "Second difference"
]

for n, (err, h_opt) in enumerate(zip(errs, h_opts)):
    c = f"C{n}"
    err_t = truncation_errs[n]
    err_r = roundoff_errs[n]
    ax.loglog(hs, err, "+-", c=c, label=labels[n])
    ax.loglog(hs, err_t, "--", c=c)
    ax.loglog(hs, err_r, "-.", c=c)
    ax.axvline(h_opt, ls=":", c=c)
ax.legend()
ax.set(xlabel='$h$', ylabel='err', ylim=[1e-12, 1]);
../_images/48390b20d3530375400f9d883d144323543c025da2c633402835806130ba56be.png

These techniques can be implemented fairly efficiently for images using scipy.ndimage.correlate1d():

from scipy.ndimage import correlate1d
from phys_581.testing import Functions

L = 10.0
Nx = 100
h = dx = L/Nx

# Note: it is important to use proper lattice points - excluding one endpoint - if
# wrapping.

x = np.arange(Nx) * dx
f = Functions(L=L, eta=0.5)

df = correlate1d(f(x), weights=[-1/2/dx, 0, 1/2/dx], mode='wrap')

assert np.allclose(
    correlate1d(f(x), weights=[-1/2/dx, 0, 1/2/dx], mode='wrap'),
    (f(x+h) - f(x-h))/2/h)

assert np.allclose(
    correlate1d(f(x), weights=[0, -1/dx, 1/dx], mode='wrap'),
    (f(x+h) - f(x))/h)

assert np.allclose(
    correlate1d(f(x), weights=[1/dx**2, -2/dx**2, 1/dx**2], mode='wrap'),
    (f(x+h) + f(x-h) - 2*f(x))/h**2)

Application: Drums#

A neat application of finite differences is to estimate the normal modes of a drum.

Consider a thin 2D drum-head streched taught over a frame outlining some specified region of the \(x\)-\(y\) plane. Small-amplitude oscillations of the drum can be described by their height \(f(x,y,t)\) and will satisfy the 2D wave equation:

\[\begin{gather*} \pdiff[2]{}{t}f(x,y,t) = c^2\nabla^2 f(x,y,t) \end{gather*}\]

where \(c\) is the speed of sound in the material which will depend on the tension and mass density of the drumhead (assumed to be uniform). The solutions must satisfy the boundary condition \(f(x,y,t) = 0\) wherever the point \((x,y)\) lies on (or outside) of this region. (Dirichlet boundary conditions.)

The strategy will be to formulate a matrix representation of the laplacian \(\nabla^2\) and then look for eigenfunctions \(\nabla^2 f_n = -\frac{\omega^2_n}{c^2} f_n\). From these eigenfunctions we can form the time-dependent solutions:

\[\begin{gather*} f(x,y,t) = \Re\left(\sum_{n} a_n e^{\I\omega_n t} f_n(x,y)\right) \end{gather*}\]

where \(a_n\) are arbitrary complex coefficients.

The trick is to note that we can trivially implement Dirichlet boundary conditions by simply discarding points outside of the region that vibrates. Thus, we first form a finite difference approximation to \(\nabla^2\) on a regular grid using a tensor product of two 1-dimensional finite difference matrices, then we simply select the rows and columns that correspond to vibrating points. The resulting matrix is Hermitian, and we simply find the eigenvalues and eigenvectors which give us the normal modes and frequencies:

Hide code cell source

 
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
from scipy import special

L = 2.2  # Size of box
R = 1.0  # Radius of drum head
N = 32  # Number of points
h = L / N  # Lattice spacing

# The abscissa are adjusted so that zero is in the center
x = y = np.arange(N) * h - L / 2.0 + h / 2.0
X = x[:, None]
Y = y[None, :]

# Here is the matrix D2 expressing Dirichlet boundary conditions
# Note that np.eye() returns an "identity" matrix with ones on
# the k'th diagonal
d2 = (np.eye(N, k=1) + np.eye(N, k=-1) - 2 * np.eye(N)) / h**2

# Here we form the tensor product using Einstein notation
D2x = np.einsum('ab,ij->aibj', d2, np.eye(N))
D2y = np.einsum('ab,ij->iajb', d2, np.eye(N))
D2 = (D2x + D2y).reshape((N**2, ) * 2)

# Get the indices into the 1D ravelled array where the membrane can fluctuate:
inds = np.where((X**2 + Y**2 < R**2).ravel())[0]

# Now restrict the Laplacian to these indices:
D2_ = D2[inds, :][:, inds]

# Make sure it is Hermitian
assert np.allclose(D2_, D2_.T)

# Find modes by diagonalizing the negative of the matrix
# so that the modes are sorted by energy
Es, Vs = np.linalg.eigh(-D2_)

# Plot several modes
modes = 10
#fig, axs = plt.subplots(2, modes, figsize=(modes, 4))  # Make figure wide enough.
fig, axs = plt.subplots(1, modes, figsize=(modes, 4))  # Make figure wide enough.
for n in range(modes):
    ax = axs[n]
    f = np.zeros((N, N))     # Full grid of zeros
    f.flat[inds] = Vs[:, n]  # Insert active points into f
    f_ = abs(f).max()        # Normalize colors so white is zero
    
    ax.contourf(
        x,
        y,           # Filled contour plot.
        f.T,         # Transpose because this follows MATLAB's conventions
        15,          # Number of contours
        cmap='bwr',  # Divergent colormap to show both + and - clearly
        vmin=-f_,    # These need to be set so that 0 is in the middle of the map 
        vmax=f_)     #   and appears white in the plots
    ax.set(aspect=1, # Make circles circular
           title=f"{Es[n]/Es[0]:.2f}")  # Energy as title
    ax.set_axis_off()  # Turn off ticks, numbers, etc.
../_images/d2d89b3863d98c6c94cea561c3453dabbb2478cf578bfb8b710448c56e35a69e.png

For more examples see Physics 571: Drums.

Warning

This trick of enforcing the boundary conditions by simply discarding lattice points outside of the zone of vibration is specific to Dirichlet boundary conditions: it depends crucially on the property that truncation is equivalent to setting the function to zero in the discarded regimes.

We will continue our exploration computing derivatives with Fourier Techniques.

Contents