The original equation is

$$\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial^2 u}{\partial \phi^2}$$

and the finite difference version used is

$$U^t_{i,j} = {c^2 \Delta t^2}\left[\frac{1}{2r_i\Delta r }(U^{t-1}_{i+1,j} - U^{t-1}_{i-1,j}) + \frac{1}{\Delta r^2}(U^{t-1}_{i-1,j} + U^{t-1}_{i+1,j} - 2U^{t-1}_{i,j}) + \frac{1}{r_i^2 \Delta \phi^2}{}(U^{t-1}_{i,j-1} + U^{t-1}_{i,j+1} - 2U^{t-1}_{i,j})\right] - U^{t-2}_{i,j} + 2U^{t-1}_{i,j}$$

with $c^2=0.25$, $\Delta r=0.1$, $\Delta \phi=0.1256$, $\Delta t=0.05$. The total radius is 5.

R = np.linespace(0, radius, 50)
phi = np.linspace(0, 2*np.pi, 50)
R, phi = np.meshgrid(R, phi)
x = R*np.cos(phi)
y = R*np.sin(phi)

thus dr=5./50=0.1, and dPhi=2*np.pi/50=0.1256.

My initial conditions looks like enter image description here

I create r[50][50], phi[50][50], and the same space for x, y ($x=r\cos(\phi)$, $y=r\sin(\phi)$).

I don't know why, but I get strange result.

Also I made a model for a Cartesian system, and it works.

  • 2
    $\begingroup$ I get strange result. Your question not very informative. $\endgroup$
    – sebix
    Commented Jul 25, 2015 at 19:13
  • 2
    $\begingroup$ I don't see how your equation can solve the wave equation, it does not show any iteration in time (or eigenvalue extraction). It is better if you write the equation in the form you are using to iterate. Also, I don't completely understand your second term, does that mean that you have dissipation in your system? $\endgroup$
    – nicoguaro
    Commented Jul 30, 2015 at 13:35
  • 1
    $\begingroup$ On top of what @nicoguaro says, you are missing a 2 as the coefficient for the $u_{i,j}$ term in the $u_{rr}$ term. $\endgroup$ Commented Jul 30, 2015 at 15:50
  • 1
    $\begingroup$ I suppose that the boundary conditions are $u(r=5)=0$. I also suppose that when $r$ tends to 0, the derivatives are 0 (you have a local extremum), what makes the bracket terms to vanish. You also should consider that for the angular variable you have periodicity, i.e., $u(0)=u(2\pi)$, $u_\phi(0)=u_\phi(2\pi)$, what should be reflected in the difference equations for the first and last value of $j$. As a final comment, you should reconsider the use of that colormap in the surface. $\endgroup$
    – nicoguaro
    Commented Jul 30, 2015 at 21:40
  • 1
    $\begingroup$ Why is $\Delta t^2$ on the bottom? Shouldn't it be on the top? I.e., smaller timesteps mean smaller adjustments, etc. The units are also wrong. Also, please give more details about your output than "I get strange result". $\endgroup$
    – Kirill
    Commented Jul 31, 2015 at 0:25

1 Answer 1


I think that you have two problems:

  • your timestep is too big (the CFL condition is not satisfied); and
  • you are not updating the values for $\phi=2\pi$.

I am not sure about the CFL condition for polar coordinates, but I believed that it should be $$2c\frac{\Delta t}{\Delta r \Delta \phi} \leq C_\max \enspace ,$$

what implies that $\Delta t < \frac{\Delta r \Delta \phi}{2c}$.

Regarding the values for $\phi=2\pi$, you need to update for $\phi=0$ and $\phi=2\pi$ although you are not looping until the last angular value.

These things are corrected in the code below (and it works). One thing that is not working terribly good is the solution for $r=0$... but I have no ideas now on how solve it. But I think that it needs to be solved as a particular case.

I attach some images for the mode $(1,1)$ of vibration (since I am not good with the animation thing in Python).

import numpy as np
from scipy.special import jv, jn_zeros
from mayavi import mlab

#%% Parameters
Nr = 50
N_phi = 50
N_steps = 1000
radius = 5.
c = 0.5
dphi = 2*np.pi/N_phi
dr = 5./Nr
dt = 0.005  
if dt< dr*dphi/2/c:
    # The maximum value for dt is dr*dphi/(2*c)
    dt = dr*dphi/(4*c)

#%% Initial conditions
r = np.linspace(0, radius, Nr)
phi = np.linspace(0, 2*np.pi, N_phi)
R, phi = np.meshgrid(r, phi)
X = R*np.cos(phi)
Y = R*np.sin(phi)
kth_zero = jn_zeros(1, 1)
Z = np.cos(phi) * jv(1, kth_zero*R/radius)
T = np.zeros((N_steps, Nr, N_phi))
T[0, :, :] = Z.T
T[1, :, :] = Z.T

#%% Stepping
k1 = c*dt**2/dr**2
for t in range(2, N_steps):
    for i in range(0, Nr-1):
        for j in range(0, N_phi-1):
            ri = max(r[i], 0.5*dr)  # To avoid the singularity at r=0
            k2 = c*dt**2/(2*ri*dr)
            k3 = c*dt**2/(dphi*ri)**2
            T[t, i, j] = 2*T[t-1, i, j] - T[t-2, i, j] \
            + k1*(T[t-1, i+1, j] - 2*T[t-1, i, j] + T[t-1, i-1, j])\
            + k2*(T[t-1, i+1, j] - T[t-1, i-1, j])\
            + k3*(T[t-1, i, j+1] - 2*T[t-1, i, j] + T[t-1, i, j-1])

        T[t, i, -1] = T[t, i, 0]  # Update the values for phi=2*pi

surf = mlab.mesh(X.T, Y.T, 10*T[999], colormap='RdYlBu')

Timestep: 0 enter image description here Timestep: 199 enter image description here Timestep: 399 enter image description here Timestep: 599 enter image description here Timestep: 799 enter image description here Timestep: 999 enter image description here

PS: Please, generate videos and share it in the comments (I want to see). And, use a different colormap, the one you used is not so good.

  • 2
    $\begingroup$ You are brilliant, absolutely brilliant. Amazing answer. Everything is working. I didn't take into account  3 things: 1. The time 2. updating the values 3. transposing X,Y,Z Thanks so much, amazing $\endgroup$ Commented Aug 2, 2015 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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