On désigne par I₀ un intervalle de ú non réduit à un point et t₀ un point fixé dans I₀ ; on se donne une fonction f définie et continue sur I₀ ´ ú ^m à valeurs dans ú ^m, et on cherche à trouver une fonction y continue et dérivable sur l’intervalle I₀, à valeurs dans ú ^m, telle que :

" t Î I₀, y¢ (t) = f(t,y(t)) (¤)

y(t₀) = y₀ [$]

Ce problème (E) s’appelle un problème de Cauchy pour le système différentiel (¤) ; la condition [$] s’appelle une condition de Cauchy. Une fonction y qui vérifie ces équations est appelée une intégrale de ce système différentiel.

Dans ce TP nous nous focaliserons sur les intervalles I₀ de la forme [t₀,t₀+T] qui seront une subdivision d’un intervalle [a,b]. Dans de nombreux exemples physiques, la variable t représente le temps ; l’instant t₀ est alors appelé instant initial et la condition [$] est appelé condition initiale.

L’équation (¤) est dite du premier ordre car elle ne fait intervenir qu’une dérivée première de la fonction y.

On désire trouver une approximation de la fonction y. On choisit des points x_i (º t_i) Î [a,b]. On cherche alors une approximation y_i de y(x_i). Il existe 2 cas :

• Le calcul de y_n fait intervenir les valeurs de x_n-1, y_n-1 et les données. Ce sont les méthodes à 1 pas.
• Le calcul de y_n fait intervenir les valeurs de :

Pour approcher (E), on choisit un ensemble fini de points x_i : a = x₀ < x₁ < ... < x_N = b, et on calcule les valeurs approchées y_i de y(x_i). En fait, pour plus de commodité dans les calculs, on imposera aux points x_i d’être équidistants. On définit le pas de l’approximation : h = (b – a) / N, alors x_i = a + i*h i = 0, 1, 2, ..., N

En x_n, la solution y de (E) vérifie : y(x_n+1) = y(x_n) + h.f(x_n,y(x_n)) + h.e (h) avec e ® 0 quand h® 0.

La méthode d’Euler explicite suggérée par cette relation consiste à approcher y(x_n) par la suite y_n donnée par :

La fonction eulerexpl calcule une solution approchée du problème de Cauchy défini plus haut par la méthode d’Euler explicite. Pour approcher ce système, nous considérons l’intervalle [a,b] qui va être subdiviser en N intervalles dont le pas est noté h. f est la fonction vérifiant de système différentiel (¤) 

Soit y0 une condition initiale du problème, l’instruction [X,y]=eulerexpl(‘f’,a,b,h,y0) renvoie le vecteur X contenant les points de la subdivision, ainsi que y – la solution approchée cherchée.

 

Interprétation graphique de cette méthode :

Le point A₀ est connu : A₀(a,a ).

On peut donc calculer la tangente à la courbe y(x) en ce point :

la pente est donnée par y¢ (a) = f(a,a ). L’approximation consiste à

confondre la courbe et sa tangente. On remplace le point inconnu

par A₁ point de même abscisse, mais situé sur la tangente T₀.

De proche en proche on construira A₁, A₂, …, A_N.

Pour gagner du temps, on peut remplacer le " feval " par la fonction elle-même soit en écrivant directement son expression calculée en (X(i),y(i,:)), soit en la définissant dans la fonction " eulerexpl " par un " inline " par exemple :

La figure 13 nous présente diverses représentations graphiques de solutions approchées d’équations différentielles du 1^er ordre par la méthode d’Euler explicite ainsi que la représentation de la véritable solution obtenu par calcul sur papier. Les systèmes sont les suivants :

J Le troisième système illustre une décharge de condensateur : un circuit électrique comprenant un condensateur de capacité C = 4.10^-4 F et une résistance ohmique R = 1000 W est fermé à l’instant t = 0 où la charge du condensateur est Q = 10 C. La décharge du condensateur est décrite par l’équation différentielle :

que l’on va réécrire sous la forme :  

J Enfin, le dernier exemple est le problème du parachutiste ou d’une chute en milieu résistant : un parachutiste de masse m = 70 kg est lancé à l’instant t = 0 avec une vitesse v₀ = 4 m/s, on cherche sa vitesse à chaque instant sachant que la résistance k = 10 N /m de l’air est proportionnelle à sa vitesse (on pose g = 9.81 N/m²).

En considérant le principe fondamental de la dynamique, l’équation différentielle est :

Voyons voir quelle peut être l’erreur maximale commise par cette méthode :

[X,y]=eulerexpl2(0,5,0.01,1); %Pour le 1^er exemple :

" [X,y]=eulerexpl2(pi/2,6*pi,0.01,1); %Pour le second exemple

" [X,y]=eulerexpl2(0,2,0.1,10); %Pour le 3^ème exemple

" [X,y]=eulerexpl2(0,40,0.1,4); %Pour le 4^ème exemple

Observons alors l’influence de h sur la représentation graphique de la solution approchée – c’est ce que traite la figure 14 – pour le système suivant :

[E1,E2,E3,E4] = 1.1973 _ 0.1444 _ 0.0147 _ 0.0015 %Erreur maximale dans chaque cas

On considère la relation : y(x_n) = y(x_n+1) – h y¢ (x_n+1) + he (h) avec e (h)® 0 quand h® 0. Un raisonnement analogue au précédent conduit à approcher y(x_n) par la suite y_n solution de :

Cette méthode est moins simple d’utilisation que la première puisqu’elle fait appel en son sein à la méthode du point fixe ; l’approximation y_n étant connue, il n’est pas clair que ce système définisse y_n+1 de manière unique – nous l’admettrons cependant (condition précisée au § 1.2 du TP).

La fonction eulerimpl3 calcule une solution approchée du problème de Cauchy par la méthode d’Euler implicite. Pour approcher ce système, nous considérons l’intervalle [a,b] qui va être subdiviser en N intervalles dont le pas est noté h. Dans le calcul de la solution, la fonction fait appel à la méthode du point fixe dont une fonction notée pointfix lui est associée.

Soient y0 une condition initiale du problème, eps la précision demandée sur 2 y_i – y_i-12 ₂ et nmax le nombre maximale d’itérations intervenant dans la boucle du point fixe , l’instruction [X,y,erreurpf]= eulerimpl3 (a,b,h,y0,eps,nmax) renvoie le vecteur X contenant les points de la subdivision, y – la solution approchée cherchée et erreurpf : erreur à chaque pas sur la solution approchée de Y=G(Y) par la méthode du point fixe..

%eps : précision demandée sur ½ y_i – y_i-1½

%nmax : nombre maximal d’itérations pour la méthode du point fixe

%erreurpf : on a soit erreur £ eps ou i ³ nmax

X=a:h:b; %Vecteur contenant la subdivision de [a,b]

N=length(X);

g=inline(‘2*y’,‘x’,‘y’); %par exemple

y(1,:)=y0; %le 1^er itéré correspond à la condition initiale

for n=1:N-1,

y1=y(n,:); %calcul du i^ème itéré en faisant intervenir la boucle du

end;

y2=y3+h*feval(nameg,X(n+1),y1); %|y2-y1|® eps ou que i=nmax

erreurpf=norm(y2-y1,2);

end

Y=y2; %le (n+1)^ème itéré est y2

La figure 15 reprend les mêmes exemples que la figure 13. Les résultats sont assez identiques. Nous pouvons cependant noter une défaillance de cette méthode pour le cas 2 avec le système différentiel :

 En effet, sur l’intervalle [p /2,p [ la solution approchée coïncide très bien avec la véritable solution, en revanche dès qu’on se situe dans un intervalle pour x ³ p , la solution n’existe pas. Cette méthode accepte peut être mal le passage de ln y = ln sin(x) + c à y = d*sin(x) ou bien la définition même de la fonction G(x,y) = y*cotan(x). Pour preuve, voyons la feuille de commande suivante :

" length(y1) % y1 : solution approchée par eulerimpl3 du système précédent

ans = Les valeurs tombent brutalement à - ¥ …

-Inf i=158

ans =

-Inf

ans =

-Inf

ans =

-Inf

ans =

-Inf

ans =

-Inf

" y1(1000:1002)

ans =

NaN …pour devenir inexistantes (Not-a-Number) : il y a une forme indéterminée

NaN

" [E1,E2,E3,E4] %Ei correspondant au cas i

" load l %Pour mémoire avec la méthode d’Euler explicite nous avions trouvé :

On constate que cette méthode est un peu plus précise qu’Euler explicite. Nous pouvons faire également le même genre de remarque quant à l’influence du pas – la figure 16 nous le montre. Cette figure permet aussi de constater que la méthode implicite et la méthode explicite encadre de chaque côté la courbe de la véritable solution. Ici – la méthode implicite est une approximation par excès (courbe au dessus de la courbe " vraie ") alors que l’explicite est une approximation par défaut. Mais leur rôle peut s’inverser suivant le système différentiel à résoudre.

Pour q Î [0,1] le schéma est le suivant :

Pour q = 0, on retrouve la méthode d’Euler explicite, et pour q = 1, la méthode d’Euler implicite. Pour q > 0, on a une équation de point fixe à résoudre.

La fonction teta calcule une solution approchée du problème de Cauchy par la q - méthode. Pour approcher le système de Cauchy, nous considérons l’intervalle [a,b] qui va être subdiviser en N intervalles dont le pas est noté h. Dans le calcul de la solution, la fonction fait appel à la méthode du point fixe dont une fonction notée pointfixe lui est associée.

Soient t une valeur de 0£ q £ 1, y0 une condition initiale du problème, eps la précision demandée sur 2 y_i – y_i-12 ₂ et nmax le nombre maximal d’itérations intervenant dans la boucle du point fixe , l’instruction [X,y]= teta (a,b,h,y0,t,eps,nmax) renvoie le vecteur X contenant les points de la subdivision, ainsi que y – la solution approchée cherchée.

%fonction donnant une solution approchée d'une EDO par la q -méthode

function [X,y] = teta (a,b,h,y0,t,eps,nmax); %0£ t£ 1

X=a:h:b; %vecteur contenant les points de équidistants de la subdivision

N=length(X);

g=inline('2*y','x','y'); %g : fonction vérifiant le syst. diff.

y(1,:)=y0; %1^er itéré y’(x)=g(x,y(x))

for n=1:N-1, %calcul du (n+1)^ème itéré

y1=y(n,:);

[Y,erreurpf]=pointfixe(g,h,X,n,y1,t,eps,nmax); %h:pas de la subdivision

y(n+1,:)=Y;

end;

En guise d’illustration nous avons comparer, d’abord, les temps de calculs, ainsi que le nombre d’opérations effectuées par les différentes méthodes vues jusqu’ici dans ce TP. Nous pouvons voir que la q -méthode est beaucoup plus longue que les autres du fait du nombre d’opérations à effectuer. Mais cet inconvénient devient vite un avantage si on recherche la précision du calcul. A la fin de cette session de Matlab, nous comparons les erreurs maximales commises entre la q –méthode et les formules d’Euler (E1 et E3 : elles coincident parfaitement – cf. figure 17) et entre la q -méthode et la véritable solution (E2 et E4).

Le système à résoudre est :

On enlèverait la valeur absolue dans le calcul d’erreur, on verrait que E2 = - E4, i.e. que la q - méthode enveloppe la vraie solution.

Voyons pour finir l’influence du q sur le calcul de la solution du système précédent : cf. figure 18.

Nous pouvons affirmer que la q -méthode permet une approximation plus fine que les méthodes d’Euler. En effet, il faut d’abord remarquer que la vraie solution du problème se trouve dans l’espace délimité par les solutions des 2 méthodes d’Euler. D’autre part, suivant la valeur de q , la solution obtenue par la q -méthode balaye cet espace, donc il existe un q (» 0.5) pour lequel cette solution coïncide avec la solution vraie.

Si la matrice A est strictement triangulaire inférieure, i.e. si i£ j Þ a_ij = 0, la résolution du système est immédiate : on a y_n,1 = y_n, on obtient y_n,2 à partir de y_n et y_n,1, y_n,3 à partir de y_n, y_n,1, y_n,2 etc... On dit dans ce cas que la méthode de Runge-Kutta est explicite : c’est le cas des méthodes données aux exemples 1, 2 & 3 du TP.

Si la matrice A est triangulaire inférieure, i.e. si i<j Þ a_ij = 0, la résolution se ramène à résoudre successivement q fois une équation à une inconnue, la méthode de Runge-Kutta est dite alors semi-implicite ; numériquement la résolution est plus facile que dans le cas général où on a un système de q équations à q inconnues. Dans le cas général, on parle de méthode de Runge-Kutta implicite.

La fonction rk44 calcule une solution approchée du problème de Cauchy par la méthode de Runge-Kutta à l’ordre 4, c’est la méthode classique. Pour approcher le système différentiel, nous considérons l’intervalle [a,b] qui va être subdiviser en N intervalles dont le pas constant est noté h.

Soient y0 une condition initiale et R le tableau donné en énoncé dans le TP qui nous livrera la matrice A et les 2 vecteurs B et C , l’instruction [X,Y]=rk44(‘f’,a,b,y0,h,R) nous livre, en plus du vecteur X, la solution Y cherchée. f est la fonction vérifiant le système différentiel considéré.

Le tableau est donné dans le TP – mais précisons la suite, qui est déterminée d’après le tableau :

Nous supposons l’équation différentielle linéaire : y¢ =2y, y(0)=1 sur [0,1]. Le pas de discrétisation est choisi constant – noté h : nous avons les résultats suivant :

" a=0;b=1;h=0.1;y0=1; %définition des paramètres

" flops(0);tic,[x,y1]=eulerexpl2(a,b,h,y0);toc,disp(flops); %EULER EXPLICITE

" %flops nous donne le nombre d'opérations effectuées par la fonction

" flops(0);tic,[x,y2]=eulerimpl3(a,b,h,y0,0.00001,1000);toc,disp(flops); %EULER IMPLICTE

" flops(0);tic,[x,y3]=teta(a,b,h,y0,0.5,0.00001,1000);toc,disp(flops); %q -METHODE (q =0.5)

" flops(0);tic,[x,y4]=rk2('f',a,b,y0,h);toc,disp(flops); %Euler modifié

" flops(0);tic,[x,y5]=rk2('f',a,b,y0,h);toc,disp(flops); %Heun

" flops(0);tic,[x,y6]=rk2('f',a,b,y0,h);toc,disp(flops); %RK2 (a =1.5)

" flops(0);tic,[x,y7]=rk44('f',a,b,y0,h,R);toc,disp(flops); %RK4 CLASSIQUE

" E1=0;for i=1:length(x), if e1(i)>E1, E1=e1(i);end;end; % euler explicite

" E2=0;for i=1:length(x), if e2(i)>E2, E2=e2(i);end;end; % euler implicite

" E3=0;for i=1:length(x), if e3(i)>E3, E3=e3(i);end;end; % q -méthode avec q = 0.5

" E4=0;for i=1:length(x), if e4(i)>E4, E4=e4(i);end;end; % RK2 avec a = 0.5

" E5=0;for i=1:length(x), if e5(i)>E5, E5=e5(i);end;end; % RK2 avec a = 1

" E6=0;for i=1:length(x), if e6(i)>E6, E6=e6(i);end;end; % RK2 avec a = 1.5

" E7=0;for i=1:length(x), if e7(i)>E7, E7=e7(i);end;end; % RK4

"

" [E1 E2 E3 E4 E5 E6 E7] % h = 0.1

Au niveau du temps de calcul, Euler explicite et Runge – Kutta à l’ordre 2 sont très rapide d’autant plus que Runge – Kutta est assez précise dans l’approximation de la solution. Nous pouvons d’ailleurs constater que dans la méthode de RK2 le résultat est indépendant de a . Nous concluons également que la méthode RK4 est la méthode la plus précise pour un coût pas trop élevé comparé à Euler implicite et la q -méthode.

Reprenons ces calculs pour un pas de h = 0.01 :

" [E1 E2 E3 E4 E5 E6 E7] % h = 0.01

" [E1 E2 E3 E4 E5 E6 E7] % h = 0.0001

Il modélise des réactions chimiques oscillantes. La variable y est définie dans ú ², c’est donc un couple [y₁,y₂] dont chaque composante dépend de x : on écrit : y(x) = [y₁(x),y₂(x)]. Le problème à résoudre possède donc 2 équations et 2 conditions initiales :

Nous allons définir une fonction brusselator(x,y) définie de ú ² dans ú ²  et calculons une solution approchée du système par 4 méthodes : Euler explicite (eulerexpld2), Euler implicite (eulerimpld), Runge-Kutta classique ( rk44d) et enfin par un solveur Matlab nommé ODE45.

La fonction eulerexpld2 résout le système posé précédemment par la méthode d’Euler explicite. L’instruction [X,y1,y2]=eulerexpld2(‘brusselator’,0,20,0.01,[1 3]) nous livre le vecteur X contenant les points de la subdivision de [a,b] dont le pas h constant est pris égal à 0.01 ainsi que les 2 composantes du vecteur y – solution approchée du problème. y0 = [1 3] est la condition initiale du problème.

La fonction eulerimpld résout le problème par la méthode d’Euler implicite. Soient eps l’erreur relative entre 2 termes consécutifs dans la méthode du point fixe et nmax le nombre maximal d’itérations, l’instruction [X,y1,y2]=eulerimpld(0,20,0.01,[1 3],eps,nmax) nous livre le vecteur X contenant les points de la subdivision de [a,b] dont le pas h constant vaut également 0.01 ainsi que les 2 composantes du vecteur y – solution approchée du problème.

y0 = [1 3] est la condition initiale du problème. C’est un vecteur à 2 dimensions car le mouvement est contenu dans un plan.

 Algorithme de la fonction eulerimpld :

La fonction rk44d résout le problème par la méthode de Runge-Kutta classique. Soit R le tableau contenant la matrice A et les 2 vecteurs B et C définis dans le paragraphe sur les méthodes de Runge-Kutta, l’instruction [X,y1,y2]=rk44d(‘brusselator’,0,20,[1 3],0.01,R) nous livre le vecteur X contenant les points de la subdivision de [a,b] dont le pas h constant vaut également 0.01 ainsi que les 2 composantes du vecteur y – solution approchée du problème.

Nous avons choisi l’ordre 4 pour cette méthode car d’après les conclusions déjà établies, RK4 nous fournit une assez bonne approximation de la solution.

La figure 19 affiche les représentations graphiques des composantes de y obtenues par les diverses méthodes précédentes. Nous avons opté pour un pas de 0.01 car d’une part les temps de calculs ne sont pas trop long et d’autre part les résultats sont assez acceptables (cf. page suivante). La fonction ode45 de Matlab calcule automatiquement le pas – nous pouvons donc constater d’après la forme des courbes que le pas choisi est proche du pas calculé par ode45.

" flops(0);tic,eulerexpld2('brusselator',0,20,0.01,[1 3]);toc,disp(flops) %Euler explicite

" flops(0);tic,eulerimpld(0,20,0.01,[1 3],0.0001,100);toc,disp(flops) %Euler implicite

" flops(0);tic,rk44d('brusselator',0,20,[1 3],0.01,R);toc,disp(flops) %Runge-Kutta classique

" flops(0);tic,[t,y]=ode45('brusselator',[0,20],[1 3]);toc,disp(flops) %ode45

Regardons le pas d’ode45 :

" length(t)

ans =

213

" t(15)-t(14),t(100)-t(99),t(213)-t(212)

ans =

0.2212

ans =

0.1383

ans =

0.1643

Finalement, le pas d’ode45 n’est pas réellement proche de h=0.01. Nous pouvons alors penser que pour l’approximation de la solution de ce problème, il ne suffit que d’une grossière subdivision de [0 20] puisque apparemment h_ode45 << 0.01. En traçant la courbe, paramétrée par les 2 composantes de la solution trouvée par quelques solveurs, nous constatons cependant que pour h=0.01, on se rapproche bien de la solution d’ode45.

C’est ce qu’indique la figure suivante :

Le mouvement d’un corps céleste de masse négligeable dans un champ de gravité de 2 corps en rotation circulaire l’un autour de l’autre, est décrit par les équations figurant sur le TP. Ce problème modélise la trajectoire d’un vaisseau spatial allant de la Terre à la Lune . Nous sommes en présence d’un problème dont l’équation peut s’écrire dans sa globalité sous la forme : Y² = f(Y,Y¢ ). Il faut donc se ramener à un problème du type Z¢ = g(Z). Pour cela nous allons écrire les vecteurs Z et Z¢ de la manière suivante :

Nous venons de définir une fonction celeste qui associe à un vecteur de ú ⁴ un autre vecteur de ú ⁴. Il faut également précisé que si Z¢ =celeste(Z), il faut le comprendre comme Z¢ (x)=celeste(Z(x)). Son algorithme est le suivant :

La résolution du système peut se faire grâce au solveur ode45 dont la syntaxe est la suivante :

" xn=11.124340337266085135; %borne sup de l’intervalle [a,b]

" plot(z(:,1),z(:,2)) %courbe paramétrée

Le résultat nous est fourni par la figure 20. Nous voyons clairement que l’orbite est périodique : la fusée par de la Terre en direction de la Lune, tourne autour d’elle 2 fois de suite mais pendant ces 2 tours le satellite naturel de la Terre se déplace – c’est ce qu’indique les 2 positions marquées " Lune " sur la courbe – puis regagne la Terre au même point car la durée de rotation de la Terre est beaucoup plus longue que celle de la Lune du fait de son diamètre.

Nous pouvons également utiliser d’autres solveurs déjà vus dans ce TP : Euler explicite, Runge-Kutta classique, par exemple. Contrairement à ode45 qui calcule automatiquement le pas, nous devons procéder par approches successives pour obtenir une orbite périodique. Les figures 20 bis et 20 ter, montre l’influence du pas, donc du nombre de points de la subdivision [0 x_n] sur l’orbite cherchée. Nous sommes amenés à constater que plus h diminue, plus on se rapproche d’une trajectoire fermée.

Nous avons modifier nos fonctions de bases pour essayer d’améliorer la rapidité des calculs.

Affichons simplement les algorithmes :

Nous voyons nettement que ode45 est beaucoup plus rapide que les 2 autres fonctions. De plus, ode45 calcule un pas non constant, ce qui permet une meilleure approximation des solutions, puisque cette fonction peut rapprocher les points de la subdivision ou les éloigner suivant les variations de la fonction.