Physics Stack Exchange is a question and answer site for active researchers, academics and students of physics. It's 100% free.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I want to discretize the wave equation

$$\frac{1}{c^2}\frac{\partial^2\psi\left(\vec{r},t\right)}{\partial t^2}=\triangle\psi\left(\vec{r},t\right)$$

in polar coordinates. I find the following relations:

$$\vec{r}=r\mathbf{e}_r$$ $$\triangle=\frac{1}{r}\partial_r r \partial_r+\frac{1}{r^2}\partial_\varphi^2$$


$$r=i\Delta r,\qquad\varphi=j\Delta\varphi,\qquad t=n\Delta t$$


$$\psi\left(r,\varphi,t\right)=\psi\left(i\Delta r,j\Delta \varphi,n\Delta t\right)\equiv\psi_{i,j}^{\left(n\right)}$$

Using the forward difference in $r$ and central differences in $\varphi$ and $t$, I get:

$$\psi_{i,j}^{\left(n+1\right)}=2\psi_{i,j}^{\left(n\right)}-\psi_{i,j}^{\left(n-1\right)}+\frac{c^2 \Delta t^2}{\Delta r^2}\left[\psi_{i+2,j}^{\left(n\right)}+\frac{1}{i}\psi_{i+2,j}^{\left(n\right)}-2\psi_{i+1,j}^{\left(n\right)}-\frac{1}{i}\psi_{i+1,j}^{\left(n\right)}+\psi_{i,j}^{\left(n\right)}\right]+\frac{c^2\Delta t^2}{i^2\Delta r^2\Delta \varphi^2}\left[\psi_{i,j+1}^{\left(n\right)}+\psi_{i,j-1}^{\left(n\right)}-2\psi_{i,j}^{\left(n\right)}\right]$$

Now this is divergent at $r=0$, did I something wrong?

And how can I set restrictions to the steps $\Delta r$, $\Delta \varphi$ and $\Delta t$?

share|cite|improve this question
up vote 6 down vote accepted

If you look at the Laplacian: $$ \nabla^2=\frac{1}{r}\,\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\right)+\frac1{r^2}\frac{\partial^2}{\partial\phi^2} $$ you can clearly see that this diverges at $r=0$ so discretization of this should also diverge.

There are three solutions to remedying the divergent feature that I can think of:

  1. Choose a problem such that you don't need to worry about $r=0$ (e.g., $r_{min}=0.1$)
  2. Use cell-centered values for $r$ (e.g., $r_{min}=\frac12dr$)
  3. Use two meshes a Cartesian one for $r<1$ and a polar one for $r\geq1$, merging them at $r\simeq1$ (requires at least 1 cell of overlap).

I've never used option #3, so I cannot speak anything about it. Option #1 works in certain cases, but not for all. Option #2 I have worked with and have seen in other hydrodynamics codes that I've used, so it might be your best bet.

As far as your discretization goes, I think it might be more clear to use $r_i$ and $\phi_j$ for the positions and $dr$ and $d\phi$ for the cell widths. Using this, $$ \frac{1}{r}\,\frac{\partial}{\partial r}\left(r\frac{\partial\psi}{\partial r}\right)=\frac{1}{r_i}\left(r_{i+1/2}\frac{\psi^n_{i+1,j}-\psi^n_{i,j}}{dr^2}-r_{i-1/2}\frac{\psi^n_{i,j}-\psi^n_{i-1,j}}{dr^2}\right) \\ \frac{1}{r^2}\frac{\partial^2\psi}{\partial\phi^2}=\frac{1}{r_i^2}\left(\frac{\psi^n_{i,j+1}-2\psi^n_{i,j}+\psi^n_{i,j-1}}{d\phi^2}\right) $$ where $r_{i+1/2}=\frac12\left(r_{i+1}+r_i\right)$.

Your choice for $dr$ and $d\phi$ are up to you. It might be easier to have $dr=d\phi$, but it might not be worth enforcing that behavior. The timestep, however, is limited by the Courant-Friedrichs-Levy condition. For the wave equation, your requirement is such that $$ c^2\frac{dt}{\min(dr^2,d\phi^2)}\leq\frac12 $$ Formally, the constant can be 1 instead of $\frac12$, but for multidimensional systems, it is appropriate to use $1/n_{dim}$. For time-stepping purposes, this is going to be the slow point of your code. It might be worth looking into implicit methods instead of the explicit method you are currently working with.

share|cite|improve this answer
thanks for your answer! but I have some questions: 1) Does Option #2 ensure continuity at $r=0$? 2) Why does the CFL-condition not contain the wave speed? And why the maximum of $dr^2$ and $d\phi^2$? 3) Actually, I have to simulate a drum. It's a problem with polar symmetry, could I leave out the dependence of $\phi$? – Andy May 15 '14 at 21:24
(1) it will so long as your scheme will conserve it. (2) it should contain $c^2$, I just forgot it. You want the max to ensure the diffusion is appropriately captured. (3) you might be able to, but I cannot tell you for sure. – Kyle Kanos May 15 '14 at 21:28
@Andy: Actually, scratch that. You do want the minimum because the CFL conditions tells you how far something can be advected/diffused out of the cell before it is unphysical. Using the minimum of the two ensures that the physics is contained in the cell. – Kyle Kanos May 16 '14 at 1:05

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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