A silica bead of radius is held in water ( ) at by an optical trap of stiffness

PHYSICS · Multi-part · Fermi University Physics · Mathematical Methods for Physics · Probability, statistics and stochastic methods · Numerical and computational mathematics

Error, Chance and the Machine, Ch Partial differential equations: stability, convergence and solvers

Problem

A silica bead of radius a=0.50 μma = 0.50\,\mu\mathrm{m} is held in water (η=1.002×10−3 Pa s\eta = 1.002\times10^{-3}\,\mathrm{Pa\,s}) at T=293 KT = 293\,\mathrm{K} by an optical trap of stiffness ktrap=1.0 pN μm−1k_{\rm trap} = 1.0\,\mathrm{pN\,\mu m^{-1}}, which pulls it back along one axis with the force −ktrapx-k_{\rm trap}x. Use kB=1.3806×10−23 J K−1k_{\rm B} = 1.3806\times10^{-23}\,\mathrm{J\,K^{-1}}. For times long against the velocity relaxation time (Chapter 9) the position obeys dx=−xτtrap dt+2D dW,τtrap=γktrap,\mathrm{d}x = -\frac{x}{\tau_{\rm trap}}\,\mathrm{d}t + \sqrt{2D}\,\mathrm{d}W,\qquad \tau_{\rm trap} = \frac{\gamma}{k_{\rm trap}}, with γ=6πηa\gamma = 6\pi\eta a and D=kBT/γD = k_{\rm B}T/\gamma, and the density p(x,t)p(x,t) of many beads obeys ∂p∂t=D ∂2p∂x2+1τtrap ∂∂x(xp).\frac{\partial p}{\partial t} = D\,\frac{\partial^2p}{\partial x^2} + \frac{1}{\tau_{\rm trap}}\,\frac{\partial}{\partial x}\bigl(xp\bigr). The stationary density is the Boltzmann density ∝e−ktrapx2/2kBT\propto e^{-k_{\rm trap}x^2/2k_{\rm B}T}. We find its variance by four routes and compare them. (Here ktrapk_{\rm trap} is a stiffness and not a wavenumber.) ((a)) Find γ\gamma, DD, τtrap\tau_{\rm trap} and the standard deviation σ\sigma of the Boltzmann density. ((b)) Simulate xx by the Euler-Maruyama method of Chapter 12 with step h=1.0 msh = 1.0\,\mathrm{ms}. Write the step as xn+1=ϕxn+2Dh Znx_{n+1} = \phi x_n + \sqrt{2Dh}\,Z_n, find ϕ\phi, and predict the stationary variance of the simulated positions as a multiple of σ2\sigma^2. ((c)) Run 4×1054\times10^5 steps from x0=0x_0 = 0, after discarding the first 200200. (The solution quotes NumPy's PCG64 with seed 1.) The record is a sequence of the kind solved in Chapter 8, so give the exact τint\tau_{\rm int} of xnx_n in steps. Estimate τint\tau_{\rm int} of xn2x_n^2 from the run by the window rule of Chapter 11, and quote x2‾/σ2\overline{x^2}/\sigma^2 with its standard error. How many standard errors is it from the prediction in (b), and from 11? What would the standard error be if the correlation between steps were ignored? ((d)) Solve the Fokker-Planck equation on −320 nm≤x≤320 nm-320\,\mathrm{nm} \le x \le 320\,\mathrm{nm} with Δx=5 nm\Delta x = 5\,\mathrm{nm} and p=0p = 0 at the two ends, so that there are 127127 unknowns. At node jj, where xj=jΔxx_j = j\Delta x, difference the drift term centrally as (xj+1pj+1−xj−1pj−1)/(2Δx τtrap)(x_{j+1}p_{j+1} - x_{j-1}p_{j-1})/(2\Delta x\,\tau_{\rm trap}). (a) Find the explicit step limit of the diffusion term and the number of such steps in τtrap\tau_{\rm trap}. Find also the time Δx/∣A∣\Delta x/|A| in which the drift speed ∣A∣=∣x∣/τtrap|A| = |x|/\tau_{\rm trap} at x=320 nmx = 320\,\mathrm{nm} carries density across one cell, and say which of the two limits is the tighter. (b) Use Crank-Nicolson with Δt=0.5 ms\Delta t = 0.5\,\mathrm{ms}. Find α=DΔt/Δx2\alpha = D\Delta t/\Delta x^2. Write the three coefficients of row jj of the tridiagonal system for pn+1\mathbf{p}^{n+1}, find the sum of the sizes of the two off-diagonal ones, and show that the matrix is not symmetric but is diagonally dominant, so that the solver of Chapter 2 needs no pivoting. (c) Central differencing of a drift term gives no spurious oscillation when ∣A∣Δx/D<2|A|\Delta x/D < 2 (quoted). Check this at x=320 nmx = 320\,\mathrm{nm}. (d) Start from a Gaussian of mean 120 nm120\,\mathrm{nm} and standard deviation 25 nm25\,\mathrm{nm}, not from a single occupied cell, and explain why. Run 1919 steps (9.5 ms9.5\,\mathrm{ms}) and 190190 steps (95 ms95\,\mathrm{ms}), normalising pp so that ∫p dx=1\int p\,\mathrm{d}x = 1 before taking moments. Give the variance as a multiple of σ2\sigma^2 at both times, and compare the first with the exact value from the moment equations of Chapter 9. ((e)) Sample the Boltzmann density by random-walk Metropolis (Chapter 11), with proposal y=x+σξy = x + \sigma\xi, ξ\xi a standard normal number. Predict the acceptance rate. Run 2×1052\times10^5 steps from x=0x = 0 after discarding 10001000 (PCG64, seed 1), estimate τint\tau_{\rm int} of x2x^2 by the same window rule, and quote x2‾/σ2\overline{x^2}/\sigma^2 with its standard error. ((f)) Compare the four variances as multiples of σ2\sigma^2. For each of the Euler-Maruyama, Crank-Nicolson and Metropolis results, say what limits its accuracy (the time step, the grid or the run length) and whether halving hh, halving Δt\Delta t or doubling the run would change it.

Browse the problem library