1.1.5 Diffusion-Reaction. Case when both ends insulated

Solve

\[ \frac{\partial u}{\partial t}=\alpha \frac{\partial ^{2}u}{\partial x^{2}}-\beta u\qquad 0<x<\pi ,t>0 \] with, \(a>0,\beta >0\) \begin{align*} \frac{\partial u\left ( 0,t\right ) }{\partial x} & =0\\ \frac{\partial u\left ( \pi ,t\right ) }{\partial x} & =0 \end{align*}

And initial conditions \[ u\left ( x,0\right ) =x \]

Solution

Let \(u=X\left ( x\right ) T\left ( t\right ) \). Substituting into the PDE gives\[ T^{\prime }X=\alpha X^{\prime \prime }T-\beta XT \] Dividing by \(XT\neq 0\) gives\begin{align*} \frac{T^{\prime }}{T} & =\frac{\alpha X^{\prime \prime }}{X}-\beta \\ \frac{T^{\prime }}{T}+\beta & =\frac{\alpha X^{\prime \prime }}{X}\\ \frac{T^{\prime }}{\alpha T}+\frac{\beta }{\alpha } & =\frac{X^{\prime \prime }}{X} \end{align*}

Since each sides depends on different variable and both are equal, they must be equal to same constant, say \(-\lambda \)\[ \frac{T^{\prime }}{\alpha T}+\frac{\beta }{\alpha }=\frac{X^{\prime \prime }}{X}=-\lambda \] This gives two ODE’s to solve\begin{align} X^{\prime \prime }+\lambda X & =0\tag{1}\\ X^{\prime }\left ( 0\right ) & =0\nonumber \\ X^{\prime }\left ( \pi \right ) & =0\nonumber \end{align}

And\begin{align} \frac{T^{\prime }}{\alpha T}+\frac{\beta }{\alpha } & =-\lambda \tag{2}\\ T^{\prime }+\beta T & =-\lambda \alpha T\nonumber \\ T^{\prime }+T\left ( \lambda \alpha +\beta \right ) & =0\nonumber \end{align}

Starting with (1).  

Assuming \(\lambda <0\), the solution is \begin{align*} X\left ( x\right ) & =A\cosh \left ( \sqrt{\lambda }x\right ) +B\sinh \left ( \sqrt{\lambda }x\right ) \\ X^{\prime } & =\sqrt{\lambda }A\sinh \left ( \sqrt{\lambda }x\right ) +\sqrt{\lambda }B\cosh \left ( \sqrt{\lambda }x\right ) \end{align*}

Applying first B.C. gives \(0=\sqrt{\lambda }B\), hence \(B=0\). Therefore the solution becomes \(X\left ( x\right ) =A\cosh \left ( \sqrt{\lambda }x\right ) \) and \(X^{\prime }=\sqrt{\lambda }A\sinh \left ( \sqrt{\lambda }x\right ) \). Applying second B.C. gives \[ 0=\sqrt{\lambda }A\sinh \left ( \sqrt{\lambda }\pi \right ) \] But \(\sinh \left ( \sqrt{\lambda }\pi \right ) =0\) only when its argument is zero, which is not the case here. This means \(A=0\), leading to trivial solution. Therefore \(\lambda <0\) is not eigenvalue.

Assuming \(\lambda =0\). The solution is \(X=Ax+B\). Hence \(X^{\prime }=A\). Applying first B.C. Gives \(A=0\). Hence solution is \(X=B\). Second boundary condition gives no additional information. Therefore  \(X=A_{0}\) is the solution with \(\lambda =0\).

Assuming \(\lambda >0\). The solution is\begin{align*} X & =A\cos \left ( \sqrt{\lambda }x\right ) +B\sin \left ( \sqrt{\lambda }x\right ) \\ X^{\prime } & =-\sqrt{\lambda }A\sin \left ( \sqrt{\lambda }x\right ) +\sqrt{\lambda }B\cos \left ( \sqrt{\lambda }x\right ) \end{align*}

Applying first B.C. gives \(0=\sqrt{\lambda }B\), hence \(B=0\). and the solution becomes \(X=A\cos \left ( \sqrt{\lambda }x\right ) ,X^{\prime }=-\sqrt{\lambda }A\sin \left ( \sqrt{\lambda }x\right ) \). Applying second B.C. gives\[ 0=-\sqrt{\lambda }A\sin \left ( \sqrt{\lambda }\pi \right ) \] Therefore \begin{align*} \sin \left ( \sqrt{\lambda }\pi \right ) & =0\\ \sqrt{\lambda }\pi & =n\pi \qquad n=1,2,3,\cdots \\ \lambda & =n^{2}\qquad n=1,2,3,\cdots \end{align*}

Therefore the space solution is\begin{equation} X\left ( x\right ) =A_{0}+\sum _{n=1}^{\infty }A_{n}\cos \left ( nx\right ) \tag{3} \end{equation} Now the time ODE (2) is solved.\[ T^{\prime }+T\left ( \lambda \alpha +\beta \right ) =0 \] Integrating factor is \(e^{\int \lambda \alpha +\beta dt}=e^{\left ( \lambda \alpha +\beta \right ) t}\). Hence   \(\frac{d}{dt}\left ( Te^{\left ( \lambda \alpha +\beta \right ) t}\right ) =0\) or \(Te^{\left ( \lambda \alpha +\beta \right ) t}=c\) where \(c\) is constant. Therefore\[ T\left ( t\right ) =ce^{-\left ( \lambda \alpha +\beta \right ) t}\] For \(\lambda =0\), the solution is \[ T_{0}\left ( t\right ) =c_{0}e^{-\beta t}\] and for \(\lambda >0\), the solution is \[ T_{n}\left ( t\right ) =c_{n}e^{-\left ( n^{2}\alpha +\beta \right ) t}\] Hence the time domain solution is\begin{equation} T\left ( t\right ) =c_{0}e^{-\beta t}+\sum _{n=1}^{\infty }c_{n}e^{-\left ( n^{2}\alpha +\beta \right ) t} \tag{4} \end{equation} Combining (3,4) the solution is\begin{align} u_{n} & =X_{n}T_{n}\nonumber \\ u\left ( x,t\right ) & =\sum _{n=0}^{\infty }X_{n}T_{n}\nonumber \\ & =A_{0}+c_{0}e^{-\beta t}+\sum _{n=1}^{\infty }b_{n}\cos \left ( nx\right ) e^{-\left ( n^{2}\alpha +\beta \right ) t} \tag{5} \end{align}

where both constants from space and time eigenfunctions are combined into one constant \(b_{n}\) above in the sum. Now initial conditions are used to find coefficients. At \(t=0\)\begin{align*} u\left ( x,0\right ) & =A_{0}+c_{0}+\sum _{n=1}^{\infty }b_{n}\cos \left ( nx\right ) \\ x & =b_{0}+\sum _{n=1}^{\infty }b_{n}\cos \left ( nx\right ) \end{align*}

Where \(b_{0}=A_{0}+c_{0}\). Multiplying both sides by \(\cos \left ( mx\right ) \) and integrating gives\[ \int _{0}^{\pi }x\cos mxdx=\int _{0}^{\pi }\left [ b_{0}\cos mx+\sum _{n=1}^{\infty }b_{n}\cos mx\cos \left ( nx\right ) \right ] dx \] For \(m=0\), all terms in the sum in RHS vanish, and only \(b_{0}\) left\begin{align*} \int _{0}^{\pi }xdx & =\int _{0}^{\pi }b_{0}dx\\ & =b_{0}\pi \end{align*}

Hence\begin{equation} b_{0}=\frac{1}{\pi }\int _{0}^{\pi }xdx \tag{6} \end{equation} For \(n>0\)\begin{align*} \int _{0}^{\pi }x\cos mxdx & =\int _{0}^{\pi }\left [ b_{0}\cos mx+\sum _{n=1}^{\infty }b_{n}\cos mx\cos \left ( nx\right ) \right ] dx\\ & =\int _{0}^{\pi }b_{0}\cos mxdx+\sum _{n=1}^{\infty }b_{n}\int _{0}^{\pi }\cos mx\cos \left ( nx\right ) dx \end{align*}

But \(\int _{0}^{\pi }b_{0}\cos mxdx=b_{0}\int _{0}^{\pi }\cos mxdx=0\), and all terms in the sum vanish except for \(n=m\). The above becomes\begin{align} \int _{0}^{\pi }x\cos mxdx & =b_{n}\frac{\pi }{2}\nonumber \\ b_{n} & =\frac{2}{\pi }\int _{0}^{\pi }x\cos mxdx \tag{7} \end{align}

From (6), where this is the case for \(n=0\)\begin{align*} b_{0} & =\frac{1}{\pi }\left ( \frac{x^{2}}{2}\right ) _{0}^{\pi }\\ & =\frac{1}{2\pi }\left ( \pi ^{2}\right ) \\ & =\frac{\pi }{2} \end{align*}

But \(b_{0}=A_{0}+c_{0}\), Therefore \begin{align*} A_{0}+c_{0} & =\frac{\pi }{2}\\ A_{0} & =\left ( \frac{\pi }{2}-c_{0}\right ) \end{align*}

From (7), where this is the for the case \(n>0\) and using integration by part, \(\int _{0}^{\pi }x\cos nxdx\), let \(u=x,dv=\cos nx,\rightarrow du=1,v=\frac{\sin nx}{n}\), therefore\begin{align*} \int _{0}^{\pi }x\cos nxdx & =\left ( x\frac{\sin nx}{n}\right ) _{0}^{\pi }-\frac{1}{n}\int _{0}^{\pi }\sin nxdx\\ & =\frac{1}{n}\left ( \pi \sin n\pi -0\right ) +\frac{1}{n}\left ( \frac{\cos nx}{n}\right ) _{0}^{\pi }\\ & =0+\frac{1}{n^{2}}\left ( \cos nx\right ) _{0}^{\pi }\\ & =\frac{1}{n^{2}}\left ( \cos n\pi -1\right ) \end{align*}

Hence \begin{align*} b_{n} & =\frac{2}{\pi }\int _{0}^{\pi }x\cos mxdx\\ & =\frac{2}{\pi }\frac{\left ( -1^{n}-1\right ) }{n^{2}} \end{align*}

Therefore the final solution is, from (5), by combing all above results, becomes\begin{align*} u\left ( x,t\right ) & =A_{0}+c_{0}e^{-\beta t}+\sum _{n=1}^{\infty }b_{n}\cos \left ( nx\right ) e^{-\left ( n^{2}\alpha +\beta \right ) t}\\ & =\left ( \frac{\pi }{2}-c_{0}\right ) +c_{0}e^{-\beta t}+\sum _{n=1}^{\infty }\frac{2}{\pi }\frac{\left ( -1^{n}-1\right ) }{n^{2}}\cos \left ( nx\right ) e^{-\left ( n^{2}\alpha +\beta \right ) t}\\ & =\frac{\pi }{2}+c_{0}\left ( e^{-\beta t}-1\right ) +\frac{2}{\pi }\sum _{n=1}^{\infty }\frac{\left ( -1^{n}-1\right ) }{n^{2}}\cos \left ( nx\right ) e^{-\left ( n^{2}\alpha +\beta \right ) t} \end{align*}

Solution is not unique, since there is unknown \(c_{0}\). To find \(c_{0}\), the solvability condition for \(\nabla ^{2}u=0\) with Neumann B.C. is used, which says that total flux must be zero at steady state, which is the case here, since flux is zero as given (insulated). Since solvability condition is satisfied, then energy is conserved. Solution at \(t=\infty \) from above is\[ u\left ( x,\infty \right ) =\frac{1}{2}\pi -c_{0}\] Since energy is conserved then comparing it to energy at initial conditions\[ \int _{0}^{\pi }\rho cu\left ( x,0\right ) dx=\int _{0}^{\pi }\rho cu\left ( x,\infty \right ) dx \] But \(u\left ( x,0\right ) =x\), hence\begin{align*} \int _{0}^{\pi }xdx & =\int _{0}^{\pi }\left ( \frac{\pi }{2}-c_{0}\right ) dx\\ \frac{\pi ^{2}}{2} & =\left ( \frac{1}{2}\pi -c_{0}\right ) \pi \\ \frac{\pi }{2} & =\frac{\pi }{2}-c_{0}\\ c_{0} & =0 \end{align*}

Hence the final solution is\[ u\left ( x,t\right ) =\frac{\pi }{2}+\frac{2}{\pi }\sum _{n=1}^{\infty }\frac{\left ( -1^{n}-1\right ) }{n^{2}}\cos \left ( nx\right ) e^{-\left ( n^{2}\alpha +\beta \right ) t}\] On the behavior of coefficients. From the final solution, the transient behavior is exponentially damped harmonics (since \(n^{2}\alpha +\beta \) is positive) with coefficient on the harmonic of the form \[ \frac{\left ( -1^{n}-1\right ) }{n^{2}}=-2\left \{ 1,0,\frac{1}{9},0,\frac{1}{25},0,\frac{1}{36},\cdots \right \} \] In other words, it has \(\frac{1}{n^{2}}\) behavior. This is expected from Fourier series theory. The coefficients are determined by initial conditions, which is \(f\left ( x\right ) =x\). The cosine Fourier series of \(f\left ( x\right ) \) is made by first making an even extension, then making periodic extension. When doing so, one can see that \(f^{\prime }\left ( x\right ) \) has a jump discontinuity at \(x=0\). From Fourier series theory, if a function \(f\left ( x\right ) \) has a derivative such that \(f^{\prime }\left ( x\right ) \) has jump discontinuity, then the Fourier series of \(f\left ( x\right ) \) itself will have \(\frac{1}{n^{2}}\)behavior.

Example

Using \(\alpha =2,\beta =1\), here is animation for 3.5 seconds

Source code is