L'analyse numérique est une discipline des mathématiques. Elle s'intéresse tant aux fondements théoriques qu'à la mise en pratique des méthodes permettant de résoudre, par des calculs purement numériques, des problèmes d'analyse mathématique.

Définition: "L'analyse numérique" est l'étude des algorithmes permettant de résoudre les problèmes de mathématiques continues (distinguées des mathématiques discrètes).

Cela signifie qu'elle s'occupe principalement de répondre numériquement à des questions à variable réelle ou complexe comme l'algèbre linéaire numérique sur les champs réels ou complexes, la recherche de solutions numériques d'équations différentielles et d'autres problèmes liés survenant dans les sciences physiques et l'ingénierie.

Certains problèmes de mathématique continue peuvent être résolus de façon exacte par un algorithme. Ces algorithmes sont appelés alors "méthodes directes". Des exemples sont l’élimination de Gauss-Jordan pour la résolution d'un système d'équations linéaires et l'algorithme du simplexe en programmation linéaire (voir plus loin).

Cependant, aucune méthode directe n'est connue pour certains problèmes (et il est même démontré que pour une classe de problèmes dits "NP complets", il n'existe aucun algorithme fini de calcul direct en temps polynômial). Dans de tels cas, il est parfois possible d'utiliser une méthode itérative pour tenter de déterminer une approximation de la solution. Une telle méthode démarre depuis une valeur devinée ou estimée grossièrement et trouve des approximations successives qui devraient converger vers la solution sous certaines conditions. Même quand une méthode directe existe cependant, une méthode itérative peut être préférable car elle est souvent plus efficace et même souvent plus stable (notamment elle permet le plus souvent de corriger des erreurs mineures dans les calculs intermédiaires).

L’utilisation de l’analyse numérique est grandement facilitée par les ordinateurs. L’accroissement de la disponibilité et de la puissance des ordinateurs depuis la seconde moitié du 20ème siècle a permis l’application de l’analyse numérique dans de nombreux domaines scientifiques, techniques et économiques, avec souvent des effets révolutionnaires.

Lors de simulations numériques de systèmes physiques, les conditions initiales sont très importantes dans la résolution d'équations différentielles (voir les différents chapitres du site ou apparaissent des effets chaotiques). Le fait que nous ne puissions les connaîtres avec exactitude fait que les résultats des calculs ne peuvent jamais êtres parfaitement exactes (nous le savons très bien pour la météo qui en est l'exemple connu le plus flagrant). Cet effet, est une conséquence des résultats de la physique fondamentale (basée sur des mathématiques pures) qui démontre que l'on ne peut connaître parfaitement un système en y effectuant des mesures puisqu'elles perturbent directement ce dernier (principe d'incertitude de Heisenberg) et elles font l'objet de la théorie du Chaos (classique ou quantique).

Avec les nouveaux outils informatiques à disposition en ce début de 21ème siècle, il est devenu pratique et passionnant de connaître les méthodes numériques afin de s'amuser dans certains logiciels (OpenGL, 3D Studio Max, Maple, Matlab, Mathematica, Comsol, etc.) à simuler sous forme graphique 2D ou 3D des systèmes physiques.

Remarques:

R1. Beaucoup de méthodes numériques utilisées en informatique se basent sur des raisonnements qui ont déjà été étudiés dans d'autres sections de ce site et sur lesquelles nous ne reviendrons pas.

R2. Ce chapitre étant à la limite entre l'ingénierie et la mathématique appliquée, nous avons décidé de donner parfois des exemples d'applications des outils développés.

Définition: Un "algorithme" est une suite finie de règles à appliquer dans un ordre déterminé à un nombre fini de données pour arriver, en un nombre fini d’étapes (dont la quantité, ou réciproquement le temps d'exécution est définie par le terme "coût"), à un certain résultat, et cela indépendamment du type de données

Les algorithmes sont intégrés dans des calculateurs par l'intermédaire de "programmes".

Définition: Un "programme" est la réalisation (l'implémentation) d'un algorithme au moyen d'un langage donné (sur une architecture donnée). Il s'agit de la mise en oeuvre du principe.

Axiomes de la programmation (anecdotique) :

A1. Plus nous écrivons de code, plus nous produirons d'erreurs

A2. Il n'existe pas de programmes sans de possibles erreurs (du au programme lui-même, à l'électronique sous-jacente ou le plus souvent à l'utilisateur même)

Basiquement voici la démarche minimale à suivre lors du développement d'un produit informatique (au niveau du code):

M1. Auditer les besoins présents et anticiper les besoins futurs

M2. Définir les objectifs

M3. Calculer la faisabilité, le risque d'erreur, la durée nécessaire au développement

M4. Créer les algortithmes en langage formel (cela comprend la gestion des erreurs!)

M5. Optimiser la complexité et contrôler les algorithmes

M6. Choix d'une stratégie de développement (modulable ou autre)

M7. Traduire les algorithmes dans la technologie choisie (ce choix est un sujet assez sensible...)

Remarque: Il faudrait dans l'étape (7) respecter les normes de nommage et de représentation du code ainsi que les conventions (traditions) de fréquence d'insertions des commentaires.

M7. Tests (maintenance préventive)

Remarque: Le déboguage (la gestion des erreurs) d'un programme et les tests de fonctionnement doivent prendre autant, voir plus, de temps que le développement du programme lui-même.

Lors du développement d'un programme scientifique, il peut être intéressant, voire même rigoureux de s'attarder la notion de "complexité" de ce dernier. Sans aller trop loin, voyons un peu de quoi il s'agit :

COMPLEXITÉ

La complexité d'un algorithme est la mesure du nombre d'opérations fondamentales qu'il effectue sur une jeu de données (le nombre de fois qu'est exécuté une ligne de code). La complexité est donc exprimée comme une fonction de la taille du jeu de données (du nombre de lignes de code).

Les hypothèses permettant un calcul de cette complexité sont :

H1. (les quatre opérations fondamentales ont le même coût)

H2. Un accès mémoire une opération arithmétique

H3. Un contrôle de comparaison une opération arithmétique

H4. Un seul processeur

Nous notons  l'ensemble des données de taille n et  le coût (en temps) de l'algorithme sur la donnée ou le jeu de données de d de taille n.

- La "complexité au meilleur" est donnée par la fonction:

  (1)

C'est le plus petit temps qu'aura exécuté l'algorithme sur un jeu de données (de lignes de code) de taille fixée, ici à n dont le coût (la durée) d'exécution est . C'est une borne inférieure de la complexité de l'algorithme sur un jeu de données de taille n.

- La "complexité au pire":

  (2)

C'est le plus grand temps qu'aura à exécuté l'algorithme sur un jeu de données de taille fixée, ici à n. Il s'agit donc d'un maximum, et l'algorithme finira toujours avant d'avoir effectué  opérations. Cependant, cette complexité peut ne pas refléter le comportement usuel de l'algorithme, le pire cas pouvant ne se produire que très rarement, mais il n'est pas rare que le cas moyen soit aussi mauvais que le pire cas.

- La "complexité moyenne":

  (3)

Il s'agit de la moyenne des complexités de l'algorithme sur des jeux de données de taille n (en toute rigueur, il faut bien évidemment tenir compte de la probabilité d'apparition de chacun des jeux de données). Cette moyenne reflète le comportement général de l'algorithme si les cas extrêmes sont rares ou si la complexité varie peu en fonction des données. Cependant, la complexité en pratique sur un jeu de données particulier peut être nettement plus important que la complexité moyenne, dans ce cas la complexité moyenne ne donnera pas une bonne indication du comportement de l'algorithme.

En pratique, nous ne nous intéressons qu'à la complexité au pire, et à la complexité moyenne

Définition: Un algorithme est dit "algorithme optimal" si sa complexité est la complexité minimale parmi les algorithmes de sa classe.

Comme nous l'avons fait sous-entendre précédemment, nous nous intéressons quasi-exclusivement à la complexité en temps des algorithmes. Il est parfois intéressant de s'intéresser à d'autres de leurs caractéristiques, comme la complexité en espace (taille de l'espace mémoire utilisé), la largeur de bande passante requise, etc.

Pour que le résultat de l'analyse d'un algorithme soit pertinent, il faut avoir un modèle de la machine sur laquelle l'algorithme sera implémenté (sous forme de programme). Nous prenons habituellement comme référence, la machine à accès aléatoire (RAM)" et à processeur unique, où les instructions sont exécutées l'une après l'autre, sans opérations simultanées et sans processus stochastiques (contrairement au possibles machines quantiques futures).

Les algorithmes usuels peuvent êtres classés en un certain nombre de grandes classes de complexité dont l'ordre O varie d'une certain manière:

- Les algorithmes sub-linéaires dont la complexité est en général de l'ordre

- Les algorithmes linéaires en complexité  et ceux en complexité en  sont considérés comme rapides.

- Les algorithmes polynomiaux en  pour  sont considérés comme lents, sans parler des algorithmes exponentiels (dont la complexité est supérieur à tout en n) que l'on s'accorde à dire impraticables dès que la taille des données est supérieur à quelques dizaines d'unités.

Remarque: Une bonne complexité est de type pour . Une mauvaise complexité est de type ou .

Exemples:

E1. Evaluation d'un polynôme :

  (4)

L'évaluation directe de la valeur de conduit à une complexité

  (5)

Nous devons à Horner un algorithme plus performant qui utilise une factorisation du polynôme sous la forme :

  (6)

Nous pouvons facilement montrer que cette factorisation maintient le même nombre d'additions . La complexité qui en découle est donc . Le gain est incontestablement important. De plus, cette factorisation évite tout calcul de puissance. mais réduit le nombre de multiplication à

E2. Soit A et B deux matrices carrées de dimensiosn n. Les principales opérations sur ces matrices ont les complexités suivantes (nous laisserons le soin au lecteur de vérifier car c'est vraiment trivial) :

- Lecture (itérations) :

- Calcul de la trace :

- Addition : tel que

- Multiplication : tel que

- Déterminant (par la méthode directe de Cramer). Nous renvoyons le lecteur au chapitre d'algèbre linéaire pour le détail des méthodes de calcul du déterminant d'un matrice. Nous pouvons alors montrer que la complexité d'un déterminant d'ordre n est n multiplications, n-1 additions plus n fois la complexité d'un déterminant d'ordre n-1. Par cumul, nous arrivons à :

En faisant l'hypothèse que l'ordinateur utilisé effectue une opération élémentaire en secondes (ce qui est déjà une bonne machine). Nous obtenons alors les temps de calculs suivants pour plusieurs valeurs de n :

n	5	10	15	20	50
t

  (7)

d'où la nécessité de faire un calcul de complexité avant de mettre l'algorithme en route (à moins que vous ne travailliez exclusivement pour les générations futures, à condition qu'il y en ait encore…).

Finalement, pour résumer un peu, nous distinguons quelques types de complexités classiques : O(1) indépendant de la taille de la donnée, O(log(n)), complexité logarithmique, O(n) complexité linéaire, O(nlog(n)), complexité quasi-linéaire, O(n²), complexité quadratique, O(n³), complexité cubique, O(n³) complexité exponentielle, O(n!), complexité factorielle

NP-COMPLÉTUDE

Nous allons introduire pour la culture générale le concept de NP-complétude, c'est-à-dire que nous allons tenter de définir sans trop de formalisme (comme à l'habitude) la classe des problèmes dit "NP-complets". Ces problèmes sont ceux pour lesquels personne à l'heure actuelle ne connaît d'algorithme efficace (i.e. seulement des algorithmes exponentiels). Ce sont des problèmes vraiment difficiles par opposition à ceux pour lesquels nous connaissons des algorithmes de complexité polynomiale.

Définitions:

D1. Les problèmes de "classe P" sont des bons problèmes dans le sens où le calcul de leur solution est faisable dans un temps raisonnabl avec un algorithme à complexité polynomiale. Plus formellement, ce sont les problèmes pour lesquels nous pouvons construire une machine déterministes (voir la remarque après les définitions) dont le temps d'exécution est de complexité polynomiale (le sigle "P" signifiant "Polynomial Time").

D2. Les problèmes de "classe E" sont des problèmes dans le sens où le calcul de leur solution est faisable dans un temps expontentiel par nature de type . Plus formellement, ce sont les problèmes pour lesquels nous pouvons construire une machine déterministes dont le temps d'exécution est de complexité exponentielle (le sigle "E" signifiant "Expontential Time").

D3. Les problèmes de la "classe NP" sont ceux pour lesquels nous pouvons construire une machine de Turing non déterministe (voir la remarque après les définitions) dont le temps d'exécution est de complexité polynomiale (le sigle "NP" provient de "Nondeterministic Polynomial time" et non de "Non Polynomial").

Remarque: Contrairement aux machines déterministes qui exécutent une séquence d'instructions bien déterminée, les machines non déterministes ont la remarquable capacité de toujours choisir la meilleur séquence d'instructions qui mène à la bonne réponse lorsque celle-ci existe. Il va sans dire qu'une telle machine ne peut pas exister autrement que dans l'esprit d'un théoricien (à moins qu'avec les ordinateurs quantiques…)! Néanmoins, comme nous le verrons par la suite, ce concept abstrait n'est pas sans intérêt et constitue en fait la base de toute la théorie de la NP-complétude.

Il importe de remarquer à ce stade de la discussion que la classe P est incluse dans la classe NP, nous écrivons , car si nous pouvons construire une machine déterministe pour résoudre efficacement (en un temps au pire polynomial) un problème, nous pouvons certainement (du moins dans notre esprit) en construire une non déterministe qui résout aussi efficacement le même problème. Par ailleurs, il ne faut pas croire que ce concept de divination optimale qu'est la machine déterministe permet de tout résoudre puisqu'il existe en informatique théorique d'autres types de problèmes n'appartenant pas à la classe NP qui sont les problèmes indécidables.

Pour savoir si un problème donné appartient ou non à la classe NP, il s'agit simplement de l'exprimer sous la forme dont la réponse est soit OUI, soit NON. Le problème appartient alors à la classe NP si par définition, nous arrivons à démontrer à l'aide d'un algorithme de complexité polynomiale que n'importe quelle instance OUI du problème est bel et bien correcte. Nous n'avons pas à nous préoccuper des instances NON du problème puisque la machine non déterministe, par définition, prend toujours la bonne décision (lorsque celle-ci existe).

Par exemple, le problème consistant à trouver un cycle hamiltonien (cycle qui passe une et une seule fois par tous les sommets du graphe - voir chapitre de théorie des graphes) dans un graphe appartient à NP puisque, étant donné un cycle, il est trivial de vérifier en temps linéaire qu'il contient bien une et une seule fois chaque sommet.

Un autre exemple de problème difficile des mathématiques est la factorisation d'un entier en produit de facteurs premiers. Nous ne savons pas à ce jour s'il existe un algorithme polynomial qui réussisse cette opération. Autrement dit, nous ne savons pas si ce problème est dans P. En revanche, étant donné des nombres premiers il est trivial de vérifier que : ce problème est donc dans NP. Il semblerait (nous n'avons pas vérifié ce résultat et ni la possibilité de le faire) que la complexité du meilleur algorithme de factorisation en nombre premier soit du type :

  (8)

il resterait donc du travail à faire (si quelque'un peut nous fournir les détails qui ont amené à ce résultat, nous sommes preneurs).

Remarque: Si un problème est dans NP, alors il existera un algorithme en temps exponentiel pour le résoudre mais le contraire n'est pas toujours vrai (il faut donc être prudent).

Parmi l'ensemble des problèmes NP, il en existe un sous-ensemble qui contient les problèmes les plus difficiles : nous les appelons les problèmes "NP-complet" N.P.C. Ainsi, un problème NP-complet possède la propriété que tout problème dans NP peut être transformé en celui-ci en temps polynomial.

La raison essentielle pour laquelle les problèmes NPC sont les plus difficiles parmi les problèmes de NP est que ces premiers peuvent toujours servir comme des sous-routines pour solutionner ces derniers. Cette réduction à un ou des sous-routines assez facilement faisable puisque réalisable, si elle existe, en temps polynomial. Un problème NPC est donc complet en ce sens qu'il contient l'essentiel de la complexité des problèmes appartenant à NP, et qu'une solution polynomiale à ce problème implique une solution polynomiale à tous les problèmes NP.

Autrement formulé, les problèmes NPC ont un complexité exponentielle et ils ont tous la même classe de complexité (modulo les polynômes).

Finalement, ce qu'il importe de bien comprendre et de retenir de tout cette théorie, son idée maîtresse, est que si nous trouvons un jour un algorithme de complexité polynomiale pour un seul de ces problèmes vraiment difficiles que sont les problèmes NPC, alors d'un seul coup NP devient égal à P et tous les problèmes difficiles deviennent faciles !

Pour résumer, être dans P, c'est trouver une solution en un temps polynomial, tandis qu'être dans NP, c'est prouver en un temps polynomial qu'une proposition de réponse est une solution du problème. Ainsi, tout problème qui est dans P se trouve dans NP. Un champ de recherche majeur des mathématiques actuelles est l'investigation de la réciproque : a-t-on P=NP? Autrement dit, peut-on trouver en un temps polynomial ce que l'on peut prouver en temps polynomial?

Remarque: Ce problème est si important en informatique qu'il fait partie (arbitrairement) des 7 problèmes du millénaire, dont la résolution est primée 1 million de dollars par le Clay Mathematic Institute.

Passons maintenant à l'étude de quelques applications types des méthodes numériques dont il est très souvent fait usage dans l'industrie. Nous irons du plus simpleau plus compliqué et sans oublier que beaucoup de méthodes ne se trouvant pas dans ce chapitre peuvent parfois être trouvées dans d'autres section du site!

PARTIE ENTIÈRE

Le plus grand entier inférieur ou égal à un nombre réel x est par , qui se lit "partie entière de x".

Ainsi, le nombre M est entier si et seulement si [M]=M. De même, le naturel A est divisible dans l'ensemble des naturels par le naturel B si et seulement si:

  (9)

Nous notons aussi {x} pour désigner la partie fraction de x; on a ainsi:

  (10)

Soit . Alors nous avons les propriétés suivantes :

P1. , ,  où

P2. , lorsque

P3. , si

P4.

P5.  si ,  si

P6. si

P7. Si , alors [x / a] représente le nombre d'entiers positifs inférieurs ou égaux à x qui sont divisibles par a.

Démonstrations:

La première partie de P1 est simplement la définition de [x] sous forme algébrique. Les deux autres parties sont des réarrangements de la première partie. Dans ce cas, nous pouvons écrire

  (11)

où .

Pour P2, la somme est vide pour  et, dans ce cas, on adopte la convention selon laquelle la somme vaut 0. Alors, pour , la somme compte le nombre d'entiers positifs n qui sont plus petits ou égaux à x. Ce nombre est évidemment [x].

La démonstration de P3 sera supposée évidente.

Pour prouver P4, nous écrivons :

,   (12)

où n et m sont des entiers et où  et . Alors:

  (13)

En écrivant , où , nous avons

  (14)

où .

Il s'ensuit que:

0 si -1 si   (15)

et on obtient la démonstration P5.

Pour démontrer P6, nous éccrivons :

  (16)

où , et :

  (17)

où . Nous obtenons ainsi:

  (18)

puisque . Par ailleurs:

  (19)

et nous avons ainsi le résultat.

Pour la dernière partie, nous observons que, si  sont tous les entiers positifs  qui sont divisibles par a, il suffit de prouver que . Puisque  , alors:

  (20)

c'est-à-dire :

  (21)

soit le résultat attendu.

C.Q.F.D.

Remarque: La méthode d'arrondi de valeurs réelles est donnée dans le chapitre d'économétrie.

ALGORITHME D'HÉRON

Soit à calculer la racine carrée:

  (22)

Il existe un algorithme dit "algorithme d'Héron" qui permet de calculer la valeur de cette racine carrée.

Démonstration:

  (23)

Nous obtenons alors la relation:

  (24)

C.Q.F.D.

Exemple:

Soit à calculer :

  (25)

Nous prenons :

Itération	xi /2	A/2xi	xi+1	Ecart
1	5	0.5	5.50	~2.3
2	2.750	0.90	3.659 090 909	~0.49
3	1.82954	1.3664	3.196 005 083	~0.033
4	1.59800	1.5644	3.162 455 624	~0.0002
5	1.58122	1.5810	3.162 277 665	~0.5 10-8

Dans le cas de la racine cubique, la démonstration est semblable et nous obtenons:

  (27)

ALGORITHME D'ARCHIMÈDE

Le calcul de la constante universelle "pi" notée  est très certainement le plus grand intérêt de l'algorithmique puisque l'on retrouve cette constante un peu partout en physique et mathématique (nous pouvons vous conseiller un très bon ouvrage sur le sujet).

Nous rappelons que nous n'en avons pas donné la valeur ni en géométrie, ni dans les autres sections de ce site jusqu'à maintenant. Nous allons donc nous attacher à cette tâche.

Nous définissons définit en géométrie le nombre   dit "pi", quelque soit le système métrique utilisé (ce qui fait son universalité), comme le rapport de la moitié du périmètre d'un cercle par son rayon tel que:

  (28)

Nous devons le premier algorithme du calcul de cette constante à Archimède (287-212 av. J.-C.) dont voici la démonstration.

Démonstration:

Soit un n-polygone inscrit dans un cercle:

  (29)

Le principe de l’algorithme d’Archimède est le suivant: Soit le périmètre d’un polygone régulier de n côtés inscrit dans un cercle de rayon 1/2. Archimède arrive par induction que:

  (30)

Nous avons pour périmètre d'un n-polygone:

 et   (31)

Avec :

  (32)

Donc:

  (33)

Il suffit d'un ordinateur ensuite et de plusieurs itérations pour évaluer avec une bonne précision la valeur de .  Evidemment, on utilise l'algorithme d'Héron pour calculer la racine…

C.Q.F.D.

Remarque: Il existe un très grand nombre d'algorithmes pur calculer . Celle présentée ci-dessus, sans être la plus esthétique, est historiquement la première.

CALCUL DU NOMBRE D'EULER

Parmi la constante , il existe d'autres constantes mathématiques importantes qu'il faut pouvoir générer à l'ordinateur (de nos jours les valeurs constantes sont stockées telles quelles et ne sont plus recalculées systématiquement). Parmi celles-ci, se trouve le "nombre d'Euler" noté e (cf. chapitre d'Analyse Fonctionnelle). Voyons comment calculer ce nombre:

Soit la série de Taylor, pour une fonction indéfiniment dérivable f donnée par (cf chapitre sur les Suites Et Séries):

  (34)

Comme (cf. chapitre de Calcul Différentiel Et Intégral):

  (35)

nous avons:

  (36)

Donc en résumé:

  (37)

Cette relation donne un algorithme pour calculer l'exponentielle  à un ordre n donnée de précision.

Remarque: Pour diminuer la complexité de cet algorithme, la factorielle peut être calculée avec la formule exposée ci-après.

CALCUL DE LA FACTORIELLE (FORMULE DE STIRLING)

Evidemment, la factorielle pourrait être calculée avec une simple itération. Cependant, ce genre de méthode génére un algorithme à complexité exponentielle ce qui n'est pas le mieux. Il existe alors une autre méthode :

Soit, la définition de la factorielle :

  (38)

Et d'après les propriétés des logarithmes :

  (39)

Si n est très grand (mais alors très…) alors:

  (40)

Lorsque , la limite inférieure est négligeable et alors:

  (41)

Après une petite simplification élémentaire, nous obtenons:

  (42)

Cette dernière relation est utile si l'on suppose bien évidemment que la constante d'euler est une valeur stockée dans la machine...

SYSTEMES D'ÉQUATIONS LINÉAIRES

Il existe de nombreuses méthodes de résolution de systèmes d'équations linéaires. La plupart d'entre elles ont été mises au point pour traiter des systèmes particuliers. Nous en étudierons une, appelée la "méthode de réduction de Gauss", qui est bien adaptée à la résolution des petits systèmes d'équations (jusqu'à 50 inconnues).

Remarques:

R1. La validité de certaines des opérations que nous allons effectuer ici pour résoudre les systèmes linéaires se démontrent dans le chapitre traitant de l'algèbre linéaire. Au fait, pour être bref, le tout fait appel à des espaces vectoriels dont les vecteurs-colonnes sont linéairement indépendants.

R2. Les systèmes admettent une solution si et seulement si (rappel) le rang de la matrice augmentée est inférieur ou égal au nombre d'équations (cf. chapitre d'Algèbre Linéaire).

UNE ÉQUATION À UNE INCONNUE

Considérons le cas le plus simple: celui d'une équation à une inconnue:

  (43)

a et b sont les coefficients de l'équation et b en est l'inconnue. Résoudre cette équation, c'est déterminer x en fonction de a et b. Si a est différent de 0 alors:

  (44)

est la solution de l'équations. Si a est nul et si b est différent de 0 alors l'équation n'admet pas de solutions. Si a et b sont nuls, alors l'équation possède une infinité de solutions.

DEUX ÉQUATIONS À DEUX INCONNUES

Un système (linéaire) de deux équations à deux inconnues s'écrit:

  (45)

 sont les coefficients du système d'équations,  et  en sont les inconnues.

Remarque: Les notations usitées ci-dessus n'ont rien à voir avec le calcul tensoriel.

Pour résoudre le système, nous procédons comme suit:

A l'aide de manipulations algébriques (addition ou soustraction des différentes égalités entre elles - manipulations autorisées par l'indépendance linéaire des vecteurs-lignes) nous transformons le système en un autre donné par :

  (46)

Ensuite, nous résolvons l'équation , ce qui donne :

  (47)

Nous en déduisons:

  (48)

La transformation entre les deux systèmes:

  (49)

s'effectue simplement en multipliant chaque coefficient de la première égalité par  et en soustrayant les résultats obtenus des coefficients correspondants de la seconde égalité. Dans ce cas, l'élément  est appelé "pivot".

TROIS ÉQUATIONS À TROIS INCONNUES

Examinons encore le cas des systèmes de trois équations à trois inconnues:

  (50)

Nous pouvons par la suite des opérations élémentaires (cf. chapitre d'Algèbre Linéaire) réduire le système linéaire précédent en le système suivant:

  (51)

en ensuite résoudre l'équation:

  (52)

puis la deuxième:

  (53)

et enfin:

  (54)

Revenons à la transformation des systèmes. Elle s'effectue en deux étapes :

1. Dans la première, nous choisissons  comme pivot et nous éliminons les coefficients  et  de la manière suivante: il faut multiplier chaque coefficient de la première égalité par  et soustraire les résultats obtenus de la deuxième égalité, ainsi  devient nul. De même, en multipliant les coefficients de la première équation par  et en soustrayant les résultats obtenus de la troisième égalité,  disparaît. Le système d'équation s'écrivant alors:

  (55)

2. La deuxième étape consiste à traiter le système de deux équations à deux inconnues formé des deuxième et troisième égalités, et ce, en choisissant  comme pivot. Cette méthode de résolution peut paraître compliquée, mais elle a l'avantage de pouvoir être généralisée et être appliquée à la résolution de systèmes de n équations n inconnues.

N ÉQUATIONS À N INCONNUES

Pour simplifier l'écriture, les coefficients seront toujours notés  et non pas , etc. lors de chaque étape du calcul.

Soit le système linéaire (nous pourrions très bien le représenter sour la forme d'une matrice augmentée afin d'alléger les écritures) :

  (56)

Il faut choisir  comme pivot pour éliminer . Ensuite, l'élimination de  s'effectue en prenant  comme pivot. Le dernier pivot à considérer est bien évidemment , il permet d'éliminer . Le système prend alors la forme:

  (57)

En résolvant d'abord la dernière équation, puis l'avant dernière et ainsi de suite jusqu'à la première.

Cette méthode, appelée "méthode de résolution de Gauss" ou encore "méthode du pivot" doit cependant être affinée pour éviter les pivots de valeur 0. L'astuce consiste à permuter l'ordre dans lequel sont écrites les équations pour choisir comme pivot le coefficient dont la valeur absolue est la plus grande. Ainsi, dans la première colonne, le meilleur pivot est l'élément  tel que:

  (58)

Il est amené à  par échange des première et j-ème lignes. L'élimination du reste de la première colonne peut alors être effectuée. Ensuite, on recommence avec les n-1 équations qui restent.

POLYNÔMES

L'ensemble des polynômes et à coefficients réels a été étudié dans le chapitre d'analyse fonctionnelle en détails.. Nous allons ici traiter de l'aspect numérique de quelques problèmes liés aux polynômes.

Mis à part l'addition et la soustraction de polynômes que nous supposerons comme triviaux (optimisation de la complexité mis à part comme le schéma de Horner), nous allons voir comment multiplier et diviser deux polynômes.

Voyons d'abord comment multiplier deux polynômes :

Soit :

  (59)

alors :

  (60)

avec :

  (61)

C'était simple….

Un tout petit peu plus difficile maintenant : la division euclidienne des polynômes (cf. chapitre de Calcul Algébrique).

Reprenons :

  (62)

mais en imposant cette fois-ci .

La division s'écrira donc nous le savons :

  (63)

avec ou sinon .

Il est connu d'avance que nous avons bien évidemment :

et deg(r(x))<m.

Nous avons donc par définition q(x) qui est le quotient de la division et r(x) le reste de la division euclidienne de f(x) par g(x).

Rien ne nous interdit donc de poser de la manière la plus générale qui soit :

  (64)

Exemple:

Soit :

  (65)

donc :

  (66)

Nous avons donc :

  (67)

Ensuite :

  (68)

et enfin :

  (69)

 

Donc de manière générale :

  (70)

comme :

  (71)

le premier reste est donc :

  (72)

Ensuite :

  (73)

Le deuxième reste est alors :

  (74)

etc.

Nous continuons jusqu'à ce que .

RÉGRESSIONs ET INTERPOLATIONS

Les régressions (ou "interpolations") sont des outils très utiles aux statisticiens, ingénieurs, informaticiens souhaitant établir un loi de corrélation entre deux (ou plus) variables dans un contexte d'études et d'analyse ou d'extrapolation.

Il existe un grand nombre de méthodes d'interpolation : de la simple résolution d'équations du premier degré (lorsque uniquement deux points d'un mesure sont connus) aux équations permettant d'obtenir à partir d'un grand nombre de points des informations essentielles à l'établissement d'une loi (ou fonction) de régression linéaire, polynomiale ou encore logistique.

RÉGRESSION LINÉAIRE

Nous présentons ici deux algorithmes (méthodes) utiles et connus dans les sciences expérimentales (nous en avons déjà parlé lors de notre étude des statistiques). L'objectif ici est de chercher à exprimer la relation linéaire entre deux variables x et y indépendantes :

- x est la variable indépendante ou "explicative". Les valeurs de x sont fixées par l'expérimentateur et sont supposées connues sans erreur

- y est la variable dépendante ou "expliquée" (exemple : réponse de l'analyseur). Les valeurs de y sont entachées d'une erreur de mesure. L'un des buts de la régression sera précisément d'estimer cette erreur.

Nous cherchons une relation de la forme:

  (75)

C'est l'équation d'une droite, d'où le terme de "régression linéaire".

DROITE DE RÉGRESSION

Il existe aussi une autre manière commune de faire une régression linéaire du type :

  (76)

qui consiste à se baser sur les propriétés de la covariance et de l'espérance (cf. chapitre de Statistiques) et très uitilisée en finance.

Soit x, y deux variables dont l'une dépend de l'autre (souvent c'est y qui dépend de x) nous avons selon la propriété de linéarité de la covariance qui est rappelons-le :

  (77)

la relation suivante :

  (78)

Il vient donc pour la pente (nous réutiliserons cette relation lors de l'étude du rendement d'un portefeuille selon le modèle de Sharpe dans le chapitre d'économétrie) :

  (79)

et pour l'ordonnée à l'origine nous utilisons les propriétés de l'espérance démontrées dans le chapitre de statistiques :

  (80)

ce qui donne b sous la forme :

  (81)

Remarque: C'est la méthode utilisée par MS Excel lors de l'utilisation de la fonction PREVISION().

MÉTHODE DES MOINDRES CARRÉS

Du fait de l'erreur sur y, les points expérimentaux, de coordonnées , ne se situent pas exactement sur la droite théorique. Il faut donc trouver l'équation de la droite expirementale qui passe le plus près possible de ces points.

La "méthode des moindres carrés" consiste à chercher les valeurs des paramètres a, b qui rendent minimale la somme des carrés des écarts résiduels (: Sum of Squared Residuals ) entre les valeurs observées  et les valeurs calculées théoriques de :

  (82)

où n est le nombre de points et:

  (83)

d'où:

  (84)

Cette relation fait apparaître la somme des carrés des écarts comme une fonction des paramètres a,b. Lorsque cette fonction est minimale, les dérivées par rapport à ses paramètres s'annulent:

  (85)

Remarque: Cette méthode de recherche de minimum (optimisation) est nommée "méthode des multiplicateurs de Lagrange" dans le monde de l'économétrie. Dans notre exemple est la grandeur scalaire qui fait office de multiptlicateur de Lagrange.

Soit après simplification :

  (86)

Le système ci-dessus est dit appelé "système des équations normales". Il admet pour solutions (après quelques substitutions et simplifications triviales) l'expression de la pente et de l'ordonnée à l'origine de l'équation recherchée :

  (87)

Remarque: C'est la méthode utilisée par MS Excel lors de l'utilisation de la fonction REGRESSION().

RÉGRESSION LOGISTIQUE

Bien souvent, les données statistiques disponibles sont relatives à des caractères qualitatifs. Or, comme nous allons le voir, les méthodes d'inférence traditionnelles ne permettent pas de modéliser et d'étudier ce type caractères. Des méthodes spécifiques doivent être utilisées tenant compte par exemple de l'absence de continuité des variables traitées ou de l'absence d'ordre naturel entre les modalités que peut prendre la caractère qualitatif. Ce sont ces méthodes spécifiques les plus usuelles qui seront l'objet du texte qui suit.

Comme nous l'avons vu plus haut, la régression linéaire simple a donc pour but de modéliser la relation entre une variable dépendante quantitative et une variable explicative quantitative.

Lorsque la "variable de classe" Y à expliquer est binaire (oui-non, présence-absence,0-1,etc.) nous approchons dans un premier temps celle-ci par une fonction de probabilité , qui nous donne à l'opposé la probabilité d'appartenir à la classe , que nous nommerons "régression logistique" ou encore "régression logit" (très souvent utilisée dans les cadre des réseaux de neurones formels). Ensuite, dans une deuxième étape, nous définissons pour un cas binaire une valeur "cutoff". Par exemple, si nous prenons un cutoff de 0.5 alors les cas pour lesquels  appartiendront à la classe 1 (et inversement dans le cas contraire).

Remarques:

R1. Au fait, la régression logistique n'est qu'une simple loi de distribution de probabilités.

R2. Il n'est évidemment pas possible d'appliquer systématiquement la régression logistique à n'importe quel type d'échantillon de données! Parfois il faut chercher ailleurs…

Nous utiliserons ce modèle, quand il est adapté au contexte, au lieu de la régression linéaire. Effectivement, si nous utilisons une régression linéaire les valeurs prédictives peuvent devenir négatives ou supérieures à 1… ce qui est gênant pour une variable dichotomique…

Remarque: Lorsque le nombre de modalités est égal à 2, nous parlons de "variable dichotomique" (oui-non) ou d'un "modèle dichotomique", s'il est supérieur à 2, nous parlons de "variables polytomiques" (satisfait-non satisfait-émerveillé).

Considérons par exemple la variable dichotomique : fin des études. Celle-ci prend deux modalités (a fini,a poursuivi). L'âge est une variable explicative de cette variable et nous cherchons à modéliser la probabilité d'avoir terminé ses études en fonction de l'âge.

Exemple:

Pour construire le graphique suivant, nous avons calculé et représenté en ordonnées, pour des jeunes d'âge différent, le pourcentage de ceux qui ont arrêté leurs études.

  (88)

Mais comment obtient-t-on pareil graphique avec une variable dichotomique??? Au fait c'est simple. Imaginez un échantillon de 100 individus. Pour ces 100 individus supposez pour un âge donné que 70% "a fini" et 30% "a terminé". Eh bien la courbe représente simplement la proportion des deux classes pour l'âge donné. Il est même parfois indiquée les grosseurs des classes avec des cercles sur tout la longueur des asymptotes horizontales pour bien signifier qu'il s'agit d'une variable dichotomique.

Les points sont distribués selon une courbe en S (une sigmoïde) : il y a deux asymptotes horizontales car la proportion est comprise entre 0 et 1. Nous voyons imédiatement qu'un modèle linéaire serait manifestement inadapté.

Cette courbe évoqueront pour certains, à juste titre, une courbe cumulative représentant une fonction de répartition (d'une loi normale par exemple, mais d'autres lois continues ont la même allure). Ainsi, pour ajuster une courbe à cette représentative, nous pourrions nous orienter vers la fonction de répartition d'une loi normale, et au lieu d'estimer les paramètres a et b de la régression linéaire, nous pourrions estimer les paramètres  de la loi normale (qui est très similaire à la loi logistique comme nous le démontrerons plus loin). Nous parlons alors d'un "modèle Probit".

La loi qui nous intéresse cependant est donc la loi logistique. Contrairement à la loi normale, nous savons calculer l'expression de sa fonction de répartition dichotomique (probabilité cumulée) qui est du type (c'est son premier avantage!):

  (89)

pour une variable de prédiction x où P est donc la probabilité d'avoir un 1.

S'il y a plusieurs variable prédictives nous avons alors :

  (90)

Lorsque nous optons pour cette fonction de répartition de la loi logistique, nous obtenons le donc modèle de régression logistique ou "modèle Logit".

Ainsi, nous estimerons la probabilité cumulée d'avoir fini ses études pour un individu d'âge x par (il existe plusieurs manières d'écrire cette loi suivant les habitudes et le contexte) :

  (91)

il en découle la fonction de distribution :

  (92)

Nous pouvons calculer aussi l'espérance de la fonction de distribution en appliquant ce qui a déjà été vu au chapitre de statistique mais une partie de cette intégrale ne peut être résolue que numériquement par contre… :

  (93)

et :

  (94)

Ainsi, nous voyons que si nous posions :

  (95)

Nous retombons sur une fonction de répartition ayant parfaitement les mêmes caractéristiques qu'une loi normale (moyenne nulle et variance unitaire).

Exemple:

La fonction sigmoïde (de répartition) est présentée ci-dessous pour  :

  (96)

Les paramètres a, b sont ajustés selon le principe du maximum de vraisemblance (cf. chapitre de Statistiques). De plus, ces paramètres doivent généralement être ajustés de manière itérative, à l'aide d'un programme auquel nous fournissons des valeurs initiales, et qui optimise ces valeurs de manière récurrente (nous n'entrerons pas dans ces détails qui dépassent le cadre théorique de ce site à ce jour).

La dernière relation peut par ailleurs être transformée de la façon suivante :

  (97)

d'où :

  (98)

Ce que certains écrivent aussi :

  (99)

Le résultat de cette dernière transformation est appelé "logit". Il est égal au logarithme de "l'odds" P/1-P.

Donc lorsque les coefficients a et b ont été déterminés, l'expression précédente permet de déterminer P connaissant x facilement (il s'agit de résoudre une équation linéaire). Par ailleurs, puisque x est une variable dichotomique les coefficients sont très facilement interprétables.

Remarque: L'odds est également appelé "cote" par analogie à la cote des chevaux au tiercé. Par exemple, si un étudiant a 3 chances sur 4 d'être reçu, contre 1 chance sur 4 d'être collé, sa cote est de 3 contre un 1, soit un odds=3.

Revenons un peu sur l'odds car il est possible d'introduire la notion de fonction logistique en faisant la démarche inverse de celle présentée ci-dessus (soit de commencer par la définition de l'odds pour aller jusqu'au logit) et ceci peut parfois même s'avérer plus pédagogique.

Supposons que nous connaissons la taille (hauteur) d'une personne pour prédire si la personne est un homme ou une femme. Nous pouvons donc parler de probabilité d'être un homme ou une femme. Imaginons que la probabilité d'être un homme pour une hauteur donnée est 90%. Alors l'odds d'être un homme est :

  (100)

Dans notre exemple, l'odds sera de 0.90/0.10 soit. Maintenant, la probabilité d'être une femme sera donc de 0.10/0.90 soit 0.11. Cette asymétrie des valeurs est peu parlante parce que l'odds d'être un homme devrait être l'opposé de l'odds d'être une femme idéalement. Nous résolvons justement cette asymétrie à l'aide du logarithme naturel. Ainsi, nous avons :

 et    (101)

Ainsi, le logit (logarithme de l'odds) est exactement l'opposé de celui d'être une femme de par la propriété du logarithme:

  (102)

Exemple:

Imaginons qu'une banque souhaite faire un scoring des ses débiteurs. Comme elle a trois succursales (la banque) elle a les données (fictives…) suivantes à disposition :

- 1^ère succursale :

- 2^ème succursale :

- 3^ème succursale :

Nous pouvons voir que la proportion totale des bons débiteurs est de . Quand le crédit est inférieur à 27'500, la proportion de bons débiteurs est de . Quand le montant des crédits est inférieur à 28'000 la proportion de bons débiteurs est de . Quand le montant des crédits est inférieur à 28'500, la proportion de bons débiteurs est de , et pour des montants inférieurs à 30'000 la proportion est de .

Nous poserons pour la régression logistique que  est un bon risque de crédit et  est un mauvais risque. Ensuite, nous créons le tableau suivant qui est un récapitulatif des données des trois succursales :

Ce qui donne graphiquement en Kilo-Francs :

  (107)

Une fois ceci fait, nous utilisons la transformation en logit :

  (108)

Ce qui donne :

Une régression linéaire par la méthode des moindres carrés donne :

  (110)

avec pour équation :

  (111)

La fonction logistique avec sa représentation vient alors immédiatement :

  (112)

INTERPOLATION POLYNÔMIALE

Il existe de nombreuses techniques d'interpolation de pôlynomes plus ou moins complexes et élaborées. Nous nous proposons ici de présenter quelques unes dans l'odre croissant de difficulté et de puissance d'application.

COURBES DE BÉZIER (SPLINE)

L'ingénieur russe Pierre Bézier (Peugot), aux débuts de la Conception Assistée par Ordinateur (C.A.O), dans les années 60, a donné un moyen de définir des courbes et des surfaces à partir de points. Ceci permet la manipulation directe, géométrique, des courbes sans avoir à donner d'équation à la machine!!

Le thème des Courbes de Bézier est une notion à multiples facettes, vraiment très riche, au croisement de nombreux domaines mathématiques très divers : Analyse, Cinématique, Géométrie Différentielle, Géométrie Affine, Géométrie Projective, Géométrie Fractale, Probabilités, ...

Les Courbes de Bézier sont par ailleurs devenues incontournables dans leurs applications concrètes dans l'industrie, l'infographie, ...

Voilà l'approche mathématique de cette technique:

D'abord, nous savons que l'équation d'une droite que nous noterons dans le domaine (par respect de tradition) M joignant deux points  est:

  (113)

Ce qui est juste puisque lorsque  nous sommes en A et lorsque  nous sommes en B. Donc  et le point M parcoure tout le segment [AB]. Par construction, si A et B étaient des masses physique égales à l'unité,  représente le barycentre (centre de gravité) du système pour un t donné.

Par définition, le segment [AB] est par définition la "courbe de Bézier de degré 1" avec points de contrôle A et B et les Polynômes 1-t et t sont les "polynômes de Bernstein de degré 1".

Construisons maintenant une courbe paramétrée en rajoutant une 2ème étape à ce qui précède:

  (114)

 

1ère étape :

- Soit  le barycentre de (A,1-t) et (B, t) et où  décrit [AB].
- Soit  le barycentre de (B,1-t) (C, t) et où  décrit [BC].

2ème étape :

- Soit M(t) le barycentre de (,1-t) (,t).

Par construction, M(t) se situe donc à la même proportion du segment  que  par rapport au segment [AB] ou par rapport au segment [BC].

La courbe obtenue est alors l'enveloppe des segments : en tout point M, la tangente à la courbe est donc le segment .

M(t) décrit alors une Courbe de Bézier de degré 2, qui, par construction
commence en A et se finit en C, et a pour tangentes [AB] en A et [BC] en C.

C'est en fait un arc de parabole (que nous pourrions noter très logiquement [ABC] ) :

  (115)

Par le même schéma, nous pouvons définir une courbe de Bézier de n points  soit . C'est ce que nous appelons "l'algorithme de Casteljau". Ainsi, soit:

  (116)

Nous avons:

  (117)

La récurrence se terminant pour:

  (118)

Ainsi, pour  nous avons trivialement:

  (119)

Soit:

  (120)

Ainsi, nous avons forcément avec deux points l'équation d'une droite.

Considérons maintenant M(t) la courbe de Bézier d'ordre 3, nous avons donc les points définis par:

  (121)

Nous avons par la relation de récurrence:

  (122)

où nous avons éliminé les termes contenant des point non déterminés.

Nous avons donc:

  (123)

Il vient alors:

  (124)

Et donc:

  (125)

Soit sous forme matricielle:

  (126)

Par le même raisonnement, nous avons pour une courbe de Bézier d'ordre 4:

  (127)

Et là nous pouvons creuser un peu les coefficients en comparant les coefficients de Bézier des courbes d'ordre 2, 3 et 4.

D'abord, reprenons la courbe de Bézier précédente:

  (128)

Nous remarquons d'abord aisément la proportionnalité suivante:

  (129)

et si on regarde de plus près les coefficients nous remarquons que nous avons aussi:

  (130)

Il ne s'agit ni plus ni moins que du triangle de Pascal!! Et donc les constantes sont simplement les coefficients binomiaux (cf. chapitre de Calcul Algébrique) donnés par pour l'ordre n dans notre cas par:

  (131)

Ainsi, les polynômes de Bernstein sont donnés par:

  (132)

et finalement les courbes de Bernstein par d'ordre n :

  (133)

Au fait, si nous avions noté plus haut la somme sous la forme suivante:

  (134)

Nous aurions alors les polynômes de Bernstein qui seraient donnés (ce qui est plus respectueux des traditions….):

  (135)

C'est une relation très pratique car elle permet de calculer facilement et rapidement le polynôme correspondant à une courbe de Bézier d'ordre n.

Nous avons alors:

  (136)

Remarque: Une courbe de Bézier est totalement modifiée dès qu'on déplace un point de contrôle. Nous disons alors que que la méthode de Bézier est une "méthode globale".

Un exemple très connu des courbes de Bézier d'ordre 3 est l'outil Plume des produits phares Adobe Photoshop ou Adobe Illustrator. Effectivement ces deux outils créent une succession de courbes de Bézier d'ordre 3 dont le point  est défini après coup à la souris en utilisant des poignées ("torseurs" dans le langage de spécialistes Adobe…). Voici un exemple prix d'un de ces logiciels fait avec un tracé à la plume de 5 points (soit 4 splines):

  (137)

Tant que l'utilisateur ne bouge pas les poignées de points alors tous les points sont alignés sur la droite. Nous avons alors l'impression d'avoir une spline d'ordre 2.

MÉTHODE D'EULER

Il s'agit là de la méthode numérique la plus simple (elle est triviale et dans l'idée assez proche de la méthode de Newton même si leur objectif n'est pas le même). En fait elle ne fournit une approximation (au sens très large du terme) d'une fonction f(x) dont la dérivée est connue.

Ici les points choisis sont équidistants, c'est-à-dire (h étant le "pas"). Nous notons la valeur exacte et , la valeur approchée.

Il y a plusieurs méthodes pour procéder (comme souvent) :

1. Graphiquement :

Nous nous déplaçons d'un pas en en suivant le vecteur de pente f(x,y)

Par construction nous savons (cf. chapitre de Calcul Différentiel Et Intégral) que (nous adoptons une notation un peu particulière dans ce contexte) :

  (138)

qui correspond donc bien à la pente (non instantané bien sûr!) en de la "courbe intégrale" passant par ce point. D'où :

  (139)

2. Analytiquement :

Nous remplaçons dans la dernière relation par . Nous obtenons alors :

  (140)

appelée "équation aux différences pour la méthode d'Euler".

L'application en est triviale et ne nécessite pas d'exemples particuliers.

POLYNÔME DE COLLOCATION

Soit une fonction connue sous forme explicite ou sous forme tabulée, et supposons qu'un certain nombre de valeurs :

  (141)

en sont données. Les points sont appelés les "points d'appui".

"Interpoler" f signifie estimer les valeurs de f pour des abscisses situées entre et , c'est-à-dire dans l'intervalle d'interpolation, par une fonction approximante , qui vérifie les "conditions de collocations" (rien à voir avec votre colocataire !) :

  (142)

La fonction p s'appelle "fonction de collocation" sur les . Lorsque p est un polynôme, nous parlons de "polynôme de collocation" ou de "polynôme d'interpolation".

"Extrapoler" f signifie approcher f(x) par p(x) pour des abscisses situées "hors" de l'intervalle d'interpolation.

Remarque: Il va sans dire que l'interpolation est un outil très important pour tous les chercheurs, statisticiens et autres.

Quand nous connaissons un polynôme de degré n en n+1 points, nous pouvons donc connaître par une méthode simple (mais pas très rapide – mais il existe plusieurs méthodes) complètement ce polynôme.

Pour déterminer le polynôme, nous allons utiliser les résultats exposés précédemment lors de notre étude des système d'équations linéaires. Le désavantage de la méthode présentée ici est qu'il faut deviner à quel type de polynôme nous avons affaire et savoir quels sont les bons points qu'il faut choisir…

Un exemple particulier devrait suffire à la compréhension de cette méthode, la généralisation en étant assez simple (voir plus loin).

Soit un polynôme du second degré :

  (143)

et nous avons connaissances des points suivants (dont vous remarquerez l'ingéniosité des points choisis par les auteurs de ces lignes…) :

  (144)

Nous en déduisons donc le système d'équations :

  (145)

Système qui une fois résolu dans les règles de l'art nous donne :

  (146)

Voyons le cas général :

Théorème : Soient des points d'appui, avec si . Alors il existe un polynôme de degré inférieur ou égal à , et un seul, tel que pour .

Démonstration:

Posons :

  (147)

Les conditions de collocation :

  (148)

s'écrivent donc :

  (149)

Il s'agit d'un système de équations à inconnues.

Sont déterminant s'écrit :

  (150)

relation que nous appelons "déterminant de Vandermonde". Nous savons que si le système à une solution, le déterminant du système doit être non nul (voir chapitre d'algèbre linéaire).

Montrons par l'exemple (en reprenant un polynôme du même degré que celui que nous avons utilisé plus haut) que le déterminant se calcule selon la relation suivante précédente (le lecteur généralisera par récurrence) :

Donc dans le cas , nous considérons le déterminant :

  (151)

qui correspond dont au système (pour rappel) :

  (152)

Calculons ce déterminant suivant la colonne 1 (en faisant usage des cofacteurs comme démontré dans le chapitre d'lgèbre linéaire) :

  (153)

Ce dernier polynôme peut s'écrire :

  (154)

Ce qui s'écrit :

  (155)

Comme les sont l'énoncé de notre problème tous différents tels que alors le système a une solution unique. Ce qui prouve qu'il existe toujours alors un polynôme d'interpolation.

C.Q.F.D.

RECHERCHE DES RACINES

Bien des équations rencontrées en pratique ou en théorie ne peuvent pas être résolues exactement par des méthodes formelles ou analytiques. En conséquence, seule une solution numérique approchée peut être obtenue en un nombre fini d'opération.

Evartiste Galois a démontré, en particulier, que l'équation  ne possède pas de solution algébrique (sauf accident...) si  est un polynôme de degré supérieur à 4.

Il existe un grand nombre d'algorithmes permettant de calculer les racines de l'équation  avec une précision théorique arbitraire. Nous n'en verrons que les principaux. Attention, la mise en œuvre de tels algorithmes nécessite toujours une connaissance approximative de la valeur cherchée et celle du comportement de la fonction près de la racine. Un tableau des valeurs de la fonction et sa représentation graphique permettent souvent d'acquérir ces connaissances préliminaires.

Si l'équation à résoudre est mise sous la forme , nous traçons les courbes représentant g et h. Les racines de l'équation  étant données par les abscisses des points d'intersection des deux courbes.

Remarque: Avant de résoudre numériquement l'équation , il faut vérifier que la fonction f choisie. Il faut par exemple, que la fonction f soit strictement monotone au voisinage de la racine , lorsque la méthode de Newton est appliquée. Il est souvent utile, voire indispensable, de déterminer un intervalle  tel que:

- f est continue sur  ou

-

-  unique,

MÉTHODE DES PARTIES PROPORTIONNELLES

La mise en œuvre, sur calculatrice, de cette méthode est particulièrement simple. Les conditions à vérifier étant seulement:

- f est continue

- f est monotones dans un voisinage de la racine

Dans un petit intervalle, nous pouvons remplacer une courbe par un segment de droite. Il y a plusieurs situations possible mais en voici une particulière généralisable facilement à n'importe quoi:

  (156)

Sur cette figure nous tirons à l'aide des théorèmes de Thalès (cf. chapitre de Géométrie Euclidienne):

  (157)

d'où :

  (158)

Si , nous pouvons négliger f(a) au dénominateur et il vient :

  (159)

L'algorithme consiste donc à réaliser les étapes suivantes :

1. Choisir a et b, calculer f(a) et f(b) 

2. Déterminer . Si  est assez petit, nous arrêtons le calcul et affichons x et f(x)

3. Sinon nous procédons comme suit:

- Nous remplaçons b par a et f(b) par f(a) 

- nous remplaçons a par x et f(a)  par f(x)

- nous retournons au point (2)

MÉTHODE DE LA BISSECTION

La condition préalable à satisfaire pour cette méthode est de trouver un intervalle  tel que:

1.  f(x) est continue sur

2.

Il faut encore fixer  qui est définit comme la borne supérieure de l'erreur admissible.

La méthode consiste à appliquer successivement les 4 étapes suivantes:

1. Calcul de

2. Evaluation de  f(x)

3. Si  alors le travail est terminé, il faut afficher x et  f(x)

4. Sinon on procède comme suit:

- on remplace a par x si

- on remplace b par x si  ou

- on retourne en (1)

L'étape (3) impose la condition  pour l'arrêt des calculs. Il est parfois préférable de choisir un autre critère de fin de calcul. Celui-ci impose à la solution calculée d'être confinée dans un intervalle de longueur  contenant . Ce test s'énonce comme suit:

3'. Si , le travail est terminé et  est affiché. Il est bien sûr évident que

MÉTHODE DE LA SÉCANTE (REGULA FALSI)

  (160)

Les conditions préalables sont les suivantes:

Il faut déterminer un intervalle  tel que:

1.  f(x) est continue sur

2.

Si  est le point de coordonnées , alors les points  sont alignés sur la sécante. La proportion suivante (Thalès) est donc vraie:

  (161)

nous en déduisons:

  (162)

La méthode consiste à appliquer successivement les étapes suivantes:

1. Calcul de

2. Evaluation de

3. Si , le travail est terminé. Il faut afficher

4. Sinon on procède comme suit:

- on remplace  par  si

- on remplace  par  si  ou

- on retourne en (1)

La condition (3) peut être remplacée par la condition:

3'. Si , alors le travail est terminé et on affiche

Remarque: Si l'intervalle [a,b] contient plusieurs racines, cette méthode converge vers l'une d'entre elles; toutes les autres sont perdues.

MÉTHODE DE NEWTON

Considérons la figure suivante:

  (163)

Si  est une approximation de la racine , nous remarquons que  en est une meilleure.  est l'intersection de la tangente à la courbe en  et de l'axe des abscisses.  est encore une meilleure approximation de ,  est obtenu de la même manière que  mais à partir de .

Le méthode de Newton consiste en la formalisation de cette constatation géométrique.

Pour utiliser cette technique, rappelons que si nous prenons une fonction f qui est dérivable en , alors nous pouvons la réécrire sous la forme:

  (164)

où  est la dérivée de f en  et  est une fonction qui tend vers 0 comme  pour  lorsque x tend vers  (c'est un terme correctif qui sous-tend la suite des termes du développement de Taylor).

En appliquant ce résultat à la résolution de , nous obtenons:

  (165)

La fonction  empêche la résolution de cette équation par rapport à . En négligeant le terme , l'équation se réécrit:

  (166)

et se résout aisément par rapport à :

  (167)

Mais  ne satisfait pas, en générale, l'égalité . Mais comme nous l'avons déjà souligné,  est plus petit que  si la fonction f satisfait à certaines conditions.

La méthode de Newton consiste à remplacer l'équation:

  (168)

par:

  (169)

et à résoudre itérativement cette équation.

Les conditions suivantes sont suffisantes pour assurer la convergence de la méthode:

Dans un intervalle  comprenant  et  il faut que:

1. La fonction soit deux fois dérivable

2. La dérivée f ' ne s'annule pas (monotonie)

3. La deuxième dérivée soit continue et ne s'annule pas (pas de point d'inflexion)

Remarque: Il suffit souvent de vérifier les conditions (1) et (2) pour que le processus soit convergent.

La condition (2) est évident, en effet si  alors l'itération peut conduire à une erreur de calcul (singularité).

La condition (3) est moins évidente, mais le graphique suivant illustre un cas de non-convergence. Dans ce cas, le processus à une boucle calculant alternativement  et .

  (170)

Si la fonction f est donnée analytiquement, sa dérivée peut-être déterminée analytiquement. Mais dans bien des cas, il est utile, voire indispensable de remplacer  par le quotient différentiel:

  (171)

où h doit être choisi suffisamment petit pour que la différence:

  (172)

soit elle aussi suffisamment petite.

L'itération s'écrit alors:

  (173)

Convergence de la méthode:

Si la méthode de résolution est convergent, l'écart entre  et  diminue à chaque itération. Ceci est assuré, par exemple, si l'intervalle  contenant , voit sa longueur diminuer à chaque étape. La méthode de Newton est intéressante car la convergence est quadratique:

  (174)

alors que la convergence des autres méthodes est linéaire:

  (175)

Considérons, par exemple, la méthode de la bissection vue précédemment. A chaque itération la longueur de l'intervalle  diminue de moitié. Ceci nous assure que l'écart  est réduit de moitié à chaque étape du calcul:

  (176)

Pour démontrer la convergence quadratique de la méthode de Newton, il faut utiliser les développements limités de f et f ' au voisinage de :

  (177)

Mais:

  (178)

donc:

  (179)

En soustrayant  à gauche et à droite de l'égalité et en mettant les deux termes du seconde membre au même dénominateur, il vient:

  (180)

et dès que  est assez petit, le dénominateur peut être simplifié.

  (181)

ce qui montre bien que la convergence est quadratique.

AIRES ET SOMMES DE RIEMANN

Considérons la figure suivante :

  (182)

Nous désirons calculer l'aire comprise entre l'axe x, la courbe de f et les droites d'équations  et . Nous supposons dans ce cas que la fonction f est à valeurs positives:

  (183)

Ce problème, dans sa généralité, est difficile voire impossible à résoudre analytiquement. Voici donc quelques méthodes numériques permettant le calcul approché de cette aire.

MÉTHODE DES RECTANGLES

Nous subdivisons l'intervalle  en n sous-intervalles dont les bornes sont . Les longueurs de ces sous intervalles sont . Nous construisons les rectangles dont les côtés sont  et .

  (184)

L'aire de ces rectangles vaut:

  (185)

Si les  sont suffisamment petits,  est une bonne approximation de l'aire cherchée. Nous pouvons recommencer cet exercice en choisissant  et  comme côtés des rectangles. Nous obtenons alors:

  (186)

La figure correspondante est la suivante:

  (187)

Encore une fois, l'aire de ces rectangles approche l'aire cherchée. Afin de simplifier la programmation, il est utile de choisir des intervalles de longueur identique:

  (188)

Si nous avons n rectangles, h vaut alors . Les aires  et  deviennent:

  (189)

MÉTHODE DES TRAPÈZES

Afin d'augment la précision des calculs, il est possible de calculer:

  (190)

Dans le cas où tous les intervalles sont de longueur égale,  vaut:

  (191)

Il existe une foule d'autres méthodes permettant la résolution de ce problème (dont la méthode de Monte-Carlo que nous verrons plus loin).

Dans le cas où la fonction f n'est pas à valeurs positives, nous ne parlons plus d'aire mais de "somme de Riemann". Les sommes à calculer sont alors:

  (192)

et:

  (193)

Tous les calculs doivent êtres conduits de la même manière, mais les résultats peuvent être positifs, négatifs ou nuls.

PROGRAMMATION LINÉAIRE

L'objectif de la programmation linéaire est de trouver la valeur optimale d'une fonction linéaire sous un système d'équations d'inégalités de contraintes linéaires. La fonction à optimiser est baptisée "fonction économique" (utilisée en économie dans le cadre d'optimisations) et on la résout en utilisant une méthode dite "méthode simplexe" (voir plus loin) dont la représentation graphique consiste en un "polygone des contraintes". 

Remarques:

R1. La programmation linéaire est beaucoup utilisée (pour ne citer qu les cas les plus connus) dans la logistique, la finance d'entreprise ou encore aussi en théorie de la décision lorsque nous devons résoudre un jeu à stratégie mixte (voir le chapitre de théorie de la décision et des jeux pour voir un exemple pratique).

R2. Dans le cadre de résolution de problèmes où interviennent des produits de deux variables nous parlons alors logiquement "programmation quadratique". C'est typiquement le cas en économétrie dans la modélisation des portefeuilles (cf. chapitre d'Econométrie).

R3. La programmation quadratique et linéaire sont réunies dans l'étude générale de ce que nous appellons la "recherche opérationnelle".

La recherche opérationnelle à pour domaine l'étude de l'optimisation de processus quels qu'ils soient. Il existe de nombreux algorithmes s'inspirant des problèmes du type exposés lors de notre étude de la programmation linéaire. Nous nous attarderons en particulier sur l'algorithme plus utilisé qui est "l'algorithme du simplexe".

Lorsqu'on peut modéliser un problème sous forme d'une fonction économique à maximiser dans le respect de certaines contraintes, alors on est typiquement dans le cadre de la programmation linéaire.

Soit une fonction économique Z telle que:

  (194)

où les  sont des variables qui influent sur la valeur de Z, et les  les poids respectifs de ces variables modélisant l'importance relative de chacune de ces variables sur la valeur de la fonction économique.

Les contraintes relatives aux variables s'expriment par le système linéaire suivant:

  (195)

Sous forme générale et matricielle ce genre de problème s'écrit:

  (196)

Voyons un exemple qui consiste à résoudre le problème simple suivant:

Une usine fabrique 2 pièces P1 et P2 usinées dans deux ateliers A1 et A2. Les temps d'usinage sont pour P1 de 3 heures dans l'atelier A1 et de 6 heures dans l'atelier A2 et pour P2 de 4 heures dans l'atelier A1 et de 3 heures dans l'atelier A2.

Le temps de disponibilité hebdomadaire de l'atelier A1 est de 160 heures et celui de l'atelier A2 de 180 heures.

La marge bénéficiaire est de 1'200.- pour une pièce P1 et 1'000.- pour une pièce P2.

Quelle production de chaque type doit-on fabriquer pour maximiser la marge hebdomadaire?

Le problème peut se formaliser de la façon suivante:

  (197)

La fonction économique étant:

  (198)

à maximiser.

Résolution graphique du problème (ou méthode du "polygone des contraintes"): les contraintes économiques et de signe sont représentées graphiquement par des demi-plans. Les solutions, si elles existent appartiennent donc à cet ensemble appelé "région des solutions admissibles" :

  (199)

Remarque: Dans le cas général, pour ceux qui aiment le vocabulaire des mathématiciens..., la donnée d'une contrainte linéaire correspond géométriquement à la donnée d'un demi-espace d'un espace à n dimensions (n étant le nombre de variables). Dans les cas élémentaires, l'ensemble des points de l'espace qui vérifient toutes les contraintes est un convexe limité par des portions d'hyperplan (voir le cas 2 variables, facile à illustrer). Si la fonction de coût est linéaire, l'extrèmum est à un sommet (facile à voir). L'algorithme du simplex de base (voir plus loin) part d'un sommet et va au sommet d'à côté qui maximise localement le coût. Et recommence tant que c'est possible.

Pour trouver les coordonnées des sommets, on peut utiliser le graphique si les points sont faciles à déterminer.

Il s'agit donc de chercher à l'intérieur de ce domaine (connexe), le couple maximisant la fonction économique.

Or, l'équation Z est représentée par une droite de pente constante (-1.2) dont tous les points fournissent la même valeur Z pour la fonction économique.

En particulier, la droite passe par l'origine et donne une valeur nulle à la fonction économique. Pour augmenter la valeur de Z et donc la fonction économique, il suffit d'éloigner de l'origine (dans le quart de plan) la droite de pente -1.2.

Pour respecter les contraintes, cette droite sera déplacée, jusqu'à l'extrême limite où il n'y aura plus qu'un point d'intersection (éventuellement un segment) avec la région des solutions admissibles.

  (200)

La solution optimale se trouve donc nécessairement sur le pourtour de la région des solutions admissibles et les parallèles formées par la translation de la fonction économique s'appellent les "droites isoquantes" ou "droites isocoûts"...

Voyons maintenant comment résoudre ce problème de manière analytique avant de passer à la partie théorique.

Nous avons donc le "système canonique" :

  (201)

avec :

  (202)

Nous introduisons d'abord des "variables d'écart" afin de transformer les 2 inégalités par des égalités. Le système d'équations devient alors une "forme standard" :

 

  (203)

Remarque: Il y a autant de variables d'écart que d'inéquations!

La situation peut se résumer dans le tableau suivant (nous omettons la représentation des variables d'écart dans le tableau-matrice qui ne servent qu'à égaliser les équations) :

			Contraintes
			Total
	3	4	160
	6	3	180
Fonction économique	1'200	1'000

Nous déterminons maintenant le pivot (voir plus loin la méthode du pivot), pour cela nous choisissons la colonne où le coefficient économique est le plus grand. Ici c'est la colonne 1.

Ensuite, nous effectuons les procédures suivantes :

1. Le pivot est remplacé par son inverse

2. On divise les éléments de la ligne du pivot (pivot exclu) par le pivot

3. On divise les éléments de la colonne du pivot (pivot exclu) par le pivot mais on change leur signe ensuite

4. Pour les autres éléments de la première ligne : élément de la ligne 1 diminué de l'élément correspondant sur la ligne de pivot multiplié par 3/6 (rapport des valeurs dans la colonne de pivot)

Nous obtenons dès lors :

			Contraintes
			Total
Fonction économique

Ce qui donne :

			Contraintes
			Total
	0.5	2.5	70
	0.166	0.5	30
Fonction économique	-200	400

Nous n'atteignons la solution optimale que lorsque tous les éléments de la marge sont négatifs ou nuls. Il faut donc continuer (car il reste 500 dans la colonne ) ... ici, on atteint déjà l'optimum au troisième tableau, mais ce n'est pas une généralité (le pivot est 2.5 cette fois). On recommence dans les opérations :

			Contraintes
			Total
Fonction économique

Ce qui donne :

			Contraintes
			Total
	-0.2	0.4	28
	0.266	-0.2	16
Fonction économique	-120	-160

Le processus est terminé car tous les termes de la fonction économique sont négatifs. Le programme optimum est donc de et pour un résultat de :

  (209)

ALGORITHME DU SIMPLEXE

Pour mettre en œuvre cet algorithme, nous dveons poser le problème sous un forme "standard" et introduire la notion de "programme de base" qui est l'expression algébrique correspondant à la notion de "point extrême du polyèdre des programmes admissibles" étudiée lors de la programmation linéaire (noté ci-après P.L.). En effet, nous verrons que la solution d'un problème de type P.L. si elle existe, peut toujours être obtenue en un programme de base. La méthode du simplexe va donc consister à trouver un premier programme de base puis à construire une suite de programmes de base améliorant constamment la fonction économique et donc conduisant à l'optimum.

Un problème de P.L. est donc mis sous sa "forme standard" s'il implique la recherche du minimum de la fonction objectif sous des contraintes ayant la forme d'équation linéaires et de conditions de non négativité des variables, c'est-à-dire s'il se pose sous la forme que nous avons vu lors de notre étude de la programmation linéaire:

  (210)

C'est-à-dire aussi, en utilisant des notation matricielles:

  (211)

où les matrices  correspondent, respectivement, aux coefficients des niveaux d'activité dans la fonction objectif, aux coefficients techniques des activités et aux seconde membres des contraintes.

Nous allons voir maintenant comme un problème général de P.L. peut toujours être ramené à une forme standard. La notion de "variable d'écart" est essentielle pour effectuer cette "réduction".

Chercher le maximum d'une fonction f(x) revient à chercher le minimum de la fonction de signe opposé -f(x) . D'autre part une contrainte qui se présente comme une inéquation:

  (212)

peut être remplacée par l'équation:

  (213)

impliquant une variable supplémentaire, , appelée donc "variable d'écart", et soumise à la contrainte de non-négativité, .

Bien évidemment, dans un cas contraire tel où le système est du type:

  (214)

Nous écrirons:

  (215)

impliquant donc également une variable supplémentaire et soumise à la contrainte de non-négativité, .

Ce travail de mise en forme standard nous permet donc de nous retrouver un système d'équation linéaires à résoudre (nous avons vu précédemment sur le site comme résoudre ce genre de système avec l'algorithme du pivot).

La matrice A qui représente les composantes du système d'équations peut s'exprimer dans différentes variantes en fonction de la base vectorielle choisie (voir le chapitre d'analyse vectorielle dans la section d'algèbre). Nous allons introduire la notion de "forme canonique utilisable" associée au choix d'une base et nous montrerons que cette reformulation du système de contraintes va nous permettre de progresser vers l'optimum.

La matrice A peut, après introduction des variables d'écart se décompenser en deux sous-matrices , une contenant les variables initiales D et l'autre comportant les variables d'écart B tel que:

  (216)

Remarque: Les variables d'écart sont des variables et non des constantes!! Il convient dans un système où les variables sont au nombre de n tel qu'une équation du système s'écrirait:

  (217)

d'ajouter une variable d'écart tel que:

   (218)

où .

et sur chaque ligne m, la variable d'écart ajoutée doit être différentes de celles déjà insérées dans le système. C'est la raisons pour laquelle nous pouvons décomposer la matrice en deux-sous matrices.

Les colonnes de la matrice B sont bien évidemment, par définition de la méthode, des colonnes unités, linéairement indépendantes. Ces colonnes forment une base de l'espace vectoriel des colonnes à m éléments (ou dimensions) – le nombre de lignes du système. Nous appelons B la "matrice de la base".

Les variables associées aux composantes colonnes de la matrice B seront dès maintenant appelées les "variables de bases". Dans ce cas, les variables de base sont donc essentiellement les variables d'écart . Les variables associées aux colonnes de la matrice D seront appelées les "variables hors-base"; il s'agit des variables .

Remarque: Rappelons que dans l'expression de la fonction économique, seule les variables hors-base apparaissent.

En résumé, tout P.L. une fois mis sous forme standard est tel que:

- il existe une sous-matrice carrée de la matrice A des coefficients techniques, qui est appelée matrice de base et qui est égale à la matrice carrée unité I de dimension  (effectivement il y autant de variables d'écart que de lignes dans le système d'équations original – au nombre de m - et autant de colonnes puisque chaque variable d'écart à un indice différent).

- les variables de base associées n'apparaissent pas dans l'expression de la fonction économiques

- le second membre des contraintes est constitué d'élément non négatifs

Nous disons que le problème est mis sous "forme canonique utilisable associée à la base B, correspondant aux variables de base ".

Remarque: Nous pouvons intervertir les matrices (et donc changer de base canonique) B et D (ce qui revient à dire que la matrice des variables de base devient la matrices des variables hors-base et inversement).

Il est maintenant commode d'introduire les notations suivantes:

  (219)

qui sont respectivement le vecteurs des variables de base et le vecteur des variables hors-base.

Ainsi, le système d'équations décrivant les contraintes peut s'écrire indifféremment:

  (220)

ou bien aussi:

  (221)

Si la matrice B est une matrice de base, elle est non singulière et admet donc une matrice inverse . En multipliant cette équation, à gauche et à droite, par  nous obtenons:

  (222)

Le système d'équations aura alors été mis sous une forme résolue en .

Pour obtenir une forme canonique utilisable associée à la base B, correspondant aux variables de base, il ne reste plus qu'à éliminer les variables de base de l'expression de la fonction économique.

Ecrivons cette fonction en séparant les vecteurs  et , nous obtenons:

  (223)

Nous pouvons alors facilement exprimer  en fonction de . En utilisant le système d'équations mis sous forme résolue en , nous avons dans un premier temps:

  (224)

que nous substituons dans la fonction économique, pour obtenir:

  (225)

Nous regroupons les termes en  et nous avons:

  (226)

Nous avons alors toujours système d'équations mais ne comportant plus d'inégalités mais au contraire des égalités !!! Reste plus qu'à démontrer que la solution de ce système dit "programme base" par la méthode du pivot est optimal.

Définition: nous appelons coût réduit de l'activité hors base j, le coefficient correspondant  de la ligne .

Soit un problème de programmation linéaire sous forme standard:

  (227)

La matrice A à m lignes (autant qu'il y a de contraintes) et n colonnes, avec . Si nous sélectionnons m variables de base et si nous annulons les  variables hors base, la matrice A:

  (228)

et le système se réduit à:

  (229)

La matrice B est de dimension . Si elle définit une base, elle admet une matrice inverse . Une solution du système est donc:

  (230)

Si l'expression  est non négative , nous avons une solution admissible qui vérifie les contraintes et que nous appellerons un programme de base:

  (231)

Le problème de programmation linéaire, s'écrit aussi sous la forme suivante, que nous appelons "forme canonique utilisable associée au programme de base":

  (232)

Définition: Nous appelons "coût réduit" de l'activité hors base j, le coefficient correspondant  de chaque ligne j de l'expression

A partir des développements effectués précédemment nous pouvons énoncer le résultat suivant:

Proposition 1: si dans la forme canonique utilisable associée à un programme de base, tous les coûts réduits sont  alors le programme de base est optimal.

Proposition 2: La solution d'un problème de P.L., si elle existe, peut toujours être trouvée en un programme de base.

La démonstration précise de ce résultat est assez délicate. Nous pouvons cependant en avoir une intuition en considérant, une fois de plus, la notion de coût réduit.

En effet, pour un programme de base donné, considérons la forme canonique utilisable associée à la base. Sur cette forme nous pouvons vérifier que, ou bien le programme de base est optimal (tous les coûts réduits sont ), ou bien que le programme de base peut être amélioré et remplacé par un nouveau programme de base donnant à z une valeur plus petite (un coût réduit est négatif et la variable hors-base associée peut être augmentée jusqu'à ce qu'une ancienne variable de base s'annule). Comme il y a un nombre fini de programmes de base (au plus égal au nombre ), la solution de P.L. se trouve nécessairement en un programme de base.

MÉTHODE DE MONTE-CARLO

La méthode de Monte-Carlo est un moyen très efficace de contourner les problèmes mathématiques et physiques les plus complexes. Elle trouve ses applications dans des domaines variés dont voici quelques exemples:

- Aiguille de Buffon

- Problèmes de neutronique liés à la bombe atomique

- Calcul d'intégrale ou d'espérance de variables aléatoires

- Résolution d'équations elliptiques ou paraboliques

- Résolution de systèmes linéaires

- Résolution de problèmes d'optimisation (recherche opérationnelle, gestion de projets)

- Finance

Il existe 2 types de problèmes qui peuvent être traités par la méthode de Monte-Carlo: les problèmes probabilistes, qui ont un comportement aléatoire, et les problèmes déterministes, qui n'en ont pas.

Pour ce qui est du cas probabiliste, il consiste à observer le comportement d'une série de nombres aléatoires qui simule le fonctionnement du problème réel et à en tirer les solutions.

Pour le cas déterministe, le système étudié est complètement défini et on peut en principe prévoir son évolution, mais certains paramètres du problème peuvent être traités comme s'il s'agissait de variables aléatoires. Le problème déterministe devient alors probabiliste et ré solvable de façon numérique. On parle alors d'estimation de "Monte-Carlo" ou d'une approche de "MC élaborée".

La dénomination de méthode de "Monte-Carlo" date des alentours de 1944. Des chercheurs isolés ont cependant utilisé bien avant des méthodes statistiques du même genre: par exemple, Hall pour la détermination expérimentale de la vitesse de la lumière (1873), ou Kelvin dans une discussion de l'équation de Boltzmann (1901), mais la véritable utilisation des méthodes de Monte-Carlo commença avec les recherches sur la bombe atomique.

Au cours de l'immédiat après-guerre, Von Neumann, Fermi et Ulam avertirent le public scientifique des possibilités d'application de la méthode de Monte-Carlo (par exemple, pour l'approximation des valeurs propres de l'équation de Schrödinger). L'étude systématique en fut faite par Harris et Hermann Khan en 1948. Après une éclipse due à une utilisation trop intensive pendant les années 1950, la méthode de Monte-Carlo est revenue en faveur pour de nombreux problèmes: en sciences physiques, en sciences économiques, pour des prévisions électorales, etc., bref, partout où il est fructueux d'employer des procédés de simulation.

Le mieux pour comprendre la méthode de Monte-Carlo c'est d'abord d'avoir un très bon générateur de nombres aléatoire (ce qui est très difficile). Prenons comme exemple le générateur de Maple:

rand();

restart;rand();

 

Nous voyons donc que la fonction par défaut de générateur de nombres aléatoires de Maple est à utiliser avec la plus grand prudence puisque un ré-initialisation du système suffit à retrouver des valeurs aléatoires… égales. Cependant il existe des libraires spécialisées dans Maple tel que:

restart;readlib(randomize):randomize():rand();

> restart;readlib(randomize):randomize():rand();

Epreuve à priori réussie (au fait, il nous faudrait faire un beaucoup plus grand nombre d'essais afin de bien vérifier que le générateur ne suive pas une loi de distribution connue… ce qui n'est malheureusement jamais le cas).

Une fois le générateur crée et testé, nous pouvons voir quelques résultat de la méthode de Monte-Carlo. Ainsi, dans le calcul des intégrales, celle-ci s'avère très utile et très rapide en terme de vitesse de convergence.

CALCUL D'UNE INTÉGRALE

Soit à calculer l'intégrale d'une fonction f définie et positive sur l'intervalle :

  (233)

Soit :

  (234)

Nous considérons le rectangle englobant de la fonction sur défini par les points  . Nous tirons un grand nombre N de points aux hasard dans ce rectangle. Pour chaque point, nous testons s'il est au-dessous de la courbe. Soit F la proportion de points situés au-dessous, nous avons:

  (235)

L'algorithme Maple est donné par:

intmonte:=proc(f,a,b,N)
local i,al,bl,m,F,aleaabs,aleaord,estaudessous;
m:=round(max(a,b)*10^4);
al:=round(a*10^4);
bl:=round(b*10^4);
aleaabs:=rand(al..bl);
aleaord:=rand(0..m);
F:=0;
for i from 1 to N do
     estaudessous:=(f(aleaabs()/10^4)-aleaord()/10^4)>=0;
     if estaudessous then
          F:=F+1;
     fi
od:
RETURN((b-a)*max(a,b)*F/N)
end:

Remarque: Pour appeler cette procédure il suffit d'écrire >intmonte(f,a,b,N) mais en remplaçant le premier argument passé en paramètre par l'expression d'une fonction et les autres arguments par des valeurs numériques bien évidemment.

CALCUL DE PI

Pour le calcul de  le principe est le même et constitue donc à utiliser la proportion du nombre de points dans un quartier de cercle (cela permet de simplifier l'algorithme en se restreignant aux coordonnées strictement positives) inscrit dans un carré relativement au nombre de points total (pour tester si un point est à l'extérieur du cercle, nous utilisons bien évidemment le théorème de Pythagore) tel que:

  (236)

L'algorithme Maple est donné par:

estalinterieur:=proc(x,y) x^2+y^2<1 end:
calculepi:=proc(N)
local i,F,abs,ord,alea,erreur,result;
alea:=rand(-10^4..10^4);
F:=0;
for i from 1 to N do
     abs:=alea()/10^4;ord:=alea()/10^4;
       if estalinterieur(abs,ord) then
            F:=F+1;
       fi
od;
RETURN(4*F/N)
end:

DICHOTOMIE

La dichotomie consiste pour un objet de taille N à exécuter un algorithme de façon à réduire la recherche à un objet de taille . On répète l'algorithme de réduction sur ce dernier objet. Ce type d'algorithme est souvent implémenté de manière récursive. Lorsque cette technique est utilisable, elle conduit à un algorithme très efficace et très lisible.

Un exemple simple est la recherche de la racine d'une fonction continue (nous avons déjà étudié différentes méthodes plus haut pour résoudre ce genre de problèmes).

Supposons qu'une fonction soit croissante et continue sur un intervalle . Nous avons donc . Nous  calculons . Si  alors la racine est dans l'intervalle  sinon elle est dans l'intervalle . Nous a donc ramené le problème à une taille inférieure. Nous arrêterons l'algorithme quand la précision sera suffisante.

L'algorithme Maple est donné par:

zero:=proc(f,a,b,pre)
local M;
M:=f((a+b)/2);
if abs(M)<pre then 
     RETURN((a+b)/2)
elif M>0 then
     zero(f,a,(a+b)/2,pre)
else zero(f,(a+b)/2,b,pre)
     fi
end:

et ce ne sont que quelques exemples auxquels la méthode est applicable!!