So it is well known that the heat equation for a rod with perfectly insulated ends at $x= 0$ and $x = l$ is given by: $u_{t} = k u_{xx}$ It has the following BCs:

$\frac{\partial u}{\partial x}|_{(0,t)}=0 \ \hbox{ and } \frac{\partial u}{\partial x}|_{(l,t)}=0$

and IC :

$\frac{\partial u}{\partial x}|_{(x,0)}= h(x)$

Just wondering how the BCs for a spherically symmetric sphere with radius $r=a$ that is perfectly insulated around the boundary would look like. Should it look like this? If we only consider radial heat flow, the heat equation would be:

$u_t = k (u_{rr}+ \frac{2}{r} u_{r})$.

By substituting $v = ur$, we have: $v_t = k v_{rr}$

with the BCs given by: $v_{r}(0,t) = v_{r} (a,t) = 0$

and IC given by: $v = r h(r)$.

  • $\begingroup$ What do you mean by $v=ur$? $\endgroup$
    – K.defaoite
    Dec 1, 2020 at 16:18
  • $\begingroup$ i added a clarification in the question $\endgroup$
    – Christoph
    Dec 3, 2020 at 2:49
  • $\begingroup$ Are you sure you want $\partial_r[v](0,t)=\partial_r[v](a,t)=0$ and not $\partial_r[u](0,t)=\partial_r[u](a,t)=0$ ? The former will (I think) lead to complex, non-physical solutions. $\endgroup$
    – K.defaoite
    Dec 3, 2020 at 17:12
  • $\begingroup$ Updated with Neumann boundary conditions on $u$. $\endgroup$
    – K.defaoite
    Dec 4, 2020 at 0:17
  • $\begingroup$ You do not have a boundary at $r=0$. The boundary conditions will be $\partial_ru(a,t)=0$ and $|u(0,t)|$ bounded. $\endgroup$
    – Artem
    Dec 4, 2020 at 1:11

1 Answer 1


SO the heat equation is $$\partial_t u=\lambda\nabla^2u$$ With $\lambda$ being some diffusivity constant. We can use the Laplacian in spherical coordinates and remove the $\theta,\phi$ dependency to get $$\partial_t u=\frac{\lambda}{r^2}\partial_r\left(r^2\partial_ru\right)$$ Let's assume the equation is separable. $$\frac{G'(t)}{G(t)}=\frac{\lambda}{r^2F(r)}\partial_r\left(r^2F'(r)\right)=-c$$ Let's focus on the spacial part. Make the substitution $F(r)=f(r)/r$ so $$F'(r)=\frac{r\cdot f'(r)-f(r)}{r^2}$$ Which leads us to $$f''(r)=\frac{-c}{\lambda}f(r)$$ $$f(r)=a_1\sin\left(\sqrt{\frac{c}{\lambda}}r\right)+a_2\cos\left(\sqrt{\frac{c}{\lambda}}r\right)$$ Imposing the boundary conditions, namely $f(0)=f(a)=0$ we require $a_2=0$ and we get constants $$c_n=\lambda \frac{n^2\pi^2}{a^2}$$ Which gives our spacial eigenfunctions $$F_n(r)=A_n\operatorname{sinc}\left(\frac{n\pi}{a}r\right)=A_nj_0\left(\frac{n\pi}{a}r\right)=A_n\sqrt{\frac{a}{2nr}}J_{1/2}\left(\frac{n\pi}{a}r\right)$$ Using Bessel functions (this will be important later)

The time component follows: $$\frac{G'(t)}{G(t)}=-c=-\lambda \frac{n^2\pi^2}{a^2}$$ $$G_n(t)=B_n\exp\left(-\lambda \frac{n^2\pi^2}{a^2}t\right)$$

So we can write our solution as (using $A_n\sqrt{\frac{a}{2n}}\cdot B_n= C_n$) $$u(r,t)=\sum_{n=1}^\infty C_n\frac{J_{1/2}\left(\frac{n\pi}{a}r\right)}{\sqrt{r}}\exp\left(-\lambda \frac{n^2\pi^2}{a^2}t\right)$$ Initial conditions:

We can see that $$\sqrt{r}\cdot u(r,0)=\sum_{n=1}^\infty C_nJ_{1/2}\left(\frac{n\pi}{a}r\right)$$ This is a perfect time to use Fourier-Bessel series. Let $\alpha_{1/2,n}$ be the $n$th root of $J_{1/2}$, which in this case is the rather simple $\alpha_{1/2,n}=n\pi$. So $$\sqrt{r}u_0(r)=\sum_{n=1}^\infty C_nJ_{1/2}\left(\frac{\alpha_{1/2,n}}{a}r\right)$$ The literature then tells us that $$C_n=\frac{\int_0^a r^{3/2}u_0(r)J_{1/2}\left(\frac{\alpha_{1/2,n}}{a}r\right)\mathrm{d}r}{(a/2)J_{3/2}(\alpha_{1/2,n})^2}$$

EXAMPLE: Let $u_0(r)=(r^2-1)\log(r)$. The first few coefficients are $$\{C_n\}=\{0.624023,0.386661,0.10893,...\}$$ Here is a plot of the initial data $u_0(r)$ as well as its first five partial series- (r^2-1)log(r) series

And the Mathematica code used (so you can try yourself):

Subscript[u, 0][r_] := (r^2 - 1) Log[r]
c[n_] := 2*
   r^(3/2) Subscript[u, 0][r] BesselJ[1/2, n*\[Pi]*r], {r, 0, 1}]/
  BesselJ[3/2, n*\[Pi]]^2
Subscript[u, 0][r_, N_] := 
 Sum[c[n]*(Sqrt[2] Sin[n*\[Pi]*r])/(\[Pi]*r*Sqrt[n]), {n, 1, N}]
Plot[{Subscript[u, 0][r], Subscript[u, 0][r, 1], 
  Subscript[u, 0][r, 2], Subscript[u, 0][r, 3], Subscript[u, 0][r, 4],
   Subscript[u, 0][r, 5]}, {r, 0.01, 1}, 
 PlotLegends -> {"\!\(\*SubscriptBox[\(u\), \(0\)]\)(r)", 
   "\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;1)", 
   "\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;2)", 
   "\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;3)", 
   "\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;4)", 
   "\!\(\*SubscriptBox[\(u\), \(0\)]\)(r;5)"}]

Here $u_0(r;N)$ is the $N$th partial series. Another example when $u_0(r)=\sqrt{1-r}$ :


EDIT: What if we have Neumann, not Dirichlet, boundary conditions? That is, what if we want $$\partial_r[u](0,t)=\partial_r[u](a,t)=0$$ Well, let's return to our spacial part of our solution: $$F(r)=a_1\frac{\sin\left(\sqrt{\frac{c}{\lambda}}r\right)}{r}+a_2\frac{\cos\left(\sqrt{\frac{c}{\lambda}}r\right)}{r}$$ The derivative is $$F'(r)=\frac{\cos\left(\sqrt{\frac{c}{\lambda}}r\right)\left(a_1\sqrt{\frac{c}{\lambda}}r-a_2\right)-\sin\left(\sqrt{\frac{c}{\lambda}}r\right)\left(a_1+a_2\sqrt{\frac{c}{\lambda}}r\right)}{r^2}$$ In order to get $F'(0)=0$ we require $a_2=0$. Then our expression for the derivative simplifies to $$F'(r)=a_1\frac{\sqrt{\frac{c}{\lambda}} r\cos\left(\sqrt{\frac{c}{\lambda}}r\right)-\sin\left(\sqrt{\frac{c}{\lambda}}r\right)}{r^2}$$ What we're trying to do is find values of $c$ such that $F'(a)=0$ given $\lambda,a$. Unfortunately, the roots of $F(r)$ are not easily obtainable in closed form. Letting $\beta_{\lambda,n}$ be the $n$th positive root of $F'(r)$ given the constant $\lambda$ (we can find these via numerical algorithms) we require $$\beta_n=\sqrt{\frac{c_n}{\lambda}}a\implies c_n=\lambda \left(\frac{\beta_n}{a}\right)^2$$ So our spacial eigenfunctions are $$F_n(r)=A_n\operatorname{sinc}\left(\frac{\beta_n}{a}r\right)=A_n j_0\left(\frac{\beta_n}{a}r\right)=A_n\sqrt{\frac{\pi}{2\frac{\beta_n}{a}r}}J_{1/2}\left(\frac{\beta_n}{a}r\right)\to A_n\frac{1}{\sqrt{r}}J_{1/2}\left(\frac{\beta_n}{a}r\right)$$ Where of course we redefine $A_n$ as we go along for convenience. After this point, we follow the same process as before with the Fourier Bessel series and whatnot - just in this case we have to compute the roots $\beta_n$ beforehand.

  • $\begingroup$ I've just noticed you wanted $\partial_r[u](0,t)=\partial_r[u](a,t)=0$, not $u(0,t)=u(a,t)=0$. This makes things more difficult. I'll try to update the answer soon. $\endgroup$
    – K.defaoite
    Dec 3, 2020 at 16:40

You must log in to answer this question.

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