diff --git a/methodes-spectrales/main.py b/methodes-spectrales/main.py new file mode 100644 index 0000000000000000000000000000000000000000..1577937bd2695a2fda83eedf6eae6b40145dab72 --- /dev/null +++ b/methodes-spectrales/main.py @@ -0,0 +1,47 @@ +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()