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?

  • $\begingroup$ Did you check these answers 1, 2? $\endgroup$
    – nicoguaro
    Commented Nov 8, 2021 at 1:32
  • $\begingroup$ @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? $\endgroup$
    – Alex I
    Commented Nov 8, 2021 at 2:11
  • $\begingroup$ 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 $\endgroup$
    – Charlie S
    Commented Nov 9, 2021 at 0:44

1 Answer 1


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.

From the introduction:

In the present paper we investigate a method for treating the coordinate singularity whereby singular coordinates are redefined so that data are differentiated smoothly through the pole, and we avoid placing a grid point directly at the pole. This eliminates the need for any pole equation. Despite the simplicity of the present technique, it appears to be an effective and systematic way to treat many scalar and vector equations in cylindrical and spherical coordinates.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.