9.10 Approximation des d´ eriv´ ees
363
favorable quand f est une fonction p´ eriodique de p´ eriode b − a, auquel cas
u i+n = u i pour tout i ∈ Z. Dans le cas non p´ eriodique, le syst` eme (9.64)
doit ˆ etre compl´ et´ e par des relations aux noeuds voisins des extr´ emit´ es de
l’intervalle d’approximation. Par exemple, la d´ eriv´ ee premi` ere en x 0 peut ˆ etre
calcul´ ee en utilisant la relation
u 0 + αu 1 =
1
h
(Af 1 + Bf 2 + Cf 3 + Df 4 ),
et en imposant
A = −
3 + α + 2D
2
, B = 2 + 3D, C = −
1 − α + 6D
2
,
afin que le sch´ ema soit au moins pr´ ecis ` a l’ordre deux (voir [Lel92] pour les
relations `
a imposer dans le cas des m´ ethodes d’ordre plus ´ elev´ e).
Le Programme 79 propose une impl´ ementation MATLAB des sch´ emas aux
diff´ erences finies compactes (9.64) pour l’approximation de la d´ eriv´ ee d’une
fonction f suppos´ ee p´ eriodique sur l’intervalle [a, b[. Les param` etres d’entr´ ee
alpha, beta et gamma contiennent les coefficients du sch´ ema, a et b sont les
extr´ emit´ es de l’intervalle, f est une chaˆ ıne contenant l’expression de f et n
d´ esigne le nombre de sous-intervalles de [a, b]. En sortie les vecteurs u et x
contiennent les valeurs approch´ ees u i et les coordonn´ ees des noeuds. Remarquer qu’en posant alpha=gamma=0 et beta=1, on retrouve l’approximation
par diff´ erences finies centr´ ees (9.62).
Programme 79 - compdiff : Sch´ emas aux diff´ erences finies compactes
function [u,x] = compdiff(alpha,beta,gamma,a,b,n,f)
%COMPDIFF Sch´ ema aux diff´ erences finies compactes.
% [U,X]=COMPDIFF(ALPHA,BETA,GAMMA,A,B,N,F) calcule la d´ eriv´ ee
% premi` ere d’une fonction F sur l’intervalle ]A,B[ en utilisant un
% sch´ ema aux diff´ erences finies compactes avec les
% coefficients ALPHA, BETA et GAMMA.
h=(b-a)/(n+1); x=[a:h:b]; fx = eval(f);
A=eye(n+2)+alpha*diag(ones(n+1,1),1)+alpha*diag(ones(n+1,1),-1);
rhs=0.5*beta/h*(fx(4:n+1)-fx(2:n-1))+0.25*gamma/h*(fx(5:n+2)-fx(1:n-2));
if gamma == 0
rhs=[0.5*beta/h*(fx(3)-fx(1)), rhs, 0.5*beta/h*(fx(n+2)-fx(n))];
A(1,1:n+2)=zeros(1,n+2);
A(1,1)= 1; A(1,2)=alpha; A(1,n+1)=alpha;
rhs=[0.5*beta/h*(fx(2)-fx(n+1)), rhs];
A(n+2,1:n+2)=zeros(1,n+2);
A(n+2,n+2)=1; A(n+2,n+1)=alpha; A(n+2,2)=alpha;
rhs=[rhs, 0.5*beta/h*(fx(2)-fx(n+1))];
else
rhs=[0.5*beta/h*(fx(3)-fx(1))+0.25*gamma/h*(fx(4)-fx(n+1)), rhs];
363
favorable quand f est une fonction p´ eriodique de p´ eriode b − a, auquel cas
u i+n = u i pour tout i ∈ Z. Dans le cas non p´ eriodique, le syst` eme (9.64)
doit ˆ etre compl´ et´ e par des relations aux noeuds voisins des extr´ emit´ es de
l’intervalle d’approximation. Par exemple, la d´ eriv´ ee premi` ere en x 0 peut ˆ etre
calcul´ ee en utilisant la relation
u 0 + αu 1 =
1
h
(Af 1 + Bf 2 + Cf 3 + Df 4 ),
et en imposant
A = −
3 + α + 2D
2
, B = 2 + 3D, C = −
1 − α + 6D
2
,
afin que le sch´ ema soit au moins pr´ ecis ` a l’ordre deux (voir [Lel92] pour les
relations `
a imposer dans le cas des m´ ethodes d’ordre plus ´ elev´ e).
Le Programme 79 propose une impl´ ementation MATLAB des sch´ emas aux
diff´ erences finies compactes (9.64) pour l’approximation de la d´ eriv´ ee d’une
fonction f suppos´ ee p´ eriodique sur l’intervalle [a, b[. Les param` etres d’entr´ ee
alpha, beta et gamma contiennent les coefficients du sch´ ema, a et b sont les
extr´ emit´ es de l’intervalle, f est une chaˆ ıne contenant l’expression de f et n
d´ esigne le nombre de sous-intervalles de [a, b]. En sortie les vecteurs u et x
contiennent les valeurs approch´ ees u i et les coordonn´ ees des noeuds. Remarquer qu’en posant alpha=gamma=0 et beta=1, on retrouve l’approximation
par diff´ erences finies centr´ ees (9.62).
Programme 79 - compdiff : Sch´ emas aux diff´ erences finies compactes
function [u,x] = compdiff(alpha,beta,gamma,a,b,n,f)
%COMPDIFF Sch´ ema aux diff´ erences finies compactes.
% [U,X]=COMPDIFF(ALPHA,BETA,GAMMA,A,B,N,F) calcule la d´ eriv´ ee
% premi` ere d’une fonction F sur l’intervalle ]A,B[ en utilisant un
% sch´ ema aux diff´ erences finies compactes avec les
% coefficients ALPHA, BETA et GAMMA.
h=(b-a)/(n+1); x=[a:h:b]; fx = eval(f);
A=eye(n+2)+alpha*diag(ones(n+1,1),1)+alpha*diag(ones(n+1,1),-1);
rhs=0.5*beta/h*(fx(4:n+1)-fx(2:n-1))+0.25*gamma/h*(fx(5:n+2)-fx(1:n-2));
if gamma == 0
rhs=[0.5*beta/h*(fx(3)-fx(1)), rhs, 0.5*beta/h*(fx(n+2)-fx(n))];
A(1,1:n+2)=zeros(1,n+2);
A(1,1)= 1; A(1,2)=alpha; A(1,n+1)=alpha;
rhs=[0.5*beta/h*(fx(2)-fx(n+1)), rhs];
A(n+2,1:n+2)=zeros(1,n+2);
A(n+2,n+2)=1; A(n+2,n+1)=alpha; A(n+2,2)=alpha;
rhs=[rhs, 0.5*beta/h*(fx(2)-fx(n+1))];
else
rhs=[0.5*beta/h*(fx(3)-fx(1))+0.25*gamma/h*(fx(4)-fx(n+1)), rhs];