RAPPORT TP MEF 1D
P1&P2 Matlab – Méthode des Eléments Finis - 2016
Encadré par : M. Anass Rachid
Rochdi Réalisé par : MARTAB Rochdi NAFAE Soukaina Soukaina
Sommaire Introduction : .......................... ........................................ ........................... .......................... ........................... ........................... ........................... ......................... ........... 2 MEF P1 :............................ .......................................... ........................... ........................... ........................... .......................... ........................... ........................... .................. ..... 3 Maillage :.......................... :....................................... ........................... ........................... .......................... .......................... ........................... ........................... ................ ... 3 Page | 1 Les fonctions de forme Phi1 et Phi2 : .......................... ........................................ ........................... ........................... ......................... ........... 3 Calcule des Matrices élémentaires : .......................... ....................................... ........................... ........................... ........................... .............. 4 Calcule des int_phi : .................... ................................. ........................... ........................... .......................... .......................... ........................... ................... ..... 6 Génération des matrices élémentaires : ................................... ................................................ .......................... ......................... ............ 6 Assemblage : .......................... ........................................ ........................... ........................... ........................... .......................... ........................... ....................... ......... 7 Résolution : .......................... ........................................ ........................... .......................... ........................... ........................... ........................... ......................... ........... 8 Main : .......................... ........................................ ........................... ........................... ........................... .......................... ........................... ........................... .................... ....... 9 MEF P2.......................... ........................................ ........................... .......................... ........................... ........................... ........................... ........................... .................... ....... 10 Maillage P2 : .......................... ........................................ ........................... ........................... ........................... .......................... ........................... ..................... ....... 10 Les fonctions de formes : ......................... ....................................... ........................... ........................... ........................... .......................... ................ ... 10 Calcule des matrices élémentaires : ................... ................................. ........................... .......................... ........................... ................... ..... 11 Calcule des int_Phi : ...................... .................................... ........................... .......................... .......................... ........................... ........................... ............... 13 Générations des matrices élémentaires : ........................... ........................................ .......................... ........................... ................. ... 14 Générations da la matrice du second membre : ......................... ....................................... ........................... .................... ....... 14 Assemblage : .......................... ........................................ ........................... ........................... ........................... .......................... ........................... ..................... ....... 15 Résolution : .......................... ........................................ ........................... .......................... ........................... ........................... ........................... ....................... ......... 16 L’interface graphique .......................... ....................................... .......................... .......................... ........................... ........................... ......................... ............ 17
....................................... .......................... .......................... ........................... ........................... ......................... ............ 18 Variables d’entrés : .......................... Conditions des boites de dialogues : ....................... .................................... .......................... ........................... ........................... ............... 19 Calcule : ......................... ....................................... ........................... ........................... ........................... .......................... ........................... ........................... ................ ... 20 Initialisation : ......................... ....................................... ........................... ........................... ........................... .......................... ........................... ..................... ....... 22 Conclusion ........................... ........................................ ........................... ........................... .......................... .......................... ........................... ........................... ............... 22
Introduction : (ME F) est une méthode d’approximation L a méthode des éléments-finis (MEF) d’approximation numérique de solutions de problèmes aux limites statiques ou dynamiques tels que : – – diffusion thermique – – mécanique des milieux continus (solides et fluides) – électromagnétisme mais en en fait, absolument absolument tous les problèmes problèmes d’équations aux dérivées partielles (EDP) aux limites.
Il s’agit, comme dans dans toutes les méthodes
numériques, de trouver une approximation discrète. Pour faire bref, d’un problème différentiel aux limites linéaire, on trouve une formulation variationnel associées équivalente, dont on calcule une approximation de la solution en projetant sur un espace de dimension finie, ce qui revient à résoudre au final un système linéaire. L’appellation éléments finis vient de la décomposition du domaine d’étude en éléments : ils sont souvent représentés par un maillage. Historiquement, l’origine de la méthode peut se trouver dans les travaux de Fermat et Bernoulli ´ (1743) avec le calcul des variations, puis il faut attendre le début du XX éme siècle avec les progrès en analyse avec la méthode de Galerkin se basant sur des théorèmes de projection dans les espaces de Hilbert. En 1943 Robert Courant introduit le principe variationnel avec des fonctions de base a supports locaux ouvrant la voie a une division d’un domaine considéré en ” éléments”. Cependant ce n’est qu’avec le développement des ordinateurs que ces travaux trouvent leurs applications avec les travaux pionniers de Zienckiewiz et Argyris qui définiront la méthode en 1960. Ce qui amène le succès de la méthode et sa puissance est l’apport du calcul matriciel, introduit par un ingénieur civil anonyme. La méthode connait alors un développement fulgurant accompagne par les progrès de l’informatique. La méthode des éléments-finis est une méthode puissante basée sur une théorie mathématique rigoureuse. Aujourd’hui, les éléments-finis sont un outil majeur, incontournable en mécanique (fluides et solides, interactions, structures), et applicable dans de nombreux domaines impliquant des problèmes d’EDP aux limites comme comme par exemple en mathématiques financières ou l’électromagnétisme. l’électromagnétisme. De nombreux codes industriels (solveurs) existent et sont généralement couplés à un logiciel de CAO 1 ou Computer Aide Design (CAD) en Anglais. Citons Ansys, Abaqus, Robot, LS-dyna, Feap, Code- Aster, Aster, Cast3M et bien d’autres.
Page | 2
MEF P1 :
On considère le problème suivant :
−∗ + ∗ +∗=
u(a)=u(b)=0 avec a et b les bords du problème
Maillage : function [X, function [X, T] = MaillageP1(a, b, h) %---------------------------%Génére un maillage de type P1 % X la table des coordonnées % T la table de connectivité moyennant les indices %---------------------------n = floor((b - a)/h) + 1; X = a + h*(0:n-1)'; [~, I] = sort(X); T = [ I(1:n-1), I(2:n)];
Dans La MEF, on commence d’abord avec le avec le maillage, cette fonction nous génère le maillage de P1, elle prend comme entrer le pas et les bords et elle retourne X, la table des coordonnes, et T la table de connectivité.
Les fonctions de forme Phi1 et Phi2 : function
[y y2] = phi1_P1(x, x1, x2)
%calcule la fonction de forme phi1 telle que: % phi1(x1)=1 et phi1(x2)=0 % Ti % |---------------| % x1 x2 % phi1(x)=a*x+b et ab=[a;b] ab=[x1 1;x2 1]\[1;0]; y = polyval(ab,x); y2 = polyval(polyder(ab),x); end function
[y y2] = phi2_P1(x, x1, x2)
ab=[x1 1;x2 1]\[0;1]; y = polyval(ab,x); y2 = polyval(polyder(ab),x); end
Page | 3
Ces deux codes permettent le calcul des deux fonctions de formes phi1 et phi2 en utilisant une commande préprogrammée sur Matlab (polyval) celle-ci retourne les coefficients d'un polynôme entré en argument. Ces codes retournent également les dérives des fonctions de forme moyennant la Page | 4 commande de Matlab Polyder.
Calcule des Matrices élémentaires : function [G] = A_P1(x1,x2) xm=(x1+x2)*0.5; [y1x1 dy1x1] = phi1_P1(x1, x1, x2); [y2x1 dy2x1] = phi2_P1(x1, x1, x2); [y1x2 dy1x2] = phi1_P1(x2, x1, x2); [y2x2 dy2x2] = phi2_P1(x2, x1, x2); [y1xm dy1xm] = phi1_P1(xm, x1, x2); [y2xm dy2xm] = phi2_P1(xm, x1, x2); A = [dy1x1 dy1xm dy1x2;dy2x1 dy2xm dy2x2]; %B = [y1x1 y1xm y1x2;y2x1 y2xm y2x2]; for i=1:2 for i=1:2 for j=1:2 for j=1:2 F=[A(i,1)*A(j,1) A(i,2)*A(j,2) A(i,3)*A(j,3)] y=(x2-x1)/6*(F(1)+4*F(2)+F(3)) G(i,j)=y end end end function [G] = B_P1(x1,x2) function [G] xm=(x1+x2)*0.5; [y1x1 dy1x1] = phi1_P1(x1, x1, x2); [y2x1 dy2x1] = phi2_P1(x1, x1, x2); [y1x2 dy1x2] = phi1_P1(x2, x1, x2); [y2x2 dy2x2] = phi2_P1(x2, x1, x2); [y1xm dy1xm] = phi1_P1(xm, x1, x2); [y2xm dy2xm] = phi2_P1(xm, x1, x2); %A = [dy1x1 dy1xm dy1x2;dy2x1 dy2xm dy2x2]; B = [y1x1 y1xm y1x2;y2x1 y2xm y2x2]; for i=1:2 for i=1:2 for j=1:2 for j=1:2 F=[B(i,1)*B(j,1) B(i,2)*B(j,2) B(i,3)*B(j,3)]; y=(x2-x1)/6*(F(1)+4*F(2)+F(3)); G(i,j)=y; end end end
function [G] = G_P1(x1,x2) function [G] xm=(x1+x2)*0.5; [y1x1 dy1x1] = phi1_P1(x1, x1, x2);
[y2x1 dy2x1] = phi2_P1(x1, x1, x2); [y1x2 dy1x2] = phi1_P1(x2, x1, x2); [y2x2 dy2x2] = phi2_P1(x2, x1, x2); [y1xm dy1xm] = phi1_P1(xm, x1, x2); [y2xm dy2xm] = phi2_P1(xm, x1, x2); A = [dy1x1 dy1xm dy1x2;dy2x1 dy2xm dy2x2]; B = [y1x1 y1xm y1x2;y2x1 y2xm y2x2]; for i=1:2 for i=1:2 for j=1:2 for j=1:2 F=[A(i,1)*B(j,1) A(i,2)*B(j,2) A(i,3)*B(j,3)]; y=(x2-x1)/6*(F(1)+4*F(2)+F(3)); G(i,j)=y; end end end
L e code en dessus nous permet de calculer les trois matrices élémentaires qui sont sous la formes suivantes :
∗ ′ ∫ ∗ ′ ∫ =_1= =_1= ∫ ∗ ′ ∫ ∗ ′ ∫ ∗ ∫ ∗ =_1= =_1= ∫ ∗ ∫ ∗ ∫ ∫ _1= _1= ∫ ∫ Ces codes nous permettent de calculer les trois matrices élémentaires, on fait d’abord appel appel aux fonctions de formes Phi, puis on construit des matrices ou on stock les éléments nécessaires pour le calcul des intégrales par la suite. On utilise la quadrature de Simpson pour calculer les intégrales. Quadrature de Simpson :
∫ =∫ − + =∫ = = +4∗( )+ 6 2 )+
Page | 5
Calcule des int_phi : function y = int_phi2(f,x1,x2) function y %calcule l'intégrale sur l'élément Ti de f*phi1 %moyennant la quadrature de Simpson % Ti Page | 6 % |--------|--------| % x1 xm x2 % xm=(x2+x1)*0.5; % y=(x2x1)/6*(f(x1)*phi2_P1(x1,x1,x2)+4*f(xm)*phi2_P1(xm,x1,x2)+f(x2)*phi2_P1 (x2,x1,x2)); end function y = int_phi1(f,x1,x2) function y %calcule l'intégrale sur l'élément Ti de f*phi1 %moyennant la quadrature de Simpson % Ti % |--------|--------| % x1 xm x2 % xm=(x2+x1)*0.5; % y=(x2x1)/6*(f(x1)*phi1_P1(x1,x1,x2)+4*f(xm)*phi1_P1(xm,x1,x2)+f(x2)*phi1_P1 (x2,x1,x2)); end
Ces deux codes nous permettent de calculer les intégrales sur les éléments Ti de f*Phi.
Génération des matrices élémentaires : function [ elemKi ] = mat_elem_P1(alpha,beta,gama,X,T,i) %calcule la matrice élémentaire dans l'élément Ti % Ti % |---------------| % x1 x2 % x1=X(T(i,1)); x2=X(T(i,2)); elemKi=alpha*A_P1(x1,x2)+beta*B_P1(x1,x2)-gama*G_P1(x1,x2); end
Dans ce code, on fait appel aux matrices déjà construites, et on les assemble en les multipliant par les coefficients correspondants. 1.
Le vecteur second membre
function [ elemFi ] = SM_elem_P1(f,X,T,i) %calcule le second membre élémentaire dans l'élément Ti % Ti % |---------------| % x1 x2 %elle fait appel aux deux fonctions: % y = int_phi1(f,x1,x2) % y = int_phi2(f,x1,x2) x1=X(T(i,1)); x2=X(T(i,2)); elemFi=[int_phi1(f,x1,x2);int_phi2(f,x1,x2)]; end
On construit la matrice du second membre en faisant appel au fonctions int_phi. La matrice du second membre est :
∫ ∗ = ∫ ∗ Assemblage : function [K,F] = AssemblageP1(alpha,beta,gama,f, X, T) function [K,F] % Assemblage des matrices élémentaires "elemki" dans la matrice globale K % Assemblage des seconds membres élémentaires "elemFi" dans le second % membre global F % cette fonction fait appel aux fonctions: % 1) mat_elem_P1(alha,X,T,i) % 2) SM_elem_P1(f,X,T,i) % n = size(X,1); % nombre des noeuds t = size(T,1); % nombre des éléments K = zeros(n); F = zeros(n,1); % Assemblage de la matrice golbale K % for k=1:t for k=1:t % boucle sur les éléments elemKi = mat_elem_P1(alpha,beta,gama,X,T,k); % matrice de l'élément numéro k for i=1:2 for i=1:2 %boucle sur les numéros locaux for j=1:2 for j=1:2 %boucle sur les numéros locaux I=k+i-1; % numéros globaux dans K J=k+j-1; % numéros globaux dans K K(I,J)=K(I,J)+ elemKi(i,j); end end end
Page | 7
% Assemblage du second membre élémentaire golbal F % for k=1:t for k=1:t % boucle sur les éléments elemFi = SM_elem_P1(f,X,T,k); % SM de l'élément numéro k for i=1:2 for i=1:2 %boucle sur les numéros locaux I=k+i-1; % numéros globaux dans K F(I)=F(I)+ elemFi(i); end end end
On initialise d’abord les deux deux matrices K et F, puis on construit le nombre des nœuds et des éléments. On O n réalise l’assemblage à l’aide des boucle for, et cela pour construire la matrice global K et celle du second membre F.
Résolution : function [U] = EF_P1(alpha,beta,gama,f,a,b,h) % fonction traite l'equation alpha*u"=f sur un[a,b] avec u(a)=u(b)=0 % %Creation de la matrice globale K et du second membre globale F [X, T] = MaillageP1(a, b, h); Nn=size(X,1); %nombre des noeuds [K,F] = AssemblageP1(alpha,beta,gama,f, X, T); % % Conditions aux bords for j=1:Nn for j=1:Nn K(1,j)=0.; K(Nn,j)=0.; end K(1,1)=1.; K(Nn,Nn)=1.; % F(1)=0; F(Nn)=0; % % Résolution U=K\F; end
système K.U=F. Cette fonction « \ » nous permet de résoudre le système
Page | 8
Main : clear clc a=-pi/2; b=pi/2 ; h=0.1; alpha=1; beta=0; ue=@(x) cos(x); f=@(x) cos(x); [X, T] = MaillageP1(a, b, h); [U] = EF_P1(alpha,beta,gama, EF_P1(alpha,beta,gama,f,a,b,h); f,a,b,h); % Post-traite Post-traitement ment %%%%%% comparaison graphique figure('name' figure('name', , 'Comparaision: solutions exacte et approchée '); '); fplot(ue, [min(X), max(X)], 'b' 'b'); ); hold on on; ; plot(X,U, 'r.-' 'r.-'); ); %%%%%% Analyse d'erreur erreur=zeros(size(X,1),1); for i=1:size(X,1) i=1:size(X,1) erreur(i)=U(i)-ue(a+(i-1)*h); end figure('name' figure('name', , 'analyse erreur'); erreur'); plot(X,erreur, 'r-o' 'r-o'); ); ylabel('erreur' ylabel('erreur'); ); xlabel('noeuds' xlabel('noeuds'); );
Premièrement on choisit une fonction qui vérifie les conditions au bord, et On construit le problème, puis On fait appel à la fonction de maillage et celle de la résolution. La figure de la solution exacte et approché est tracée sur un graphe ainsi que celle de l’analyse de l’erreur .
Page | 9
MEF P2 On considère toujours le même problème, cette fois on utilise un polynôme de degré deux, ce qui nécessite 3 fonctions de forme chacune de degré deux. Et les Page | 10 matrices seront de forme 3*3.
Maillage P2 : function [X, X1, T] = MaillageP2(a, b, h) %---------------------------%Génére un maillage de type P1 % X la table des coordonnées % T la table de connectivité moyennant les indices %---------------------------n = floor((b - a)/h) + 1; X1 = a + h*(0:n-1)'; X = zeros(2*n-1,1); for j=0:(n-1) X(2*(n-1)-2*j+1) = X1(n-j); if j ~= (n-1) X(2*(n-1)-2*j) = X1(n-j)-0.5*h;
end end [~, I] = sort(X); T = [ I(1:2:(2*n-3)), I(2:2:2*(n-1)), I(3:2:(2*n-1))];
end
s ortie la table de connectivite avec l’ajout d’un Ce maillage nous donne comme sortie point intermédiaire dans chaque élément Ti, et deux tables de cordonnées X1 et X, X contient tous les nœuds tandis que X1 ne contient que les extrémités des Ti.
Les fonctions de formes : function
[y y2] y2] = phi1_P2(x, phi1_P2(x, x1, x2, x3)
%calcule la fonction de forme phi1 telle que: % phi1(x1)=1 et phi1(x2)=0 % Ti % |---------------| % x1 x2 % phi1(x)=a*x+b et ab=[a;b] ab=[x1^2 x1 1;x2^2 x2 1;x3^2 x3 1]\[1;0;0]; y = polyval(ab,x); polyval(ab,x); y2 = polyval(polyder(ab),x);
end function
[y y2] y2] = phi2_P2(x, phi2_P2(x, x1, x2, x3)
%calcule la fonction de forme phi1 telle que: % phi1(x1)=1 et phi1(x2)=0 % Ti % |---------------| % x1 x2 % phi1(x)=a*x+b et ab=[a;b] ab=[x1^2 x1 1;x2^2 x2 1;x3^2 x3 1]\[0;1;0]; y = polyval(ab,x); polyval(ab,x); y2 = polyval(polyder(ab),x);
end function
[y y2] y2] = phi3_P2(x, phi3_P2(x, x1, x2, x3)
%calcule la fonction de forme phi1 telle que: % phi1(x1)=1 et phi1(x2)=0 % Ti % |---------------| % x1 x2 % phi1(x)=a*x+b et ab=[a;b] ab=[x1^2 x1 1;x2^2 x2 1;x3^2 x3 1]\[0;0;1]; y = polyval(ab,x); polyval(ab,x); y2 = polyval(polyder(ab),x);
end
En utilisant les polynômes de degré deux, on se trouve dans l’obligation de travailler avec trois points ce qui veut dire trois fonctions de formes. On obtient donc les trois équations suivantes :
1111 =1;=0; 2222 =0;=1;3333 == 00 11 =0; 22 =0; 33 = 1
On util ise les deux comm andes de Matlab Matlab Polyval et polyder.
Calcule des matrices élémentaires : function [G] = A_P2(x1,x2,x3) xm=(x1+x3)*0.5; [y1x1 dy1x1] = phi1_P2(x1, x1, x2, x3); [y2x1 dy2x1] = phi2_P2(x1, x1, x2, x3); [y3x1 dy3x1] = phi3_P2(x1, x1, x2, x3); [y1x3 dy1x3] = phi1_P2(x3, x1, x2, x3); [y2x3 dy2x3] = phi2_P2(x3, x1, x2, x3); [y3x3 dy3x3] = phi3_P2(x3, x1, x2, x3); [y1xm dy1xm] = phi1_P2(xm, x1, x2, x3); [y2xm dy2xm] = phi2_P2(xm, x1, x2, x3); [y3xm dy3xm] = phi3_P2(xm, x1, x2, x3); A = [dy1x1 dy1xm dy1x3;dy2x1 dy2xm dy2x3;dy3x1 dy2x3;dy3x1 dy3xm dy3x3]; %B = [y1x1 y1xm y1x3;y2x1 y2xm y2x3;y3x1 y3xm y3x3]; for i=1:3 for j=1:3 F=[A(i,1)*A(j,1) F=[A(i,1)*A(j,1) A(i,2)*A(j,2) A(i,2)*A(j,2) A(i,3)*A(j,3)]; A(i,3)*A(j,3)]; y=(x3-x1)/6*(F(1)+4*F(2)+F(3)); G(i,j)=y;
Page | 11
end end end function [G] = B_P2(x1,x2,x3) xm=(x1+x3)*0.5; [y1x1 dy1x1] = phi1_P2(x1, x1, x2, x3); [y2x1 dy2x1] = phi2_P2(x1, x1, x2, x3); [y3x1 dy3x1] = phi3_P2(x1, x1, x2, x3); [y1x3 dy1x3] = phi1_P2(x3, x1, x2, x3); [y2x3 dy2x3] = phi2_P2(x3, x1, x2, x3); [y3x3 dy3x3] = phi3_P2(x3, x1, x2, x3); [y1xm dy1xm] = phi1_P2(xm, x1, x2, x3); [y2xm dy2xm] = phi2_P2(xm, x1, x2, x3); [y3xm dy3xm] = phi3_P2(xm, x1, x2, x3); %A = [dy1x1 dy1xm dy1x3;dy2x1 dy2xm dy2x3;dy3x1 dy3xm dy3x3]; B = [y1x1 y1xm y1x3;y2x1 y2xm y2x3;y3x1 y3xm y3x3]; for i=1:3 for j=1:3 F=[B(i,1)*B(j,1) F=[B(i,1)*B(j,1) B(i,2)*B(j,2) B(i,2)*B(j,2) B(i,3)*B(j,3)]; B(i,3)*B(j,3)]; y=(x3-x1)/6*(F(1)+4*F(2)+F(3)); G(i,j)=y;
end end end
function [G] = G_P2(x1,x2,x3) xm=(x1+x3)*0.5; [y1x1 dy1x1] = phi1_P2(x1, x1, x2, x3); [y2x1 dy2x1] = phi2_P2(x1, x1, x2, x3); [y3x1 dy3x1] = phi3_P2(x1, x1, x2, x3); [y1x3 dy1x3] = phi1_P2(x3, x1, x2, x3); [y2x3 dy2x3] = phi2_P2(x3, x1, x2, x3); [y3x3 dy3x3] = phi3_P2(x3, x1, x2, x3); [y1xm dy1xm] = phi1_P2(xm, x1, x2, x3); [y2xm dy2xm] = phi2_P2(xm, x1, x2, x3); [y3xm dy3xm] = phi3_P2(xm, x1, x2, x3); A = [dy1x1 dy1xm dy1x3;dy2x1 dy2xm dy2x3;dy3x1 dy2x3;dy3x1 dy3xm dy3x3]; B = [y1x1 y1xm y1x3;y2x1 y2xm y2x3;y3x1 y3xm y3x3]; for i=1:3 for j=1:3 F=[A(i,1)*B(j,1) F=[A(i,1)*B(j,1) A(i,2)*B(j,2) A(i,2)*B(j,2) A(i,3)*B(j,3)]; A(i,3)*B(j,3)]; y=(x3-x1)/6*(F(1)+4*F(2)+F(3)); G(i,j)=y;
end end end
Les nouveaux matri ces sont d e types 3*3 3*3 :
Page | 12
∫ ′ ′ ∫ ′ ′ ∫ ′ ′ =_2= ∫ ′ ′ ∫ ′ ′ ∫ ′ ′ [∫∫′ ′ ∫∫′′ ∫∫ ′′] =_2= ∫ ∫ ∫ ∫ ∫ ∫ ∫ ′ ∫ ′ ∫ ′ _2= ∫ ′ ∫ ′ ∫ ′ ∫ ′ ∫ ′ ∫ ′ Calcule des int_Phi : function y function y = int_phi1_P2(f,x1,x2,x3) %calcule l'intégrale sur l'élément Ti de f*phi1 %moyennant la quadrature de Simpson % Ti % |--------|--------| % x1 xm x2 %
xm=(x3+x1)*0.5; %
y=(x3x1)/6*(f(x1)*phi1_P2(x1,x1,x2,x3)+4*f(xm)*phi1_P2(xm,x1,x2,x3)+f(x3)*p hi1_P2(x3,x1,x2,x3)) end
function y function y = int_phi3_P2(f,x1,x2,x3) %calcule l'intégrale sur l'élément Ti de f*phi1 %moyennant la quadrature de Simpson % Ti % |--------|--------| % x1 xm x2 %
xm=(x3+x1)*0.5; %
y=(x3x1)/6*(f(x1)*phi3_P2(x1,x1,x2,x3)+4*f(xm)*phi3_P2(xm,x1,x2,x3)+f(x3)*p hi3_P2(x3,x1,x2,x3)); end
function y function y = int_phi2_P2(f,x1,x2,x3) %calcule l'intégrale sur l'élément Ti de f*phi1
Page | 13
%moyennant la quadrature de Simpson % Ti % |--------|--------| % x1 xm x2 %
xm=(x3+x1)*0.5; %
Page | 14 y=(x3x1)/6*(f(x1)*phi2_P2(x1,x1,x2,x3)+4*f(xm)*phi2_P2(xm,x1,x2,x3)+f(x3)*p hi2_P2(x3,x1,x2,x3)); end
On calcule les int_Phi. On procède de la même façon qu’auparavant
Générations des matrices élémentaires : function [ function [ elemKi ] = mat_elem_P2(alpha,beta,gama,X,T,i) %calcule la matrice élémentaire dans l'élément Ti % Ti % |---------------| % x1 x2 %
x1=X(T(i,1)); x3=X(T(i,3)); x2=(x1+x3)/2; elemKi=alpha*A_P2(x1,x2,x3)+beta*B_P2(x1,x2,x3)-gama*G_P2(x1,x2,x3); end
De la même manière qu ’en P1.
Générations da la matrice du second membre : function [ function [ elemFi ] = SM_elem_P2(f,X,T,i) %calcule le second membre élémentaire dans l'élément Ti % Ti % |---------------| % x1 x2 %elle fait appel aux deux fonctions: % y = int_phi1(f,x1, int_phi1(f,x1,x2) x2) % y = int_phi2(f,x1, int_phi2(f,x1,x2) x2)
x1=X(T(i,1)); x2=X(T(i,2)); x3=X(T(i,3)); elemFi=[int_phi1_P2(f,x1,x2,x3);int_phi2_P2(f,x1,x2,x3);int_phi3_P2(f, x1,x2,x3)]; %
end
L a matrice du second membre de P2 est :
∫ ∗ = ∫ ∗ [∫ ∗ ] Assemblage : function [K,F] function [K,F] = AssemblageP2(alpha,beta,gama,f, X, T) % % % % % % %
Assemblage des matrices élémentaires "elemki" dans la matrice globale K Assemblage des seconds membres élémentaires "elemFi" dans le second membre global F cette fonction fait appel aux fonctions: 1) mat_elem_P1 mat_elem_P1(alha,X,T,i (alha,X,T,i) ) 2) SM_elem_P1( SM_elem_P1(f,X,T,i) f,X,T,i)
n = size(X,1); t = size(T,1);
% nombre des noeuds % nombre des éléments
K = zeros(n); F = zeros(n,1); % Assemblage de la matrice golbale K
c=0; for k=1:2:(2*t-1) for k=1:2:(2*t-1) % boucle sur les éléments c=c+1; elemKi = mat_elem_P2(alpha,beta,gama,X,T,c); % matrice de l'élément numéro k for i=1:3 for i=1:3 %boucle sur les numéros locaux for j=1:3 for j=1:3 %boucle sur les numéros locaux % numéros globaux dans K I=k+i-1; J=k+j-1; % numéros globaux dans K % l'assemblage K(I,J)=K(I,J)+ elemKi(i,j); end end end % Assemblage du second membre élémentaire golbal F
c=0; for k=1:2:(2*t-1) for k=1:2:(2*t-1) c=c+1; % boucle sur les éléments elemFi = SM_elem_P2(f,X,T,c); % SM de l'élément numéro k for i=1:3 for i=1:3 %boucle sur les numéros locaux % numéros globaux dans K I=k+i-1; F(I)=F(I)+ elemFi(i); % l'assemblage end end end
Page | 15
Cette fonction permet de faire l’assemblage des matrices élémentaires elemFi ainsi que l’assemblage des elemAi de la même manière qu’en qu’en P1, le changement s’impose au niveau du nombre de nœuds n nœuds n et les incrémentations dans chaque boucle. Page | 16
Résolution : function [U] function [U] = EF_P2(alpha,beta,gama,f,a,b,h) % fonction traite l'equation alpha*u"=f sur un[a,b] avec u(a)=u(b)=0 % %Creation de la matrice globale K et du second membre globale F
[X, X1, T] = MaillageP2(a, b, h); Nn=size(X,1); %nombre des noeuds [K,F] = AssemblageP2(alpha,beta,gama,f, X, T); % % Conditions aux bords
for j=1:Nn for j=1:Nn K(1,j)=0.; K(Nn,j)=0.; end K(1,1)=1.; K(Nn,Nn)=1.; % F(1)=0; F(Nn)=0;
% % Résolution U=K\F;
end
De la même manière on procède pour résoudre l’équation et trouver U.
L interface graphique ’
Pour créer une interface sur Matlab on tape Guide sur la commande window, une interface .fig s’affiche ou on peut ajouter des boutons, des graphs, des choix etc. Pour Page | 17 chaque élément ajouté on doit modifier son tag pour faciliter la manipulation du code et éviter l’ambiguïté l’utiliser après. Une fois on enregistre le .fig un code est est générer.
Un fichier PDF s’ouvre pour une aide simplifiée à propos de l’interface : l’interface :
Page | 18
Variables d’entrés :
travers la fonction suivante suivante pour recevoir recevoir les On fait appel aux tag ajoutées à travers éléments entrées par l’utilisateur de l’interface : l’interface :
tag = str2doublegethandles.tag,′String′;
alpha = str2double(get(handles str2double(get(handles.alpha, .alpha,'String' 'String')); )); h = str2double(get(handles str2double(get(handles.h, .h,'String' 'String')); )); a = str2num(get(handles.a, str2num(get(handles.a,'String' 'String')); )); b = str2num(get(handles.b, str2num(get(handles.b,'String' 'String'));…………
Conditions des boites de dialogues : On a essayé de contrôler les données entrées par l’utilisateur par des conditions pour ne pas avoir des données erronées , par exemple sur l’intervalle l’intervalle [a, b], a doit être plus petit que b Page | 19 (a
L e script est pour appeler la boite de dialogue suivante après qu’on appuis sur le bouton Valider :
m aillage avant de Valider les données On ne peut effectuer aucun calcule ou maillage d’entrés par le bouton Valider (ce qui laisse le programme de vérifier tous les données données d’entrer s’ils sont valident sont valident avant de faire aucun calcule) Pour les boites dialogues (modaldlg1 à 8), On a créé d’autres interfaces, quand fait appel à eux selon certaines conditions, pour guider l’utilisateur.
Calcule : L ’interface nous permet de faire seulement un calcule d’approximation sans avoir une solution exacte pour calculer l’erreur si on veut ( On doit seulement décocher la Page | 20 petite case à coté de « Solution U »
Pour qu’un bouton effectue la fonction désirée, On fait appel aux codes déjà réalisés sur Matlab. Pour les boutons suivants :
L e programme garde en mémoire les valeurs d’approximation P1 et P2 dont la Matrice S : kj=floor(2*((b-a)/h)+1); S(1:(kj+1)/2,1)=U1;
… [U2] = EF_P2(alpha,beta,gama,f,a,b,h); S(:,2)=U2;
compara ison d’erreur entre l’approximation P1 & l’approximation Pour avoir une comparaison P2.
Exemple :
Page | 21
En appuient sur « Comp Graph P1&P2 » une comparaison graphique des deux graphes s’affiche, mais ce n’est pas évident de voir les erreurs entre les 2 approximation. Après Après qu’on appuie sur « Erreur P1&P2 », il est plus claire la différence d’erreurs entre les 2 approximations (avec un tableau affichant les valeurs des erreurs P1 sur la 1ere ligne et les valeurs des erreurs de P2 sur la 2eme ligne :
Initialisation : les donnée et remettre le tous comme on a ouvri l’interface pour la Pour initialiser les 1ere fois, un bouton de Reset est programmé à faire ça : Page | 22 clear; %Vider la memoir des variables clc; %Vider la fenêtre des commande si elle est activée cla; %Vider les axes d’affichage des graphes set(findobj(0,'style' set(findobj(0, 'style', ,'edit' 'edit'), ),'string' 'string', ,' '); '); %Vider les valeur d’entr set(findobj(0,'style' set(findobj(0, 'style', ,'text' 'text', ,'tag' 'tag', ,'equation' 'equation'), ),'string' 'string', ,' '); '); %vider le panel d’affichage de equation entrée
set(findobj(0,'tag' set(findobj(0, 'tag', ,'val' 'val'), ),'data' 'data',0); ,0); % Vider le tableau des valeurs set(findobj(0,'tag' set(findobj(0, 'tag', ,'erg' 'erg'), ),'data' 'data',0); ,0); % Vider le tableau des erreurs
Le reste du
code de l’interface e st
dans les les documents joints sou s format
(.m)
Conclusion
Ce premier projet sur Matlab nous a énormément aidés à comprendre la notion et l’ordonnancement l’ordonnancement des étapes de résolution en éléments finis en 1D pour P1et P2. On a pu également grâce à ce projet de travailler avec l’interfaces graphiques fournie par Matlab, et savoir comment le convertir en. Exe.
