1.4.1 Heat PDE in 2D inside a Disk, zero temperature at boundary (circumference) of disk.

   1.4.1.1 Animations

Solve for \(u\left ( r,\theta ,t\right ) \), the heat PDE \[ u_{t}=k\left ( u_{rr}+\frac{1}{r}u_{r}+\frac{1}{r^{2}}u_{\theta \theta }\right ) \] where \(0<r<a\) is the radius of the disk and \(-\pi <\theta <\pi \) and \(t>0\). Initial conditions are\[ u\left ( r,\theta ,0\right ) =f\left ( r,\theta \right ) \] Boundary conditions are\begin{align*} u\left ( a,\theta ,t\right ) & =0\\ u\left ( 0,\theta ,t\right ) & <\infty \end{align*}

Boundary conditions for \(\theta \) is that it is periodic\begin{align*} u\left ( r,-\pi ,t\right ) & =u\left ( r,\pi ,t\right ) \\ u_{\theta }\left ( r,-\pi ,t\right ) & =u_{\theta }\left ( r,\pi ,t\right ) \end{align*}

Analytical solution

Using separation of variables, let \(u=R\left ( r\right ) \Theta \left ( \theta \right ) T\left ( t\right ) \). The PDE becomes\[ \frac{1}{k}T^{\prime }R\Theta =R^{\prime \prime }\Theta T+\frac{1}{r}R^{\prime }\Theta T+\frac{1}{r^{2}}\Theta ^{\prime \prime }RT \] Dividing by \(R\Theta T\) gives\begin{equation} \frac{1}{k}\frac{T^{\prime }}{T}=\frac{R^{\prime \prime }}{R}+\frac{1}{r}\frac{R^{\prime }}{R}+\frac{1}{r^{2}}\frac{\Theta ^{\prime \prime }}{\Theta }=-\lambda ^{2} \tag{1} \end{equation} Where \(\lambda \) is the first separation constant. This gives the time ODE\begin{equation} T^{\prime }+\lambda ^{2}kT=0 \tag{2} \end{equation} And from (1)\[ \frac{R^{\prime \prime }}{R}+\frac{1}{r}\frac{R^{\prime }}{R}+\frac{1}{r^{2}}\frac{\Theta ^{\prime \prime }}{\Theta }=-\lambda ^{2}\] Multiplying by \(r^{2}\) and rearranging\[ r^{2}\frac{R^{\prime \prime }}{R}+r\frac{R^{\prime }}{R}+r^{2}\lambda ^{2}=-\frac{\Theta ^{\prime \prime }}{\Theta }=\mu ^{2}\] Where \(\mu \) is the second separation constant. This gives the \(R\) ODE as\begin{equation} r^{2}R^{\prime \prime }+rR^{\prime }+\left ( r^{2}\lambda ^{2}-\mu ^{2}\right ) R=0 \tag{3} \end{equation} And the \(\Theta \) ODE as\begin{equation} \Theta ^{\prime \prime }+\mu ^{2}\Theta =0 \tag{4} \end{equation} The eigenvalues for (4) determine the Bessel equation (3) order. Hence starting with (4), the ODE and its boundary conditions becomes\begin{align*} \Theta ^{\prime \prime }+\mu ^{2}\Theta & =0\\ \Theta \left ( -\pi \right ) & =\Theta \left ( \pi \right ) \\ \Theta ^{\prime }\left ( -\pi \right ) & =\Theta ^{\prime }\left ( \pi \right ) \end{align*}

\(\mu =0\) leads to solution \begin{align*} \Theta & =c_{1}\theta +c_{2}\\ \Theta ^{\prime } & =c_{1} \end{align*}

First BC gives\begin{align*} -c_{1}\pi +c_{2} & =c_{1}\pi +c_{2}\\ c_{1} & =0 \end{align*}

And since second BC \(\Theta ^{\prime }=c_{1}\), this implies \(\Theta =\) constant is a solution. So \(\mu =0\) is an eigenvalue, with \(\Theta _{0}\left ( \theta \right ) =\) constant being the eigenfunction.\(\mu <0\) is not possible eigenvalue, since the solution is made up of real exponentials When \(\mu >0\) the solution is\[ \Theta =A\cos \mu \theta +B\sin \mu \theta \] To satisfy the periodic boundary conditions, \(\mu \) must be an integer, and since \(\mu >0\), then \(\mu =n\) for \(n=1,2,3,\cdots \). Therefore\begin{align} \Theta _{0}\left ( \theta \right ) & =\text{constant}\qquad n=0\tag{5A}\\ \Theta _{n}\left ( \theta \right ) & =A_{n}\cos n\theta +B_{n}\sin n\theta \qquad n=1,2,3,\cdots \tag{5B} \end{align}

Since \(\mu \) is now found, then the Bessel ODE (3) is solved.  \begin{align} r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) +\left ( r^{2}\lambda ^{2}-n^{2}\right ) R\left ( r\right ) & =0\qquad n=0,1,2,3,\cdots \tag{5C}\\ R\left ( a\right ) & =0\nonumber \\ R\left ( 0\right ) & <\infty \nonumber \end{align}

Note that \(\lambda =0\) is not possible. For example, if \(\lambda =0\) then (5C) becomes Euler ODE\[ r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) +n^{2}R\left ( r\right ) =0\qquad n=0,1,2,3,\cdots \] For the case \(n=0\), the ODE becomes \(r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) =0\) whose solution is \(R\left ( r\right ) =c_{1}+c_{2}\ln \left ( r\right ) \). Since bounded at \(r=0\), the solution becomes \(R\left ( r\right ) =c_{1}\) and since \(R\left ( a\right ) =0\) then \(c_{1}=0\), leading to trivial solution. And for the case of \(n>0\), the ODE becomes \(r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) +n^{2}R\left ( r\right ) =0\) whose solution is \(R\left ( r\right ) =c_{1}r^{n}+c_{2}\frac{1}{r^{n}}\). Since bounded at \(r=0\), then \(c_{2}=0\) and the solution becomes \(R\left ( r\right ) =c_{1}r^{n}\). Using BC \(R\left ( a\right ) =0\), gives \(0=c_{1}a^{n}\). Hence \(c_{1}=0\) leading to trivial solution. This shows that \(\lambda =0\) is not possible. Back to solving (5C) now that it was shown that \(\lambda \) is not zero, the first step is to convert this to Bessel ODE in the classical form to obtain a solution. Let \(t=\lambda r\), then \(R^{\prime }\left ( r\right ) =R^{\prime }\left ( t\right ) \lambda \) and \(R^{\prime \prime }\left ( r\right ) =R^{\prime \prime }\left ( t\right ) \lambda ^{2}\). The ODE becomes\begin{align*} \frac{t^{2}}{\lambda ^{2}}\lambda ^{2}R^{\prime \prime }\left ( t\right ) +\frac{t}{\lambda }\lambda R^{\prime }\left ( t\right ) +\left ( \frac{t^{2}}{\lambda ^{2}}\lambda ^{2}-n^{2}\right ) R\left ( t\right ) & =0\\ t^{2}R^{\prime \prime }\left ( t\right ) +tR^{\prime }\left ( t\right ) +\left ( t^{2}-n^{2}\right ) R\left ( t\right ) & =0 \end{align*}

This is now in standard Bessel ODE form. This is of order \(n\). The solution is given by\[ R_{n}\left ( t\right ) =C_{n}J_{n}\left ( t\right ) +D_{n}Y_{n}\left ( t\right ) \] Where \(J_{n}\left ( t\right ) \) is the Bessel function of order \(n\) and \(Y_{n}\left ( t\right ) \) is the Bessel function of second kind of order \(n\). Converting back to \(r\) gives\[ R_{n}\left ( \lambda r\right ) =C_{n}J_{n}\left ( \lambda r\right ) +D_{n}Y_{n}\left ( \lambda r\right ) \] Since the solution is bounded at \(r=0\) and since \(Y_{n}\left ( 0\right ) \) blows up, then \(D_{n}=0\). The solution simplifies to\[ R_{n}\left ( \lambda r\right ) =C_{n}J_{n}\left ( \lambda r\right ) \] Applying the second boundary conditions, when \(r=a\) then\[ 0=C_{n}J_{n}\left ( \lambda a\right ) \] For non-trivial solution \(J_{n}\left ( \lambda a\right ) =0\). Hence \(\lambda a\) are the zeros of \(J_{n}\left ( z\right ) \). Let the positive zeros of \(J_{n}\left ( z\right ) \) be \(j_{m}\). For \(m=1,2,3,\cdots \). Therefore \[ \lambda _{nm}=\frac{j_{m}}{a}\qquad n=0,1,2,\cdots ,m=1,2,3,\cdots \] is the \(m^{th}\) eigenvalue for the \(n^{th}\) Bessel function \(J_{n}\). So there are two indices to handle in these problems. The order of the Bessel function is determined from the \(\Theta _{n}\left ( \theta \right ) \) eigenvalues, and then once the order \(n\) is fixed, the second eigenvalue \(\lambda _{nm}\) is determined from the zeros of the Bessel function. Hence the \(R\) solution is\[ R_{n}\left ( r\right ) =C_{n}J_{n}\left ( \lambda _{nm}r\right ) \qquad n=0,1,2,3,\cdots ,m=1,2,3,\cdots \] Now that \(\lambda _{nm}=\frac{j_{m}}{a}\) is known, then the time ODE (2) is solved\begin{align*} T_{nm}^{\prime }+\lambda _{nm}^{2}kT_{nm} & =0\\ T_{nm} & =e^{-k\lambda _{nm}^{2}t}\qquad n=0,1,2,3,\cdots ,m=1,2,3,\cdots \end{align*}

Hence the fundamental solution is\[ u_{nm}\left ( r,\theta ,t\right ) =\Theta _{n}\left ( \theta \right ) T_{nm}\left ( t\right ) R_{nm}\left ( r\right ) \] And the complete solution is\begin{align*} u\left ( r,\theta ,t\right ) & =\sum _{n=0}^{\infty }\Theta _{n}\left ( \theta \right ) \sum _{m=1}^{\infty }T_{nm}\left ( t\right ) R_{nm}\left ( r\right ) \\ & =\sum _{n=0}^{\infty }\left ( A_{n}\cos n\theta +B_{n}\sin n\theta \right ) \left ( \sum _{m=1}^{\infty }e^{-k\lambda _{nm}^{2}t}C_{nm}J_{n}\left ( \lambda _{nm}r\right ) \right ) \end{align*}

The above can now be written as\[ u\left ( r,\theta ,t\right ) =\sum _{n=0}^{\infty }\cos n\theta \left ( \sum _{m=1}^{\infty }e^{-k\lambda _{nm}^{2}t}A_{n}C_{nm}J_{n}\left ( \lambda _{nm}r\right ) \right ) +\sum _{n=0}^{\infty }\sin n\theta \left ( \sum _{m=1}^{\infty }e^{-k\lambda _{nm}^{2}t}B_{n}C_{nm}J_{n}\left ( \lambda _{nm}r\right ) \right ) \] Let \(A_{n}C_{nm}\) be combined into one constant \(A_{nm}\) and \(B_{n}C_{nm}\) combined into one constant \(B_{nm}.\) The above becomes\[ u\left ( r,\theta ,t\right ) =\sum _{n=0}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \cos n\theta +\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }B_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \sin n\theta \] The constants \(A_{nm},B_{nm}\) are found from initial conditions. At \(t=0\), \(u\left ( r,\theta ,0\right ) =f\left ( r,\theta \right ) \). Hence the above becomes\[ f\left ( r,\theta \right ) =\sum _{n=0}^{\infty }\sum _{m=1}^{\infty }A_{nm}J_{n}\left ( \lambda _{nm}r\right ) \cos n\theta +\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }B_{nm}J_{n}\left ( \lambda _{nm}r\right ) \sin n\theta \] Multiplying both sides by \(\cos n\theta \) and integrating simplifies to\[ \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \cos n\theta d\theta =\int _{-\pi }^{\pi }\cos ^{2}n\theta d\theta \left ( \sum _{m=1}^{\infty }A_{nm}J_{n}\left ( \lambda _{nm}r\right ) \right ) \] For \(n=0\) case\[ \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \cos n\theta d\theta =2\pi \sum _{m=1}^{\infty }A_{0m}J_{0}\left ( \lambda _{0m}r\right ) \] For \(n>0\) case\[ \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \cos n\theta d\theta =\pi \sum _{m=1}^{\infty }A_{nm}J_{n}\left ( \lambda _{nm}r\right ) \] Orthogonality is applied again. Multiplying both side by \(rJ_{n}\left ( \lambda _{nm}r\right ) \) and integrating (weight is \(r\)).

For \(n=0\) case

\begin{align} \int _{0}^{a}\left ( \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \cos n\theta d\theta \right ) rJ_{n}\left ( \lambda _{nm}r\right ) dr & =2\pi A_{nm}\int _{0}^{a}rJ_{n}^{2}\left ( \lambda _{nm}r\right ) dr\nonumber \\ A_{0m} & =\frac{\int _{0}^{a}\left ( \int _{-\pi }^{\pi }f\left ( r,\theta \right ) d\theta \right ) rJ_{0}\left ( \lambda _{0m}r\right ) dr}{2\pi \int _{0}^{a}rJ_{0}^{2}\left ( \lambda _{0m}r\right ) dr}\qquad n=0,m=1,2,3,\cdots \tag{6AA} \end{align}

For \(n>0\) case\begin{align} \int _{0}^{a}\left ( \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \cos n\theta d\theta \right ) rJ_{n}\left ( \lambda _{nm}r\right ) dr & =\pi A_{nm}\int _{0}^{a}rJ_{n}^{2}\left ( \lambda _{nm}r\right ) dr\nonumber \\ A_{nm} & =\frac{\int _{0}^{a}\left ( \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \cos \left ( n\theta \right ) d\theta \right ) rJ_{n}\left ( \lambda _{nm}r\right ) dr}{\pi \int _{0}^{a}rJ_{n}^{2}\left ( \lambda _{nm}r\right ) dr}\qquad n=1,2,\cdots ,m=1,2,3,\cdots \tag{6AB} \end{align}

Similarly, \begin{equation} B_{nm}=\frac{\int _{0}^{a}\left ( \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \sin \left ( n\theta \right ) d\theta \right ) rJ_{n}\left ( \lambda _{nm}r\right ) dr}{\pi \int _{0}^{a}rJ_{n}^{2}\left ( \lambda _{nm}r\right ) dr}\qquad n=1,2,\cdots ,m=1,2,3,\cdots \tag{6B} \end{equation} This complete the solution. Hence the solution is\[ u\left ( r,\theta ,t\right ) =\sum _{n=0}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \cos n\theta +\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \sin n\theta \] With \(A_{nm}\) given by (6AA,6AB) and \(B_{nm}\) by (6B) and \(\lambda _{nm}=\frac{j_{m}}{a}\) where \(j_{m}\) is the \(m^{th}\) zero of the \(J_{n}\left ( z\right ) \).

1.4.1.1 Animations

Example 1

This animation uses \(a=1,f\left ( r,\theta \right ) =\left ( a-r^{2}\right ) ,k=1\). It used \(n=0\cdots 6\) terms and \(m=1\cdots 6\) terms for the double sum. Animation time is 0.4 seconds.

Example 2

This animation uses \(a=1,f\left ( r,\theta \right ) =\left ( r-r^{2}\right ) \sin \theta ,k=1\). Same number of terms above the above animation. Animation time is 0.15 seconds.

Example 3

This animation uses \(a=1,f\left ( r,\theta \right ) =\left ( r-r^{2}\right ) \sin \theta \cos \theta ,k=1\). Same number of terms above the above animation. Animation time is 0.15 seconds.