# Does anyone know how to add a forcing term at the center of a cicular membrane?

I am here once again searching for wisdom, as some of you might notice I asked a related question a few days ago. And now I am struggling to add a forcing term at the center of the membrane, in order to recreate the Circular Chaldni figures, unless as at first approximation.

The idea is that with my original code see working code for arbitrary modes as initial condition if it works for any initial condition, it might work for a forcing term at the center (spoiler it is not) that is because the forcing terms will excite at certain combination of modes.

Here is my code and I will add some images.

# -*- coding: utf-8 -*-
"""
Created on Thu Nov  2 11:28:50 2023

@author: Manuel
"""

import numpy as np
from scipy.special import jv, jn_zeros
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

#%% Parameters
Nr = 50
N_phi = 50
N_steps = 100
c = 10
dphi = 2*np.pi/N_phi
dr = 5./Nr
dt = 0.001
CFL=c*dt/min(dr,dphi)
print(CFL)

#%% 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)
u = np.zeros((N_steps, Nr, N_phi))

#%% Forzante
u = np.zeros((N_steps, Nr, N_phi))
# Create an array to store the results
f = np.zeros((N_steps, Nr, N_phi))

# Create an array for the time values
t= np.linspace(0, dt * N_steps*10, N_steps)
frec=10
m=1
n=1

# Set the center element in each time step
for i in range(N_steps):
f[i, 0,0] = np.cos(2*np.pi*t[i]/N_steps)

#%% Plot de combinación inicial
fig = plt.figure()
M = f  # Or choose another time step to visualize

ax.plot_surface(X.T, Y.T, M, cmap='viridis')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
#%% Stepping
k1 = c*dt**2/dr**2
for t in range(0, N_steps-1):
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
u[t+1, i, j] = 2*u[t, i, j] - u[t-1, i, j] \
+ k1*(u[t, i+1, j] - 2*u[t, i, j] + u[t, i-1, j])\
+ k2*(u[t, i+1, j] - u[t, i-1, j])\
+ k3*(u[t, i, j+1] - 2*u[t, i, j] + u[t, i, j-1])+ f[t,i,j]
u[t+1, i, -1] = u[t+1, i, 0]  # Update the values for phi=2*pi

#%% Plot 11
fig = plt.figure()
Z = u  # Or choose another time step to visualize # Crece hasta 150 maximo en 610
# Recien en 900 empieza a cambiar el sentido
ax.plot_surface(X.T, Y.T, Z, cmap='viridis')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
#%% Forzante
fig = plt.figure()

# Inicializa la superficie 3D
Z_3d = u
surf = [ax3d.plot_surface(X.T, Y.T, Z_3d, cmap='viridis')]
# Set Z-axis limits
u_flat = u.flatten()

maximo= np.max(u_flat)
minimo=np.min(u_flat)
if abs(maximo)<=abs(minimo):
ax3d.set_zlim(-abs(minimo), abs(minimo))
else:
ax3d.set_zlim(-maximo, maximo)

ax3d.set_title("Frecuencia{} ".format(frec))
ax3d.set_xlabel('X')
ax3d.set_ylabel('Y')
ax3d.set_zlabel('Z')
ax3d.view_init(azim=70, elev=13)
# Función de actualización para la animación
def update(i):
Z_3d = u[i]  # Datos para la iteración actual
surf.remove()  # Elimina la superficie anterior
surf = ax3d.plot_surface(X.T, Y.T, Z_3d, cmap='viridis')  # Actualiza el gráfico 3D

ani = FuncAnimation(fig, update, frames=1000, interval=1, blit=False, repeat=False)
# Set the view angle (azimuth and elevation)

plt.show() I look at the final u vector, and I do not understand why it has so many 0 values, for the axis time and for the angular axis all are 0 except for the first row, that does not make any sense. I do not understan why it is not workin, the only thing that I know is that something is very wron with my code.

Once again, my mates and I have successfully deciphered how to tackle this. It's worth noting that when introducing the forcing term, altering 'i' impacts the radial position, while adjusting 'j' influences the angle. Instead of generating a Dirac forcing term, we opted for a more realistic approach by employing "if i in range(0,2)" to create a full cylinder. This is undoubtedly a superior choice, as a Dirac function isn't representative of real-world scenarios.

Please note that the cylinder can only be performed if i varies from 0 to n.

And is someone figure out how to do the same with a circular membrane that does not have fixed boundary condition let me know.

The code is the following:

import numpy as np
from scipy.special import jn_zeros
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

#%% Parameters
Nr = 50
N_phi = 50
N_steps = 1000
c = 10
dphi = 2*np.pi/N_phi
dt = 0.0001
CFL=c**2*dt/min(dr,dphi)
print(CFL)

#%% 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)

#%% MODE 01
m=0
n=1

print("la frecuencia es",frec,"modo",m,n)

#frec=257

#%% Stepping
u = np.zeros((N_steps, Nr, N_phi))
k1 = (c*dt)**2/dr**2
for t in range(0, N_steps-1):
for i in range(0, Nr-1):
for j in range(-1, 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
u[t+1, i, j] = 2*u[t, i, j] - u[t-1, i, j] \
+ k1*(u[t, i+1, j] - 2*u[t, i, j] + u[t, i-1, j])\
+ k2*(u[t, i+1, j] - u[t, i-1, j])\
+ k3*(u[t, i, j+1] - 2*u[t, i, j] + u[t, i, j-1])
if i==12 and j==29:
u[t+1,i,j]=u[t+1,i,j]+1/5*np.sin(frec*dt*t)

#%% Plot 11
fig = plt.figure()
Z = u  # Or choose another time step to visualize # Crece hasta 150 maximo en 610
# Recien en 900 empieza a cambiar el sentido
ax.plot_surface(X.T, Y.T, Z, cmap='viridis')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

#%% Animation

fig = plt.figure()

# Inicializa la superficie 3D
Z_3d = u
surf = [ax3d.plot_surface(X.T, Y.T, Z_3d, cmap='viridis')]
# Set Z-axis limits
u_flat = u.flatten()

maximo= np.max(u_flat)
minimo=np.min(u_flat)
if abs(maximo)<=abs(minimo):
ax3d.set_zlim(-maximo, maximo)
else:
ax3d.set_zlim(-maximo, maximo)

ax3d.set_zlim(-200, 200)
#ax3d.set_title("Frecuencia{} ".format(frec))
ax3d.set_xlabel('X')
ax3d.set_ylabel('Y')
ax3d.set_zlabel('Z')
ax3d.view_init(azim=70, elev=13)
# Función de actualización para la animación
def update(i):
index = i * 10
if index < len(u):
Z_3d = u[index]
surf.remove()
surf = ax3d.plot_surface(X.T, Y.T, Z_3d, cmap='viridis')

ani = FuncAnimation(fig, update, frames=range(N_steps), interval=10, blit=False, repeat=False)

plt.show()
#%% Store
# Specify the filename for the GIF