home

PDF (letter size

Contents

These are collection of PDE problems solved analytically and animated. Most of the animations where done in Mathematica, some in Maple and Matlab.

Chapter 1
Heat PDE

Animations moved to This page

Chapter 2
Wave PDE

Animations moved to This page

Chapter 3
Laplace and Poisson PDE

Animations moved to This page

Chapter 4
Burger’s PDE

Solve $$\begin{array}{}\text{(1)}& u_t+u{u}_x=D{u}_{xx}\end{array}$$

BC$$\begin{array}{rl}u(0,t)& =0t>0\\ u(L,t)& =0t>0\end{array}$$

Initial conditions$$u(x,0)=f\left(x\right)0<x<L$$

Where $D$ is the diffusion constant.

Solution

Using Cole-Hopf, let $$\begin{array}{}\text{(2)}& u(x,t)=-2D\frac{{\varphi }_x}{\varphi }\end{array}$$ where $\varphi \equiv \varphi (x,t)$. Rewriting equation (1) as

$$\begin{array}{rl}{u}_t& =D{u}_{xx}-u{u}_x\\ \text{(3)}& & ={(D{u}_x-\frac{u^2}{2})}_x\end{array}$$

Substituting (2) into (3) gives

$$\begin{array}{rl}{(-2D\frac{{\varphi }_x}{\varphi })}_t& ={[D{(-2D\frac{{\varphi }_x}{\varphi })}_x-\frac{1}{2}{(-2D\frac{{\varphi }_x}{\varphi })}^2]}_x\\ -2D{\left(\frac{{\varphi }_x}{\varphi }\right)}_t& ={(-2D^2{\left(\frac{{\varphi }_x}{\varphi }\right)}_x-2D^2{\left(\frac{{\varphi }_x}{\varphi }\right)}^2)}_x\\ \text{(4)}& -2D{\left(\frac{{\varphi }_x}{\varphi }\right)}_t& =-2D^2{({\left(\frac{{\varphi }_x}{\varphi }\right)}_x+{\left(\frac{{\varphi }_x}{\varphi }\right)}^2)}_x\end{array}$$

But $${\left(\frac{{\varphi }_x}{\varphi }\right)}_t=-\frac{1}{{\varphi }^2}{\varphi }_t{\varphi }_x+\frac{1}{\varphi }{\varphi }_{xt}$$

And

$${\left(\frac{{\varphi }_t}{\varphi }\right)}_x=-\frac{1}{{\varphi }^2}{\varphi }_x{\varphi }_t+\frac{1}{\varphi }{\varphi }_{tx}$$

Therefore ${\left(\frac{{\varphi }_x}{\varphi }\right)}_t={\left(\frac{{\varphi }_t}{\varphi }\right)}_x$. Using this in LHS of (4) gives

$$\begin{array}{}\text{(5)}& -2D{\left(\frac{{\varphi }_t}{\varphi }\right)}_x=-2D^2{({\left(\frac{{\varphi }_x}{\varphi }\right)}_x+{\left(\frac{{\varphi }_x}{\varphi }\right)}^2)}_x\end{array}$$

And $$\begin{array}{rl}{\left(\frac{{\varphi }_x}{\varphi }\right)}_x+{\left(\frac{{\varphi }_x}{\varphi }\right)}^2& =-\frac{1}{{\varphi }^2}{\varphi }_x^2+\frac{{\varphi }_{xx}}{\varphi }+\frac{{\varphi }_x^2}{{\varphi }^2}\\ & =\frac{{\varphi }_{xx}}{\varphi }\end{array}$$

Using the above in the RHS of (5) gives

$$\begin{array}{}\text{(6)}& -2D{\left(\frac{{\varphi }_t}{\varphi }\right)}_x=-2D^2{\left(\frac{{\varphi }_{xx}}{\varphi }\right)}_x\end{array}$$

Integrating both side w.r.t. $x$ gives

$$-2D\frac{{\varphi }_t}{\varphi }=-2D^2\frac{{\varphi }_{xx}}{\varphi }+G\left(t\right)$$

Where $G\left(t\right)$ is the constant of integration since $\varphi \equiv \varphi (x,t)$. The above simplifies to the heat PDE in $\varphi (x,t)$

$$\begin{array}{}\text{(6A)}& {\varphi }_t=D{\varphi }_{xx}+G\left(t\right)\varphi \end{array}$$

Let $$\begin{array}{}\text{(6B)}& \psi =\varphi e^{-\int G\left(t\right)dt}\end{array}$$

Then $$\begin{array}{rl}{\psi }_t& ={\varphi }_t{e}^{-\int G\left(t\right)dt}-\varphi G\left(t\right)e^{-\int G\left(t\right)dt}\\ & =({\varphi }_t-\varphi G\left(t\right))e^{-\int G\left(t\right)dt}\end{array}$$

But from (6A), we see that ${\varphi }_t-\varphi G\left(t\right)=D{\varphi }_{xx}$. Therefore the above becomes

$${\psi }_t=D{\varphi }_{xx}{e}^{-\int G\left(t\right)dt}$$

But from (6B), we see that ${\varphi }_{xx}{e}^{-\int G\left(t\right)dt}={\psi }_{xx}$, therefore the above becomes

$${\psi }_t=D{\psi }_{xx}$$

Which is the heat PDE. The original BC and initial conditions are now transformed to $\psi$ to solve the above. Since $u=-2D\frac{{\varphi }_x}{\varphi }$, then solving this first for $\varphi$

$$\begin{array}{rl}\frac{{\varphi }_x}{\varphi }& =-\frac{1}{2D}u\\ \frac{\partial \varphi }{\partial x}\frac{1}{\varphi }& =-\frac{1}{2D}u\\ \frac{d\varphi }{\varphi }& =-\frac{1}{2D}udx\end{array}$$

Integrating gives

$$\begin{array}{rl}\ln \varphi & =-\frac{1}{2D}{\int }_0^{x}uds+C_0\\ \varphi (x,t)& =C{e}^{-\frac{1}{2D}{\int }_0^{x}uds}\end{array}$$

Since $u=-2D\frac{{\varphi }_x}{\varphi }$, then the constant $C$ cancels out. Then we it can be set to any value as it does not affect the solution. Let $C=1$ and the above becomes

$$\begin{array}{}\text{(7)}& \varphi (x,t)=e^{-\frac{1}{2D}{\int }_0^{x}uds}\end{array}$$

(7) is now used to transform the initial conditions. When $u(x,0)=f\left(x\right)$ the above becomes

$$\begin{array}{rl}\varphi (x,0)& =e^{-\frac{1}{2D}{\int }_0^{x}u(s,0)ds}\\ & =e^{-\frac{1}{2D}{\int }_0^{x}f\left(s\right)ds}\end{array}$$

Since from (6B), $\psi (x,t)=\varphi e^{-\int G\left(t\right)dt}=\varphi e^{-{\int }_0^{t}G\left(s\right)ds}$ therefore

$$\begin{array}{rl}\psi (x,0)& =\varphi (x,0)e^0\\ & =\varphi (x,0)\\ & =e^{-\frac{1}{2D}{\int }_0^{x}f\left(s\right)ds}\end{array}$$

To transform the boundary conditions, $u=-2D\frac{{\varphi }_x}{\varphi }$ is used. When $u(0,t)=0$ then $0=-2D\frac{{\varphi }_x(0,t)}{\varphi (0,t)}$ or

$${\varphi }_x(0,t)=0$$

But $\psi =\varphi e^{-{\int }_0^{t}G\left(s\right)ds}$, then ${\psi }_x={\varphi }_x{e}^{-{\int }_0^{t}G\left(s\right)ds}$ and therefore

$$\begin{array}{rl}{\psi }_x(0,t)& ={\varphi }_x(0,t)e^{-{\int }_0^{t}G\left(s\right)ds}\\ & =0\end{array}$$

Similarly, when $u(L,t)=0$ then $0=-2D\frac{{\varphi }_x(L,t)}{\varphi (L,t)}$ or

$${\varphi }_x(L,t)=0$$

Which gives

$${\psi }_x(L,t)=0$$

Hence the heat PDE to solve is

$${\psi }_t=D{\psi }_{xx}$$

BC$$\begin{array}{rl}{\psi }_x(0,t)& =0t>0\\ {\psi }_x(L,t)& =0t>0\end{array}$$

Initial conditions$$\psi (x,0)=e^{-\frac{1}{2D}{\int }_0^{x}f\left(s\right)ds}0<x<L$$

The above heat PDE is now solved for $\psi (x,t)$. This solution is transformed back to $u(x,t)$. First using $\psi =\varphi e^{-\int G\left(t\right)dt}$ to find $\varphi (x,t)$, then using $u(x,t)=-2D\frac{{\varphi }_x}{\varphi }$, to find $u(x,t)$.

So in summary, there are two transformations needed. Going from $u(x,t)\to \varphi (x,t)$ uses Cole-Hopf. Going from $\varphi (x,t)\to \psi (x,t)$ uses $\psi (x,t)=\varphi (x,t)e^{-\int G\left(t\right)dt}.$ It is $\psi (x,t)$ which is solved as the heat PDE ${\psi }_t=D{\psi }_{xx}$ and not $\varphi (x,t)$, which is just an intermediate transformation.

4.1 Example 1

Let $L=2\pi ,0<x<2\pi ,D=\frac{1}{10},$$$u(x,0)=f\left(x\right)=\sin x$$ And boundary conditions $u(0,t)=u(2\pi ,0)=0$. From above, after carrying the forward transformation, the PDE to solve is found to be

$${\psi }_t=D{\psi }_{xx}$$

With transformed boundary conditions

$$\begin{array}{rl}{\psi }_x(0,t)& =0t>0\\ {\psi }_x(2\pi ,t)& =0t>0\end{array}$$

And transformed initial conditions

$$\psi (x,0)=e^{-\frac{1}{2D}\int \sin \left(x\right)dx}=e^{\frac{1}{2D}\cos \left(x\right)}$$

This heat PDE is standard and has known solution by separation of variables which is

$$\begin{array}{rl}\psi (x,t)& =c_0+\sum _{n=1}^{\mathrm{\infty }}{c}_n{e}^{-D{\lambda }_{n}t}\cos \left(\sqrt{{\lambda }_n}x\right)\\ {\lambda }_n& ={\left(\frac{n\pi }{L}\right)}^{2}n=1,2,3,\cdots \end{array}$$

Or, since $L=2\pi$,

$$\begin{array}{rl}\psi (x,t)& =c_0+\sum _{n=1}^{\mathrm{\infty }}{c}_n{e}^{-D\frac{n^2}{4}t}\cos \left(\frac{n}{2}x\right)\\ {\lambda }_n& =\frac{n^2}{4}n=1,2,3,\cdots \end{array}$$

Where $$\begin{array}{rl}{c}_0& =\frac{1}{L}{\int }_0^L\psi (x,0)dx\\ & =\frac{1}{2\pi }{\int }_0^{2\pi }{e}^{\frac{\cos \left(x\right)}{2D}}dx\\ & =\mathrm{BesselI}(0,\frac{1}{2D})\end{array}$$

And

$$\begin{array}{rl}{c}_n& =\frac{2}{L}{\int }_0^L\psi (x,0)\cos \left(\sqrt{{\lambda }_n}x\right)dx\\ & =\frac{1}{\pi }{\int }_0^{2\pi }{e}^{\frac{\cos \left(x\right)}{2D}}\cos \left(\frac{n}{2}x\right)dx\end{array}$$

The after integral has no closed form solution. Hence the solution is

$$\psi (x,t)=\mathrm{BesselI}(0,\frac{1}{2D})+\frac{1}{\pi }\sum _{n=1}^{\mathrm{\infty }}({\int }_0^{2\pi }{e}^{\frac{\cos \left(x\right)}{2D}}\cos \left(\frac{n}{2}x\right)dx)e^{-D\frac{n^2}{4}t}\cos \left(\frac{n}{2}x\right)$$

But $\psi =\varphi e^{-\int G\left(t\right)dt}$, therefore

$$\varphi (x,t)=\psi (x,t)e^{\int G\left(t\right)dt}$$

But I do not know what $G\left(t\right)$ is. This was used during the forward transformation only and was eliminated. So how to find $\varphi (x,t)$? Need to find $\varphi (x,t)$ to be able to find $u(x,t)$ from $u(x,t)=-2D\frac{{\varphi }_x}{\varphi }$.

Chapter 5
Misc. PDE’s

5.1 FitzHugh-Nagumo in 2D

This was solved numerically in Matlab.

5.1.1 Example 1

The equations to solve are the following on the unit square in 2D.$$\begin{array}{rl}\frac{\partial v(x,y,t)}{\partial t}& =D{\mathrm{\nabla }}^{2}v+(a-v)(v-1)v-w+I\\ \frac{\partial w(x,y,t)}{\partial t}& =ϵ(v-\gamma w)\end{array}$$

Using $a=0.1,\gamma =2,ϵ=0.005,D=5\times {10}^{-5},I=0,$ hence the PDE’s are$$\begin{array}{rl}\frac{\partial v(x,y,t)}{\partial t}& =(5\times {10}^{-5}){\mathrm{\nabla }}^{2}v+(0.1-v)(v-1)v-w\\ \frac{\partial w(x,y,t)}{\partial t}& =0.005(v-2w)\end{array}$$

Initial conditions, $t=0$$$\begin{array}{rl}v(x,y,0)& =\exp (-100(x^2+y^2))\\ w(x,y,0)& =0\end{array}$$

Boundary conditions are homogeneous Neumann for $v$. (I solved this numerically, fractional step method. ADI for the diffusion solve).

5.1.2 Example 2

The equations to solve are the following on the unit square in 2D.$$\begin{array}{rl}\frac{\partial v(x,y,t)}{\partial t}& =D{\mathrm{\nabla }}^{2}v+(a-v)(v-1)v-w+I\\ \frac{\partial w(x,y,t)}{\partial t}& =ϵ(v-\gamma w)\end{array}$$

Using $a=0.1,\gamma =2,ϵ=0.005,D=5\times {10}^{-5},I=0,$ hence the PDE’s are$$\begin{array}{rl}\frac{\partial v(x,y,t)}{\partial t}& =(5\times {10}^{-5}){\mathrm{\nabla }}^{2}v+(0.1-v)(v-1)v-w\\ \frac{\partial w(x,y,t)}{\partial t}& =0.005(v-2w)\end{array}$$

Initial conditions, $t=0$$$\begin{array}{rl}v(x,y,0)& =1-2x\\ w(x,y,0)& =0.05y\end{array}$$

Boundary conditions are homogeneous Neumann for $v$. (I solved this numerically, fractional step method. ADI for the diffusion solve). See my Matlab web page for source code.

Chapter 6
Appendix

6.1 Summary table

Heat PDE $\frac{\partial u}{\partial t}=k\frac{{\partial }^{2}u}{\partial x^2}$ in $1D$ (in a rod)

Heat PDE $\frac{\partial u}{\partial t}=\alpha \frac{{\partial }^{2}u}{\partial x^2}-\beta u$ in $1D$ (in a rod) with $\alpha ,\beta >0$ for $0<x<\pi$

(TO DO) Heat PDE for periodic conditions $u(-L)=u\left(L\right)$ and $\frac{\partial u(-L)}{\partial x}=\frac{\partial u\left(L\right)}{\partial x}$

$${\lambda }_n={\left(\frac{n\pi }{L}\right)}^2,n=1,2,3,\cdots$$$$u(x,t)=\overbrace{a_0}+\overbrace{\sum _{n=1}^{\mathrm{\infty }}{A}_n\cos \left(\sqrt{{\lambda }_n}x\right)e^{-k{\lambda }_{n}t}+\sum _{n=1}^{\mathrm{\infty }}{B}_n\sin \left(\sqrt{{\lambda }_n}x\right)e^{-k{\lambda }_{n}t}}$$

6.2 Using Mathematica to obtain the eigenvalues and eigenfunctions for heat PDE in 1D

6.2.1 $u(0)=0,u(L)=0$

For eigenvalue

ClearAll[y,x,L]; 
op={-y''[x]  ,DirichletCondition[y[x]==0,x==0], 
DirichletCondition[y[x]==0,x==L]}; 
(*or simply*) 
op={-y''[x],DirichletCondition[y[x]==0,True]}; 
eig=DEigenvalues[op,y[x],{x,0,L},5]
 

$$\{\frac{{\pi }^2}{L^2},\frac{4{\pi }^2}{L^2},\frac{9{\pi }^2}{L^2},\frac{16{\pi }^2}{L^2},\frac{25{\pi }^2}{L^2}\}$$

FindSequenceFunction[eig,n]
 

$$\frac{{\pi }^2{n}^2}{L^2}$$

For eigenfunctions

ClearAll[y,x,L]; 
op={-y''[x]  ,DirichletCondition[y[x]==0,x==0], 
DirichletCondition[y[x]==0,x==L]}; 
eigf=Last@DEigensystem[op,y[x],{x,0,L},5]
 

$$\{\sin \left(\frac{\pi x}{L}\right),\sin \left(\frac{2\pi x}{L}\right),\sin \left(\frac{3\pi x}{L}\right),\sin \left(\frac{4\pi x}{L}\right),\sin \left(\frac{5\pi x}{L}\right)\}$$

FullSimplify[FindSequenceFunction[eigf,n]];
 

$$\sin \left(\frac{\pi nx}{L}\right)$$

6.2.2 $u^{\prime }(0)=0,u^{\prime }(L)=0$

For eigenvalue

ClearAll[y,x,L]; 
op={-y''[x]+NeumannValue[0,True]}; 
eig=DEigenvalues[op,y[x],{x,0,L},5]
 

$$\{0,\frac{{\pi }^2}{L^2},\frac{4{\pi }^2}{L^2},\frac{9{\pi }^2}{L^2},\frac{16{\pi }^2}{L^2}\}$$

FindSequenceFunction[eig,n]
 

$$\frac{{\pi }^2(n-1{)}^2}{L^2}$$

ClearAll[y,x,L]; 
op={-y''[x]+NeumannValue[0,True]}; 
eigf=Last@DEigensystem[op,y[x],{x,0,L},5]
 

$$\{1,\cos \left(\frac{\pi x}{L}\right),\cos \left(\frac{2\pi x}{L}\right),\cos \left(\frac{3\pi x}{L}\right),\cos \left(\frac{4\pi x}{L}\right)\}$$

FullSimplify[FindSequenceFunction[eigf, n]];
 

$$\cos \left(\frac{\pi (n-1)x}{L}\right)$$

6.2.3 $u^{\prime }(0)=0,u(L)=0$

For eigenvalue

ClearAll[y,x,L]; 
op={-y''[x]+NeumannValue[0,x==0],DirichletCondition[y[x]==0,x==L]}; 
eig=DEigenvalues[op,y[x],{x,0,L},6] 
\begin{MMAinline} 
 
\[ 
\left\{\frac{\pi^2}{4 L^2},\frac{9 \pi^2}{4 L^2},\frac{25 \pi^2}{4 L^2}% 
,\frac{49 \pi^2}{4 L^2},\frac{81 \pi^2}{4 L^2},\frac{121 \pi^2}{4 L^2}\right\} 
\] 
 
\begin{MMAinline} 
Simplify[FindSequenceFunction[eig,n]]
 

$$\frac{{\pi }^2(1-2n{)}^2}{4L^2}$$

For eigenfunctions

eigf=Last@DEigensystem[op,y[x],{x,0,L},7]
 

$$\{\cos \left(\frac{\pi x}{2L}\right),\cos \left(\frac{3\pi x}{2L}\right),\cos \left(\frac{5\pi x}{2L}\right),\cos \left(\frac{7\pi x}{2L}\right),\cos \left(\frac{9\pi x}{2L}\right)\}$$

FullSimplify[FindSequenceFunction[eigf,n]];
 

$$\cos \left(\frac{\pi (1-2n)x}{2L}\right)$$

6.2.4 $u(0)=0,u^{\prime }(L)=0$

Eigenvalue are the same as above

ClearAll[y,x,L]; 
op={-y''[x]+NeumannValue[0,x==L],DirichletCondition[y[x]==0,x==0]}; 
eig=DEigenvalues[op,y[x],{x,0,L},6]
 

$$\{\frac{{\pi }^2}{4L^2},\frac{9{\pi }^2}{4L^2},\frac{25{\pi }^2}{4L^2},\frac{49{\pi }^2}{4L^2},\frac{81{\pi }^2}{4L^2},\frac{121{\pi }^2}{4L^2}\}$$

FindSequenceFunction[eig,n]
 

$$\frac{{\pi }^2(1-2n{)}^2}{4L^2}$$

For eigenfunctions

eigf=Last@DEigensystem[op,y[x],{x,0,L},6]
 

$$\{\sin \left(\frac{\pi x}{2L}\right),\sin \left(\frac{3\pi x}{2L}\right),\sin \left(\frac{5\pi x}{2L}\right),\sin \left(\frac{7\pi x}{2L}\right),\sin \left(\frac{9\pi x}{2L}\right),\sin \left(\frac{11\pi x}{2L}\right)\}$$

FullSimplify[FindSequenceFunction[eigf, n]];
 

$$-\sin \left(\frac{\pi (1-2n)x}{2L}\right)$$

6.2.5 $u(0)=0,u\left(L\right)+u^{\prime }(L)=0$

For eigenvalue

(*can only find them numerially*) 
ClearAll[y,x]; 
op={-y''[x]+NeumannValue[y[x],x==1],DirichletCondition[y[x]==0,x==0]}; 
eig=DEigenvalues[op,y[x],{x,0,1},6]//N 
(* {4.11586,24.1393,63.6591,122.889,201.851,300.55}*) 
NSolve[Tan[Sqrt[lam]]+Sqrt[lam]==0&& 0<lam<130,lam] 
(*{{lam->4.11586},{lam->24.1393},{lam->63.6591},{lam->122.889}}*)
 

For eigenfunctions

ClearAll[y,x]; 
op={-y''[x]+NeumannValue[y[x],x==1],DirichletCondition[y[x]==0,x==0]}; 
eig=Last@DEigensystem[op,y[x],{x,0,1},6]//N
 

$$\{\sin (2.02876x),\sin (4.91318x),\sin (7.97867x),\sin (11.0855x),\sin (14.2074x),\sin (17.3364x)\}$$

6.2.6 $u\left(0\right)+u^{\prime }(0)=0,u^{\prime }(L)=0$

For eigenvalue

 
(*can only find them numerially. Notice the sign difference now. *) 
ClearAll[y,x]; 
op={-y''[x]+NeumannValue[-y[x],x==0],DirichletCondition[y[x]==0,x==1]}; 
eig=DEigenvalues[op,y[x],{x,0,1},6]//N 
(* {0.,20.1908,59.6814,118.914,197.924,296.774}*) 
NSolve[Tan[Sqrt[lam]]-Sqrt[lam]==0&& 0<=lam<130,lam] 
(*{{lam->0.},{lam->20.1907},{lam->59.6795},{lam->118.9},{lam->197.858}% 
,{lam->296.554}}*)
 

For eigenfunctions, mathematica only gives plots, so not shown.

6.2.7 $u(-L)=0,u\left(L\right)=0$

For eigenvalue

ClearAll[y,x,L]; 
op={-y''[x],DirichletCondition[y[x]==0,x==-L],DirichletCondition[y[x]==0,x==L]}% 
; 
eig=DEigenvalues[op,y[x],{x,-L,L},6]
 

$$\{\frac{{\pi }^2}{4L^2},\frac{{\pi }^2}{L^2},\frac{9{\pi }^2}{4L^2},\frac{4{\pi }^2}{L^2},\frac{25{\pi }^2}{4L^2},\frac{9{\pi }^2}{L^2}\}$$