1
$\begingroup$

I am currently employing a centered finite difference method to simulate the excitation of a circular membrane by a forcing term in Python, and the code is functioning as intended. To further enhance the simulation, I am considering the incorporation of a 3D Fourier transform. This would allow me to analyze the response of the circular membrane at different frequencies, such as when the frequency is set to 257, to observe the excitation of specific normal modes.

Could you provide guidance on calculating and plotting the 3D Fourier transform in Python for this purpose?

This is the array where the wave data is being stored

u = np.zeros((N_steps, Nr, N_phi))

and this is my code

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 = 10000 radius = 1 c = 10 dphi = 2*np.pi/N_phi dr = radius/Nr dt = 0.0001 CFL=c**2*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) #%% MODE 01 m=0 n=1 frec= jn_zeros(m, n)[n-1]*c/radius 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() ax = fig.add_subplot(111, projection='3d') Z = u[999] # 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() 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(-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) # Function that animates def update(i): index = i * 10 if index < len(u): Z_3d = u[index] surf[0].remove() surf[0] = 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()

I tried this

#%% Sección de FFT from scipy.signal import find_peaks FS= np.fft.fftn(u) # Calculate the magnitude spectrum along each axis fouriertime = np.abs(FS.sum(axis=(1, 2))) plt.plot(fouriertime) plt.title('Transforamda temporal'), plt.xlabel('Time Step') subfouriertime =fouriertime[:300] # Find peaks in the subset peaks, _ = find_peaks(subfouriertime) # Plot the spectrum with identified peaks plt.figure() plt.plot(subfouriertime, label='Temporal Transform') plt.plot(peaks, subfouriertime[peaks], 'x', label='Peaks > 0.2e7') plt.xlabel('Time Step') plt.ylabel('Magnitude') plt.legend() plt.show()

But as you see, I do not get a peak at 24.04 hz which is the frequency of mode 01.

I know that me code is not the best but I am new with these simulations.

This is the result of the Fourier transform

Fourier tranform

And if someone is wondering what is the result of the wave propagation:

Mode 01

Improve this question
$\endgroup$
7
  • $\begingroup$ Hello and welcome to Signal Processing SE. Although I believe your question is better suited for the Physics SE, or the Computational Science SE, I’m pretty sure there’s people here that could provide help. It would be extremely helpful though, if you could add some plots from your code and describe them along with what “seems to be” wrong and what you expected as a “correct” result, so that people won’t have to run your code to “understand” your problem. Furthermore, documenting you code (variables especially) to make a connection to a formulation would also be very helpful. $\endgroup$ Commented Nov 19, 2023 at 12:12
  • $\begingroup$ Hi, I asked at Stackoverflow and they suggested that I have to ask here. $\endgroup$ Commented Nov 19, 2023 at 13:59
  • $\begingroup$ I don’t have time to go through your code, but there seems to be a very strong DC component (frequency boin #0). Try removing the mean(s) from your signal before taking its FFT and see if that helps. $\endgroup$ Commented Nov 19, 2023 at 16:49
  • 1
    $\begingroup$ By removing the mean from your signal… sorry but I’m not sure how else to describe this… for a 1d signal it would be something like y= x-mean(x) $\endgroup$ Commented Nov 19, 2023 at 18:05
  • 1
    $\begingroup$ So, do you see how a 1D DFT is extended into a 2D DFT? $\endgroup$ Commented Nov 20, 2023 at 4:12

1 Answer 1

Reset to default
1
$\begingroup$

Do not expect to see anything interesting in the Fourier transform spectra, whether 3D or 1D spectra, of forced oscillations with resonant frequency. Any natural oscillations that may (and should) arise at the moment when the external force is exerted will be concealed by ever increasing resonant oscillation as there is no damping in your membrane problem.

This code shows what happens in time domain. You'll easily understand some other minor changes within loops. I also renamed freq into ω to make clear that the frequency parameter used with external force and natural mode oscillation functions is ANGULAR FREQUENCY measured in radian per second.

#%% Parameters Nr = 50 N_phi = 50 N_steps = 10000 radius = 1 c = 10 dphi = 2*np.pi/N_phi dr = radius/Nr dt = 0.0001 CFL=c*dt/min(dr,dphi*0.5*dr) print(CFL) if (CFL>=1): print("stop calculation, quit") quit() #%% 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) #%% MODE 01 m=0 n=1 omega= jn_zeros(m, n)[n-1]*c/radius print("la frecuencia es ",omega,"modo es [",m,n,"]") #%% Stepping u = np.zeros((N_steps, Nr, N_phi)) k1 = (c*dt)**2/dr**2 #u[0,12,29] = 1 for t in range(0, N_steps-1): for i in range(0, Nr-1): ri = max(r[i], 0.5*dr) # excluding i=0 to avoid the singularity k2 = (c*dt)**2/(2*ri*dr) for j in range(-1, N_phi-1): 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(omega*dt*t) uu = u[:,12,29] plt.figure('uu') plt.plot(uu) plt.show() uuFourier = np.fft.fft(uu) plt.figure('FFT(uu)') plt.plot(np.abs(uuFourier)) plt.show()

Notice that I corrected the Courant–Friedrichs–Lewy condition of the OP code: the code presented in the question travels across SE sites with its wrong formula for CFL. Also, while the special treatment of stencils at the coordinate origin helps the OP "avoid the singularity", this naive approach does nothing to eliminate excessively tiny steps (and consequently unnecessary computations) in this vicinity. The consistent approach, which eliminates the clustering of nodes near the origin, is presented in A Second-order Finite Difference Scheme For The Wave Equation on a Reduced Polar Grid (https://www.math.arizona.edu/~leonk/papers/polarFD7.pdf) by B. Holman, L. Kunyanski.

Back to the OP problem. In the beginning, the code prints

la frecuencia es 24.04825557695773 modo es [ 0 1 ]

in the end, the code shows the time-domain and FFT plots of the membrane motion at the point where the external force is exerted:

u_at_12,29

the FFT plot, zoomed-in, is:

FFT(u_at_12,29)zoomed

As the FFT bin index is an integer entity, we see the maximum value at an integer value of 4 Hz. Precise conversion of angular frequency [rad/s, oscillation function parameter] into temporal frequency [Hz, or s-1, the domain of FFT] gives the value of 24.048[rad/s]/2π = 3.827[Hz] as the frequency where the continuous FT spectrum of membrane motion has its maximum. To enhance the precision of numerical computation of pole/zero values, you should have more samples, that is, N_steps value should be increased.

The membrane motion at the other location (for example, uu = u[:,12,4], that is, at the point placed opposite) demonstrates a similar behavior, with only a starting time lag during which the excitation travels to this other point from the point where the force is exerted:

u_at_12,4

To see the very complex pattern of membrane motion, one can use a shock excitation, when the motion is started with only one point pushed away at a starting distance from its equilibrium position before it is let go the next moment. The code for this exercise:

#%% Parameters Nr = 50 N_phi = 50 N_steps = 10000 radius = 1 c = 1 dphi = 2*np.pi/N_phi dr = radius/Nr dt = 0.001 CFL=c**2*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) #%% MODE 01 m=0 n=1 omega = jn_zeros(m, n)[n-1]*c/radius print("la frecuencia es ", omega, "radian/s; ", "modo es ",m,n) #%% Stepping u = np.zeros((N_steps, Nr, N_phi)) k1 = (c*dt)**2/dr**2 u[0,12,29] = 1 for t in range(0, N_steps-1): for i in range(0, Nr-1): ri = max(r[i], 0.5*dr) # To avoid the singularity at r=0 k2 = (c*dt)**2/(2*ri*dr) for j in range(-1, N_phi-1): 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(omega*dt*t) uu = u[:,12,4] plt.figure('bang') plt.plot(uu) plt.show() uuFourier = np.fft.fft(uu) plt.figure('FFT(bang)') plt.plot(np.abs(uuFourier)) plt.show()

the plots:

bang

and it is a very complex behavior for useful analysis.

In reality, to be able to analyze a frequency spectrum of the membrane response to external force, you should scan the frequency of external force. You may want to examine a frequency response at all the membrane points, and the source data array for FFT analysis becomes 4D: u[extfreq, t, r, phi]. As you see, there is no magic in computing FFT for this data array: it is still a 1D FFT with extfreq, r, and phi parameters fixed and scanned over all the values you are interested in.

There exists a different problem of computing a Fourier transform of the membrane shape (vertical displacement values) when the external force of given frequency is exerted to membrane at a single point. The domain of this 2D Fourier transform is a 2D wave number k ([k_r, k_phi]), the parameters are the external frequency value and the time moment when we record the membrane shape. And still, there is no 3D Fourier transform in this task, only a 2D Fourier transform is computed.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Before you proceed with PDE research and computations, verify, if the CFL condition of your code is correct. If it is in error, this error was copy-pasted into my "exercises". If this is the case, please:1) Edit both your question and my answer, to secure that this error never migrates further on. Be responsible, always correct the noticed errors in the code of your posts. 2) Verify the other formulas containing velocity c, and if these are in error, correct the errors and compute with your original parameter values c=10, dt=0.0001. $\endgroup$ Commented Nov 25, 2023 at 5:49
  • 1
    $\begingroup$ The code is working without problem. But i need to study a bit more to understand the cfl condition.No sooner I understand that I will perform the corresponding modifications. $\endgroup$ Commented Nov 26, 2023 at 11:25
  • $\begingroup$ The fact that such an in depth and comprehensive answer only gets 1 upvote is frustrating to me. I wish active members of this community were better at recognizing efforts such as these. $\endgroup$ Commented Dec 9, 2023 at 14:53
  • $\begingroup$ The community is absolutely right, the question does not belong in the DSP realm. I'm regretting that I've got involved in this question here, but as this happened, I'm obliged to correct the noticed errors. For better or for worse, the SE posts are progressively invading web search results, and we are responsible for any misinformation (all the more disinformation) that may be communicated to visitors through our posts. $\endgroup$ Commented Dec 10, 2023 at 3:30

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.