Programmation en Matlab

Simulation de la propagation d'une onde 1D

Vibration d'une corde

Quelques éléments de programmation pour le TP 8: Simulation de la propagation d'une onde monodimensionnelle.

Initialisation: paramètres et constantes physiques

clear all;close all % Pour bien commencer...
tic % Si on veut mesurer par la suite le temps d'exécution (avec toc)
c=1;% célérité de l'onde

deltax=0.001;
deltat=0.001;
ga=(c*deltat/deltax)^2;

L=1; % Longueur (de la corde par exemple)
T=1; % Temps max de la simulation

% Grille de calcul, spaciale et emporelle
xx=[0:deltax:L];nx=length(xx);
tt=[0:deltat:T];nt=length(tt);

Initialisation: position et vitesse

Pour correctement initialiser notre onde, on a besoin des deux conditions au bord du système (2): position initiale via la fonction h, et vitesse initiale via la fonction h_tilde (ht dans le code).
Courbe représentant une Impulsion

x0=L/2; % Centre de l'impulsion
r=L/5;  % "rayon" de l'impulsion

h=exp(-1./abs((xx-x0).^2-r^2));
H=zeros(size(xx));H(round(nx/3):round(2*nx/3))=1; % filtre pour imposer h(x)=0 si |x-x0|>r
h=h.*H; % On applique le filtre à h
h=h/max(h); % enfin on normalise h (ainsi max(h)=1)

plot(xx,h); % s'il nous prend l'envie de tracer h...

Par la suite on peut bien sûr modifier la fonction h: avec 2 impulsions h=h1+h2, ou encore un signal rectangulaire (simulation de l'attaque d'une corde de piano), ...

En ce qui concerne la vitesse initiale proposée dans le TP, c'est plus simple: vitesse initiale nulle:

ht=zeros(1,nx); % aucune vitesse initiale

vitesse initiale qu'on pourra bien sûr aussi s'ameser à modifier à loisir ensuite…

Initialisation de la matrice

Il s'agit de la discrétisation des conditions au bord de (2):

%% Initialisation de la matrice:
% Pour initialiser une matrice en Matlab, on la remplit de zeros
F=zeros(nx,nt);
%% Initialisation temporelle
% Position initiale (donnée par la fonction impulsion h)
F(:,1)=h(:);
% Vitesse initiale (donnée par la fonction ht)
F(:,2)=F(:,1)+deltat*ht(:);
% spatiale
F(1,:)=0;
F(nx,:)=0;

Equation des ondes discrétisée

Enfin, on s'attaque à la propagation à proprement parler: la discrétisation (3) de l'équation des ondes (1):

for j=2:nt-1
  for i=2:nx-1
    F(i,j+1)=-F(i,j-1)+ga*F(i+1,j)+ga*F(i-1,j)-2*(1-ga)*F(i,j);
  end
end

ou bien mieux en Matlab (comparer les deux méthodes et leur temps de calcul avec tic et toc):

for j=2:nt-1
  F(2:nx-1,j+1)=-F(2:nx-1,j-1)+ga*F(3:nx,j)+ga*F(1:nx-2,j)+2*(1-ga)*F(2:nx-1,j);
end

Affichage

On peut tout afficher en une seule fois, toutes les valeurs (amplitudes et en fonction du temps) étant dans la matrice F.
L'affichage du contenu d'une matrice peut se faire par exemple avec pcolor:

figure % Nouvelle figure
pcolor(F);
shading interp; % On lisse la représentation par interpolation des données

On peut aussi représenter l'évolution de l'onde en fonction du temps, tel un film, avec une succession des représentations des positions à chaque instant:

figure % Nouvelle figure
for j=1:nt
  plot(xx,F(:,j));
  % et on fixe l'échelle: identique à chaque tracer !
  axis([0 1 1.2*min(min(F)) 1.2*max(max(F))]);
  pause(0.01) % et on se laisse un peu de temps pour voir
end