# Auteur : Clément de la Salle # Agrégation de physique, ENS de Lyon, 2019-2020 import numpy as np import matplotlib.pyplot as plt from cmath import * from matplotlib.animation import FuncAnimation from time import time from matplotlib.widgets import Slider, Button, RadioButtons # plt.rcParams['animation.ffmpeg_path'] = 'C:\\Program Files\\ffmpeg-20190709-5a481b1-win64-static\\bin\\ffmpeg.exe' def sans_dispersion(f) : '''Renvoie la vitesse de phase en fonction de la fréquence dans le cas sans dispersion''' return 1 def Klein_Gordon(f, omega_c = .99, c = 1) : '''Vitesse de phase en fonction de la fréquence pour Klein-Gordon avec une fréquence de coupure omega_c et une célérité c : k^2 = (omega^2 - omega_c^2) / c^2''' if f ** 2 > omega_c ** 2 : return c / (1 - omega_c**2 / (f**2)) ** .5 else : return 0 # Si omega < omega_c, l'onde ne se propage pas... On traduit ça par une vitesse de phase nulle (partie réelle de omega / k) def gravite(f, g = 1) : if abs(f) < 1e-12 : return 0 return g / f def groupe_constant(f, omega0, k0, vg = 1) : if f - omega0 + vg * k0 != 0 : return f * vg / (f - omega0 + vg * k0) else : return 0 class Onde : def __init__(self, X, Y, dispersion, *fargs) : ''' X : tableau des valeurs en abscisses (généralement np.linspace) Y : tableau des valeurs en ordonnées (deux exemples tout en bas du programme : gaussien et battements) dispersion : fonction renvoyant la vitesse de phase correspondant à une fréquence f (deux exemples : sans_dispersion et Klein-Gordon) si cette fonction prend d'autres arguments, on peut les ajouter à la suite (*fargs) ''' self.dispersion = dispersion # Vitesse de phase en fonction de la fréquence self.args_disp = fargs # Arguments supplémentaires pour la relation de dispersion (ex: omega_c et c pour Klein-Gordon self.Y = Y # Tableau des ordonnées self.X = X # Tableau des abscisses self.dx = (self.X.max() - self.X.min()) / (len(self.X) - 1) # Pas des abscisses # On calcul le spectre du paquet en entrée a = np.fft.ifftshift(self.Y) A = np.fft.fft(a) B = np.real(self.dx * np.fft.fftshift(A)) # on effectue un fftshift pour positionner la frequence zero au centre self.spectre = B # calcul des frequences avec fftfreq n = self.X.size freq = np.fft.fftfreq(n, d = self.dx) f = np.fft.fftshift(freq) self.spectre = np.concatenate((np.array([f]).transpose(), np.array([self.spectre]).transpose()), axis = 1) # Première colonne de spectre : fréquences # Deuxième colonne de spectre : amplitude associée def evol(self, Nmax, dt) : ''' Cette fonction calcul l'avolution du paquet d'onde sur Nmax itération avec un pas de temps dt. la résolution est entièrement analytique et ne fait appel à aucune intégration numérique, resserrer ce pas de temps ne rendra pas meilleur le calcul... Par contre l'animation sera plus fluide :) ''' hist = [] # Va contenir tous les tableaux des ordonnées à chaque instant histX = [] # Va contenir tous les tableaux des abscisses à chaque instant t = 0 vitesses = np.zeros((len(self.spectre))) # Va contenir les vitesses de phases associées aux fréquences du paquet for k, f in enumerate(self.spectre[:, 0]) : # On est obligé d'utiliser une boucle (et pas simplement des opérations d'array) à cause des cas particuliers (comme pour Klein-Gordon où seules certaines fréquences se propagent) vitesses[k] = self.dispersion(f, *self.args_disp) vg = vitesses.mean() # C'est à peu près la vitesse de groupe... Ce sera juste utile pour l'animation print(vg) debut_temps = time() division = 10 # Nombre de division du chargement ticks_chargement = np.linspace(0, Nmax, division) # On divise l'échelle en divisions curseur = 1 # On va repérer où l'on en est dans le chargement for k in range(Nmax) : profil = np.zeros((len(self.X))) # Pour chaque instant k, on va recréer le profil en sommant les cosinus du spectre (TF inverse) for [f, amplitude], v in zip(self.spectre, vitesses) : if abs(v) > 1e-10 : # Dans certains cas (comme Klein-Gordon par exemple), l'onde ne se propage pas... Par convention, on renvoie un vitesse nulle. Il faut donc exclure les fréquences qui ne se propagent pas. profil += amplitude * np.cos(f * (t - self.X / v)) # On reconstitue le paquet en sommant les contributions de chaque fréquence, en prenant en compte qu'elle ne se propagent pas toutes à la même vitesse if k > ticks_chargement[curseur] : print(''.join(['#' for k in range(curseur)]) + ''.join([' ' for k in range(division - curseur)]) + '| ' + str(round(time() - debut_temps, 1))) curseur += 1 t += dt hist.append(profil) histX.append(self.X) self.X = self.X + vg * dt self.hist = hist self.histX = histX def it_anim(self, k) : '''Fonction d'itération pour FuncAnimation''' self.line.set_ydata(self.hist[self.i]) self.line.set_xdata(self.histX[self.i]) self.i += 1 self.ax.set_xlim(self.histX[self.i][0], self.histX[self.i][-1]) # J'ai choisi de laisser déplacer la fenêtre avec le paquet d'onde... Certains peuvent préférer voir le paquet évoluer sur ton son espace, il suffit alors de retirer cette ligne return self.line, def animer(self, dt = 10) : '''Lance l'animation de l'onde en attendant un temps dt (en ms) entre chaque frame''' self.i = 0 self.dt = dt self.fig = plt.figure() self.ax = self.fig.add_subplot(111) self.line, = self.ax.plot(self.histX[0], self.hist[0]) mini = np.array(self.hist).min() maxi = np.array(self.hist).max() self.ax.set_ylim(mini - .1 * (maxi - mini), maxi + .1 * (maxi - mini)) # mini = np.array(self.histX).min() # maxi = np.array(self.histX).max() # self.ax.set_xlim(mini - .1 * (maxi - mini), maxi + .1 * (maxi - mini)) self.ax.set_xlim(self.histX[self.i][0], self.histX[self.i][-1]) self.ax.grid(True) case_start = plt.axes([0.55, 0.1, 0.1, 0.075]) start = Button(case_start, 'Start') start.on_clicked(self.start) plt.show() def start(self, event) : ani = FuncAnimation(self.fig, self.it_anim, interval = self.dt, frames = len(self.hist) - 2, repeat = False, blit = False) # ani.save('C:\\Users\\cleme\\Documents\\ENS\\Agreg\\Codes Pyhtons\\dispersion.mp4', fps=80, extra_args=['-vcodec', 'libx264'], savefig_kwargs={'facecolor':'white'}) plt.show() ## Petits exemples # Signal gaussien avec Klein-Grodon : # On voit les hautes fréquences se déplacer plus vite # On voit également des petits défauts devant le paquet, ceci est du au calcul de FFT qui n'est pas parfait """ T = np.linspace(-100, 100, 1000) # Tableau de l'abscisse (espace) x = np.exp(-1000 * T**2) # Tableau des ordonnées (fonction Gausienne) A = Onde(T, x, Klein_Gordon, .99, 1) # Création de l'onde constituée du paquet précédemment créé et de la relation de dispersion de Klein-Gordon (avec comme pulsation de coupure omega_c = 0.99 et comme célérité c = 1 : k^2 = (omega^2 - omega_c^2) / c^2 A.evol(200, .5) # Calcul de l'évolution temporelle : 200 itérations avec un pas de temps de 0.5 # la résolution est entièrement analytique et ne fait appel à aucune intégration numérique, resserrer ce pas de temps ne rendra pas meilleur le calcul... Par contre l'animation sera plus fluide :) A.animer() # Lance l'animation (possibilité de donner en argument le temps à attendre entre deux images de l'animation, par défaut 10 ms) """ # Battements avec Klein-Gordon : # On voit l'enveloppe se déplacer moins vite que la phase """ T2 = np.linspace(-200, 200, 1000) x2 = np.cos(2 * pi * T2) + np.cos(2.01 * pi * T2) B = Onde(T2, x2, Klein_Gordon, .99, 1) B.evol(400, .2) B.animer() """ # On jette un caillou dans l'eau """ C = Onde(T, x, gravite) C.evol(200, .5) C.animer() """ T = np.linspace(-10, 10, 1000) x = np.exp(-50 * T**2) * np.cos(50 * T) D = Onde(T, x, Klein_Gordon, 2, 1) D.evol(500, .5) D.animer() """ E = Onde(T, x, Klein_Gordon, .99, 1) E.evol(200, .5) E.animer() """