302
8 Approximation numérique des problèmes aux limites
Programme 8.4. newmarkwave : méthode de Newmark pour l’équation des
ondes
function [xh , uh ]= n e w m arkwa ve( xspan , tspan , nstep , param ,...
c ,u0 , v0 ,g ,f , varargin )
% N E W M ARK WAVE résout l ’ équation des ondes avec la
% méthode de Newmark .
% [ XH , UH ]= N E W M ARKW AVE( XSPAN , TSPAN , NSTEP , PARAM ,C ,...
% U0 , V0 ,G ,F ) résout l ’ équation des ondes
% D ^2 U / DT ^2 - C D ^2 U/ DX ^2 = F
% dans ] XSPAN (1) ,XSPAN (2)[ x ] TSPAN (1) , TSPAN (2)[ en
% u t i lisa nt la méthode de Newmark avec les c o n d it ions
% i n i tial es U (X ,0)= U0( X) , DU / DX(X ,0)= V0 (X ) et les
% c o n di tions de D i r i chlet U(X , T )= G(X , T) pour X= XSPAN (1)
% et X= XSPAN (2). C est une c o n st ante positive .
% NSTEP (1) est le nombre de pas d ’ i n t é grati on en
% espace , NSTEP (2) est le nombre de pas d ’ i n t é grat ion
% en temps . PARAM (1)= ZETA et PARAM (2)= THETA .
% U0( X) , V0 (X ) , G(X , T) et F(x , T) sont des f o n ct ions
% inline , anonymes ou définies par un M - file .
% XH contient les noeuds de d i s c réti sat ion.
% UH contient la solution n u m é riqu e au temps TSPAN (2).}
% [ XH , UH ]= N E W M ARKW AVE( XSPAN , TSPAN , NSTEP , PARAM ,C ,...
% U0 , V0 ,G ,F ,P1 , P2 ,...) passe les p a r a mètr es
% s u p p lé ment air es P1 , P2 ,... aux f o n ctio ns U0 , V0 ,G , F.
h = ( xspan (2) -xspan (1))/ nstep (1);
dt = ( tspan (2) -tspan (1))/ nstep (2);
zeta = param (1); theta = param (2);
N = nstep (1)+1;
e = ones (N ,1); D = spdiags ([ e -2* e e ] ,[ -1 ,0 ,1] ,N ,N );
I = speye (N ); lambda = dt/ h;
A = I -c * lambda ^2* zeta * D;
An = I+ c* lambda ^2*(0.5 - zeta )*D ;
A (1 ,:) = 0; A (1 ,1) = 1; A(N ,:) = 0; A (N ,N ) = 1;
xh = ( linspace ( xspan (1) , xspan (2) ,N )) ’;
fn = feval (f , xh , tspan (1) , varargin {:});
un = feval ( u0 ,xh , varargin {:});
vn = feval ( v0 ,xh , varargin {:});
[L , U ]= lu( A );
alpha = dt ^2* zeta ; beta = dt ^2*(0.5 - zeta );
theta1 = 1 - theta ;
for t = tspan (1)+ dt: dt : tspan (2)
fn1 = feval (f , xh ,t , varargin {:});
rhs = An* un+ dt *I * vn+ alpha * fn1 + beta * fn ;
temp = feval (g ,[ xspan (1) , xspan (2)] ,t , varargin {:});
rhs ([1 ,N ]) = temp ;
uh = L\ rhs;
uh = U\ uh ;
v = vn + dt *((1 -theta )*( c* D* un/ h ^2+ fn )+...
theta *(c *D * uh/ h ^2+ fn1 ));
fn = fn1;
un = uh ;
vn = v ;
end
Comme alternative au schéma de Newmark, on peut considérer le
schéma saute-mouton
u
n+1
j
− 2u
n
j + u
n−1
j
= c
Δt
Δx
2
(u
n
j+1 − 2u
n
j + u
n
j−1 ),
(8.79)
8 Approximation numérique des problèmes aux limites
Programme 8.4. newmarkwave : méthode de Newmark pour l’équation des
ondes
function [xh , uh ]= n e w m arkwa ve( xspan , tspan , nstep , param ,...
c ,u0 , v0 ,g ,f , varargin )
% N E W M ARK WAVE résout l ’ équation des ondes avec la
% méthode de Newmark .
% [ XH , UH ]= N E W M ARKW AVE( XSPAN , TSPAN , NSTEP , PARAM ,C ,...
% U0 , V0 ,G ,F ) résout l ’ équation des ondes
% D ^2 U / DT ^2 - C D ^2 U/ DX ^2 = F
% dans ] XSPAN (1) ,XSPAN (2)[ x ] TSPAN (1) , TSPAN (2)[ en
% u t i lisa nt la méthode de Newmark avec les c o n d it ions
% i n i tial es U (X ,0)= U0( X) , DU / DX(X ,0)= V0 (X ) et les
% c o n di tions de D i r i chlet U(X , T )= G(X , T) pour X= XSPAN (1)
% et X= XSPAN (2). C est une c o n st ante positive .
% NSTEP (1) est le nombre de pas d ’ i n t é grati on en
% espace , NSTEP (2) est le nombre de pas d ’ i n t é grat ion
% en temps . PARAM (1)= ZETA et PARAM (2)= THETA .
% U0( X) , V0 (X ) , G(X , T) et F(x , T) sont des f o n ct ions
% inline , anonymes ou définies par un M - file .
% XH contient les noeuds de d i s c réti sat ion.
% UH contient la solution n u m é riqu e au temps TSPAN (2).}
% [ XH , UH ]= N E W M ARKW AVE( XSPAN , TSPAN , NSTEP , PARAM ,C ,...
% U0 , V0 ,G ,F ,P1 , P2 ,...) passe les p a r a mètr es
% s u p p lé ment air es P1 , P2 ,... aux f o n ctio ns U0 , V0 ,G , F.
h = ( xspan (2) -xspan (1))/ nstep (1);
dt = ( tspan (2) -tspan (1))/ nstep (2);
zeta = param (1); theta = param (2);
N = nstep (1)+1;
e = ones (N ,1); D = spdiags ([ e -2* e e ] ,[ -1 ,0 ,1] ,N ,N );
I = speye (N ); lambda = dt/ h;
A = I -c * lambda ^2* zeta * D;
An = I+ c* lambda ^2*(0.5 - zeta )*D ;
A (1 ,:) = 0; A (1 ,1) = 1; A(N ,:) = 0; A (N ,N ) = 1;
xh = ( linspace ( xspan (1) , xspan (2) ,N )) ’;
fn = feval (f , xh , tspan (1) , varargin {:});
un = feval ( u0 ,xh , varargin {:});
vn = feval ( v0 ,xh , varargin {:});
[L , U ]= lu( A );
alpha = dt ^2* zeta ; beta = dt ^2*(0.5 - zeta );
theta1 = 1 - theta ;
for t = tspan (1)+ dt: dt : tspan (2)
fn1 = feval (f , xh ,t , varargin {:});
rhs = An* un+ dt *I * vn+ alpha * fn1 + beta * fn ;
temp = feval (g ,[ xspan (1) , xspan (2)] ,t , varargin {:});
rhs ([1 ,N ]) = temp ;
uh = L\ rhs;
uh = U\ uh ;
v = vn + dt *((1 -theta )*( c* D* un/ h ^2+ fn )+...
theta *(c *D * uh/ h ^2+ fn1 ));
fn = fn1;
un = uh ;
vn = v ;
end
Comme alternative au schéma de Newmark, on peut considérer le
schéma saute-mouton
u
n+1
j
− 2u
n
j + u
n−1
j
= c
Δt
Δx
2
(u
n
j+1 − 2u
n
j + u
n
j−1 ),
(8.79)