MÃ©thodes numeriques pour la physique - CINaM

M´ ethodes numeriques pour la physique

Licences de Physique Faculté des Sciences de Luminy

Hubert Klein ([email protected]) CiNaM UPR CNRS 3118 case 913 faculté des Sciences de Luminy 13288 Marseille cedex 9

Table des mati` eres Introduction

4

Repr´ esentation des nombres dans un ordinateur

6

I

9

Tri et mod´ elisation de donn´ ees

1 Algorithmes de tri de 1.1 Tri à bulle . . . . . 1.2 Tri par insertion . 1.3 Tri shell . . . . . . 1.4 Tri rapide . . . . .

. . . .

10 11 12 13 15

2 Cacul de quelques grandeurs statistiques usuelles 2.1 Moyenne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 variance, écart-type . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Médian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17 17 18 18

3 R´ egression lin´ eaire par la m´ ethode des moindres carr´ es 3.1 Définition du problème . . . . . . . . . . . . . . . . . . . . 3.2 Implémentation . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Exemple d’application : ajustement d’une fonction linéaire 3.3.1 prise en compte de l’incertitude sur les mesures . .

19 19 19 20 20

II

donn´ ees . . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . y .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . . . . . . = a.x + b . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

M´ ethodes num´ eriques

4 Recherche des racines d’une fonction 4.1 Encadrement . . . . . . . . . . . . . 4.2 Dichotomie . . . . . . . . . . . . . . 4.2.1 Convergence de la méthode . 4.3 Méthode de Newton - Raphson . . . 4.3.1 Convergence de la méthode .

23 . . . . . 2

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

24 24 26 26 28 28

` TABLE DES MATIERES

3

5 Minimisation ou maximisation d’une 5.1 Recherche par la section d’or (golden 5.1.1 Principe . . . . . . . . . . . . 5.1.2 Convergence . . . . . . . . .

fonction 31 section search) . . . . . . . . . . . . . . . . 31 . . . . . . . . . . . . . . . . . . . . . . . . . 32 . . . . . . . . . . . . . . . . . . . . . . . . . 35

6 Interpolation 6.1 Interpolation linéaire . . . . . . . . . . 6.2 Interpolation quadratique . . . . . . . 6.3 Généralisation ` a l’ordre n . . . . . . . 6.3.1 Du bon usage de l’interpolation

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

36 36 38 39 40

7 Int´ egration de fonctions 7.1 Intégration par la méthode des rectanges . . . 7.1.1 Incertitude de la méthode . . . . . . . 7.2 Intégration par la méthode des trapèzes . . . 7.2.1 Incertitude de la méthode . . . . . . . 7.3 Intégration par la méthode de Simpson . . . . 7.3.1 Incertitude de la méthode . . . . . . . 7.4 Méthode de Romberg . . . . . . . . . . . . . 7.5 Détermination d’un pas optimal d’intégration

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

41 41 43 43 44 44 45 45 46

8 D´ erivation num´ erique 8.1 Formule d’Euler . . . . . . . 8.1.1 Erreur intrinsèque . 8.2 Dérivation au demi-pas . . . 8.2.1 Stabilité numérique

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

48 48 48 49 49

. . . . . . . . .

51 52 53 54 55 55 56 57 58 59

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

9 R´ esolution d’´ equations diff´ erentielles lin´ eaires 9.1 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.1.1 Précision et stabilité de la méthode . . . . . . . . . . . . . 9.1.2 Exemple d’application : la loi de décroissance radioactive 9.2 Méthode du saut de grenouille . . . . . . . . . . . . . . . . . . . 9.2.1 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.3 Algorithme Verlet-vitesses . . . . . . . . . . . . . . . . . . . . . . 9.4 Méthode de Runge Kutta ` a l’ordre 2 . . . . . . . . . . . . . . . . 9.4.1 Stabilité de la méthode . . . . . . . . . . . . . . . . . . . 9.5 Méthode Runge Kutta ` a l’ordre 4 . . . . . . . . . . . . . . . . . . A R´ ef´ erences bibliographiques et informatiques

M´ ethodes Informatiques pour la physique

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

62

H. Klein

Introduction Ce document n’a pas pour but d’être un cours d’informatique, encore moins un cours d’analyse numérique. L’idée en est plutˆ ot d’utiliser la puissance de calcul d’un ordinateur pour traiter des problèmes physique très variés. Pourquoi utiliser un ordinateur ? Pour gagner du temps par exemple lors d’un calcul fastidieux que l’on pourrait réaliser à la main, mais aussi dans des cas o` u le problème ne peut pas être traité avec un crayon et une feuille de papier (certaines équations n’ont par exemple pas de solutions analytiques). Nous aborderons au cours de ce cours les méthodes classiques de calcul numérique qui nous permettrons d’intégrer des fonctions, de calculer leurs racines, d’intégrer des systèmes d’équations différentielles. . .Dans une deuxième partie nous aborderons l’usage d’un ordinateur pour traiter des données. Cet aspect statistique du calcul numérique est particulièrement utilisé par les expérimentateur, notamment en physique o` u les volumes de données peuvent être particulièrement important et les calculs répétitifs. L’utilisation d’un ordinateur permet dans ces cas de gagner beaucoup de temps. Les calculs numériques sont apparus bien avant l’apparition des ordinateurs. Au XVII siècle par exemple, John Napier (1550-1616) introduisit les logarithmes (Mirifici logarithmorum canoni descriptio, Edimbourg,1614). Les multiplications pouvait être transformée en addition, un calcul de racine carrée en division par deux. C’était déjà un gain de temps appréciable . . .De nombreux contemporains de Napier réalisèrent des prouesses en calculs numériques (calculs des orbites planétaires par exemple), et tout ceci ` a la main. Outre le temps requis par de tels calculs, les erreurs étaient monnaie courante. Aussi, l’idée d’automatiser ces calculs fastidieux apparut rapidement. Le baron Gaspard de Prony était chargé pendant le premier empire de mettre en place des tables pour le calcul des impˆ ots fonciers (dont les règles étaient déjà complexes à l’époque). Il eut l’idée de fractionner ce travail énorme en plusieurs partie. Des mathématiciens décomposèrent tous les calculs nécessaires en opérations élémentaires, des gestionnaires organisèrent le travail et collectèrent les résultats. Enfin, les calculs élémentaires furent assurés par des opérateurs pour lesquels la qualification requise était de savoir faire des additions. Cette approche est exactement celle d’un programmeur devant résoudre un problème ` a l’aide d’un ordinateur : sa tâche est d’écrire un algorithme, c’est à dire de décomposer un problème en une série de tâches élémentaires compréhensibles par un ordinateur (additions, comparaisons . . .), et de définir l’ordre d’exécution de ces tâches. Dans ce cas, les opérateurs humains sont remplacés par un calculateur beaucoup plus rapide. La mécanisation du calcul suivit, notamment grˆ ace à un anglais, Charles Babbage au début du XIX siècle. Il eut l’idée d’associer cette idée de décomposition d’un problème en tâches simples à une machine ` a calculer mécanique ressemblant à la machine de Pascal et à un système 4


5

de lecture de cartes perforées utilisées sur les métiers à tisser Jacquard. Sa machine à différences était l’analogue mécanique de nos calculateurs électroniques. Cependant, son idée n’aboutit pas, car il se heurta a de très importantes difficultés techniques. Mais l’idée était lancée . . . Cette idée fut concrétisée en 1890 par Herman Hollerith qui répondit à un appel d’offre du gouvernement américain pour traiter les données du recensement national. Le traitement des données du recensement précédent avait nécessité sept ans de travail. Hollerith proposa une machine qui lisait les données sur des cartes perforées, et le résultat fut obtenu en six semaines. Un gain de temps appréciable. Hollerith fonda alors une société, la Tabulating Machine Company qui devint en 1924 la International Business Machines, plus connue sous le sigle IBM . . . Au cours de la deuxième guerre mondiale, les efforts de recherche dans ce domaine se sont intensifiés : les armées avaient besoin de gros volumes de calcul (établissement de tables de portée de tir pour l’artillerie, décryptage des communications ennemies . . .). Les premiers calculateurs électroniques ont vu le jour très peu de temps après la guerre (l’ENIAC en 1946 aux Etats-Unis). L’invention du transistor en 1951 par John Bardeen, William Shockley et Walter Brattain a ensuite ouvert la voie aux calculateurs électroniques modernes. Leurs performances n’ont depuis jamais cessé de croˆıtre.


H. Klein

Repr´ esentation des nombres dans un ordinateur On parle souvent de précision du résultat au cours de calculs numériques. Cette précision n’est jamais une précision absolue. En effet, outre le fait qu’un algorithme puisse introduire une imprécision dans les résultats, les nombres eux-mêmes sont représentés avec une certaine précision dans la mémoire d’un ordinateur.

Repr´ esentation des nombres entiers Dans la mémoire d’un ordinateur, les nombres ont une représentation binaire, c’est à dire une suite de 0 et de 1. Aisi pour les nombres entiers on obtient la correspondance suivante entre représentation décimale (première ligne) et binaire (deuxième ligne) 0 1 2 3 4 5 6 7 8 9 10 ... 0 1 10 11 100 101 110 111 1000 1001 1010 . . . Une case contenant em mémoire un 0 ou un 1, est appelée un bit. Dans la majorité des ordinateurs, ces bits sont regroupés par 8, les octets ou bytes. Ces octets constituent alors un emplacement correspondant ` a une adresse dans la mémoire. Le tableau précédent s’écrirait en fait 0 1 2 3 4 5 ... 00000000 00000001 00000010 00000011 00000100 00000101 . . . Avec un octet on peut représenter les nombres entiers 0 . . . 28 − 1, c’st à dire de 0 à 255. Afin d’avoir une dynamique de représentation plus importante, on groupe des octets en mots, la taille de ces mots dépendant du nombre de bits qu’un ordinateur sait manipuler en une opération. Sur les ordinateurs personnels actuels, la taille de ces mots est de 32 bits soit 4 octets. Dans un mot de quatre octets, on peut coder les entiers allant de 0 à 232 − 1 soit de 0 à 4 294 967 296. Pour pouvoir représenter les nombres négatifs, on change de convention de codage : un nombre négatif est le complément ` a deux de sa valeur positive. Par exemple 00000001 = 1 et 11111110=-1 (au lieu de 254). Aisi de 00000000=0 à 01111111=127 on code des nombres positifs, et de 11111110=-1 ` a 11111111=-128 des nombres négatifs. Sur 32 bits, on peut coder des nombres entiers signés allant de −231 à 231 − 1 6


7

Repr´ esentation des nombres r´ eels Parler de√nombre réels est abusif, on ne peut en effet traiter que des nombres rationels, par exemple 2 = 1.4142213562 . . .. La restriction est encore plus grande puisque l’on ne peut représenter qu’un nombre fini de chiffres : 1/3 devient par exemple 0.33333333. Si l’on supprime la virgule, ce nombre peut être représenté√comme un entier, il suffit ensuite d’indiquer la position de la virgule par une puissance de 10. 2 = 1.4142213562 peut s’écrire 14142213562.10−9 ou encore 0, 14142213562.101 . L’avantage de la deuxième écriture est que l’expression des chiffres composant le nombre est comprise entre 0 et 1, chaque chiffre correspondant à une puissance négative de 10. On pourra donc représenter un nombre avec 3 papramètres : le signe, la position de la virgule (ou l’exposant) et les chiffres (ou la mantisse). La décomposition que nous avons fait en base 10 peut sans problème être transposée en base 2. Aisi un nombre réel codé sur 4 octets se répartit de la manière suivante – 1 bit de signe – 8 bits pour l’exposant – 23 bits pour la mantisse, comme le premier chiffre sera toujours zéro, on peut l’omettre et gagner ainsi un bit significatif Cette convention appelée notation en virgule flottante permet de coder les nombres réels de ±1, 175494.10−38 ` a ±3, 402823.1038 avec 7 chiffres significatifs. Il est aussi possible de représenter un nombre réels sur 64 bits (8 octets) en double précision avec 15 chiffres significatifs de ±2, 225074.10−308 ` a 1, 797693.10308 . Il convient cependant de se rappeler que cette représentation des nombres réels n’est pas exacte, c’st une approximation.

Cons´ equences La conséquence de ces représentations, est que les calculs se feront toujours avec une certaine précision intrinsèque (par exemple 7 chiffres significatifs). Imaginons que nous voulions réaliser le calcul suivant avec des nombres réels codés sur 4 octets si+1 = si + 10−9 avec s0 = 1, si l’on répète le calcul 109 fois, le résultat donnera 1... La raison en est que pour réaliser l’addition la machine dispose d’une précision de 10−7 , donc 1+10−9 = 1. Le même calcul conduit en double précision donnerait un résultat correct (la précision est alors de 15 chiffres). Il faut donc extrêmement prudent lorsque l’on réalise des opérations entre des nombres qui peuvent prendre des valeurs très différentes. Il est aussi prudent de se rappeler que lorsque l’on fait des calculs itératifs (ce qui est très fréquent en calcul numérique), les erreurs peuvent s’ajouter les unes aux autres pour finalement conduire à des résultats dont l’erreur n’est pas négligeable, voire ` a des résultats aberrants. Un deuxième cas classqiue d’erreur est le suivant : on compare un nombre réel if(a == 0.0) then .... Le résultat d’une telle comparaison est aléatoire à cause des arrondis ! Imaginons que a soit le résultat d’un long calcul, va-t-il valoir 0 ou 0, 1234567.10−35 ? Il est prudent de remplacer M´ ethodes Informatiques pour la physique

H. Klein


8

la comparaison précédente par if(|a| < ǫ), qui se traduirait par : si |a| est inférieur à la précision alors . . . Il faut cependant se poser la question de la précision adaptée à un calcul : si l’on est capable représenter un nombre avec une précision de 10−7 , cela ne représente pas une précision absolue mais relative. On pourra par exemple par exemple représenter x ≈ 1 avec cette précision, mais x ≈ 1020 sera représenté avec une précision de 10−7 × 1020 = 1013 . . . D’une mani` ere g´ en´ erale, lorsque l’on fait des calculs num´ eriques, il est utile de v´ erifier le calcul sur des cas o` u l’on connait la solution afin de v´ erifier le fonctionnement correct d’un programme.


H. Klein

Premi` ere partie

Tri et mod´ elisation de donn´ ees

9

Chapitre 1

Algorithmes de tri de donn´ ees Ce chapitre n’a a priori pas grand chose à voir avec le calcul numérique, c’est cependant une étape nécessaire et préalable ` a des manipulations plus sophistiquées ( par exemple des calculs statistiques) . C’est de plus un domaine o` u de très nombreux algorithmes ont été proposés, et nous allons en détailler quelques uns. Un peu d’algorithmique ne peut pas faire de mal . . . Pour effectuer un tri, nous allons placer les données dans un tableau, nous supposerons dans la suite que ce tableau comport N éléments. Le but sera de réduire au minimum le nombre d’opérations pour trier ces N éléments. Cette optimisation est souvent indispensable, lorsque le nombre de données est important, o` u lorsqu’on est amené a` répéter le tri un grand nombre de fois. Le nombre d’opérations nécessaires pour trier un tableau de N éléments, n − op est variable en fonction des algorithmes :

N log2 N ≤ n − op ≤ N 2

Dans ce chapitre nous allons aborder deux types de méthodes : – simples, mais peu performantes : tri à bulle ou par insertion – moins simples mais plus performantes : tri rapide ou shell 10

´ CHAPITRE 1. ALGORITHMES DE TRI DE DONNEES

1.1

11

Tri ` a bulle

Le principe du tri ` a bulle est des plus simples : supposons que l’on souhaite classer un tableau de nombres entiers par ordre croissant, on va alors parcourir le tableau contenant les N nombres. Le plus petit nombre sera rangé dans la preière case du tableau. Il ne reste plus qu’à trier N − 1 éléments pour lesquels on répète la même opération. Une illustration de l’algorithme 1 est donnée dans la figure 1.1 Algorithme 1 Tri d’un tableau par tri ` a bulle Requiert: un tableau tab[N] de N éléments Fournit: le tableau trié par ordre croissant pour i = 0 à N − 1 faire tampon ← i pour j = i + 1 ` a N − 1 faire si tab[j] ≤ tab[tampon] alors tampon ← j fin si fin pour ttampon ← tab[i] tab[i] ← tab[tampon] tab[tampon] ← ttampon fin pour renvoyer tab[]

Fig. 1.1 – Tri d’un tableau par la méthode du tri à bulle. Cet algorithme requiert en moyenne N 2 opérations pour trier un tableau de N éléments. Pour cette raison cette méthode de tri, très simple à implémenter n’est utilisée que pour de petits ensemble de données. M´ ethodes Informatiques pour la physique

H. Klein


1.2

12

Tri par insertion

Cette méthode de tri est bien connue des joueurs de cartes : on commence par trier les deux premières cartes l’une par rapport ` a l’autre, puis on trie la troisième par rapport aux deux premières, et ainsi de suite jusqu’à ce que toutes les cartes soient triées. Une illustration de l’algorithme 2 est donnée ` a la figure 1.2. Algorithme 2 Tri d’un tableau par insertion Requiert: un tableau tab[N] de N éléments Fournit: le tableau trié par ordre croissant pour i = 1 à N faire tampon ← tab[i] j ←i−1 tant que j ≥ 0 ET tab[j] > tampon faire tab[j + 1] ← tab[j] j ←j−1 fin tant que tab[j + 1] ← tampon fin pour renvoyer tab[]

Fig. 1.2 – Tri d’un tableau par la méthode du tri par insertion. Cet algorithme demande également N 2 opérations pour trier un tableau de N éléments, ce n’est donc pas a priori une méthode très performante pour de grands ensembles car la recherche de l’emplacement o` u ranger un élément dans l’ensemble trié nécessite éventuellement de le comparer ` a tous les autres éléments.Le grand intérêt de cette méthode est que dans le cas o` u l’on veut insérer un élément dans un ensemble déjà trié, il est simplement nécessaire de parcourir une fois l’ensemble déjà trié pour positionner le nouvel élément. C’est la raison pour laquelle cette méthode est particulièrement adaptée au tri incrémental. Considérons par exemple un système de réservation dans un hotel : tous les clients sont affichés par ordre alphabétique. Si la liste est mise ` a jour lorsque de nouveaux clients arrivent, la mise en place d’un nouveau M´ ethodes Informatiques pour la physique

H. Klein


13

client dans la liste ne nécessitera qu’un simple décalage des données.

1.3

Tri shell

Le tri shell est une variante plus performante du tri par insertion. Cet algorithme appliqué ` a un tableau rangé au hasard nécesite environ N 1.25 opérations. Cependant pour N > 50, le tri rapide est plus performant. L’idée de base est la suivante : soit un ensemble de nombres n1 . . . n1 6. Nous allons d’abord trier par insertion 8 groupes de 2 nombres (n1 , n9 ), (n2 , n10 ), . . . , (n8 , n16 ). Puis 4 groupes de 4 nombres (n1 , n5 , n9 ), (n2 , n6 , n10 , n14 ), . . . , (n4 , n8 , n12 , n16 ). Puis 2 groupes de 8 nombres (n1 , n3 , n5 , n7 , n9 , n11 , n13 , n15 ), (n2 , n4 , n6 , n8 , n10), n12 , n14 , n16 , et enfin un groupe de 16 nombres. A priori, seule la dernière étape est nécessaire. L’intérêt des étapes préliminaires est d’effectuer un filtrage qui va rapprocher tous les éléments de leurs places finales. Ainsi, à la dernière étape, seule quelques comparaisons sont nécessaires pour que chaque élément trouve sa place dans l’ensemble trié. La performance de l’algorithme dépend de l’incrément utilisé entre nombres. Dans notre cas, nous démarrons avec les couples (n1 , n9 ) . . ., puis (n1 , n5 . . .. L’incrément sera ici 8, 4, 2, 1, cette série n’est pas un bon choix. Il y a eu de nombreux travaux sur ces suites d’incréments, et il en ressort que le meilleur compromis est de la forme 3k − 1 , . . . , 40, 13, 4, 1 3 cette série peut-être générée par la récurrence suivante i1 = 1, ik+1 = 3ik + 1, k = 1, 2, . . . Il a été démontré qu’en utilisant cette série d’incréments, l’algorithme nécessite N 1.5 opérations dans le pire des cas. L’algorithme 3 est illustré par la figure 1.3.

Fig. 1.3 – Tri d’un tableau par la méthode du tri shell.


H. Klein


14

Algorithme 3 Tri d’un tableau par la méthode du tri shell Requiert: un tableau tab[N] de N éléments Fournit: le tableau trié par ordre croissant inc ← 1 r´ ep´ eter inc ← inc ∗ 3 inc ← inc + 1 tant que inc ≤ N r´ ep´ eter inc ← inc/3 pour i = inc ` a N faire tampon ← tab[i] j←i tant que tab[j − inc] > tampon faire tab[j] ← tab[j − inc] j ← j − inc si j < inc alors sortir fin si fin tant que tab[j] ← tampon fin pour tant que inc > 1 renvoyer tab[]


H. Klein


1.4

15

Tri rapide

Pour de grands ensemble de données , cet algorithme est l’algorithme le plus rapide sur la plupart des ordinateurs. Le principe de cet algorithme (4) est de partitionner le tableau à trier par rapport à une valeur pivot. Tout les éléments dont la valeur est supérieure au pivot seront placés après le pivot (haut de la figure 1.4, les éléments inférieurs avant. Cette opération est ensuite répétée sur les deux sous-ensembles de valeurs jusqu’à ce que le tableau soit complètement trié (bas de la figure 1.4). La valeur choisie pour le pivot a peu d’importance, on peut très bien choisir une valeur au hasard du tableau, mais pour plus de clarté, dans la figure 1.4, le médian des valeurs a été choisi comme pivot pour chaque opération. Algorithme 4 Tri d’un tableau par la méthode du tri rapide. Dans cet algorithme, le pivot est le dernier élément du tableau. Requiert: un tableau tab[N] de N éléments, initialement debut = 1, f in = N − 1 Fournit: le tableau trié par ordre croissant si debut < f in alors pivot ← tab[N − 1] i ← debut − 1 j ← f in r´ ep´ eter r´ ep´ eter i ←i+1 tant que tab[i] < pivot r´ ep´ eter j ←j−1 tant que tab[j] ≥ pivot tampon ← tab[i] tab[i] = tab[j] tab[j] ← tampon le pivot est maintenant dans tab[i] tant que i < j appel recursif 1 : debut, f in = i − 1 appel recursif 2 : debut = i + 1, f in fin si renvoyer tab[]


H. Klein


16

Fig. 1.4 – Tri d’un tableau par la méthode du tri rapide. En haut : on a choisi le pivot 28, premier tri du tableau. En bas : Visualisation des opérations de tris succesives.


H. Klein

Chapitre 2

Cacul de quelques grandeurs statistiques usuelles Ce chapitre ne se veut pas exhaustif, nous allons simplement aborder le calcul des grandeurs statistiques les plus couramment utilisées par les physiciens, comme la moyenne. Nous verrons notamment que, si certains calculs sont très simples à réaliser, ils peuvent néanmoins receler certains pièges. Quand un distribution de données a tendance à être groupé autour d’une valeur particulière, il peut alors être utile de le caractériser par des nombres reliés à ses moments, c.a.d. la somme de puissances entières de ses valeurs.

2.1

Moyenne

Le calcul d’une valeur moyenne x est certainement la manière la plus simple de caractériser une distribution de N données xi n 1 X xi x= N i=1 l’information obtenue est assez pauvre : on sait maintenant que les données sont réparties plus ou moins régulièrement autour de la valeur moyenne x. Le calcul est élémentaire ` a réaliser ` a l’aide d’un ordinateur, mais il peut poser des problèmes assez sérieux. Il convient bien sˆ ur de s’assurer que le type de représentation des nombres est compatibles avec les données utilisées. Lors du calcul de la moyenne, il faut réaliser la somme de toutes les données. Ceci peut rapidement conduire à un dépassement de capacité. La solution peut être de faire la division par n pour chaque élément avant de réaliser la somme. Mais, outre le rallongement du temps de calcul, cela peut conduire à des erreurs d’arrondis, du fait du nombre d’opérations réalisées. Il ne faudra donc pas oublier de tester les d´ epassements de capacit´ e lors de ce calcul. L’utilisation d’une valeur moyenne pour estimer une distribution de données n’est pas forcément judicieux : imaginons une distribution très dispersée de valeurs, dans ce cas la valeur moyenne ne sera pas représentative. Il peut alors être plus judicieux d’utiliser le médian (voir 2.3) 17

CHAPITRE 2. CACUL DE QUELQUES GRANDEURS STATISTIQUES USUELLES

2.2

18

variance, ´ ecart-type

Lorsque l’on a caractérisé la valeur centrale d’une distribution, on caractérise souvent sa ”largeur” ou sa dispersion autour de cette valeur. La mesure la plus courante est appelée variance. V ar(x1 . . . xN ) =

N 1 X (xi − x)2 N − 1 i=1

ou écart-type lorsque l’on considère sa racine carrée σ(x1 . . . xN ) =

q

V ar(x1 . . . xN )

La variance est une estimation de l’écart moyen (au carré) entre les valeurs xi et la moyenne x

2.3

M´ edian

Le médian d’une distribution de probabilité p(x) est la valeur xmed pour laquelle les valeurs supérieures ou inférieures de x sont équiprobables Z

xmed

−∞

1 p(x)dx = = 2

Z

∞

p(x)dx

xmed

Il est assez simple de déterminer le médian d’une distribution x1 . . . xN , en déterminant la valeur xi qui a un nombre équivalent de valeurs au dessus et au dessous d’elle, si le nombre d’éléments de la distribution est pair. Dans le cas d’un nombre d’éléments impair, on calculera la valeur moyenne des deux éléments centraux. xmed = x(N +1)/2 , N pair 1 xmed = (xN/2 + x(N/2)+1 ), N impair 2 Pour trouver un median, on peut classer toutes les valeurs d’une distribution par ordre croissant ou décroissant, et d’utiliser ensuite l’équation ci-dessus. Dans le cas o` u une distribution est dispersée, le médian peut être une meilleure estimation que la moyenne.


H. Klein

Chapitre 3

R´ egression lin´ eaire par la m´ ethode des moindres carr´ es 3.1

D´ efinition du probl` eme

Nous allons maintenant nous intéresser un problème qui se pose souvent aux expérimentateurs. A l’issue d’une expérience, un expérimentateur a souvent enregistré des donn´ ees relatives par exemple à l’observation d’un phénomène physique. Ces données sont toujours entachées de bruit (le plus souvent dˆ u aux instruments utilisés). L’expérimentateur souhaite néanmoins comparer ses données à un mod` ele. Nous ne nous intéresserons ici qu’au cas o` u le modèle est connu. L’expérimentateur souhaite donc uniquement s’assurer que le phénomène observé suit bien la loi o` u le modèle qu’il connaˆıt. Il souhaitera également déterminer les valeurs numériques des paramètres intervenant dans ce modèle, et qui ne sont pas connues a ` priori.

3.2

Impl´ ementation

L’expérimentateur est en possession de n données (xi , yi ). La loi connue s’exprime comme y = f (x, p1 , p2 . . . pm ). Cette loi dépend des paramètres p qui ne sont pas nécessairement connus. Pour vérifier si les données se superposent correctement à la loi f , il suffit de calculer les écarts yi − f (x, pj ) en faisant varier la valeur de pj jusqu’à ce que l’écart soit minimal (o` u nul si l’accord est parfait). Il faut bien évidemment réaliser ce calcul pour tous les points expérimentaux. Les statistiques nous offrent une manière élégante de calculer cela. On introduit pour cela une variable dépendant des paramètres pj que l’on appelle χ2 (prononcer”ki-2”). Cette variable est définie comme χ2 (pj ) =

n X i=1

(f (xi , {pj }) − yi )2

on va ensuite chercher l’ensemble de paramètres pj qui va minimiser χ2 . 19

´ ´ ´ ´20 CHAPITRE 3. REGRESSION LINEAIRE PAR LA METHODE DES MOINDRES CARRES

3.3

Exemple d’application : ajustement d’une fonction lin´ eaire y = a.x + b

Dans ce cas précis, χ2 s’exprimera χ2 =

n X i=1

(axi + b − yi )2

on cherche à minimiser cette grandeur, on va donc chercher a` satisfaire les conditions n ∂χ2 X = 2xi (axi + b − yi ) = 0 ∂a i=1 n ∂χ2 X = 2(axi + b − yi ) = 0 ∂b i=1

on en déduit que

a=

n

Pn Pn i=1 yi i=1 xi i=1 xi yi − Pn Pn 2 n i=1 xi − ( i=1 xi )2

Pn

2 Pn y − Pn x y Pn x i i=1 i i=1 i i i=1 xi Pi=1 P b= n ni=1 x2i − ( ni=1 xi )2 deux paramètres a et b qui minimisent le χ2 ,

Pn

On trouve ainsi les c’est à dire les paramètres qui fournissent le meilleur accord entre les données expérimentales (xi , yi ) et la fonction y. La valeur numérique de χ2 va nous renseigner sur cet accord. Une valeur nulle indiquera un accord parfait entre les données et le modèle, plus la valeur sera grande, moins bon sera l’accord. C’est ensuite à l’utilisateur de cette méthode de déterminer quelles valeurs peuvent être satisfaisantes ...

3.3.1

prise en compte de l’incertitude sur les mesures

Dans l’introduction, nous avons parlé de données, et de bruit. Un expérimentateur rigoureux va toujours tenter de caractériser ce bruit qui va se traduire par une erreur sur chaque mesure, que l’on notera σi . Dans la suite, nous considérerons que cette erreur ne concerne que les points yi , ce qui signifie que l’on connaˆıt exactement les valeurs xi . Nous avons donc une incertitude σi pour chaque couple de mesures (xi , yi ). Nous exprimerons alors χ2 comme X axi + b − yi 2 2 χ = σi on va de nouveau chercher ` a le minimiser en posant

et

∂χ2 X 2xi (axi + b − yi ) =0 = ∂a σi2


∂χ2 X 2(axi + b − yi ) =0 = ∂b σi2 H. Klein


alors, en posant Sxy =

X xi y i

σi2

, S = S11

∆ = SSxx − (Sx )2 ) on obtient

SSxy − Sx Sy ∆ Sxx Sy − Sx Sxy b= ∆ ce sont à nouveau les paramètres qui fournissent le meilleur accord entre données et modèle en tenant compte des erreurs de mesures. Et l’on peut maintenant calculer l’erreur sur ces paramètres S σa = ∆ Sxx σb = ∆ Les statistiques nous permettent également de calculer un facteur de confiance sur l’opération qui vient d’être réalisée : le facteur de confiance pondéré (wheighted en anglais ? ? ?) a=

v u u Rw = t P

χ2 f (xi ,{pj }) σi2

par exemple Rw = 0.1 signifie que f est éloignée des données de 10 % en unités de σ.

Fig. 3.1 – Exemple d’application de la régression linéaire à un ensemble de données. Chaque couple de mesures (xi , yi ) est accompagné d’une incertitude sur yi , σi , représentée comme une barre verticale. Dans l’exemple de la figure 3.1, est représenté un ensemble de données qui suivent visiblement un modèle linéaire. L’application de la régression linéaire non pondérée donne les résultats suivants : – a = 1.25 M´ ethodes Informatiques pour la physique

H. Klein


– b = 0.5 le tracé de cette fonction sur la même figure, montre que l’accord entre données expérimentales et régression est excellent. L’application de la régression pondérée donne les résultats suivants : – a = 1.25, σa = 0.19 – b = 0.5, σb = 0.42 – Rw est quant ` a lui très proche de zéro. on voit donc que si la régression est tout a` fait pertinente sur cet ensemble de données, l’incertitude sur le paramètre b est loin d’être négligeable.


H. Klein

Deuxi` eme partie

M´ ethodes num´ eriques

23

Chapitre 4

Recherche des racines d’une fonction Nous nous intéressons ici ` a un des problèmes les plus basiques : résoudre numériquement une équation du type f (x) = 0 Dans le cadre de ce document, nous n’allons nous intéresser qu’au cas o` u il y a une seule variable indépendante x, c’est ` a dire au problème ` a une dimension. Excepté le cas des problèmes linéaires, la recherche des racines de la fonction f (x) va se faire par itérations. Il faut pour cela partir d’une solution approchée, même grossière, et bˆ atir un algorithme qui va “améliorer” la solution jusqu’à satisfaire un critère de convergence. Nous allons nous intéresser à deux types d’algorithmes : – méthode dichotomique – méthode utilisant la dérivée Mais avant d’aborder ces algorithmes, nous allons étudier des algorithmes simples d’encadrement permettant de définir des solutions approchées de f (x) = 0.

4.1

Encadrement

On dira qu’une racine d’une fonction f (x) est encadrée par l’intervalle (a, b) si f (a) et f (b) ont des signes opposés. Si la fonction est continue, alors il existe au moins une racine de f dans cet intervalle (thèorème de la valeur intermédiaire). Si la fonction n’est pas continue, on peut à la place d’une racine, trouver une discontinuité qui traverse le zéro dans l’intervalle (a, b). D’un point de vue numérique, cela va être équivalent à une racine car le comportement va être le même, c.a.d. que l’on va traverser le zéro entre deux nombres flottants de la représentation en précision finie de la machine. Seules les fonctions du type f (x) =

1 x−c

vont présenter des singularités mais pas de racines. Dans ce cas, certains algorithmes pourront converger vers c, mais dans ce cas il n’y a (heureusement) pas de confusion possible avec une racine, puisqu’une évaluation de |f (x)| va donner un résultat très grand plutˆ ot que proche de zéro. Il est donc important d’avoir une idée de la variation de f (x) dans l’intervalle étudié pour 24

CHAPITRE 4. RECHERCHE DES RACINES D’UNE FONCTION

25

choisir un algorithme adapté, d’o` u la nécessité d’avoir des algorithmes simples permettant de déterminer les intervalles contenant une racine pour une fonction quelconque. L’algorithme 5 permet, partant d’un intervalle de départ (x1 , x2 ), de l’élargir jusqu’à ce qu’il contienne une racine, les valeurs renvoyées sont alors le nouvel intervalle (x1 , x2 ) contenant une racine. Ce type de procédure est une bonne manière de commencer une recherche de racine, et a de bonnes chances de fonctionner si la fonction étudiée a des signes opposés à la limite x → ±∞.

Algorithme 5 Encadrement d’une fonction Requiert: un intervalle (x1 , x2 ), un nombre d’essais N Requiert: f (x) définie dans l’intervalle (x1 , x2 ) étudié, x1 6= x2 Fournit: l’intervalle x1 , x2 contenant une racine f1 = f (x1 ) f2 = f (x2 ) α ← 1.6 (par exemple) pour i = 1 à N faire si f1 < f2 alors Renvoyer x1 et x2 fin si si |f1 | < |f2 | alors f1 = f (x1 + α(x1 − x2 )) sinon f2 = f (x2 + α(x2 − x1 ) fin si fin pour

L’algorithme 6 permet lui, ` a l’inverse, de réduire l’intervalle de départ. On va simplement subdiviser l’intervalle de départ (x1 , x2 ) en n intervalles, et renvoyer les sous-intervalles contenant une racine sous forme de deux tableaux. A l’aide de ces algorithmes, nous sommes maintenant capables de trouver des intervalles contenant les racines d’une fonction, c’est un pré-requis pour estimer plus précisemment la valeur d’une racine particulière d’une fonction. M´ ethodes Informatiques pour la physique

H. Klein


26

Algorithme 6 Réduction de l’encadrement d’une fonction Requiert: L’intervalle (x1 , x2 ) fourni contient une racine de f Requiert: nombre maximum de racines recherchées N Fournit: nr le nombre de racines trouvés, x1 [] et x2 [] deux tableaux contenant les bornes des intervales contenant les racines nr ← 0 dx ← (x2 − x1 )/n x ← x1 f p = f (x) pour i = 0 à N faire x ← x + dx f c = f (x) si f p ∗ f c ≤ 0 alors nr ← nr + 1 x1 [nr] ← x − dx x2 [nr] ← x fin si fp = fc fin pour renvoyer nr

4.2

Dichotomie

Le principe est extrêmement simple, il est schématisé sur la figure 4.1. On choisit un intervalle (a, b) encadrant une racine de la fonction f étudiée. On divise cet intervalle en 2 par un point c. On évalue f (a).f (c), si f (a).f (c) < 0, alors le nouvel intervalle devient (a, c), sinon (c, b). Et l’on recommence ainsi jusqu’à ce que l’intervalle satisfasse le critère de convergence. Cette méthode a de bonnes chances de fonctionner si la fonction étudiée a des signes opposés à la limite x → ±∞. L’algorithme 4.2 est un algorithme simple permettant de déterminer une valeur approchée de la racine d’une fonction dans le cas o` u f (x) traverse le zéro dans l’intervalle (x1 , x2 ) étudié.

4.2.1

Convergence de la m´ ethode

Dans cette algorithme, après chaque itération, l’intervalle (x1 , x2 ) a diminué d’un facteur 2. Après n itérations, la racine est contenu dans un intervalle de taille ǫn , avec ǫn ǫn+1 = 2 Donc, pour parvenir ` a la précision ǫ désirée en partant de ǫ0 = x2 − x1 il faudra ǫ0 n = log2 ǫ itérations. Le cas étudié ici o` u l’on divise l’intervalle par 2 à chaque itération présente une convergence lin´ eaire On voit en effet que pour ǫn+1 = α × ǫm n avec α < 1 M´ ethodes Informatiques pour la physique

H. Klein


27

Fig. 4.1 – principe de la méthode de recherche de racines par dichotomie. L’intervalle initial ( 1 , 2 ) est divisé en 2 par 3 , puis l’intervalle devient ( 1 , 3) Algorithme 7 Recherche d’une racine par dichotomie Requiert: l’intervalle d’étude (x1 , x2 ) et la précision désirée ǫ Fournit: l’intervalle encadrant la racine ` a la précision près f 1 ← f (x1 ) f 2 ← f (x2 ) r´ ep´ eter f temp ← f ((x1 + x2 )/2) si f 1 ∗ f temp < 0 alors x2 ← (x1 + x2 )/2 sinon x1 ← (x1 + x2 )/2 fin si tant que |x2 − x1 | > ǫ renvoyer x1 et x2

– la convergence est linéaire si m = 1 – la convergence est superlinéaire si m > 1 Pour obtenir des résultats cohérents avec un algorithme de ce type, il faut faire attention au critère de convergence ǫ. Il est en effet nécessaire de se rappeler qu’un ordinateur utilise un nombre fini de bits pour représenter un nombre en virgule flottante. Si, analytiquement la fonction f passe par zéro, il est possible que la valeur calculée ne vale jamais zéro. Si la racine recherchée est proche de 1, demander une précision absolue de 10−8 est raisonnable, mais ne l’est plus si la racine est proche de 1020 par exemple. On pourrait fournir une précision relative (fractionnaire), mais dans ce cas le calcul devient problématique si la racine est proche de zéro. . .D’une manière générale, il est préférable de fournir une tolérance absolue pour que les itérations se poursuivent jusqu’à ce que l’intervalle de recherche deviennent plus petit que cette tolérance en unités absolues. La valeur ǫm (|x1 | + |x2 |)/2 o` u ǫm est la précision de représentation de la machine, et (x1 , x2 ) l’intervalle initial est un bon choix. Il faut cependant toujours se poser M´ ethodes Informatiques pour la physique

H. Klein


28

la question de ce que représente une tolérance raisonnable lorsque la racine est proche de zéro.

4.3

M´ ethode de Newton - Raphson

Cette méthode est la plus utilisée pour la recherche de racines dans les problèmes à une dimension. Elle requiert cependant l’évaluation de f (x) et de f ′ (x). Le principe est le suivant : on prend la tangente de f (x) au point xi (solution approchée), et on la prolonge jusqu’à l’axe des abscisses (fig. 4.2). On utilise alors cette nouvelle abscisse comme point xi+1 , et l’on recommence tant que xi+1 − xi ne satisfait pas le critère de convergence.

Fig. 4.2 – Principe de la méthode de Newton-Raphson : on part de la solution approchée 1 , l’extrapolation de la tangente ` a la courbe pour cette abscisse fournit la nouvelle solution approchée 2. Cela peut se traduire algébriquement grˆ ace aux développements de Taylor δ2 ′′ f (x) . . . 2 Si on néglige les termes d’ordre supérieurs a` 1, écrire f (x+δ) = 0 revient à écrire f (x)+δf ′ (x) ≈ 0 d’o` u f (x) δ=− ′ f (x) Lorsque l’on se trouve loin d’une racine de la fonction étudiée, les termes du développement d’ordre supérieur ` a 1 ne peuvent pas être négligés, aussi la méthode peut donner des résultats aberrants. Par exemple si la solution approchée de départ est loin d’une racine, l’intervalle iniial de recherche peut être proche d’un extrémum local (fig.4.3), ce qui signifie que f ′ (x) → 0. Cela sera fatal dans le cadre de cette méthode. . . L’algorithme 4.3, présente une implémentation de cette méthode. f (x + δ) ≈ f (x) + δf ′ (x) +

4.3.1

Convergence de la m´ ethode

Pour ǫ petit f (x + ǫ) ≈ f (x) + ǫf ′ (x) + ǫ2 M´ ethodes Informatiques pour la physique

f ′′ (x) + ... 2 H. Klein


29

Fig. 4.3 – Cas particulier o` u l’on trouve un extrémum local. Dans ce cas la méthode de NewtonRaphson diverge Algorithme 8 Recherche de racines par la méthode de Newton-Raphson Requiert: valeur de départ x, précision désirée ǫ Fournit: valeur x ` a la précision désirée si f ′ (x) = 0 alors Sortie du programme fin si r´ ep´ eter xn ← x x ← x − f (x)/f ′ (x) tant que |x − xn | ≥ ǫ ET f ′ (x) = 0 renvoyer x f ′ (x + ǫ) ≈ f ′ (x) + ǫf ′′ (x) + . . . d’après la formule de Newton-Paphson xi+1 = xi −

f (xi ) f (xi ) ⇔ ǫi+1 = ǫi − ′ ′ f (xi ) f (xi )

Quand une solution approchée xi diffère de la solution vraie α par une quantité ǫi = α − xi , on a f ′′ (xi ) f (α) = f (xi + ǫi ) = f (xi ) + ǫi f ′ (xi ) + ǫ2i + ... 2 f (α) = 0 ⇔ f (xi + ǫi ) = f (xi ) + (α − xi )f ′ (xi ) + ǫ2i donc

f ′′ (xi ) + ... = 0 2

′′ f (xi ) 2 f (xi ) − x ≃0 +α + ǫ i i f ′ (xi ) 2f ′ (xi )

| M´ ethodes Informatiques pour la physique

{z

=−xi+1

}

H. Klein


d’o` u ǫi+1 = α − xi+1 ≃ −ǫ2i

30

f ′′ (xi ) 2f ′ (xi )

on voit donc que la convergence est quadratique, la méthode est plus efficace que la dichotomie o` u la convergence n’est que linéaire. Ceci suppose bien sˆ ur que l’on sache évaluer f ′ (x) en tout point, et que f ′ (α) 6= 0. Si la forme analytique de f ′ (x) n’est pas connue, il faut calculer numériquement cette quantité par f (x + dx) − f (x) f ′ (x) ≈ dx on√doit donc évaluer f (x) 2 fois ` a chaque itération, l’ordre de convergence va donc passer de 2 à 2. De plus, si l’intervalle dx est trop petit, il y a un risque que les erreurs d’arrondi fassent diverger le calcul. A l’opposé, si l’intervalle dx est grand, la convergence ne sera que linéaire. En conclusion on peut donc dire que la méthode de Newton-Raphson n’est pas une méthode intéressante si l’on ne peut déterminer de forme analytique de f ′ (x).


H. Klein

Chapitre 5

Minimisation ou maximisation d’une fonction Le problème peut-être posé simplement : considérons une fonction f dépendant d’une ou plusieurs variables indépendantes. Nous voulons trouver les valeurs de ces variables pour lesquelles f prend une valeur maximale ou minimale. La recherche d’un maximum ou d’un minimum peuvent être réalisée par un programme identique en travaillant simplement sur f ou 1/f . Un extrémum d’une fonction peut être soit local (minimum dans un intervalle donné)soit global (minimum réel) comme on le voit sur la figure 5.1

Fig. 5.1 – Les extrémums d’une fonction peuvent être locaux ou globaux C’est une des raisons pour lesquelles la recherche d’un extrémum est souvent une tâche délicate. Ce chapitre va être extrêmement succint, puisque nous n’aborderons qu’une méthode de minimisation d’une fonction présentant une variable indépendante (probléme 1D).

5.1

Recherche par la section d’or (golden section search)

Nous avons déjà vu dans le chapitre 4 que pour obtenir une valeur approchée d’une racine d’une fonction il était suffisant de l’encadrer dans un intervalle (a, b). Ceci est insuffisant pour 31

CHAPITRE 5. MINIMISATION OU MAXIMISATION D’UNE FONCTION

32

caractériser un minimum. Nous aurons besoins dans ce cas de trois points (a, b, c). Un intervalle tel que a < b < c et f (b) < f (a), f (b) < f (c) caractérise un minimum dans l’intervalle (a, c). En d’autres termes pour encadrer un minimum nous avons besoin d’un triplet de point tel que le point central présente une valeur de f inférieure à celles des bornes de l’intervalle.

5.1.1

Principe

Le principe de cette recherche est résumé sur la figure 5.2. Le minimum est encadré à l’origine par ( 1 , 3 , 2 ). La fonction est évaluée en 4 , qui remplace 2 . Ensuite en 5 qui remplace 1 puis en 6 qui remplace 4 . La règle est de garder un point central pour lequel la valeur de f est inférieure à celles des bornes. Après cette série d’itérations, le minimum est alors encadré par ( 5 , 3 , 6 ). Le principe est analogue ` a celui utilisé lors d’une recherche de racine par dichotomie (cf. chapitre 4.2)

Fig. 5.2 – Encadrements succesif d’un minimum d’une fonction f . Le minimum est encadré à l’origine par ( 1 , 3 , 2 ). Après 3 évaluations de f le minimum est finalement encadré par ( 5 , 3 , 6 ). Le seul point délicat, est de définir une méthode pour choisir un nouveau point d’encadrement dans l’intervalle (a, b, c) initial. Supposons que b soit une fraction ω du segment [ac], alors c−b b−a =ω⇔ =1−ω c−a c−a Si l’on choisit un nouveau point x eloigné d’une fraction z par rapport à b x−b =z c−a le nouveau segment d’encadrement aura alors une longueur ω + z par rapport au segment initial, ou de longueur 1 − ω. Si nous voulons minimiser le pire cas (attitude plutˆ ot saine...) nous allons choisir z pour que ces deux longueurs soient équivalentes. 1 − ω = ω + z ⇔ z = 1 − 2ω M´ ethodes Informatiques pour la physique

H. Klein


33

Fig. 5.3 – Comment choisir un nouveau point√pour réduire l’encadrement du minimum. La position optimale est ω = (b − a)/(c − a) = (3 − 5)/2. Le nouveau point x doit être choisi pour respecter cette proportion par rapport ` a (a, b) ou (b, c). On montre alors simplement que le nouveau point x est le symétrique de b dans son intervalle d’origine, c.a.d. que |b − a| = |x − c|. Cela signifie que le point x se trouve dans le plus grand des segments [ab], [bc]. Reste maintenant ` a définir la position du point x à l’intérieur de ce segment. La valeur ω provient d’une étape précédente de calcul, et si l’on suppose qu’elle est optimale, alors z doit être choisie de la même manière. Cette similarité d’échelle implique que x doit être situé à la même fraction du segment [bc] (si [bc] est le segment le plus long) que b l’était du segment [ac]. Cela conduit ` a la relation z =ω 1−ω en combinant cette équation avec la définition de z on arrive à l’équation quadratique √ 3− 5 ω − 3ω = 0 ⇔ ω = ≈ 0.38197 2 2

En d’autres termes, ce résultat signifie que l’intervalle d’encadrement (a, b, c) a son point central b est situé a une distance telle que b−a = ω ≈ 0.38197 c−a ou c−b = 1 − ω ≈ 0.61803 c−a ces fractions sont celles de la moyenne d’or ou de la section d’or supposées avoir certaines propriétés esthétiques par les pythagoriciens du monde antique. C’est la raison pour laquelle cette méthode optimale de détermination du minimum d’une fonction est appelée recherche par section d’or (golden search section). L’algorithme 9 décrit une impémentation de cette méthode. M´ ethodes Informatiques pour la physique

H. Klein


34

Algorithme 9 Recherche d’un minimum par la methode golden section search Requiert: un intervalle de départ (a, b, c), une tolérance ǫ Fournit: la meilleure estimation du minimum R ← 0.38197 C ←1−R x0 ← a x3 ← c si (c − b) > (b − a) alors x1 ← b x2 = (c − b) ∗ R + b sinon x2 ← b x1 ← a + (b − a) ∗ R fin si f1 ← f (x1 ) f2 ← f (x2 ) tant que |x3 − x0 | > ǫ ∗ (|x1 | + |x2 |) faire si f2 < f1 alors x0 ← x1 ; x1 ← x2 ; x2 ← C ∗ x1 + R ∗ x3 f1 ← f2 ; f2 ← f (x2 ) sinon x3 ← x2 ; x2 ← x1 ; x1 ← C ∗ x2 + R ∗ x0 f2 ← f1 ; f1 ← f (x1 ) fin si fin tant que si f1 < f2 alors renvoyer x1 sinon renvoyer x2 fin si


H. Klein


5.1.2

35

Convergence

Cette méthode garantit que chaque nouvelle évaluation de la fonction va encadrer le minimum dans un intervalle qui vaut 0.61803 fois le précédent, en d’autres termes ǫi+1 = 0.61803ǫi Cela signifie donc que la convergence est linéaire, mais moins rapide que dans le cas de la méthode de recherche de racines par dichotomie o` u le nouvel intervalle valait la moitié du précédent. La méthode abordée dans ce chapitre n’est bien évidemment pas la seule. Il existe des méthodes pouvant être plus performantes, ou mieux adaptées à des cas particuliers, et également un très vaste champ d’investigation concernant les problèmes à N dimensions que nous n’avons pas abordés. Je ne saurais trop vous conseiller la lecture du chapitre 10 des Numerical Recipes 1 pour vous cultiver sur ce sujet. . .

1

Numerical Recipes in C, W.H. Press et al., Cambridge University Press,2002. Le même livre existe pour d’autres language de programmation : Fortran et C++. M´ ethodes Informatiques pour la physique

H. Klein

Chapitre 6

Interpolation L’interpolation est un concept fondamental en calcul numérique : on évalue une fonction f en x0 , x1 . . . xn . Que faire si l’on veut connaˆıtre la valeur de f (x) si x0 < x < x1 ? Il n’y a pas d’évaluation de f pour cette valeur, on a donc recours a l’interpolation qui consiste en une estimation d’une valeur inconnue s’appuyant sur des valeurs connues qui l’encadrent. Ce problème se pose fréquemment pour l’expérimentateur scientifique qui réalise des mesures d’une grandeur en un certain nombre de points, et qui a parfois besoin de connaˆıtre une approximation de cette grandeur en dehors des valeurs mesurées. Une autre utilisation courante de l’interpolation est le lissage de courbe. Cela rend le résultat graphiquement plus satisfaisant, mais la courbe obtenue n’est pas plus précise que la série f (x0 ), f (x1 ) . . . f (xn ), mais peut permettre d’avoir une estimation d’une valeur x non comprise dans la série. Ceci à la condition que l’on soit capable d’estimer l’erreur commise lors du processus d’interpolation. De plus, l’interpolation est à la base de nombreuses méthodes numériques d’intégration de fonctions. Ces techniques sont basées sur les polynômes de Lagrange. On peut en effet montrer que pour une fonction f définie sur un intervalle [a, b] pour n + 1 valeurs distinctes, il existe un polynôme d’ordre n tel que Pn (x) = f (xi ) pour i = 0 . . . n. Nous détaillerons dans ce chapitre les calculs pour des polynômes d’ordre 1 et 2, puis nous généraliserons à l’ordre n.

6.1

Interpolation lin´ eaire

Il est bien connu qu’il n’y a qu’une droite passant par deux points (xi , fi ) et (xi+1 , fi+1 ), comme cela est schématisé sur la figure 6.1. Cette droite peut-être décrite par un polynôme Pi (x) de premier degré en x : fi+1 − fi (x − xi ) Pi (x) = fi + xi+1 − xi si l’on suppose que f a une dérivée seconde, l’écart entre f et Pi (et donc l’erreur commise lors d’une interpolation peut-être estimée par max |f (x) − Pi (x)| ≤

(xi+1 − xi )2 max |f ′′ (x)| 2 36

CHAPITRE 6. INTERPOLATION

37

pour xi ≤ xi+1 . Si l’écart entre deux évaluations successives de la fonction (le pas d’échantillongage) vaut xi+1 − xi = h alors l’erreur s’exprime max |f (x) − Pi (x)| ≤

h2 max |f ′′ (x)| 2

Fig. 6.1 – La fonction f est définie aux points xi , entre ces points, il est possible de l’interpoler par un polynôme d’ordre 1 passant par deux points consécutifs. exemple : soit la fonction f (x) = sin(x) évaluée (tabulée) en n + 1 points dans l’intervalle [0; π2 ] avec une précision désirée de 10−m . Dans ce cas le pas d’échantillonage h = π/2n, et l’on souhaite 2 π d2 sin(x) 1 −m max |f (x) − P (x)| ≤ 10 ⇔ | ≤ 10−m × max | 2n 2 dx2 d2 sin(x) = − sin(x) dx2 donc max |

d2 sin(x) |=1⇔ dx2

π 2n

2

×

1 ≤ 10−m 2

d’o` u

π2 8n2 π2 m −m m 2 .10 ≤ 10 ⇔ ≥ 10 ⇔ n ≥ 8n2 π2 8 ce qui conduit finalement ` a π n ≥ √ .10m/2 2 2

Concrètement, cela signifie que si l’on souhaite pouvoir interpoler un point de la fonction f avec une précision de 10−m , il va être nécessaire d’évaluer la fonction f en n points. Le tableau 6.1 résume ceci pour quelques valeurs de m. Ceci vous montre qu’augmenter la précision peut couter extrêmement cher en temps de calcul. Il pourrait être astucieux de faire varier le pas d’échantillonage h de la fonction f , de fa¸con à suivre le rayon de courbure local de la fonction : là o` u il est grand choisir h n’a pas besoin d’être grand, et inversement. On pourrait construire un algorithme réalisant cette opération en M´ ethodes Informatiques pour la physique

H. Klein


38

précision 10−3 10−6 10−8

n ∼ 35 ∼ 1000 ∼ 11000

Tab. 6.1 – Nombre d’évaluations n nécessaires pour parvenir à une précision donnée dans le cas d’une interpolation linéaire. évaluant la dérivée seconde de la fonction en chaque point. Néanmoins, le calcul d’une dérivée seconde n’est pas toujours facile ` a réaliser, il vaut donc mieux utiliser une méthode plus précise, c.a.d. utiliser un polynôme d’ordre plus élevé.

6.2

Interpolation quadratique

Le principe est identique ` a l’interpolation linéaire : il existe un polynôme unique d’ordre 2 Pi+1 qui passe par 3 points décrivant notre fonction f (voir figure 6.2) (x − xi+2 )(x − xi ) (x − xi+1 )(x − xi+2 ) + fi+1 (xi − xi+1 )(xi − xi+2 ) (xi+1 − xi+2 )(xi+1 − xi ) (x − xi )(x − xi+1 ) +fi+2 (xi+2 − xi )(xi+2 − xi+1 )

Pi+1 (x) = fi

Fig. 6.2 – La fonction f est définie aux points xi , entre ces points, il est possible de l’interpoler par un polynôme d’ordre 2 passant par trois points consécutifs. Cette expression est simple lorsque le pas d’échantillonage est constant. Si ce pas vaut h alors x = xi+1 + u.h. Si u = −1 x = xi , si u = 0 x = xi+1 et si u = 1 x = xi+2 . L’équation précédente peut alors s’écrire plus simplement Pi+1 (u) = fi+1 + M´ ethodes Informatiques pour la physique

fi+2 − 2fi+1 + fi 2 fi+2 − fi u+ u 2 2 H. Klein


39

L’erreur est maintenant majorée par le terme d’ordre 3 du développement en série de Taylor de f (x) au point xi+1 : h3 (3) max |f (x) − Pi+1 (x)| ≤ |f (x)| max 6 Exemple : Si l’on reprend l’exemple précédent de f (x) = sin(x) évaluée en n + 1 points, avec une précision de 10−m , on obtient l’inégalité π n≥ √ .10m/3 236 Le tableau 6.2 présente le nombre d’évaluations nécessaires de la fonction pour une précision donnée. Il est clair que cette méthode est plus performante dans le cas de la fonction étudiée que l’interpolation d’ordre 1. précision 10−3 10−6 10−8

n ∼9 ∼ 90 ∼ 400

Tab. 6.2 – Nombre d’évaluations n nécessaires pour parvenir à une précision donnée dans le cas d’une interpolation quadratique.

6.3

G´ en´ eralisation ` a l’ordre n

Un polynôme de Lagrange d’ordre n s’exprime Pn (x) =

n X i=0

fi

(x − x0 ) . . . (x − xi−1 )(x − xi+1 ) . . . (x − xn ) (xi − x0 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )

ou plus simplement Pn (x) =

n X

Li (x)f (xi )

i=0

avec Li (x) =

n Y (x − xj ) j=0

j6=i

(xi − xj )

L’erreur est dans ce cas majorée par max |f (x) − Pn (x)| ≤ max |f (n+1) (x)|.

(xn − x0 )n+1 (n + 1)!

Ceci suppose que la dérivée d’ordre (n + 1) de f existe et est continue. L’algorithme 10 permet de calculer la valeur d’un polynôme de Lagrange d’ordre n en un point, connaissant la fonction f en n + 1 points. M´ ethodes Informatiques pour la physique

H. Klein


40

Exemple : si l’on reprend l’exemple de la fonction f (x) = sin(x) définie en n + 1 points sur l’intervalle [0; π/2], et que l’on veuille interpoler une valeur avec une précision de 10−m alors −m

10

6.3.1

≤

π 2n

n+1

.

1 (n + 1)!

Du bon usage de l’interpolation

Il convient de faire attention ` a l’ordre du polynôme utilisé. En effet si l’on considère la fonction f (x) = sin(10x), l’erreur |f (n+1) |max ≈ 10n peut augmenter très rapidement avec l’ordre du polynôme. On peut établir la règle de prudence suivante : – si f (x) varie lentement, elle pourra être approximée correctement par un polynôme de degré élevé, c.a.d. que l’erreur diminuera avec le degré du polynôme – si la dérivée de f varie brutalement ou est discontinue, la fonction sera mieux approximée par un polynôme de degré faible, c.a.d. que l’erreur peut augmenter rapidement avec l’ordre du polynôme. D’une manière générale, l’utilisation de polynôme d’ordre n > 4, 5 donne rarement de bons résultats numériques. Algorithme 10 Interpolation par un polynôme de Lagrange d’ordre n Requiert: une fonction f définie en x0 . . . xn , l’ordre du polynôme n, la valeur xi nter o` u l’on souhaite approximer la fonction Fournit: l’approximation de f au point souhaité lagr ← 0 pour i = 0 . . . n faire li ← 1 pour j = 0 . . . n faire si j 6= i alors li ← li ∗ (xinter − xj )/(xi − xj ) fin si fin pour lagr ← lagr + li ∗ f (xi ) fin pour renvoyer lagr


H. Klein

Chapitre 7

Int´ egration de fonctions Dans ce chapitre, nous n’allons pas aborder la détermination des primitives des fonctions (ceci concerne le calcul formel, et des logiciels comme MAPLE, MATHEMATHICA...), mais le calcul numérique d’une intégrale Z

b

f (x)dx

a

qui revient à calculer l’aire sous la courbe f (x) entre a et b. Nous allons aborder 3 méthodes : – intégration par la méthode des rectanges – intégration par la méthode des trapèzes – intégration par la méthode de Simpson ces méthodes sont basées sur l’interpolation d’ordre 0,1,2 respectivement. Ces méthodes sont généralisables à l’ordre n. Nous terminerons ce chapitre par une discussion sur la recherche d’un pas optimal d’intégration.

7.1

Int´ egration par la m´ ethode des rectanges

Soit une fonction f (x) dépendant d’une variable indépendante. Nous souhaitons calculer I=

Z

xn

f (x)dx

x0

numériquement, cela revient ` a calculer la somme I=

n−1 X i=0

fi (xi+1 − xi )

Pour cela, nous allons tabuler f (x) en n + 1 points f0 , f1 . . . fn dans l’intervalle d’intégration [x0 , xn ]. On peut approximer l’aire sous la courbe f (x) par une série de rectangles comme représenté sur la figure 7.1, ce qui mathématiquement correspond à l’interpolation de f par un pˆ olynome d’ordre 0 (une constante). Si la distribution des points est uniforme, ce qui revient à xi+1 − xi = h, l’intégration revient ` a calculer I=h

n−1 X i=0

41

fi

´ CHAPITRE 7. INTEGRATION DE FONCTIONS

42

Fig. 7.1 – Découpage de l’aire sous la courbe f par une série de rectangle. Il y a cependant 2 possibilités pour réaliser ce calcul , comme cela est montré à la figure 7.2. On peut en effet majorer ou minorer la fonction selon que l’on calcule respectivement I=h

n−1 X

fi+1

i=0

ou I=h

n−1 X

fi

i=0

ceci dans le cas o` u f ′ (x) > 0. Ceci pose deux problèmes. Tout d’abord les résultats ne seront pas identiques suivant que l’on majore ou minore f . On peut éventuellement résoudre le problème en réalisant les calculs dans les deux cas, puis en faisant une moyenne des deux résultats. Ensuite, pour les mêmes raisons, on n’obtiendra pas les mêmes résultats selon que l’on intègre de a → b ou de b → a. Cette limitation doit être prise en compte lors de calculs utilisant cette méthode.

Fig. 7.2 – Il y a deux manières de délimiter l’aire sous la courbe par une série de rectangle : fi (xi+1 − xi ) (minoration) ou fi+1 (xi+1 − xi ) (majoration)


H. Klein


7.1.1

43

Incertitude de la m´ ethode

Nous avons vu que cette méthode revient à approximer la fonction f par un polynôme d’ordre 0. Dans ce cas, l’incertitude sur chaque point fi s’exprime ǫ ≤ h|f ′ (x)|max , donc pour chaque intervalle h d’intégration h2 ′ |f (x)|max ǫ≤ 2 L’incertitude totale pour n pas d’intégration sera donc ǫ ≤ n.

h2 ′ |f (x)|max 2

Cette méthode est extrêment simple a` mettre en œuvre, et ne nécessite pas de gros calculs. Cependant si l’on souhaite conserver une incertitude acceptable, il est généralement nécessaire d’utiliser un pas d’échantillonage petit.

7.2

Int´ egration par la m´ ethode des trap` ezes

Le principe de cette méthode est analogue à celui de l’intégration par une série de rectangles. Dans cette méthode on cherche ` a augmenter la précision du calcul tout en évitant la disymétrie qui apparaˆıt dans la définition de l’intégration par rectangles. Nous allons donc ici aussi tabuler f (x) en n + 1 points f0 , f1 . . . fn dans l’intervalle d’intégration [x0 , xn ]. On va ensuite interpoler linéairement (polynôme d’ordre 1, c.a.d. une droite) pour approximer f (x) entre deux points tabulés. On va donc remplacer les arcs xi , xi+1 par leurs cordes (figure 7.3). L’intégration va

Fig. 7.3 – Découpage de l’aire sous la courbe f par une série de trapèzes. alors se faire en sommant l’aire des trapèzes, soit I=

n−1 X i=0

fi+1 + fi (xi+1 − xi ) 2

si l’échantillonage des points est régulier, c.a.d. si xi+1 − xi = h alors l’intégrale s’exprimera I=

n−1 X h (f0 + fn ) + h. fi 2 i=1

en effet chaque trapèze aura chacun de ces côtés en commun avec les trapèzes adjacents, à l’exception du premier et du dernier qui n’ont qu’un côté en commun. M´ ethodes Informatiques pour la physique

H. Klein


7.2.1

44


Nous avons vu que cette méthode revient à approximer la fonction f par un polynôme d’ordre 1. Dans ce cas, l’incertitude sur chaque point fi s’exprime h2 ′′ |f (x)|max 2

ǫ≤ donc pour chaque intervalle h d’intégration ǫ≤

h3 ′′ |f (x)|max 4

L’incertitude totale pour n pas d’intégration sera donc ǫ ≤ n.

h3 ′′ |f (x)|max 4

La précision est supérieure ` a celle de l’intégration par rectangles, le terme d’erreur étant maintenant le terme du deuxième ordre du développement de f en série de Taylor.

7.3

Int´ egration par la m´ ethode de Simpson

Comme les autres méthodes, cette méthode est basée sur une distribution fi de la fonction f (x) que l’on interpole. La méthode de Simpson correspond au cas de l’interpolation parabolique (polynôme d’ordre 2 passant par 3 points fi ). Le polynôme d’interpolation Pi (x)d’ordre 2 a été défini au chapitre 6.2, et l’on peut donc exprimer Z

xi+1

f (x)dx =

Z

xi+1

xi−1

xi−1

Pi (x)dx = h

Z

1 −1

Pi (u)du

en posant x = xi+1 + u.h, donc si u = −1 x = xi , si u = 0 x = xi+1 , et si u = 1 x = xi+2 . On peut alors calculer Z xi+1 h f (x)dx = (fi−1 + 4fi + fi+1 ) 3 xi−1 l’expression de l’intégrale sur un intervalle de calcul de 2n + 1 points s’exprime alors I=

2n hX (fi−1 + 4fi + fi+1 ) 3 i=1

que l’on peut réécrire plus commodemment

I=





n−1 n−1 X X  h (f1 + fn ) + 4.f + 2.fi  i   3 i>1 i>1 pair


impair

H. Klein


7.3.1

45


Nous avons vu que cette méthode revient à approximer la fonction f par un polynôme d’ordre 2. Dans ce cas, l’incertitude sur chaque point fi est majorée par le terme d’ordre 3 du développement en série de Taylor de f h3 (3) |f (x)|max 6 donc pour chaque intervalle h d’intégration ǫ≤

h4 (3) |f (x)|max 12 L’incertitude totale pour n pas d’intégration sera donc ǫ≤

ǫ ≤ n.

7.4

h4 (3) |f (x)|max 12

M´ ethode de Romberg

La méthode de Romberg génère un tableau triangulaire d’estimations numériques de l’intégrale en applicant l’extrapolation de Richardson à la méthode des trapèzes. Cela permet d’améliorer l’ordre de convergence de la méthode des trapèzes en divisant l’intervalle d’étude et en formant des combinaisons judicieuses des différents termes calculés. Nous n’allons pas détailler cette méthode en analyse numérique, mais juste la définir. Le principe est le suivant : si l’on veut trouver une estimation de f (x) intégrée de a à b on calcule le premier terme R(0, 0) par la méthode des trapèzes, puis les termes suivants en somant des nouveaux termes et en faisant des combinaisons avec les termes calculés précédemment : R(0, 0) =

1 2 (b

R(n, 0) =

P2n−1 1 k=1 2 R(n − 1, 0) + h 1 m 4m −1 (4 R(n, m − 1) −

R(n, m) =

− a)(f (a) + f (b)) f (a + (2k − 1)h)

R(n − 1, m − 1))

avec h = (b − a)/2n , n ≥ 1, m ≥ 1. On utilise alors les termes diagonaux R(n, n) comme estimation de l’intégrale avec une erreur en O(h2n+1 . Si vous développez ` a la main les équations précédentes, vous constaterez que les valeurs de R situées dans la première colonne, correspondent à la méthode des trapèzes pour nu intervalle divisé en 2 puis en 3 . . .. La deuxième colonne correspond aux mêmes calculs avec la méthode de Simpson, et ainsi de suite. Généralement, pour implémenter cette méthode, on poursuit le calcul jusqu’à ce que (critère de convergence) la différence entre 2 valeurs successives d’une ligne soit inférieure ` a une précision souhaitée. Les premiers termes R(n, 0) sont calculés par la méthode des trapèzes, avec un nombre de pas qui augmentent à chaque itération, les termes suivants sont estimés ` a partir des calculs précédents. Les valeurs numériques suivantes correspondent à l’estimation de Z

1.25

cos(x2 )dx = 0.9774

0


H. Klein


n=1 n=2 n=3 n=4 n=5

m=1 6.301851e-01 8.930121e-01 9.406354e-01 9.568902e-01 9.643337e-01

46

m=2

m=3

m=4

m=5

9.806211e-01 9.565098e-01 9.623084e-01 9.668149e-01

9.549024e-01 9.626950e-01 9.671153e-01

9.628187e-01 9.671855e-01

9.672026e-01

le calcul a été arrêté lorsque m = 5. Vous voyez cependant que l’on est encore loin d’une estimation correcte : la valeur exacte est 0.9774 soit une erreur de ǫ = 1.02e − 1.

7.5

D´ etermination d’un pas optimal d’int´ egration

La question ici est : lorsque l’on cherche à calculer In =

Z

a

f (x)dx

b

quel est le pas optimal h = xi+1 −xi permettant de calculer In avec une précision ǫ donnée. Nous avons jusqu’à présent considéré le pas h comme une constante. Reprenons le cas de l’intégration par trapèze. Si l’on souhaite déterminer In à la précision ǫ, il suffit de choisir un pas h0 , de calculer I0 , puis de diminuer le pas en h1 , calculer I1 ... Ceci tant que |In+1 − In | ≥ ǫ. La méthode est simple, mais les calculs sont répétitifs et potentiellement longs. Il existe un moyen de réduire le temps de calcul en introduisant une relation de récurrence permettant d’exprimer In+1 en fonction de In . Nous allons détailler ce calcul, toujours dans le cas d’une intégration par trapèzes.

1 On pose h0 = b − a, et on calcule I0 =

h0 [f (a) + f (b)] 2

2 On pose h1 = h0 /2 et on calcule I0 2

I1 = = =

z

}|

{

h1 h1 [f (a) + f (a + h1 ) + f (a + h1 ) + f (b)] = [f (a) + f (b)] +h1 f (a + h1 ) 2 2 1 I0 + h1 f (a + h1 ) 2 1 [I0 + h0 (f (a + h1 ))] 2

3 On pose h2 = h1 /2 = h0 /4 et on calcule I2 =

h2 [f (a) + f (a + h2 ) + f (a + h2 ) + f (a + 2h2 ) + f (a + 2h2 ) + f (a + 3h2 ) + f (a + 3h2 ) + f (b)] 2


H. Klein


=

h2 [f (a) + f (b)] +h2 [f (a + h2 ) + f (a + 2h2 ) +f (a + 3h2 )] {z } | |2 {z } f (a+h1 )

h1 [f (a)+f (b)] 4

= =

47

{z

|

I1 2

1 I1 + h2 [f (a + h2 ) + f (a + 3h2 )] 2 1 {I1 + h1 [f (a + h2 ) + f (a + 3h2 )]} 2

}

4 On pose h3 = h2 /2 = h0 /8 h3 [f (a) + f (b) + 2f (a + h3 ) + 2f (a + 2h3 ) + 2f (a + 3h3 ) + 2f (a + 4h3 ) 2 +2f (a + 5h3 ) + 2f (a + 6h3 ) + 2f (a + 7h3 )] 1 I2 + h3 [f (a + h3 ) + f (a + 3h3 ) + f (a + 5h3 ) + f (a + 7h3 )] = 2 1 = {I2 + h2 [. . .]} 2

I3 =

On peut alors généraliser en In =



2n−1 X



1 f (a + (2i − 1)hn ) In−1 + hn−1 2 i=1impair

Cette relation permet de calculer économiquement le terme n connaissant le terme précédent. D’une manière générale, un programme réalisant une intégration sera divisé de la manière suivante : – une fonction permettant de réaliser un pas d’intégration – une fonction réalisant le calcul de l’intégrale en utilisant des relations de récurrence du type de celle que nous venons de détailler pour fournir un résultat avec une précision connue.


H. Klein

Chapitre 8

D´ erivation num´ erique Nous allons nous intéresser dans ce chapitre au calcul de la dérivée d’une fonction échantillonée, c.a.d d’une fonction dont seulement certaines valeurs sont connues. La définition mathématique de la dérivée f ′ (x) d’une fonction f (x) définie sur un intervalle I avec x ∈ I est lim

f ′ (x) =ǫ → 0 ou

f (x + ǫ) − f (x) ǫ

f (x) − f (x − ǫ) ǫ selon que l’on réalise une estimation ` a droite ou à gauche respectivement. Comme l’on calcule la limite pour ǫ → 0, les deux expressions sont strictement équivalentes. lim

f ′ (x) =ǫ → 0

8.1

Formule d’Euler

En calcul numérique en revanche, si les mêmes définitions sont utilisées, la situation est différente puisque l’on va maintenant travailler sur des différences finies, c.a.d que ǫ sera égal à une constante non nulle, ǫ = h. Les deux définitions précédentes ne sont alors plus équivalentes. Il faudra utiliser f (x + h) − f (x) f ′ (x) = h ou f (x) − f (x − h) f ′ (x) = h

8.1.1

Erreur intrins` eque

Si l’on réalise un développement en série de Taylor de f (x + h) on obtient f (x + h) = f (x) + h.f ′ (x) +

h2 ′′ .f (x) + . . . 2

il apparaˆıt donc que l’estimation de la dérivée f ′ (x) =

f (x + h) − f (x) h 48

´ ´ CHAPITRE 8. DERIVATION NUMERIQUE

49

présente une erreur correspondant au terme d’ordre 2 du dévellopement de Taylor, cette erreur est un terme h2 .|f ′′ (x)|. On voit donc que cette méthode n’est pas très précise, elle donne de plus des résultats différents suivant que l’on réalise une estimation à gauche ou à droite. Pour ces raisons, cette méthode est peu employée. Nous allons maintenant voir qu’il est possible de l’améliorer à peu de frais.

8.2

D´ erivation au demi-pas

En conservant la même démarche, nous allons refaire en calcul en symétrisant le problème c.a.d. que pour estimer la dérivée au point x, nous allons utiliser les valeurs de la fonction f (x + h) et f (x − h). Si l’on écrit les développements limités de f (x) dans ces deux situations f (x + h) = f (x) + h.f ′ (x) +

h3 h2 ′′ .f (x) + f (3) + . . . 2 6

f (x − h) = f (x) − h.f ′ (x) +

h2 ′′ h3 .f (x) − f (3) + . . . 2 6

on en déduit que f (x + h) − f (x − h) = 2h.f ′ (x) + 2

h3 (3) f + ... 6

et il apparaˆıt alors que f ′ (x) ≃

f (x + h) − f (x − h) 2h 2

avec une erreur correspondant au terme h3 |f (3) (x)|. Nous sommes maintenant en possession d’une méthode symétrique, avec un terme d’erreur à l’ordre 2.

8.2.1

Stabilit´ e num´ erique

La seule question qu’il convient de se poser ici est : comment doit-on choisir le pas h utilisé dans le calcul ? La réponse intuitive est que ce pas doit être le plus petit possible, nous allons voir que c’est faux. Quant on discrétse une fonction continue, on commet une erreur numérique ǫ=

∆f f

o` u ∆f est l’écart entre la valeur réelle de la fonction et la valeur tabuléee. Pour le calcul au demi-pas, l’incertitude numérique pour les évaluations de f est Inumerique =

ǫ ǫ|f (x + h)| + ǫ|f (x − h)| = |f (x)| 2h h

L’incertitude totale comprenant l’erreur intrinsèque au calcul de la dérivée s’exprime Itotal = M´ ethodes Informatiques pour la physique

ǫ h2 |f (x)| + |f (3) (x)| h 3 H. Klein

´ ´ CHAPITRE 8. DERIVATION NUMERIQUE

50

L’incertitude optimale se définit telle que ∂Itotal =0 ∂h ce qui conduit à un pas de calcul optimal h=


3ǫf 2f (3)

1 3

H. Klein

Chapitre 9

R´ esolution d’´ equations diff´ erentielles lin´ eaires C’est un problème très fréquent en physique (phénomènes transitoires dans un circuit RLC par exemple), chimie (cinétique des réactions chimiques), astronomie . . .Le problème à l’ordre 1 peut-être résolu par les méthodes d’intégration du chapitre 7 dy = r(x) ⇔ y(x) = dx

Z

r(x)dx

Dans ce chapitre, nous considèrerons 3 types de bases d’équations différentielles respectivement décroissantes, croissantes ou oscillantes : dy + αy ⇔ y = y0 exp(−αx) dx dy − αy ⇔ y = y0 exp(αx) dx dy ± iωy ⇔ y = y0 exp(±iωx) dx et nous aborderons plusieurs méthodes numériques de résolution de ce type d’équation, en mettant l’accent sur le fait que toutes les méthodes ne sont pas adaptées à tous les types d’équation. Les problèmes physiques sont souvent d’ordre supérieur, comme dans l’équation dy d2 y +α = r(x) dx dx on peut facilement se ramener au problème à l’ordre 1 en réécrivant l’équation du deuxième ordre précédente comme un système d’équation du premier ordre dy dx dz dx

= z(x) = r(x) − αz(x)

51

´ ´ ´ ´ CHAPITRE 9. RESOLUTION D’EQUATIONS DIFFERENTIELLES LINEAIRES

52

système que nous pouvons résoudre par intégration. On peut généraliser à une équation différentielle d’ordre N , que l’on récrit comme un système de N équations fi d’ordre 1 dyi (x) = f (x, y1 . . . yN ) dx o` u les fonctions fi sont connues. Pour résoudre le système, il suffit d’intégrer les fonctions fi . Pour résoudre numériquement ce système, on part d’un point donné (x0 , y0 ) représentant les conditions initiales du problème. On choisit ensuite un pas d’intégration fini h = ∆x comme pour les méthodes du chapitre 7, et l’on intègre pas à pas avec xn = xn−1 + h yn = yn−1 + ∆y dans ce chapitre, nous allons détailler différentes méthodes pour évaluer les pas finis ∆y.

9.1

M´ ethode d’Euler

Soit

dyi = fi (x, y) dx

(1)

La formule la plus simple que nous puissions utiliser est la formule d’Euler. C’est à dire que dy et dx sont remplacés par des accroissements finis ∆y et ∆x, et l’on réécrit ensuite les équations en multipliant par ∆x, ela revient ` a écrire, pour un accroissement h y(x + h) = y(x) + hy(x) + O(h2 ) le terme O(h2 ) étant le terme d’erreur (` a l’ordre 2) On peut alors réécrire l’équation (1) en négligeant les termes du deuxième ordre yi+1 = yi + hf (xi , yi ) dans le cas o` u l’on ne considère qu’un échantillonage de x x0 . . . xn , avec h = xi+1 − xi . Si l’on connait les valeurs initiales {x0 , y0 }, cette méthode permet de calculer toutes les valeurs successives yi , comme on le voit sur la figure 9.1 Avec une méthode de ce type, on fait avancer la solution de xi à xi + h = xi+1 à chaque pas d’intégration. Nous avons cependant vu au chapitre 6 que cette méthode est peu performante, et surtout asymétrique. Cela peut poser quelques problèmes : si par exemple vous tracez une trajectoire dépendant d’une équation différentielle entre t = 0 et t, vous obtiendrez des résultats différents en réalisant les mêmes calculs entre t et t = 0... Cette méthode est donc ` a éviter dans la majorité des cas. Il est utile de l’introduire, pour permettre la compréhension des méthodes d’intégration des équations différentielles, et elle reste cependant utilisable dans certains cas simples. M´ ethodes Informatiques pour la physique

H. Klein


53

Fig. 9.1 – Méthode d’Euler. Dans cette méthode d’intégration d’une équation différentielle, la dérivée au point de départ xn de chaque intervalle est extrapolée pour trouver la valeur suivante en xi+1 de la fonction.

9.1.1

Pr´ ecision et stabilit´ e de la m´ ethode

Quelle est la précision de la méthode d’Euler ? Pour y répondre, nous allons considérer le dévelopement en série de Taylor de la fonction y(x) autour de la valeur xi yi+1

dy = yi + h dx

h2 + 2 i

d2 y dx2

!

+ ... i

si l’on soustrait cette expression ` a l’expression de yi+1 défini au chapitre 9.1 on obtient yi+1

dy = yi + h dx

h2 + 2 i

d2 y dx2

!

dy + . . . ≈ yi − hf (yi , xi ) = yi + h dx i

i

Nous voyons donc que le terme en h est correctement exprimé dans l’approximation de la méthode d’Euler, mais que le terme d’ordre supérieur est faux. Ceci nous permet de décrire la méthode d’Euler comme une méthode du premier ordre. Le terme d’erreur parfois appelé erreur d’arrondi est ici le terme en h2 , on peut réécrire l’équation de la manière suivante y(x + h) = y(x) + hy(x) + O(h2 ) le terme O(h2 ) étant le terme d’erreur (` a l’ordre 2) Nous venons de voir que la méthode d’Euler est une méthode du premier ordre. Il reste à se poser une question importante : quelle est sa stabilité ? Supposons qu’à un moment du calcul, la solution numérique diffère de la vraie solution d’une quantité δy. Nous écrivons

yi+1 + δyi+1 = yi + δyi − h f (yi , xi ) +

∂f ∂y

i

δyi

o` u le terme entre crochets est le développement de f en série de taylor, en soustrayant ceci à l’expression de yi+1 , on obtient δyi+1 M´ ethodes Informatiques pour la physique

∂f = 1−h ∂y

δyi = gδyi

i

H. Klein


54

si g a une magnitude supérieure ` a 1, alors δyi va avoir tendance à augnmenter quand i augmente. Il y a alors un risque que cette erreur devienne prédominante devant la solution recherchée. On dira que la méthode d’Euler est stable si |g| ≤ 1 ou −1 ≤ 1 − h

∂f ≤1 ∂y

h étant positif par définition, la deuxième inégalité impose donc que la dérivée de f par rapport à y soit elle aussi positive. La première inégalité introduit une restriction sur δx h ≤ 2/

∂f ∂y

Si cette dérivée est complexe, il faut faire attention au calcul de |g|. Dans ces cas il peut être plus facile de rechercher des solutions ` a la condition |g2 | ≤ 1. Dans le cas d’une équation d’oscillation du type dy ± iωy = 0 ⇔ y = yo exp(±iωx) dx cette inégalité conduira ` a 1 + h2 ω 2 ≤ 1 qui est impossible ` a satisfaire pour des ω et h réels. De même pour une équation du type dy − αy = 0 ⇔ y = yo exp(αx) dx cette analyse n’a pas de sens. Dans ce cas il serait plus judicieux de s’intéresser à l’erreur relative, et imposer que cette erreur dans la solution n’augmente pas. décroissance h ≤ 2/α

9.1.2

croissance instable

oscillation instable

Exemple d’application : la loi de d´ ecroissance radioactive

Pour illustrer ce propos, voici un exemple dans un cas extrêmement simple : la loi de décroissance radioactive, o` u nous sommes capables de donner une solution analytique exacte très facilement. Une substance radioactive est caractérisée par une constante de décroissance radioactive k. On peut écrire la loi N (t) = N0 exp(−k.t) permettant d’exprimer le nombre de noyaux radioactifs a` un instant t, N (t) en fonction du nombre de noyaux N0 ` a l’instant t = 0 et de la constante de décroissance radioactive. Cette loi correspond à l’équation différentielle dN dN + k.N = 0 ⇔ = −k.N = f (t) dt dt M´ ethodes Informatiques pour la physique

H. Klein


55

Nous allons résoudre numériquement cette équation avec k = 100s−1 et N0 = 1 en utilisant la méthode d’Euler. Nous avons montré précédemment que cette méthode était stable si le pas d’intégration h (ici un temps) respectait le critère h ≤ 2/α, dans ce cas précis α = k. Il est donc nécessaire dans ce cas de choisir un pas h ≤ 2.10−2 s pour satisfaire le critère de stabilité. La figure 9.2 représente les résultats de la simulation pour des pas d’intégration inférieurs ou supérieurs à ce critère, ainsi que la solution exacte en trait plein. 1 dt = 1 ms -100*t e

dt=30ms dt=1ms

4

0,8

N /N0

N / N0

0,6

2

0,4

0 0,2

0

0

0,05 Temps (s)

0,1

-2

0 Temps (s)

Fig. 9.2 – Application de la méthode d’Euler : simulation de l’évolution temporelle du nombre de noyaux d’une source radioactive (constante de décroissance radioactive 100s−1 . Le critère de stabilité de la méthode impose un pas d’intégration inférieur à 0.02 s. Courbe de gauche : le critère est respecté, le calcul numérique donne des résultats corrects (la solution exacte est tracée en trait plein). Lorsque ce critère n’est pas respecté (courbe de droite) le résultat diverge. Un simple examen de ce graphe permet de comprendre qu’il est important de choisir un pas d’intégration satisfaisant le critère de stabilité de la méthode, sous peine de voir le calcul diverger rapidement . . .On voit donc qu’une méthode peut donner de bons résultats dans certains cas, et pas dans d’autres. Il convient donc de se poser ces questions avant de se lancer dans des calculs complexes. . .

9.2

M´ ethode du saut de grenouille

L’idée ici est d’améliorer la méthode d’Euler en symétrisant le problème. On va simplement écrire yi+1 = yi−1 + 2h.f (xi , yi ) ici, le terme d’erreur est en h3 /6, on parlera donc d’une méthode d’ordre2.

9.2.1

Stabilit´ e

Si l’on tient le même raisonnement que pour la méthode d’Euler, on peut exprimer la variation de la solution par rapport ` a la solution vraie pour une équation décroissante : δyi+1 = δyi−1 − 2h M´ ethodes Informatiques pour la physique

∂f δyi ∂y H. Klein


56

on va écrire δyi = gδyi−1 et δyi+1 = g2 δyi−1 ce qui conduit à g2 = 1 − 2h. dont les solutions s’écrivent ∂f g = h. ± ∂y

s

∂f g ∂y

∂f 1 + h. ∂y

2

on a toujours, pour les deux racines g1 .g2 = −1 et g1 6= g2 ce qui implique que l’une des deux racines est supérieure ` a 1. La méthode est donc instable dans ce cas, ainsi que dans le cas d’une équation croissante. Il y a une exception dans le cas ou ∂f ∂y est purement imaginaire, le terme sous la racine dans l’expression de g peut alors être négatif, il est alors possible que g1 .g2 < 1. C’est le cas pour une équation oscillante o` u ∂f ∂y = ±iω. Dans ce cas, l’algorithme est stable si le terme sous la racine dans l’expression de g est positif, c’est à dire si h ≤ 1/ω décroissance instable

9.3

croissance instable

oscillation h ≤ 1/ω

Algorithme Verlet-vitesses

Les algorithmes de Verlet sont très utilisés pour intégrer les équations du mouvement de Newton et calculer les trajectoires de particules en dynamique moléculaire (et dans les jeux vidéos). L’algorithme du Verlet réduit les erreurs par rapport à une méthode calculant la position au temps t + h uniquement ` a partir de la position aux temps t et t − h. C’est une méthode à l’ordre 3 (erreur en O(h4 ). Le principe est d’écrire le développement limité de la position dans les deux directions temporelles : x(t + h) = x(t) + v(t).h +

b(t).h3 a(t).h2 + + O(h4 ) 2 6

a(t).h2 b(t).h3 − + O(h4 ) 2 6 o` u h est le pas d’intégration, x est la position, v la vitesse, a l’accélération, et b la dérivée troisième de la position. En sommant les deux expressions précédentes, on obtient x(t − h) = x(t) − v(t).h +

x(t + h) = 2.x(t) − x(t − h) + a(t).h2 + O(h4 ) Une variante de cet algorithme, que nous détaillons ici est l’algorithme Verlet vitesses, basé sur le même principe pais dans lequel on inclue explicitement la vitesse : 1 x(t + h) = x(t) + v(t).h + a(t).h2 2 a(t) + a(t + h) .h 2 on peut montrer que l’erreur est du même ordre que pour l’algorithme du Verlet. On l’implémente généralement de la manière suivante : v(t + h) = v(t) +


H. Klein


1. on calcule

57

1 x(t + h) = x(t) + v(t).h + a(t).h2 2

2. on calcule v(t + h/2) = v(t) + a(t).h/2 3. à partir du potentiel d’interaction on calcule a(t + h) 4. et finalement on calcule v(t + h) = v(t + h/2) +

a(t + h).h 2

Cet algorithme suppose que l’accélération a(t + h) ne dépend que de x(t + h), et pas de v(t + h).

9.4

M´ ethode de Runge Kutta ` a l’ordre 2

Nous avons vu deux méthodes : l’une est stable pour des équations décroissantes, l’autre pour les équations oscillantes en respectant certaines conditions. Est-il possible de combiner ces deux conditions de stabilité dans une seule méthode ? La réponse est oui, avec la méthode de Runge Kutta à l’ordre 2. Le principe est schématisé sur la figure 9.3. Les valeurs des dérivées des yi au point initial sont utilisées pour estimer les yi au centre du pas d’intégration, puis les valeurs de ces dérivées sont utilisées pour faire avancer la solution jusqu’au point xi+1 .

Fig. 9.3 – Méthode de Runge Kutta d’ordre 2. Dans cette méthode d’intégration d’une équation différentielle, on obtient une précision au deuxième ordre en utilisant la dérivée au point de départ xi pour trouver un point intermédiaire. On utilise ensuite la dérivée en ce point pour trouver la valeur suivante de la fonction. On calcule la dérivée de y(x) au point xi k1 = M´ ethodes Informatiques pour la physique

dy dx

= h.f (xi , yi )

x=xi H. Klein


58

puis on calcule la dérivée au point xi + h/2 k2 =

dy dx

x=xi + h 2

= h.f (xi +

k1 h , yi + ) 2 2

on écrit alors que yi+1 = yi + k2 + O(h3 ) c’est une méthode du deuxième ordre, l’erreur correspond donc au terme suivant d’ordre 3. On voit que cette méthode requiert deux évaluations de f (x) à chaque pas d’intégration. Contrairement ` a la méthode d’Euler, cette méthode est symétrique, c.a.d que l’intégration de l’équation différentielle de x0 → xi ou de xi → x0 conduira au même résultat.

9.4.1

Stabilit´ e de la m´ ethode

On exprime comme précédemment l’écart à la vrai solution δyi+1 = δyi

"

1 ∂f ∂f + h 1−h ∂y 2 ∂y

2 #

ce qui conduit en posant que le terme entre crochets est inférieur à 1 aux conditions de stabilité suivantes décroissance h ≤ 2/α

croissance instable

oscillation 1 + 14 (hω)4 ≤ 1

nous avons ici combiné les critères de stabilité des méthodes d’Euler et de la grenouille au sein d’une seule méthode d’ordre 2. Dans le cas d’une équation oscillante, on peut considérer la méthode comme instable, mais en pratique si hω < 1 les erreurs générées sont négligeables devant le vrai résultat, et la méthode est utilisable. Exemple d’application Reprenons l’exemple présenté pour la méthode d’Euler, dN dN + k.N = 0 ⇔ = −k.N = f (t) dt dt Nous allons résoudre numériquement cette équation avec k = 100s−1 en utilisant la méthode RK2. Nous avons montré précédemment que cette méthode était stable si le pas d’intégration h (ici un temps) respectait le critère h ≤ 2/α, dans ce cas précis α = k. Il est donc nécessaire dans ce cas de choisir un pas h ≤ 2.10−2 s pour satisfaire le critère de stabilité. La figure 9.4 représente les résultats de la simulation pour des pas d’intégration inférieurs ou supérieurs à ce critère, ainsi que la solution exacte en trait plein. A nouveau,un simple examen du graphe permet de comprendre qu’il est important de choisir un pas d’intégration satisfaisant le critère de stabilité de la méthode, sous peine de voir le calcul diverger rapidement. On notera cependant que le comportement est différent de celui observé avec la méthode d’Euler. M´ ethodes Informatiques pour la physique

H. Klein


59

2 -2

pas de temps : 2,1.10 -3

N / N0

pas de temps 10

s

s

1

0

0

0,1 Temps (s)

0,2

Fig. 9.4 – Application de la méthode Runge Kutta d’ordre 2 : simulation de l’évolution temporelle du nombre de noyaux d’une source radioactive (constante de décroissance radioactive 100s−1 . Le critère de stabilité de la méthode impose un pas d’intégration inférieur à 0.02 s. Lorsque ce critère n’est pas respecté, le résultat diverge par rapport à la solution exacte tracée en trait plein.

9.5

M´ ethode Runge Kutta ` a l’ordre 4

Dans la méthode RK2, on utilise l’estimation des yi au demi-pas d’intégration pour le calcul des dérivées. La méthode de Runge-Kutta d’ordre 4 tente de fournir une meilleure estimation à l’aide d’une moyenne de 4 estimations, ce principe est schématisé sur la figure 9.5. La première estimation est celle d’Euler : on fait un demi-pas de longueur h/2 par la méthode d’Euler pour évaluer les valeurs yi au demi-pas. Cela permet de calculer de nouvelles valeurs des dérivées. On revient ensuite au point d’origine et l’on utilise ces dérivées pour refaire un demi-pas suivant la méthode d’Euler et ainsi obtenir de nouvelles valeurs des yi : on ne trouvera pas les mêmes valeurs qu’à la première étape puisque l’on n’utilise pas les mêmes valeurs des dérivées (n’oublions pas qu’il ne s’agit que d’estimations). Cela permet d’obtenir une troisième estimation des dérivées. Cette estimation est utilisée pour faire un pas h complet suivant la méthode d’Euler, toujours ` a partir du même point de départ, cela conduit à une quatrième estimation des dérivées des yi . Il ne reste plus alors qu’à faire une moyenne pondérée de ces quatres estimations en affectant un poids double aux deux estimations intermédiaires, et faire un pas complet de longueur h pour faire avancer la solution jusqu’au point xi+1 . L’algorithme 11 propose une implémentation de cette méthode. Cette méthode requiert 4 évaluations de f (x, y) par pas d’intégration h : une au point de départ de l’intervalle, 2 au demi-intervalle, et une à la fin de l’intervalle : k1 = h.f ()xi , yi ) k1 h k2 = h.f (xi + , yi + ) 2 2 M´ ethodes Informatiques pour la physique

H. Klein


60

Fig. 9.5 – Méthode Runge Kutta d’ordre 4. Dans cette méthode d’intégration d’une équation différentielle, on obtient une précision au quatrième ordre en évaluant 4 fois la dérivée : une fois au point de départ, deux fois au milieu de l’intervalle, et une fois en un point final estimé. h k2 ), yi + ) 2 2 = h.f (xi + h, yi + k3 )

k3 = h.f (xi + k4

et l’on calcule alors la valeur de y au point xn+1 en sommant ces quatres termes yi+1 = yi +

k1 k2 k3 k4 + + + + O(h5 ) 6 3 3 6

le terme d’erreur est alors d’ordre 5. Cette méthode peut paraitre ressembler à un bricolage, mais c’est une des méthodes les plus utilisées en physique : elle est en effet stable et efficace dans beaucoup de cas. Pour des raisons de lourdeur des calculs, nous ne détaillerons pas les critères de stabilité de cette méthode, qui est cependant la plus robuste de toute celles que nous venons d’aborder.


H. Klein


61

Algorithme 11 Résolution d’un système d’équations différentielles du premier ordre par la méthode Runge-Kutta d’ordre 4 Requiert: un pas d’intégration h, les n valeurs yi du début d’intervalle x dans un tableau, une fonction deriv permettant le calcul des dérivées Fournit: les valeurs fi au point x + h dh ← h/2 d1 ← deriv(x, y[], n) pour i = 1 à n faire ytemp [i] = y[i] + d1 [i] ∗ dh fin pour d2 ←deriv(x + dh, ytemp [], n) pour i = 1 à n faire ytemp [i] = y[i] + d2 [i] ∗ dh fin pour d3 ← deriv(x + dh, ytemp [], n) pour i = 1 à n faire ytemp [i] = y[i] + d3 [i] ∗ h fin pour d4 ← deriv(x + h, ytemp[] , n) pour i = 1 à n faire y[i] ← y[i] + h d16[i] + d23[i] + d33[i] + d46[i] fin pour renvoyer y


H. Klein

Annexe A

R´ ef´ erences bibliographiques et informatiques La littérature couvrant ce domaine est importante, je ne citerai ici que quelques références qui m’ont servi ` a la réalisation de ce document. Libre à vous d’en utiliser d’autres. 1. Numerical Recipes in C, second edition, Cambridge University Press, 1992 W.H. Press, S.A. Teukolski, W.T. Vetterling, B.P. Flannery La bible pour l’utilisateur ! Le livre peut être téléchargé au format PDF à l’adresse http ://libwww..lanl.gov/numerical/bookcpdf.html 2. Physique numérique, Ed. Vuibert, 1998 P. Depondt Malheureusement épuisé 3. Méthodes de calcul numérique, Ed. Hermès, 2001 J.P. Nougier Vous trouverez ici quelques références d’outils informatiques libres. Cela signifie que ces outils sont gratuits, et que vous disposez de leurs codes sources. Libre à vous de les utiliser, de lire les sources pour les comprendre . . .Ces outils développés sur platte-forme UNIX, fonctionnent aussi sous Windows (sauf Linux bien sˆ ur). 1. Un système d’exploitation stable, fonctionnant sur architecture PC : Linux. Voir par exemple le site de la distribution Mandrake http ://www.linux-mandrake.com ou http ://www.linux.org 2. La référence des compilateurs C : gcc, il compile aussi le C++, l’ObjC . . . http ://www.gnu.org/software/gcc 3. Besoin d’une librairie de calcul numérique ? Pensez à la GNU Scientific Library http ://www.gnu.org/software/gsl 4. Un logiciel de tracé de courbe très complet : Gnuplot, il fonctionne même sous Windows http ://www.gnuplot.info 62

List of Algorithmes 1 2 3 4 5 6 7 8 9 10 11

Tri d’un tableau par tri ` a bulle . . . . . . . . . . . . . . . . . . . . . . . . . . . . Tri d’un tableau par insertion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Tri d’un tableau par la méthode du tri shell . . . . . . . . . . . . . . . . . . . . . Tri d’un tableau par la méthode du tri rapide. Dans cet algorithme, le pivot est le dernier élément du tableau. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Encadrement d’une fonction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Réduction de l’encadrement d’une fonction . . . . . . . . . . . . . . . . . . . . . Recherche d’une racine par dichotomie . . . . . . . . . . . . . . . . . . . . . . . . Recherche de racines par la méthode de Newton-Raphson . . . . . . . . . . . . . Recherche d’un minimum par la methode golden section search . . . . . . . . . . Interpolation par un polynôme de Lagrange d’ordre n . . . . . . . . . . . . . . . Résolution d’un système d’équations différentielles du premier ordre par la méthode Runge-Kutta d’ordre 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

11 12 14 15 25 26 27 29 34 40 61