Méthode numériques équations différentielles

République Algérienne Démocratique et Populaire Ecole Normale Supérieurs d’Enseignement Technique -ORANDépartement de Génie Electrique

Magistère 1ere année en électrotechnique Option : Analyse et Commande des Machines Electrique

Monographie de Méthode numérique

Méthode de résolution des équations différentielles ODE

Mr BOUZID Allal El Moubarek

BELAIDI Sous la direction de : Mr. BELAIDI

Anne universitaire 2008-2009

SOMMAIRE Introduction générale

01

Chapitre I : Résolution des ODE Introduction

02

I.1 Définition

02

I.2 Rappel des solutions EDO simple

04

I.3 Le problème de Cauchy

04

I.4 Principe général des méthodes numériques

05

I.5 Propriétés des méthodes numériques

05

I.6 Les principales méthodes numériques

07

Chapitre II Méthodes a un pas Introduction

09

II.1 Méthodes d'Euler explicite et implicite

09

II.2 Méthode d'Euler amélioré

15

II.3 Méthode d'Euler-Cauchy

15

II.4 Méthode de Crank – Nicholson

18

II.5 Méthode de Heun

19

II.6.Méthode de Runge-Kutta

20

II.6.1 Méthode de Rungue Kutta 2

20

II.6.2 Méthode de Runge et Kutta explicite a un point

21

II.6.3 Méthodes RK implicites à un point

22

II.6.4 Méthodes RK explicites à 2 points intermédiaires

22

II.6.5 Méthodes RK explicites à 3 points intermédiaires

22

II.6.6 Formule générale des méthodes RK explicites

24

II.7. Méthode de Burlish – Stoer

26

II.8 Comparaison entre quelques méthodes

27

II.9 Comparaison entre toutes les méthodes

29

II.10. La méthode de Runge-Kutta-Fehlberg

32

Conclusion

38

Chapitre III Méthode à pas multiples Introduction

40

III-1 Méthode Adams–Bashforth–Moulton Method

40

III.2.Méthode de Hamming

44

III.3.Comparaison des méthodes (RK4 ; Hamming, ABM)

46

III.4 Méthode de Milne Simpson

51

III.5. Méthodes de Gear

52

III.5.1 Méthode de Gear à 2 pas, implicite, d'ordre 2

52


52


52

III.6 problème de la valeur marginale

52

III.6.1 Méthode de Tir

53

III.7 Méthode des différences finis (Finite Difference)

Bibliographie

56

61

Introduction

Introduction

Introduction On regroupe sous le terme générique de « méthodes numériques », toutes les techniques de calcul qui permettent de résoudre de manière exacte ou, le plus souvent, de manière approchée un problème donné. Le concept de calcul est assez vaste et doit être pris au sens large. Il peut s'agir de déterminer l'inconnue d'une équation, de calculer la valeur d'une fonction en un point ou sur un intervalle, d'intégrer une fonction, d'inverser une matrice, etc. Bien que la mise en équation d'un problème et sa résolution passent naturellement par les Mathématiques, les problématiques sous-jacentes concernent des disciplines aussi variées que la Physique, l'Astrophysique, la Biologie, la Médecine, l'Economie, etc. Il existe ainsi une grande variété de problèmes possibles avec pour chacun d'eux, des méthodes très spécifiques. De fait, le nombre total de méthodes numériques dont nous disposons à l'heure actuelle est vraisemblablement gigantesque. L’ordinateur est aujourd’hui un outil incontournable pour simuler et modéliser les systèmes, mais il faut encore savoir exprimer nos problèmes en langage formalisé des mathématiques pures. Nous sommes habitués à résoudre les problèmes de façon analytique, alors que l’ordinateur ne travaille que sur des suites de nombres. On verra dès lors qu’il existe souvent plusieurs approches pour résoudre un même problème, ce qui conduit à des algorithmes différents. Un algorithme, pour être utile, doit satisfaire un certain nombre de conditions. Il doit être : • Rapide : le nombre d’opérations de calcul pour arriver au résultat escompté doit être aussi réduit que possible. • Précis : l’algorithme doit savoir contenir les effets des erreurs qui sont inhérentes à tout calcul numérique. Ces erreurs peuvent être dues à la modélisation, à la représentation sur ordinateur ou encore à la troncature. • Souple : l’algorithme doit être facilement transposable à des problèmes différents. La résolution numérique d’équations différentielles est très souvent nécessaire, faute de l’existence de solutions analytiques.

-1-

Chapitre I Définition des équations différentielles

CHAPITRE I

Définition des équations différentielles

Introduction Les équations différentielles ordinaires (EDO) apparaissent très souvent dans la modélisation de la physique et des sciences de l'ingénieur. Trouver la solution d'une EDO ou d'un système d'EDO est ainsi un problème courant, souvent difficile ou impossible à résoudre de façon analytique. Il est alors nécessaire de recourir à des méthodes numériques pour résoudre ces équations différentielles. I.I Définitions : Équation d'un (premier ordre) de la forme y' = f (x, y) exprime le taux de changement de la variable dépendante y en ce qui concerne un changement de la variable indépendante X comme fonction f (x, y) de la variable indépendante X et de la variable dépendante Y. Une solution de l'équation y' ; = f (x, y) est une fonction, que par abus de langue nous dénotons comme y (x), tel que y' ; (x) = f (x, y (x)). Le théorème d'existence et d'unicité des équations ordinaires du premier ordre déclare que si la fonction f est continue et remplit la condition de Lipschitz sur une certaine région fermée U×V, avec x dans U et y dans V, contenant le point intérieur (x0, y0), alors là existe une solution unique y (x) et un sous-intervalle U0 contenant x0 et contenu dans U tels que y (x0) = y0 et (x, y (x)) est contenu dans U×V pour x dans U0. La condition que y (x0) = y0 s'appelle l'état initial. Si la fonction f est différentiable dans la région fermée, alors la condition de Lipschitz est automatiquement satisfaite. En outre, si la fonction f est Cn, c.-à-d. f a les dérivés continus jusqu'à la commande n, pour n > 0, dans la région fermée, alors la solution est Cn+1 dans un voisinage de x0. Notez que si l'état de Lipschitz n'est pas satisfaisant, alors il n'existe pas une solution ou il peut exister plus d'une solution. Si la fonction f (x, y) est indépendant de la variable dépendante y, alors la solution est y (x) = y0 + ∫ x0 f (t), où l'intégrale peut être rapprochée par une des techniques numériques dans l'intégration numérique de section. Si la fonction f (x, y) est indépendant de la variable indépendante X, l'équation y' ; =f (y) serait une équation autonome. Les méthodes données ci-dessous s'appliquent aux équations autonomes et non autonomes mais dans le cas autonome, l'utilisateur peut vouloir éditer les dossiers et supprimer des références à la variable indépendante partout où la variable indépendante est mise en référence par le F. Système de premier ordre des équations - équation évoluée : La variable dépendante y peut être une quantité y de vecteur dans ce cas nous avons un système des équations du premier ordre. Le théorème d'existence et d'unicité indiqué ci-dessus généralise immédiatement pour inclure la situation dans laquelle y est un vecteur dimensionnel de n. Une équation y de nième ordre (n) = f (x, y, y' ; , y (2),…, y (n-1)) peut être transformé à un système des équations du premier ordre en définissant le z= de quantité de vecteur (z1,…, Zn) par z1=y, z2=y' ; , z3=y (2),…, zn=y (n-1), puis z' ; 1=z2, z' ; 2=z3,…, z' ; n-1=zn, et z' ; n=f (x, z1,…, Zn).

ENSET Oran

-2-

CHAPITRE I


Tandis qu'une équation du second degré peut être transformé à un système de premier ordre comme décrit ci-dessus mais parce que les équations du second degré sont omniprésentes dans des méthodes spéciales de physique et de technologie ont été développés pour les résoudre. Soit une fonction numérique notée y(x) définie sur un intervalle de R et de classe Cp (continûment dérivable d’ordre p). On appelle équation différentielle d’ordre p une équation de la forme F(x, y, y′, y′′, . . . y(p)) = 0 (1.1) Ou y′ représente la dérivée première par rapport x, y′′ la dérivée seconde, etc... Plus généralement, on appelle système différentiel un ensemble d’équations différentielles reliant une variable x et un certain nombre de fonction yi(x) ainsi que leurs dérivées. L’ordre du système différentiel correspond à l’ordre de dérivation le plus élevé parmi l’ensemble des fonctions. On appelle solution de l’équation différentielle (1.1) toute fonction y(x) de classe Cp qui vérifie l’équation (1.1). On appelle forme canonique d’une équation différentielle une expression du type y(p) = f(x, y, y′, y′′, . . . y(p−1)) (1.2)

ENSET Oran

-3-

CHAPITRE I


I.2 RAPPEL - SOLUTIONS D'EDO SIMPLES Nous rappelons ici les principes de résolution des équations différentielles simples. Dans ce qui suit, la variable indépendante « t » sera notée et la variable dépendante « x(t) ».

ENSET Oran

-4-

CHAPITRE I


I.3 LE PROBLEME DE CAUCHY On appelle problème de Cauchy ou problème à la valeur initiale le problème qui consiste à trouver une fonction y(x) définie sur l'intervalle [a; b] telle que :

y ' ( x)  f ( x, y ( x))

;

x  [a, b]

y ( a )  y0 Si la fonction f est continue et vérifie une condition de Lipschitz par rapport à la deuxième variable alors le problème admet une solution unique. On dit que le problème est bien posé. Attention : un problème bien posé ne signifie pas qu'il peut se résoudre numériquement ! I.4 PRINCIPE GENERAL DES METHODES NUMERIQUES Pour obtenir une approximation numérique de la solution y(x) sur l'intervalle [a; b], nous allons estimer la valeur de cette fonction en un nombre fini de points xi, pour i = 0; 1; :::; n, constituants les noeuds du maillage. La solution numérique discrète obtenue aux points xi est notée : yi = y (xi). L'écart entre deux abscisses, noté h, est appelé pas de discrétisation. Ce pas, dans les méthodes les plus simples, est constant, mais il peut être judicieux de travailler avec un pas variable hi = xi – xi-1. Le choix du maillage et de la répartition des nœuds, peut s'avérer crucial. Les techniques de résolution des EDO sont basées sur :  l'approximation géométrique de la fonction  les formules d'intégration numérique (rectangle, trapèze, Simpson...)  les développements de Taylor au voisinage de xi I.5 PROPRIETES DES METHODES NUMERIQUES Plusieurs notions mathématiques sont introduites lors de la résolution d'EDO au moyen de leurs équivalents discrétisés. Les trois principales sont la convergence, la stabilité et la consistance, permettant de relier la solution exacte des équations continues à la solution exacte des équations discrétisées et à la solution numérique obtenue. A ces propriétés, il convient d'ajouter la notion de précision ainsi que des aspects informatiques comme la facilité de mise en oeuvre, les coûts CPU et mémoire. Consistance d'une méthode : La consistance est la propriété qui assure que la solution exacte de l'équation discrétisée tende vers la solution exacte de l'équation continue lorsque le pas de discrétisation h tend vers 0.

ENSET Oran

-5-

CHAPITRE I


Stabilité d'une méthode : C'est la propriété qui assure que la différence entre la solution numérique obtenue et la solution exacte des équations discrétisées reste bornée. Le stabilité indique si l'erreur augmente ou non au cours du calcul. Une méthode peut être stable sous condition (elle sera dite conditionnellement stable) ou toujours stable (elle sera dite inconditionnellement stable). Ordre de précision d'une méthode ~

L'erreur de troncature " est définie comme la différence entre la solution exacte l'approximation numérique obtenue yn, soit :  n

y et



 y ( xn )  yn  (h p ) . L'ordre de

précision de la méthode est donné par l'entier p. Convergence et taux de convergence d'une méthode Une méthode est convergente si, lorsque le pas de discrétisation tend vers 0, la solution numérique tend vers la solution exacte de l'équation continue. Résultat théorique : Une méthode stable et consistante est convergente.

ENSET Oran

-6-

CHAPITRE I


I.6 LES PRINCIPALES METHODES NUMERIQUES : Les principales méthodes de résolution numérique des EDO sont séparées en deux grands types : 

Les méthodes à un pas

Pour ces méthodes, le calcul de la valeur discrète yn+1 au noeud tn+1 fait intervenir la valeur yn obtenue à l'abscisse précédente. Les principales méthodes sont : - Méthodes d'Euler explicite et implicite. - Méthode d'Euler amélioré. - Méthode d'Euler-Cauchy. - Méthode de Crank-Nicholson. - Méthodes de Runge et Kutta. - Runge-Kutta-Fehlberg Method RKF45.



Les méthodes à pas multiples

Pour ces méthodes, le calcul de la valeur discrète yn+1 au noeud tn+1 fait intervenir plusieurs valeurs yn; yn-1; yn-2;…. obtenues aux abscisses précédentes. Les principales méthodes sont : - Méthode de Nystrom ou saute-mouton. - Méthode d'Adams-Bashforth-Moulton. - Méthode de Gear. - Méthode de Hamming. - Méthode Milne-Simpson. - Méthode Linear Shooting. - Méthode Finite-Difference.

ENSET Oran

-7-

Chapitre II Méthode de résolution des ODE a un pas

CHAPITRE II

Méthode de résolution des ODE a un pas

Introduction Le but de ce chapitre est de montrer comment programmer facilement et efficacement les schémas classiques de résolution numérique des équations différentielles à l'aide de Matlab. On met en évidence leurs ordres de convergence respectifs sur un exemple simple. II.1 METHODES A UN PAS  La formulation générale des méthodes à un pas explicite est :

y 0 donné y n 1  y n   ( xn , yn , h) Où la fonction  définit la méthode utilisée.  La formulation générale des méthodes à un pas implicite est :

y 0 donné y n 1  y n   (x n , yn , yn 1 , h) L'obtention de la solution à chaque abscisse nécessite la résolution d'une équation. Ces méthodes sont obtenues en intégrant l'équation différentielle et en utilisant des formules d'intégration numérique pour le second membre. L'ordre du schéma est égal au degré du polynôme pour lequel l'intégration est exacte + 1. Résultats théoriques : 1) Si la fonction  est lipchitzienne par rapport à la deuxième variable alors les méthodes explicites sont stables. 2) Les méthodes implicites sont toujours stables. 3) Les méthodes sont consistantes si x  [a; b];  (x; y; 0) = (x; y). II.1 Méthodes d'Euler explicite et implicite En parlant des solutions numériques aux ODE, chacune commence par la méthode d'Euler, puisqu'elle est facile à comprendre et simple pour programmer. Explication de la méthode d’Euler : On part d’un point M de coordonnées (XM, YM) appartenant à la solution approchée. Le point suivant est le point E de coordonnées (XE, YE), tel que XE = XM +h et tel que le point E est sur la droite passant par M de pente F (XM, YM), On itère le procédé en partant du point (X0, Y0) et en reprenant à chaque étape le point E comme nouveau point M, autant de fois qu’il le faut pour que X décrive l’intervalle demandé. Voir la figure II.1.

ENSET ORAN

-9-

CHAPITRE II


Figure II.1. Méthode d’Euler Quoique sa basse exactitude la garde d'être employée couramment pour résoudre des odes, elle nous donne un indice au concept de base de la solution numérique pour une équation simplement et clairement. Considérons une équation de premier ordre :

y ' (t )  ay (t )  r

avec

y(0)  y0

II.1.1

Sa solution analytique est de la forme suivante :

r r y (t )  ( y0  )e at  a a

II.1.2

Ce qui peut être obtenu par l'utilisation d’une méthode conventionnelle ou la technique de transformé de Laplace. Cependant, une solution analytique exacte n’exister pas pour chaque équation ; même si elle existe, il n'est pas facile de trouver même en utilisant un ordinateur équipé des possibilités du calcul symbolique. C'est pourquoi nous devrions étudier les solutions numériques aux équations différentielles. Comment traduisons-nous l'équation en forme qui peut facilement être manipulée par ordinateur? Tout d'abord, nous devons remplacer la dérivée y’(t)= dy/dt dans l'équation par un dérivé numérique, où le pas « h » est déterminé sur les conditions d’exactitude et les contraintes de calcul.

y (t  h)  y (t )  ay (t )  r h y (t  h)  (1  ah) y (t )  hr

ENSET ORAN

- 10 -

avec y(0)  y 0

(II.1.3)

CHAPITRE II


Et résoudre cette équation différentielle point par point avec l'augmentation de t par h chaque fois, de t = 0.

y (h)  (1  ah) y (0)  hr  (1  ah) y0  hr y (2h)  (1  ah) y (h)  hr  (1  ah) 2 y0  (1  ah)hr  hr

(II.1.4)

y (3h)  (1  ah) y (2h)  hr  (1  ah)3 y0  m0 (1  ah) m hr 2

......................................................................................... C'est un ordre numérique {y (KH)}, que nous appelons une solution numérique de l’équation (II.1.1). Pour être spécifiques, on laisse les paramètres et la valeur initiale d'Eq. (II.1.1) soit a = 1, r = 1, et y0 = 0. Puis, la solution analytique (II.1.2) devient :

y (t )  1  e  at

(II.1.5)

Nous faisons un programme de MATLAB « comparaison1.m », qui emploie la méthode d'Euler pour l'équation (II.1.1), Et la solution numérique (II.1.4) avec le pas h = 0.5 et h = 0.25 sont énumérée dans le tableau II.1 et représentée dans fig. II.2. Tableau II.1 une solution numérique de l'équation (II.1.1) obtenue par la méthode de l'Euler

ENSET ORAN

- 11 -

CHAPITRE II


Les graphiques semblent nous indiquer qu'un petit pas « h » aide réduire l'erreur afin trouver la solution numérique la plus près de la solution analytique (vraie). En fait, un pas « h » trop petite fait non seulement le temps de calcul plus longtemps (proportionnel en tant que 1/h), mais a également comme conséquence des erreurs un peu plus grandes dues à l'effet approximatif accumulé. C'est pourquoi nous devrions rechercher d'autres méthodes pour diminuer les erreurs plutôt que réduisons simplement le pas « h ». La méthode d'Euler peut également être appliquée pour résoudre une équation de vecteur de premier ordre

y ' (t )  f (t , y )

avec y(t 0 )  y0

(II.1.6)

Ce qui est équivalent à une équation scalaire d'ordre élevé. L'algorithme peut être décrit près :

y k 1  yk  hf (t k , y k )

ENSET ORAN

avec y(t 0 )  y0

- 12 -

(II.1.7)

CHAPITRE II


II.1.1 Méthode explicite, d'ordre 1, dont l'algorithme est :

y 0 donné y n 1  y n  h  ( xn , yn )

(II.1.8)

La méthode peut s'interpréter de plusieurs manières : 1) Via les formules d'intégration numérique : la méthode est le résultat de l'application de la formule des rectangles basée au point xn. 2) Géométriquement : la méthode revient à remplacer localement en chaque point xn la courbe solution par sa tangente. 3) Via les développements de Taylor : la méthode provient du développement de Taylor d'ordre 1 de la fonction y au voisinage de xn.

ENSET ORAN

- 13 -

CHAPITRE II


II.1.2 Méthode implicite, d'ordre 1, dont l'algorithme est :

y 0 donné y n 1  y n  hf(x n 1 , y n 1 )

ENSET ORAN

- 14 -

(II.1.9)

CHAPITRE II


II.2 Méthode d'Euler amélioré Méthode explicite dont l'algorithme est :

y0 donné h f ( xn , y n ) 2 h  yn   f ( xn , yn )  f ( xn1 , y * n1 ) 2

y * n1  yn  yn1

(II.2.1)

Géométriquement, la méthode consiste à remplacer dans la méthode d'Euler la pente de la tangente en (xn; yn) par la valeur corrigée au milieu de l'intervalle [xn; xn+1].

II.3 Méthode d'Euler-Cauchy Cette méthode est dite aussi la méthode dite de Newton-Cauchy modifiée ou méthode du point milieu. C’est une méthode explicite dont l'algorithme est :

y0 donné h f ( xn , y n ) 2 h  yn  hf ( xn  , y * n1 ) (II.2.3) 2

y * n1  yn  yn1

ENSET ORAN

- 15 -

CHAPITRE II


La méthode du point milieu consiste à approcher, au voisinage de t0, la solution de l'équation différentielle y'(t) = f (t, y (t)), y(t0) = y0 t0 [a;b], par une parallèle à la tangente au "point milieu de la solution obtenue par la méthode d'Euler" : Le "point milieu de la solution obtenue par la méthode d'Euler" est le point de coordonnées :

(t 0  h/2, y 0  h/2f (t 0 , y 0 )) Et sa tangente a donc comme pente :

  f (t 0  h/2, y 0  h/2f (t 0 , y 0 ))

La méthode du point milieu fait passer du point (t0, y0) au point (t1, y1) avec :

t1  t 0  h y1  y 0  h *   y 0  h * f (t 0  h/2, y 0  h/2 * f (t 0 , y 0 )) qui réalise une seule étape de cette méthode et qui permet de passer d'un point d'abscisse t0 au point d'abscisse t0 + h, ces points étant situés sur le graphe de la solution approchée par cette méthode.

La méthode d’Euler est une méthode numérique peu coûteuse numériquement, mais peu précise quand on intègre sur plusieurs pas de temps. Des améliorations sont possibles des que l’on considère des points intermédiaires, ce que nous allons voir ci-dessous en considérant des méthodes dites de Runge-Kutta.

ENSET ORAN

- 16 -

CHAPITRE II


Comparaison entre Euler explicite Euler modifier et Euler Cauchy : En prenant la même équation précédente on trouve les résultats suivant. Fig (II.3) :

Dans la figure II.3 on remarque que seule la méthode d’Euler Cauchy s’approche de la solution exacte avec un pas de h= 0.2 Et pour la figure II.4 qui se trouve ci-dessous on voit bien que l’erreur est petite pour la méthode d’Euler modifier par contre la méthode d’Euler explicite prend une très grande erreur. erreur_Euler_explecite = 0.0192 erreur_Euleur_Modifier = 0.0091 erreur_euler_cauchy = 0.0013 c’es la plus petite valeur si on compare les trois erreur relative

ENSET ORAN

- 17 -

CHAPITRE II


II.4 Méthode de Crank - Nicholson Méthode implicite, d'ordre 2, dont l'algorithme est :

y0 donné y n1  yn 

h  f ( xn , yn )  f ( xn1 , yn1 ) (II.4.1) 2

Elle est obtenue en utilisant la formule d'intégration numérique des trapèzes.

ENSET ORAN

- 18 -

CHAPITRE II


II.5 Méthode de Heun Dans cette méthode (encore appelée méthode de Newton modifiée), on évalue l'incrément en deux étapes. La première étape consiste à écrire : xn 1

yn1  yn   f ( yn , xn )dx (II.5.1) xn

en estimant l'intégrale par la méthode des trapèzes, soit :

h yn1  yn  ( f ( yn , xn )  f ( yn1 , xn1 )) (II.5.2) 2

(A)

Cette formule est une formule implicite en yn+1 et elle n'est généralement pas inversible. De ce fait, on ne connaît pas encore yn+1 à ce stade. On trouve une valeur approchée de f (yn+1; xn+1) grâce à la méthode d'Euler. C'est la seconde étape. On a

pn1  yn  hf ( yn , xn )  yn1

(II.5.3)

La méthode de Heun met en jeu ce que l'on appelle classiquement un schéma de type prédicateur/correcteur. Le schéma prédicateur donne une estimation de yn+1 qui permet de prédire la valeur de la fonction en xn+1. Le schéma correcteur associé est donné par la relation (II.5.2). Globalement, on a donc :

yn1  yn 

ENSET ORAN

h  f ( yn , xn )  f ( yn  hf ( yn , xn ), xn1 ) 2

- 19 -

(II.5.4)

CHAPITRE II


II.6.Méthode de Runge-Kutta Les méthodes de Runge-Kutta sont couramment utilisées car elles allient précision, stabilité et simplicité. Toutes nécessitent plusieurs itérations pour effectuer un pas. II.6.1 Méthode de Rungue Kutta 2 La méthode de Runge-Kutta d’ordre 2 est la plus simple ; elle combine deux itérations successives de la méthode d’Euler explicite. Dans un premier temps la dérivée en (xn , yn ) est évaluée pour faire une première estimation du point suivant (noté « A » dans la figure ci II.5 dessous). L’estimation provisoire de yn+1 en ce point est ensuite utilisée pour affiner le calcul de la dérivée. Une nouvelle approximation de celle-ci est obtenue en prenant sa valeur à mi-parcours (prise au point B). C’est cette valeur de la dérivée qui sera ensuite utilisée pour estimer le prochain pas yn+1 (point C). La procédure d’intégration se résume à :

  hf ( xn , yn ) h    hf ( xn  , y n  ) 2 2 y n1  yn     (h 3 )

(II.6.1)

Même si cette méthode demande deux fois plus opérations de calcul que la méthode d’Euler explicite pour effectuer un seul pas, le résultat est plus précis et plus stable.

ENSET ORAN

- 20 -

CHAPITRE II


Figure II.5 : Principe de fonctionnement de la méthode de Runge-Kutta du second ordre pour effectuer un seul pas Une méthode explicite, d'ordre 2, peut être obtenue par l'utilisation de la formule des trapèzes. L'algorithme, noté RK2, s'écrit :

y0 donné y * n1  yn  hf ( xn , yn )

(II.6.2)

h yn1  yn  [ f ( xn , yn )  f ( xn1 , y * n1 )] 2 Elle est obtenue en utilisant la formule d'intégration numérique des trapèzes. II.6.2 Méthode de Runge et Kutta explicite a un point

y n ,1  y n 

h f ( xn , y n ) 2

yn1  yn  h(1   ) f ( xn , yn )  f ( xn 

h , y n ,1 ) 2

(II.6.3)

y0 : donné Avec α un nombre réel compris entre 0 et 1. Les valeurs de α couramment utilisées sont α=1, α=1/2 et α=3/4. Ces méthodes sont d’ordre 2.

ENSET ORAN

- 21 -

CHAPITRE II


II.6.3 Méthodes RK implicites à un point La formule de récurrence est définie par la relation :

yn1  yk  h[(1   ) f ( xn , yn )  f ( xn1 , yn1 )]

(II.6.4)

Oû θ est un nombre réel appartenant à l’intervalle] 0, 1]   

Si θ = 0, on retrouve la méthode d’Euler. Si θ = 1/2, la méthode est d’ordre 2 et s’appelle la méthode des trapèzes. Si θ ≠ 1/2, la méthode est d’ordre 1.

II.6.4 Méthodes RK explicites à 2 points intermédiaires Ces méthodes sont définies par les relations :

h f ( xn , y n ) 3 2h h yn , 2  yn  f ( xn  , yn ,1 ) 3 3 h 2h yn1  yn  ( f ( xn , yn )  3 f ( xn  , yn , 2 )) 4 3 y n ,1  y n 

Ou par

y n ,1  y n 

(II.6.5)

h f ( xn , y n ) 2

h yn , 2  yn  h( f ( xn , yn )  2 f ( xn  , yn ,1 )) 2 h h yn1  yn  ( f ( xn , yn )  4 f ( xn  , yn ,1 )  f ( xn1 , yn , 2 )) 6 2

(II.6.6)

Ces deux méthodes sont d’ordre 3. La première est parfois appelée méthode de Heun. II.6.5 Méthodes RK explicites à 3 points intermédiaires Une méthode explicite, d'ordre 4, peut être obtenue par l'utilisation de la formule de Simpson. La méthode avec un pas de h : A partir du point M, on construit les point E, RK1, RK2. « RK3 » est le milieu du segment (M, RK2). Le point suivant de la solution approchée est le point RK4, de coordonnées (XRK4, YRK4), avec : XRK2 = XE = XM + h, et qui, de plus, appartient à la droite passant par M et de pente :

ENSET ORAN

- 22 -

CHAPITRE II


1/6(F (XM, YM) + 2 F (XRK1, YRK1) + 2 F (XRK3, YRK3) + F (XRK2, YRK2)) Voir la figure ci-dessous. (II.6)

Figure.II.6 Méthode de Runge-Kutta d’ordre 4 On itère le procédé en partant du point (X0, Y0) et en reprenant à chaque étape le point E comme nouveau point M, autant de fois qu’il le faut pour que X décrive l’intervalle demandé. La méthode suivante est de loin la plus connue et utilisée. Les relations de récurrence sont les suivantes :

h f ( xn , y n ) 2 h h yn , 2  yn  f ( xn  , yn ,1 ) (II.6.7) 2 2 h h y n , 3  y n  f ( xn  , y n , 2 ) 2 2 h h h yn1  yn  ( f ( xn , yn )  2 f ( xn  , yn ,1 )  2 f ( xn  , yn , 2 )  f ( xn1 , yn , 3 )) 6 2 2 y n ,1  y n 

ENSET ORAN

- 23 -

CHAPITRE II


Méthode est d’ordre 4. La méthode de Runge-Kutta d’ordre 4 combine quatre itérations successives. Sa précision est généralement encore meilleure. La procédure d’intégration devient

y0 donné

k1  hf  xn , yn  h k   k 2  hf  xn  , yn  1  2 2  h k   k3  hf  xn  , yn  2  2 2  k 4  hf  xn  h, yn  k3  yn1  yn 

(II.6.8)

1 k1  2k 2  2k3  k4  6

II.6.6 Formule générale des méthodes RK explicites : Les méthodes de Runge-Kutta s’écrivent de manière générale

K1  h[ f ( xn  1h, yn   1,1 K1   1, 2 K 2  ...   1,n K n )] K 2  h[ f ( xn   2 h, yn   2 ,1 K1   2 , 2 K 2  ...   2 ,n K n )] ...... K n  h[ f ( xn   n h, yn   n ,1 K1   n , 2 K 2  ...   n ,n K n )]

(II.6.9)

yn1  yn  h[ 1 K1   2 K 2  ... n K n ] Les coefficients sont détermines afin que l’ordre soit le plus élevé possible. On note :  A la matrice de coefficients (αi,j),  Γ le vecteur des coefficients γi.  Θ le vecteur des coefficients θi. Quand la matrice A est triangulaire inférieure stricte, αi,j = 0 pour j ≥ i, on dit que la méthode est explicite. Si seule la partie triangulaire supérieure est nulle, αi,j = 0 pour j > i, la méthode est dite implicite ; sinon elle est totalement implicite. ENSET ORAN

- 24 -

CHAPITRE II

ENSET ORAN


- 25 -

CHAPITRE II


II.7. Méthode de Burlish - Stoer L'idée de cette méthode repose sur les trois principes suivants : la résolution de l'équation différentielle pour un accroissement de x est donnée par une fonction qui dépend de h, mais qui tend vers une limite finie (indépendante de h) ; on cherche à estimer la valeur exacte du système différentiel à intégrer en calculant pour différents pas et en prenant la limite d'un pas tendant vers zéro. La figure II.7 illustre graphiquement les différentes estimations de la valeur de la fonction solution en utilisant trois pas d'intégration différents. La limite d'un pas de temps nul correspond à la solution exacte donnée par la courbe en trait plein. Ce point est très comparable à celui de la méthode de Romberg décrite pour l'intégration. Le second principe de la méthode consiste à extrapoler cette limite non pas sur la base de développements polynomiaux mais de développements en fractions rationnelles. Le troisième principe consiste à utiliser des fonctions d'erreur qui sont paires en pas d’intégration.

Fig. II.7 {Schéma illustrant le principe de la méthode de Burlish-Stoer. Les lignes en pointille correspondent aux intégrations réalisées avec trois pas d'intégration différents. La courbe en trait plein correspond à la solution exacte.

ENSET ORAN

- 26 -

CHAPITRE II


II.8 Comparaison entre quelques méthodes Sur le programme suivant on va comparer l’allure des courbes pour quelques méthodes, que nous avons étudiées précédemment et voir la rapidité de chaque méthode Avec la fonction y’(t)= -y+1 qui a pour solution exacte : y(t) = 1-exp(-a*t) et y(0) = 0

Voici le temps de résolution pour chaque méthode : time_Euler = 0.6880 time_Heun = 1.4380 time_RK4 = 2.7660

ENSET ORAN

- 27 -

CHAPITRE II


On voit bien que la méthode d’Euler est la plus rapide parce que dans les autres méthodes y a plus de calcul pour approché la solution exacte. Pour l’allure des courbes on obtient la figure II.8 :

Figure II.8 Solution numérique d’une équation du 1er Ordre

ENSET ORAN

- 28 -

CHAPITRE II


II.9 Comparaison entre toutes les méthodes Si on prend la comparaison de toutes les méthodes vue dans ce chapitre on obtient la figure II.8 suivi du programme sous Matlab ci-dessous enregistré sous le nom de « comparaison » Programme

ENSET ORAN

- 29 -

CHAPITRE II


Suite du programme

Sur la figure II.9 on voit bien que Euler Explicite s’éloigne de la vraie solution par contre les plus précise sont RK4 et Heun donc c’est les méthodes les plus précises surtout pour RK4 erreur_RK4 = 1.1979e-007 erreur la plus petite erreur_Euler_explecite = 0.0070 erreur_Euler_Implecite = 0.0077 erreur_Euleur_Modifier = 0.0034 erreur_Heun = 2.3774e-004 erreur_euler_cauchy = 2.3774e-004 ENSET ORAN

- 30 -

CHAPITRE II


Figure II.9 : Solution graphique pour les différentes ODE Si on compare les résultats obtenu par la figure II.10 ainsi que les résultats de calcul affiché sur le Works-pace du Matlab on voit bien que la méthode de Runge – Kutta a une erreur très faible et que la solution s’approche du résultas exacte par contre pour les autres implicites on remarque qu’il y a une très grande erreurs La figure suivant II.10 nous montre les erreurs de chaque méthode par rapport à la solution exacte

ENSET ORAN

- 31 -

CHAPITRE II


FigureII.10 : Allure des différentes méthodes de résolution des ODE II.10. La méthode de Runge-Kutta-Fehlberg Une façon pour garantir l'exactitude dans la solution d'un système d’équation différentielle est de résoudre le problème en employant deux fois les tailles du « pas » h et h/2 et de comparer les réponses aux points de maille correspondant à la taille du « pas » le plus grand. Mais ceci exige une quantité significative de calcul pour la taille du pas la plus petite et doit être répété si on le détermine que l'accord n'est pas assez bon. La méthode de Runge-Kutta-Fehlberg (RKF45 dénoté) est à sens unique pour essayer de résoudre ce problème. Elle a un procédé pour déterminer si la taille du « pas » appropriée h est employée. À chaque étape, deux approximations différentes pour la solution sont faites et comparées : Si les deux réponses sont en accord étroit, l'approximation est acceptée. Si les deux réponses ne sont pas d'accord sur une exactitude spécifique, la taille du « pas » est réduite. Si les réponses sont d'accord sur des chiffres plus significatifs qu'exigés, la taille du pas est augmentée. ENSET ORAN

- 32 -

CHAPITRE II


Chaque étape exige l'utilisation des six valeurs suivantes :

k1  hf ( xn , yn ), 1 1 k 2  hf ( xn  h, yn  k1), 4 4 3 3 9 k 3  hf ( xn  h, yn  k1  k 2 ), 8 32 32 k 4  hf ( xn 

12 1932 7200 7296 h, y n  k1  k2  k 3), 13 2197 2197 2197

k 5  hf ( xn  h, yn 

439 3680 845 k1  8k 2  k3  k 4), 216 513 4104 (A)

1 8 3544 1859 11 k 6  hf ( xn  h, yn  k1  2k 2  k3  k 4  k 5). 2 27 2565 4104 40 Alors une approximation à la solution du système d’équation différentielle est faite pour suivre une méthode de Runge-Kutta d'ordre 4 :

yn 1  yn 

25 1408 2197 1 k1  k3  k 4  k5 216 2565 4101 5

(B)

Là où les quatre valeurs de fonction f1, f3, f4, et f5 sont employées. Notez que le f2 n'est pas employé dans la formule (B). Une meilleure valeur pour la solution est déterminée suivre une méthode de Runge-Kutta d'ordre 5 :

z n 1  yn 

ENSET ORAN

16 6656 28,561 9 2 k1  k3  k 4  k 5  k 6. 135 12,825 56,430 50 55

- 33 -

(C)

CHAPITRE II


La taille du pas optimale SH peut être déterminée en multipliant les temps scalaires de s avec la taille du pas courant H. Le s scalaire est :

 tolh s    2 z n 1  yn 1

14

  tolh   0.84  z y n 1   n 1

   

14

(D)

Là où le tol est la tolérance spécifique de commande d'erreur. Exemple 1 : résolution d’une équation avec la méthode RKF45 et RK4

y '  1  y 2 avec y (0)  0 [0, 1.4] Tableau II.2 : RKF45

La dérivation de la formule (D) peut être trouvée dans les livres avancés sur l'analyse numérique. Il est important d'apprendre qu'une taille fixe du « pas » n'est pas la meilleure stratégie quoiqu'elle donne une table gentille apparaissante des valeurs. Si les valeurs nécessaires ne sont pas dans la table, l'interpolation polynôme devrait être employée.

ENSET ORAN

- 34 -

CHAPITRE II


Tableau II.3 : RK4

Comparaison entre RKF45 et RK4 pour la résolution d’une équation différentielle. Un programme RKF45 a été employé avec le tol de valeur = 2 × 10-5 pour la tolérance de commande d'erreur. Il a changé la taille du pas automatiquement et a produit les 10 approximations à la solution dans le tableau II.2. Un programme RK4 a été employé avec la taille du pas a priori de h = 0.1, ce qui a exigé de l'ordinateur de produire de 14 approximations aux points équidistants dans le tableau .II.3. Les approximations au bon point final sont :

y (1.4)  y10  5.7985045

et

y(1.4)  y14  5.7919748

et

E14  0.0059089

Et l’erreur est de :

E10  0.0006208

Pour les méthodes RKF45 et RK4, respectivement. La méthode RKF45 a une erreur très petite

ENSET ORAN

- 35 -

CHAPITRE II

ENSET ORAN


- 36 -

CHAPITRE II

ENSET ORAN


- 37 -

CHAPITRE II


Conclusion Les intégrateurs implicites sont souvent moins précis, plus compliqués à utiliser. Mais ils sont plus stables que les intégrateurs explicites : ils divergent moins rapidement. En d’autres termes : moins précis à court terme, plus précis à long termes. Pour les méthodes RKF45 et RK4, respectivement. La méthode RKF45 a une erreur très petite Donc on peut dire que : •

Runge-Kutta d’ordre 4 :  Souvent la réponse par défaut  Précision est généralement encore meilleure  Bon rapport qualité/prix  Mais pas une réponse universelle !

•

Euler :  Beaucoup de défauts  Déconseillée pour presque tout  Mais tellement rapide à implémenter

•

Heun :  Fait la moyenne entre 2 pentes ; les méthodes de Runge-Kutta font la moyenne de m pentes …

La précision dépend de l’ordre  Euler : ordre 1  Heun : ordre 2  Runge-Kutta : ordre 4  Cependant utilisation d’ordre supérieur n’est pas très intéressant …

ENSET ORAN

- 38 -

Chapitre III Méthode de résolution des ODE à pas multiples

CHAPITRE III

Méthode de résolution des ODE à pas multiples

III Méthode à pas multiples Introduction Un intégrateur “intelligent” possède une procédure de contrôle de la méthode de convergence, c’est à dire un moyen d’estimer l’erreur commise par le calcul sur un pas díntégration et la possibilité de choisir en conséquence un nouveau pas si le système différentiel aborde une région ou la fonction prend des valeurs plus importantes. Le calcul de cette estimation entraîne un surcoût de calcul, qu’il convient de bien gérer afin de minimiser cet effort supplémentaire. L’idée la plus simple pour estimer cette erreur consiste à calculer la solution donnée par un algorithme (de Runge-Kutta d’ordre 4 par exemple) pour deux pas díntégration différents, h et 2h. Soit y(x + 2h) la solution exacte à x + 2h et y(x + h) la solution exacte à x + h, on a  

y(x + 2h) = y1 + (2h)5+ O(h6) y(x + 2h) = y2 + 2(h5) + O(h6)

Où  est une fonction qui reste constante sur l’intervalle x, x + 2h à l’ordre h5. La première équation correspond à une intégration avec un pas égal à 2h tandis que la seconde correspond à deux intégrations successives avec un pas de h. La différence  = y2 − y1 fournit une estimation de l’erreur commise avec un pas d’intégration h. III-1 Méthode Adams–Bashforth–Moulton Method La méthode d'Adams-Bashforth-Moulton (ABM) se compose de deux étapes. La première étape est de rapprocher f (t, y) par le polynôme (de Lagrange) de degré 4 assortissant les quatre points :

t k 3 , f k 3 , t k 2 , f n2 , t k 1 , f n1 , t , f  Méthodes basées sur des techniques d'intégration numérique, dont la formulation générale est : p

yk 1  yk  h   j f (t k  j , y k  j ) j  1

et substituez le polynôme dans la forme intégrale d'équation pour obtenir une évaluation prévue de yk+1. h

pk 1  yk   l3 (t )dt  yk  0

h (9 f k 3  37 f k 2  59 f k 1  55 f k ) 24

La deuxième étape est de répéter le même travail avec les quatre points mis à jour :

ENSET ORAN

- 40 -

(III.1.1a)

CHAPITRE III


t k 2 , f k 2 , t k 1 , f k 1 , t k , f k , t k 1 , f k 1 

(f k 1  f (t k 1 , pk 1 ))

pour obtenir une évaluation corrigée de yk+1. h

ck 1  yk   l '3 (t )dt  yk  0

h ( f k 3  f n1  19 f k  9 f k 1 ) 24

(III.1.1b)

Alternativement, nous écrivons à Taylor l'expansion de série de yk+1 au sujet de tk et de celui du yk au sujet de tk+1 avec

et remplacez les premiers, deuxièmes, et troisième dérivés par leurs approximations de différence.

ENSET ORAN

- 41 -

CHAPITRE III


De ces équations et dans la prétention cela f k 1 facteur prédictif/correcteur avec : (4)

 f k(4)  K , nous pouvons écrire les erreurs de

251 5 (4) 251 5 h fk  Kh 720 720 19 5 (4) 19 EC ,k 1  yk 1  p k 1  h f k 1  Kh 5 720 720 EP ,k 1  yk 1  p k 1 

(III.1.4a) (III.1.4b)

Nous ne pouvons pas encore employer ces formules pour estimer les erreurs de facteur prédictif/correcteur, puisque K est inconnu. Mais, de la différence entre ces deux formules

EP ,k 1  EC ,k 1  ck 1  p k 1 

270 5 270 270 Kh  E P,k 1  E C,k 1 720 251 19

Nous pouvons obtenir les formules pratiques pour estimer les erreurs avec

251 (ck 1  pk 1 ) 270 19 EC ,k 1  yk 1 c k 1   (ck 1  pk 1 ) 270 EP ,k 1  yk 1  p k 1 

ENSET ORAN

- 42 -

(III.1.6a) (III.1.6b)

(III.1.5)

CHAPITRE III


Ces formules nous donnent de grandes évaluations de la façon dont valeurs prévues/corrigées étroites sont à la valeur vraie et ainsi peuvent être employées pour les améliorer aussi bien que pour ajuster l'étape-taille.

251 (ck  pk )  mk 1 270 19  ck 1  (ck 1  pk 1 )  y k 1 270

pk 1  pk 1  c k 1

(III.1.7a) (III.1.7b)

On s'attend à ce que ces formules de modification récompensent nos efforts que nous avons fait pour les dériver. La méthode d'Adams-Bashforth-Moulton (ABM) avec les formules de modification peut être décrite par Eqs. (III.1.1a), (III.1.1b), et (III.1.7a), (III.1.7b) récapitulé ci-dessous et est moulé dans ode_ABM courant de MATLAB le « () ». On s'attend à ce qu'ainsi cet arrangement a besoin de seulement deux évaluations de fonction (appels) par itération, tout en ayant une erreur de troncation d'O (h5) et fonctionne mieux que les méthodes discutées jusqu'ici. Il est mis en application par la routine intégrée de MATLAB « ode113 () » avec beaucoup de techniques sophistiquées additionnelles.

h (9 f n3  37 f n2  59 f n1  55 f n ) 24 251 Modificateur : m n 1  p n 1  (cn  p n ) 270 h Correcteur : c n 1  yn  ( f n2  5 f n1  19 f n  9 f ( xn1 , mn1 )) 24 19 y n 1  cn1  (cn 1  pn1 ) 270 Prédicteur :

ENSET ORAN

p n 1  yn 

- 43 -

(III.1.8a) (III.1.8b) (III.1.8c) (III.1.8d)

CHAPITRE III


Programme de ABM

III.2.Méthode de Hamming Dans cette section, nous présentons juste l'algorithme de la méthode de Hamming [H-1] récapitulée dans la boîte ci-dessus et le « ode_Ham courant correspondant () », qui est une autre méthode multi pas de pré enregistreur/correcteur comme la méthode d'Adams-Bashforth Moulton (ABM). Cet arrangement a besoin également de seulement deux évaluations de fonction (appels) par itération, tout en ayant l'erreur d'O (h5) et est ainsi comparable à la méthode d'ABM discutée dans la section précédente.

ENSET ORAN

- 44 -

CHAPITRE III


Méthode de Hamming avec des formules de modification

4h (2 f k 2  f k 1  2 f k ) (III.1.9a) 3 112 Modificateur : m k 1  p k 1  (ck  p k ) (III.1.9b) 121 1 Correcteur : c k 1  9 yk  y k 2  3h( f k 1  2 f k  f ( xk 1 , mk 1 )) 8 9 y k 1  ck 1  (ck 1  pk 1 ) (III.1.9d) 121 Prédicteur :

p k 1  yk 3 

Programme de Hamming :

ENSET ORAN

- 45 -

(III.1.9c)

CHAPITRE III


III.3.Comparaison des méthodes (RK4 ; Hamming, ABM) Les facteurs principaux à considérer dans l'évaluation/comparant différentes méthodes numériques sont l'exactitude de la solution numérique et son temps de calcul. Dans cette section, nous comparerons les routines : “ode_RK4()”, “ode_ABM()”, “ode_Ham()”, “ode23()”, “ode45()”, and “ode113()” . Il est important de noter que l'évaluation/comparaison des méthodes numériques n'est pas aussi simple parce que leurs exécutions peuvent dépendre de la caractéristique du problème actuel. Il devrait également noter qu'il y a autre facteurs à considérer, comme la stabilité, polyvalence, preuve contre l'erreur d'exécution, et ainsi de suite. Ces points sont considérés dans la plupart des routines de fonction intégrée de MATLAB. La première chose que nous allons faire est de valider l'efficacité des modificateurs (Eqs. (II.1.8b, d) et (III.1.9b, d)) dans l'ABM (Adams-Bashforth-Moulton) méthode et la méthode de Hamming et RK4. Pour ce travail, nous écrivons et lançons le programme « nm643_1.m »

ENSET ORAN

- 46 -

CHAPITRE III


Pour obtenir les résultats représentés dans fig. III.1 pour l'équation

y' (t)  - y(t)  1 avec y(0)  0

(6.4.10)

Solution exacte

y (t )  1  e  t

La Fig. III.1 nous montre un fait intéressant ce, bien que la méthode d'ABM et la méthode de Hamming, même sans modificateur, sont théoriquement prévus avoir une meilleure exactitude que la méthode RK4, ils s'avèrent fonctionner mieux que RK4 seulement avec modificateurs. On remarque que l’erreur relative est plus grande pour la méthode de RK4 par rapport aux autres méthodes. Tableau III-1 : Résultats d'appliquer plusieurs routines pour résoudre une équation simple err =1.0e004 * time =

ode_RK4 0.0061

ode_ABM 0.0050

ode_Ham 0.0049

ode23 0.1275

ode45 0.0062

ode113 0.0294

0.3280

0.2030

0.1720

0.4840

0.1880

0.3280

Figure III.1 Solutions numériques et leurs erreurs pour l'équation y’(t) = -y(t)+1

ENSET ORAN

- 47 -

CHAPITRE III


Naturellement, ce n'est pas toujours le cas, comme illustré dans fig. III.2, ce que nous avons obtenu en appliquant les mêmes routines pour résoudre une autre équation : y' (t)  y(t)  1 avec y(0)  0

(6.4.11)

Là où la véritable solution analytique est : y(t)  e t - 1

(6.4.12)

Il est remarquable que, sans modificateurs, la méthode d'ABM semble être meilleure que Méthode de Hamming ; Cependant, avec les modificateurs, c'est l'autre manière autour du mineurs ils courent un cou et une course de cou. N'importe qui verra que les méthodes de pré enregistreur/correcteur telles que la méthode d'ABM (l'ode_ABM ()) et la méthode de Hamming (ode_Ham ()) nous donne une meilleure solution numérique avec moins d'erreur et de temps plus court de calcul que les routines intégrées de MATLAB « ode23 () », « ode45 () », et « ode113 () » aussi bien que la méthode RK4 (ode_RK4 ()), comme énuméré dans le tableau III.2. Mais, une conclusion générale ne devrait pas être déduite juste d'un exemple. Tableau III-2 : Résultats d'appliquer plusieurs routines pour résoudre une équation simple

err =1.0e-004 *

ode_RK4 0.0925

ode_ABM 0.0203

ode_Ham 0.0179

ode23 0.4770

ode45 0.0422

ode113 0.1249

time =

0.05

0.03

0.03

0.07

0.05

0.05

ENSET ORAN

- 48 -

CHAPITRE III


Figure III.2 Solutions numériques et leurs erreurs pour l'équation y’(t) = y(t)+1

ENSET ORAN

- 49 -

CHAPITRE III

ENSET ORAN


- 50 -

CHAPITRE III


III.4 Méthode de Milne Simpson :

Prédicteur : Correcteur :

ENSET ORAN

4h (2 f k 2  f k 1  2 f k ) 3 h y k 1  yk 1  ( f k 1  4 f1  f k 1 ) 3

p k1  yk 3 

- 51 -

CHAPITRE III


III.5. Méthodes de Gear Les méthodes de Gear ne sont pas construites à partir de techniques d'intégration numérique mais directement à partir de polynôme d'interpolation passant par p points (xi+1; yi+1), (xi; yi),...,(xi-p+1; yi+1). III.5.1 Méthode de Gear à 2 pas, implicite, d'ordre 2 : y0 donné y1 calculé avec une méthode à un pas yn1 

4 1 2h y n  yn1  f ( xn1 , yn1 ) 3 3 3

III.5.2 Méthode de Gear à 3 pas, implicite, d'ordre 3 : y0 donné y1 , y 2 calculé avec une méthode à un pas yn1 

18 9 2 2h yn  yn1  yn2  f ( xn1 , y n1 ) 11 11 11 11

III.5.3 Méthode de Gear à 4 pas, implicite, d'ordre 4 : y0 donné y1 , y 2 , y 3 calculé avec une méthode à un pas yn1 

48 36 16 3 12h yn  y n1  y n 2  y n 3  f ( xn1 , y n1 ) 25 25 25 25 25

III.6 problème de la valeur marginale Un problème de la valeur marginale (BVP) est une équation de Nième-ordre avec certaines des valeurs de la variable dépendante X (t) et son dérivé spécifiques au temps initial t0 et d'autres spécifiques au temps final tf. [BVP]N : x (N) (t)  f (t, x(t), x' (t), x (2) (t), . . . , x (N -1) (t)) Avec les valeurs limites x(t1 )  x 10 , x(t 2 )  x 21 , . . . , x (N-1) (t N )  x N, N-1

ENSET ORAN

- 52 -

(6.6.1)

CHAPITRE III


Dans certains cas, quelques relations entre les valeurs initiales et les valeurs finales peuvent être données comme état de frontière mélangé au lieu de valeurs initiales/finales spécifiques. Cette section couvre la méthode de tir et la méthode de différence finie cela peut être employé pour résoudre un BVP de second ordre :

[BVP]2 : x " (t)  f (t, x(t)

avec x(t 0 )  x 0 (t f )  x f

(III.6.2)

III.6.1 Méthode de Tir L'idée de cette méthode est d'assumer la valeur de x’ (t0), résolvez alors le différentiel d’équation (IVP) avec l'état initial [x (t0) x’ (t0)] et subsistance ajustant la valeur de x’ (t0) et de résoudre l'IVP répétitivement jusqu'à la valeur finale X (tf) de la solution assortit, le xf donné de valeur avec assez d'exactitude. Il est semblable à l'ajustement de l'angle de mettre le feu à un canon de sorte que la coquille frappe par la suite la cible et à cela est pourquoi cette méthode est appelée la méthode de tir. Ceci peut être regardé comme problème non linéaire d'équation, si nous considérons x’ (t0) comme variable indépendante et la différence entre la valeur finale en résultant X (tf) et un xf désiré comme (non-adaptation) fonction de x’ (t0). Ainsi l'arrangement de solution peut être systématisé en employant la méthode sécante et est présenté dans la routine « bvp2_shoot de MATLAB () ». Par exemple, considèrent un BVP comprenant l'équation de second ordre

x" (t )  2x 2 (t)  4tx(t)x' (t)

ENSET ORAN

1 1 avec x(0)  , x(1)  4 3

- 53 -

(III.6.3)

CHAPITRE III


Programme de la méthode du TIR

La solution x(t) et ses dérivés x’(t) sont comme suite :

x(t) 

1 4 - t2

et

x' (t) 

2t  2tx 2 (t) 2 (4 - t )

(III.6.4)

Notez que cette équation du second ordre peut être écrite sous forme d'équation d'état est :

  x'1 (t )   x 2 (t )   x' (t ) 2 x 2 (t )  4tx (t ) x (t )  2   1  1 2

ENSET ORAN

avec

- 54 -

 x 1 ( 0)   x 0  1 / 4   x ( 0)    x  1 / 3   2   f 

(III.6.5)

CHAPITRE III


Afin d'appliquer la méthode de tir, nous avons placé la conjecture initiale de x2(0)=x’ (0) avec :

dx0[1]  x 2 (0) 

x f  x0 t f  t0

(III.6.6)

et résolvez l'équation d'état avec l'état initial [x1(0) x2(0)=dx0[1]]. Puis, selon le signe de la différence « e (1) » entre la valeur finale « x1 (1) » de la solution et du xf final de valeur de cible, nous rendons la prochaine conjecture « dx0 [2] » plus grande/plus petite que la conjecture initiale dx0 [1] et résolvons l'équation d'état encore avec l'état initial [x1 (0) dx0 [2]]. Nous pouvons commencer vers le haut la méthode sécante avec les deux valeurs initiales dx0 [1] et dx0 [2] et répétez l'itération jusqu'à la différence (l'erreur) « e(k) » devient suffisamment petit. Pour ce travail, nous composons le programme de MATLAB “extir.m”, ce qui emploie la routine “bvp2_shoot ()” pour obtenir la solution numérique et la compare à la véritable solution analytique.

Programme de la fonction

La figure III.3 prouve que la solution numérique obtient plus près de la véritable solution analytique après chaque rond de l'ajustement. Avec une erreur estimée a err = 3.5166e-012

ENSET ORAN

- 55 -

CHAPITRE III


Figure III.3 La solution d'un BVP obtenu en employant la méthode de tir. III.7 Méthode des différences finis (Finite Difference) L'idée de cette méthode est de diviser l'intervalle entier [t0, tf] en segments de N de largeur h= (tf-t0)N, et rapprochez les premiers et deuxièmes dérivés dans l’équations pour chaque point de grille par les formules de différence centrale. Ceci mène à un système de tri diagonal des équations en ce qui concerne (n-1) variables [xi= x(t0+ih) = 1…………..N-1]. Cependant, afin que ce système d’équations soit résolu facilement, il devrait être linéaire, en impliquant que ses coefficients peuvent ne pas contenir aucune limite de x Par exemple, on considère un BVP comprenant l'équation linéaire de second ordre

x" (t )  a1 (t ) x' (t )  a0 (t ) x(t )  u (t ) avec x(t 0 )  x0 , x(t f )  xf

(AB)

Selon la méthode de différence finie, nous divisons l'intervalle de solution [t0 tf] dans N, dans des segments de N et convertir l'équation pour chaque point de grille ti = t0 + ih, dans une équation à différences

ENSET ORAN

- 56 -

CHAPITRE III


xi 1  2 xi  xi 1 xi 1  xi 1  a  a0 i xi  ui 1i h2 2h ( 2  ha1i ) xi 1  (4  2h 2 a0 i ) xi  (2  ha1i ) xi 1  2h 2ui Puis, tenant compte de l'état de frontière cela x0 = x(t0) et xN = x(tf), nous rassemblons le tout les (N-1) équations pour construire un système tri diagonal avec des équations

Ceci peut être résolu efficacement en employant « trid » courant de MATLAB le « () », qui est consacré à un système de tri diagonal des équations linéaires. Le procédé entier de la méthode de différence finie pour résoudre une équation linéaire de second ordre avec des états de frontière est moulé dans la routine « bvp2_fdf de MATLAB () ». Cette routine est conçue pour accepter les deux coefficients a1 et a0 et l'entrée u de côté droit de l’équation AB, en tant que ses trois premiers entrés arguments, où n'importe lequel de ces trois arguments entrés peut être indiqué comme nom de fonction au cas où la limite correspondante ne serait pas une valeur numérique, mais une fonction du temps t. Nous faisons le programme « dofdf » pour employer cette routine pour résoudre le BVP de second ordre.

2 2 x" (t )  x' (t )  2 x(t )  0 avec x(1)  5, x(2)  3 t t

ENSET ORAN

- 57 -

(BB)

CHAPITRE III


Programme

Pour obtenir le résultat représenté dans fig. III.4 et, en plus, pour employer la commande symbolique de calcul « dsolve () » et « subs () » d'obtenir la solution analytique

x (t )  t 

4 t2

et substituez le vecteur de temps dans la solution analytique pour obtenir ses valeurs numériques pour le contrôle.

ENSET ORAN

- 58 -

CHAPITRE III


Figure III.4 solution d'un BVP obtenu en employant la méthode de différence finie.

ENSET ORAN

- 59 -

CHAPITRE III


On note à propos de la méthode de tir et de la méthode de différence finie : 1- Tandis que la méthode de tir s'applique à BVP linéaire/non linéaire, la méthode de différence finie convient à BVP linéaire. Cependant, nous pouvons également appliquer la méthode de différence finie d'une façon itérative pour résoudre BVP non linéaire 2- Les deux méthodes peuvent être modifiées pour résoudre BVP avec des états de frontière mélangés

ENSET ORAN

- 60 -

Bibliographie

Bibliographie

Bibliographie Burden, Richard L., and Fairs, J. Douglas, Numerical Analysis, 7th ed., Brooks/ Cole Thomson, Pacific Grove, CA, 2001. Canale, Raymond, and Chapra, Steven, Numerical Methods for Engineers: with Software and Programming Applications, McGraw-Hill, New York, 2002. Fausett, L. V., Applied Numerical Analysis Using MATLAB, Prentice-Hall, Upper Saddle River, NJ, 1999. Hamming, R. W., Numerical Methods for Scientists and Engineers, 2nd ed., McGraw-Hill, New York, 1973. Kreyszig, Erwin, Advanced Engineering Mathematics, 8th ed., John Wiley & Sons, New York, 1999. Lindfield, G. R., and Peny, J. E. T., Numerical Methods Using MATLAB, 8th ed., Prentice-Hall, Upper Saddle River, NJ, 2000. Luenberger, D. G., Linear and Nonlinear Programming, 2nd ed., Addison-Wesley Publishing Company, Reading, MA, 1984. Mathews, J. H., and Fink, K. D., Numerical Methods Using MATLAB, PrenticeHall, Upper Saddle River, NJ, 1999. Maron, Melvin J., Numerical Analysis, Macmillan, Inc., New York, 1982. Oppenheim, Alan V., and Schafer, Ronald W., Discrete-Time Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1989. Peaceman, D. W., and Rachford, H. H., The numerical solution of parabolic and elliptic differential equations, J. Soc. Ind. Appl. Math. 3, 28–41 (1955). Schilling, R. J., and Harris, S. L., Applied Numerical Methods for Engineers Using MATLAB and C, Brooks/Cole, Pacific Grove, CA, 2000. Stoer, J., and Bulirsch, R., Introduction to Numerical Analysis, Springer-Verlag, New York, 1980.

ENSET ORAN

- 61 -

Bibliographie M. Crouzeix, A. L. Mignot. Analyse numérique des équations différentielles. Collections mathématiques appliquées pour la maîtrise. Masson, Paris 1984. K. Dekker, J.-G. Verwer. Stability of Runge-Kutta methods for stiff nonlinear differential equations. CWI Monographs. North-Holland, Amsterdam 1984. E. Hairer, S. P. N_rsett, G. Wanner. Solving ordinary differential equations. I. Springer-Verlag, Berlin, second edition 1993. Nonstiff problems. E. Hairer, G. Wanner. Solving ordinary differential equations. II, volume 14 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition 1996. Stiff and differential-algebraic problems.

ENSET ORAN

- 62 -

Méthode numériques équations différentielles

Recommend Documents