[Math] Stability of the BTCS scheme for the heat equation in a disk

finite differencesnumerical methodspartial differential equations

Consider the $1$-D heat equation:
$$
u_t = a \Delta u = au_{xx} \\
u(0,t) = u(1,t) = 0 \\
u(x,0) = u_0(x)
$$
where $a > 0$ is constant and $u_0$ is given. It is a classic result that the implicit finite difference method BTCS is unconditionally stable.

Now, for $u = u(r)$, consider the analogous problem in the unit disk:
$$
u_t = a \Delta u = \dfrac{a}{r}\left(ru_r\right)_r = au_{rr} + \dfrac{a}{r}u_r \\
u(1,t) = 0 \\
u_r(0,t) = 0 \\
u(r,0) = u_0(r)
$$

If we assume we have a uniform mesh $0 \leq j \leq N$, $0 \leq i \leq M$, $h=1/N$, $k = 1/M$ with $U_j^i = u(jh,ik)$ then the BTCS scheme is
$$
\dfrac{U_0^{i+1} – U_0^i}{k} = 4a\dfrac{U_1^{i+1} – U_0^{i+1}}{h^2} \\
\dfrac{U_j^{i+1} – U_j^i}{k} = a\dfrac{U_{j+1}^{i+1} – U_j^{i+1} + U_{j-1}^{i+1}}{h^2} + a\dfrac{U_{j+1}^{i+1} – U_{j-1}^{i+1}}{2jh^2}, 1 \leq j \leq N
$$

I assume this is also unconditionally stable, but we cannot apply von Neumann analysis because the coefficients are no longer constant (and have a singularity at $r=0$). Could anyone confirm that this is in fact unconditionally stable and provide a reference or an outline how to show it? Perhaps by using an energy method in a weighted $L^2$-norm?

Best Answer

When you do the von Neumann analysis, I guess you end up with an amplification factor $G=U^{i+1}/U^i$ which depends both on $d=ka/h^2$ and $j$ (and $\theta$). So for stability you have to find $d$ which satisfies $|G(d,j,\theta)|\leq 1$, for $0<\theta<2\pi$ and $0<j<N$.

The equation you are solving is not the same as the 1D equation, so it's not obvious if it will also be unconditionally stable.