From b7661507e85c3a4c1449374537f345b92097c9ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Laure=CE=B7t?= Date: Sat, 10 Jun 2023 20:51:59 +0200 Subject: [PATCH] init --- TP1/Affichage/SetClear.m | 13 ++++++ TP1/Affichage/SetCtrl.m | 64 +++++++++++++++++++++++++++++ TP1/Affichage/SetWait.m | 28 +++++++++++++ TP1/Affichage/plot_uh.m | 59 ++++++++++++++++++++++++++ TP1/diffusivity.m | 45 ++++++++++++++++++++ TP1/forcing.m | 27 ++++++++++++ TP1/laplacian.m | 39 ++++++++++++++++++ TP2/advection.m | 89 ++++++++++++++++++++++++++++++++++++++++ TP2/reference.m | 45 ++++++++++++++++++++ 9 files changed, 409 insertions(+) create mode 100755 TP1/Affichage/SetClear.m create mode 100755 TP1/Affichage/SetCtrl.m create mode 100755 TP1/Affichage/SetWait.m create mode 100755 TP1/Affichage/plot_uh.m create mode 100755 TP1/diffusivity.m create mode 100755 TP1/forcing.m create mode 100755 TP1/laplacian.m create mode 100755 TP2/advection.m create mode 100755 TP2/reference.m diff --git a/TP1/Affichage/SetClear.m b/TP1/Affichage/SetClear.m new file mode 100755 index 0000000..0a4cee9 --- /dev/null +++ b/TP1/Affichage/SetClear.m @@ -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 \ No newline at end of file diff --git a/TP1/Affichage/SetCtrl.m b/TP1/Affichage/SetCtrl.m new file mode 100755 index 0000000..83201db --- /dev/null +++ b/TP1/Affichage/SetCtrl.m @@ -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 ); diff --git a/TP1/Affichage/SetWait.m b/TP1/Affichage/SetWait.m new file mode 100755 index 0000000..3a08f33 --- /dev/null +++ b/TP1/Affichage/SetWait.m @@ -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 diff --git a/TP1/Affichage/plot_uh.m b/TP1/Affichage/plot_uh.m new file mode 100755 index 0000000..2c771c4 --- /dev/null +++ b/TP1/Affichage/plot_uh.m @@ -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 diff --git a/TP1/diffusivity.m b/TP1/diffusivity.m new file mode 100755 index 0000000..5370ccb --- /dev/null +++ b/TP1/diffusivity.m @@ -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 diff --git a/TP1/forcing.m b/TP1/forcing.m new file mode 100755 index 0000000..d3ddd98 --- /dev/null +++ b/TP1/forcing.m @@ -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 diff --git a/TP1/laplacian.m b/TP1/laplacian.m new file mode 100755 index 0000000..66a2ce7 --- /dev/null +++ b/TP1/laplacian.m @@ -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 diff --git a/TP2/advection.m b/TP2/advection.m new file mode 100755 index 0000000..998c511 --- /dev/null +++ b/TP2/advection.m @@ -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 diff --git a/TP2/reference.m b/TP2/reference.m new file mode 100755 index 0000000..f401b45 --- /dev/null +++ b/TP2/reference.m @@ -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