Heat equation/Solution to the 2-D Heat Equation in Cylindrical Coordinates

Definition
We start by changing the Laplacian operator in the 2-D heat equation from rectangular to cylindrical coordinates by the following definition: $$D:=(0,a) \times (0,b)~.$$ By changing the coordinate system, we arrive at the following nonhomogeneous PDE for the heat equation: $$u_t=k \left [ \frac{1}{r} \left ( u_r + ru_{rr} \right ) + \frac{1}{r^2} u_{\theta \theta} \right ] + h(r,\theta,t), \text{ where } (r,\theta) \in D, t \in (0,\infty)~.$$ We choose for the example the Robin boundary conditions and initial conditions as follows: $$ \begin{cases} \left\vert u(0,\theta,t) \right\vert < \infty \\ \alpha_1u(a,\theta,t)+\beta_1u_r(a,\theta,t)=0 \\ \alpha_2u(r,0,t)-\beta_2u_\theta(r,0,t)=0 \\ \alpha_3u(r,b,t)+\beta_3u_\theta(r,b,t)=0 \\ u(r,\theta,0)=f(r,\theta) \end{cases} $$

Separate Variables
$$u(r,\theta,t)=R(r)\Theta(\theta)T(t)$$ $$\Rightarrow R\Theta T'=k \left [ \frac{1}{r}(R'\Theta T+rR\Theta T) + \frac{1}{r^2} R\ThetaT \right ]$$ $$\Rightarrow \frac{T'}{kT} = \frac{1}{r} \left ( \frac{R'}{R}+ r \frac{R}{R} \right ) + \frac{1}{r^2} \frac{\Theta}{\Theta}$$ This means that a separation constant can be found that both sides will equal. Let's define it to be $$\lambda^2~.$$ This yields: $${\color{Blue}T'+k\lambda^2 T=0}$$ and multiplying the other side by $$r^2$$ yields: $$r \left ( \frac{R'}{R}+r\frac{R}{R} \right ) + \lambda^2r^2=-\frac{\Theta}{\Theta}$$ After defining another separation constant $$\mu^2$$, it yields: $${\color{Blue}\Theta''+\mu^2\Theta=0}$$ Multiplying the other side by R yields: $${\color{Blue}r^2R''+rR'+\left ( \lambda^2r^2-\mu^2 \right ) R = 0}$$ We now have separate differential equations for each variable.

Translate Boundary Conditions
$$\left\vert R(0)\Theta(\theta)T(t) \right\vert < \infty \Rightarrow \left\vert R(0) \right\vert < \infty$$ $$\alpha_1R(a)\Theta(\theta)T(t)+\beta_1R'(a)\Theta(\theta)T(t)=0 \Rightarrow \alpha_1R(a)+\beta_1R'(a)=0$$ $$\alpha_2R(r)\Theta(0)T(t)+\beta_2R(r)\Theta'(0)T(t)=0 \Rightarrow \alpha_2\Theta(0)-\beta_2\Theta'(0)=0$$ $$\alpha_3R(r)\Theta(b)T(t)+\beta_2R(r)\Theta'(b)T(t)=0 \Rightarrow \alpha_3\Theta(b)+\beta_3\Theta'(b)=0$$

Solve SLPs
$$ \left. \begin{align} \Theta''+\mu^2\Theta=0 \\ \alpha_2\Theta(0)-\beta_2\Theta'(0)=0\\ \alpha_3\Theta(b)+\beta_3\Theta'(b)=0 \end{align} \right \} \Rightarrow \begin{align} & \text{Eigenvalues } \mu_m \text{ are solutions to } (\alpha_2\alpha_3-\beta_2\beta_3\mu^2)\sin (\mu B)+(\alpha_2\beta_3+\alpha_3\beta_2)\mu\cos(\mu B)=0\\ & \Theta_m(\theta)=\beta_2\mu_m\cos(\mu_m \theta )+\alpha_2\sin (\mu_m \theta ), m=0,1,2,\cdots \end{align} $$ The SLP for $$R$$ is a singular Bessel type, whose eigenvalues $$\lambda_{mn}$$ depends on $$\mu_m$$ and are non-negative solutions to the following equation: $$\left ( \alpha_1 a + \beta_1 \mu_m \right ) J_{ \mu_m }(\lambda a) - \beta_1 a \lambda J_{\mu_m+1}(\lambda a) = 0$$ and the eigenfunction is: $$R_{mn}(r)=J_{\mu_m}(\lambda_{mn} r), \quad m,n=0,1,2,\cdots$$ where $$J_\nu(x)$$ is the Bessel function of the first kind of order $$\nu$$.

Solve Time Equation
$$T'+k\lambda_{mn}^2 T=0$$ $$\Rightarrow T_{mn}(t)=C_{mn}e^{-k\lambda_{mn}^2t}$$

Step 2: Satisfy Initial Condition
Let's define the solution as an infinite sum: $$u(r,\theta,t):=\sum_{m,n=0}^\infty T_{mn}(t)R_{mn}(r)\Theta_m(\theta)~.$$ With the initial condition: $$ \begin{align} f(r,\theta)&=u(r,\theta,0) \\ &=\sum_{m,n=0}^\infty T_{mn}(0)R_{mn}(r)\Theta_m(\theta) \\ &=\sum_{m,n=0}^\infty C_{mn}R_{mn}(r)\Theta_m(\theta) \end{align} $$ where $$C_{mn}=\frac{\int\limits_0^b \int\limits_0^a f(r,\theta)R_{mn}(r)\Theta_m(\theta) r dr d\theta }{\int\limits_0^a R_{mn}^2(r) r dr \int\limits_0^b \Theta_m^2(\theta) d\theta}~.$$ The weight function in the inner product $$w(r)=r$$ in integrals involving the Bessel functions. The Bessel functions $$R_m$$ are orthogonal relative to the "weighted" scalar product $$< f,g >_w:=\int\limits_0^a f(r)g(r)r dr~.$$

Step 3: Solve Non-homogeneous Equation
Solving the non-homogeneous equation involves defining the following functions: $$u(r,\theta,t):=\sum_{m,n=0}^\infty T_{mn}(t)R_{mn}(r)\Theta_m(\theta)$$ $$h(r,\theta,t):=\sum_{m,n=0}^\infty H_{mn}(t)R_{mn}(r)\Theta_m(\theta), \quad H_{mn}(t)=\frac{\int\limits_0^b \int\limits_0^a h(r,\theta,t) R_{mn}(r) \Theta_m(\theta) r dr d\theta}{\int\limits_0^a R_{mn}^2(r) r dr \int\limits_0^b \Theta_m^2(\theta) d\theta}$$ Substitute the new definitions into the non-homogeneous equations: $$ \begin{align} \sum T_{mn}'(t)R_{mn}(r)\Theta_m(\theta) = & k \left [ \frac{1}{r} \left ( \sum T_{mn}(t)R_{mn}'(r)\Theta_m(\theta) + r \sum T_{mn}(t)R_{mn}(r)\Theta_m(\theta) \right ) + \frac{1}{r^2} \sum T_{mn}(t)R_{mn}(r)\Theta_m(\theta) \right ] \\ & + \sum H_{mn}(t)R_{mn}(r)\Theta_m(\theta) \end{align} $$ We will use the following substitutions in our equation above: $$ \begin{cases} \Theta_m''(\theta)=-\mu_m^2\Theta_m(\theta) \\ rR_{mn}''(r)+R_{mn}'(r)=\frac{1}{r} \left ( \mu_m^2 - \lambda_{mn}^2r^2 \right ) R_{mn}(r) \end{cases} $$ We can eliminate the derivatives by substituting: $$ \begin{align} \sum T_{mn}'(t)R_{mn}(r)\Theta_m(\theta) & = \sum \left \{ T_{mn}(t) k \left [ \frac{1}{r} \underbrace{\left ( R_{mn}'(r)+rR_{mn}(r) \right ) }_{\frac{1}{r} \left ( \mu_m^2-\lambda_{mn}^2r^2 \right ) R_{mn}(r)} \Theta_m(\theta) + \frac{1}{r^2} R_{mn}(r) \underbrace{\Theta_m(\theta)}_{-\mu_m^2 \Theta_m(\theta) } \right ] \right \} + \sum H_{mn}(t)R_{mn}(r)\Theta_m(\theta) \\ & = \sum \left \{ T_{mn}(t) k \left [ \frac{1}{r^2} \left ( \mu_m^2-\lambda_{mn}^2r^2 \right ) - \frac{1}{r^2} \mu_m^2 \right ] \right \} R_{mn}(r) \Theta_m(\theta) + \sum H_{mn}(t)R_{mn}(r)\Theta_m(\theta) \\ & = \sum \left [ \left ( -k\lambda_{mn}^2 \right ) T_{mn}(t)+H_{mn}(t) \right ] R_{mn}(r) \Theta_m(\theta) \end{align} $$ From the linear independence of $$R_{mn} \otimes \Theta_m$$, it follows that: $$T_{mn}'(t) + k\lambda_{mn}^2 T_{mn}(t) = H_{mn}(t)~.$$ This first-order ODE can be solved with the following integration factor: $$\mu(t) = e^{k\lambda_{mn}^2t} $$ Thus, the equation becomes: $$\left [ e^{ k \lambda_{mn}^2 t } T_{mn}(t) \right ]' = e^{ k \lambda_{mn}^2 t} H_{mn}(t)$$ $$\Rightarrow T_{mn}(t) = e^{-k \lambda_{mn}^2 t} \int\limits_0^t e^{k \lambda_{mn}^2 t} H_{mn}(s) ds + C_{mn} e^{-k \lambda_{mn}^2 t}$$ We satisfy the initial condition: $$ \begin{align} u(r,\theta,0) & = f(r,\theta) \\ & = \sum_{m,n=0}^\infty T_{mn}(0)R_{mn}(r)\Theta_m(\theta) \\ & = \sum_{m,n=0}^\infty C_{mn}R_{mn}(r)\Theta_m(\theta), \quad C_{mn}=\frac{\int\limits_0^b \int\limits_0^a f(r,\theta) R_{mn}(r) \Theta_m(\theta) r dr d\theta}{\int\limits_0^a R_{mn}^2(r) r dr \int\limits_0^b \Theta_m^2(\theta) d\theta} \end{align} $$