246
7 Equations différentielles ordinaires
Programme 7.8. newmark : méthode de Newmark
function [t ,u ]= newmark ( odefun , tspan ,y0 , Nh , param ,...
varargin )
% NEWMARK résout une équation d i f f ér ent iell e du second
% ordre avec la méthode de Newmark
% [T ,Y ]= NEWMARK ( ODEFUN , TSPAN , Y0 , NH , PARAM ) avec TSPAN =
% [ T0 TF ] intègre le système d ’ é q u atio ns différen -
% tielles y ’ ’=f (t ,y ,y ’) du temps T0 au temps TF avec
% la c o n dit ion initiale Y0 =(y ( t0 ) ,y ’( t0) en u t i lisan t
% la méthode de Newmark sur une grille de NH
% i n t e rva lles é q u i dis trib ués.
% PARAM contient les p a r a mèt res zeta et theta .
% La fonction ODEFUN (T , Y) doit r e t ourn er un vecteur
% c o n t enant les é v a l uatio ns de f (t ,y ) et de même
% d i m e nsion que Y . Chaque ligne de la solution Y
% c o r r espo nd à un temps contenu dans le vecteur
% colonne T.
tt= linspace ( tspan (1) ,tspan (2) , Nh +1);
y = y0 (:); u= y . ’;
global glob_h glob_t glob_y g l o b_ odefu n;
global g l o b _zeta g l o b_ theta g l o b _var argi n glob_fn ;
glob_h =( tspan (2) -tspan (1))/ Nh;
glob_y = y; g l o b _ode fun= odefun ;
g l o b_ze ta = param (1); g l o b _th eta = param (2);
g l o b _var argi n= varargin ;
if ( exist ( ’ O C T A VE _VE RSIO N’) )
o_ver = O C T A VE _VER SION;
version = str2num ([ o_ver (1) , o_ver (3) , o_ver (5)]);
end
if ( ~ exist ( ’ O C T A V E_VE RSIO N’ ) | version >= 320 )
options = optimset ;
options . Display = ’ off ’;
options . TolFun =1.e -12;
options . M a x F un Evals =10000;
end
glob_fn = feval ( odefun , tt (1) , glob_y , varargin {:});
for glob_t = tt (2: end)
if ( exist ( ’ O C T A VE_ VERS ION’ ) & version < 320 )
w = fsolve ( ’ n e w m arkfu n’ , glob_y );
else
w = fsolve ( @( w ) n e w m arkf un(w ) , glob_y , options );
end
glob_fn = feval ( odefun , glob_t ,w , varargin {:});
u = [ u; w . ’]; glob_y = w ;
end
t = tt;
clear glob_h glob_t glob_y g l o b _o defun;
clear g l o b_z eta g l o b _the ta g l o b _v arar gin glob_fn ;
end
function z= n e w m ar kfun( w)
global glob_h glob_t glob_y g l o b _ode fun;
global g l o b_ zeta g l o b _the ta g l o b _v arar gin glob_fn ;
fn1= feval ( glob_odefun , glob_t ,w , g l o b _var argi n{:});
z (1)= w (1) - glob_y (1) - glob_h * glob_y (2) -...
glob_h ^2*( g l o b_zet a* fn1 +(0.5 - g l o b_ze ta)* glob_fn );
z (2)= w (2) - glob_y (2) -...
glob_h *((1 -g l o b _thet a)* glob_fn + g l o b_ theta* fn1 );
end