390
10 R´ esolution num´ erique des ´ equations diff´ erentielles ordinaires
les p + 1 coefficients a i ; le vecteur colonne b qui contient les p + 2 coefficients
b i ; le pas de discr´ etisation h ; le vecteur des donn´ ees initiales u0 aux instants
t0 ; les macros fun et dfun contiennent les fonctions f et ∂f/∂y. Si la m´ ethode multi-pas est implicite, on doit fournir une tol´ erance tol et un nombre
maximal d’it´ erations itmax. Ces deux param` etres contrˆ olent la convergence
de l’algorithme de Newton utilis´ e pour r´ esoudre l’´ equation non lin´ eaire (10.46)
associ´ ee ` a la m´ ethode multi-pas. En sortie, le code renvoie les vecteurs u et t
qui contiennent la solution calcul´ ee aux instants t.
Programme 80 - multistep : M´ ethodes multi-pas lin´ eaires
function [t,u]=multistep(a,b,tf,t0,u0,h,fun,dfun,tol,itmax)
%MULTISTEP M´ ethode multi-pas.
% [T,U]=MULTISTEP(A,B,TF,T0,U0,H,FUN,DFUN,TOL,ITMAX) r´ esout le
% probl` eme de Cauchy Y’=FUN(T,Y) pour T dans ]T0,TF[ en utilisant une m´ ethode
% multi-pas avec les coefficients A et B. H est le pas de temps. TOL
% est la tol´ erance des it´ erations de point fixe quand la m´ ethode
% choisie est implicite.
y = u0; t = t0; f = eval (fun); p = length(a) - 1; u = u0;
nt = fix((tf - t0 (1) )/h);
for k = 1:nt
lu=length(u);
G=a’*u(lu:-1:lu-p)+ h*b(2:p+2)’*f(lu:-1:lu-p);
lt=length(t0);
t0=[t0; t0(lt)+h];
unew=u(lu);
t=t0(lt+1); err=tol+1; it=0;
while err>tol & it<=itmax
y=unew;
den=1-h*b(1)*eval(dfun);
fnew=eval(fun);
if den == 0
it=itmax+1;
else
it=it+1;
unew=unew-(unew-G-h*b(1)* fnew)/den;
err=abs(unew-y);
end
end
u=[u; unew]; f=[f; fnew];
end
t=t0;
return
Dans les prochaines sections, nous examinons quelques familles de m´ ethodes
multi-pas.