Code matlab

MarcStephane - Modifié le 18 sept. 2023 à 10:47

Bonjour à vous tous et je souhaite à tous santé dans ces moments difficile
j'ai fait un programme dans lequel je simule sur matlab la propagation d'onde dans un milieu contenant des defauts en 2D;
le résultat me plait bien, mais je voudrais que quelqu'un m'aide pour limiter la dispersion (je jaune) qui apparait lors de la simulation

 
close<span> all </span><span>% Pour bien commencer...</span>
tic <span>% Si on veut mesurer par la suite le temps d'exécution (avec alors toc)</span>
c=10; <span>% célérité de l'onde</span>
 
<span>% Domaine et grille de calcul</span>
Lx=100;Ly=100;
T=10;
 
dx=Lx/200;
dy=Ly/200;
dt=sqrt<span>(</span>dx^2+dy^2<span>)</span>/<span>(</span>2*c<span>)</span>;
 
gax=<span>(</span>c*dt<span>)</span>^2/<span>(</span>dx^2<span>)</span>;
gay=<span>(</span>c*dt<span>)</span>^2/<span>(</span>dy^2<span>)</span>;
 
<span>% Grille de calcul et nombre de points</span>
xx=0:dx:Lx ;nx=length<span>(</span>xx<span>)</span>;
yy=0:dy:Ly ;ny=length<span>(</span>yy<span>)</span>;
tt=0:dt:T ;nt=length<span>(</span>tt<span>)</span>;
<span>%     F(i,j+1)=-F(i,j-1)+2*F(i,j)+(F(i-1,j)+F(i+1,j)-2*F(i,j))*(2*b+3*d*(F(i+1,j)-F(i-1,j))/deltax)*((F(i+1,j)-F(i-1,j))/deltax-1);</span>
<span>% Fréquence: 2Hz, de la source</span>
nu=2;
<span>% Pour initialiser une matrice en Matlab, on la remplit de zeros</span>
u=zeros<span>(</span>nx,ny,nt<span>)</span>;
 
<span>for</span> i=2:nx-1
    <span>for</span> j=2:ny-1
        u<span>(</span>i,j,2<span>)</span>=u<span>(</span>i,j,1<span>)</span>+<span>...</span>
            gax/2*<span>(</span>u<span>(</span>i-1,j,1<span>)</span>+u<span>(</span>i+1,j,1<span>)</span>-2*u<span>(</span>i,j,1<span>)</span><span>)</span>+<span>...</span>
            gay/2*<span>(</span>u<span>(</span>i,j-1,1<span>)</span>+u<span>(</span>i,j+1,1<span>)</span>-2*u<span>(</span>i,j,1<span>)</span><span>)</span>;
<span>    end</span>
<span>end</span>
<span>for</span> k=3:nt-1
    <span>for</span> i=2:nx-1
        <span>for</span> j=2:ny-1
            u<span>(</span>i,j,k+1<span>)</span>=2*u<span>(</span>i,j,k<span>)</span>-u<span>(</span>i,j,k-1<span>)</span>+<span>...</span>
                gax*<span>(</span>u<span>(</span>i-1,j,k<span>)</span>+u<span>(</span>i+1,j,k<span>)</span>-2*u<span>(</span>i,j,k<span>)</span><span>)</span>+<span>...</span>
                gay*<span>(</span>u<span>(</span>i,j-1,k<span>)</span>+u<span>(</span>i,j+1,k<span>)</span>-2*u<span>(</span>i,j,k<span>)</span><span>)</span>;
<span>        end</span>
<span>    end</span>
    <span>% et on impose la source à chaque instant:</span>
      sx=nx/3; 
      sy=ny;
      <span>% de calcul de la source</span>
    u<span>(</span>sx,sy,k+1<span>)</span>= u<span>(</span>sx,sy,k<span>)</span>+ sin<span>(</span>2*pi*nu*k*dt<span>)</span>;
   u<span>(</span>50,160,k<span>)</span>=0;
<span>% %    u(100,70,k)=0;</span>
<span>%     u(160,60,k)=0;</span>
<span>% et on impose aussi ici les conditions aux limites: </span>
      <span>% conditions sur les bords du domaine </span>
      <span>% et sur l'éventelle obstacle </span>
      <span>% (fente simple ou double dans le TP)</span>
      u<span>(</span>1,:,k+1<span>)</span>=u<span>(</span>2,:,k<span>)</span>;
u<span>(</span>13,:,k+1<span>)</span>=u<span>(</span>4,:,k<span>)</span>;
u<span>(</span>:,1,k+1<span>)</span>=u<span>(</span>:,2,k<span>)</span>;
u<span>(</span>:,6,k+1<span>)</span>=u<span>(</span>:,17,k<span>)</span>;
<span>end</span>
figure<span>(</span>1<span>)</span>;clf;whitebg<span>(</span><span>'w'</span><span>)</span>;
<span>for</span> k=1:nt
   pcolor<span>(</span>u<span>(</span>:,:,k<span>)</span><span>)</span>;
    <span>% Les axes n'ont ici aucun sens physique, on les retire</span>
    axis<span> on</span>;
    <span>% On fixe l'échelle des couleurs</span>
    caxis<span>(</span><span>[</span>-0.1 0.8<span>]</span><span>)</span>;
    <span>% On lisse la représentation par interpolation des données</span>
    shading interp;
    <span>% Et on laisse un peu de temps entre chaque affichage</span>
    pause<span>(</span>0.1<span>)</span>
<span>end</span>
A voir également: