90 lines
1.8 KiB
Matlab
Executable file
90 lines
1.8 KiB
Matlab
Executable file
function advection(scheme,Nt,Nx)
|
|
% Script calculant une approximation de la solution du problème
|
|
% d'advection linéaire 1D.
|
|
%
|
|
% Inputs
|
|
% ------
|
|
%
|
|
% scheme : schéma numérique à utiliser.
|
|
%
|
|
% Nt : nombre de pas de temps.
|
|
%
|
|
% Nx: nombre de pas d'espace.
|
|
%
|
|
% Exemple: advection('LaxWendroff',200,500);
|
|
|
|
|
|
% Cadre experimental
|
|
a=1.; % vitesse
|
|
L=3; % longueur du domaine spatial.
|
|
T=1; % longueur de la fenêtre temporelle.
|
|
ic=0; % condition initiale : 0 -> porte, sinon densité gaussienne
|
|
|
|
% Définition de la grille
|
|
dx=L/(Nx+1);
|
|
dt=T/(Nt+1);
|
|
xx=0:Nx+1;
|
|
xx=xx';
|
|
%dt=dx/a;
|
|
|
|
% Nombre de Courant
|
|
lambda=a*dt/dx
|
|
|
|
% Condition initiale
|
|
u0=zeros(Nx+2,1);
|
|
u0=reference(ic,lambda,Nx,dx,0);
|
|
rmse=zeros(Nt+1,1);
|
|
|
|
% Boucle temporelle
|
|
uh=u0;
|
|
uh_last = uh;
|
|
|
|
% calcul de A en fonction du scheme
|
|
e = ones(Nx+2, 1);
|
|
switch scheme
|
|
case 'explicite'
|
|
A = spdiags([(1-lambda)*e lambda*e], [0 -1], Nx+2, Nx+2);
|
|
case 'implicite'
|
|
A = spdiags([-lambda/2*e e lambda/2*e], [-1 0 1], Nx+2, Nx+2);
|
|
case 'LaxWendroff'
|
|
A = spdiags([(lambda/2+lambda^2/2)*e (1-lambda^2)*e (-lambda/2+lambda^2/2)*e], [-1 0 1], Nx+2, Nx+2);
|
|
end
|
|
|
|
for k=1:Nt+1
|
|
|
|
% intérieur du domaine
|
|
switch scheme
|
|
case 'explicite'
|
|
uh = A*uh;
|
|
case 'implicite'
|
|
uh = A\uh_last;
|
|
uh_last = uh;
|
|
case 'LaxWendroff'
|
|
uh = A*uh;
|
|
end
|
|
|
|
% Conditions aux limites
|
|
uh(1) = 0;
|
|
uh(length(uh)) = 0;
|
|
|
|
%Erreur RMS
|
|
uref=reference(ic,lambda,Nx,dx,k);
|
|
rmse(k)=norm(uh-uref,2)/norm(uref,2);
|
|
|
|
% Affichage de la solution
|
|
figure(1)
|
|
plot(dx*xx,uh,'b-',dx*xx,uref,'r-');
|
|
axis([0 L -1 max(abs(u0))+1]);
|
|
legend('Solution numerique','Solution de reference');
|
|
xlabel('Domaine spatial')
|
|
ylabel('u')
|
|
pause(0.1);
|
|
end
|
|
|
|
% Affichage de l'erreur RMS
|
|
figure(2)
|
|
plot(rmse);
|
|
legend('Erreur RMS')
|
|
|
|
end
|