ENSEEIHT/TP-equation-derivees-partielles
1
0
Fork 0
Code Issues Pull requests Packages Projects Releases Wiki Activity

init

Browse source
This commit is contained in:
Laureηt 2023-06-10 20:51:59 +02:00
commit b7661507e8
Signed by: Laurent
SSH key fingerprint: SHA256:kZEpW8cMJ54PDeCvOhzreNr4FSh6R13CMGH/POoO8DI
9 changed files with 409 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

13
TP1/Affichage/SetClear.m Executable file
View file

@ -0,0 +1,13 @@
global PosWait;
global SizeWait;
global WaitBar;
present = isempty(PosWait)*isempty(SizeWait)*isempty(WaitBar);
if (present==0)
clear Hfig KcFrame WaitFrame black fig fy4 lightgreen orange red whitepaper
clear KcBar PosKc SizeKc WaitLabel blue fx4 lightblue magenta pink white
close(10)
end
PosWait = ''; SizeWait = ''; WaitBar = '';
clear PosWait SizeWait WaitBar present

64
TP1/Affichage/SetCtrl.m Executable file
View file

@ -0,0 +1,64 @@
global PosWait;
global SizeWait;
global WaitBar;
global PosKc;
global SizeKc;
global KcBar;
fig = 10;
figure(fig), clf
Hfig = gcf;
set( Hfig, 'Position', [1000, 600, 380, 100], 'Resize', 'off' );
% Define some default colors :
orange = [0.8 0.7 0.1];
magenta = [1.0 0.0 1.0];
pink = [0.8 0.5 0.7];
blue = [0.2 0.2 1.0];
lightblue = [0.5 0.7 0.8];
lightgreen = [0.2 0.8 0.3];
white = [1.0 1.0 1.0];
red = [1.0 0.0 0.0];
whitepaper = [0.9 0.9 0.6];
black = [0.0 0.0 0.0];
fx4 = 35; fy4 = 35;
WaitLabel = uicontrol( Hfig, ...
'Style', 'text', ...
'Position', [fx4+15 fy4+20 280 15], ...
'BackgroundColor', lightblue, ...
'ForegroundColor', blue, ...
'Visible', 'on', ...
'FontSize', [10], ...
'FontWeight', 'normal', ...
'HorizontalAlignment', 'center', ...
'String', 'Wait please ...' );
WaitFrame = uicontrol( Hfig, ...
'Style', 'frame', ...
'Position', [fx4+15 fy4+5 280 15], ...
'Visible', 'on', ...
'BackgroundColor', white );
PosWait = [fx4+16 fy4+6 1 13];
SizeWait = 278;
WaitBar = uicontrol( Hfig, ...
'Style', 'frame', ...
'Position', PosWait, ...
'Visible', 'off', ...
'ForegroundColor', red, ...
'BackgroundColor', red );
KcFrame = uicontrol( Hfig, ...
'Style', 'frame', ...
'Position', [fx4+15 fy4-10 280 15], ...
'Visible', 'on', ...
'BackgroundColor', white );
PosKc = [fx4+16 fy4-9 1 13];
SizeKc = 278;
KcBar = uicontrol( Hfig, ...
'Style', 'frame', ...
'Position', PosKc, ...
'Visible', 'off', ...
'ForegroundColor', blue, ...
'BackgroundColor', blue );

28
TP1/Affichage/SetWait.m Executable file
View file

@ -0,0 +1,28 @@
function [] = SetWait( jvp, n, Kc, KcExp )
global PosWait;
global SizeWait;
global WaitBar;
global PosKc;
global SizeKc;
global KcBar;
present = isempty(PosWait)*isempty(SizeWait)*isempty(WaitBar);
if (present~=0)
return
end
PosWait(3) = max( fix( (jvp/n)*SizeWait ), 1 );
set( WaitBar, 'Position', PosWait, 'Visible', 'on' ); drawnow;
if (nargin > 2)
present = isempty(PosKc)*isempty(SizeKc)*isempty(KcBar);
if (present~=0)
return
end
PosKc(3) = max( fix( (Kc/KcExp)*SizeKc ), 1 );
set( KcBar, 'Position', PosKc, 'Visible', 'on' ); drawnow;
end

59
TP1/Affichage/plot_uh.m Executable file
View file

@ -0,0 +1,59 @@
function nbfig_out=plot_uh(uh,dx1,dx2,N1,N2,nbfig)
%
% Cette fonction affiche la solution de l'EDP
%
% Inputs
% ------
%
% uh : solution numérique de l'EDP.
%
% dx1 : pas d'espace dans la direction x1.
%
% dx2 : pas d'espace dans la direction x2.
%
% N1 : nombre de points de grille dans la direction x1.
%
% N2 : nombre de points de grilles dans la direction x2.
%
% nbfig : identifiant de la figure
%
% Outputs:
% -------
%
% nbfig_out : identifiant de la figure
%
% Construction de la grille pour l'affichage.
x1min=0.;
x1max=dx1*(N1+1);
x2min=0.;
x2max=dx2*(N2+1);
if(x1max==0)
disp('Erreur : L1=0')
return
else
rxy=x2max/x1max;
end
x1=x1min:dx1:x1max;
x2=x2min:dx2:x2max;
[X,Y]=meshgrid(x1,x2);
[m,n]=size(X);
% Mise au format 2D de la solution.
Z=zeros(N2+2,N1+2);
Z(2:N2+1,2:N1+1)=reshape(uh,N2,N1);
% Affichage
nbfig_out=nbfig;
figure(nbfig_out)
s=surf(X,Y,Z,'FaceAlpha',0.5);
s.EdgeColor='none';
pbaspect([1 rxy 1]);
colorbar;
xlabel('x_1')
ylabel('x_2')
zlabel('u')
end

45
TP1/diffusivity.m Executable file
View file

@ -0,0 +1,45 @@
function [uhref]=diffusivity(nu,L1,L2,N1,N2)
%
% Cette fonction résoud le problème du calcul d'une solution du Laplacien anisotrope.
%
% Inputs
% ------
%
% nu: [nu1;nu2] valeurs des paramètres de diffusivité.
%
% L1 : longeur du domaine dans la direction x1.
%
% L2 : longueur du domaine dans la direction x2.
%
% N1 : nombre de points de grille dans la direction x1.
%
% N2 : nombre de points de grille dans la direction x2.
%
%
% Outputs:
% -------
%
% uhref : vecteur de taille N1*N2 contenant une approximation de la solution
%
% Ajout du repertoire Affichage à l'environnement
addpath('Affichage');
% Construction de la grille
dx1 = L1/(N1+1);
dx2 = L2/(N2+1);
% Calcul de la solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcul de la matrice du systeme
L = laplacian(nu,dx1,dx2,N1,N2);
% Calcul de l'approximation de la solution de l'EDP
b = forcing(nu, dx1, dx2, N1, N2);
uhref = L\b;
%Affichage de la solution approximee
fig_ref = plot_uh(uhref, dx1, dx2, N1, N2, 2); drawnow;
end

27
TP1/forcing.m Executable file
View file

@ -0,0 +1,27 @@
function b=forcing(nu, dx1, dx2, N1, N2)
%
% Cette fonction construit le vecteur de forçage de l'EDP
%
% Inputs
% ------
%
% nu : nu=[nu1;nu2], coefficients de diffusivité dans les directions x1 et x2.
%
% dx1 : pas d'espace dans la direction x1.
%
% dx2 : pas d'espace dans la direction x2.
%
% N1 : nombre de points de grille dans la direction x1.
%
% N2 : nombre de points de grilles dans la direction x2.
%
% Outputs:
% -------
%
% b : vecteur de forçage (dimension N1N2)
%
%
b = -ones(N1*N2, 1);
end

39
TP1/laplacian.m Executable file
View file

@ -0,0 +1,39 @@
function L = laplacian(nu,dx1,dx2,N1,N2)
%
% Cette fonction construit la matrice de l'opérateur Laplacien 2D anisotrope
%
% Inputs
% ------
%
% nu : nu=[nu1;nu2], coefficients de diffusivité dans les dierctions x1 et x2.
%
% dx1 : pas d'espace dans la direction x1.
%
% dx2 : pas d'espace dans la direction x2.
%
% N1 : nombre de points de grille dans la direction x1.
%
% N2 : nombre de points de grilles dans la direction x2.
%
% Outputs:
% -------
%
% L : Matrice de l'opérateur Laplacien (dimension N1N2 x N1N2)
%
%
g1 = -nu(1)/dx1^2;
g2 = -nu(2)/dx2^2;
a = -2*(g1+g2);
e = ones(N1*N2, 1);
e_inf = e;
e_inf([N1:N1:N1*N2-1]) = 0;
e_sup = e;
e_sup([N1+1:N1:N1*N2-1]) = 0;
L = spdiags([g1*e g2*e_inf a*e g2*e_sup g1*e], [-N2 -1 0 1 N2], N1*N2, N1*N2);
end

89
TP2/advection.m Executable file
View file

@ -0,0 +1,89 @@
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

45
TP2/reference.m Executable file
View file

@ -0,0 +1,45 @@
function uref=reference(ic,lambda,Nx,dx,kt)
% Evaluation de la solution de référence en les points du maillage.
% uref(x)=u0(x-a*t)
%
% Inputs :
% ------
%
% ic : choix de la condition initiale: 0 -> fonction porte, sinon densité gaussienne.
%
% lambda: nombre de Courant.
%
% Nx : Nombre de pas d'espace.
%
% dx : Pas d'espace.
%
% kt : Indice du pas de temps courant.
%
% Output :
% ------
%
% uref : solution évaluée en les noeuds du maillage.
%
%*****************************************
uref=zeros(Nx+2,1);
if(ic==0)
% Fonction porte
xi=0.25;
xj=1;
xx=0:Nx+2;
tmp=find((dx*xx>=xi+dx*lambda*kt)&(dx*xx<=xj+dx*lambda*kt));
uref(tmp)=2;
else
% Densité gaussienne
gaussian=@(x,m,var)exp(-(x-m)^2/(2*var^2))/(var);
n0=1;
var0=0.1;
for j=2:Nx+1
uref(j)=(1/5)*gaussian((j-lambda*kt)*dx,n0,var0);
end
end
end