9.5: Solution of the Diffusion Equation

Homogeneous Boundary Conditions

We consider one dimensional diffusion in a pipe of length \(L\), and solve the diffusion equation for the concentration \(u(x, t)\), \[\label{eq:1}u_t=Du_{xx},\quad 0\leq x\leq L,\quad t>0.\]

Both initial and boundary conditions are required for a unique solution. That is, we assume the initial concentration distribution in the pipe is given by \[\label{eq:2}u(x,0)=f(x),\quad 0\leq x\leq L.\]

Furthermore, we assume that boundary conditions are given at the ends of the pipes. When the concentration value is specified at the boundaries, the boundary conditions are called Dirichlet boundary conditions. As the simplest example, we assume here homogeneous Dirichlet boundary conditions, that is zero concentration of dye at the ends of the pipe, which could occur if the ends of the pipe open up into large reservoirs of clear solution, \[\label{eq:3}u(0,t)=0,\quad u(L,t)=0,\quad t>0.\]

We will later also discuss inhomogeneous Dirichlet boundary conditions and homogeneous Neumann boundary conditions, for which the derivative of the concentration is specified to be zero at the boundaries. Note that if \(f(x)\) is identically zero, then the trivial solution \(u(x, t) = 0\) satisfies the differential equation and the initial and boundary conditions and is therefore the unique solution of the problem. In what follows, we will assume that \(f(x)\) is not identically zero so that we need to find a solution different than the trivial solution.

The solution method we use is called separation of variables. We assume that \(u(x, t)\) can be written as a product of two other functions, one dependent only on position \(x\) and the other dependent only on time \(t\). That is, we make the ansatz \[\label{eq:4}u(x,t)=X(x)T(t).\]

Whether this ansatz will succeed depends on whether the solution indeed has this form. Substituting \(\eqref{eq:4}\) into \(\eqref{eq:1}\), we obtain \[XT'=DX''T,\nonumber\] which we rewrite by separating the \(x\) and \(t\) dependence to opposite sides of the equation:

\[\frac{X''}{X}=\frac{1}{D}\frac{T'}{T}.\nonumber\]

The left hand side of this equation is independent of \(t\) and the right hand side is independent of \(x\). Both sides of this equation are therefore independent of both \(x\) and \(t\) and equal to a constant. Introducing \(−\lambda\) as the separation constant, we have \[\frac{X''}{X}=\frac{1}{D}\frac{T'}{T}=-\lambda ,\nonumber\] and we obtain the two ordinary differential equations \[\label{eq:5}X''+\lambda X=0,\quad T'+\lambda DT=0.\]

Because of the boundary conditions, we must first consider the equation for \(X(x)\). To solve, we need to determine the boundary conditions at \(x = 0\) and \(x = L\). Now, from \(\eqref{eq:3}\) and \(\eqref{eq:4}\), \[u(0,t)=X(0)T(t)=0,\quad t>0.\nonumber\]

Since \(T(t)\) is not identically zero for all \(t\) (which would result in the trivial solution for \(u\)), we must have \(X(0) = 0\). Similarly, the boundary condition at \(x = L\) requires \(X(L) = 0\). We therefore consider the two-point boundary value problem \[\label{eq:6}X''+\lambda X=0,\quad X(0)=X(L)=0.\]

The equation given by \(\eqref{eq:6}\) is called an ode eigenvalue problem. Clearly, the trivial solution \(X(x) = 0\) is a solution. Nontrivial solutions exist only for discrete values of \(\lambda\). These discrete values of \(\lambda\) and the corresponding functions \(X(x)\) are called the eigenvalues and eigenfunctions of the differential equation.

Since the form of the general solution of the ode depends on the sign of \(\lambda\), we consider in turn the cases \(\lambda > 0\), \(\lambda < 0\) and \(\lambda = 0\). For \(\lambda > 0\), we write \(\lambda =\mu^2\) and determine the general solution of \[X''+\mu^2X=0\nonumber\] to be \[X(x)=A\cos\mu x+B\sin\mu x.\nonumber\]

Applying the boundary condition at \(x = 0\), we find \(A = 0\). The boundary condition at \(x = L\) then yields \[B\sin\mu L=0.\nonumber\]

The solution \(B = 0\) results in the trivial solution for \(u\) and can be ruled out. Therefore, we must have \[\sin\mu L=0,\nonumber\] which is an equation for \(\mu\). The solutions are \[\mu =n\pi /L,\nonumber\] where \(n\) is an integer. We have thus determined the eigenvalues \(\lambda =\mu^2 > 0\) to be \[\label{eq:7}\lambda_n=(n\pi /L)^2,\quad n=1,2,3,\ldots ,\] with corresponding eigenfunctions \[\label{eq:8}X_n=\sin (n\pi x/L).\]

For \(\lambda < 0\), we write \(\lambda = −\mu^2\) and determine the general solution of \[X''-\mu^2X=0\nonumber\] to be \[X(x)=A\cosh \mu x+B\sinh \mu x,\nonumber\] where we have previously introduced the hyperbolic sine and cosine functions in §4.4. Applying the boundary condition at \(x = 0\), we find \(A = 0\). The boundary condition at \(x = L\) then yields \[B\sinh\mu L=0,\nonumber\] which for \(\mu\neq 0\) has only the solution \(B = 0\). Therefore, there is no nontrivial solution for \(u\) with \(\lambda < 0\). Finally, for \(\lambda = 0\), we have \[X''=0,\nonumber\] with general solution \[X(x)=A+Bx.\nonumber\]

The boundary condition at \(x = 0\) and \(x = L\) yields \(A = B = 0\) so again there is no nontrivial solution for \(u\) with \(\lambda = 0\).

We now turn to the equation for \(T(t)\). The equation corresponding to the eigenvalue \(\lambda_n\), using \(\eqref{eq:7}\), is given by \[T'+\left(n^2\pi^2 D/L^2\right)T=0,\nonumber\] which has solution proportional to \[\label{eq:9}T_n=e^{-n^2\pi ^2Dt/L^2}.\]

Therefore, multiplying the solutions given by \(\eqref{eq:8}\) and \(\eqref{eq:9}\), we conclude that the functions \[\label{eq:10}u_n(x,t)=\sin (n\pi x/L)e^{-n^2\pi^2 Dt/L^2}\] satisfy the pde given by \(\eqref{eq:1}\) and the boundary conditions given by \(\eqref{eq:3}\) for every positive integer \(n\).

The principle of linear superposition for homogeneous linear differential equations then states that the general solution to \(\eqref{eq:1}\) and \(\eqref{eq:3}\) is given by \[\begin{align}u(x,t)&=\sum\limits_{n=1}^\infty b_nu_n(x,t)\nonumber \\ &=\sum\limits_{n=1}^\infty b_n\sin (n\pi x/L)e^{-n^2\pi^2 Dt/L^2}.\label{eq:11}\end{align}\]

The final solution step is to satisfy the initial conditions given by \(\eqref{eq:2}\). At \(t = 0\), we have \[\label{eq:12}f(x)=\sum\limits_{n=1}^\infty b_n\sin (n\pi x/L).\]

We immediately recognize \(\eqref{eq:12}\) as a Fourier sine series (9.4.4) for an odd function \(f(x)\) with period \(2L\). Equation \(\eqref{eq:12}\) is a periodic extension of our original \(f(x)\) defined on \(0 ≤ x ≤ L\), and is an odd function because of the boundary condition \(f(0) = 0\). From our solution for the coefficients of a Fourier sine series (9.4.3), we determine \[\label{eq:13}b_n=\frac{2}{L}\int_0^L f(x)\sin\frac{n\pi x}{L}dx.\]

Thus the solution to the diffusion equation with homogeneous Dirichlet boundary conditions defined by \(\eqref{eq:1}\), \(\eqref{eq:2}\) and \(\eqref{eq:3}\) is given by \(\eqref{eq:11}\) with the \(b_n\) coefficients computed from \(\eqref{eq:13}\).

Example \(\PageIndex{1}\)

Determine the concentration of a dye in a pipe of length \(L\), where the dye has total initial mass \(M_0\) and is initially concentrated at the center of the pipe, and the ends of the pipe are held at zero concentration.

Solution

The governing equation for concentration is the diffusion equation. We model the initial concentration of the dye by a delta-function centered at \(x = L/2\), that is, \(u(x, 0) = f(x) = M_0\delta (x − L/2)\). Therefore, from \(\eqref{eq:13}\), \[\begin{aligned}b_n&=\frac{2}{L}\int_0^L M_0\delta (x-\frac{L}{2})\sin\frac{n\pi x}{L}dx \\ &=\frac{2M_0}{L}\sin (n\pi /2) \\ &=\left\{\begin{array}{ll} 2M_0/L &\text{if }n=1,5,9,\ldots ; \\ -2M_0/L&\text{if }n=3,7,11,\ldots ; \\ 0&\text{if }n=2,4,6,\ldots .\end{array}\right.\end{aligned}\]

With \(b_n\) determined, the solution for \(u(x, t)\) given by \(\eqref{eq:11}\) can be written as \[u(x,t)=\frac{2M_0}{L}\sum\limits_{n=0}^\infty (-1)^n\sin\left(\frac{(2n+1)\pi x}{L}\right)e^{-(2n+1)^2\pi^2 Dt/L^2}.\nonumber\]

When \(t ≫ L^2/D\), the leading-order term in the series is a good approximation and is given by \[u(x,t)\approx\frac{2M_0}{L}\sin (\pi x/L)e^{-\pi ^2Dt/L^2}.\nonumber\]

The mass of the dye in the pipe is decreasing in time, diffusing into the reservoirs located at both ends. The total mass in the pipe at time \(t\) can be found from \[M(t)=\int_0^L u(x,t)dx,\nonumber\] and since \[\int_0^L\sin (\pi x/L)dx=\frac{2L}{\pi},\nonumber\] we have for large times \[M(t)=\frac{4M_0}{\pi}e^{-\pi^2 Dt/L^2}.\nonumber\]

Inhomogeneous Boundary Conditions

Consider a diffusion problem where one end of the pipe has dye of concentration held constant at \(C_1\) and the other held constant at \(C_2\), which could occur if the ends of the pipe had large reservoirs of fluid with different concentrations of dye. With \(u(x, t)\) the concentration of dye, the boundary conditions are given by \[u(0,t)=C_1,\quad u(L,t)=C_2,\quad t>0.\nonumber\]

The concentration \(u(x, t)\) satisfies the diffusion equation with diffusivity \(D\):

\[u_t=Du_{xx}.\nonumber\]

If we try to solve this problem directly using separation of variables, we will run into trouble. Applying the inhomogeneous boundary condition at \(x = 0\) directly to the ansatz \(u(x, t) = X(x)T(t)\) results in \[u(0,t)=X(0)T(t)=C_1;\nonumber\] so that \[X(0)=C_1 /T(t).\nonumber\] However, our separation of variables ansatz assumes \(X(x)\) to be independent of \(t\)! We therefore say that inhomogeneous boundary conditions are not separable.

The proper way to solve a problem with inhomogeneous boundary conditions is to transform it into another problem with homogeneous boundary conditions. As \(t\to\infty\), we assume that a stationary concentration distribution \(v(x)\) will attain, independent of \(t\). Since \(v(x)\) must satisfy the diffusion equation, we have \[v''(x)=0,\quad 0\leq x\leq L,\nonumber\] with general solution \[v(x)=A+Bx.\nonumber\]

Since \(v(x)\) must satisfy the same boundary conditions of \(u(x, t)\), we have \(v(0) = C_1\) and \(v(L) = C_2\), and we determine \(A = C_1\) and \(B = (C_2 − C_1)/L\).

We now express \(u(x, t)\) as the sum of the known asymptotic stationary concentration distribution \(v(x)\) and an unknown transient concentration distribution \(w(x, t)\):

\[u(x,t)=v(x)+w(x,t).\nonumber\]

Substituting into the diffusion equation, we obtain \[\frac{\partial}{\partial t}(v(x)+w(x,t))=D\frac{\partial ^2}{\partial x^2}(v(x)+w(x,t))\nonumber\] \[\text{or}\nonumber\] \[w_t=Dw_{xx},\nonumber\] since \(v_t = 0\) and \(v_{xx} = 0\). The boundary conditions satisfied by \(w\) are \[\begin{aligned} w(0,t)&=u(0,t)-v(0)=0, \\ w(L,t)&=u(L,t)-v(L)=0,\end{aligned}\] so that \(w\) is observed to satisfy homogeneous boundary conditions. If the initial conditions are given by \(u(x, 0) = f(x)\), then the initial conditions for \(w\) are \[\begin{aligned} w(x,0)&=u(x,0)-v(x) \\ &=f(x)-v(x).\end{aligned}\]

The resulting equations may then be solved for \(w(x, t)\) using the technique for homogeneous boundary conditions, and \(u(x, t)\) subsequently determined.

Pipe with Closed Ends

There is no diffusion of dye through the ends of a sealed pipe. Accordingly, the mass flux of dye through the pipe ends, given by (9.1.1), is zero so that the boundary conditions on the dye concentration \(u(x, t)\) becomes \[\label{eq:14}u_x(0,t)=0,\quad u_x(L,t)=0,\quad t>0,\] which are known as homogeneous Neumann boundary conditions. Again, we apply the method of separation of variables and as before, we obtain the two ordinary differential equations given by \(\eqref{eq:5}\). Considering first the equation for \(X(x)\), the appropriate boundary conditions are now on the first derivative of \(X(x)\), and we must solve \[\label{eq:15}X''+\lambda X=0,\quad X'(0)=X'(L)=0.\]

Again, we consider in turn the cases \(\lambda > 0,\: \lambda < 0\) and \(\lambda = 0\). For \(\lambda > 0\), we write \(\lambda = \mu^2\) and determine the general solution of \(\eqref{eq:15}\) to be \[X(x)=A\cos\mu x+B\sin\mu x,\nonumber\] so that taking the derivative \[X'(x)=-\mu A\sin \mu x+\mu B\cos\mu x.\nonumber\]

Applying the boundary condition \(X' (0) = 0\), we find \(B = 0\). The boundary condition at \(x = L\) then yields \[-\mu A\sin\mu L=0.\nonumber\]

The solution \(A = 0\) results in the trivial solution for \(u\) and can be ruled out. Therefore, we must have \[\sin\mu L=0,\nonumber\] with solutions \[\mu =n\pi /L,\nonumber\] where \(n\) is an integer. We have thus determined the eigenvalues \(\lambda =\mu^2 > 0\) to be \[\label{eq:16}\lambda_n=(n\pi /L)^2,\quad n=1,2,4,\ldots ,\] with corresponding eigenfunctions \[\label{eq:17}X_n=\cos (n\pi x/L).\]

For \(\lambda < 0\), we write \(\lambda = −\mu^2\) and determine the general solution of \(\eqref{eq:15}\) to be \[X(x)+A\cosh\mu x+B\sinh\mu x,\nonumber\] so that taking the derivative \[X'(x)=\mu A\sinh\mu x+\mu B\cosh\mu x.\nonumber\]

Applying the boundary condition \(X' (0) = 0\) yields \(B = 0\). The boundary condition \(X' (L) = 0\) then yields \[\mu A\sinh\mu L=0,\nonumber\] which for \(\mu\neq 0\) has only the solution \(A = 0\). Therefore, there is no nontrivial solution for \(u\) with \(\lambda < 0\). Finally, for \(\lambda = 0\), the general solution of \(\eqref{eq:15}\) is \[X(x)=A+Bx,\nonumber\] so that taking the derivative \[X'(x)=B.\nonumber\]

The boundary condition \(X' (0) = 0\) yields \(B = 0\); \(X' (L) = 0\) is then trivially satisfied. Therefore, we have an additional eigenvalue and eigenfunction given by \[\lambda_0=0,\quad X_0(x)=1,\nonumber\] which can be seen as extending the formula obtained for eigenvalues and eigenvectors for positive \(\lambda\) given by \(\eqref{eq:16}\) and \(\eqref{eq:17}\) to \(n = 0\).

We now turn to the equation for \(T(t)\). The equation corresponding to eigenvalue \(\lambda_n\), using \(\eqref{eq:16}\), is given by \[T'+\left(n^2\pi^2D/L^2\right)T=0,\nonumber\] which has solution proportional to \[\label{eq:18}T_n=e^{-n^2\pi^2 Dt/L^2},\] valid for \(n = 0,\: 1,\: 2,\ldots\). Therefore, multiplying the solutions given by \(\eqref{eq:17}\) and \(\eqref{eq:18}\), we conclude that the functions \[\label{eq:19}u_n(x,t)=\cos (n\pi x/L)e^{-n^2\pi^2 Dt/L^2}\] satisfy the pde given by \(\eqref{eq:1}\) and the boundary conditions given by \(\eqref{eq:14}\) for every nonnegative integer \(n\).

The principle of linear superposition then yields the general solution as \[\begin{align}u(x,t)&=\sum\limits_{n=0}^\infty c_nu_n(x,t) \nonumber \\ &=\frac{a_0}{2}+\sum\limits_{n=1}^\infty a_n\cos (n\pi x/L)e^{-n^2\pi^2 Dt/L^2},\label{eq:20}\end{align}\] where we have redefined the constants so that \(c_0 = a_0/2\) and \(c_n = a_n\), \(n = 1,\: 2,\: 3,\ldots\). The final solution step is to satisfy the initial conditions given by \(\eqref{eq:2}\). At \(t = 0\), we have \[\label{eq:21}f(x)=\frac{a_0}{2}+\sum\limits_{n=1}^\infty a_n\cos (n\pi x/L),\] which we recognize as a Fourier cosine series (9.4.2) for an even function \(f(x)\) with period \(2L\). We have obtained a cosine series for the periodic extension of \(f(x)\) because of the boundary condition \(f'(0) = 0\), which is satisfied by an even function with continuous first derivative. From our solution (9.4.1) for the coefficients of a Fourier cosine series, we determine \[\label{eq:22}a_n=\frac{2}{L}\int_0^L f(x)\cos\frac{n\pi x}{L}dx.\]

Thus the solution to the diffusion equation with homogenous Neumann boundary conditions defined by \(\eqref{eq:1}\), \(\eqref{eq:2}\) and \(\eqref{eq:14}\) is given by \(\eqref{eq:20}\) with the coefficients computed from \(\eqref{eq:22}\).

Example \(\PageIndex{2}\)

Determine the concentration of a dye in a pipe of length \(L\), where the dye has total initial mass \(M_0\) and is initially concentrated at the center of the pipe, and the ends of the pipe are sealed

Solution

Again we model the initial concentration of the dye by a delta-function centered at \(x = L/2\). From \(\eqref{eq:22}\), \[\begin{aligned}a_n&=\frac{2}{L}\int_0^L M_0\delta \left(x-\frac{L}{2}\right) \cos\frac{n\pi x}{L}dx \\ &=\frac{2M_0}{L}\cos (n\pi /2) \\ &=\left\{\begin{array}{ll}2M_0/L &\text{if }n=0,4,8,\ldots ; \\ -2M_0/L&\text{if }n=2,6,10,\ldots ; \\ 0&\text{if }n=1,3,5\ldots .\end{array}\right.\end{aligned}\]

The first two terms in the series for \(u(x, t)\) are given by \[u(x,t)=\frac{M_0}{L}\left[1-2\cos (2\pi x/L)e^{-4\pi^2 Dt/L^2}+\ldots\right].\nonumber\]

Notice that as \(t\to\infty\), \(u(x, t) → M_0/L\): the dye mass is conserved in the pipe (since the pipe ends are sealed) and eventually diffused uniformly throughout the pipe of length \(L\).