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

Question

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.

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')
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.