--- jupytext: formats: ipynb,md:myst text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.14.4 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- ```{code-cell} :tags: [hide-cell] import mmf_setup;mmf_setup.nbinit() import logging;logging.getLogger('matplotlib').setLevel(logging.CRITICAL) %matplotlib inline import numpy as np, matplotlib.pyplot as plt ``` (sec:Derivatives)= Derivatives =========== ```{code-cell} :tags: [margin, hide-input] 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)); ``` ## 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: :::{margin} The truncation error terms here are [Lagrange form of the remainder][] and are exact for some $x^* \in [x, x+h]$ and valid if $f(x)$ is appropriately differentiable. [Lagrange form of the remainder]: **To Do:** make this precise. Here $x^*$ is a point somewhere in the interval $[x-h, x+h]$, as follows from the mean value theorem, $\bar{f}$ is a typical value of $f$ at the evaluated points, and $\epsilon$ is the machine precision. The round-off error is estimated by assuming that the relative error of all computations is accurate to machine precision $\epsilon$. The code below shows that taking the maximum absolute values of the various derivatives works nicely. ::: \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*} :::{admonition} 5-point stencil for $f''(x)$. :class: dropdown 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*} ::: ```{code-cell} :tags: [margin] 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() ``` 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*} :::::{doit} 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 {cite}`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$. ::: ::::: ```{code-cell} :tags: [hide-input] 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]); ``` These techniques can be implemented fairly efficiently for images using {func}`scipy.ndimage.correlate1d`: ```{code-cell} 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) ``` (sec:Drums)= ## 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: ```{code-cell} :tags: [hide-input] %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. ``` 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 {ref}`sec:FourierTechniques`. [Physics 571: Drums]: