Méthode d'Euler implicite

1. Définition et stabilité

La méthode d'Euler est présentée dans Méthode d'Euler explicite. On reprend ici les mêmes notations. La méthode d'Euler implicite consiste à chercher la valeur approchée à l'instant t_n+1 avec la relation suivante :

$$y_{n+1}=y_n+hf(t_{n+1},y_{n+1})\text{(1)}$$

Cette méthode consiste donc à prendre la dérivée à la fin de l'intervalle [t_n,t_n+1] au lieu de la prendre au début.

La valeur y_n+1 est obtenue par résolution d'une équation. Dans le cas d'un système différentiel, le nombre d'équation à résoudre est égal aux nombres de fonctions inconnues.

Pour étudier la stabilité, on considère l'équation linéaire suivante :

$$\frac{dy}{dt}=qy\left(t\right)\text{(2)}$$

Pour la condition initiale y(0)=1, la solution est :

$$y\left(t\right)=exp\left(qt\right)\text{(3)}$$

Pour cette équation, la méthode d'Euler implicite s'écrit :

$$y_{n+1}=y_n+hq{y}_{n+1}\text{(4)}$$

Les valeurs approchées obtenues par cette méthode sont donc :

$$y_n=\frac{y_0}{{(1-hq)}^n}\text{(5)}$$

On pose z=hq. Le domaine de stabilité est défini comme le domaine du plan complexe où la suite précédente est bornée, c'est-à-dire :

$$|1-z|\ge 1\text{(6)}$$

Le domaine de stabilité est donc tout le plan complexe à l'exclusion du disque de centre z=1 et de rayon 1. Le demi espace Re(z)<0 est donc entièrement dans le domaine de stabilité. Si la partie réelle de q est négative, la solution approchée à toujours le même comportement que la solution exacte, c'est-à-dire qu'elle tend vers 0 lorsque t tend vers l'infini.

2. Itération de Newton-Raphson

L'équation (1) est linéaire pour les systèmes différentiels linéaires. Dans le cas général, il faut résoudre un système d'équations algébriques non linéaires. Celui peut être résolu par la méthode de Newton-Raphson décrite ci-dessous.

Pour préciser les calculs à effectuer, nous écrivons explicitement les équations différentielles. Soit P le nombre d'équations. L'équation différentielle k est :

$$\frac{d{y}^k}{dt}=f^k(t,y^1\left(t\right),\dots ,y^i\left(t\right),\dots ,y^P\left(t\right))\text{(7)}$$

L'indice désignant la fonction inconnue est donc en exposant. La méthode d'Euler implicite s'écrit alors, pour $1\le k\le P$:

$$y_{n+1}^k=y_n^k+h{f}^k(t_{n+1},y_{n+1}^1,\dots ,y_{n+1}^i,\dots ,y_{n+1}^P)\text{(8)}$$

Posons

$$G^k(y_{n+1}^1,\dots ,y_{n+1}^i,\dots ,y_{n+1}^P)=y_n^k+h{f}^k(t_{n+1},y_{n+1}^1,\dots ,y_{n+1}^i,\dots ,y_{n+1}^P)-y_{n+1}^k\text{(9)}$$

Le système à résoudre pour déterminer les $y_{n+1}^k$ s'écrit alors :

$$G^k(y_{n+1}^1,\dots ,y_{n+1}^i,\dots ,y_{n+1}^P)=0\text{(10)}$$

Soit Y_n+1 la matrice colonne contenant les $y_{n+1}^k$. Une première estimation de Y_n+1 est donnée par Y_n. Pour simplifier la notation, on note Y cette estimation. La méthode de Newton-Raphson consiste à considérer le développement de Taylor à l'ordre 1 suivant :

$$G^k(Y+\delta Y)=G^k\left(Y\right)+\sum _{j=1}^P\frac{\partial G^k}{\partial y^j}\delta y^j\text{(11)}$$

L'estimation suivante est donnée par $Y\text{'}=Y+\delta Y$ avec :

$$\sum _{j=1}^P\frac{\partial G^k}{\partial y^j}\delta y^j=-G^k\left(Y\right)\text{(12)}$$

Il faut donc résoudre un système linéaire pour calculer les $\delta Y$. Soit J le jacobien du système différentiel évalué en Y, c'est-à-dire la matrice carrée PxP définie par :

$$J_{ij}=\frac{\partial f^i}{\partial y^j}\left(Y\right)\text{(13)}$$

Le système à résoudre s'écrit sous forme matricielle :

$$(hJ-I_d)\delta Y=-G\text{(14)}$$

Ce système est résolu avec la décomposition LU puis Y' est calculé. Si $\left|\delta Y\right|$ est supérieure à une tolérance, l'itération est reprise.

3. Programmation avec Scilab

Pour tester le fonctionnement, on considère l'équation différentielle de l'oscillateur harmonique, qui se met sous la forme d'un système différentiel du premier ordre :

$$\begin{array}{ccc}\hfill & \frac{d{y}_1}{dt}=y_2\hfill & \hfill \text{(15)}\\ \hfill & \frac{d{y}_2}{dt}=-4{\pi }^2{y}_1\hfill & \hfill \text{(16)}\end{array}$$

Le système différentiel est défini par une fonction :

function [dy]=oscillateur(t,y),
 dy(1)=y(2),
 dy(2)=-4*%pi^2*y(1),
endfunction
   

Rappelons tout d'abord les fonctions réalisant la méthode d'Euler explicite :

function [y]=pasEuler(systeme,ne,h,t,yn),
   dy = systeme(t,yn);
   for e=1:ne,
    y(e) = yn(e)+h*dy(e)
   end;
endfunction

function [mt,my]=euler(yi,T,N,systeme,ne),
   h = T/N;
   my = zeros(ne,N+1);
   mt = zeros(1,N+1);
   my(:,1)=yi;
   mt(1)=0;
   for k=1:N,
    my(:,k+1)=pasEuler(systeme,ne,h,mt(k),my(:,k));
    mt(k+1) = mt(k)+h;
   end;
endfunction
   

Pour la méthode implicite, on doit aussi définir le jacobien :

function [J]=jacobienOscillateur(t,y),
   J(1,1) = 0; J(1,2)=1;
   J(2,1)=-4*%pi^2; J(2,2)=0;
endfunction   
   

La fonction suivante effectue le calcul de y_n+1 à partir de y_n :

function [y]=pasEulerImplicite(systeme,jacobien,Id,P,h,t,yn,tol),
   stop = 0; it=0; maxit=100;
   y = yn;
   t1 = t+h;
   while stop==0,
    it = it+1;
    dy = systeme(t1,y);
    G = yn+h*dy-y;
    [deltaY,ker] = linsolve(h*jacobien(t1,y)-Id,G);
    y = y+deltaY;
    if it>maxit then stop=1; end;
    stop = 1;
    for k=1:P,
     if abs(deltaY(k))>tol then stop=0; end;
    end;
   end;
endfunction
   

La fonction suivante effectue le calcul complet, avec une condition initiale (instant 0) définie dans une matrice colonne.

function [mt,my]=eulerImplicite(yi,T,N,systeme,jacobien,P,tol),
   Id = eye(P,P);
   h = T/N;
   my = zeros(P,N+1); 
   mt = zeros(1,N+1);
   my(:,1)=yi;
   mt(1)=0;
   for k=1:N,
    my(:,k+1)=pasEulerImplicite(systeme,jacobien,Id,P,h,mt(k),my(:,k),tol);
    mt(k+1) = mt(k)+h;
   end;
endfunction
   

Voyons le fonctionnement avec le système défini plus haut :

yi=[1;0];
P=2;
N=5000;
T=20;
tol=1e-5;
[mt,my] = eulerImplicite(yi,T,N,oscillateur,jacobienOscillateur,P,tol);
   

figA=scf();
plot2d(mt,my(1,:));
xtitle('','t','y(1)');
   

figA.pdf

La solution exacte est une sinusoïde d'amplitude 1. Comme prévu par l'étude de stabilité, l'amplitude d'oscillation décroît alors qu'elle augmente avec la méthode d'Euler explicite. En réduisant le pas d'un facteur 2, on réduit l'erreur du même facteur (méthode d'ordre 1) :

yi=[1;0];
N=10000;
[mt,my] = eulerImplicite(yi,T,N,oscillateur,jacobienOscillateur,P,tol);
   

figB=scf();
plot2d(mt,my(1,:));
xtitle('','t','y(1)');
   

figB.pdf

Bien entendu, la méthode d'Euler est en pratique inutilisable en raison du coût très élevé en calcul pour obtenir une erreur faible. On voit néanmoins sur cet exemple une propriété vérifiée ausi sur des méthodes d'ordre plus élevé : l'application de la méthode à un système conservatif revient à introduire une dissipation dans le système. Voyons ce qu'il en est pour un système dissipatif, l'oscillateur amorti. Le calcul est aussi conduit avec la méthode d'Euler explicite.

function [dy]=amorti(t,y),
 dy(1)=y(2),
 dy(2)=-4*%pi^2*y(1)-0.1*y(2),
endfunction

function [J]=jacobienAmorti(t,y),
 J(1,1)=0; J(1,2)=1;
 J(2,1)=-4*%pi^2; J(2,2)=-0.1;
endfunction

P=2;
N=10000;
T=40;
tol=1e-5;
[mt,my] = eulerImplicite(yi,T,N,amorti,jacobienAmorti,P,tol);
function y=f(t),
 y=exp(-0.05*t)*(cos(2*%pi*t)+0.05/(2*%pi)*sin(2*%pi*t))
endfunction
[mt1,my1] = euler(yi,T,N,amorti,P);
   

figC=scf();
plot2d(mt,my(1,:),style=color('red'));
plot2d(mt1,my1(1,:),style=color('blue'))
fplot2d(mt,f)
xtitle('','t','y(1)');
   

figC.pdf

Pour comparaison, la solution exacte est tracée (en noir) et la solution obtenue par la méthode d'Euler explicite (en bleu). Sur cet exemple de système dissipatif, la méthode d'Euler implicite conduit à une erreur plus faible que la méthode explicite, et surtout à un comportement beaucoup plus réaliste.