466
12 Probl` emes transitoires paraboliques et hyperboliques
Programme 88 - thetameth : θ-sch´ ema pour l’´ equation de la chaleur
function [u,x] = thetameth(I,n,u0,f,bc,nu,theta)
%THETAMETH Theta-sch´ ema.
% [U,X]=THETAMETH(I,N,U0,F,BC,NU,THETA) r´ esout l’´ equation de la chaleur
% avec le THETA-sch´ ema.
nx=n(1); h=(I(2)-I(1))/nx; x=[I(1):h:I(2)]; t=I(3);
uold=(eval(u0))’; nt=n(2); k=(I(4)-I(3))/nt; e=ones(nx+1,1);
K=spdiags([(h/(6*k)-nu*theta/h)*e, (2*h/(3*k)+2*nu*theta/h)*e, ...
(h/(6*k)-nu*theta/h)*e],-1:1,nx+1,nx+1);
B=spdiags([(h/(6*k)+nu*(1-theta)/h)*e, (2*h/(3*k)-nu*2*(1-theta)/h)*e, ...
(h/(6*k)+nu*(1-theta)/h)*e],-1:1,nx+1,nx+1);
K(1,1)=1;
K(1,2)=0;
B(1,1)= 0;
B(1,2)=0;
K(nx+1,nx+1)=1; K(nx+1,nx)=0; B(nx+1,nx+1)=0; B(nx+1,nx)=0;
[L,U]=lu(K);
t=I(3);
x=[I(1)+h:h:I(2)-h];
fold=(eval(f))’;
fold=h*fold;
fold=[bc(1); fold; bc(2)];
for time=I(3)+k:k:I(4)
t=time;
fnew=(eval(f))’; fnew=h*fnew; fnew=[bc(1); fnew; bc(2)];
b=theta*fnew+(1-theta)*fold+B*uold;
y=L\b; u=U\y; uold=u;
end
x=[I(1):h:I(2)];
return
Exemple 12.1 V´ erifions la pr´ ecision en temps du θ-sch´ ema appliqu´ e ` a l’´ equation de
la chaleur (12.1) sur le domaine espace-temps ]0, 1[×]0, 1[, avec un f choisi de sorte
que la solution exacte soit u = sin(2πx) cos(2πt). On fixe le pas d’espace h = 1/500
et on prend un pas de temps ∆t ´ egal `
a (10k)
−1 , avec k = 1, . . . , 4. On utilise des
´ el´ ement finis affines par morceaux pour la discr´ etisation en espace. La Figure 12.1
montre la convergence du sch´ ema d’Euler r´ etrograde (θ = 1, trait plein) du sch´ ema
de Crank-Nicolson (θ = 1/2, trait discontinu) en norme L
2 (]0, 1[) (´ evalu´ ee au temps
t = 1), quand ∆t tend vers z´ ero. Comme pr´ evu, la m´ ethode de Crank-Nicolson est
beaucoup plus pr´ ecise que celle d’Euler.
•