Chapitre 1: Introduction au Calcul Scientifique PDF

Notations Ensembles usuels en mathématiques On désigne généralement les ensemble les plus usuels par une lettre à double barre : N l’ensemble des entiers naturels N∗ l’ensemble des entiers strictement positifs Z l’ensemble des entiers relatifs (positifs, négatifs ou nuls) Z∗ l’ensemble des entiers...

Notations Ensembles usuels en mathématiques On désigne généralement les ensemble les plus usuels par une lettre à double barre : N l’ensemble des entiers naturels N∗ l’ensemble des entiers strictement positifs Z l’ensemble des entiers relatifs (positifs, négatifs ou nuls) Z∗ l’ensemble des entiers 6= 0 ³ ´ p Q l’ensemble des nombres rationnels q, p ∈ Z, q ∈ Z∗ R l’ensemble des réels R∗ l’ensemble des réels autres que 0 C l’ensemble des nombres complexes Rn [x] l’espace vectoriel des polynômes en x de degré inférieur ou égal à n Intervalles Inégalité Notation ensembliste Représentations graphique a b a b a ≤x ≤b [a, b] a b a b a > 1 + 1 2 2 ② les instructions sans chevrons dans une boite grisée sont des bouts de code à écrire dans un fichier 1 print(’Coucou!’) 6 © G. Faccanoni  Introduction au calcul scientifique On peut définir le CALCUL SCIENTIFIQUE comme la discipline qui permet de reproduire sur un ordinateur un phénomène ou un processus décrit par un modèle mathématique. A NALYSE MATHÉMATIQUE P HÉNOMÈNE PHYSIQUE , ÉCO - M ODÈLE MATHÉMATIQUE Bien posé NOMIQUE , BIOLOGIQUE... Mise en équations : Bien conditionné Observation expérimentale équations différentielles, inté- Propriétés de la solution (stabi- Modèle conceptuel grales, stochastiques... lité,...) Solutions analytiques C ALCULS P ROGRAMMATION Langage de programmation (C, A NALYSE NUMÉRIQUE C++, Fortran, Java, Python, Mat- Méthodes de discrétisation lab, Scilab, Octave...) Analyse des algorithmes (rapi- Structure des données dité, précision, souplesse) Implémentation de l’algorithme Estimation des erreurs Optimisation P OST P ROCESSING Visualisation C ALCUL SCIENTIFIQUE Analyse des résultats L’ordinateur est aujourd’hui un outil incontournable pour simuler et modéliser des systèmes complexes, mais il faut encore savoir exprimer nos problèmes (physiques, économiques, biologiques...) en langage formalisé des mathématiques pures sous la forme d’équations mathématiques (différentielles, intégrales...). Nous sommes habitués à résoudre les pro- blèmes de façon analytique, alors que l’ordinateur ne travaille que sur des suites de nombres. On verra qu’il existe souvent plusieurs approches pour résoudre un même problème, ce qui conduit à des algorithmes différents. Un des objectifs de ce cours est de fournir des bases rigoureuses pour développer quelques algorithmes utiles dans la résolution de problèmes en mathématique, économie, physique... 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, aux données, à la représentation sur ordinateur ou encore à la troncature) ; souple : l’algorithme doit être facilement transposable à des problèmes différents. Le choix et l’optimisation des algorithmes numériques mis en pratique sont absolument cruciaux tant pour les calculs de type industriel souvent très répétitifs et devant donc pouvoir être exécutés en un temps très court, que pour les cal- culs de référence pour lesquels la seule limite est la patience de celui qui les fait. Par exemple, en fluidodynamique, en laissant tourner une station de travail pendant quelques jours, les numériciens résolvent des systèmes frisant le milliard d’inconnues. L’expérience montre qu’entre une approche numérique standard et une approche soigneusement réfléchie 7 Introduction au calcul scientifique Jeudi 27 août 2015 et optimisée un gain de temps de calcul d’un facteur 100, voire davantage, est souvent observé. Il est clair qu’on peut pas- ser ainsi, grâce à cet effort, d’un calcul totalement déraisonnable à un calcul parfaitement banal : tout l’enjeu de l’analyse numériques est là ! C’est dire l’importance pour tous scientifique de bien connaître ces méthodes, leurs avantages et leurs limites. p Exemple Calcul de A Sur ordinateur, l’addition de deux entiers peut se faire de façon exacte mais non le calcul d’une racine carrée. On procède alors par approximations successives jusqu’à converger vers la solution souhaitée. Il existe pour cela divers algorithmes. Le suivant est connu depuis l’antiquité (mais ce n’est pas celui que les ordinateurs utilisent). Soit A un nombre réel positif dont on cherche la racine carrée. Désignons par x 0 la première estimation de cette racine (généralement le plus grand entier dont le carré est inférieur à A ; par exemple, si A = 178, alors x 0 = 13 car 132 = 169 < 178 et 142 = 196 > 178) et par ε0 l’erreur associée : p A = x 0 + ε0. Cherchons une approximation de ε0. On a A = (x 0 + ε0 )2 = x 02 + 2x 0 ε0 + ε20. Supposons que l’erreur soit petite face à x 0 , ce qui permet de négliger le terme en ε20 : A ≃ x 02 + 2x 0 ε0. Remplaçons l’erreur ε0 par un ε′0 , qui en est une approximation, de telle sorte que A = x 02 + 2x 0 ε′0. On en déduit que A − x 02 ε′0 = 2x 0 donc la quantité 1 A µ ¶ x 1 ≡ x 0 + ε′0 =+ x0 2 x0 constitue une meilleure approximation de la racine que x 0 (sous réserve que le développement soit convergent). De plus, rien ne nous empêche de recommencer les calculs avec x 1 , puis x 2 , etc., jusqu’à ce que la précision de la machine ne permette plus de distinguer le résultat final de la véritable solution. On peut donc définir une suite, qui à partir d’une estimation initiale x 0 devrait en principe converger vers la solution recherchée. Cette suite est 1 A µ ¶ x k+1 = + xk , x 0 > 0. 2 xk L’algorithme du calcul de la racine carrée devient donc p 1. Démarrer avec une première approximation x 0 > 0 de A. 2. À chaque itération k, calculer la nouvelle approximation x k+1 = 12 xA + x k. ³ ´ k A−x 2 3. Calculer l’erreur associée ε′k+1 = 2x k+1. k+1 4. Tant que l’erreur est supérieure à un seuil fixé, recommencer au point 2 Le tableau ci-dessous illustre quelques itérations de cet algorithme pour le cas où A = 5 : k xk ε′k 0 2.0000000000 0.2360679775 1 2.2500000000 0.0139320225 2 2.2361111111 0.0000431336 3 2.2360679779 0.0000000004 4 2.2360679775 0.0000000000 On voit que l’algorithme converge très rapidement et permet donc d’estimer la racine carrée d’un nombre moyennant un nombre li- mité d’opérations élémentaires (additions, soustractions, divisions, multiplications). Il reste encore à savoir si cet algorithme converge toujours et à déterminer la rapidité de sa convergence. L’analyse numérique est une discipline proche des mathématiques appliquées, qui a pour objectif de répondre à ces questions de façon rigoureuse. Les erreurs Le simple fait d’utiliser un ordinateur pour représenter des nombres réels induit des erreurs. Par conséquent, plutôt que de tenter d’éliminer les erreurs, il vaut mieux chercher à contrôler leur effet. Généralement, on peut identifier plusieurs niveaux d’erreur dans l’approximation et la résolution d’un problème physique. 8 © G. Faccanoni Jeudi 27 août 2015 Au niveau le plus élevé, on trouve l’erreur qui provient du fait qu’on a réduit la réalité physique à un modèle mathéma- tique. De telles erreurs limitent l’application du modèle mathématique à certaines situations et ne sont pas dans le champ du contrôle du Calcul Scientifique. On ne peut généralement pas donner la solution explicite d’un modèle mathématique (qu’il soit exprimé par une in- tégrale, une équation algébrique ou différentielle, un système linéaire ou non linéaire). La résolution par des algorithmes numériques entraîne immanquablement l’introduction et la propagation d’erreurs d’arrondi. De plus, il est souvent né- cessaire d’introduire d’autres erreurs liées au fait qu’un ordinateur ne peut effectuer que de manière approximative des calculs impliquant un nombre infini d’opérations arithmétiques. Par exemple, le calcul de la somme d’une série ne pourra être accompli qu’en procédant à une troncature convenable. On doit donc définir un problème numérique, dont la solution diffère de la solution mathématique exacte d’une erreur, appelée erreur de troncature. La somme des erreurs d’arrondis et de troncature constitue l’erreur de calcul. L’erreur de calcul absolue est la différence entre x, la solution exacte du modèle mathématique, et x̃, la solution obtenue à la fin de la résolution numérique, tandis que (si x 6= 0) l’erreur de calcul relative est définie par l’erreur de calcul absolue divisé par x. Le calcul numérique consiste généralement à approcher le modèle mathématique en faisant intervenir un paramètre de discrétisation, que nous noterons h et que nous supposerons posi- tif. Si, quand h tend vers 0, la solution du calcul numérique tend vers celle du modèle mathématique, nous dirons que le calcul numérique est convergent. Si de plus, l’erreur (absolue ou relative) peut être majorée par une fonction de C h p où C est indépendante de h et où p est un nombre positif, nous dirons que la méthode est convergente d’ordre p. Quand, en plus d’un majorant, on dispose d’un minorant C 1 h p (C 1 étant une autre constante (≤ C ) indépendante de h et p), on peut remplacer le symbole ≤ par ≃. Remarque Pourquoi ne pas faire toujours confiance à une calculatrice Voici un exemple d’un calcul pour lequel une calculatrice donne un résultat faux mais qu’un élève sait faire (en un temps raisonnable). p Calculer sin((1 + 2)200 π). Une calculatrice p (parpexemple celle de Google) renvoie un résultat pbidon, alors que n’importe p 200 quel élève compétent voit que (1+ 2)200 +(1− 2)200 est un entier pair, si bien que p sin((1+ 2) 200 pπ) = − sin((1− 2) π) est extrêmement petit, et 1 on peut même en faire une estimation de tête : sin((1 − 2)200 π) < ((1 − 2)200 π) < 2200 π, or 210 = 1024 ≥ 1000 = 103 donc 1 p 220×10 π ≤ 10−60 π ; de tête on peut dire que sin((1 + 2)200 π) ∈]0; −4 × 10−60 [. La valeur correcte est environ −8.75 × 10−77 , donc cette estimation n’est pas terrible, mais c’est tout de même mieux que les 0.97 retournés par Google. Source : http://www.madore.org/~david/weblog/2014-06.html#d.2014-06-04.2205 © G. Faccanoni 9  1. Résolution d’équations non linéaires Recherche des solutions de l’équation non linéaire f (x) = 0 où f est une fonction donnée Un des problèmes classiques en mathématiques appliquées est celui de la recherche des valeurs pour lesquelles une fonction donnée s’annule. Dans certains cas bien particuliers, comme pour les fonctions x 7→ x + 1, x 7→ cos(2x) ou encore x 7→ x 2 − 2x + 1, le problème est simple car il existe pour ces fonctions des formules qui donnent les zéros explicitement. Toutefois, pour la plupart des fonctions f : R → R il n’est pas possible de résoudre l’équation f (x) = 0 explicitement et il faut recourir à des méthodes numériques. Ainsi par exemple une brève étude de la fonction f (x) = cos(x) − x montre qu’elle possède un zéro à proximité de 0.7 mais ce zéro ne s’exprime pas au moyen de fonctions usuelles et pour en obtenir une valeur approchée il faut recourir à des méthodes numériques. Plusieurs méthodes existent et elles différent pas leur vitesse de convergence et par leur robustesse. Lorsqu’il s’agit de calculer les zéros d’une seule fonction, la vitesse de la méthode utilisée n’est souvent pas cruciale. Cependant, dans certains applications il est nécessaire de calculer les zéros de plusieurs milliers de fonctions et la vitesse devient alors un élément stratégique. Par exemple, si on veut représenter graphiquement l’ensemble de points (x, y) du plan pour lequel x 2 sin(y) + e x+y − 7 = 0, on peut procéder comme suit : pour une valeur donnée de x, on cherche l’ensemble des valeurs de y pour lesquelles x 2 sin(y)+e x+y −7 = 0, i.e. l’ensemble des zéros de la fonction f (y) = x 2 sin(y)+e x+y −7 = 0, et on représente tous les couples obtenus. On choisit une nouvelle valeur pour x et on répète l’opération. Pour obtenir une courbe suffisamment précise, il faut procéder à la recherche des zéros d’un très grand nombre de fonction et il est alors préférable de disposer d’une méthode rapide. Soit f : R → R une fonction continue donnée dont on veut évaluer numériquement un ou plusieurs zéros xb, c’est-à-dire qu’on cherche tous les xb tels que f (b x ) = 0. Les méthodes numériques pour approcher xb consistent à : ① localiser grossièrement le (ou les) zéro(s) de f en procédant à l’étude du graphe de f et/ou à des évaluations qui sont souvent de type graphique ; on note x 0 cette solution grossière ; ② construire, à partir de x 0 , une suite x 1 , x 2 , x 3 ,... telle que limk→∞ x k = xb où f (b x ) = 0. On dit alors que la méthode est convergente. Définition Méthode itérative à deux niveaux On appelle méthode itérative à deux niveaux un procédé de calcul de la forme x k+1 = G(x k ), k = 0, 1, 2,... dans lequel on part d’une valeur donnée x 0 pour calculer x 1 , puis à l’aide de x 1 on calcul x 2 etc. La formule même est dite formule de récurrence. Le procédé est appelé convergent si x k tend vers un nombre fini lorsque k tend vers +∞. Il est bien évident qu’une méthode itérative n’est utile que s’il y a convergence vers les valeurs cherchées. On peut parfaitement envisager des méthodes itératives multi-niveaux, comme par exemples les schémas à trois niveaux dans lesquels on part de deux valeurs données x 0 et x 1 pour calculer x 2 , puis à l’aide de x 1 et x 2 on calcule x 3 etc. Définition Ordre de convergence Soit p un entier positif. On dit qu’une méthode (à deux niveaux) convergente est d’ordre p s’il existe une constante C telle que x − x k |p. x − x k+1 | ≤ C |b |b ou encore x k+1 − xb lim = C. k→∞ (x k − xb)p Si p = 1 (et C < 1) on parle de convergence linéaire, si p = 2 on parle de convergence quadratique. 11 1. Résolution d’équations non linéaires Jeudi 27 août 2015 1.1. Étape ① : localisation des zéros Définition x ) = 0. Il est dit simple si f ′ (b Soit f : R → R une fonction. On dit que xb est un zéro de f si f (b x ) 6= 0, multiple sinon. p Si f est de classe C avec p ∈ N, on dit que xb est un zéro de multiplicité p si ( f (i ) (b x) = 0 ∀i = 0,... , p − 1 (p) f (b x ) 6= 0. Pour localiser grossièrement le (ou les) zéro(s) de f on va d’abord étudier de la fonction f , puis on va essayer d’utiliser un corollaire du théorème des valeurs intermédiaires et le théorème de la bijection afin de trouver un intervalle qui contient un et un seul zéro. Théorème des valeurs intermédiaires Soit f une fonction continue sur un intervalle I = [a; b] de R. Alors f atteint toutes les valeurs intermédiaires entre f (a) et f (b). Autrement dit : ⋆ si f (a) ≤ f (b) alors pour tout d ∈ [ f (a), f (b)] il existe c ∈ [a; b] tel que f (c) = d ; ⋆ si f (a) ≥ f (b) alors pour tout d ∈ [ f (b), f (a)] il existe c ∈ [a; b] tel que f (c) = d. Ce théorème donne alors le corollaire immédiat suivant. Corollaire des zéros d’une fonction continue Soit une fonction continue f : [a, b] → R. Si f (a) · f (b) < 0, alors il existe (au moins un) xb ∈]a, b[ tel que f (b x ) = 0. Ce théorème garantit juste l’existence d’un zéro. Pour l’unicité on essayera d’appliquer le théorème de la bijection dont l’énoncé est rappelé ci-dessous. Théorème de la bijection Soit f une fonction continue et strictement monotone sur un intervalle I de R, alors f induit une bijection de I dans f (I ). De plus, sa bijection réciproque est continue sur I , monotone sur I et de même sens de variation que f. 1.2. Étape ② : construction d’une suite convergente Ayant encadré les zéros de f , la construction de suites qui convergent vers ces zéros peut se faire à l’aide de plusieurs méthodes numériques. Ci-dessous on décrit les méthodes les plus connues et on étudie leurs propriétés (convergence locale vs globale, vitesse de convergence, etc.). 1.2.1. Méthodes de dichotomie (ou bissection), de Lagrange (ou Regula falsi ) Dans les méthodes de dichotomie et de L AGRANGE, à chaque pas d’itération on divise en deux un intervalle donné et on choisit le sous-intervalle où f change de signe. Concrètement, soit f : [a; b] → R une fonction strictement monotone sur un intervalle [a, b]. On suppose que l’équation f (x) = 0 n’a qu’une et une seule solution dans cet intervalle. On se propose de déterminer cette valeur avec une précision donnée. Soit [a 0 , b 0 ] un intervalle dans lequel f (a 0 ) f (b 0 ) < 0 et soit c 0 ∈]a 0 , b 0 [. Si f (a 0 ) f (c 0 ) < 0, alors la racine appartient à l’intervalle [a 0 , c 0 ] et on reprend le procédé avec a 1 = a 0 et b 1 = c 0. Sinon, c’est-à-dire si f (a 0 ) f (c 0 ) ≥ 0 on pose a 1 = c 0 et b 1 = b 0. On construit ainsi une suite d’intervalles emboîtés [a k , b k ]. Les suites a k et b k sont adjacentes 1 et convergent vers xb. Définition Méthodes de dichotomie et de L AGRANGE Soit deux points a 0 et b 0 (avec a 0 < b 0 ) d’images par f de signe contraire (i.e. f (a 0 )· f (b 0 ) < 0). En partant de I 0 = [a 0 , b 0 ], les méthodes de dichotomie et de L AGRANGE produisent une suite de sous-intervalles I k = [a k , b k ], k ≥ 0, avec I k ⊂ I k−1 pour k ≥ 1 et tels que f (a k ) · f (b k ) < 0. ⋆ Dans la méthode de dichotomie, on découpe l’intervalle [a k ; b k ] en deux intervalles de même longueur, i.e. on divise 1. Deux suites (u n )n∈N et (v n )n∈N sont adjacentes si ⋆ (u n ) est croissante, ⋆ (v n ) est décroissante, ⋆ limn (u n − v n ) = 0. Si deux suites sont adjacentes, elles convergent et ont la même limite. 12 © G. Faccanoni Jeudi 27 août 2015 1. Résolution d’équations non linéaires [a k ; b k ] en [a k ; c k ] et [c k ; b k ] où c k est ak + bk ck =. 2 ⋆ Dans la méthode de L AGRANGE, plutôt que de diviser l’intervalle [a k ; b k ] en deux intervalles de même longueur, on découpe [a k ; b k ] en [a k ; c k ] et [c k ; b k ] où c k est l’abscisse du point d’intersection de la droite passant par (a k , f (a k )) et (b k , f (b k )) et l’axe des abscisses, autrement dit c k est solution de l’équation f (b k ) − f (a k ) (c − b k ) + f (b k ) = 0 bk − ak qui est bk − ak a k f (b k ) − b k f (a k ) ck = bk − f (b k ) =. f (b k ) − f (a k ) f (b k ) − f (a k ) Dans les deux cas, pour l’itération suivante, on pose soit [a k+1 ; b k+1 ] = [a k ; c k ] soit [a k+1 ; b k+1 ] = [c k ; b k ] de sorte à ce que f (a k+1 )· f (b k+1 ) < 0. La suite (c k )k∈N converge vers xb puisque la longueur de ces intervalles tend vers 0 quand k tend vers +∞. Soit ε l’erreur maximale qu’on peut commettre, les algorithmes s’écrivent alors comme suit : D ICHOTOMIE : L AGRANGE : Require: a, b > a, ε, f : [a, b] → R Require: a, b > a, ε, f : [a, b] → R k ←0 k ←0 ak ← a ak ← a bk ← b bk ← b ak + bk bk − ak xk ← xk ← ak − f (a k ) 2 f (b k ) − f (a k ) while b k − a k > ε or | f (x k )| > ε do while b k − a k > ε or | f (x k )| > ε do if f (a k ) f (x k ) < 0 then if f (a k ) f (x k ) < 0 then a k+1 ← a k a k+1 ← a k b k+1 ← x k b k+1 ← x k else else a k+1 ← x k a k+1 ← x k b k+1 ← b k b k+1 ← b k end if end if a k+1 + b k+1 b k+1 − a k+1 x k+1 ← x k+1 ← a k+1 − f (a k+1 ) 2 f (b k+1 ) − f (a k+1 ) k ← k +1 k ← k +1 end while end while On n’est pas obligé de stoker tous les intervalles et les itérées, on peut gagner de la mémoire en les écrasant à chaque étape : D ICHOTOMIE : L AGRANGE : Require: a, b > a, ε, f : [a, b] → R Require: a, b > a, ε, f : [a, b] → R a +b b−a x← x ←a− f (a) 2 f (b) − f (a) while b − a > ε or | f (x)| > ε do while b − a > ε or | f (x)| > ε do if f (a) f (x) < 0 then if f (a) f (x) < 0 then b←x b←x else else a←x a←x end if end if a +b b−a x← x ←a− f (a) 2 f (b) − f (a) end while end while Remarque Avec la méthode de la dichotomie, les itérations s’achèvent à la m-ème étape quand |x m − xb| ≤ |I m | < ε, où ε est une tolérance fixée et |I m | désigne la longueur de l’intervalle I m. Clairement I k = b−a 2k , donc pour avoir une erreur |x m − xb| < ε, © G. Faccanoni 13  1. Résolution d’équations non linéaires Jeudi 27 août 2015 on doit prendre le plus petit m qui vérifie µ ¶ b−a ln(b − a) − ln(2) m ≥ log2 =. ε ln(2) Notons que cette inégalité est générale : elle ne dépend pas du choix de la fonction f. Exemple Fond d’investissement Un compte d’épargne donne un taux T ∈ [0; 1] d’intérêt par an avec un virement annuel des intérêts sur le compte. Cela signifie que si le premier janvier 2014 on met v euros sur ce compte, à la fin de la n-ème année (i.e. au 31 décembre 2014 + n) on en retire v(1 + T )n euros. On décide alors d’ajouter au début de chaque année encore v euros. Cela signifie que si on verse ces v euros le premier janvier 2014 + m avec 0 < m < n, au 31 décembre 2014 + n ajoutent v(1 + T )n−m euros. Si à la fin de la n-ème année on en retire un capital de M > v euros, quel est le taux d’intérêt annuel de cet investissement ? À la fin de la n-ème année, le capital versé au premier janvier 2014 est devenu v(1 + T %)n , celui versé au premier janvier 2015 est devenu v(1 + T %)n−1... par conséquente, le capital final M est relié au taux d’intérêt annuel T par la relation n X (1 + T )n − 1 1+T ¡ ¢ M =v (1 + T )k = v =v (1 + T )n − 1. k=1 (1 + T ) − 1 T On en déduit que T est racine de l’équation algébrique non linéaire f (T ) = 0 où 1+T ¡ ¢ f (T ) = v (1 + T )n − 1 − M. T Étudions la fonction f : ⋆ f (T ) > 0 pour tout T > 0, ⋆ limT →0+ f (T ) = nv − M < (n − 1)v, limT →+∞ f (T ) = +∞, ¡ ¢ ⋆ f ′ (T ) = v2 1 + (1 + T )n (T n − 1) > 0 pour tout T > 0 (comparer le graphe de −1/(1 + T )n et de nT − 1) T En étudiant la fonction f on voit que, comme nv < M dès que n > 1, elle admet un unique zéro dans l’intervalle ]0, +∞[ (on peut même prouver qu’elle admet un unique zéro dans l’intervalle ]0, M [). Supposons que v = 1000 ¤ et qu’après 5 ans M est égal à 6000 ¤. En étudiant la fonction f on voit qu’elle admet un unique zéro dans l’intervalle ]0.01, 0.1[. Si on applique la méthode de la dichotomie avec ε = 10−12 , après 37 itérations de la méthode de dichotomie, la méthode converge vers 0.061402411536. On conclut ainsi que le taux d’intérêt T est approximativement égal à 6.14%. Ce résultat a été obtenu en faisant appel à la fonction dichotomie définie à la page 24 comme suit : 2 a = 0.01 # T=1% 3 b = 0.1 # T=10% 4 tol = 1.0e-12 5 maxITER = 50 6 def f(x): 7 return 1000.*(1+x)/x*((1+x)**5-1.)-6000. 8 print dichotomie(f,a,b,tol,maxITER) Exemple Soit f (x) = −39 − 43x + 39x 2 − 5x 3. On cherche a estimer x ∈ [1; 5] tel que f (x) = 0. y D ICHOTOMIE f (5) = 2 f (x) = −39 + (−43 + (39 − 5x)x)x f (3) = 1 f (2.5) = 0.3984375 1 2 f (2) = −0.1875 2.5 3 5 x f (1) = −1 I 0 = [1; 5] I 1 = [1; 3] I 2 = [2; 3] I 3 = [2; 2.5] 14 © G. Faccanoni  Jeudi 27 août 2015 1. Résolution d’équations non linéaires y L AGRANGE f (5) = 2 f (x) = −39 + (−43 + (39 − 5x)x)x f (2.3̄) = 0.197530824 1 2.11 f (2.11) = −0.06002 5 x 2.3̄ f (1) = −1 I 0 = [1; 5] I 1 = [1; 2.3̄] I 2 = [2.11; 2.3̄] La méthode de dichotomie est simple mais elle ne garantit pas une réduction monotone de l’erreur d’une itération à l’autre : tout ce dont on est assuré, c’est que la longueur de l’intervalle de recherche est divisée par deux à chaque étape. Par conséquent, si le seul critère d’arrêt est le contrôle de la longueur de I k , on risque de rejeter de bonnes approximations de xb. En fait, cette méthode ne prend pas suffisamment en compte le comportement réel de f. Il est par exemple frappant que la méthode ne converge pas en une seule itération quand f est linéaire (à moins que le zéro xb ne soit le milieu de l’intervalle de recherche initial). 1.2.2. Méthode de la sécante Soit f : R → R une fonction continue et soit xb ∈ [a, b] un zéro de f. La méthode de la sécante est une variante de la méthode de L AGRANGE dans laquelle on ne demande plus à ce que le zéro soit entre a k et b k. Pour calculer x k+1 on prend l’intersection de l’axe des abscisses avec la droite passant par les points (x k , f (x k )) et (x k−1 , f (x k−1 )), i.e. on cherche x solution du système linéaire ( f (x )− f (x ) y = xk −x k−1 (x − x k ) + f (x k ), k k−1 y = 0, ce qui donne x k − x k−1 x = xk − f (x k ). f (x k ) − f (x k−1 ) Définition Méthode de la Sécante Il s’agit d’une méthode à trois niveaux : approcher les zéros de f se ramène à calculer la limite de la la suite récurrente    x 0 donné,   x 1 donné,   x k − x k−1  x k+1 = x k − f (x k ), f (x k ) − f (x k−1 ) p 1+ 5 Cette méthode a ordre de convergence 2. 1.2.3. Méthodes de point fixe En s’amusant avec une calculatrice de poche ou avec le code python ci-dessous 1 import math 2 x = 1 3 for i in range (1,100): 4 −−−→x = math.cos(x) 5 −−−→print("x[{0}]={1}".format(i,x)) on peut vérifier qu’en partant de la valeur 1 et en appuyant plusieurs fois de suite sur la touche «cosinus», on obtient cette suite de valeurs : x 0 = 1, © G. Faccanoni 15 1. Résolution d’équations non linéaires Jeudi 27 août 2015 x 1 = cos(x 0 ) = 0.540302305868, x 2 = cos(x 1 ) = 0.857553215846, x 3 = cos(x 2 ) = 0.654289790498,... x 55 = 0.739085133171,... x 100 = 0.739085133215 qui tend vers la valeur 0.73908513.... En effet, on a par construction x k+1 = cos(x k ) pour k = 0, 1,... (avec x 0 = 1). Si cette suite converge, sa limite ℓ satisfait l’équation cos(ℓ) = ℓ. Pour cette raison, ℓ est appelé point fixe de la fonction cosinus. Définition Point fixe Soit ϕ : R → R une fonction. Si xb ∈ R est tel que ϕ(b x ) = xb, on dit que xb est un point fixe de ϕ (l’image de xb par ϕ est lui-même). On peut se demander comment exploiter cette procédure pour calculer les zéros d’une fonction donnée. Remarquons qu’on peut voir ℓ comme un point fixe du cosinus, ou encore comme un zéro de la fonction f (x) = x − cos(x). La méthode proposée fournit donc un moyen de calculer les zéros de f. Précisons ce principe : soit f : [a, b] → R la fonction dont on cherche le zéro. Il est toujours possible de transformer le problème (Pb-1) “chercher x tel que f (x) = 0” en un problème équivalent (i.e. admettant les mêmes solutions) (Pb-2) “chercher x tel que x − ϕ(x) = 0”. Pour que les deux problèmes soient équivalent, la fonction auxiliaire ϕ : [a, b] → R doit être choisie de manière à ce que ϕ(bx ) = xb si et seulement si f (b x ) = 0 dans [a; b] (on dit alors que le problème (Pb-2) est consistant avec le problème (Pb-1)). Clairement, il existe une infinité de manières pour opérer cette transformation. Par exemple, on peut poser ϕ(x) = x − f (x) ou plus généralement ϕ(x) = x + γ f (x) avec γ ∈ R∗ quelconque. On peut même remplacer γ par une fonction de x pour autant qu’elle ne s’annule pas. Définition Méthode de point fixe Supposons que xb ∈ R soit un point fixe de ϕ. La méthode de point fixe consiste en la construction d’une suite (x k )k∈N définie par récurrence comme suit : ( x0 donné, x k+1 = ϕ(x k ) ∀k ∈ N. Naturellement une telle suite n’est pas forcement convergente. Par contre, si elle converge, c’est-à-dire si la suite x k a une limite que nous notons ℓ, et si ϕ est continue, alors cette limite est nécessairement un point fixe de ϕ puisque µ ¶ ℓ = lim x k+1 = lim ϕ(x k ) = ϕ lim x k = ϕ(ℓ). k→∞ k→∞ k→∞ On utilise alors l’algorithme itératif suivant pour construire la suite (comme l’ordinateur ne peux pas construire une infinité de termes, on calcul les premiers termes de la suite et on s’arrête dès que la différence entre deux éléments de la suite est inférieure à une tolérance ε > 0 donnée) : Require: x 0 , ε, ϕ : [a, b] → R k ←0 x 1 ← x 0 + 2ε while |x k+1 − x k | > ε do x k+1 ← ϕ(x k ) k ← k +1 end while On va maintenant s’intéresser à la convergence de la suite construite par une méthode de point fixe. Théorème Convergence (globale) des itérations de point fixe Considérons une fonction ϕ : [a; b] → R. On se donne x 0 ∈ [a; b] et on considère la suite x k+1 = ϕ(x k ) pour k ≥ 0. Si les deux conditions suivantes sont satisfaites : 1. condition de stabilité : ϕ(x) ∈ [a, b] pour tout x ∈ [a, b] 16 © G. Faccanoni Jeudi 27 août 2015 1. Résolution d’équations non linéaires y y b b x1 f x3 x5 x7 x6 x4 x6 x5 x4 x2 x3 f x2 x1 ϕ([a; b]) ⊂ [a; b] ϕ([a; b]) ⊂ [a; b] 0 ≤ ϕ′ (x) < 1 −1 < ϕ′ (x) ≤ 0 convergence a a convergence a x0 x1 x2 x 3 x 4x 5 b x a x0 x2 x 4x 6 x 5 x 3 b x x1 F IGURE 1.1.: Interprétation géométrique du théorème de point fixe 2. condition de contraction stricte : il existe K ∈ [0; 1[ tel que |ϕ(x) − ϕ(y)| ≤ K |x − y| pour tout x, y ∈ [a, b] alors ⋆ ϕ est continue, ⋆ ϕ a un et un seul point fixe x b dans [a, b], ⋆ la suite x k+1 = ϕ(x k ) converge vers x b pour tout choix de x 0 dans [a, b]. Démonstration. Continuité La condition de contraction stricte implique que ϕ est continue puisque, si on prend une suite (y k )k∈N ∈ [a, b] qui converge vers un élément x de [a, b], alors nous avons |ϕ(x)−ϕ(y n )| ≤ K |x − y n | et par suite limk→∞ ϕ(y k ) = ϕ(x). Existence Commençons par prouver l’existence d’un point fixe de ϕ. La fonction g (x) = ϕ(x) − x est continue dans [a, b] et, grâce à la condition de stabilité, on a g (a) = ϕ(a) − a ≥ 0 et g (b) = ϕ(b) − b ≤ 0. En appliquant le théorème des valeurs intermédiaires, on en déduit que g a au moins un zéro dans [a, b], i.e. ϕ a au moins un point fixe dans [a, b]. Unicité L’unicité du point fixe découle de la condition de contraction stricte. En effet, si on avait deux points fixes distincts xb1 et xb2 , alors |b x 1 ) − ϕ(b x 1 − xb2 | = |ϕ(b x 2 )| ≤ K |b x 1 − xb2 | < |b x 1 − xb2 | ce qui est impossible. Convergence Prouvons à présent que la suite x k converge vers l’unique point fixe xb quand k tend vers +∞ pour toute donnée initiale x 0 ∈ [a; b]. On a 0 ≤ |x k+1 − xb| = |ϕ(x k ) − ϕ(b x )| ≤ K |x k − xb| où K < 1 est la constante de contraction. En itérant k + 1 fois cette relation on obtient |x k+1 − xb| ≤ K k+1 |x 0 − xb|, i.e., pour tout k ≥ 0 |x k+1 − xb| ≤ K k+1. |x 0 − xb| En passant à la limite quand k tend vers +∞ on obtient |x k+1 − xb| tend vers zéro. Il est important de disposer d’un critère pratique assurant qu’une fonction ϕ est contractante stricte. Pour cela, rappelons quelques définitions. Théorème Si ϕ : [a; b] → [a; b] est de classe C 1 ([a, b]) et si |ϕ′ (x)| < 1 pour tout x ∈ [a, b], alors la condition de contraction stricte est satisfaite avec K = max|ϕ′ (x)|. [a;b] © G. Faccanoni 17 1. Résolution d’équations non linéaires Jeudi 27 août 2015 Démonstration. Considérons la fonction affine g qui transforme l’intervalle [0; 1] dans l’intervalle [x; y] : g : [0; 1] → [x; y] t 7→ x + t (x − y) Alors Zy Z1 Z1 ϕ(x) − ϕ(y) = ϕ′ (ζ) dζ = ϕ′ (g (t ))g ′ (t ) dt = (x − y) ϕ′ (x + t (x − y)) dt x 0 0 et donc ¯ Z1 ¯ Z ¯ ¯ ¯ ′ ¯ ¯ ¯ 1¯ ′ ¯ ¯ ¯ ¯ϕ(x) − ϕ(y)¯ = ¯(x − y) ϕ (x + t (x − y)) dt ¯ ≤ ¯(x − y)¯ ¯ϕ (x + t (x − y))¯ dt ≤ K ¯(x − y)¯. ¯ ¯ 0 0 Le théorème de convergence globale assure la convergence, avec un ordre 1, de la suite (x k )k∈N vers le point fixe xb pour tout choix d’une valeur initiale x 0 ∈ [a; b]. Mais en pratique, il est souvent difficile de déterminer à priori l’intervalle [a; b] ; dans ce cas, le résultat de convergence locale suivant peut être utile. Théorème d’O STROWSKI ou de convergence (locale) des itérations de point fixe Soit [a; b] ∈ R et ϕ : [a; b] → [a; b] une application de classe C 1 ([a; b]). Soit xb ∈ [a; b] un point fixe de ϕ. On peut distinguer trois cas : ① Soit |ϕ′ (b x )| < 1. Par continuité de ϕ′ il existe un intervalle [b x − ε; xb + ε] ⊂ [a; b] sur lequel |ϕ′ (b x )| < K < 1, donc ϕ est contractante stricte sur cet intervalle. On a nécessairement ϕ([b x −ε; xb +ε]) ⊂ [bx −ε; xb +ε] et par conséquent la suite (x k )k∈N converge vers xb pour tout x 0 ∈ [bx − ε; xb + ε]. On dit que xb est un point fixe attractif. De plus, ⋆ si 0 < ϕ′ (b x ) < 1 la suite converge de façon monotone, c’est-à-dire l’erreur x k − xb garde un signe constant quand k varie ; ⋆ si −1 < ϕ′ (b x ) < 0 la suite converge de façon oscillante, c’est-à-dire l’erreur x k − xb change de signe selon la parité de k. ② Si |ϕ′ (b x )| > 1, alors il n’existe aucun intervalle [b x − ε; xb + ε] ⊂ [a; b] tel que la suite (x k )k∈N converge vers xb pour tout x 0 ∈ [b x − ε; xb + ε] à l’exception du cas x 0 = xb. On dit que xb est un point fixe répulsif. De plus, ⋆ si ϕ′ (bx ) > 1 la suite diverge de façon monotone, ⋆ si ϕ′ (bx ) < −1 la suite diverge en oscillant. ③ Si |ϕ′ (b x )| = 1 on ne peut en général tirer aucune conclusion : selon le problème considéré, il peut y avoir conver- gence ou divergence. Exemple ⋆ La fonction ϕ(x) = cos(x) vérifie toutes les hypothèses du théorème d’O STROWSKI : elle est de classe C ∞ (R) et |ϕ′ (x b)| = | sin(xb)| ≃ 0.67 < 1, donc il existe par continuité un intervalle [c, d ] qui contient x b tel que |ϕ ′ (x b )| < 1 pour x ∈ [c, d ]. p p ⋆ La fonction ϕ(x) = x 2 − 1 possède deux points fixes x b1 p= (1 + 5)/2 et xb2 = (1 − 5)/2 mais ne vérifie l’hypothèse du théorème d’O STROWSKI pour aucun d’eux puisque |ϕ(xb1,2 )| = |(1 ± 5)/2| > 1. Les itérations de point fixe ne convergent pas. ⋆ La fonction φ(x) = x − x 3 admet x b = 0 comme point fixe. On a φ′ (xb) = 1 et x k → xb pour tout x 0 ∈ [−1; 1] car ⋆ si x 0 = ±1 alors x k = x b pour tout k ≥ 1, ⋆ si x 0 ∈]−1, 1[ alors on montre que x k ∈]−1, 1[ pour tout k ≥ 1. De plus, la suite est monotone décroissante si 0 < x 0 < 1, monotone croissante si −1 < x 0 < 0 donc elle converge vers ℓ ∈ [−1; 1]. Les uniques candidats limites sont les solutions de l’équation ℓ = φ(ℓ) et par conséquente x k → 0. ⋆ La fonction φ(x) = x + x 3 admet aussi x b = 0 comme point fixe. À nouveau φ′ (xb) = 1 mais dans ce cas la suite diverge pour tout choix de x 0 6= 0. Le théorème d’O STROWSKI dit que, de manière générale, la méthode de point fixe ne converge pas pour des valeurs arbi- traires de x 0 , mais seulement pour des valeurs suffisamment proches de xb, c’est-à-dire appartenant à un certain voisinage de xb. Au premier abord, cette condition semble inutilisable : elle signifie en effet que pour calculer xb (qui est inconnu), on devrait partir d’une valeur assez proche de xb. En pratique, on peut obtenir une valeur initiale x 0 en effectuant quelques itérations de la méthode de dichotomie ou en examinant le graphe de f. Si x 0 est convenablement choisi alors la méthode de point fixe converge. Proposition Calcul de l’ordre de convergence d’une méthode de point fixe Soit xb un point fixe d’une fonction ϕ ∈ C p+1 pour un entier p ≥ 1 dans un intervalle [a; b] contenant xb. Si ϕ(i ) (bx ) = 0 pour 1 ≤ i ≤ p et ϕ(p+1) (b x ) 6= 0, alors la méthode de point fixe associée à la fonction d’itération ϕ est d’ordre p + 1 et x k+1 − xb ϕ(p+1) (b x) lim =. k→+∞ (x k − x b)p+1 (p + 1)! 18 © G. Faccanoni Jeudi 27 août 2015 1. Résolution d’équations non linéaires Démonstration. Écrivons le développement de TAYLOR avec le reste de L AGRANGE de ϕ en x = xb : X p ϕ(i ) (b x) ϕ(p+1) (ξ) x) + ϕ(x) = ϕ(b (x − xb)i + (x − xb)p+1 i =1 i! (p + 1)! x ) = xb et ϕ(i ) (b où ξ est entre x et xb. Comme ϕ(b x ) = 0 pour 1 ≤ i ≤ p, cela se simplifie et on a ϕ(p+1) (ξ) ϕ(x) = xb + (x − xb)p+1. (p + 1)! En évaluant l’expression ainsi trouvée en x k et sachant que ϕ(x k ) = x k+1 , on a alors ϕ(p+1) (ξ) x k+1 − xb = (x k − xb)p+1. (p + 1)! Lorsque k → +∞, x k tend vers xb et donc ξ, qui se trouve entre x k et xb, tend vers xb aussi. Alors x k+1 − xb ϕ(p+1) (ξ) ϕ(p+1) (b x) lim = lim =. k→+∞ (x k − x b) p+1 k→+∞ (p + 1)! (p + 1)! ϕ(p+1) (b x) Pour un ordre p fixé, la convergence de la suite vers xb est d’autant plus rapide que (p+1)! est petit. Méthodes de point fixe particulièrement connues Soit f : [a, b] → R une fonction continue (continûment dérivable pour la méthode de la corde 2 et la méthode de N EW- TON ) et soit x b un zéro de f. Supposons que l’on connaisse une valeur x 0 proche de xb. Approcher les zéros de f se ramène au problème de la détermination des points fixes de la fonction ϕ, ce qui se fait en construisant la suite récurrente ( x 0 donné, x k+1 = ϕ(x k ). Pour choisir la fonction ϕ il est nécessaire de prendre en compte les informations données par les valeurs de f et, éven- tuellement, par sa dérivée f ′ ou par une approximation convenable de celle-ci (si f est différentiable). Écrivons pour cela le développement de TAYLOR de f en xb au premier ordre : f (b x − x) f ′ (ξ) où ξ est entre xb et x. Le problème x ) = f (x)+(b “chercher xb tel que f (b x ) = 0” devient alors x − x) f ′ (ξ) = 0”. “chercher xb tel que f (x) + (b Cette équation conduit à la méthode itérative suivante : “pour tout k ≥ 0, étant donné x k , déterminer x k+1 en résolvant l’équation f (x k )+(x k+1 − x k )q k = 0, où q k est égal à f ′ (ξk ) (ou en est une approximation) avec ξk un point entre x k et x k+1.” La méthode qu’on vient de décrire revient à chercher l’intersection entre l’axe des x et la droite de pente q k passant par le point (x k , f (x k )), ce qui s’écrit sous la forme d’une méthode de point fixe avec f (x k ) x k+1 = ϕ(x k ) ≡ x k − , k ≥ 0. qk Considérons maintenant quatre choix particuliers de q k et donc de ϕ qui définissent des méthodes célèbres : b−a b−a Méthode de la Corde 1 : qk = =⇒ ϕ(x) = x − f (x) f (b) − f (a) f (b) − f (a) f (x) Méthode de la Corde 2 : q k = f ′ (x 0 ) =⇒ ϕ(x) = x − ′ f (x 0 ) f (x) Méthode de N EWTON : q k = f ′ (x k ) =⇒ ϕ(x) = x − ′ f (x) Proposition Si la méthode de la corde converge, elle converge à l’ordre 1 ; si la méthode de N EWTON converge, elle converge à l’ordre 2 si la racine est simple, à l’ordre 1 sinon. © G. Faccanoni 19 1. Résolution d’équations non linéaires Jeudi 27 août 2015 Démonstration. f (b)− f (a) Méthodes de la Corde qk = (méthode 1) ou q k = f ′ (x 0 ) (méthode 2). Si f ′ (b b−a x ) = 0 alors ϕ′ (b x ) = 1 et on ne peut pas ′ assurer la convergence de la méthode. Autrement, la condition |ϕ (b x )| < 1 revient à demander que 0 < f ′ (b x )/q k < 2. Ainsi la pente de la corde doit avoir le même signe que f ′ (bx ) et, pour la méthode 1, l’intervalle de recherche [a; b] doit être tel quel f (b) − f (a) b−a 1, alors la méthode n’est plus du second ordre. En effet, f (x) = (x − xb)m h(x) où h est une fonction telle que h(b x ) 6= 0. On a alors f (x) (x − xb)h(x) ϕ(x) = 1 − = 1− , f ′ (x) mh(x) + (x − xb)h ′ (x) ¡ ¢ h(x) m(m − 1)h(x) + 2(x − xb)h ′ (x) + (x − xb)2 h ′′ (x) 1 ϕ′ (x) = ¡ ¢2 , ϕ′ (b x) = 1 −. mh(x) + (x − xb)h ′ (x) m Si la valeur de m est connue a priori, on peut retrouver la convergence quadratique en modifiant la méthode de N EWTON comme suit : f (x) ϕ(x) = x − m ′. f (x) Attention À noter que même si la méthode de N EWTON permet en général d’obtenir une convergence quadratique, un mauvais choix de la valeur initiale peut provoquer la divergence de cette méthode (notamment si la courbe représentative de f présente au point d’abscisse x 0 un tangente à peu près horizontale). D’où l’importance d’une étude préalable soignée de la fonction f (cette étude est d’ailleurs nécessaire pour toute méthode de point fixe). Remarque Interprétation géométrique de la méthode de N EWTON et des méthodes de la corde Soit f : R → R une fonction continûment dérivable et soit xb un zéro simple de f , c’est-à-dire f (b x ) = 0 et f ′ (b x ) 6= 0. Suppo- sons que l’on connaisse une valeur x k proche de xb. Pour calculer x k+1 on prend l’intersection de l’axe des abscisses avec la droite tangente au graphe de f passant par le point (x k , f (x k )), i.e. on cherche x solution du système linéaire ( y = f ′ (x k )(x − x k ) + f (x k ), y = 0. On obtient f (x k ) x = xk − f ′ (x k ) ce qui correspond à la méthode de N EWTON. 20 © G. Faccanoni Jeudi 27 août 2015 1. Résolution d’équations non linéaires y f y = f ′ (x 1 )(x − x 1 ) + f (x 1 ) xb x2 x1 x0 x y = f ′ (x 0 )(x − x 0 ) + f (x 0 ) Soit f : R → R une fonction continue et soit xb ∈ [a, b] un zéro de f. Cette fois-ci, pour calculer x k+1 on prend l’intersection de l’axe des abscisses avec la droite passant par le point (x k , f (x k )) et parallèle à la droite passant par les points (a, f (a)) et (b, f (b)), i.e. on cherche x solution du système linéaire ( f (b)− f (a) y= b−a (x − x k ) + f (x k ), y = 0, ce qui donne b−a x = xk − f (x k ). f (b) − f (a) Il s’agit de la méthode de la corde 1. Cette méthode permet d’éviter qu’à chaque itération on ait à évaluer f ′ (x k ) car on f (b)− f (a) remplace f ′ (x k ) par b−a. Une variante de la méthode de la corde consiste à calculer x k+1 comme l’intersection entre l’axe des abscisses et la droite passant par le point (x k , f (x k )) et parallèle à la droite tangente au graphe de f passant par le point (x 0 , f (x 0 )), i.e. on cherche x solution du système linéaire ( y = f ′ (x 0 )(x − x k ) + f (x k ), y = 0, ce qui donne f (x k ) x = xk − f ′ (x 0 ) Dans cette variante on remplace f ′ (x k ) par f ′ (x 0 ). y y f f y = f ′ (x 0 )(x − x 0 ) + f (x 0 ) f (b)− f (a) y = b−a (x − x 0 ) + f (x 0 ) f (b)− f (a) y = b−a (x − x 1 ) + f (x 1 ) xb a x1 x2 b x1 x2 x0 x xb x0 x y = f ′ (x 0 )(x − x 1 ) + f (x 1 ) Exemple On se trouve en possession d’une calculatrice qui ne sait effectuer que les opérations addition, soustraction et multiplication. Lorsque a > 0 est donné, on veut calculer sa valeur réciproque 1/a. Le problème peut être ramené à résoudre l’équation x = 1/a ce qui équivaut à chercher le zéro de la fonction f : R+ ∗ →R 1 x 7→ − a x Selon la formule de N EWTON on a x k+1 = (1 + a)x k + x k2 , © G. Faccanoni 21 1. Résolution d’équations non linéaires Jeudi 27 août 2015 une récurrence qui ne requiert pas de divisions. Pour a = 7 et partant de x 0 = 0.2 par exemple, on trouve x 1 = 0,12, x 2 = 0,139 2,x 3 = 0,142 763 520 0, x 4 = 0,142 857 081 5, etc. Cette suite converge vers 1/7 ≃ 0,142 857 142 857. Exemple Comparaison des méthodes de N EWTON pour différentes formulation de la fonction initiale Dans R∗ + on veut résoudre l’équation x = e 1/x. (1.1) En transformant l’équation donnée de différentes manières, on arrive à différentes formules de récurrence : 1. L’équation (1.1) équivaut à chercher le zéro de la fonction f : R∗ + →R x 7→ x − e 1/x En utilisant la méthode de N EWTON on trouve la formule itérative f (x ) x − e 1/xk x − e 1/xk x k+1 = x k − ′ k = x k − k 1/x = x k − x k2 k2. f (x k ) 1+ e k x + e 1/xk k x k2 2. Si on pose y = 1/x, alors on a l’équivalence x = e 1/x ⇐⇒ y = e −y , donc la solution x de l’équation (1.1) est la réciproque du zéro de la fonction g : R∗ + →R y 7→ 1 − ye y En utilisant la méthode de N EWTON on trouve la formule itérative g (y ) 1 − yk e yk e −y k − y k y k+1 = y k − ′ k = y k − = yk + g (y k ) −(1 + y k )e y k 1 + yk et x k = 1/y k. 3. L’équation (1.1) est encore équivalente à chercher le zéro de la fonction h : R∗ + →R x 7→ 1 − x ln(x) En utilisant la méthode de N EWTON on trouve la formule itérative h(x ) 1 − x k ln(x k ) 1 + xk x k+1 = x k − ′ k = x k + =. h (x k ) 1 + ln(x k ) 1 + ln(x k ) La représentation graphique de f montre qu’il n’existe qu’une seule racine. Comme f (1.7) f (1.9) < 0, elle se trouve dans l’intervalle [1.7; 1.9]. En partant de x 0 = 1.8 on trouve les suites suivantes : Formule 1 Formule 2 Formule 3 x1 = 1,762 878 141 2 1.7418849724 1.7634610883 x2 = 1.7632228030 1.7751466845 1.7632228446 x3 = 1.7632228344 1.7564077294 1.7632228344 La solution est x ≃ 1,763 222 834 35. Critères d’arrêt Supposons que (x n )n∈N soit une suite qui converge vers xb zéro de la fonction f. Nous avons le choix entre deux types de critères d’arrêt pour interrompre le processus itératif d’approximation de xb : ceux basés sur le résidu et ceux basés sur l’incrément. Nous désignerons par ε une tolérance fixée pour le calcul approché de xb et par e n = xb − x n l’erreur absolue. Nous supposerons de plus f continûment différentiable dans un voisinage de la racine. Contrôle du résidu : les itérations s’achèvent dès que | f (x n )| < ε. Il y a des situations pour lesquelles ce test s’avère trop restrictif ou, au contraire, trop optimiste. ⋆ si | f ′ (b x )| ≃ 1 alors |e n | ≃ ε : le test donne donc une indication satisfaisante de l’erreur ; ⋆ si | f ′ (b x )| ≪ 1, le test n’est pas bien adapté car |e n | peut être assez grand par rapport à ε (voir la figure 1.2 à droite) ; ⋆ si enfin | f ′ (b x )| ≫ 1 alors |e n | ≪ ε et le test est trop restrictif (voir la figure 1.2 à gauche). Contrôle de l’incrément : les itérations s’achèvent dès que |x n+1 − x n | < ε. Soit (x n )n∈N la suite produite par la méthode de point fixe x n+1 = ϕ(x n ). Comme xb = ϕ(b x ) et x n+1 = ϕ(x n ), si on développe au premier ordre on sait qu’il existe ξn ∈ I xb,xn tel que x ) − ϕ(x n ) = ϕ′ (ξn )(b e n+1 = xb − x n+1 = ϕ(b x − x n ) = ϕ′ (ξn )e n 22 © G. Faccanoni Jeudi 27 août 2015 1. Résolution d’équations non linéaires y y f f (x k ) f xk xk f (x k ) xb xb ek x ek x F IGURE 1.2.: Deux situations pour lesquelles le résidu e k = x k − xb est un mauvais estimateur d’erreur : | f ′ (x)| ≫ 1 (à gauche), | f ′ (x)| ≪ 1 (à droite), pour x dans un voisinage de xb. où I xb,xn est l’intervalle d’extrémités xb et x k. En utilisant l’identité x − x n+1 ) + (x n+1 − x n ) = e n+1 + (x n+1 − x n ) = ϕ′ (ξn )e n + (x n+1 − x n ), e n = (b on en déduit que x n+1 − x n en =. 1 − ϕ′ (ξn ) Par conséquent, ce critère fournit un estimateur d’erreur satisfaisant si ϕ′ (x) ≃ 0 dans un voisinage de xb. C’est le cas notamment des méthodes d’ordre 2, dont la méthode de N EWTON. Cette estimation devient d’autant moins bonne que ϕ′ s’approche de 1. Notons d’ailleurs que si la méthode de point fixe converge avec K < 1 et si on considère le critère d’arrêt |x n+1 − x n | < ε alors ε |e n | = |x n − xb| ≤ ≤ 2ε. 1−K En effet, il suffit de considérer les inégalités suivantes : |e n+1 | = |x n+1 − xb| = |ϕ(x n ) − ϕ(b x )| ≤ K |x n − xb| = K |e n |, |e n+1 | = |x n+1 − xb| = |x n+1 − x n + x n − xb| ≥ |x n − xb| − |x n+1 − x n | = |e n | − |x n+1 − x n | > |e n | − ε © G. Faccanoni 23  1. Résolution d’équations non linéaires Jeudi 27 août 2015 ✴✴✴✴✴✴✴✴✴✴✴✴✴ Codes Python ✴✴✴✴✴✴✴✴✴✴✴✴ dichotomie, lagrange, newton et point_fix sont quatre fonctions (informatiques) qui renvoient la valeur approchée du zéro d’une fonction (mathématique) f. En paramètre elles reçoivent f, la fonction dont on cherche la racine, a et b sont les extrémités de l’intervalle de recherche pour les méthodes de dichotomie et de L AGRANGE, x_init est la donnée initiale pour les méthodes de N EWTON et de point fixe, maxITER est le nombre maximal d’itérations et tol est la tolérance. Méthodes numériques. 1 #!/usr/bin/python 2 #-*- coding: Utf-8 -*- 3 4 import math, sys 5 6 def dichotomie(f,a,b,tol,maxITER): 7 −−−→fa = f(a) 8 −−−→if abs(fa)tol)) and (ktol) and (ktol) and (k