Wave propagation and interference phenomenom simulation with Matlab

Young's double-slit experiment


Résultats, animation et observations

L'animation ci-dessous montre les résultats de la simulation détaillée par la suite:



La fenêtre graphique est divisée en 3 figures: (voir mise en forme avancée et subplots),
  • à gauche: la simulation complète, propagation de l'onde, dont le passage au travers des deux fentes
  • bas à droite: onde se propageant après le passage dans la double fente
  • haut droite: écran d'observation (en haut des deux autres graphiques), sur lequel on observe des figures d'interférence. La courbe en rouge est l'observation théorique (voir plus bas)

Principes des calculs numériques et de la simulation - programme Matlab

Les calculs numériques de la simulation reposent sur la méthode des différences finies, voir par exemple cete page ou aussi introduction à la programmation et la simulation en Matlab, et plus précisément TP 9: Simulation bidimensionnelle de la propagation d'une onde, Diffraction par une simple ou double fente

Les calculs, simulations, et animations sont réalisés ici avec Matlab, à partir du programme: 
% Nouveau départ... propre...
clear all;close all
tic
c=10;

% Grille de calcul
Lx=100;
Ly=100;
T=15;

deltax=Lx/200;
deltay=Ly/200;
deltat=sqrt(deltax^2+deltay^2)/(2*c)
deltat=0.9*deltat

gax=(c*deltat)^2/(deltax^2);
gay=(c*deltat)^2/(deltay^2);

xx=[0:deltax:Lx];nx=length(xx);
yy=[0:deltay:Ly];ny=length(yy);
tt=[0:deltat:T];nt=length(tt);

Reflection_parois=0; % si 1, reflection sur les bords


Largeurx=1;Nlgx=round(Largeurx/Lx*nx);
posfx=50;nfx1=round(posfx/Lx*nx);nfx2=nfx1+Nlgx;

LargeurFy=2;NFy=round(LargeurFy/Ly*ny);
posF1y=Ly/2-LargeurFy/2;posF2y=Ly/2+LargeurFy/2;

cF1y=round(posF1y/Ly*ny);cF1y=96;
cF2y=round(posF2y/Ly*ny);cF2y=104;
nF1y1=cF1y-NFy/2;nF1y2=nF1y1+NFy/2;
nF2y1=cF2y-NFy/2;nF2y2=nF2y1+NFy/2;

Fentex=[nfx1:nfx2];
Fentey=[1:nF1y1,nF1y2:nF2y1,nF2y2:ny];

% Position de la source
sx=round(nx/4);sy=round((nF1y2+nF2y1)/2);

% Fréquence de la source: 2 Hertz
sfrq=2;

Mat=zeros(nx,ny,nt);
% Indices de Mat, hors bords:
i=[2:nx-1];j=[2:ny-1];

Mat(i,j,2)=Mat(i,j,1)+...
    gax/2*(Mat(i-1,j,1)+Mat(i+1,j,1)-2*Mat(i,j,1))+...
    gay/2*(Mat(i,j-1,1)+Mat(i,j+1,1)-2*Mat(i,j,1));


for k=3:nt-1
     Mat(i,j,k+1)=2*Mat(i,j,k)-Mat(i,j,k-1)+...
                gax*(Mat(i-1,j,k)+Mat(i+1,j,k)-2*Mat(i,j,k))+...
                gay*(Mat(i,j-1,k)+Mat(i,j+1,k)-2*Mat(i,j,k));
   
    Mat(sx,sy,k+1)=sin(2*pi*sfrq*k*deltat);
    
    if (Reflection_parois==0)
        Mat(1,:,k+1)=Mat(2,:,k);
        Mat(nx,:,k+1)=Mat(nx-1,:,k);
        Mat(:,1,k+1)=Mat(:,2,k);
        Mat(:,ny,k+1)=Mat(:,ny-1,k);
    end
    
    Mat(Fentex,Fentey,k+1)=0;
end
m=max(max(max(Mat)))
% Affichage de la double fente
Mat(Fentex,Fentey,:)=1;

% 3 zones d'affichage:
%  3 subplots: 
%    - gauche: whole simulation: wave propgation + double-slit
%    - bas droite: onde après les fentes
%    - haut droite: écran, montrant les interférences
fig=figure(1);clf;whitebg('w');
set(fig,'position',[0 0 640 480]);

% Initialisation de la variable Film
%  qui contiendra toutes les images de l'animation
Film=[];

ecranx=200;
Interf=squeeze(Mat(ecranx,:,:));
MatEcran=Mat;
MatEcran(ecranx,:,:)=5;
MatEcranSup=MatEcran(101:end,:,:);

% Résultats théoriques:, see below 
lbd=c/sfrq; % longueur d'onde
k=2*pi/lbd; % nombre d'onde
D=Lx/2;
D=(ecranx-max(Fentex))/200*Lx
d=(nF2y1-nF1y2)/200*Ly;
w=k*d/D;


for kk=1:1:nt
    subplot('position',[0.05 0.2 0.45 0.6])
    pcolor(squeeze(MatEcran(:,:,kk)));
    shading interp,axis off
    caxis([-0.4 0.4])
    
    subplot('position',[0.55 0.05 0.4 0.4])
    pcolor(squeeze(MatEcranSup(:,:,kk)));
    shading interp
    caxis([-0.003 0.003]);
    axis off
    
    subplot('position',[0.55 0.55 0.4 0.4])
    axis off;hold on
    % La courbe théorique:
    plot(xx,(.75e-6)*(cos(w*(xx-49))).^2,'-r','linewidth',2);

    plot(xx,Interf(:,kk).^2)
    axis([xx(30) xx(end-30) 0 1e-6])
    title('{\bf \''Ecran}','fontsize',16,'interpreter','Latex')
    
    %Film=[Film getframe(fig)];
    pause(0.05)
end


%movie2avi(Film,'Film.avi','compression','none','fps',15)

Détails théoriques et mathématiques: phénomène d'interference



\[\psset{unit=1cm,arrowsize=5pt} \begin{pspicture}(-.8,-3)(4.6,3) \psline[linewidth=1.5pt](0,-3)(0,-.4) \psline[linewidth=1.5pt](0,-.2)(0,.2) \psline[linewidth=1.5pt](0,.4)(0,3) \psline[linecolor=blue]{<->}(-.4,-.3)(-.4,.3) \rput[r](-.6,0){\blue$d$} \psline[linecolor=blue]{<->}(0,-2)(4,-2) \rput[r](2,-2.3){\blue$D$} \psline[linewidth=1.5pt](4,-3)(4,3) \psline(0,-.3)(4,1.5)(0,.3) \rput(2,1.15){$d_1$}\rput(2,.3){$d_2$} \psline[linestyle=dashed](-.25,0)(4.5,0) \psline[linecolor=blue]{->}(4.2,0)(4.2,1.5) \rput[l](4.3,.8){\blue$x$} \end{pspicture}\]


L'amplitude $A(x)$ à un point de l'écran est la somme des amplitudes des deux ondes provenant des sources secondaires cohérentes (de même fréquence):
\[A(x)(x)=A_1(x)+A_2(x)=I_0e^{i(\omega t-kd_1)}+I_0e^{i(\omega t-kd_2)}\]

où $I_0$ est l'intensité lumineuse, $\omega$ est la pulsation, $t$ est le temps d'observation (par rapport à une origine temporelle quelconque, mais de toute façon sans intérêt pour la suite), $k$ est le nombre d'onde, et $d_1$ et $d_2$ les distances telles que représentées sur l'illustration précédente.
On peut aussi écrire cette somme selon
\[A(x)=I_0e^{i\omega t}\left( e^{-ikd_1}+e^{-ikd_2}\rp\]

et donc, en factorisant par l'angle moitié comme il l'est courant avec des exponentielles complexes:
\[\begin{array}{ll} A(x) &=I_0e^{i\omega t}e^{-ik(d_1+d_2)/2}\Bigl[ e^{ik(d_2-d_1)/2}+e^{ik(d_2-d_1)/2}\Bigr]\\[.8em] &=I_0e^{i\omega t}e^{-ik(d_1+d_2)/2}\times 2\cos\left( k\dfrac{d_2-d_1}{2}\right) \enar\]


Ensuite, l'intensité lumineuse est proportionnelle au carré du module de l'amplitude:
\[I(x)=|A(x)|^2 =4I_0^2\cos^2\left( k\dfrac{d_2-d_1}{2}\right) \]


$\Delta=d_2-d_1$ est appelée différence de phase; l'expression précédente montre ainsi que l'intensité dépend seulement de la différence de trjat parcouru and chaque onde.
Géométriquement, avec les notations précédentes et le théorème de Pythagore,
\[\begin{array}{ll} d_1&=\sqrt{\left( x-\dfrac{d}{2}\rp^2+D^2} \\[1em] d_2&=\sqrt{\left( x+\dfrac{d}{2}\rp^2+D^2} \enar\]


En supposant que la distance $D$ depuis les fentes jusqu'à l'écran d'observation est grande comparée à celle $d$ séparant les deux fentes, on a avec un développement limité au 1er ordre,
\[\begin{array}{ll} d_1&=\sqrt{x^2-dx+\dfrac{d^2}{4}+D^2}\\[1em] &=D\sqrt{1+\dfrac{x^2}{D^2}-\dfrac{d}{D^2}x+\dfrac{d^2}{4D^2}}\\[1em] &\simeq D\left( 1+\dfrac{x^2}{2D^2}-\dfrac{d}{2D^2}x\rp\\ &\simeq D+\dfrac{x^2}{2D}-\dfrac{d}{2D}x \enar\]

et de même pour $d_2$:
\[ d_2\simeq D+\dfrac{x^2}{2D}+\dfrac{d}{2D}x \]


Finallement, au 1er ordre,
\[d_2-d_1=\dfrac{d}{D}x\]

et alors, l'intensité sur l'écran est donnée par
\[I(x)=4I_0^2\cos^2\lp\dfrac{kd}{D}x\rp\]


et est donc proportionnelle à:
\[\psset{xunit=1cm,yunit=2cm,arrowsize=7pt} \begin{pspicture}(-5,-.3)(5.5,1.4) \psline[linewidth=1pt]{->}(-5,0)(5.5,0) \psline[linewidth=1pt]{->}(0,-.2)(0,1.4) \psplot[plotpoints=200,linecolor=red,linewidth=1.5pt]{-5}{5}{x 180 mul 3.14 div cos 2 exp} \end{pspicture}\]

Cette expression est celle codée dans la simulation Matlab et peut être observée dans l'animation en haut de la page (courbe sinusoïdale rouge).

Voir aussi:
Top Lien


LongPage: h2: 1 - h3: 3