Added May 26, 2019.
Problem 7.9.4 (a) from Richard Haberman Applied Partial Differential Equations, 4th edition.
Solve Heat PDE \(u_t = k \nabla ^2 u \) inside cylinder with radius \(a\) and height \(H\) with initial conditions \(u(r,\theta ,z,0)=f(r,z)\) independent of \(\theta \) if the boundary conditions are \(u(r,\theta ,0,t)=0\), \(u(r,\theta ,H,t)=0\), \(u(a,\theta ,z,t)=0\).
Since it says independent of \(\theta \), will use the PDE as \begin {align*} u_t = k \left (u_{rr} + \frac {1}{r} u_r + u_{zz} \right ) \end {align*}
Instead of the full Laplacian \begin {align*} u_t = k \left (u_{rr} + \frac {1}{r} u_r + \frac {1}{r^2} u_{\theta \theta } + u_{zz} \right ) \end {align*}
Mathematica ✗
ClearAll["Global`*"]; lap = Laplacian[u[r, z, t], {r, theta, z}, "Cylindrical"]; bc = {u[r, 0, t] == 0, u[r, H, t] == 0, u[a, z,t] == 0}; ic = u[r, z, 0] == f[r,z]; sol = AbsoluteTiming[TimeConstrained[DSolve[{D[u[r,z,t],t]==k*lap, bc, ic}, u[r, z,t], {r, z,t}, Assumptions -> {a > 0, r < a, H > 0}], 60*10]];
Failed
Maple ✓
restart; lap:=VectorCalculus:-Laplacian(u(r,z,t),'cylindrical'[r,theta,z]); bc := u(r,0,t)=0,u(r,H,t)=0, u(a,z,t)=0; ic := u(r,z,0) = f(r,z); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([diff(u(r,z,t),t) = k*lap,bc,ic],u(r,z,t)) assuming a>0,r<a,H>0,k>0),output='realtime'));
Failed to convert to latex
Hand solution
Solve \(u_{t}=k\nabla ^{2}u\) inside cylinder with radius a and height \(H\) with initial conditions \(u(r,z,0)=f(r,z)\) independent of \(\theta \) if the boundary conditions are \(u(r,0,t)=0,u(r,H,t)=0,u(a,z,t)=0\)
Let \(u\left ( r,z,t\right ) =T\left ( t\right ) \Phi \left ( r,z\right ) \). Substituting into the PDE gives\begin {align*} T^{\prime }\left ( t\right ) \Phi \left ( r,z\right ) & =k\left ( T\left ( t\right ) \nabla ^{2}\Phi \right ) \\ \frac {T^{\prime }}{kT} & =\frac {\nabla ^{2}\Phi }{\Phi \left ( r,z\right ) }=-\lambda \end {align*}
Where \(\lambda \) is the separation constant assumed positive. This gives \(T^{\prime }+\lambda kT=0\) with solution \(T\left ( t\right ) =Ce^{-k\lambda t}\). And\begin {equation} \nabla ^{2}\Phi +\lambda \Phi \left ( r,z\right ) =0\tag {1} \end {equation} Where \(\nabla ^{2}\Phi \left ( r,z\right ) =\Phi _{rr}+\frac {1}{r}\Phi _{r}+\Phi _{zz}\) in this case, the Laplacian in cylindrical with no \(\theta \) dependency. Let \(\Phi \left ( r,z\right ) =R\left ( r\right ) Z\left ( z\right ) \). Substituting in (1) gives\begin {align} R^{\prime \prime }Z+\frac {1}{r}R^{\prime }+Z^{\prime \prime }+\lambda RZ & =0\nonumber \\ \frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}+\lambda & =-\frac {Z^{\prime \prime }}{Z}=\upsilon \tag {2} \end {align}
Where \(\upsilon \) is the second separation constant assumed positive. This gives \begin {align*} Z^{\prime \prime }+\upsilon Z & =0\\ Z\left ( z\right ) & =A\cos \left ( \sqrt {\upsilon }z\right ) +B\sin \left ( \sqrt {\upsilon }z\right ) \end {align*}
At \(z=0,Z=0\), hence \(A=0\) and the solution becomes \(Z\left ( z\right ) =B\sin \left ( \sqrt {\upsilon }z\right ) \). At \(z=H,Z=0\), hence for non-trivial solution we want \(\sqrt {\upsilon }H=n\pi ,n=1,2,3,\cdots \). Therefore\[ \upsilon _{n}=\left ( \frac {n\pi }{H}\right ) ^{2}\qquad n=1,2,3,\cdots \] And the corresponding eigenfunctions\[ Z_{n}\left ( z\right ) =\sin \left ( \frac {n\pi }{H}z\right ) \] From (2) the radial equation becomes\begin {align*} \frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}+\lambda & =\left ( \frac {n\pi }{H}\right ) ^{2}\\ R^{\prime \prime }+\frac {1}{r}R^{\prime }+\left ( \lambda -\left ( \frac {n\pi }{H}\right ) ^{2}\right ) R & =0 \end {align*}
This is Bessel ODE. The solution is \[ R_{n}\left ( r\right ) =C_{1}J_{0}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}r\right ) +C_{2}Y_{0}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}r\right ) \] Since \(Y_{0}\) blows up at \(r=0\), it is discarded leaving \(R_{n}\left ( r\right ) =C_{1}J_{0}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}r\right ) \). At \(r=a,R_{n}\left ( a\right ) =0\). For non-trivial solution we want \(J_{0}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}a\right ) =0\). Therefore \(\sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}a\) are the zeros of Bessel function \(J_{0}\left ( x\right ) \). This allows us to determine all possible values of \(\sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}\) (since \(a\) is given constant). Let the \(m^{th}\) zero of \(J_{0}\left ( x\right ) \) be called \(\Lambda _{m}\). Hence \begin {align*} \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}a & =\Lambda _{m}\qquad m=1,2,3,\cdots \\ \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}} & =\frac {\Lambda _{m}}{a}\\ \lambda _{nm} & =\frac {\Lambda _{m}^{2}}{a^{2}}+\left ( \frac {n\pi }{H}\right ) ^{2}\qquad n=1,2,3,\cdots ,m=1,2,3,\cdots \end {align*}
For each \(n\), there are \(m\) values of \(\lambda _{nm}\). The corresponding radial eigenfunction is \[ R_{m}\left ( r\right ) =J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \qquad m=1,2,3,\cdots \] Therefore the complete solution is \begin {align*} u\left ( r,z,t\right ) & =\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\lambda _{nm}t}\sin \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \\ & =\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\left ( \frac {\Lambda _{m}^{2}}{a^{2}}+\left ( \frac {n\pi }{H}\right ) ^{2}\right ) t}\sin \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \end {align*}
What is left is to determine \(A_{nm}\). This is done using initial conditions by using orthogonality. At \(t=0\)\[ f\left ( r,z\right ) =\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}\sin \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \] Multiplying both sides by \(rJ_{0}\left ( \frac {\Lambda _{m^{\prime }}}{a}r\right ) \) and integrating\begin {align*} \int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m^{\prime }}}{a}r\right ) rdr & =\int _{0}^{a}\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}\sin \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) J_{0}\left ( \frac {\Lambda _{m^{\prime }}}{a}r\right ) rdr\\ \int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr & =\sum _{n=1}^{\infty }A_{nm}\sin \left ( \frac {n\pi }{H}z\right ) \int _{0}^{a}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr \end {align*}
Multiplying both sides by \(\sin \left ( \frac {n^{\prime }\pi }{H}z\right ) \) and integrating\begin {align*} \int _{0}^{H}\left ( \int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr\right ) \sin \left ( \frac {n^{\prime }\pi }{H}z\right ) dz & =\int _{0}^{H}\sum _{n=1}^{\infty }A_{nm}\sin \left ( \frac {n\pi }{H}z\right ) \sin \left ( \frac {n^{\prime }\pi }{H}z\right ) \left ( \int _{0}^{a}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr\right ) dz\\ \int _{0}^{H}\int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \sin \left ( \frac {n\pi }{H}z\right ) rdrdz & =A_{nm}\int _{0}^{H}\int _{0}^{a}\sin ^{2}\left ( \frac {n\pi }{H}z\right ) J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz\\ A_{nm} & =\frac {\int _{0}^{H}\int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \sin \left ( \frac {n\pi }{H}z\right ) rdrdz}{\int _{0}^{H}\int _{0}^{a}\sin ^{2}\left ( \frac {n\pi }{H}z\right ) J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz} \end {align*}
Hence the final solution is\[ u\left ( r,z,t\right ) =\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }\frac {\int _{0}^{H}\int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \sin \left ( \frac {n\pi }{H}z\right ) rdrdz}{\int _{0}^{H}\int _{0}^{a}\sin ^{2}\left ( \frac {n\pi }{H}z\right ) J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz}e^{-k\left ( \frac {\Lambda _{m}^{2}}{a^{2}}+\left ( \frac {n\pi }{H}\right ) ^{2}\right ) t}\sin \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \] Where \(\Lambda _{m}\) is the \(m^{th}\) zero of \(J_{0}\left ( x\right ) \). To verify the solution, it is compared to numerical solution, using the following values \(a=1,H=3,k=\frac {1}{100},f\left ( r,z\right ) =\left ( a-r\right ) \sin \left ( \frac {z}{H}\pi \right ) \). The summation was taken up to \(n=10,m=10\). This animation only looks at cross section of the cylinder in the middle. The height indicates the amount of heat. As time passes the cylinder cools down. It runs for 15 seconds. (The animation will only show on the HTML version, not the PDE version of this report).
The following is the source code used to generate the above
____________________________________________________________________________________
Added May 26, 2019.
Problem 7.9.4 (b) from Richard Haberman Applied Partial Differential Equations, 4th edition.
Solve Heat PDE \(u_t = k \nabla ^2 u \) inside cylinder with radius \(a\) and height \(H\) with initial conditions \(u(r,z,0)=f(r,z)\) independent of \(\theta \), subject to boundary conditions \(u_z(r,0,t)=0\), \(u_z(r,H,t)=0\), \(u_r(a,z,t)=0\).
Since it says independent of \(\theta \), will use the PDE as \begin {align*} u_t = k \left (u_{rr} + \frac {1}{r} u_r + u_{zz} \right ) \end {align*}
Instead of the full Laplacian \begin {align*} u_t = k \left (u_{rr} + \frac {1}{r} u_r + \frac {1}{r^2} u_{\theta \theta } + u_{zz} \right ) \end {align*}
Mathematica ✗
ClearAll["Global`*"]; lap = Laplacian[u[r, z, t], {r, theta, z}, "Cylindrical"]; bc = {Derivative[0,1,0][u][r,0, t] == 0, Derivative[0,1,0][u][r, H, t] == 0, Derivative[1,0,0][u][a, z,t] == 0}; ic = u[r,z, 0] == f[r,z]; sol = AbsoluteTiming[TimeConstrained[DSolve[{D[u[r,z,t],t]==k*lap, bc, ic}, u[r, z,t], {r, z,t}, Assumptions -> {a > 0, r < a, H > 0}], 60*10]];
Failed
Maple ✗
restart; lap:=VectorCalculus:-Laplacian(u(r,z,t),'cylindrical'[r,theta,z]); bc:=eval(diff(u(r,z,t),z),z=0)=0,eval(diff(u(r,z,t),z),z=H)=0, eval(diff(u(r,z,t),r),r=a)=0; ic := u(r,z,0) = f(r,z); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([diff(u(r,z,t),t) = k*lap,bc,ic],u(r,z,t)) assuming a>0,r<a,H>0,k>0),output='realtime'));
sol=()
Hand solution
Solve \(u_{t}=k\nabla ^{2}u\) inside cylinder with radius a and height \(H\) with initial conditions \(u(r,z,0)=f(r,z)\) independent of \(\theta \) if the boundary conditions are \(u_{z}(r,0,t)=0,u_{z}(r,H,t)=0,u_{r}(a,z,t)=0\)
Let \(u\left ( r,z,t\right ) =T\left ( t\right ) \Phi \left ( r,z\right ) \). Substituting into the PDE gives\begin {align*} T^{\prime }\left ( t\right ) \Phi \left ( r,z\right ) & =k\left ( T\left ( t\right ) \nabla ^{2}\Phi \right ) \\ \frac {T^{\prime }}{kT} & =\frac {\nabla ^{2}\Phi }{\Phi \left ( r,z\right ) }=-\lambda \end {align*}
Where \(\lambda \) is the separation constant assumed positive. This gives \(T^{\prime }+\lambda kT=0\) with solution \(T\left ( t\right ) =Ce^{-k\lambda t}\). And\begin {equation} \nabla ^{2}\Phi +\lambda \Phi \left ( r,z\right ) =0 \tag {1} \end {equation} Where \(\nabla ^{2}\Phi \left ( r,z\right ) =\Phi _{rr}+\frac {1}{r}\Phi _{r}+\Phi _{zz}\) in this case, the Laplacian in cylindrical with no \(\theta \) dependency. Let \(\Phi \left ( r,z\right ) =R\left ( r\right ) Z\left ( z\right ) \). Substituting in (1) gives\begin {align} R^{\prime \prime }Z+\frac {1}{r}R^{\prime }+Z^{\prime \prime }+\lambda RZ & =0\nonumber \\ \frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}+\lambda & =-\frac {Z^{\prime \prime }}{Z}=\upsilon \tag {2} \end {align}
Where \(\upsilon \) is the second separation. \(\ \)Then \(Z^{\prime \prime }+\upsilon Z=0\)
case \(\upsilon =0\)
This gives\begin {align*} Z\left ( z\right ) & =Az+B\\ Z^{\prime } & =A \end {align*}
At \(z=0,Z^{\prime }=0\), hence \(A=0\). The solution becomes \(Z\left ( z\right ) =B\) and \(Z^{\prime }\left ( z\right ) =0\). Which satisfies the boundary conditions at \(z=H\). Hence \(\upsilon =0\) is eigenvalue with \(Z_{0}=1\) as eigenfunction.
case \(\upsilon >0\) \begin {align*} Z^{\prime \prime }+\upsilon Z & =0\\ Z\left ( z\right ) & =A\cos \left ( \sqrt {\upsilon }z\right ) +B\sin \left ( \sqrt {\upsilon }z\right ) \\ Z^{\prime }\left ( z\right ) & =-A\sqrt {\upsilon }\sin \left ( \sqrt {\upsilon }z\right ) +B\sqrt {\upsilon }\cos \left ( \sqrt {\upsilon }z\right ) \end {align*}
At \(z=0,Z^{\prime }=0\), hence \(B=0\) and the solution becomes \(Z\left ( z\right ) =A\cos \left ( \sqrt {\upsilon }z\right ) \) and \(Z^{\prime }\left ( z\right ) =-A\sqrt {\upsilon }\sin \left ( \sqrt {\upsilon }z\right ) \). At \(z=H,Z_{t}=0\), hence for non-trivial solution we want \(\sqrt {\upsilon }H=n\pi ,n=1,2,3,\cdots \). Therefore\[ \upsilon _{n}=\left ( \frac {n\pi }{H}\right ) ^{2}\qquad n=1,2,3,\cdots \] And the corresponding eigenfunctions\[ Z_{n}\left ( z\right ) =\cos \left ( \frac {n\pi }{H}z\right ) \] From (2) the radial equation becomes\[ \frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}+\lambda =v \]
case \(\upsilon =0\)
The above becomes\[ R_{0}^{\prime \prime }+\frac {1}{r}R_{0}^{\prime }+R_{0}\lambda _{0}=0 \] This is Bessel ODE whose solution is \[ R_{0}\left ( r\right ) =c_{1}J_{0}\left ( \sqrt {\lambda _{0}}r\right ) +c_{2}Y_{0}\left ( \sqrt {\lambda _{0}}r\right ) \] Since \(Y_{0}\) blows at \(r=0\) the solution becomes\begin {align*} R_{0}\left ( r\right ) & =c_{1}J_{0}\left ( \sqrt {\lambda _{0}}r\right ) \\ R_{0}^{\prime }\left ( r\right ) & =c_{1}J_{0}^{\prime }\left ( \sqrt {\lambda _{0}}r\right ) \\ & =-c_{1}\sqrt {\lambda _{0}}J_{1}\left ( \sqrt {\lambda _{0}}r\right ) \end {align*}
At \(r=a\) then \(R_{0}^{\prime }\left ( a\right ) =0\). For nontrivial solution we want \(J_{1}\left ( \sqrt {\lambda }a\right ) =0\). Hence \(\sqrt {\lambda }a\) are zeros of \(J_{1}\left ( x\right ) \). This determines \(\lambda \). Let \(\Lambda _{m}\) be the zeros of \(J_{1}\left ( x\right ) \), therefore\begin {align*} \sqrt {\lambda _{0}}a & =\Lambda _{m}\qquad n=0,m=1,2,3,\cdots \\ \lambda _{0,m} & =\left ( \frac {\Lambda _{m}}{a}\right ) ^{2}\qquad m=1,2,3,\cdots \end {align*}
Hence the radial eigenfunction for \(\upsilon =0\) is given by \[ R_{0,m}\left ( r\right ) =c_{1}J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \qquad m=1,2,3,\cdots \] case \(\upsilon >0\)\begin {align*} \frac {R^{\prime \prime }}{R}+\frac {1}{r}\frac {R^{\prime }}{R}+\lambda & =\left ( \frac {n\pi }{H}\right ) ^{2}\\ R^{\prime \prime }+\frac {1}{r}R^{\prime }+\left ( \lambda -\left ( \frac {n\pi }{H}\right ) ^{2}\right ) R & =0 \end {align*}
This is Bessel ODE. The solution is \[ R_{n}\left ( r\right ) =C_{1}J_{0}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}r\right ) +C_{2}Y_{0}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}r\right ) \] Since \(Y_{0}\) blows up at \(r=0\), it is discarded leaving \(R_{n}\left ( r\right ) =C_{1}J_{0}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}r\right ) \). Hence \(R_{n}^{\prime }\left ( r\right ) =C_{1}J_{0}^{\prime }\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}r\right ) =-C_{1}\sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}J_{1}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}r\right ) \)
At \(r=a,R_{n}^{\prime }\left ( a\right ) =0\). For non-trivial solution we want \(J_{1}\left ( \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}a\right ) =0\). Therefore \(\sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}a\) are the zeros of Bessel function \(J_{1}\left ( x\right ) \). This allows us to determine all possible values of \(\sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}\) (since \(a\) is given constant). Let the \(m^{th}\) zero of \(J_{1}\left ( x\right ) \) be called \(\Lambda _{m}\). Hence \begin {align*} \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}}a & =\Lambda _{m}\qquad m=1,2,3,\cdots \\ \sqrt {\lambda -\left ( \frac {n\pi }{H}\right ) ^{2}} & =\frac {\Lambda _{m}}{a}\\ \lambda _{nm} & =\frac {\Lambda _{m}^{2}}{a^{2}}+\left ( \frac {n\pi }{H}\right ) ^{2}\qquad n=1,2,3,\cdots ,m=1,2,3,\cdots \end {align*}
For each \(n\), there are \(m\) values of \(\lambda _{nm}\). The corresponding radial eigenfunction is \[ R_{m}\left ( r\right ) =J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \qquad m=1,2,3,\cdots \] Now we look at the time solution again. \(T\left ( t\right ) =Ce^{-k\lambda t}\). For \(n=0\) this gives \(T_{0,m}\left ( t\right ) =C_{0,m}e^{-k\lambda _{0,m}t}=C_{0,m}e^{-k\left ( \frac {\Lambda _{m}}{a}\right ) ^{2}t}\). For \(n>0\) the solution becomes \[ T_{n,m}\left ( t\right ) =C_{n,m}e^{-k\lambda _{n,m}t}=Ce^{-k\left ( \frac {\Lambda _{m}^{2}}{a^{2}}+\left ( \frac {n\pi }{H}\right ) ^{2}\right ) t}\] Therefore the complete solution is \begin {align*} u\left ( r,z,t\right ) & =\sum _{m=1}^{\infty }B_{m}R_{0,m}\left ( r\right ) T_{0,m}\left ( t\right ) +\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\lambda _{nm}t}\cos \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \\ & =\sum _{m=1}^{\infty }B_{m}J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) e^{-k\left ( \frac {\Lambda _{m}}{a}\right ) ^{2}t}+\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\left ( \frac {\Lambda _{m}^{2}}{a^{2}}+\left ( \frac {n\pi }{H}\right ) ^{2}\right ) t}\cos \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \end {align*}
What is left is to determine \(A_{nm}\) and \(B_{m}\). This is done using initial conditions by using orthogonality. At \(t=0\)\[ f\left ( r,z\right ) =\sum _{m=1}^{\infty }B_{m}J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) +\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}\cos \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \] For \(n=0\)\begin {align*} f\left ( r,z\right ) & =\sum _{m=1}^{\infty }B_{m}J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \\ \int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m^{\prime }}}{a}r\right ) rdr & =\int _{0}^{a}\sum _{m=1}^{\infty }B_{m}J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) J_{0}\left ( \frac {\Lambda _{m^{\prime }}}{a}r\right ) rdr\\ \int _{0}^{a}f\left ( r,0\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr & =\int _{0}^{a}B_{m}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr\\ B_{m} & =\frac {\int _{0}^{a}f\left ( r,0\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr}{\int _{0}^{a}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr}\qquad m=1,2,3,\cdots \end {align*}
Integrating again both sides over \(z\) gives (we can think of \(\cos \left ( \frac {n\pi }{H}z\right ) =1\) with \(n=0\) as the eigenfunction in this case).
\begin {align*} \int _{0}^{H}\int _{0}^{a}f\left ( r,0\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz & =B_{m}\int _{0}^{H}\int _{0}^{a}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz\\ B_{m} & =\frac {\int _{0}^{H}\int _{0}^{a}f\left ( r,0\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz}{\int _{0}^{H}\int _{0}^{a}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz} \end {align*}
For \(n>0\)\[ f\left ( r,z\right ) =\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}\cos \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \] Multiplying both sides by \(rJ_{0}\left ( \frac {\Lambda _{m^{\prime }}}{a}r\right ) \) and integrating\begin {align*} \int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m^{\prime }}}{a}r\right ) rdr & =\int _{0}^{a}\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}\cos \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) J_{0}\left ( \frac {\Lambda _{m^{\prime }}}{a}r\right ) rdr\\ \int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr & =\sum _{n=1}^{\infty }A_{nm}\cos \left ( \frac {n\pi }{H}z\right ) \int _{0}^{a}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr \end {align*}
Multiplying both sides by \(\cos \left ( \frac {n^{\prime }\pi }{H}z\right ) \) and integrating\begin {align*} \int _{0}^{H}\left ( \int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr\right ) \cos \left ( \frac {n^{\prime }\pi }{H}z\right ) dz & =\int _{0}^{H}\sum _{n=1}^{\infty }A_{nm}\sin \left ( \frac {n\pi }{H}z\right ) \cos \left ( \frac {n^{\prime }\pi }{H}z\right ) \left ( \int _{0}^{a}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdr\right ) dz\\ \int _{0}^{H}\int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \cos \left ( \frac {n\pi }{H}z\right ) rdrdz & =A_{nm}\int _{0}^{H}\int _{0}^{a}\cos ^{2}\left ( \frac {n\pi }{H}z\right ) J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz\\ A_{nm} & =\frac {\int _{0}^{H}\int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \cos \left ( \frac {n\pi }{H}z\right ) rdrdz}{\int _{0}^{H}\int _{0}^{a}\cos ^{2}\left ( \frac {n\pi }{H}z\right ) J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz} \end {align*}
Hence the final solution is\[ u\left ( r,z,t\right ) =\sum _{m=1}^{\infty }B_{m}J_{0}\left ( \left ( \frac {\Lambda _{m}}{a}\right ) ^{2}r\right ) e^{-k\left ( \frac {\Lambda _{m}}{a}\right ) ^{2}t}+\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\left ( \frac {\Lambda _{m}^{2}}{a^{2}}+\left ( \frac {n\pi }{H}\right ) ^{2}\right ) t}\cos \left ( \frac {n\pi }{H}z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \]
With
\begin {align*} A_{nm} & =\frac {\int _{0}^{H}\int _{0}^{a}f\left ( r,z\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) \cos \left ( \frac {n\pi }{H}z\right ) rdrdz}{\int _{0}^{H}\int _{0}^{a}\cos ^{2}\left ( \frac {n\pi }{H}z\right ) J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz}\\ B_{m} & =\frac {\int _{0}^{H}\int _{0}^{a}f\left ( r,0\right ) J_{0}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz}{\int _{0}^{H}\int _{0}^{a}J_{0}^{2}\left ( \frac {\Lambda _{m}}{a}r\right ) rdrdz} \end {align*}
Where \(\Lambda _{m}\) is the \(m^{th}\) zero of \(J_{1}\left ( x\right ) \).
To verify the solution, it is compared to numerical solution, using the following values \(a=1,H=3,k=\frac {1}{100},f\left ( r,z\right ) =\left ( a-r\right ) \sin \left ( \frac {z}{H}\pi \right ) \). The summation was taken up to \(n=10,m=10\). This animation only looks at cross section of the cylinder in the middle. The height indicates the amount of heat. As time passes the initial temperature averages inside (cylinder is insulated). It runs for 15 seconds. (The animation will only show on the HTML version, not the PDE version of this report).
The following is the source code used to generate the above
____________________________________________________________________________________
Added May 26, 2019.
Problem 7.9.4 (c) from Richard Haberman Applied Partial Differential Equations, 4th edition.
Solve Heat PDE \(u_t = k \nabla ^2 u \) inside cylinder with radius \(a\) and height \(H\) with initial conditions \(u(r,z,0)=f(r,z)\) subject to boundary conditions \(u(r,0,t)=0\), \(u(r,H,t)=0\), \(u_r(a,z,t)=0\).
Since it says independent of \(\theta \), will use the PDE as \begin {align*} u_t = k \left (u_{rr} + \frac {1}{r} u_r + u_{zz} \right ) \end {align*}
Instead of full Laplacian \begin {align*} u_t = k \left (u_{rr} + \frac {1}{r} u_r + \frac {1}{r^2} u_{\theta \theta } + u_{zz} \right ) \end {align*}
Mathematica ✗
ClearAll["Global`*"]; lap = Laplacian[u[r, z, t], {r, theta, z}, "Cylindrical"]; bc = {u[r, 0, t] == 0, u[r, H, t] == 0, Derivative[1,0,0][u][a, z,t] == 0}; ic = u[r, z, 0] == f[r,z]; sol = AbsoluteTiming[TimeConstrained[DSolve[{D[u[r,z,t],t]==k*lap, bc, ic}, u[r, z,t], {r, z,t}, Assumptions -> {a > 0, r < a, H > 0}], 60*10]];
Failed
Maple ✗
restart; lap:=VectorCalculus:-Laplacian(u(r,z,t),'cylindrical'[r,theta,z]); bc := u(r,0,t)=0,u(r,H,t)=0, eval(diff(u(r,z,t),r),r=a)=0; ic := u(r,z,0) = f(r,z); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([diff(u(r,z,t),t) = k*lap,bc,ic],u(r,z,t)) assuming a>0,r<a,H>0,k>0),output='realtime'));
sol=()
____________________________________________________________________________________
Added May 26, 2019.
Problem 7.9.3 (a) from Richard Haberman Applied Partial Differential Equations, 4th edition.
Solve Heat PDE \(u_t = k \nabla ^2 u \) inside quarter circular cylinder \(0<\theta <\frac {\pi }{2}\) with radius \(a\) and height \(H\) with initial conditions \(u(r,\theta ,z,0)=f(r\theta ,z)\) subject to boundary conditions \(u(r,\theta ,0,t)=0\), \(u(r,\theta ,H,t)=0\), \(u(r,0,z,t)=0\), \(u(r,\frac {\pi }{2},z,t)=0\), \(u(a,\theta ,z,t)=0\).
\begin {align*} u_t = k \left (u_{rr} + \frac {1}{r} u_r + \frac {1}{r^2} u_{\theta \theta } + u_{zz} \right ) \end {align*}
Mathematica ✗
ClearAll["Global`*"]; lap = Laplacian[u[r, theta, z, t], {r, theta, z}, "Cylindrical"]; bc = {u[r, theta, 0, t] == 0, u[r, theta, H, t] == 0, u[r, 0, z,t] == 0,u[r,Pi/2,z,t]==0,u[a,Pi,z,t]==0}; ic = u[r, theta, z, 0] == f[r,theta,z]; sol = AbsoluteTiming[TimeConstrained[DSolve[{D[u[r,theta,z,t],t]==k*lap, bc, ic}, u[r, theta, z,t], {r, theta, z,t}, Assumptions -> {a > 0, r < a, H > 0,theta>0,theta<Pi/2}], 60*10]];
Failed
Maple ✗
restart; lap:=VectorCalculus:-Laplacian(u(r,theta,z,t),'cylindrical'[r,theta,z]); bc := u(r,theta,0,t)=0,u(r,theta,H,t)=0, u(r,0,z,t)=0, u(r,Pi/2,z,t)=0,u(a,Pi,z,t)=0; ic := u(r,theta,z,0) = f(r,theta,z); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([diff(u(r,theta,z,t),t) = k*lap,bc,ic],u(r,theta,z,t)) assuming a>0,r<a,H>0,theta>0,theta<Pi/2,k>0),output='realtime'));
sol=()
____________________________________________________________________________________
Added May 26, 2019.
Problem 7.9.3 (b) from Richard Haberman Applied Partial Differential Equations, 4th edition.
Solve Heat PDE \(u_t = k \nabla ^2 u \) inside quarter circular cylinder \(0<\theta <\frac {\pi }{2}\) with radius \(a\) and height \(H\) with initial conditions \(u(r,\theta ,z,0)=f(r\theta ,z)\) subject to boundary conditions \(u_z(r,\theta ,0,t)=0\), \(u_z(r,\theta ,H,t)=0\), \(u_\theta (r,0,z,t)=0\), \(u_\theta (r,\frac {\pi }{2},z,t)=0\), \(u_r(a,\theta ,z,t)=0\).
\begin {align*} u_t = k\left ( u_{rr} + \frac {1}{r} u_r + \frac {1}{r^2} u_{\theta \theta } + u_{zz} \right ) \end {align*}
Mathematica ✗
ClearAll["Global`*"]; lap = Laplacian[u[r, theta, z, t], {r, theta, z}, "Cylindrical"]; bc = {Derivative[0,0,1,0][u][r, theta, 0, t] == 0, Derivative[0,0,1,0][u][r, theta, H, t] == 0, Derivative[0,1,0,0][u][r, 0, z,t] == 0,Derivative[0,1,0,0][u][r,Pi/2,z,t]==0,Derivative[1,0,0,0][u][a,theta,z,t]==0}; ic = u[r, theta, z, 0] == f[r,theta,z]; sol = AbsoluteTiming[TimeConstrained[DSolve[{D[u[r,theta,z,t],t]==k*lap, bc}, u[r, theta, z,t], {r, theta, z,t}, Assumptions -> {a > 0, r < a, H > 0,theta>0,theta<Pi/2}], 60*10]];
Failed
Maple ✗
restart; lap:=VectorCalculus:-Laplacian(u(r,theta,z,t),'cylindrical'[r,theta,z]); bc:=eval(diff(u(r,theta,z,t),z),z=0)=0,eval(diff(u(r,theta,z,t),z),z=H)=0,eval(diff(u(r,theta,z,t),theta),theta=0)=0,eval(diff(u(r,theta,z,t),theta),theta=Pi/2)=0,eval(diff(u(r,theta,z,t),r),r=a)=0; ic := u(r,theta,z,0) = f(r,theta,z); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([diff(u(r,theta,z,t),t) = k*lap,bc,ic],u(r,theta,z,t)) assuming a>0,r<a,H>0,theta>0,theta<Pi/2,k>0),output='realtime'));
sol=()
____________________________________________________________________________________
Added May 26, 2019.
Problem 7.9.3 (c) from Richard Haberman Applied Partial Differential Equations, 4th edition.
Solve Heat PDE \(u_t = k \nabla ^2 u \) inside quarter circular cylinder \(0<\theta <\frac {\pi }{2}\) with radius \(a\) and height \(H\) with initial conditions \(u(r,\theta ,z,0)=f(r\theta ,z)\) subject to boundary conditions \(u(r,\theta ,0,t)=0\), \(u(r,\theta ,H,t)=0\), \(u_\theta (r,0,z,t)=0\), \(u(r,\frac {\pi }{2},z,t)=0\), \(u_r(a,\theta ,z,t)=0\).
\begin {align*} u_t = k \left ( u_{rr} + \frac {1}{r} u_r + \frac {1}{r^2} u_{\theta \theta } + u_{zz} \right ) \end {align*}
Mathematica ✗
ClearAll["Global`*"]; lap = Laplacian[u[r, theta, z, t], {r, theta, z}, "Cylindrical"]; bc = {u[r, theta, 0, t] == 0, u[r, theta, H, t] == 0, u[r, 0, z,t] == 0,u[r,Pi/2,z,t]==0,u[a,Pi,z,t]==0}; ic = u[r, theta, z, 0] == f[r,theta,z]; sol = AbsoluteTiming[TimeConstrained[DSolve[{D[u[r,theta,z,t],t]==k*lap, bc}, u[r, theta, z,t], {r, theta, z,t}, Assumptions -> {a > 0, r < a, H > 0,theta>0,theta<Pi/2}], 60*10]];
Failed
Maple ✗
restart; lap:=VectorCalculus:-Laplacian(u(r,theta,z,t),'cylindrical'[r,theta,z]); bc := u(r,theta,0,t)=0,u(r,theta,H,t)=0, eval(diff(u(r,theta,z,t),theta),theta=0)=0, u(r,Pi/2,z,t)=0,eval(diff(u(r,theta,z,t),r),r=a)=0; ic := u(r,theta,z,0) = f(r,theta,z); cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([diff(u(r,theta,z,t),t) = k*lap,bc,ic],u(r,theta,z,t)) assuming a>0,r<a,H>0,theta>0,theta<Pi/2,k>0),output='realtime'));
sol=()