# CFL condition in polar coordinates.

Is there a Courant–Friedrichs–Lewy (CFL) type condition when dealing with the wave equation numerically in polar coordinates? For Example here it is mentioned for the wave equation but there isn't a explanation as to where it comes from. Any references on this subject would be greatly appreciated.

On Cartesian grids, the stability condition can be derived by inserting complex exponential solutions $\exp(i (k_x x + k_y y + k_z z + \omega t))$. From this you obtain a discrete dispersion relation $f(k_x,k_y,k_z,\omega)=0$ and the stability condition is the condition under which this equation cannot have complex $\omega$ solution (or at least no solutions with $\Im(\omega)<0$).

The reason this works at all is the happy coincidence that $\exp(i (k_x x + k_y y + k_z z + \omega t))$ is also an eigenfunction of the discrete derivative operators on a cartesian grid - it is an exact solution of the discretized equations.

In other coordinate systems, exact discrete solutions are not known, and deriving exact stability conditions is much more complicated.

That being said, there are ways of doing it. This thesis derives an exact stability condition for discretized Maxwell's equations in terms of (I oversimplify a little) eigenvalues of the discrete curl operator (eq 2.113), making no assumption on the coordinate system.

There is also this paper (eq. 31), but the derivation is not exact: they assume eigensolutions of the form $\exp(i (k_r r + k_{\phi} \phi + k_z z + \omega t))$ which is only valid in the far-from-the-origin limit.

I am speaking from memory here, but I'm sure Taflove's book "Computational Electrodynamics: The Finite-Difference Time-Domain Method" (not sure which edition) has a discussion of the stability of a discretized wave equation in cylindrical coordinates. They do not calculate it exactly, but they give numerical results that indicate the stability condition depends on how close to the origin you get.

One more thing: the formula in your link $$c^2\frac{dt}{\min(dr^2,d\phi^2)}\leq\frac12$$ cannot possibly be correct since $dr$ and $d\phi$ have different units.

• Very Helpful. Thanks – Harley Aug 1 '16 at 13:03

Have a look here: CFD-online-Wiki and references.

Let assume we are using standard structured Cartesian meshes (Cartesian meshes can be used, obviously, also with non-Cartesian coordinates) and a very basic numerical method (e.g. a stable and consistent first order finite-difference scheme). In general, the CFL stability condition is not known. For the simplest case, scalar equation in one dimension, it has shown to be something like $$| \Delta t | < \frac{h_{min}}{c_{max}}$$ where $\Delta t$ is the time-step (or the step of your integration), $h_{min}$ is the minimum mesh size of your numerical domain, $c_{max}$ is the maximum 'signal speed' (eigenvalue of the Jacobian of the flux tensor...). If $d$ is the dimension of your numerical domain, then the following form is commonly used $$| \Delta t | < \frac{1}{d}\frac{h_{min}}{c_{max}}.$$ Then the CFL in the link is ok, at least formally. Please, don't confuse numerics with physics. Here, the space-dimensions are dimensionless.

If you change coordinates, then you change the PDE (then also $c_{max}$) and the numerical domain ($h_{min}$). If you do things properly, the CFL restriction takes the same form, i.e. $$| \Delta t | < \frac{1}{d}\frac{h_{min}}{c_{max}}.$$

The problem at the origin ($r=0$, if $r$ is the radius of the polar coordinates) is that the Jacobian of the coordinate transformation between polar and Cartesian coordinates is singular (e.g. look at the angular velocity in $r=0$, for a generic constant velocity initial condition).

The general idea is that: the domain of dependece of the PDE must be contained in the domain of dependence of the discretization scheme.