# Finite difference methods in cylindrical and spherical co-ordinate systems

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?

• What is the issue? E.g., for polar coordinates you have some operator in terms of $\partial{}/\partial{r}$, $\partial{}/\partial{\theta}$ and higher derivatives; and your domain is rectangular in $r,\theta$ coordinates, with grid points on $r$=const and $\theta$=const lines. Otherwise the same stuff as in Cartesian coordinates. – Maxim Umansky May 30 '20 at 4:21

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:

• A nice answer. Could you please provide some sources to your work / any other source that I can refer to for this stuff? – Siddharth Bachoti May 30 '20 at 15:36
• Also, it seems to me that the common method of taking a constant grid spacing like in case of cartesian co-ordinates ($\Delta x = \Delta y = \Delta z$) is not possible in curvilineaar co-ordinates. – Siddharth Bachoti May 30 '20 at 15:40
• My own work is just an application of these concepts to quantum mechanics, i.e. nothing that I would primarily refer to for these methods. Googling "Lebedev quadrature" is a good start into the topic. One ineresting paper I encountered by this is Beentjes, Quadrature on a spherical surface. -- note that many useful point distributions stem from quadrature rules (just as in 1D). – davidhigh May 31 '20 at 10:10