11.3 M´ ethode de Galerkin
445
11.3.7 Probl` emes d’impl´ ementation
Dans cette section, on propose une impl´ ementation de la m´ ethode des ´ el´ ements
finis affines par morceaux (k = 1) pour le probl` eme aux limites (11.37) avec
des conditions aux limites de Dirichlet non homog` enes.
Voici la liste des param` etres d’entr´ ee du Programme 82 : Nx est le nombre
de sous-intervalles ; I est l’intervalle [a, b], alpha, beta, gamma et f les macros correspondant aux coefficients de l’´ equation, bc=[ua,ub] est un vecteur
contenant les conditions aux limites de Dirichlet pour u en x = a et x = b.
Le param` etre stabfun est une chaˆ ıne de caract` eres optionnelle permettant `
a
l’utilisateur de choisir la m´ ethode de viscosit´ e artificielle pour les probl` emes
du type de ceux consid´ er´ es ` a la Section 11.3.6.
Programme 82 - ellfem : El´ ements finis affines 1D
function [uh,x] = ellfem(Nx,I,alpha,beta,gamma,f,bc,stabfun)
%ELLFEM Solveur ´ el´ ements finis.
% [UH,X]=ELLFEM(NX,I,ALPHA,BETA,GAMMA,F,BC) r´ esout le probl` eme aux
% limites:
%
-ALPHA*U’’+BETA*U’+GAMMA=F in (I(1),I(2))
%
U(I(1))=BC(1),
U(I(2))=BC(2)
% avec des ´ el´ ements finis affines. ALPHA, BETA, GAMMA et F peuvent
% ˆ etre des fonctions inline.
% [UH,X]=ELLFEM(NX,I,ALPHA,BETA,GAMMA,F,BC,STABFUN) r´ esout le
% probl` eme avec des ´ el´ ements finis stabilis´ es:
% STABFUN=’upwind’ pour la m´ ethode d´ ecentr´ ee;
% STABFUN=’sgvisc’ pour le fitting exponentiel.
%
a=I(1); b=I(2); h=(b-a)/Nx; x=[a+h/2:h:b-h/2];
valpha=feval(alpha,x); vbeta=feval(beta,x);
vgamma=feval(gamma,x); vf=feval(f,x);
rhs=0.5*h*(vf(1:Nx-1)+vf(2:Nx));
if nargin == 8
[Afe,rhsbc]=femmatr(Nx,h,valpha,vbeta,vgamma,stabfun);
else
[Afe,rhsbc]=femmatr(Nx,h,valpha,vbeta,vgamma);
end
[L,U,P]=lu(Afe);
rhs(1)=rhs(1)-bc(1)*(-valpha(1)/h-vbeta(1)/2+h*vgamma(1)/3+rhsbc(1));
rhs(Nx-1)=rhs(Nx-1)-bc(2)*(-valpha(Nx)/h+vbeta(Nx)/2+...
h*vgamma(Nx)/3+rhsbc(2));
rhs=P*rhs’;
z=L\rhs; w=U\z;
uh=[bc(1), w’, bc(2)]; x=[a:h:b];
return