# How do you handle the singularity in polar or cylindrical coordinates?

Governing equations in polar or cylindrical coordinates often have terms with $$\frac{1}{r}$$ involved. At $$r = 0$$, such terms blow up to become a "singularity." The Cartesian version of such governing equations will not have these "singularities" at $$(x,y)=(0,0)$$. Therefore, the singularity at $$r = 0$$ must be numerical.

How do you handle the singularity at $$r = 0$$ in polar or cylindrical coordinates via e.g., finite difference?

To break this down:

• How are the radial coordinates discretized and differentiated?
• Is there is a transformation of variables of $$(r,\theta)$$ involved?
• How does any of this prevent blow-up at $$r = 0$$?

Let's assume the general case where there is no (anti-) symmetry to exploit and the governing equations rely on $$r$$ and $$\theta$$ explicitly, so the polar problem is at least two-dimensional. The governing partial differential equations could be the heat equation, Navier-Stokes equations, or Laplace's equation.

Some related answers are:

Neumann Boundary Condition at r=0 in Polar Coordinates (Numerical BCs)

Finite difference methods in cylindrical and spherical co-ordinate systems

• Using some small nonzero $r_{min}$ instead of $r=0$ is an easy solution acceptable for many applications. Aug 2, 2022 at 1:03
• So, the radial mesh would be something like {$r_{min}$, $r_{min} +\Delta r$, $+ ...$, $+ r_{max}$} or would it be something like {$-r_{max}$, $+...$, $-r_{min}$, $r_{min}$, $+ ...$, $+ r_{max}$}? If it is the former, does $r_{min}$ require a "boundary" condition? What would that be? If it is the latter, is there a transformation of variables involved (I'm thinking of needing to change the definition of $\theta$)? Aug 2, 2022 at 4:06
• A related question on Math.SE Aug 2, 2022 at 12:54
• One could use the grid, e.g., $\epsilon, \Delta r, 2 \Delta r, ..., r_{max}$. At $r_{min}=\epsilon$ the boundary condition would enforce Dirichlet or Neumann (or more complex) conditions that you'd use with $r_{min}$=0. Aug 2, 2022 at 14:27
• The problem you're having is called the "coordinate singularity" and is, among other places, addressed in this paper doi.org/10.1006/jcph.1999.6382
– user20857
Aug 2, 2022 at 17:03

Let me show this specifically for the finite element discretization of the Laplace equation: $$-\frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial}{\partial r} u(r,\theta)\right) - \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} u(r,\theta) = f(r,\theta).$$ The finite element method is based on the weak formulation, which you obtain by multiplying the equation by a test function $$\varphi(r,\theta)$$ and integrating over the domain, then integrate by parts. You'd think that that leads to $$\int_0^R \int_0^{2\pi} \varphi(r,\theta) \left[ -\frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial}{\partial r} u(r,\theta)\right) - \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} u(r,\theta) \right] \, d\theta dr = \int_0^R \int_0^{2\pi} \varphi(r,\theta) f(r,\theta) \, d\theta dr,$$ followed by integration by parts.
But you don't have to do it that way -- the above approach simply assumes that $$r$$ and $$\theta$$ are independent variables. Instead, you use the proper area element, which is $$2\pi r \, dr \, d\theta$$ and you get the following form instead: $$\int_0^R \int_0^{2\pi} \varphi(r,\theta) \left[ -\frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial}{\partial r} u(r,\theta)\right) - \frac{1}{r^2} \frac{\partial^2}{\partial \theta^2} u(r,\theta) \right] \, 2\pi r \, d\theta dr \\ = \int_0^R \int_0^{2\pi} \varphi(r,\theta) f(r,\theta) \, 2\pi r \, d\theta dr.$$ The $$2\pi$$ cancel on the two sides, and so what you have left is this: $$\int_0^R \int_0^{2\pi} \varphi(r,\theta) \left[ -\frac{\partial}{\partial r} \left(r \frac{\partial}{\partial r} u(r,\theta)\right) - \frac{1}{r} \frac{\partial^2}{\partial \theta^2} u(r,\theta) \right] \, d\theta dr = \int_0^R \int_0^{2\pi} \varphi(r,\theta) f(r,\theta) \, r \, d\theta dr.$$ After integration by parts, you then have the following (omitting boundary terms for simplicity): $$\int_0^R \int_0^{2\pi} r \frac{\partial \varphi(r,\theta)}{\partial r} \frac{\partial u(r,\theta)}{\partial r} + \frac{1}{r} \frac{\partial \varphi(r,\theta)}{\partial \theta} \frac{\partial u(r,\theta)}{\partial \theta} \, d\theta dr = \int_0^R \int_0^{2\pi} \varphi(r,\theta) f(r,\theta) \, r \, d\theta dr.$$ As you can see, at least from in front of the the $$r$$ derivatives, the singular weight has disappeared. If you choose appropriate quadrature rules, then the quadrature points for the second term will never lie at $$r=0$$, and so the still singular weight will not matter in the $$\theta$$-derivative term.
• How will the $1/r$ singular weight be made irrelevant after integration? To generalize, the method to treat the singularity is to somehow remove the $r = 0$ meshpoint from consideration? Aug 2, 2022 at 4:29
• I've never seen things done this way before. Could you point me to what a "proper" area element is? It seems to be what you'd get when you transform $dx dy$ into polar.
• @Nachiket Yes, exactly. It's just a coordinate transform. The integral is over a domain in $R^2$, and the integral should be the same whether you express it in $x-y$ coordinates or $r-\theta$ coordinates. The area element that includes the factor $2\pi r$ does that. Aug 3, 2022 at 2:40
• @Nachiket You can, but the question is whether you are able to integrate the other part by parts in a convenient way so that the resulting bilinear form becomes symmetric and positive definite. You will have noticed that that is the case with the $2\pi r$ weight but that that would not have been the case when using no weight at all. Aug 5, 2022 at 16:22