methodes spectrales MVP

methodes-spectrales/main.py

+from numpy import *
+import matplotlib.pyplot as plt
+
+
+def run_simulation(L, u0, dt, timesteps, nu=1):
+    N = len(u0)
+    # compute a vector storing the frequencies for the fourier transform, it already takes the interval [0, L]
+    # into account, no need to divde with / L later in the computations
+    k = fft.fftfreq(N, d=L/(N-1))
+    solution = zeros((timesteps+1, N))
+    solution[0] = u0
+
+    previous_u2k = fft.fft(u0 ** 2)
+    current_uk = fft.fft(u0)
+    # f_L (k), this is constant so we only compute it once.
+    # note that we don't need to divide by /L because fftfreq already does it :)
+    fL = (2*pi * k) ** 2 - nu * (2 * pi * k) ** 4
+
+    for i in range(timesteps):
+        # update solution (note that we don't divide by L!)
+        current_u2k = fft.fft(solution[i] ** 2)
+        next_uk = (1 + dt/2 * fL) / (1 - dt/2 * fL) * current_uk - 1j * pi * k * (3 * current_u2k - previous_u2k) / (2 - dt * fL) * dt
+        # store after inverse FFT, ignoring its imaginary part
+        solution[i+1] = real(fft.ifft(next_uk))
+        # update values for the next iteration
+        current_uk = next_uk
+        previous_u2k = current_u2k
+
+    return solution
+
+
+# define a few constants
+L = 10
+N = 1024
+dx = L/(N-1)
+x = linspace(0, L, N)
+time_values = linspace(0, 200, int(200 / 0.05)+1)
+# run the simulation
+s = run_simulation(L, cos(2*pi/L * x) + 0.1 * cos(4*pi/L * x), 0.05, len(time_values)-1)
+# pick a nice color for the plot
+colormap = plt.get_cmap("jet")
+# print the values
+plt.pcolormesh(x, time_values, s, cmap=colormap)
+# print the colorbar on the right
+plt.colorbar()
+# open a window with the plot :)
+plt.show()