#### 1.1.5 Diﬀusion-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 diﬀerent 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 ﬁrst 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 ﬁrst 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 ﬁrst 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$$X\left ( x\right ) =A_{0}+\sum _{n=1}^{\infty }A_{n}\cos \left ( nx\right ) \tag{3}$$ 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$$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}$$ 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 ﬁnd coeﬃcients. 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$$b_{0}=\frac{1}{\pi }\int _{0}^{\pi }xdx \tag{6}$$ 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 ﬁnal 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 ﬁnd $$c_{0}$$, the solvability condition for $$\nabla ^{2}u=0$$ with Neumann B.C. is used, which says that total ﬂux must be zero at steady state, which is the case here, since ﬂux is zero as given (insulated). Since solvability condition is satisﬁed, 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 ﬁnal 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 coeﬃcients. From the ﬁnal solution, the transient behavior is exponentially damped harmonics (since $$n^{2}\alpha +\beta$$ is positive) with coeﬃcient 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 coeﬃcients are determined by initial conditions, which is $$f\left ( x\right ) =x$$. The cosine Fourier series of $$f\left ( x\right )$$ is made by ﬁrst 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