Finite difference methods in cylindrical and spherical co-ordinate systems

Question

I am quite familiar with finite difference schemes in cartesian coordinates. The key point here is that every point in the cartesian grid is treated equally as the spacing between consecutive points is same.

I want to know how one would perform finite-differencing in cylindrical (or even spherical) systems. I believe my main confusion is with angular differentiation. If we take a 2D cylindrical (polar) system, one way to divide the grid would be to make concentric circles (of $\Delta r$ spacing). For the angular spacing, we can draw radially outgoing rays, each of angular width $\Delta\phi$.

With $O(h^2)$ central differencing, the Laplacian, for example, can be given by:

$$ \nabla^2 f = 0$$ $$\Rightarrow \frac{f_{i+1,j} + f_{i-1,j} - 2f_{i,j}}{\Delta r^2} + \frac{1}{i\Delta r}\frac{f_{i+1,j} - f_{i-1,j}}{\Delta r} + \frac{1}{(i\Delta r)^2}\frac{f_{i,j+1} + f_{i,j-1} - 2f_{i,j}}{\Delta\phi^2} = 0$$

But in such a gridding scheme, as we increase $i$ (and hence $r$), the distance between two points on the same concentric circle will keep increasing. Is this how finite differencing works in cylindrical coordinates? Is such a scheme stable or does it become unstable for large $r$?

Are there better methods of finite differencing?

Answer 1

There are basically two methods: You can disrectize the angular part via grid points, or you can discretize it via basis expansion. I will focus on spherical symmetry here, the cylindrical case is quite similar.

In the basis expansion approach, one applies the ansatz $$ v(x,t) = \sum_{klm} a_{klm}(t)\,R_{klm}(r)\, Y_{lm}(\theta,\phi) $$ This is inserted into the problem to obtain equations for the radial function $R_{klm}(r)$, which are usually coupled in the angular indices $lm$. These equations are then solved by usual one-dimensional finite differences.

The other variant is to use an anguar grid. Here you can, in principle, use any grid you come up with. Standard choices are an equally spaced grid in the polar coordinate (x-y-plane), and a Legendre grid in the azimutal coordinate (z-direction). (The reason for the Legendre grid is the Jacobi determinant $r^2 sin\theta$ that arises in the transformation, by using Legendre grid points the $sin$ basically drops out).

The previous ansatz uses a product grid for the angles $\theta$ and $\phi$. Other, more sophisticated approaches try to reduce the required grid points by using specialized non-product grids. The grid points obtained in Lebedev quadrature is one popular alternative.

Here is a picture from a work of mine which shows (a) the Lebedev points and (b) the product grid: