# How can I solve wave equation for circular membrane in polar coordinates?

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

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.

• I get strange result. Your question not very informative. – sebix Jul 25 '15 at 19:13
• 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? – nicoguaro Jul 30 '15 at 13:35
• 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. – NoseKnowsAll Jul 30 '15 at 15:50
• 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. – nicoguaro Jul 30 '15 at 21:40
• 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". – Kirill Jul 31 '15 at 0:25

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
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
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')
mlab.show()


Timestep: 0 Timestep: 199 Timestep: 399 Timestep: 599 Timestep: 799 Timestep: 999

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.

• 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 – Henry Kolesnik Aug 2 '15 at 21:27