# How to solve heat equation in spherical coordinates with finite differences?

I have a problem dealing with heat transfer which is spherically symmetrical. I was thinking it should be possible to solve this as a 1d problem in spherical coordinates using the radius only.

Normally, the 1d heat equation can be solved as:

$$u_n^{k+1} = u_n^{k} + \frac{\alpha \Delta t}{\Delta x^2} \left(u_{n+1}^{k} - 2 u_{n}^{k} + u_{n-1}^{k} \right)$$

with forward difference in time and central difference in space, where $$u_n^k$$ is the value at the n-th grid position at the k-th time step.

Now, in spherical coordinates, the Laplacian is $$\nabla^2 = \frac{\partial^2}{\partial r^2} + \frac{2}{r} \frac{\partial}{\partial r} + \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} + \frac{\cos \theta}{r^2 \sin \theta} \frac{\partial}{\partial \theta} + \frac{1}{r^2 \sin^2 \theta} \frac{\partial^2}{\partial \phi^2}$$

The derivatives wrt anything other than r are zero because of symmetry, so only the first two terms are needed. This gives a finite difference scheme along the lines of:

$$u_n^{k+1} = u_n^{k} + \frac{\alpha \Delta t}{\Delta r^2} \left(u_{n+1}^{k} - 2 u_{n}^{k} + u_{n-1}^{k} \right) + \frac{2}{r} \frac{\alpha \Delta t}{\Delta r} \frac{1}{2} \left(u_{n+1}^{k} - u_{n-1}^{k} \right)$$

(Also maybe the method for calculating the first derivative could be simply a forward difference like $$\left( u_{n+1}^{k} - u_{n}^{k} \right)$$ instead).

What is the right version of the CFL condition here? Would is be determined by the second derivative term, the first derivative, or (somehow) both at the same time?

Should the grid spacing $$\Delta r$$ vary with $$r$$? (perhaps as proportional to r)

Lastly, what is the right way to handle the r=0 boundary?

• Did you check these answers 1, 2? Nov 8, 2021 at 1:32
• @nicoguaro I hadn't seen those, thanks. I don't get the derivation of CFL for polar coordinates in the second answer, esp because there is both a first and second derivative term wrt r. Same CFL for both? Also, how to handle r=0, and what grid spacing to use? Nov 8, 2021 at 2:11
• L'Hopital's rule suggests that $\lim_{r\to 0} \frac{1}{r} \frac{\partial u}{\partial r} = \frac{\partial^2 u}{\partial^2 r}$. Perhaps you can use that on your r=0 boundary Nov 9, 2021 at 0:44

Deriving the CFL condition was handled by @nicoguaro in the comments. The coordinate singularity problem at $$r = 0$$ was addressed in this paper: "Numerical Treatment of Polar Coordinate Singularities" Mohseni + Colonius JCP (2000) and is available here.