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
radius = 1
c = 10
dphi = 2*np.pi/N_phi
dr = 5./Nr
dt = 0.001
CFL=c*dt/min(dr,dphi)
print(CFL)
#%% 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)
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()
ax = fig.add_subplot(111, projection='3d')
M = f[0] # 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()
ax = fig.add_subplot(111, projection='3d')
Z = u[90] # 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()
ax3d = fig.add_subplot(111, projection='3d')
# Inicializa la superficie 3D
Z_3d = u[0]
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[0].remove() # Elimina la superficie anterior
surf[0] = 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.
