C'est une des équations les plus connues de la physique mathématique, que nous devons à Jean le Rond D'Alembert en 1747. Elle décrit la propagation d'une onde dans le vide ou un milieu matériel quelconque.
Nous avons rencontré cette équation dans la page de TangenteX.com consacrée aux équations de Maxwell et, comme promis, je vous en propose ici une méthode de résolution numérique.
Soit A un champ scalaire ou vectoriel non stationnaire. L'équation de D'Alembert s'exprime alors sous cette forme :
\( \Delta \: A - \frac{1}{c^2} \frac{\partial^2 A}{\partial t^2} = 0 \)
La grandeur c représente la célérité de l'onde. C'est une équation aux dérivées partielles hyperbolique d'ordre 2 en temps et en espace. A chaque variation touchant l'espace, exprimée par le laplacien, correspond une variation temporelle proportionnelle et opposée de la grandeur A.
L'équation de D'Alembert est une équation différentielle linéaire, ce qui est une caractéristique essentielle de la physique des ondes qui sont modélisables par cette équation. Cela signifie que n'importe quelle combinaison linéaire des solutions est aussi une solution. On désigne d'ailleurs cette équation comme l'équation des ondes linéaire et homogène, car il existe d'autres équations d'ondes, comme l'équation des télégraphistes, qui ne sont pas linéaires.
Une multitude de modélisations physiques débouchent sur l'équation de propagation des ondes élaborée par D'Alembert. Citons par exemple le mouvement d'une corde de guitare, ou la vibration d'une membrane de tambour ou encore d'une corde. Mais on peut aussi l'obtenir à partir des équations de Maxwell...
Partons de l'équations de Maxwell-Faraday en absence de charges électriques : \( \vec{rot} \: \vec{E} = -\frac{\partial \vec{B}}{\partial t} \).
Calculons \( \vec{rot} \: \vec{rot}\: \vec{E} = \vec{rot}(-\frac{\partial \vec{B}}{\partial t}) \). Or, parce que les variables sont indépendantes, je peux inverser l'ordre des dérivations et écrire \( \vec{rot}(-\frac{\partial \vec{B}}{\partial t}) = -\frac{\partial}{\partial t}\vec{rot}\vec{B} \).
L'équation de Maxwell-Ampère en l'absence de charges électriques me dit que \( \vec{rot} \: \vec{B} = \mu_0 \epsilon_0 \frac{\partial \vec{E}}{\partial t} \), d'où je peux écrire que \( \vec{rot} \: \vec{rot} \: \vec{E} = -\mu_0 \epsilon_0 \frac{\partial}{\partial t}\frac{\partial \vec{E}}{\partial t} \), soit encore :
\( \vec{rot} \: \vec{rot} \: \vec{E} = -\mu_0 \epsilon_0 \frac{\partial^2 \vec{E}}{\partial t^2} \).
En utilisant l'égalité vectorielle \( \vec{rot} \: \vec{rot} = \vec{grad} \: div - \Delta \), j'ai \( -\mu_0 \epsilon_0 \frac{\partial^2 \vec{E}}{\partial t^2} = \vec{grad}(div \: \vec{E}) - \Delta \: \vec{E} \). L'équation de Maxwell-Gauss me dit que \( div \: \vec{E} = 0 \) et donc que \( -\mu_0 \epsilon_0 \frac{\partial^2 \vec{E}}{\partial t^2} = -\Delta \: \vec{E} \). Ce qui me donne finalement l'équation :
\( \Delta \: \vec{E} \: - \:\mu_0 \epsilon_0 \frac{\partial^2 \vec{E}}{\partial t^2} = \vec{0} \).
L'équation de D'Alembert admet des solutions analytiques connues, qui sont de la forme \( A_i(x_i, t) = A_+(x_i - ct) + A_-(x_i + ct) \), où \( x_i \) représente les différentes coordonnées de l'espace. Les fonctions \( A_+ \) et \( A_- \) sont définies par les conditions initiales du problème. Elles représentent deux ondes progressives de vitesse finie respective -c et +c, qui se propagent dans deux directions opposées.
Vous vous souvenez sans doute, depuis votre cours sur l'analyse de Fourier, que n'importe quelle fonction (ou presque) peut être décomposée en sommes de sinus. Il est donc envisageable de proposer comme solution une fonction \( A(x_i,t) = a\sin(kx_i - \omega t) \). Si vous avez un peu de courage, vous pourrez vérifier que cette fonction est bien une solution de l' équation de D'Alembert. Vous pouvez aussi chercher d'autres fonctions solutions, par exemple du coté des gaussiennes...
Pour simplifier notre étude, je vais me limiter à l'équation de D'Alembert unidimensionnelle, projetée sur l'axe Ox. Je reprends les notations traditionnelles et je l'écris \( \frac{\partial^2 u(x,t)}{\partial t^2} - c^2 \frac{\partial^2 u(x,t)}{\partial x^2} = 0 \), u étant une fonction scalaire. Elle décrit par exemple le mouvement d'une corde vibrante. c est bien sur la célérité de l'onde.
Je suis en présence de ce qu'on appelle un problème de Cauchy homogène, que je résume ainsi :
\( \left\{ \begin{array}{lr} \frac{\partial^2 u(x,t)}{\partial t^2} - c^2 \frac{\partial^2 u(x,t)}{\partial x^2} = 0 \; avec \; t \gt 0 , \: x \in \mathbb{R} \\ u(x,0) = g(x) \\ \frac{\partial u(x,0)}{\partial x} = h(x) \; avec \; x \in \mathbb{R} \\ \end{array} \right. \)
Les conditions initiales du problèmes sont fixées par les fonctions h(x) et g(x) de classe \( C^1 \) et \( C^2 \) respectivement. On démontre (vous le verrez en cours) que ce problème admet une solution unique que l'on nomme "formule d'Alembert".
Vous aurez remarqué que le problème de Cauchy ne fixe pas de conditions aux limites. La solution se propage à l'infini. Dans de nombreux cas, par exemple dans le cas d'une corde vibrante ou d'ondes électromagnétiques dans une cavité, le domaine de variation de x n'est plus \( \mathbb{R} \) mais un intervalle borné, que l'on normalise généralement à [0,1]. A l'énoncé du problème de Cauchy, je vais donc ajouter les conditions aux limites, dites "conditions de bord de Dirichlet homogènes", ce qui me donne :
\( \left\{ \begin{array}{lr} \frac{\partial^2 u(x,t)}{\partial t^2} - c^2 \frac{\partial^2 u(x,t)}{\partial x^2} = 0 \; avec \; t \gt 0 , \: x \in ]0,1[ \\ u(x,0) = g(x) \\ \frac{\partial u(x,0)}{\partial x} = h(x) \; avec \; x \in ]0,1[ \\ u(0,t) = u(1,t) = 0 \; avec \; t \gt 0 \end{array} \right. \)
Ici, les conditions de bord sont nulles (Dirichlet homogène), mais on peut imaginer qu'elles ne le soient pas.
Je vous renvoie à la page décrivant la méthode des différences finies pour les aspects techniques des schémas aux différences finies.
Dans le cas présent de l'équation de D'Alembert, nous allons discrétiser les deux membres de l'équation qui sont de deuxième ordre en temps et en espace. Le pas spatial du maillage est \( \Delta x\) et son pas temporel \( \Delta t\). L'indice spatial est j, variant entre 0 et Nx, et l'indice temporel n variant entre 0 et Nt, de sorte que je puisse écrire \( u_j^n = u(x_j, t^n) \).
Pour le membre de gauche, j'obtiens, en appliquant le schéma centré d'Euler explicite :
\( \frac{\partial^2 u(x_j,t_n)}{\partial t^2} \approx \frac{u_j^{n+1} - 2 u_j^n + u_j^{n-1}}{\Delta t^2} \)
De même, pour le membre de droite :
\( \frac{\partial^2 u(x_j,t_n)}{\partial x^2} \approx \frac{u_{j+1}^n - 2 u_j^n + u_{j-1}^n}{\Delta x^2} \)
J'obtiens donc la forme discrétisée de l'équation de D'Alembert, sous forme explicite, et en posant \( \alpha = c\frac{\Delta t}{\Delta x}\) :
\( u_j^{n+1} = 2 u_j^n - u_j^{n-1} + \alpha^2 (u_{j+1}^n - 2 u_j^n + u_{j-1}^n) \)
Cette équation discrète peut s'écrire sous forme matricielle \( U^{n+1} = \alpha * A*U^n + 2*U^n - U^{n-1} \) avec A la matrice tridiagonale de coefficients 1,-2,1. Nous résoudrons le problème sous sa forme matricielle, bien adaptée à Python.
Le nombre \( \alpha = c\frac{\Delta t}{\Delta x} \) est important s'agissant de la stabilité de notre schéma numérique explicite. Il définit en effet le nombre CFL (Courant-Friedrich-Levy) qui doit respecter certaines règles pour garantir la stabilité du schéma.
En particulier, il doit être impérativement inférieur ou égal à 1, en dimension 1, ce qui impose des contraintes sur le choix du pas spatial et du pas temporel, la valeur de c étant par définition fixe.
Cette contrainte peut être interprétée physiquement. En effet, on conçoit que le pas d'espace \( \Delta x \) doive être supérieur à la distance parcourue par l'onde de célérité c pendant un temps \( \Delta t \), ce que nous impose la condition \( c\frac{\Delta t}{\Delta x} \le 1 \) ou encore \( c\Delta t \le \Delta x \) .
Je ne vais pas détailler la définition des maillages spatial et temporel, ni la séquence d'affichage, qui sont très classiques. Attention toutefois au respect de la condition CFL dans le calcul du pas temporel. Dans le code, j'ai choisi CFL = 1 et donc dt = dx/c. Mais voyons les quelques points particuliers de ce script.
Elle fait l'objet de la fonction MatDiag:
def MatDiag(n):
M = zeros((n+1,n+1),float)
for i in range(1,n):
M[i,i-1], M[i,i], M[i,i+1] = 1,-2,1 # coefficients diagonaux
M = Alpha*M
return M
Le paramètre Alpha est décrit dans le corps principal du script:
Alpha = (c*dt/dx)**2
Elles nous permettent de calculer les fonctions g(x,0) et h(x,0) du problème, en notant que h(x,0) est la dérivée de g(x,0). Pour g(x,0), j'ai choisi une fonction sinus, sachant que ce choix est arbitraire. Il est assez courant de choisir dans ce cas une fonction gaussienne.
Voici la fonction g(x,0), qui retourne le vecteur U0 qui permettra d'initialiser le schéma pour le temps t = 0
def Uzero():
U0 = sin(2*pi*X/4.0)
return U0
Puis la fonction dérivée h(x,0), qui retourne le vecteur U1, qui permettra d'initialiser le schéma pour le temps t = 1. Notez que je calcule la dérivée à partir du schéma aux DF défini plus haut. Notez aussi que je définis les conditions de bord U1[0] = U1[Nx] = 0.0 :
def dUzero():
U1 = zeros(Nx+1)
U0 = Uzero()
for i in range(1,Nx):
U1[i] = U0[i] + (Alpha/2.)*(U0[i+1] - 2*U0[i] + U0[i-1])
U1[0] = U1[Nx] = 0.0
return U1
Les intructions de construction de la matrice tridiagonale A et d'initialisation des conditions initiales :
A = MatDiag(Nx)
U[:,0] = Uzero()
U[:,1] = dUzero()
puis finalement le calcul du système matriciel
for t in range(1,Nt):
U[:,t+1] = dot(A,U[:,t]) + 2*U[:,t] - U[:,t-1]
Sur le domaine considéré, le script nous fournit l'évolution de l'onde en fonction du temps suivante:
Cette simulation d'onde se prête particulièrement à l'animation. Sur cet exemple, je n'ai pas voulu me lancer dans l'utilisation du module Visual de VPython, on verra ça plus tard. Il existe heureusement dans matplotlib un module permettant de faire quelque chose de sympa. Donc, au lieu de tracer une surface ondulante mais statique pour représenter l'évolution de notre onde dans le temps, j'ai modifié le script précédent pour tracer dans un même référentiel l'évolution de la courbe.
Les modifications apportées au script précédent ne concernent pas le calcul mais l'affichage. Tout d'abord, il faut importer la librairie qui va bien :
import matplotlib.animation as animation
puis remplacer les instructions d'affichage habituelles par celles-ci:
fig = plt.figure(figsize=(14,8))
ax = fig.add_subplot(111)
ax.set_xlabel('abscisse', fontsize = 16)
ax.set_ylabel('amplitude', fontsize = 16)
plt.ylim(-1.0,1.0)
# préparation de l'animation
lines = []
plt.hold("off")
for i in range(Nt):
p, = plt.plot(X,U[:,i],'k')
lines.append([p])
ani = animation.ArtistAnimation(fig,lines,interval=10,blit=True)
Le résultat est assez amusant...
Les scripts Python étudiés dans cette page sont disponibles dans le package OEMAlembert.zip :
Contenu et design par Dominique Lefebvre - tangenteX.com novembre 2016 -- Pour me joindre par ou PhysiqueX ou 
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.