Source Latex
\documentclass[12pt]{article}
%\usepackage{french}
\usepackage[french]{babel}
\usepackage{amsmath}
\usepackage{amsfonts}\usepackage{amssymb}
\usepackage[utf8]{inputenc}
\usepackage{a4wide}
\usepackage{graphicx}
\usepackage{epsf}
\usepackage{pst-all}
\usepackage{ifthen}
\usepackage{calc}
\usepackage{pst-all}
\usepackage{hyperref}
\hypersetup{
pdfauthor={Yoann Morel},
pdfsubject={Matlab},
pdftitle={Simulation monodimensionnelle de la propagation d'une onde},
pdfkeywords={onde, propagation, équation des ondes, 1D, monodimensionnel, Matlab, simulation, programmation, modélisation, mathématiques}
}
\hypersetup{
colorlinks = true,
linkcolor = red,
anchorcolor = red,
citecolor = blue,
filecolor = red,
urlcolor = red
}
\voffset=-1.cm
% Raccourcis diverses:
\newcommand{\nwc}{\newcommand}
\nwc{\dsp}{\displaystyle}
\nwc{\bge}{\begin{equation}}\nwc{\ene}{\end{equation}}
\nwc{\bgar}{\begin{array}}\nwc{\enar}{\end{array}}
\nwc{\bgit}{\begin{itemize}}\nwc{\enit}{\end{itemize}}
\nwc{\bgmp}{\begin{minipage}}\nwc{\enmp}{\end{minipage}}
\nwc{\ct}{\centerline}
\nwc{\V}[1]{\overrightarrow{#1}}
\nwc{\la}{\left\{}\nwc{\ra}{\right\}}
\nwc{\lp}{\left(}\nwc{\rp}{\right)}
\nwc{\lb}{\left[}\nwc{\rb}{\right]}
\nwc{\bgsk}{\bigskip}
\nwc{\vsp}{\vspace{0.1cm}}
\nwc{\vspd}{\vspace{0.2cm}}
\nwc{\vspt}{\vspace{0.3cm}}
\nwc{\vspq}{\vspace{0.4cm}}
\newcounter{ntp}%[section]
\setcounter{ntp}{1}
\newenvironment{TP}[1]{%
\setcounter{equation}{0}\setcounter{subsection}{0}
\section*{ \underline{TP \arabic{ntp} :} #1}
\addcontentsline{toc}{section}{\underline{TP \arabic{ntp} :}\ #1}
\stepcounter{ntp}}{}
\nwc{\bgTP}{\begin{TP}}\nwc{\enTP}{\end{TP}}
\def\N{{\rm I\kern-.1567em N}} % Doppel-N
\def\No{\N_0} % Doppel-N unten 0
\def\R{{\rm I\kern-.1567em R}} % Doppel R
\def\C{{\rm C\kern-4.7pt % Doppel C
\vrule height 7.7pt width 0.4pt depth -0.5pt \phantom {.}}}
\def\Z{{\sf Z\kern-4.5pt Z}} % Doppel Z
% grec et autres
\nwc{\lbd}{\lambda}
\nwc{\vphi}{\varphi}
\nwc{\ga}{\gamma}
\def\epsi{\varepsilon}
\def\tht{\theta}
\nwc{\tm}{\times}
\nwc{\ul}{\underline}
\nwc{\scp}[1]{\scriptsize{#1}}
\newcommand{\TITLE}{Simulation monodimensionnelle de la propagation d'une onde}
\newcommand{\TITLEfoot}{Simulation de la propagation d'une onde 1D}
\author{Y. Morel}
\date{}
\headheight=0cm%15pt
\textheight=26cm
\topmargin=-1.8cm
\footskip=.8cm
\textwidth=17.6cm
\oddsidemargin=-1.cm
%\parindent=0.2cm
%\pagestyle{headings}
\usepackage{fancyhdr}
\usepackage{lastpage}
\newenvironment{centerpage}{\vspace*{\fill}}{
\protect\vspace*{\fill}}
\begin{document}
\pagestyle{fancyplain}
\setlength{\headheight}{0cm}
\renewcommand{\headrulewidth}{0pt}
\renewcommand{\footrulewidth}{0.1pt}
\lhead{}\chead{}\rhead{}
\lfoot{Y. Morel \ \url{https://xymaths.fr/Matlab/}}
\rfoot{\TITLEfoot\ - \thepage/\pageref{LastPage}}
\cfoot{}%\TITLEfoot}
\vspace*{1cm}
\nwc{\TPtop}{}
\newlength{\lgr}
\newlength{\TPtopLarge}
\newlength{\EXX}
\newlength{\Tlgr}
\newlength{\TopHeight}
\newlength{\TopTitleHeight}
\newcommand{\ChapTitle}[3]{
\renewcommand{\TPtop}{{\LARGE\bf #1} {\Huge\bf#2}}
\settowidth{\lgr}{\LARGE\bf #1}
\settowidth{\TPtopLarge}{\TPtop}
\settoheight{\EXX}{\LARGE X}
\settowidth{\Tlgr}{\LARGE T}
\setlength{\Tlgr}{0.5\Tlgr}
\settototalheight{\TopHeight}{\fbox{\bgmp{\linewidth}\begin{center}{\LARGE\bf#3}\end{center}\enmp}}
%%
\noindent
\bgmp{\linewidth}
\psline(\lgr,\EXX)(1.1\lgr,\EXX)
\psline(\TPtopLarge,\EXX)(\linewidth,\EXX)%
(\linewidth,-1.\TopHeight)(\Tlgr,-1.\TopHeight)%
(\Tlgr,-0.05)
\TPtop
\vspace{-0.5\EXX}
\begin{center}
{\LARGE\bf#3}
\end{center}
\enmp%
\vspace{0.5cm}
}
\setcounter{ntp}{8}
\ChapTitle{TP}{\thentp}{Simulation monodimensionnelle\linebreak de la propagation d'une onde}
\addcontentsline{toc}{section}{TP \thentp: Simulation monodimensionnelle de la propagation d'une onde}
La propagation d'une onde, acoustique ou électromagnétique, est régie
par l'équation de Helmohltz, pour $t>0$ et $x\in D$, le domaine
géométrique étudié:
\begin{equation}\label{Helmohltz}
\frac{\partial^2 f}{\partial t^2}(x,t) - c^2\frac{\partial^2
f}{\partial x^2}(x,t) =0\,,
\end{equation}
où $c$ désigne la célérité de l'onde dans le milieu considéré
($c=3.10^8 \ \textrm{m.s}^{-1}$ pour la lumière dans le vide, par
exemple), et la fonction $f$ désigne l'amplitude de l'onde au point
$x$ et à l'instant $t$.
\bigskip
Cette équation n'est pas suffisante pour décrire complètement la
propagation de l'onde; il faut la compléter de conditions initiales
(état et vitesse de l'onde à l'instant initial $t=0$), ainsi que des
conditions au bord (comportement de l'onde lorsqu'elle rencontre les
limites physiques du domaine étudié).
Ces conditions peuvent se formuler suivant:
\begin{equation}\label{Donnees}\left\{
\begin{array}{ll}
f(x,0)&=h(x) \ \,,\ \forall x\in D \\
\dsp \frac{\partial f}{\partial t}(x,0)&=\tilde{h}(x) \ \,,\ \forall x\in D \\
f(x,t)&=g(t) \ \,,\ \forall t>0\ \, \mbox{ et, } \forall x\in\partial D\,,
\end{array}\right.
\end{equation}
les fonctions $h$, $\tilde{h}$ et $g$ étant des données du problème.
\bigskip
Le but de cette partie est de résoudre numériquement sous Matlab la
propagation d'une onde, i.e. de résoudre numériquement le jeu
d'équations (\ref{Helmohltz}) et (\ref{Donnees}), dans le cas simple
1-D, c'est-à-dire lorsque le domaine $D$ est un intervalle : $D=]a,b[$.
\vsp
La cas 1-D correspond par exemple à la simulation du mouvement d'une
corde (de piano ou de guitare par exemple), ou de l'amplitude d'un
champ électrique ou magnétique en polarisation TE ou TM.
\subsection*{Discrétisation de l'équation de Helmohltz (\ref{Helmohltz})}
\addcontentsline{toc}{subsection}{Discrétisation de l'équation des ondes}
On cherche une solution approchée de (\ref{Helmohltz}) en 1-D:
le domaine géométrique dans lequel l'onde peut se propager est
le segment $[0\,,L]$.
On s'interesse alors à l'onde $f(x,t)$ pour $x\in[0\,,L]$ et $t\in
[0\,,T]$.
\paragraph{Discrétisation de (\ref{Helmohltz}): } Pour approximer
l'onde solution de (\ref{Helmohltz}), nous allons utilisé la méthode
numérique dite des différences finies.
On se donne tout d'abord pour cela une ``grille'' régulière de calcul,
c'est-à-dire une discrétisation de l'ensemble $[0\,,L]\times[0\,,T]$:
soit les deux suites
\[ \{x_k\} =\{ k\,\delta_x\}_{k=0\dots L/\delta_x}\ \,, \mbox{ et, }
\{t_k\} =\{ k\,\delta_t\}_{k=0\dots T/\delta_t} \,,
\]
où $\delta_x$ et $\delta_t$ désignent respectivement le pas de la
grille en $x$ (espace) et en $t$ (temps).
\bigskip
Cette grille étant donnée, on cherche à calculer les valeurs de $f$
aux noeuds de la grille. On notera
\[ f_{i,j} := f(x_i,t_j)\,.
\]
La deuxième étape est l'approximation des dérivées partielles de $f$ à
l'aide des éléments $f_{i,j}$. Cela s'effectue à l'aide de la formule
de Taylor, par exemple pour la dérivée selon $x$:
\[ f(x+\delta_x,t) = f(x,t)+\delta_x \frac{\partial f}{\partial
x}(x,t) + \frac{\delta_x^2}{2} \frac{\partial^2 f}{\partial
x^2}(x,t) + O (\delta_x^3) \,,
\]
ou, en d'autres termes
\[ f_{i+1,j} = f_{i,j} + \delta_x \frac{\partial f_{i,j}}{\partial x}
+ \frac{\delta_x^2}{2} \frac{\partial^2 f_{i,j}}{\partial
x^2} + O (\delta_x^3) \,.
\]
On peut de cette manière approximer les dérivées secondes de $f$ par
rapport à $x$ et $t$ (en considérant aussi par exemple le
développement de Taylor de $f(x-\delta_x,t)$ et en sommant les
relations).
On obtient, tous calculs fait:
\[ \frac{\partial^2 f_{i,j}}{\partial t^2} \sim
\frac{1}{\delta_t^2}\left[
f_{i,j+1} + f_{i,j-1} - 2 f_{i,j}\right] \,,
\]
et,
\[\frac{\partial^2 f_{i,j}}{\partial x^2} \sim
\frac{1}{\delta_x^2}\left[
f_{i+1,j} + f_{i-1,j} - 2 f_{i,j}\right] \,.
\]
En combinant ces deux expressions avec l'équation de Helmohltz
(\ref{Helmohltz}), on obtient la relation de récurrence:
\begin{equation}\label{rec}
f_{i,j+1} = - f_{i,j-1} + \gamma f_{i+1,j} +\gamma f_{i-1,j}
-2 (1-\gamma) f_{i,j} \,,
\end{equation}
où $\gamma = c^2 \delta_t^2/\delta_x^2$.
\bigskip
Il ne reste plus alors qu'à initialiser la suite $f_{i,j}$ au moyen
des équations (\ref{Donnees}).
Pour commencer, on pourra prendre $g=0$. Cette condition modélise une
paroi infiniment dure; l'onde y sera totalement réfléchie.
La fonction $h(x)$ modélise l'état initial de l'onde, à l'instant
$t=0$.
On pourra utiliser la fonction $h$, infiniment régulière:
\[ h(x) = \left\{
\begin{array}{ll}
\displaystyle A \exp\left( -\frac{1}{|(x-x_0)^2-r^2|} \right) \ , \mbox{ pour } |x-x_0|\leq r \\
0 \ , \mbox{ ailleurs}\,,
\end{array}\right.
\]
où $A=\exp(1/r^2)$ permet de normaliser $h$, et où $x_0\in [0,L]$ est
le centre de la ``cloche'' $h(x)$ et $r>0$ sa largeur.
La fonction $\tilde{h}$ modélise l'''impulsion'', ou la vitesse,
initiale donnée à l'onde. On pourra choisir pour commencer
$\tilde{h}=0$.
\subsection*{\underline{\bf Travail à réaliser}}
\paragraph{Initialisation}
On prendra (pour commencer) les valeurs numériques suivantes:
$c=1\textrm{m.s}^{-1}$, $\delta_x=\delta_t=10^{-3} \textrm{m}$, soit
$\gamma =1$, et $L=1 \textrm{m}$, $T=1 \textrm{s}$.
\bigskip
On définit ainsi la grille de calcul:
$ (x_i,t_j)\in [0:\delta_x:L]\times[0:\delta_t:T]$.
La solution $f(x,t)$, approximée aux point $(x_i,t_j)$, sera alors
stockée dans une matrice rectangulaire $F_{i,j}$.
\bigskip
On utilisera la fonction $h(x)$ définie précédemment pour définir
l'état initial de l'onde, avec $x_0=L/2$ et $r=L/5$, tandis que l'on
imposera $\tilde{h}(x)=0$ (vitesse initiale nulle).
L'initialisation s'implémente alors simplement suivant:
\[ F(:,1)=h(:) \ , \mbox{ et, } F(:,2)=F(:,1)+\delta_t h_t(:) \,,
\]
où $h_t(x)$ est la vitesse initiale de l'onde (on pourra prendre pour
commencer $h_t=0$).
\bigskip
La condition au bord: $f(x=0,t)=f(x=L,t)\,,\forall t$ se traduit quant
à elle par:
\[ F(1,:)=0 \ ,\mbox{ et, } F(end,:)=0 \,.
\]
\paragraph{Récurrence}
Il s'agit ici d'implémenter la relation de récurrence (\ref{rec}).
Le calcul peut se faire à l'aide de deux boucles imbriquées: à chaque
pas de temps (i.e. pour chaque $j$), on calcul tout les $F_{i,j}$
récursivement.
\paragraph{Visualisation du résultat}
Le comportement complet de l'onde est contenu dans la matrice $F$:
chaque colonne contient l'amplitude de l'onde à un instant donné,
chaque ligne correspondant à une abscisse.
On peut simplement visualiser la propagation de l'onde en utilisant la
commande {\it imagesc(F)}.
L'amplitude de l'onde au $\textrm{n}^{\mbox{ième}}$ pas de temps
s'affiche alors simplement par la commande
{\it plot($[0:\delta_x:L],F(:,n)$)}.
\bigskip
Comme il s'agit d'un problème d'évolution, il peut aussi être
interessant de représenter ces tracés successivement.
On pourra, par exemple, regrouper ces graphiques successifs dans une
video à l'aide des commandes {\it getframe}, et {\it avifile}.
\paragraph{Suite de l'étude...}
On peut reprendre l'étude précédente, et modifier les conditions
initiales et la condition au bord.
On pourra, par exemple, s'interesser au cas ou deux ondes existent
initialement, et qui vont donner lieu à des phénomènes d'interférences
(constructives ou destructives).
\bigskip
La condition au bord pourra elle aussi être modifiée en la condition
dite de Neumann:
\[ \frac{\partial f}{\partial x}(x,t)=0 \ , \mbox{ pour }
x=0\ \mbox{ et } x=L \,.
\]
\label{LastPage}
\end{document}
Télécharger le fichier source 
