I have to solve the diffusion equation for a solid rod, which is in a bath of 100 degrees. At $t=0$ it is moved to a second bath, where the water temperature is 0 degrees.

I prepare the diffusion equation:

\begin{equation} u_t-\alpha\Delta=0. \end{equation}

Boundary and initial conditions:

\begin{equation} u(r,\theta,t)=0 \ \ \ \ 0<\phi<\pi, t>0 \\ u(r,\phi,0)=100 \ \ \ \ 0<r<R, 0<\phi<2\pi. \end{equation}

Since the function $u$ is of three variables $u=R\Theta T$, I prepare the given partial differential equation with the differential operator on a circle (the cross section of the rod):

\begin{equation} T_tR\Theta-\alpha\bigg[R_{rr}\Theta T+\frac{1}{r}R_r\Theta T+\frac{1}{r^2}\Theta_{\theta\theta}RT\bigg]=0 \end{equation}

Divide by $R\Theta T$ and obtain

\begin{equation} \frac{T_t}{T}-\alpha\bigg[\frac{R_{rr}}{R}+\frac{1}{r}\frac{R_r}{R}+\frac{1}{r^2}\frac{\Theta_{\theta\theta}}{\Theta}\bigg]=0 \end{equation}

Forming one ODE and one PDE:

\begin{equation} \frac{T_t}{T}=-\lambda^2 \\ \bigg[\frac{R_{rr}}{R}+\frac{1}{r}\frac{R_r}{R}+\frac{1}{r^2}\frac{\Theta_{\theta\theta}}{\Theta}\bigg]=\frac{\lambda_n^2}{\alpha} \end{equation}

The first gives $T(t)=Ce^{-\lambda^2t}$, the second partial differential equation gives two ordinary differential equations:

\begin{equation} r^2\frac{R_{rr}}{R}+r\frac{R_r}{R}+\frac{\Theta_{\theta\theta}}{\Theta}=\frac{\lambda_n^2}{\alpha}r^2 \end{equation}

which can be further split into two ODEs which must equal some constant $\mu^2$, which must be negative to give Bessel solutions:

\begin{equation} r^2\frac{R_{rr}}{R}+r\frac{R_r}{R}-\frac{\lambda_n^2}{\alpha}r^2=\frac{\Theta_{\theta\theta}}{\Theta}=-\mu_k^2 \end{equation}

this gives the Bessel equation:

\begin{equation} R_{rr}+\frac{1}{r}R_r+\bigg(\frac{\lambda_n^2}{\alpha}-\frac{\mu_k^2}{r^2}\bigg)R=0 \end{equation}

and the second-order ODE:

\begin{equation} \Theta_{\theta\theta}+\mu_k^2\Theta=0 \end{equation}

The latter has the solution $\Theta(\theta)=A\sin\mu_k\theta+B\cos\mu_k\theta$.

The former is the Bessel equation which gives the Bessel functions:

\begin{equation} R(r)=\bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_n r}{\sqrt{\alpha}}\bigg)+bY_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg) \end{equation}

So if this procedure is right, the general solution would be (ignoring the Bessel function of the second kind:

\begin{equation} u(r,\theta,t)=\sum_{n=1}^\infty \bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg)\bigg(\big(A\sin\mu_k\theta+B\cos\mu_k\theta\big)\big(e^{-\lambda_n^2t})\bigg) \end{equation}

Then we set the angular part equal to some constant, which is merged with the Bessel constant $a$ and get

\begin{equation} u(r,\theta,t)=\sum_{n=1}^\infty \bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg)\big(e^{-\lambda_n^2t})\bigg) \end{equation}

Using ICs to find the constants, the zeros of the Bessel are given by $\alpha_{n,k}$ and since the core will still be at 100 degrees at t=0, while the surface will be 0, then $R(0)=100$ and $R(R)=0$, and with the radial function given by:

\begin{equation} R(r)=aJ_{\lambda n}(i\sqrt{\mu}r) \end{equation}

we solve: \begin{equation} \alpha=aJ_{\lambda n}(i\sqrt{\mu}R) \end{equation}


\begin{equation} 100=aJ_{\lambda n}(i\sqrt{\mu}\alpha) \end{equation}

which gives:

\begin{equation} u(r,t)=\frac{\alpha J_{\lambda_n}(i\sqrt{\mu}r)}{J_{\lambda_n}\big(\frac{100R}{a\alpha}\big)}e^{-\lambda^2t} \end{equation}

The plot of the time-dependent evolution looks like this (with x=r) and x=0=R (surface) and t goes from 0 to 100. (ignore minus signs)

enter image description here

Does this entire procedure appear reasonable? Thanks

  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Xander Henderson
    Commented May 23, 2022 at 22:17

1 Answer 1


First let $u=X(\mathbf{x}) T(t)$. The equation reads $T_t X - \alpha \Delta(X) T = 0$. Take the cylindrical Laplacian and kill the $z$ and $\theta$ derivatives from symmetry considerations, so that $X(\mathbf{x})=R(\| \mathbf{x} \|)$; let $r=\| \mathbf{x} \|$.

Now you have

$$T_t R - \alpha T \left ( R_{rr} + r^{-1} R_r \right ) = 0.$$

Doing the usual separation of variables trick, you conclude that $\frac{T_t}{T}=\alpha \frac{R_{rr} + r^{-1} R_r}{R}=\lambda_n$ for the same constant $\lambda_n$. It follows that up to a constant factor that we'll sort out later, $T=e^{\lambda_n t}$.

Now in the radial part we have

$$R_{rr} + r^{-1} R_r - \frac{\lambda_n}{\alpha} R = 0.$$

This is actually an adjusted Bessel equation itself. Multiply through by $r^2$:

$$r^2 R_{rr} + r R_r + \left ( -\frac{\lambda_n}{\alpha} r^2-0^2 \right ) R = 0.$$

That $-\frac{\lambda_n}{\alpha}$ wouldn't be there in the ordinary Bessel equation. The fix comes because the first two terms are invariant under a linear change of variable, so we can replace $s=\sqrt{-\lambda_n/\alpha} r$ (dropping the $n$ as an abuse of notation) and then we have

$$s^2 R_{ss} + s R_s + \left ( s^2 - 0^2 \right ) R = 0.$$

This gives $R=c_{1,n} J_0(s) + c_{2,n} Y_0(s)$, Bessel functions of order zero. For nonsingularity at the origin we force $c_{2,n}=0$.

So now the question is about what $c_{1,n}$ is. The big catch is that in order for $J_0(\sqrt{-\lambda_n/\alpha} r)$ to have zeros, we must have $\sqrt{-\lambda_n/\alpha}$ real. (This also makes intuitive sense as it implies exponential decay.)

As a result you choose $\lambda_n$ so that $\sqrt{-\lambda_n/\alpha} R$ is a zero of $J_0$ (where $R$ is now the radius of the rod). We can pick the zero of interest to be $j_{0,n}$ which is where $n$ enters into the picture.

Now we solve that for the eigenvalue, and get $\lambda_n=-\frac{\alpha j_{0,n}^2}{R^2}$ which gives

$$u=\sum_n c_n J_0 \left ( j_{0,n} \frac{r}{R} \right ) e^{-\frac{\alpha j_{0,n}^2}{R^2} t}.$$

Now it remains to find the $c_n$ to match the initial profile.

  • $\begingroup$ Thanks Ian, I will study this on Monday! Have a good weekend. PS: Upscore my post! $\endgroup$ Commented May 21, 2022 at 15:41
  • 1
    $\begingroup$ @Luthier415Hz I fixed the mistake, it was partially my fault; you can in fact change the radial eigenvalue while keeping the angular eigenvalue frozen. You were moving them in lockstep and that confused me when it shouldn't have. $\endgroup$
    – Ian
    Commented May 21, 2022 at 15:48
  • $\begingroup$ No problem. I think that we have reached some knowledge here anyhow. $\endgroup$ Commented May 21, 2022 at 16:17

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .