Lecture 6 Notes: Solving the Diffusion Equation

18.354J Nonlinear Dynamics II: Continuum Systems 6 Lecture 6 Solving the diffusion equation We have shown, through two different arguments, that the density of random walkers on a one dimensional lattice obeys the diffusion equation, ∂n ∂2n = D 2. ∂t ∂x (103) This description is valid whenever examining the dynamics of large quantities of random walkers on scales much larger than the lattice spacing. As the next step, it is important to understand how to solve this equation, as the same mathematical problem will arise later on in our studies of fluid motion. For instance we would like to know the solution to the above equation, subject to the initial condition n(x, t = 0) = n0 (x) and the boundary conditions that n vanishes at ±∞. There are two basic techniques for solving this, each of which relies on a different method for representing the solution. In both cases, the central idea is that since the equation is linear, it is possible to “break down” any initial state into a linear combination of simpler problems. By solving the simpler problems explicitly it is then possible to reconstruct the general solution. 6.1 Fourier method This method relies on the fact that it is possible to express n(x, t) in a basis of plane waves3 , i.e. Z ∞ 1 n(x, t) = eikx n̂(k, t)dk. (104) 2π −∞ As a complement to (104), the Fourier coefficients for a given distribution are found using the Fourier transform Z ∞ n̂(k, t) = e−ikx n(x, t)dx, (105) −∞ and we define n̂0 (k) to be the Fourier coefficients of the initial condition n0 (x). The strength of this approach is that it is simple to solve the diffusion equation for a single plane wave. For example, integrating the diffusion equation in the following manner Z ∞ Z ∞ ∂n −ikx ∂2n e dx = D 2 e−ikx dx (106) ∂x −∞ ∂t −∞ gives ∂n̂(k, t) = (ik)2 Dn̂(k, t) = −k 2 Dn̂(k, t). ∂t Given the initial condition n0 (x) → n̂0 (k) the solution of (107) is 2 n̂(k, t) = n̂0 (k)e−Dk t . (107) (108) 3 An intuitive way of thinking is to note that a plane wave can be written as eikx = cos(kx) + i sin(kx). The Fourier transform thus can be interpreted as expressing a (complex) function as a superposition of real and complex sinusoidal waves. The Fourier coefficient n̂(k, t) describes how large the contribution of a wave with wave vector k ∝ 1/λ, with λ the wavelength, is in this superposition. The solution of the original problem is therefore Z ∞ 1 2 n(x, t) = n̂0 (k)eikx−Dk t dk. 2π −∞ (109) We note that the high-wavenumber components (which correspond to sharp gradients) are rapidly damped, emphasising the smoothing property of diffusion. Example Consider a distribution that is initially Gaussian (normal) about the point x = 0 at time t = 0, with standard deviation σ, n(x, 0) = √ x2 1 2πσ 2 e− 2σ2 . (110) Note that in the limit where σ → 0, this distribution corresponds to the Dirac deltafunction δ(x), a function which is localized at zero. We come back to this limit again at the end of the example. The Fourier transform of the above initial distribution is  Z ∞ Z ∞  x2 2 1 1 +ikx − −ikx− x 2 2 2σ dx = √ n̂0 (k) = √ dx. (111) e e 2σ 2πσ 2 −∞ 2πσ 2 −∞ Completing the square for the exponent i

x2 1 1 h 2 2 2 4 2 2 x + ikσ + k σ . + ikx = x + 2σ ikx = 2σ 2 2σ 2 2σ 2 (112) enables (111) to be rewritten as k2 σ 2 e− 2 n̂0 (k) = √ 2πσ 2 Z ∞ e− (x+ikσ 2 )2 2σ 2 dx. (113) −∞ To calculate the above integral, which involves a complex integrand, we use the Cauchy integral formula. It states that for a complex function f (z) ∈ C, z ∈ C, integration along a closed path in the complex plane is zero, provided that f (z) has no poles inside the path: I f (z)dz = 0 (114) Introducing the substitution z = x + ikσ 2 , dz = dx, the integral (113) can be rewritten as k2 σ 2 e− 2 n̂0 (k) = √ lim 2πσ 2 R→∞ R+ikσ 2 Z z2 e− 2σ2 dz (115) −R+ikσ 2 Let’s keep R finite for the moment. We can then think of the integral as one segment of a closed curve with rectangular shape: I z2 0 = e− 2σ2 dz Z R+ikσ 2 = e −R+ikσ 2 − z2 2σ 2 Z R dz + e R+ikσ 2 − z2 2σ 2 Z −R dz + − e R z2 2σ 2 Z −R+ikσ 2 dz + −R z2 e− 2σ2 dz (116) In the limit R → ∞, the second as well as the last integral vanish due to the exponential damping with large R, leading to Z R+ikσ 2 lim R→∞ −R+ikσ 2 − e z2 2σ 2 Z R dz = lim R→∞ −R z2 e− 2σ2 dz . (117) Inserting this into Eq. (115), we obtain k2 σ 2 e− 2 n̂0 (k) = √ 2πσ 2 Z ∞ x2 e− 2σ2 dx. (118) −∞ This means that we can basically drop the imaginary part in the original integral, Eq. (113). From 18.01 we know that Z ∞ √ 2 e−x = π, −∞ √ so by making a change of variable y = x/ 2σ 2 we have that n̂0 (k) = e− k2 σ 2 2 . (119) This result is worth keeping in mind: The Fourier transform of a Gaussian is a Gaussian. To obtain the full solution of the diffusion equation in real space, we have to insert n̂0 (k) into (104), Z ∞ 1 2 2 2 n(x, t) = eikx−Dk t−k σ /2 dk. (120) 2π −∞ We can do a bit of rearranging to get n(x, t) = 1 2π ∞ Z 2 2 eikx−k (Dt+σ /2) dk . (121) −∞ As before we complete the square for the exponent, )     ( 2 σ2 σ2 ix x2 2 k Dt + − ikx = Dt + k− + , 2 2 2(Dt + σ 2 /2) 4 (Dt + σ 2 /2)2 (122) so that the integral becomes − n(x, t) = e x2 4(Dt+σ 2 /2) 2π Z ∞ h −(Dt+σ 2 /2) k− e ix 2(Dt+σ 2 /2) i2 dk. (123) −∞ This is essentially the same integral that we had before, so we drop the imaginary part and change the integration variable, giving the result n(x, t) = p e − x2 4(Dt+σ 2 /2) 4π(Dt + σ 2 /2) . (124) This is the solution of the diffusion equation starting from a Gaussian distribution at time t = 0. Note: • Introducing σ̃ 2 = 2(Dt+σ 2 /2), the solution is √a Gaussian with a standard deviation σ̃, i.e. the width of the solution grows like Dt in time. Similarly, the amplitude decreases like √1Dt . • Let us substitute td = σ 2 /(2D). The solution then can be written as − x2 e 4D(t+td ) n(x, t) = p . 4πD(t + td ) (125) 2 σ Remember that in the limit td = 2D → 0 the initial condition (110) corresponds to a Dirac delta function. Thus an initially Gaussian distribution of particles that is diffusing may be viewed as having originated from a delta function a time td ago. Indeed, it can be shown that diffusion will cause any form of particle distribution initially localised about zero to eventually look like a Gaussian. 6.2 Green’s function method This method relies on another trick for representing the solution, that is somewhat more intuitive. Now, instead of representing n in a basis of plane wave states, we will express it as a basis of states which are localised in position. This is done by using the so-called Dirac delta function, denoted δ(x − x0 ). You should think of this of a large spike of unit area that is centered exactly at the position x0 . The definition of δ is that given any function4 f (x), Z ∞ f (x0 ) δ(x − x0 )dx0 = f (x). (126) −∞ We can represent the initial distribution of particles n(x, 0) = n0 (x) as a superposition of δ-functions Z ∞ (127) n0 (x) = n0 (x0 ) δ(x − x0 ) dx0 . −∞ This formula decomposes n0 into a continuous series of “spikes”. The idea is to then understand how each spike individually evolves and then superimpose the evolution of each 4 Intuitively, one can obtain the Dirac δ-function from the normalized Gaussian (110) by letting σ → 0. Derivatives of order n of the δ-function, denoted by δ (n) , can be defined by partial integration Z ∞ Z ∞ 0 (n) 0 0 n f (x ) δ (x − x ) dx = (−1) f (n) (x0 ) δ(x − x0 ) dx0 . −∞ −∞ The Fourier transformation of the δ-function is given by Z ∞ δ̂(k) = e−ikx δ(x) dx = 1. −∞ Applying the inverse transformation yields a useful integral representation of the Dirac δ-function Z ∞ Z ∞ 1 1 δ(x) = eikx δ̂(k) dk = eikx dk. 2π −∞ 2π −∞ spike to find the final density distribution. We define the Green’s function G(x − x0 , t) so that G(x − x0 , 0) = δ(x − x0 ), and Z ∞ n(x, t) = G(x − x0 , t) n0 (x0 ) dx0 . (128) −∞ Plugging this into the diffusion equation we see that Z ∞ Z ∞ 0 ∂ 2 G(x − x0 , t) 0 0 ∂G(x − x , t) 0 n0 (x ) n0 (x0 ) dx = D dx . ∂t ∂x2 −∞ −∞ (129) Thus G(x − x0 , t) obeys the diffusion equation and we have reduced the problem to the mathematics of solving the diffusion equation for the localised initial condition δ(x − x0 ). There are many ways of solving this problem. The one in most textbooks is to actually use the Fourier decomposition of δ(x − x0 ) and solve the equation in Fourier space, and then transform back into real space. This is an advisable procedure as the Fourier transform of δ is very simple5 . We will advocate another procedure however, that is more elegant and uses an idea that we will return to later in our studies of fluids. The idea is to use dimensional analysis to determine the solution. If we look at the diffusion equation ∂n ∂2n =D 2 ∂t ∂x (130) we see that roughly ‘∂/∂t’ ∼ ‘∂ 2 /∂x2 ’. (I’ve written this in quotes because there is a sense in which this equality is meaningless.) What I mean by it is that if you have a function n which obeys a diffusion equation, taking a single time derivative of the function gives a number of about the same size as when you take two spatial derivatives. This means √ that the characteristic length scale over which n varies is of order t. Now, since the initial distribution δ √ is perfectly localised, we expect that at time t, G(x − x0 ) will have a characteristic width t. Thus, we guess a (so-called) similarity solution   x − x0 0 √ . (131) G(x − x , t) = A(t) F t The time dependence of A(t) is determined by the conservation of particles. Since   Z ∞ Z ∞ √ Z ∞ x n dx = A(t) F √ dx = A(t) t F (y)dy (132) t −∞ −∞ −∞ √ must be constant in time (we have changed variables from x to y = x/ t), we see that A0 A(t) = √ t (133) for some constant A0 . Now let’s just plug in √ A0 G(x, t) = √ F (x/ t) t 5 Note that in the limit td = σ 2 /(2D) → 0 the initial condition (110) approaches a Dirac delta function, so we have already ‘solved’the problem: equation (125) tells us that an initially Gaussian distribution of particles that is diffusing may be viewed as having originated from a delta function a time td ago. into the diffusion equation. This gives us the following ordinary differential equation for F (y)   1 1 1 0 1 − F − yF = 3 DF 00 . (134) 3 2 2 t2 t2 Cancelling out the time factors and integrating this equation once gives 1 − F y = DF 0 . 2 This equation can be immediately integrated to give F (y) = F0 e−y (135) 2 /4D , and thus F0 (x−x0 )2 G(x − x0 , t) = √ e− 4Dt , t (136) √ R where the constant F0 = 1/ 4πD is determined by requiring that dx G = 1. Mean square displacement We would like to determine hx(t)2 i for a collection of particles starting at x0 = 0. Since G(x − x0 , t) is the solution of the diffusion equation with initial condition δ(x − x0 ), we can compute the mean square displacement from Z ∞ 2 x(t) = dx x2 G(x, t) −∞ Z ∞ x2 x2 = dx √ e− 4Dt 4πDt −∞ √ Z ∞ α 2 dx x2 e−αx (137) = √ π −∞ where α = 1/(4Dt). To evaluate the integral, note that r Z ∞ Z ∞ π d d π 2 −αx2 −αx2 . dx x e =− dx e =− = dα −∞ dα α 2α3/2 −∞ (138) which then gives 1 = 2Dt. (139) 2α We have thus recovered the fundamental result that the mean square displacement of Brownian particles grows linearly in time. x(t)2 = 6.3 Zero-flux solution: Sedimentation Consider spherical particles diffusing under the effect of a constant drift velocity u in one dimension, described by the conservation law ∂n ∂ = − Jx ∂t ∂x with current Jx = un − D ∂ n. ∂x (140) (141) Interpreting x > 0 as the distance from the bottom of a vessel and assuming reflective boundaries at x = 0, we may think of u arising from the effects of gravity u= −gm∗ 6πηa (142) where g is the gravitational acceleration, a and m∗ > 0 denote radius and effective bouyant mass6 of the particles, and η the viscosity of the fluid. The stationary zero-current solution, satisfying Jx = 0, is obtained by integrating un − D ∂ n = 0, ∂x (143) yielding an exponentially decaying density profile n(x) = Ce−x/λ (144) with characteristic sedimentation length λ=− 6πηaD D > 0. = gm∗ u (145) Sutherland and Einstein showed in 1905 that the diffusion constant D of a small particle moving in a fluid is given by kT D= , (146) 6πηa where k is the Boltzmann’s constant and T the temperature (measured on a Kelvin scale), implying that kT λ= . (147) gm∗ Note that particle shape (and mass density) enter through the buoyant mass m∗ .