11.3 M´ ethode de Galerkin
447
if ax>xlim
bp=x/(exp(x)-1); bn=x+bp; return
else
ii=1; fp=1.;fn=1.; df=1.; s=1.;
while abs(df)>eps
ii=ii+1; s=-s; df=df*x/ii;
fp=fp+df; fn=fn+s*df; bp=1./fp; bn=1./fn;
end
return
end
return
Le Programme 86 calcule la matrice de raideur A fe . Les coefficients α, β, γ et le
terme de force f sont remplac´ es par des fonctions constantes par morceaux sur
chaque sous-intervalle du maillage. Les autres int´ egrales de (11.37), impliquant
les fonctions de base et leurs d´ eriv´ ees, sont calcul´ ees de mani` ere exacte.
Programme 86 - femmatr : Construction de la matrice de raideur
function [Afe,rhsbc] = femmatr(Nx,h,alpha,beta,gamma,stabfun)
%FEMMATR Matrice de raideur et second membre.
dd=(alpha(1:Nx-1)+alpha(2:Nx))/h;
dc=-(beta(2:Nx)-beta(1:Nx-1))/2;
dr=h*(gamma(1:Nx-1)+gamma(2:Nx))/3;
ld=-alpha(2:Nx-1)/h; lc=-beta(2:Nx-1)/2; lr=h*gamma(2:Nx-1)/6;
ud=-alpha(2:Nx-1)/h; uc=beta(2:Nx-1)/2; ur=h*gamma(2:Nx-1)/6;
Kd=spdiags([[ld 0]’,dd’,[0 ud]’],-1:1,Nx-1,Nx-1);
Kc=spdiags([[lc 0]’,dc’,[0 uc]’],-1:1,Nx-1,Nx-1);
Kr=spdiags([[lr 0]’,dr’,[0 ur]’],-1:1,Nx-1,Nx-1);
Afe=Kd+Kc+Kr;
if nargin == 6
Ks,rhsbc]=feval(stabfun,Nx,h,alpha,beta); Afe = Afe + Ks;
else
rhsbc = [0, 0];
end
return
La norme H
1 de l’erreur peut ˆ etre calcul´ ee ` a l’aide du Programme 87. On doit
pour cela fournir les macros u et ux contenant la solution exacte u et sa d´ eriv´ ee
u
. La solution num´ erique est stock´ ee dans le vecteur uh, le vecteur coord
contient les coordonn´ ees de la grille et h le pas du maillage. Les int´ egrales
intervenant dans le calcul de la norme H
1 de l’erreur sont ´ evalu´ ees avec la
formule de Simpson composite (8.17).