TABLE DES MATIERES

INTRODUCTION 1
1 METHODOLOGIE DE TRAITEMENT NUMERIQUE DES PROBLEMES SCIEN-
TIFIQUES, CONCEPT DE BASE 2
1.1 NOTIONS DE BASE EN CALCUL NUMERIQUE . . . . . . . . . . . . . . . . . . 2
1.1.1 Utilisation des réels . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Utilisation des fonctions . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.3 Discrétisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.4 Itérations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.5 Erreurs d’arrondis et de troncature . . . . . . . . . . . . . . . . . . . 3
1.2 PROBLEMES ET METHODES INSTABLES . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 Problèmes instables ou mal conditionnés . . . . . . . . . . . . . . . . 3
1.2.2 Méthodes instables . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 METHODOLOGIE DE TRAITEMENT NUMERIQUE DES PROBLEMES SCIENTI-
FIQUES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 Le problème posé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.2 La méthode de résolution . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.3 L’algorithme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.4 La programmation informatique . . . . . . . . . . . . . . . . . . . . . 4
1.3.5 Le traitement machine . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3.6 Interprétation des résultats . . . . . . . . . . . . . . . . . . . . . . . 5
2 INTERPOLATION ET APPROXIMATION DES FONCTIONS 6
2.1 GENERALITES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.2 Types d’interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 INTERPOLATION POLYNOMIALE . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 Formule d’interpolation de Lagrange . . . . . . . . . . . . . . . . . . 7
2.2.2 Formule d’interpolation de Newton . . . . . . . . . . . . . . . . . . . 10
2.2.3 Erreurs d’interpolation polynomiale . . . . . . . . . . . . . . . . . . . 10
2.2.4 Interpolation spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.5 Autres types d’interpolation polynomiale . . . . . . . . . . . . . . . . 12
2.3 INTERPOLATION TRIGONOMETRIQUE, DOUBLE ET INVERSE . . . . . . . . . 12
2.3.1 Interpolation trigonométrique . . . . . . . . . . . . . . . . . . . . . . 12
2.3.2 Interpolation double . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.3 Interpolation inverse . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 APPROXIMATION PAR LA METHODE DES MOINDRES CARRES . . . . . . . . 13
2.4.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4.2 Résolution du problème de l’approximation par les moindres carrés . 14

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

TABLE DES MATIERES II

3 INTEGRATION ET DIFFERENTIATION 16
3.1 DERIVATION NUMERIQUE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.1.1 Formules classiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.1.2 Formules plus précises . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 INTEGRATIONS NUMERIQUES- FORMULES DE NEWTON-COTES . . . . . . . . 17
3.2.1 Méthode des trapèzes : . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2.2 Algorithme : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.3 Méthode de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 INTEGRALES MIXTE ET DOUBLE . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3.1 Intégrale mixte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3.2 Interpolation double . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4 RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUA-
TIONS INTEGRALES 21
4.1 GENERALITES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1.1 Définition et classification . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1.2 Equations différentielles . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1.3 Position du problème Numérique . . . . . . . . . . . . . . . . . . . . 21
4.2 METHODES POUR LES PROBLEMES A VALEURS INITIALES . . . . . . . . . . . 22
4.2.1 Méthode de Picard . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.2.2 Séries de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2.3 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2.4 Les méthodes de Runge et Kutta . . . . . . . . . . . . . . . . . . . . 24
4.2.5 Méthodes à pas multiples ou d’Adams Bashforth . . . . . . . . . . . . 26
4.2.6 Méthode de Prédicteur-Correcteur . . . . . . . . . . . . . . . . . . . . 27
4.3 METHODE POUR LES PROBLEMES AUX VALEURS AUX LIMITES . . . . . . . . 28
4.3.1 Méthode de tir . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.3.2 Méthode des fonctions complémentaires ou de superposition . . . . . 29
4.3.3 Méthode des différences finies . . . . . . . . . . . . . . . . . . . . . . 30
4.4 LES EQUATIONS INTEGRALES . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.4.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.4.2 Résolution numérique des équations intégrales linéaires . . . . . . . . 31
5 SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES
33
5.1 METHODES DIRECTES POUR LES SYSTEMES D’EQUATIONS LINEAIRES . . . . 33
5.1.1 Rappel de la méthode de Gauss . . . . . . . . . . . . . . . . . . . . . 33
5.1.2 Méthode de Jordan . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.1.3 Méthode de la décomposition triangulaire . . . . . . . . . . . . . . . 34
5.1.4 Méthode de la racine carré de Cholesky . . . . . . . . . . . . . . . . . 36
5.1.5 Fragmentation des systèmes volumineux . . . . . . . . . . . . . . . . 36
5.1.6 Systèmes à coefficients complexes . . . . . . . . . . . . . . . . . . . . 37
5.1.7 METHODES ITERATIVES . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.2 RAPPELS SUR LES METHODES DE NEWTON ET DE BISSECTION POUR LES
EQUATIONS ALGEBRIQUES NON LINEAIRES . . . . . . . . . . . . . . . . . . . 38
5.2.1 Méthode de Newton-Raphson(1960) . . . . . . . . . . . . . . . . . . . 38
5.2.2 Méthode de Dichotomie ou Méthode de Bissection . . . . . . . . . . . 41
5.3 CALCUL DE TOUTES LES RACINES D’UNE EQUATION POLYNOMIALE NON LI-
NEAIRE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.3.1 Opération sur les polynômes . . . . . . . . . . . . . . . . . . . . . . . 41
5.3.2 Calcul d’une racine réelle, passage à l’équation de degré n − 1 . . . . 43
5.3.3 Calcul de deux racines et passage à l’équation de degré n − 2 . . . . . 44
5.3.4 Formule d’itération de Newton sur l’axe réel . . . . . . . . . . . . . . 44
5.3.5 Méthode de Bairstow . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.3.6 Calcul simultané de toutes les racines . . . . . . . . . . . . . . . . . . 46

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

TABLE DES MATIERES III

5.4 SYSTEMES D’EQUATIONS NON LINEAIRES . . . . . . . . . . . . . . . . . . . . 48

6 EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABI-
LITE 50
6.1 NOTIONS GENERALES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.1.1 Définitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.1.2 Quelques équations aux dérivées partielles en Physique et en Mécanique 50
6.2 DISCRETISATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.2.1 L’équation de Poisson dans le plan . . . . . . . . . . . . . . . . . . . 51
6.2.2 Discrétisation de l’équation de la chaleur . . . . . . . . . . . . . . . . 53
6.2.3 Discrétisation de l’équation des ondes . . . . . . . . . . . . . . . . . . 55
7 RESOLUTION NUMERIQUE DES EQUATIONS DE NAVIER STOKES ET
DE BURGERS 59
7.1 EQUATION DE NAVIER STOKES . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.1.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.1.2 Discrétisation de l’équation d’advection-diffusion . . . . . . . . . . . . 61
7.2 DISCRETISATION DE L’EQUATION DE NAVIER-STOCKES DANS LE PLAN . . . . 62
7.3 EQUATIONS DE BURGERS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
7.3.1 Présentation de l’équation . . . . . . . . . . . . . . . . . . . . . . . . 62
7.3.2 Schémas de discrétisation de l’équation de Burgers sans terme de viscosité 63
7.3.3 Schéma de discrétisation de l’équation de Burgers visqueux . . . . . . 64
EXERCICES 65

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 INTRODUCTION

Le but de ce cours est de permettre aux scientifiques en général, aux physiciens en particu-
luer d’élaborer des méthodes de calcul numériques utilisables par l’ordinateur pour résoudre
les problèmes de recherche et de développement parmi lesquels :
– Le traitement d’un grand nombre de données issues des mesures expérimentales,
– La résolution de grands systèmes d’équations à plusieurs inconnues,
– La résolution des systèmes d’équations algébriques linéaires, des équations différentielles
et intégrales ainsi que des équations aux dérivées partielles.
Bien qu’il existe sur le marché des logiciels conçus pour traiter la quasi-totalité de des
problèmes su-mentionés, il se trouve que la plupart des problèmes de recherche sont com-
plexes et exigent du chercheur ou de l’ingénieur qu’il écrive lui-même son programme en
tenant compte de tous les critères et conditions particulièrs. De plus, faire fonctionner ces
logiciels demande au chercheur ou à l’ingénieur de comprendre correctement les différentes
étapes à suivre afin d’introduire les données dans le logiciel et recueillir les résultats dudit
calcul.
Le cours comprend essentiellement six chapitres :
– Interpréation et approximation des fonctions
– Intégration et dérivation numérique,
– Méthodes de résolution numérique des équations différentielles,
– Méthodes de résolution numérique des systémes d’équations algébriques linéaires et non
linéaires,
– Méthodes de résolution numérique des équations aux dérivées partielles,
– Introduction à la méthode des élements finis.
Chaque partie du cours furnira d’une part les fondements théoriques des méthodes nu-
mériques et les algorithmes qui seront traduits dans les langages de programmation tels que
FORTRAN, PASCAL et C. Des exemples de problèmes seront présentés à titre d’illustration.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 Chapitre I

METHODOLOGIE DE TRAITEMENT NUMERIQUE

DES PROBLEMES SCIENTIFIQUES, CONCEPT DE
BASE

1.1 NOTIONS DE BASE EN CALCUL NUMERIQUE

1.1.1 Utilisation des réels

Le terme calcul numérique se rapporte à l’utilisation des calculateurs scientifiques pour ré-
soudre les problèmes faisant intervenir les nombres entiers et réels. La plupart de ces nombres
parmi lesquels tous les nombres entiers, s’expriment par une suite finie de chiffres.√ D’autres
par contre nécessitent pour leur représentation une suite infinie de chiffres (π, e, 2, 103
, . . .).
Les ordinateurs quant à eux ne peuvent représenté les réels et les entiers que par une suite
finie, ce qui entraîne des erreurs d’arrondis. Chaque nombre de cette suite est appelé digit.
La longueur de la suite finie est un élément qui entre dans les qualités de performance et
de précision de l’ordinateur. On peut améliorer la précision en travaillant en double ou en
quadruple précision (représentation en binaire ou en octale).

1.1.2 Utilisation des fonctions

Il y a deux manières de représenter les fonctions dans l’ordinateur : par un tableau de
valeurs ou par un sous-programme qui calcule les valeurs de la fonction à des points choisis.

1.1.3 Discrétisation
En calcul numérique, certains problèmes sont par nature continus. Ils font intervenir des
opérations telles que la dérivation et l’intégration. L’ordinateur ne pouvant résoudre de tels
problèmes, on doit les remplacer par des problèmes discrets en utilisant par exemple les
différence finies, les éléments finis, les méthodes particulaires ou spectrales : on parle alors
de discrétisation.

1.1.4 Itérations
Certains problèmes numériques sont souvent formulés en terme de processus successifs ou
d’une suite de calcul, le résultat d’un processus étant lié à celui du processus précédent. On
parle alors d’itération.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

METHODOLOGIE DE TRAITEMENT NUMERIQUE DES PROBLEMES SCIENTIFIQUES, CONCEPT DE BASE 3

1.1.5 Erreurs d’arrondis et de troncature

Les erreurs de calcul numérique associées aux ordinateurs peuvent être classées en deux
catégories : les erreurs dues à la représentation limitée des réels et celles dues au manque
de ressources de l’ordinateur. La première catégorie est celle des erreurs d’arrondis qui aug-
mentent avec les opérations arithmétiques. La deuxième catégorie est celle des erreurs de
troncature. Elles surviennent par exemple lors de l’approximation des fonctions par des fonc-
tions rationnelles. Généralement, l’on remarque que la réduction des erreurs de troncature
entraîne l’augmentation des erreurs d’arrondis.

1.2 PROBLEMES ET METHODES INSTABLES

1.2.1 Problèmes instables ou mal conditionnés

Une difficulté courante et source des erreurs graves en calcul numérique réside dans l’ex-
trême sensibilité de la solution aux légères variations des paramètres du problème. On pour-
rait parler du chaos numérique. Cette difficulté peut survenir dans n’importe quel type de
calcul numérique : soustractions, systèmes d’équations, équations non linéaires, équations
différentielles. Par exemple, le système suivant est instable :
½
x+y = 3
0.499x + 1.001y = 1.5
Sa solution est (1; 1). Si dans la deuxième équation, nous remplaçons 0.499 par 0.5 on obtient
comme solution (3; 0), ce qui est très éloigné de la première. Ce genre de problème est dit in-
stable ou mal conditionné. Si des variations légères des paramètres d’un problème entraînent
des changements légers de la solution, on dit que le problème est stable ou bien conditionné
(bien posé). Il faut donc noter que lorsque le problème est instable, les erreurs d’arrondis
et de troncature ainsi que d’autres approximations peuvent conduire à des solutions très
éloignées de la solution exacte. Notons également qu’un problème peut être instable pour
une méthode numérique et être stable pour une autre méthode. Il faudra donc au cours des
calculs numériques faire des tests de la stabilité du problème ou de la méthode.

1.2.2 Méthodes instables

A cause des opérations arithmétiques et des erreurs d’arrondis subséquentes, les méthodes
de calcul numérique peuvent conduire à des faux résultats. Une méthode qui est bonne,
mais qui souvent, à cause des données d’un problème ou de la manière que ces données
sont transmises, produit des résultats d’une pauvre précision est dite instable. Une méthode
sera d’autant plus stable que les solutions qu’elle donne d’un problème sont proches de la
solution exacte. C’est pourquoi il faut toujours analyser les erreurs d’une méthode de calcul
numérique.

1.3 METHODOLOGIE DE TRAITEMENT NUMERIQUE DES PROBLEMES SCIEN-

TIFIQUES.

Le traitement d’un problème de calcul numérique comprend en général six phases dont
la disposition est la suivante :

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

METHODOLOGIE DE TRAITEMENT NUMERIQUE DES PROBLEMES SCIENTIFIQUES, CONCEPT DE BASE 4

1.3.1 Le problème posé

Le problème doit :
- Être posé clairement
- Avoir un sens
- Etre bien posé afin d’éviter les instabilités

1.3.2 La méthode de résolution

C’est la description du procédé permettant d’obtenir la solution du problème. Elle peut
être faite par une formulation mathématique ou par des directives de type littéral.

1.3.3 L’algorithme
C’est la décomposition en un nombre fini d’opérations élémentaires de la méthode choisie.
L’algorithme peut être complétée par un organigramme. L’organigramme est une étape très
importante du calcul numérique.
- Il confirme ou infirme la réalité du problème posé.
- Il justifie la faisabilité de la méthode proposée.
- Il donne directement accès au programme informatique
- Il permet d’exercer un contrôle sur les performances qui résulteront du traitement infor-
matique : temps d’exécution, encombrement mémoire, précision.
Optimiser un algorithme signifie :
– réduire le temps de calcul machine
– réduire la place occupée en mémoire centrale
– augmenter la précision

1.3.4 La programmation informatique

C’est la traduction de l’algorithme en un langage informatique. Les langages informatiques
les plus utilisés sont le FORTRAN, le PASCAL et le C. ce sont des langages qui évoluent
de jour en jour (FORTRAN1, 2, 3, 4, 77, 90) ; PASCAL, TURBO PASCAL et ses sous
versions. La compatibilité de ces versions est ascendante : Le FORTRAN 66 passe dans un
compilateur de FORTRAN 77.

Remarque : le langage choisi doit pouvoir satisfaire toutes les opérations figurant explici-
tement dans l’algorithme.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

METHODOLOGIE DE TRAITEMENT NUMERIQUE DES PROBLEMES SCIENTIFIQUES, CONCEPT DE BASE 5

1.3.5 Le traitement machine

C’est l’exécution de la méthode de résolution par le programme. Il s’agit ici d’utiliser le
système d’exploitation des ordinateurs : logiciel (compilateur, utilitaire, fonctions et sous
programmes requis), ressources matérielles (taille mémoire, vitesse d’exécution, . . . ).

1.3.6 Interprétation des résultats

Les résultats obtenus doivent être interprétés avec la logique scientifique car il n’existe
aucune règle permettant d’affirmer à priori qu’un programme est tout à fait correct. Il est
donc souvent nécessaire d’essayer plusieurs méthodes de résolution, plusieurs algorithmes ou
même plusieurs langages pour un même problème.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 Chapitre II

INTERPOLATION ET APPROXIMATION DES

FONCTIONS

2.1 GENERALITES

2.1.1 Position du problème

L’interpolation et l’approximation des fonctions sont des techniques numériques très uti-
lisées. Le problème est le suivant : Nous avons une grandeur y qui dépend d’une variable x.
Nous connaissons les valeurs y1 , y2 ,...,yn de y en certains points discrets x1 , x2 ,...,xn de x.
Nous voulons approximer y par une fonction analytique de la variable x, c’est à dire y = f (x)
telle que
y1 = f (x1 ) , y2 = f (x2 ) , . . . , yn = f (xn )
L’expression de f (x) pourra ainsi permettre de déterminer les valeurs de la fonction en
des points autre que ceux de la suite discrète. On parle alors de l’art de lire entre les lignes
d’une table ou d’interpolation.
Il s’agit donc de construire une expression analytique approchée d’une fonction que l’on
ne connaît que pour une suite discrète des valeurs de la variable. Un tel problème n’admet
pas de solution univoquement déterminée.

Remarques :

Ce problème est rencontré regulièrement dans des études expérimentales ou numériques.

2.1.2 Types d’interpolation

La fonction f (x) peut avoir plusieurs formes. Lorsqu’elle est un polynôme, on parle d’in-
terpolation polynomiale. Lorsque f (x) est une fonction trigonométrique, on parle d’interpo-
lation trigono- métrique. On peut aussi avoir une interpolation par les fonctions exponen-
tielles, des polynômes de Legendre, les fonctions de Bessel, etc.
Cependant en pratique, l’on choisit la forme la plus simple pour p (x), notamment la
forme polynomiale. Dans le cas où l’on sait que la fonction est périodique, il est cependant
judicieux d’utiliser l’interpolation trigonométrique. La justification de l’un ou de l’autre type
d’interpolation repose sur 2 théorèmes présentés par Weistrass (en 1885)

Théorème 1 : Toute fonction contenue dans un intervalle [ a , b ] peut être représentée dans
P
N
cet intervalle et pour chaque degré de précision par un polynôme f (x) = an xn .
n=0

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 7

Théorème 2 : Toute fonction continue de période 2π peut être représentée par une série
P
N
trigonométrique limitée de la forme f (x) = (an cos nx + bn sin nx). Les coefficients an et
n=0
bn sont déterminés par la connaissance des valeurs discrètes de la fonction sur l’intervalle
[ a , b].

2.2 INTERPOLATION POLYNOMIALE

Un polynôme de degré n est définie par ses n + 1 coefficients. Ainsi un polynôme d’in-
terpolation de degré n sera complètement déterminé par la connaissance de n + 1 valeurs
discrètes yn . Il existe plusieurs formules pour déterminer l’interpolation polynomiale.

2.2.1 Formule d’interpolation de Lagrange

L’idée de base de cette méthode est de déterminer des polynômes fk (x) de degré n qui
prennent les valeurs 1 au point xk et 0 en tout point différent de xk
½
1 si k 0 = k
⇒ fk (xk0 ) =
0 si k 0 6= k
Le polynôme d’interpolation de Lagrange aura alors la forme :

f (x) = y1 f1 (x) + y2 f2 (x) + · · · + yn fn (x)

P
x+1
⇒ f (x) = yk fk (x)
k=1

Détermination de fk (x) :

Soit par exemple un ensemble de deux points (x1 , y1 ) et (x2 , y2 ). Les polynômes fk sont
de degré 1, c’est à dire : ½
f1 (x) = a1 + b1 x
f2 (x) = a2 + b2 x
½ ½
f1 (x1 ) = 1 f2 (x1 ) = 0
P uisque et
f1 (x2 ) = 0 f2 (x2 ) = 1
Il vient
x − x2 x − x1
f1 (x) = ; f2 (x) = ⇒ y(x) = y1 f1 (x) + y2 f2 (x)
x1 − x2 x2 − x1
Si on avait un ensemble de trois (x1 , y1 ), (x2 , y2 ), et (x3 , y3 ) alors les polynômes seraient
de degré 2 et définis par :

(x − x2 ) (x − x3 ) (x − x1 ) (x − x3 ) (x − x1 ) (x − x2 )
f1 (x) = , f2 (x) = , et f3 (x) = .
(x1 − x2 ) (x1 − x3 ) (x2 − x1 ) (x2 − x3 ) (x3 − x1 ) (x3 − x2 )
Si on a un½ensemble de n + 1 points, les polynômes de degré n tel que
1 si k 0 = k
fk (xk0 ) = sont définis par :
0 si k 0 6= k

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 8

n+1
,n+1 n+1
Y Y X
fk (x) = (x − xi ) (xk − xi ) et f (x) = fk (x) yk
i=1 i=1 k=1
i6=k i6=k

Le polynôme f est unique pour un ensemble de points xi

Exemple : On donne le tableau suivant :

x -1 0 1
donc x1 = −1, x2 = 0, x3 = 1
y 4 1 2
D’après la formule de Lagrange, on a :
(x−x2 )(x−x3 ) x(x−1)
f1 (x) = (x1 −x2 )(x1 −x3 )
= 2
(x−x1 )(x−x3 ) (x+1)(x−1)
f2 (x) = (x2 −x1 )(x2 −x3 )
= 1×(−1)
= − (x2 − 1)
(x−x1 )(x−x2 ) (x+1)x
f3 (x) = (x3 −x1 )(x3 −x2 )
= 2

d’où

f (x) = y1 f1 (x) + y2 f2 (x) + y3 f3 (x)

⇒ f (x) = 2x (x − 1) − (x2 − 1) + x (x + 1)
A partir de l’expression de la fonction de Lagrange, on peut évaluer pour toute valeur xp
comprise entre deux points xk (interpolation) ou en dehors des points xk (extrapolation ) la
valeur de yp correspondante. Pour ce faire, on utilise l’algorithme suivant :

Algorithme

1◦ - Entrer les données (xk , yk , xp ) avec k = 1, 2, . . . , n + 1.

2◦ - yp ← 0
3◦ - pour k allant de 1 jusqu’à n + 1faire
ys ← 1
pour i allant de 1 jusqu’à n + 1 faire
si i 6= k
(xp −xi )
ys ← ys (x k −xi )
fin pour
yp ← yp + yk ys
fin pour
4◦ - écrire xp et yp
Fin

Traduction en PASCAL

Program inter ;
Uses crt ;
Const n =?
var xp , yp , ys : real ;
i, k : integer ;
x, y : array[ 1, . . . , n + 1] of real ;

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 9

begin
x [1] := ; x [2] := ; . . . ; x [n + 1] := ;
y [1] := ; y [2] := ; . . . ; y [n + 1] := ;
yp := 0 ;
for k = 1 to n + 1 do
begin
ys := 1 ;
for i = 1 to n + 1 do
begin
if i <> k then
begin
ys := ys × (xp − x [i])/(x [k] − x [i]) ;
end ;
yp : = yp + y [k] × ys ;
end ;
write(xp , yp ) ;
end.

Traduction en FORTRAN 77

PROGRAM inter ;
PARAMETER (n =?)
REAL xp , yp , ys
INTEGER i, k
REAL DIMENSION x (n + 1) , y (n + 1)
x (1) =
x (2) =
...
x (n + 1) =
y (1) =
y (2) =
...
y (n + 1) =
xp =
yp =
DO 10 k = 1, n + 1
ys = 1
DO 20 i = 1, n + 1
IF (i . N E . k) GO TO 30
30 ys = ys × (xp − x (i))/(x (k) − x (i))
20 CONTINUE
yp = yp + y (k) × ys
10 CONTINUE
WRITE * , xp , yp
STOP
END.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 10

2.2.2 Formule d’interpolation de Newton

a) Par les différences divisées
On appelle différence divisée première des variables xi et xj correspondant aux valeurs yi
et yj la quantité δ (xi , xj ) telle que
yi − yj
δ (xi , xj ) =
xi − xj
La différence divisée seconde est
δ (xi , xj ) − δ (xi , xk )
δ (xi , xj , xk ) = .
xj − xk
La différence divisée d’ordre m est

δ (x1 , . . . , xm−1 ) − δ (x1 , . . . , xm−2 , xm )

δ (x1 , . . . , xm ) = .
xm−1 − xm
Le polynôme d’interpolation de Newton se met sous la forme

y = y1 + δ (x1 , x2 ) (x − x1 ) + δ (x1 , x2 , x3 ) (x − x1 ) (x − x2 ) + . . . ,
+δ (x1 , x2 , . . . , xk+1 ) (x − x1 ) (x − x2 ) . . . (x − xk )
b) Par les différences progressives et régressives
La quantité :
∆yk = yk+1 − yk
est appelée différence première progressive. La différence seconde progressive est
∆2 yk = ∆yk+1 − ∆yk
et la différence progressive d’ordre α est définie par :
∆α yk = ∆α−1 yk+1 − ∆α−1 yk .
On définit également la différence régressive d’ordre α de la manière suivante :

∇α yk = ∇α−1 yk−1 − ∇α−1 yk .

Lorsque l’on a xk+1 = xk + h = x1 + kh, la formule d’interpolation par les différences
divisées se réduit à :
2
f (x) = y1 + ∆y
h
1
(x − x1 ) + ∆h2 y2!1 (x − x1 ) (x − x2 ) + · · ·
∆k+1 y1
+ hk+1 (k+1)!
(x − x1 ) (x − x2 ) . . . (x − xk+1 ) (x − xk )
Cette dernière formule est souvent appelée formule de Newton-Gregory
Si par contre on a xk+1 = xk − h = x1 − kh, alors on peut utiliser la formule de Newton
Grégory régressive en remplaçant dans la formule précédente ∆ par ∇.

2.2.3 Erreurs d’interpolation polynomiale

Dans la méthode de Lagrange, en remplaçant la valeur exacte yexp par yp on commet une
erreur de troncature qui peut être déterminée par le théorème suivant :

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 11

Théorème : Soit f (x) le polynôme de degré n qui satisfait la condition

f (xk ) = yk (k = 1, 2, . . . , n + 1) et définie dans l’intervalle [ a , b ] qui contient les points
xk .
Si f est de classe C n+1 sur [ a , b ], alors pour toute valeur xp ∈ [a , b], il existe une valeur
ξ ∈ [ a , b ] telle que :
£ ± ¤
Yexp − Yp = f (n+1) (ξ) (n + 1)! (xp − x1 ) (xp − x2 ) · · · (xp − xn+1 ) .
Cette différence représente l’erreur de troncature.
En supposant que les dérivées de la fonction f sont uniformément bornées sur [a , b]
l’erreur de troncature décroît lorsque n croît.

2.2.4 Interpolation spline

a) Problème et propriétés :
L’interpolation spline a originellement une formulation en mécanique. Elle est utilisée
(avec la méthode des éléments finis) pour résoudre numériquement les équations aux déri-
vées partielles. Son importance relève des inconvénients de l’interpolation de Lagrange. Le
problème associé à l’interpolation Spline est le suivant :
Soient x1 , x2 , . . . , xn points distincts compris dans un intervalle [a , b] et soient f1 , f2 , · · · , fn
les valeurs de la fonction en ces points, l’interpolation Spline est une interpolation cubique
qui passe par les points (xi , fi ). C’est une fonction S qui admet les propriétés suivantes :
(a) S (xi ) = fi ,
(b) S (x) , S 0 (x) , S 00 (x) sont continues sur [a , b],
(c) Dans chaque intervalle [xi , xi+1 ] pour i = 1, 2, · · · , n + 1, S (x) est un polynôme cubique
(qui peut avoir des coefficients différents dans chaque sous intervalle),
(d) Si g(x) est une autre fonction qui satisfait les conditions (a), (b) et (c) alors on doit
Rb Rb
avoir a (S 00 (x))2 dx ≤ a (g 00 (x))2 dx.
Dans la terminologie de la mécanique, ces propriétés peuvent s’énoncer de la manière
suivante :
(a) La courbe Spline doit passer par tous les nœuds ;
(b) La courbe Spline n’est pas brisée (ou ne courbe pas suivant un angle pointu) ;
(c) Entre deux nœuds, la courbe Spline est un polynôme de degré 3 ;
(d) La courbe Spline minimise l’énergie potentielle du système. Cette propriété implique
aussi que S 00 (a) = S 00 (b) = 0
b) Forme de l’interpolation Spline
Puisque S est un polynôme cubique pour xi ≤ x ≤ xi+1 , S 00 (x) est alors linéaire et la
formule d’interpolation de S 00 (x) est :
x − xi
S 00 (x) = S 00 (xi ) + [S 00 (xi+1 ) − S 00 (xi )] .
xi+1 − x1
Après intégration, on obtient
S 00 (xi+1 ) − S 00 (xi )
S 0 (x) = S 0 (xi ) + S 00 (xi ) (x − xi ) + (x − xi )2 .
2 (xi+1 − xi )

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 12

Une deuxième intégration donne :

S 00 (xi ) S 00 (xi+1 ) − S 00 (xi )

S (x) = fi + S 0 (xi ) (x − xi ) + (x − xi )2 + (x − xi )3 .
2 6 (xi+1 − xi )
A partir de ces définitions, on peut établir les relations suivantes :
· ¸
0 00 00 fi+1 − fi fi − fi−1
Si−1 hi−1 + (hi + hi−1 ) Si + Si+1 hi = 6 − .
hi hi−1
On a aussi pour i = 2, 3, · · · , n − 1
fi+1 − fi 00 hi hi
Si0 = − Si+1 − Si00
hi 6 3
où hi = xi+1 − xi et Si00 = S 00 (xi )
Les valeurs S100 et Sn00 sont déterminées par la quatrième propriété de l’interpolation Spline
à savoir
S 00 (a) = S 00 (b) = 0 ⇒ S100 = Sn00 = 0.

2.2.5 Autres types d’interpolation polynomiale

Dans les applications, d’autres renseignements peuvent être connus sur la fonction f
(par exemple la parité). Dans ce cas, le polynôme est écrit sous la forme d’une somme de
puissances paires de (si la fonction est paire) et de puissances impaires (si la fonction est
impaire).
Les coefficients du polynôme peuvent être déterminés comme dans le cas général. De
plus il peut aussi arriver qu’en plus des valeurs yk = f (xk ), nous disposons des valeurs
des dérivées yk0 = f 0 (xk ). Dans ces cas on utilise la formule d’interpolation de Hermite.
Il existe également d’autres types d’interpolation polynomiale telle que celle de Stirling
définie à partir des différences centrées sur base entière, et celle de Bessel définie à partir
des différences centrées sur base demie entière.
La formule de Stirling s’utilise surtout lorsque xp est proche d’un nœud xk tandis que
celle de Bessel est adaptée au cas où xp est proche de xk +x2 k+1 = xk+1/2 .

Remarque : L’opérateur
¡ ¢ de
¡ différence
¢ centrée sur base entière est définie par :
∆f (x) = f x + h2 − f x − h2 ⇒ ∆fk = fk+1/2 − fk−1/2

2.3 INTERPOLATION TRIGONOMETRIQUE, DOUBLE ET INVERSE

2.3.1 Interpolation trigonométrique

Quand la fonction représentée est périodique, on utilise l’interpolation trigonométrique. Ce
problème a été d’abord analysé par Gauss et ensuite par Hermite. La formule d’interpolation
de Hermite est la suivante :
n+1
,n+1 n+1
Y Y X
fk (x) = sin (x − xi ) sin (xk − xi ) et f (x) = fk (x) yk
i=1 i=1 k=1
i6=k i6=k

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 13

La différence entre la formule de Gauss et celle de Hermite réside dans le fait que chez
Gauss les arguments de la fonction sinus sont divisés par deux.

2.3.2 Interpolation double

a) Problème
Au cours de certaines études, l’on peut être confronté à un problème d’interpolation à
deux arguments, notamment dans le cas des problèmes plans où une grandeur physique
dépend des coordonnées x et y. Ce type de problème peut se résoudre de deux manières :
(i) On peut d’abord interpoler par rapport à une variable ensuite par rapport aux autres ;
(ii) On peut aussi trouver une formule d’interpolation qui utilise les différences doubles.
b) Différences doubles ou à deux voies
Soit la fonction Z = z ( x , y) dont la valeur en un point (xm , ye ) est zme .
On note ∆10 zme = ∆x zme = zm+1,e − zm,e , ∆01 zme = ∆y zme = zm,e+1 − zme .
Par exemple ∆10 = z2,1 − z11 , et ∆11 z11 = ∆2xy z11 = ∆10 (z12 − z11 ) = ∆01 z21 − ∆01 z11 .
De façon générale ∆αβ z11 = ∆α0 z1β + β∆α0 z1,β−1 + β(β−1)
2
∆α0 z1,β−2 + · · · + ∆0β z11 .
c) Formule d’interpolation double
Elle est donnée par

z (x h, y) = f (x , y) = z11 + x−x1
h
+ ∆10 z11 + y−y l
1
∆01 z11 + i
(x−x1 )(x−x2 ) 20 2(x−x1 )(y−y1 ) 11 (y−y1 )(y−y2 ) 02
+ 2!1 h2
∆ z 11 + hl
∆ z 11 + l2
∆ z 11 + · · ·
 (x−x1 )(x−x2 )+···+(x−xm−1 ) m0 
hm
∆ z11 + m(x−x1 )(x−x2h)···(x−x m−1 l
m−2 )(y−y1 )
∆(m−1)1 z11
1
+ m!  + m(m−1)(x−x1 )(x−x2 )···(x−x m−3 )(y−y1 )(y−y2 )
∆(m−2)2 z11 + · · · 
hm−2 l2
+ (y−y1 )(y−yl2m)···(y−ym−1 ) ∆0m z11

2.3.3 Interpolation inverse

Elle consiste à déterminer xp connaissant yp . On peut dans ce cas utiliser les formules
d’interpolation décrites précédemment, il suffira de remplacer y par x.

2.4 APPROXIMATION PAR LA METHODE DES MOINDRES CARRES

2.4.1 Position du problème

La méthode d’approximation par les moindres carrés est utilisée dans de nombreuses ap-
plications sous des noms différents : optimisation linéaire, méthode de régression , lissage des
courbes. Les expériences en physique peuvent contenir des milliers de mesures. L’utilisation
des méthodes d’interpolation va alors conduire à des grandes erreurs d’arrondi et des calculs
assez longs. La méthode des moindres carrés nous permet d’utiliser plusieurs fonctions pour
obtenir une approximation simple.
Par exemple, soit f1 , f2 , · · · , fm es valeurs correspondantes aux points x1 , x2 , · · · , xm et
soit f (x) = a0 + a1 x + a2 x2 + · · · + an xn un polynôme d’approximation de degré n avec
n ¿ m.
Il est clair que puisque n ¿ m il n’est pas possible d’avoir f telle que ∀k, f (xk ) = fk .
Mais on peut choisir f telle que f (xk ) ' fk . La fonction f sera bonne si la quantité (reste)

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 14

R = [f1 − f (x1 )]2 + [f2 − f (x2 )]2 + · · · + [fm − f (xm )]2

est la plus petite possible. Il s’agit donc de trouver la fonction f qui minimise le reste R.

2.4.2 Résolution du problème de l’approximation par les moindres carrés

Pour résoudre ce problème, on doit établir un système d’équations linéaires de la forme :
∂R
= 0; k = 0, 1, · · · , n.
∂ak
La solution de ce système permet alors de déterminer les coefficients ai de f .
Par exemple, considérons le cas polynomial n = 1 ⇒ f (x) = a0 + a1 x, le reste est
m
X
R= [fk − (a0 + a1 xk )]2
k=1

∂R ∂R
Pour que R soit minimale, il faut que ∂a 0
= 0 et ∂a 1
= 0.
On trouve
P 2P P P P P P
xk f k − xk xk f k m xk f k − xk f k
a0 = P P et a1 = P P .
m x2k − ( xk )2 m x2k − ( xk )2
En général on choisit f (x) sous la forme f (x) = a0 g0 (x) + a1 g1 (x) + · · · + an gn (x)
où les fonctions gk (x) peuvent être de natures différentes (polynôme rationnel, logarithme,
exponentiel, etc.)
Le système d’équations à résoudre est alors défini de la manière suivante :

 Pm

 A00 a0 + A01 a 1 + . . . + A 0n a n = g0 (xi ) fi




i=1
 A10 a0 + A11 a1 + · · · + A1n an = P g1 (xi ) fi
m

i=1

 ..

 .

 P
m


 An0 a0 + An1 a1 + · · · + Ann an = gn (xi ) fi
i=1
où m ½
X j = 0, 1, · · · , n
Akj = gk (xi ) gj (xi ) avec
k = 0, 1, · · · , n
i=1

EXERCICE 1 : Les résultats d’une expérience sont donnés dans le tableau suivant :
x 1 2 3 4 5
y 1.6 2 2.3 2.4 2.5
On se propose de déterminer l’expression y = f (x). Pour cela, on choisit a) y = a + xb ,
et b) y = a + xb + c ex . En utilisant la méthode des moindres carrés, calculer les coefficients
a, b et c.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTERPOLATION ET APPROXIMATION DES FONCTIONS 15

EXERCICE 2 : Du tableau suivant utiliser la formule d’interpolation de Lagrange pour

déterminer y pour x = 102 et x pour y = 13.5
x 93 96.2 100 104.2 108.7
y 11.38 12.80 14.7 17.07 19.91
Solution : ( xc = 102 et yc = 15, 79)
(xc = 97, 65 et yc = 13, 5).

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 Chapitre III

INTEGRATION ET DIFFERENTIATION

3.1 DERIVATION NUMERIQUE

3.1.1 Formules classiques

Soit f une fonction différentiable sur un intervalle [a, b]. Le nombre dérivé de f au point
x0 est défini par :

f (x + h) − f (x)
f 0 (x) = lim .
h→0 h
On peut également utiliser la formule

f (x + h/2) − f (x − h/2)
f 0 (x) = lim
h→0 h
De la même manière, la dérivée seconde est donnée par :

f (x + h) − 2f (x) + f (x − h)
f 00 (x) = lim
h→0 h2
Cependant, si h est trop petit, la machine commet des erreurs de calcul et d’arrondis dues
aux instabilités. Pour éviter ce genre de problèmes, il faut utiliser des formules plus précises.

3.1.2 Formules plus précises

Pour établir ces formules, de même que celles données précédemment, on utilise des po-
lynômes d’interpolation de la fonction f .
a) Dérivée à partir de l’interpolation linéaire
Dans l’interpolation linéaire, on a

x − x2 x − x1 1
f (x) = f (x1 ) + f (x2 ) + f 00 (ξ)(x − x1 )(x − x2 )
x1 − x2 x2 − x1 2
D’où

f (x2 ) − f (x1) 1 00 1 d
f 0 (x) = + f (ξ) [(x − x1 ) + (x − x2 )] + (x − x1 )(x − x2 ) [f 00 (ξ(x))]
x2 − x1 2 2 dx
L’erreur de troncature correspondant à la dérivée par interpolation linéaire est donc

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTEGRATION ET DIFFERENTIATION 17

1 1 d
E = f 00 (ξ) [(x − x1 ) + (x − x2 )] + (x − x1 )(x − x2 ) [f 00 (ξ(x))]
2 2 dx
avec ξ ∈ [x1 ; x2 ]
Si nous posons x = x1 et x2 = x1 + h, on obtient
f (x +h)−f (x
f 0 (x1 ) = 1
h
1)
avec une erreur de troncature E = 21 hf 00 (ξ) : C’est la formule
classique donnée précédemment.
b) Dérivée à partir de l’interpolation quadratique
En procédant comme dans le cas d’une interpolation linéaire, et en posant x = x2
on établit que f 0 (x) = f (x+h)−f
2h
(x−h)
avec une erreur de troncature E = 61 h2 f 000 (ξ) où
ξ ∈ [x − h; x + h].
Si l’on posait x = x1 , il viendrait le résultat suivant :
Si f est de classe C 4 sur l’intervalle [a, b] alors pour toutes valeurs de x et h telles que
a < x + 2h < b, il existe ξ ∈ [x − h; x + h] tel que

1 1
f 0 (x) = [−3f (x) + 4f (x + h) − f (x + 2h)] + h2 f 000 (ξ).
2h 3
c) Dérivée à partir de l’interpolation cubique
On obtient le résultat suivant : Si f est de classe C 6 sur l’intervalle [a, b] alors pour
toute valeur de x et h avec a < x − 2h < x + 2h < b , il existe ξ ∈ [x − 2h; x + 2h] tel
que

1 1 d £ (4) ¤
f 0 (x) = [f (x − 2h) − 8f (x − h) + 8f (x + h) − f (x + 2h)] + h4 f (ξ(x)) .
12h 30 dx

EXERCICE 1 : Soit la fonction f (x) = ex on prend h = 10−n pour n allant de 1 à 8.

Utiliser toutes les formules précédentes pour évaluer f 0 (0).

EXERCICE 2 : Etablir la formule de dérivation numérique à partir de l’interpolation qua-

dratique.

EXERCICE 3 : Au cours d’une expérience, les résultats suivants ont été obtenus :
x 1.0 1.10 1.20 1.30 1.40 1.50
y 1.6487 1.7333 1.8221 1.9155 2.0138 2.1170
Calculer les dérivées y 0 aux points x = 1.05; 1.15; 1.25; 1.35 et 1.45.

3.2 INTEGRATIONS NUMERIQUES- FORMULES DE NEWTON-COTES

Rb
L’intégrale définie I = a f (x)dx représente un nombre qui, dans certains cas favorables
peut être calculé analytiquement. Dans le cas contraire, on doit faire appel au calcul numé-
rique. Pour ce faire, la fonction f supposée continue sur [a, b] est remplacée par une formule
d’interpolation. Généralement, f est approximée par un polynôme et les formules d’intégra-
tion qui s’en déduisent sont dites formules de Newton-Cotes. Ces formules supposent que les
noeuds d’interpolation sont équidistants.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTEGRATION ET DIFFERENTIATION 18

3.2.1 Méthode des trapèzes :

Dans la méthode des trapèzes, l’interpolation de la fonction f est linéaire. En supposant
que [a, b] est court, on a
x−a x − b 1 00
f (x) = f (a) + f (b) + f (ξ)(x − a)(x − b).
b−a a−b 2
Il vient alors
b−a
I∼= [f (a) + f (b)] + E,
2
où E est l’erreur de troncature définie par
Z b
E= f 00 (ξ)(x − a)(x − b)dx.
a
Pour évaluer E, nous utilisons le théorème de la moyenne pour les intégrales.

Théorème de la moyenne : Soient f et g deux fonctions continues sur [a, b]. Si g(x) ne change
Rb Rb
pas de signe sur [a, b], alors il existe η ∈ [a, b] tel que a f (x)g(x)dx =f (η) a g(x)dx.
1 00
Puisque le produit (x−a)(x−b) reste négatif sur [a, b], nous aurons E = − 12 f (ξ)(b−a)3 .
Lorsque l’amplitude de l’intervalle [a, b] n’est pas assez petite, l’on utilise la propriété
d’additivité de l’intégrale pour faire un découpage en morceaux réduits, à savoir :
Z b n Z
X xk
f (x)dx = f (x)dx
a k=1 xk−1

où a = x0 < x1 < . . . < xn−1 < xn = b.

Prenons xk = x0 + kh, on obtient alors :
Z b
h
f (x)dx = [f (x0 ) + 2f (x1 ) + 2f (x2 ) + ... + 2f (xn−1 ) + f (xn )] + E
a 2
Pn
1 3 00
Où E = Ek , avec Ek = − 12 h f (ξk ) : C’est la formule composite de la méthode des
k=1
trapèzes.
P
n
n
X n
X f 00 (ξk )
1 3 00 nh3 k=1 nh3 00
Ek = − h f (ξk ) = − =− f (η)
k=1
12 k=1 12 n 12
h2
Où η ∈ [a, b] ; or b − a = nh d’où E − 12
(b − a)f 00 (η).
R2 dx
EXERCICE 4 : Calculer I = 1 x
pour h = 0.2 et pour h = 0.1. Comparer à la valeur
exacte.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTEGRATION ET DIFFERENTIATION 19

3.2.2 Algorithme :
Cet algorithme permet de calculer une intégrale par la méthode composite des trapèzes.
Il utilise des tableaux pour xk et fk = f (xk ).
1◦ - Définir la fonction f
2◦ - Donner a, b et n
b−a
3◦ - h ← n−1
,S ←0
◦
4 - pour k allant de 0 jusquà n faire
xk ← a + kh
fk ← f (xk )
Fin pour
5◦ - pour k allant de 1 jusquà n − 1 faire
S ← S + hfk
Fin pour
6◦ - Ecrire S

EXERCICE 5 : Ecrire un autre algorithme qui n’utilise pas les tableaux.

3.2.3 Méthode de Simpson

Dans la méthode de Simpson, la fonction f est approximée par un polynome de degré
deux qui passe par les points a, b et c = a+b2
. Cette méthode permet d’obtenir :
Z b
b−a
f (x)dx = [f (a) + 4f (c) + f (b)] + E
a 6
avec
Z
1 b 00
E= f (ξ)(x − a)(x − b)(x − c)dx
6 a
En utilisant l’additivité des intégrales on obtient :

Z b
h
f (x)dx = [f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) + 2f (x4 ) + ... + 4f (xn−1 ) + f (xn )].
a 3

EXERCICE 6 : Ecrire un algorithme pour la méthode de Simpson.

EXERCICE 7 : Etablir une formule d’intégration numérique qui utilise une interpolation
de degré 3 et qui interpole les points a, b, c = a+b
4
, et d = 3(a+b)
4
.

EXERCICE
R 8: Utiliser la méthode de Simpson pour calculer manuellement l’intégrale I =
1 4
0 1+x2
dx.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

INTEGRATION ET DIFFERENTIATION 20

3.3 INTEGRALES MIXTE ET DOUBLE

3.3.1 Intégrale mixte

Rb
C’est une intégrale de la forme I = a f (x)g(x)dx où l’expression mathématique de g(x)
est connue et où f (x) est approchée par une formule d’interpolation. Ce type d’intégrale
intervient par exemple dans le calcul des coefficients de Fourier du développement d’une
fonction : R RT
T
an = T2 0 f (x) cos(nx)dx, et bn = T2 0 f (x) sin(nx)dx. Dans ce cas, g(x) = cos(nx) ou
g(x) = cos(nx). Ce type d’intégrale intervient également
R1 dans la résolution des équations
intégrales. Par exemple l’équation f (x) = 5x + 0 (t + x)f (t)dt = 0 ou dans le calcul des
intégrales comportant des singularités.
Pour calculer ce type d’intégrale, on élabore des formules d’intégration de type Newton-
Côtes (f (x) est interpolée par un polynôme).

3.3.2 Interpolation double

R x=b R y=d
Il s’agit de calculer une intégrale de la forme I = x=a y=c Z(x, y)dxdy.
Pour ce faire, l’on peut utiliser la formule d’interpolation double de la fonction Z(x, y).

EXERCICE 9 : Etablir une formule générale pour une intégrale double en limitant l’inter-
polation à l’ordre deux.
R x=4.4 R y=2.6 dxdy
EXERCICE 10 : calculer I = x=4.0 y=2 xy
pour h = 0.2 et h = 0.3. Comparer à la
valeur exacte.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 Chapitre IV

RESOLUTION NUMERIQUE DES EQUATIONS

DIFFERENTIELLES ET EQUATIONS INTEGRALES

4.1 GENERALITES

4.1.1 Définition et classification

Une équation différentielle est une relation qui lie une fonction y à ses dérivées. S’il y a une
seule variable indépendante, on parle alors d’équation différentielle ordinaire. Dans le cas de
plusieurs variables indépendantes, on parle d’équations au dérivées partielles. Les équations
différentielles ordinaires de variables indépendante x peuvent se mettre sous la forme

F (x, y, y 0 , · · · , y (n−1) , y (n) ) = 0 (4.1)

où y (i) est la dérivée d’ordre i de y et n est l’ordre de l’équation différentielle et a ≤ x ≤ b.
Le degré d’une équation différentielle, si celle-ci peut être décrite comme un polynôme,
où les indéterminées sont des dérivées, est le degré de la dérivée de l’ordre le plus élevé. Par
exemple, l’équation (y 00 )3 + (y 0 )5 + 7y = ex est une équation différentielle de second ordre et
de degré 3. Si la fonction F est un polynôme dans laquelle les degrés de y et de ses dérivées
sont égales 0 ou à 1 et si F ne contient pas des produits de y et de ses dérivées, ou des
dérivées entre elles, alors l’équation différentielle est dite linéaire. Dans tous les autres cas,
l’équation est dite non linéaire.
La résolution d’une équation différentielle d’ordre n nécessite la connaissance préalable
de n constantes, notamment n conditions imposées à y et ses dérivées. Si l’on fixe la valeur
de y et de ses (n − 1) dérivées au point initial x = a, on a un problème aux valeurs initiales.
Si par contre les n conditions sont imposées les unes en x = a et les autres en x = b, on a
un problème aux valeurs aux limites. Origine des

4.1.2 Equations différentielles

Les équations différentielles sont d’une importance fondamentale dans toutes les disci-
plines scientifiques. Ceci est dû au fait que plusieurs lois et phénomènes ont une formulation
mathématique sous forme d’équations différentielle. Par exemple la seconde loi de New-
ton est une équation différentielle du second ordre ; la loi de la radioactivité est une
équation différentielle du premier ordre.

4.1.3 Position du problème Numérique

Considérons l’équation différentielle ordinaire du premier ordre

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 22

y 0 = f (x, y) (4.2)
Le problème est de trouver la valeur ȳ de y correspondant à x = x0 + h (où en xk =
x0 + kh, k = 1, 2, 3... et h un réel) connaissant y = y0 pour x = x0 . Dans certains cas,
un tel problème peut se résoudre en trouvant les solutions générales et particulières par
les méthodes analytiques classiques (intégration directe, facteur intégrant, transformée de
Laplace, etc.).
Lorsque aucune méthode analytique ne permet de calculer la solution générale, il faut
trouver les méthodes numériques pour approcher la valeur cherchée. En intégrant Eq.(4.2),
on obtient
Z x
y = y0 + f (x, y)dx (4.3)
x0

On a alors
Z x0 +h
ȳ = y0 + f (x, y)dx (4.4)
x0

On peut donc utiliser des méthodes numériques de calcul d’intégrales pour évaluer l’intégrale
ou faire un développement de y(x + h) en série de puissances de h.

4.2 METHODES POUR LES PROBLEMES A VALEURS INITIALES

4.2.1 Méthode de Picard

L’équation (4.3) ou (4.4) est compliquée par la présence de y dans l’intégrale. C’est une
équation intégrale. Pour résoudre une telle équation, Picard a utilisé une méthode d’approxi-
mations successives. Il obtient une première approximation y1 en remplaçant y par y0 dans
l’intégrale et obtient
Z x0 +h
y1 = y0 + f (x, y0 )dx (4.5)
x0

et l’intégrale peut être calculée par les méthodes numériques connus telles que la formule des
trapèzes ou celle de Simpson. Connaissant y1 , on obtient y2 en utilisant le même principe
Z x0 +h
y2 = y0 + f (x, y1 )dx (4.6)
x0

Ce processus peut ainsi être répété plusieurs fois. L’approximation d’ordre n étant donnée
par Z x0 +h
yn = y0 + f (x, yn−1 )dx (4.7)
x0

La méthode est acceptable lorsque les yn convergent vers une valeur fixe ȳ (valeur ap-
prochée). Cette méthode nécessite à chaque étape un calcul d’intégrale, ce qui n’est pas
toujours facile. De plus, dans la pratique, elle exige un temps de calcul long, surtout si l’on
veut trouver ȳ en plusieurs points.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 23

4.2.2 Séries de Taylor

Le développement de Taylor de y = g(x) au voisinage de (x0 , y0 ) est
1 1
y = g(x0 ) + (x − x0 )g 0 (x0 ) + (x − x0 )2 g 00 (x0 ) + (x − x0 )3 g 000 (x0 ) + · · ·
2 6
Or si
y = g(x) ⇒ y 0 = g 0 (x) = f (x, y)
il vient alors
µ ¶
00 00 df
0 ∂f ∂f ∂y ∂f ∂f
y = g (x) = f (x, y) = = + = +f
dx ∂x ∂y ∂x ∂x ∂y
d ∂ ∂
⇒ = +f
dx ∂x ∂y
De même
µ ¶µ ¶
000 ∂
000 ∂f ∂f ∂ ∂f ∂f
y = g (x) = +f+f +f
∂x ∂x ∂y ∂y ∂x ∂y
µ ¶ 2
∂ 2f ∂f ∂f ∂ 2f ∂f 2
2∂ f
= + + 2f + f + f
∂x2 ∂x ∂y ∂x∂y ∂y ∂y 2
Posons
∂f ∂f ∂ 2f ∂2f ∂ 2f
p= , q= , r= , s = et t =
∂x ∂y ∂x2 ∂x∂y ∂y 2
Soit f0 , p0 , q0 , r0 , s0 , t0 en (x0 , y0 ). Il vient alors d’après le développement en série de Taylor
avec x = x0 + h,
1 1 ¡ ¢
ȳ = y0 + hf0 + h2 (p0 + f0 q0 ) + h3 r0 + p0 q0 + 2f0 s0 + f0 q02 + f02 t0 + · · ·
2 6
Et la formule itérative s’écrit
1 1 ¡ ¢
ȳk+1 = yk + hfk + h2 (pk + fk qk ) + h3 rk + pk qk + 2fk sk + fk qk2 + fk2 tk .
2 6

4.2.3 Méthode d’Euler

Encore appelée méthode de la dérivée première, elle s’exprime par la formule

y (x + h) = y (x) + hf (x, y) .
Ici, le développement de Taylor se limite à l’ordre 2 et la formule itérative est donc :

yk+1 = yk + hf (xk , yk ) .
C’est la méthode la plus simple mais son inconvénient réside dans le fait qu’elle est souvent
très instable et exige des pas de calcul h très petits.

EXERCICE 1 : On considère l’équation y 0 = y 2 + x. Trouver (manuellement) la solution

approchée en x0 = 0.5 sachant que en x = 0, y = 1. On donne h = 0.1. Comparer à la valeur
exacte.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 24

A’ B’
a b

EXERCICE 2 : Ecrire un algorithme pour la méthode d’Euler.

4.2.4 Les méthodes de Runge et Kutta

Ce sont les méthodes élaborées par C. Runge en 1894 et améliorées par W. Kutta vers
1901. Elles utilisent les méthodes d’intégrations numériques des trapèzes et de Simpson. Ces
méthodes sont les plus utilisées car elles sont très stables.
a) Rappel sur les formules des trapèzes et de Simpson
Rb
Considérons l’intégrale a g (x)dx. Pour calculer cette intégrale sur l’intervalle [a, b] supposé
court, on peut remplacer la fonction g(x) par une interpolaion linéaire passant par les points
d’abcisses a et b. on peut également utiliser une interpolation quadratique de g(x) passant
par a, c = (a + b)/2 et b.
– La formule des trapèzes s’appuie sur une interpolation linéaire de g(x) aux points
(a, g(a)) et (b, g(b)). Pour cela g(x) est approximée par une droite affine de la forme
g(x) = αx + β. La formule d’interpolation de Lagrange donne g(x) = g(a) x−a b−a
+ g(b) x−b
a−b
.
Il vient alors Z b
b−a
g (x)dx ' [g(a) + g(b)].
a 2
Il s’agit ainsi de la surface du trapèze AA0 B 0 B.
– Dans le cas où on la fonction g(x) est approximée par un polynome du second degré,
on parle d’interpolation quadratique. Posons alors g(x) = αx2 + βx + γ et puisque g(x)
passe par les points (a, g(a)), (c = b+a2
, g(c)) et (b, g(b)), la formule d’interpolation de
Lagrange donne
(x − b)(x − c) (x − a)(x − c) (x − a)(x − b)
g(x) = g(a) + g(b) + g(c) .
(a − b)(a − c) (b − a)(b − c) (c − a)(c − b)
On en déduit la formule d’intégration de Simpson sous la forme
Z b
b−a
g (x)dx ' [g(a) + 4g(c) + g(b)].
a 6

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 25

b) Méthode de Runge-Kutta du second ordre

Elle est basée sur la formule intégrale des trapèzes. puisque
Z x0 +h
h
f (x, y)dx ' [f (x0 , y0 ) + f (x0 + h, y0 + hf (x0 , y0 ))],
x0 2

il vient y = y0 + h2 {f (x0 , y0 ) + f [x0 + h, y0 + hf (x0 , y0 )]} .

Et la formule itérative s’écrit

yk+1 = yk + h2 {f (xk , yk ) + f [xk + h, yk + hf (xk , yk )]}

Elle représente une erreur d’ordre o (h3 ).
c) Méthode de Runge du quatrième ordre
Supposons les valeurs y0 , y1 , y2 de y connues aux points x0 , x1 = x0 + h2 , x2 = x0 + h.
On a Z x0 +h
L = ȳ − y0 = f (x, y)dx
x0

En utilisant la formule d’intégration numérique de Simpson, on obtient

· ¸
h h
L= f (x0 , y0 ) + 4f (x0 + , y1 ) + f (x0 + h, y2 ) .
6 2
Les valeurs y1 et y1 étant inconnues, on utilise les approximations
h h
y1 ≈ y0 + f (x0 , y0 ) = y0 + f0 et y2 ≈ y0 + hf (x0 + h, y0 + hf0 )
2 2
Ces approximations sont faites de telle manière qu’en développant L en séries de puissance
de h, ces 3 premiers termes coïncident avec le développement en série de Taylor. Il vient alors

h
© £ ¤ ª
L≈ 6
f (x0 , y0 ) + 4f x0 + h2 , +y0 + h2 f (x0 , y0 ) + f [x0 + h, y0 + hf (x0 + h, y0 + hf (x0 , y0 ))]

Soit plus généralement yk+1 = yk + L avec

h
© £ ¤ ª
L= 6
f (xk , yk ) + 4f xk + h2 , +yk + h2 f (xk , yk ) + f [xk + h, yk + hf (xk + h, yk + hf (xk , yk ))]

Elle présente une erreur d’ordre o (h5 ). Pratiquement on fait le calcul de la manière sui-
vante : On calcule les quantités Li définies par

L1 = hf (xk , yk )
L2 = hf (xk + h, yk + L1 )
L3 = hf (x
¡ k + h, yk + L2 )¢
L4 = hf xk + h2 , yk + L21
et par la suite on a :

yk+1 = yk + 16 (L1 + 4L4 + L3 )

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 26

EXERCICE 3 : Etablir l’algorithme de cette méthode et la traduire en Pascal et en Fortran.

d) Méthode de Runge-Kutta ou de Kutta-Simpson du quatrième ordre
C’est une amélioration de la méthode de Runge. Elle présente également une erreur de
l’ordre o (h5 ) et son schéma itératif s’écrit :

yk+1 = yk + 16 (L1 + 2L2 + 2L3 + L4 )

avec

L1 = hf (x
¡ k , yk )h ¢
L2 = hf ¡xk + 2 , yk + L21 ¢
L3 = hf xk + h2 , yk + L22
L4 = hf (xk + h, yk + L3 )
C’est la méthode la plus utilisée. Il existe d’autres variantes (de la méthode ) plus élaborées
telles la méthode de Kutta Merdsen ou de Kutta-Simpson du sixième ordre.

4.2.5 Méthodes à pas multiples ou d’Adams Bashforth

a) Formule du second ordre
Soit l’équation différentielle y 0 = f (x, y), on veut connaître y (x + h) connaissant y (x) et
y (x − h). Pour cela , on fait le développement suivant :
h2 00
¡ ¢ h2 00
y (x + h) = y (x) + hy 0 + 2
y + o h3 ⇒ yk+1 = yk + hyk0 + y
2 k
Or
yk0 = f (xk , yk ) 0
et yk−1 = yk0 − hyk00 .
Il vient alors
yk0 − yk−1
0
yk00 = .
h
Finalement, on obtient
£ ¤
yk+1 = yk + h
2
3yk0 − yk−1
0
⇒ yk+1 = yk + h2 [3f (xk , yk ) − f (xk−1 , yk−1 )] .
C’est la formule d’Adams - Bashforth du second ordre présentant une erreur d’ordre o(h3 ).
Pratiquement au début des itérations, il faut connaître les valeurs de y en x0 et en x1 afin de
calculer la première valeur y2 de l’itération. Pour avoir y1 , on peut utiliser l’une des formules
approchées classiques. Par exemple en utilisant la formule d’intégration des trapèzes, on
obtient :

y1 = y0 + h2 [f (x0 , y0 ) + f [x0 + h, y0 + hf (x0 , y0 )]] .

La formule itérative peut alors débuter par k = 1, et on calcule les y2 , y3 , y4 , . . . Par
exemple,

y2 = y1 + h2 [3f (x1 , y1 ) − f (x0 , y0 )] .

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 27

EXERCICE 4 : Ecrire l’algorithme de cette méthode et la traduire en FORTRAN et en

PASCAL.
b) Formule d’ordre supérieur
Pour obtenir des formules d’ordre supérieur présentant par exemple des erreurs d’ordre
o (h4 ) ou o (h5 ), il suffit de développer à l’ordre supérieur y (x + h). Par exemple à l’ordre 3,
on a
h2 00 h3 000
¡ ¢
yk+1 = yk + hyk0 + y
2 k
+ y
6 k
+ o h4 (?)
De même à l’ordre 4, on a
h2 00 h3 000 h4 (4)
¡ ¢
yk+1 = yk + hyk0 + y
2 k
+ y
6 k
+ 24
y + o h5 (??)
A partir de la formule (?) on obtient le schéma itératif suivant :
£ ¤
yk+1 = yk + h
12
23yk0 − 16yk−1
0 0
+ 5yk−2
soit
h
yk+1 = yk + 12
[23f (xk , yk ) − 16f (xk−1 , yk−1 ) + 5f (xk−2 , yk−2 )] .
Ainsi pour calculer yk+1 , il nous faut connaître les valeurs de yk , yk−1 et yk−1 . Donc, au
début des itérations, il faut avoir y0 , y1 et y2 . Les valeurs de y1 et y2 doivent être calculées
par les méthodes classiques.
A partir de la formule (??), on établit :

h
yk+1 = yk + 24
{55f (xk , yk ) − 59f (xk−1 , yk−1 ) + 37f (xk−2 , yk−2 ) − 9f (xk−3 , yk−3 )}

Ici il nous faut connaître y0 , y1 , y2 et y3 afin de lancer le processus itératif. On utilise les
formules classiques pour déterminer les valeurs de y0 , y1 , et y2 .

4.2.6 Méthode de Prédicteur-Correcteur

Les méthodes de Runge, Runge-Kutta et de Kutta-Simpson sont fondamentalement im-
plicites car l’inconnue yk+1 est exprimée à l’aide d’une fonction contenant yk+1 . Par exemple,
dans la formule de Runge-Kutta 2 (RK2), on a normalement

yk+1 = yk + h2 [f (xk , yk ) + f (xk+1 , yk+1 )]

Afin d’avoir une expression explicite, on part des approximations pour exprimer yk+1 du
second membre de l’égalité en fonction de yk . On pourrait aussi partir des formules d’Adams-
Bashforth pour évaluer (prédire) les valeurs de yk+1 , puis insérer ces valeurs prédites dans le
second membre de la formule implicite afin d’avoir une valeur corrigée. Dans un tel processus,
la formule explicite est appelée prédicteur tandis que la formule implicite utilisée est appelée
correcteur. Par exemple une méthode de prédicteur-correcteur peut utiliser les formules
suivantes :
yk+1 = yk + h2 [3f (xk , yk ) − f (xk−1 , yk−1 )] = Prédicteur (Formule d’Adams-Bashforth du
second ordre ).
yk+1 = yk + h2 [f (xk , yk ) + f (xk+1 , yk+1 )] = Correcteur (RK2)

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 28

EXERCICE 5 : Etablir l’algorithme de la méthode de prédicteur-correcteur, utilisant comme

prédicteur la formule d’Adams-Bashforth du second ordre et comme correcteur la formule de
RK2 et traduire cet algorithme en Pascal et en Fortran.

Remarque : On peut également combiner les autres formules implicites de Runge-Kutta

à d’autres formules explicites d’Adams-Bashforth pour avoir des méthodes de prédicteur-
correcteur plus élaborées.

4.3 METHODE POUR LES PROBLEMES AUX VALEURS AUX LIMITES

Il s’agit des problèmes de la forme

¡ ¢
y (n) = f x, y, y 0 , ..., y (n−1) avec a < x < b
avec y(a) = y1 et y(b) = y2 en plus d’autres conditions sur les dérivées de y. On peut
résoudre ce type de problème par plusieurs méthodes dont la méthode de tir, la méthode des
fonctions complémentaires ou de superposition et la méthode des différences finies.

4.3.1 Méthode de tir

Elle est utilisée dans le cas des problèmes aux limites non linéaires. Il s’agit de transformer
un problème aux valeurs aux limites en un problème aux valeurs initiales par un choix adé-
quat des valeurs initiales. Ce choix se fait à l’une des bornes du domaine d’intégration et doit
être tel que la solution passe le plus près possible à l’autre borne du domaine d’intégration.
Par exemple, considérons une équation du second ordre
½ 00
y = f (x, y, y 0 )
y(a) = η1 , y(b) = η2

La méthode de tir consiste à choisir γ telle que y 0 (a) = γ et de résoudre le problème aux
valeurs initiales
½ 00
y = f (x, y, y 0 )
y(a) = η1 , y 0 (a) = γ
Le choix de γ est bon lorsque en x = b, on a y(b) ≈ η2 , sinon il faut choisir une autre
valeur pour γ.
Une procédure à utiliser pour le choix de γ est la suivante : on sait que la valeur y(b)
dépend de γ, c’est à dire y(b) = y(b, γ). Pour un bon choix de γ, on doit avoir y(b, γ) − η2 ≈
0 ⇔ f (γ) = 0 où f (γ) = y(b, γ) − η2 . Or d’après la méthode de Newton, le schéma itératif
donnant la solution de cette équation est
f (γn ) y(b, γn ) − η2 (γn − γn−1 )(y(b, γn ) − η2 )
γn+1 = γn − = γ n − = γ n − .
f 0 (γn ) y 0 (b, γn ) y(b, γn ) − y(b, γn−1 )
Ainsi en pratique, on commence par donner une valeur initiale γ0 à γ, ensuite, on choisit
la valeur suivante γ1 , les autres valeurs de γ i.e (γ2 , γ3 , . . . ) sont ensuite évaluées à partir
de la formule ci-dessus.

EXERCICE 6 : Etablir l’algorithme de la méthode de tir.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 29

4.3.2 Méthode des fonctions complémentaires ou de superposition

C’est une méthode qui transforme une équation différentielle linéaire aux valeurs aux
limites en deux problèmes différentiels aux valeurs initiales. Elle a plusieurs variantes que
nous exposerons sur deux équations différentielles de second ordre.
a) Première variante
Considérons l’équation
½
y 00 = p (x) y 0 + q (x) y + r (x)
y (a) = η1 , y (b) = η2
Supposons que :

y (x) = y1 (x) + µy2 (x) , (µ = cte 6= 0)

Il vient
½
y100 = p (x) y10 + q (x) y1 + r (x)
y200 = p (x) y20 + q (x) y2
En x = a on a
y1 (a) + µy2 (a) = η1 .
Imposons que
y1 (a) = η1 ⇒ y2 (a) = 0
En x = b, on a
η2 − y1 (b)
y1 (b) + µy2 (b) = η2 ⇒ µ =
y2 (b)
Fixons y10 (a) = 0 et y20 (a) = 1, on obtient les deux problèmes aux valeurs initiales sui-
vantes :
 ½ 00

 y1 = p (x) y10 + q (x) y1 + r (x)
 0
½ y100 (a) = η1 0 y1 (a) = 0 x ∈ ]a, b].

 y2 = p (x) y2 + q (x) y2

y2 (a) = 0 y20 (a) = 1
On résout alors numériquement ces deux problèmes aux valeurs initiales pour x ∈]a, b] et
on obtient alors les valeurs de y1 (b) et y2 (b), d’où la valeur de µ. Puisque les valeurs de y1
et y2 sont connues en tout point de ]a, b], on déduit celles de y en tout point x ∈]a, b] par la
relation y (x) = y1 (x) + µy2 (x).
b) Deuxième variante
Considérons le problème suivant

y 00 = p (x) y + q (x) (a)

y 0 (a) = α00 y (a) + α10 (b)
y 0 (b) = β00 y (b) + β10 (c)
Considérons également l’équation

y 0 = α0 (x) y + α1 (x) (d)

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 30

Où nous avons choisi α0 (x) et α1 (x) tel que y satisfait toujours l’équation (a). En diffé-
renciant l’équation (d) par rapport à x et en identifiant à l’équation (a), on obtient
½ 0
α0 + α02 = p (x) (e)
0
α1 + α0 α1 = q (x) (f )
avec α0 (a) = α00 et α1 (a) = α10 .
Les équations (e) et (f ) sont des équations différentielles du premier ordre que l’on peut
résoudre analytiquement ou numériquement pour trouver α0 (x) et α1 (x).
De l’équation (d) on a
y 0 (b) = α0 (b) y (b) + α1 (b) = β00 y (b) + β10

(
y (b) = βα10 −α1 (b)
0 (b)−β00
(g)
⇒
y 0 (b) = β00 αβ100
(b)−β10 α0 (b)
−α0 (b)
(h)
L’équation (a) est ainsi transformée en trois problèmes auxvaleurs initiales. On peut
en effet résoudre le problème y 00 = p (x) y + q (x) avec les conditions initiales (g) et (h)
ou résoudre le problème y 0 = α0 (x) y + α1 (x) avec la condition initiale (g) obtenue après
résolution de (e) et (f ).

4.3.3 Méthode des différences finies

Elle convertit le problème de résolution d’une équation différentielle à la résolution des
équations algébriques. Soit le problème :
½ 00
y = p (x) y 0 + q (x) y + r (x)
y (a) = η1 , y (b) = η2
on a : y 0 (x) = y(x+h)−y(x)
h
et y 00 (x) = y(x+h)−2y(x)+y(x−h)
h2
.
Soit x = xk , on pose

y (xk ) = yk , p (xk ) = pk , q (xk ) = qk , r (xk ) = rk .

Par substitution des expressions des dérivées dans l’équation différentielle, on aboutit à
l’équation discrète
¡ ¢
(1 − hpk ) yk+1 − 2 − hpk − h2 qk yk + yk−1 = h2 rk ; k = 0, 1, . . . , N − 1
pour k = 1, 2, .., N − 1 où h = b−a
N
. On aboutit ainsi à un système d’équations algébriques
linéaires de la forme

 y0 = η1



 (1 − hp1 ) y2 − (2 − hp1 − h2 q1 ) y1 + y0 = h2 r1

 (1 − hp2 ) y3 − (2 − hp2 − h2 q2 ) y2 + y1 = h2 r2
.. .. ..

 . . .


 2 2
 (1 − hpN −1 ) yN − (2 − hpN −1 − h qN −1 ) yN −1 + yN −2 = h rN −1

yN = η2
Ce système algébrique peut être résolu par les méthodes itératives de Jacobi ou de Gauss-
Seidel.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 31

4.4 LES EQUATIONS INTEGRALES

4.4.1 Définition
On appelle équation intégrale une équation dont la fonction inconnue se trouve dans le
symbole d’intégration, par exemple :
Z x
x
y (x) = e sin x + 2 cos (x − t)y (t) dt
0
Une équation intégrale est dite linéaire lorsque la fonction inconnue y intervient avec un
degré 1. On a aussi les équations intégro-différentielles, par exemple
Z x
0
y (x) = 3y (x) + 2x + sin (x − t)y (t) dt
0
Sous la forme générale, les équations intégrales linéaires se présentent de la manière sui-
vante
Z x
y (x) = f (x) + k (x , t)y (t) dt
0
La fonction k (x , t) est appelée noyau de l’équation. Dans les cas simples, les équations
intégrales peuvent se résoudre analytiquement en utilisant par exemple la méthode de la
fonction de Gauss, des transformées de Laplace ou de Mellin, ou par les approximations
successives. Dans le cas contraire, il faut utiliser les méthodes numériques.

4.4.2 Résolution numérique des équations intégrales linéaires

Soit l’équation intégrale linéaire
Z b
y (x) = f (x) + k (x, t) y (t) dt (i)
a
L’intégrale définie dans (i) peut être remplacée par les formules de quadrature classique
(trapèze, Simpson). On obtient alors

y (x) = f (x)+(b − a) [c1 k (x, t1 ) y (t1 ) + c2 k (x, t2 ) y (t2 ) + · · · + cn k (x, tn ) y (tn )] (ii)
où t1 , t2 ,. . . , tn sont les nœuds de subdivision de l’intervalle [ a , b] et les valeurs des
coefficients ci dépendent de la formule d’intégration utilisée. Puisque l’équation (ii) est vé-
rifiée pour tout x ∈ [ a , b], elle est également vérifiée pour x = t1 , x = t2 , · · · , x = tn . Par
conséquent, à partir de l’équation (ii), on obtient un système d’équations de la forme

y (ti ) = f (ti )+(b − a) [c1 k (ti , t1 ) y (t1 ) + c2 k (ti , t2 ) y (t2 ) + · · · + cn k (ti , tn ) y (tn )] (iii)
Pour i allant de 1 à n. Posons y (ti ) = yi et f (ti ) = fi . le système (iii) devient alors

 y1 = f1 + (b − a) [c1 k (t1 , t1 ) y1 + c2 k (t1 , t2 ) y2 + · · · + cn k (t1 , tn ) yn ]

 y2 = f2 + (b − a) [c1 k (t2 , t1 ) y1 + c2 k (t2 , t2 ) y2 + · · · + cn k (t2 , tn ) yn ]
.. .. .. ..

 . . . .

yn = fn + (b − a) [c1 k (tn , t1 ) y1 + c2 k (tn , t2 ) y2 + . . . + cn k (tn , tn ) yn ]

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DIFFERENTIELLES ET EQUATIONS INTEGRALES 32

C’est un système de n équations à n inconnues que sont les yi . La résolution de ce système

permet de déterminer les yi et en les substituant dans l’équation (ii), on obtient l’expression
de y (x), solution de l’équation (i).

Exemple : Soit l’équation

Z 1
5x 1 1
y (x) = 9
− +
9 3
(t + x) y (t)dt
0
Pour résoudre numériquement cette équation, utilisons la méthode de Simpson et consi-
dérons pour cela trois points : t1 = 0,, t2 = 0.5, et t3 = 1.
Le système algébrique à résoudre est alors défini par :

5x
y (x) = 6
− 19 + 13 × 31 × 12 [(t1 + x) y (t1 ) + 4 (t2 + x) y (t2 ) + (t3 + x) y (t3 )] (?)

Soit en remplaçant x par ti ,


 y1 = − 91 + 18 1
¡(2y2 + y3 ) ¢ (t = 0.0)
y2 = 12 + 18 y21 + 4y2 + 32 y3
5 1
(t = 0.5)

y3 = 65 + 18
1
(y1 + 6y2 + 2y3 ) (t = 1.0)
On trouve y1 = 0, y2 = 12 et y3 = 1. Il vient alors y (x) = x (en remplaçant les yi par leurs
valeurs dans l’expression (?)).
R
1 1
EXERCICE 7 : Soit l’équation intégrale y (x) = 5x
9
− 1
9
+ 3 0
(t + x) y (t)dt. Résoudre cette
équation en utilisant la formule des trapèzes sur les trois points t1 = 0, t2 = 0.5 et t3 = 1

EXERCICE 8 : Discrétiser l’équation intégo-diff ’erentielle suivante :

Z x
0
y = 2y + f (x) + (k(x, t)) y (t)dt.
0

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 Chapitre V

SYSTEMES D’EQUATIONS ALGEBRIQUES

LINEAIRES ET NON LINEAIRES

5.1 METHODES DIRECTES POUR LES SYSTEMES D’EQUATIONS LINEAIRES

5.1.1 Rappel de la méthode de Gauss

Soit le système algébrique de la forme :

 a11 x1 + a12 x2 + · · · + a1n xn = b1



 a21 x1 + a22 x2 + · · · + a2n xn = b2
a31 x1 + a32 x2 + · · · + a3n xn = b3

 .
. .. .. ..

 .
 . . .
an1 x1 + an2 x2 + · · · + ann xn = bn
La méthode d’élimination de Gauss est basée sur le principe suivant. Les inconnues sont
éliminées successivement en exprimant à partir d’une des équations une inconnue en fonction
des autres puis en substituant dans le reste des équations. On part ainsi d’un système de n
équations à n inconnues à un système de n − 1 équations à n − 1 inconnues. Le processus
se poursuit jusqu’à l’obtention d’une équation à une inconnue. On trouve cette dernière
inconnue et par un processus inverse, les autres inconnues sont calculées. A chaque étape,
l’équation qui permet d’exprimer une inconnue en fonction des autres est appelée équation
pivotale. Pour éviter les erreurs d’arrondis, l’équation pivotale doit être celle qui possède le
plus grand coefficient de l’inconnue à éliminer. L’algorithme de cette méthode peut alors se
formuler de la manière suivante :
Pour une ligne numéro i, notée Li on fait l’opération
ai1
Li ← Li − L1 ; i = 2, 3, · · · , n
a11
On obtient alors un nouveau système de la forme

 a x + a12 x2 + · · · + a1n xn = b1
 11 1


 0 + a022 x2 + · · · + a02n xn = b02
0 + a032 x2 + · · · + a03n xn = b03

 . .. .. ..
 ..
 . . .

0 + an2 x2 + · · · + ann xn = b0n
0 0

Les coefficients a0ij et b0i sont donnés par

ai1 ai1
a0ij = aij − aij et b0i = bi − b1
a11 a11

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 34

On continue le processus en considérant le système constitué des lignes 2 à n et ainsi de

suite jusqu’au système constitué des lignes n − 1 et n. Ainsi pour une étape k, c’est-à-dire
en considérant le système constitué des lignes k à n, l’algorithme suivant est utilisé :
Pour i allant de k + 1 à n faire
bi ← bi − aakk
ik
bk
Pour j allant de k à n faire
aij ← aij − aakk
ik
akj
Fin pour
Fin pour
On peut ainsi en déduire l’algorithme de triangularisation du système d’équations en
faisant varier k de 1 à n − 1.
Finalement on obtient un système algébrique triangulaire supérieur de la forme

 a11 x1 + a12 x2 + · · · + a1n xn = b1

 0 + a22 x2 + · · · + a2n xn = b2
.. .. .. ..

 . . . .

0 +0 + · · · + ann xn = bn
que l’on peut resoudre de la manière suivante :
La dernière équation donne
bn
xn =
ann
Pour la racine d’ordre k, on peut établir que

xk = (bk − ak,k+1 xk+1 − ak,k+2 xk+2 − ... − ak,n xn )/akk

ou sous forme condensée
à n
!,
X
xk = bk − ak,l xl akk
l=k+1

EXERCICE 1 : Ecrire l’algorithme complet de la méthode de Gauss. Introduire le cas où il

existe des akk = 0 et faire le tri pour que akk soit maximal (C’est la règle du pivot maximal).
Compléter cet algorithme par une procédure de résolution du système triangulaire obtenu.

5.1.2 Méthode de Jordan

Elle est analogue à la méthode de Gauss, mais ici on met aussi à zéro les termes qui se
trouvent au dessus de la diagonale principale. On obtient alors un système diagonal.

5.1.3 Méthode de la décomposition triangulaire

a) Principe
Soit le système
AX = b,

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 35

on peut de manière unique décomposer la matrice A en deux matrices triangulaires dont

l’une est triangulaire inférieure avec les éléments quelconques sur la diagonale et l’autre
triangulaire suérieure avec tous les éléments égaux à 1 sur la diagonale.
Soientt L et S ces deux matrices triangulaires, on a

A = LS
et le système d’équations devient

SX = L−1 b.
Ce dernier système est un système triangulaire que l’on peut résoudre par le schéma
précédent.
b) Décomposition de la matrice
On peut écrire L et S sous la forme
   
l11 0 0 . . . 0 1 s12 s13 . . . s1n
 l21 l22 0 . . . 0   0 1 s23 . . . s2n 
   
 .   . . 1 s3n 
L= ; S= 
 .   . . . 
 .   . . . 
ln1 ln2 ln3 . . . lnn 0 0 0 . 1
L’équation LS = A entraîne par identification les relations suivantes

½
aj1 = lj1 1≤j≤n
a1j = l11 s1j 2≤j≤n
½
aj2 = lj1 s12 + lj2 2≤j≤n
a2j = l21 s1j + l22 s2j 3≤j≤n
½
aj3 = lj1 s13 + lj2 s23 + lj3 3≤j≤n
a3j = l31 s1j + l32 s2j + l33 s3j 4≤j≤n

.. .. ..
½ . . .
aj,n−1 = lj1 s1,n−1 + lj2 s2,n−1 + ... + lj,n−2 sn−2,n−1 + lj,n−1 n−1≤j ≤n
an−1,j = ln−1,1 s1j + ln−1,2 s2,j + ... + ln−1,n−1 sn−1,j j=n

ann = ln1 s1n + ln2 s2n + ... + lnn

Avec ces équations on obtient les éléments de la première colonne de L puis ceux de la
première ligne de S, ensuite ceux de la deuxième colonne de L et ceux de la deuxième ligne
de S et ainsi de suite. On obtient alors les éléments suivants qui permettent de calculer les
éléments lij et sij . On a :

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 36

½
lj1 = aj1
½ s1,j = a1j /l11
lj2 = aj2 − lj1 s12
½ 2j = (a2j − l21 s1j ) /l22
s
lj3 = aj3 − lj1 s13 − lj2 s23
s3j = (a3j − l31 s1j − l32 s2j ) /l33
.. .. ..
 . . .
 P
n−2

 lj,n−1 = aj,n−1 − ljk sk,n−1
· k=1 ¸

 P
n−2
 sn−1,j = an−1,j − ln−1,k skj /ln−1,n−1
k=1
P
n−1
lnn = ann − lnk skn
k=1

EXERCICE 2 : Ecrire l’algorithme permettant de calculer les lij et sij .

5.1.4 Méthode de la racine carré de Cholesky

Une matrice A symétrique peut être mise sous la forme d’un produit de deux matrices
M et M 0 dont l’une est la transposée de l’autre. On a ainsi une seule matrice triangulaire
supérieure M à calculer. Ceci simplifie les calculs et diminue l’encombrement mémoire de
la machine. En utilisant le procédé précédent on établit que les éléments de M sont définis
par :
 √
 m11 = a11



 m1j = ap 1j /m11



 m22 = a22 − m212



 m = a2j − m12 m1j /m22
 . 2j


 ..

 s

P 2
p−1
 m pp = a pp − mkp

 k=1

 Á

 P
p−1

 mpj = apj − mkp mkj mpp



 k=1

 ..

 .

 q
 m = a − m2 − m2 − ... − m2
nn nn 1n 2n n−1,n

5.1.5 Fragmentation des systèmes volumineux

Quand n est très grand, il n’y a pas assez de mémoire pour stocker les n(n + 1) coeffi-
cients du système et mettre en place le programme de calcul. On fragmente alors le système
d’équation en décomposant les matrices du système en sous-matrices. Par exemple on mettra
le système AX = b sous la forme

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 37

¯ ¯¯ ¯ ¯ ¯
¯ A1 . A2 ¯ ¯ X1 ¯ ¯ B1 ¯
¯ ¯¯ ¯ ¯ ¯
¯ . . . ¯¯ . ¯ = ¯ . ¯
¯ ¯¯ ¯ ¯ ¯
¯ A3 . A4 ¯ ¯ X2 ¯ ¯ B2 ¯
Les matrices A1 et A2 ne sont pas nécessairement des matrices carrées et la décomposi-
tion par les droites en pointillés est telle que les multiplications matricielles suivantes sont
possibles.
½
A1 X1 + A2 X2 = B1
A3 X1 + A4 X2 = B2
Si A1 est régulière, alors on a

X1 = A−1
1 [B1 − A2 X2 ]

En substituant dans la deuxième équation on obtient

£ ¤ £ ¤
A3 A−1 −1
1 A2 − A4 X2 = A3 A1 B1 − B2

On peut ainsi calculer les éléments de X2 et les substituer dans l’expression de X1 pour
avoir les éléments de X1 , donc toutes les composantes de X. Cette méthode peut être géné-
ralisée en décomposant la matrice en 9, 16, 25,... sous matrices.

5.1.6 Systèmes à coefficients complexes

Soit le système AZ = b avec A = M + iN , b = c + id et Z = X + iY .
Le système algébrique peut alors se décomposer en deux sous systèmes de la forme
½
MX − NY = c
NX + MY = d
On peut alors réunir ces deux sous systèmes en un système unique, ou utiliser l’approche
du paragraphe précédent.

5.1.7 METHODES ITERATIVES

Les méthodes directes telles que la méthode d’élimination de Gauss sont appropriées pour
les systèmes de taille moyenne. Pour les systèmes de grande taille il est conseillé d’utiliser les
méthodes itératives qui sont plus rapide, plus facile à programmer et moins encombrantes.
Cependant, ces méthodes sont basées sur le calcul des suites infinies qu’il faudra tronquer
à partir d’un certain rang. D’autre part, aucune méthode itérative ne donne des résultats
assez précis. Il faut en plus que chacune des équations du système possède un coefficient
assez grand devant tous les autres. La résolution se fait alors en exprimant l’inconnu ayant
le plus grand coefficient en fonction des autres inconnues. Pour le choix des valeurs initiales
des itérations, on fixe généralement les inconnues xi = 1 ou xi = 0 pour i = 1 jusqu’à n.
Pour la condition d’arrêt de l’itération, il y a plusieurs méthodes. Elles consistent à arrêter
les itérations lorsque toutes les quantités¯ ¯
¯ ¯ ¯ x(k+1) −x(k) ¯
¯ (k+1) (k) ¯
di = ¯xi − xi ¯ < ε ou aussi di = ¯¯ (k+1) ¯¯ < ε
i i
x i
Cependant le meilleur test de convergence est d’utiliser la condition

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 38

Pn ¯ ¯
¯ (k+1) (k) ¯
¯xi − xi ¯
i=1
d= Pn ¯ ¯ < ε.
¯ (k+1) ¯
¯xi ¯
i=1

Il existe deux variantes de méthode itérative : celle de Jacobi et celle de Gauss-Seidel.

Le principe des méthodes itératives consiste à exprimer l’inconnue xi de la ligne i en
fonction des autres inconnues, c’est-à-dire
 
n
,
 X 
xi = bi − aij xj 
 aii
j=1
j6=i

Ceci suggère que les valeurs connues x1 , x2 , ..., xi−1 , xi+1 , xi+2 , ..., xn soient substituées
pour trouver xi.
La formule itérative de Jacobi est :
 
n
,
 X 
= aij xj 
(k+1) (k)
xi bi −  aii
j=1
j6=i

(k)
où xj , sont les valeurs des inconnues après k itérations.
(k+1)
Dans la formule de Gauss-Seidel, les valeurs xj pour j < i − 1 sont utilisées aussitôt
(k+1)
qu’elles sont déterminées pour calculer xi . La formule itérative de Gauss-Seidel se met
alors sous la forme :
à i−1 n
!,
X (k+1)
X (k)
x(k+1)
i
= bi − aij xj − aij xj aii
j=1 j=i+1

Notons que les deux formules sont importantes car il existe des systèmes d’équations pour
lesquels la méthode de Jacobi converge tandis que la méthode de Gauss-Seidel ne converge
pas. Dans d’autres cas c’est la méthode de Gauss-Seidel qui converge et celle de Jacobi ne
converge pas.

5.2 RAPPELS SUR LES METHODES DE NEWTON ET DE BISSECTION POUR LES

EQUATIONS ALGEBRIQUES NON LINEAIRES

5.2.1 Méthode de Newton-Raphson(1960)

La formule a été donnée par Raphson en 1698 ; mais Newton avait suggéré une méthode
plus voisine quelques mois auparavant. Il se pourrait aussi que cette formule ait été utilisée
par Héron en l’an 100 Av J. C., pour trouver la racine carrée d’un nombre. Cette méthode
est de la classe des méthodes itératives par approximation.
a) La méthode

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 39

On commence généralement par trouver une solution grossière x0 ou approchée dite ap-
proximation d’ordre zéro. Dans le cas d’un problème de physique, une telle solution peut
être suggérée par la signification du problème ou par une courbe. Généralement, on localise
un intervalle [a, b] tel que x0 appartienne à cet intervalle. Désignons la solution exacte par
x̄ = x0 + h. A l’aide de la formule de Taylor, on obtient,

f (x0 + h) ≈ f (x0 ) + hf 0 (x0 ) + O(h2 ).

Or f (x0 + h) = 0 ⇒ f (x0 ) + hf 0 (x0 ) ≈ 0 et par conséquent
f (x0 )
h≈−
f 0 (x0 )
d’où
f (x0 )
x̄ = x0 −
f 0 (x0 )
Posons x̄ = x1 , on obtient x1 = x0 − ff0(x 0)
(x0 )
. C’est une première approximation de la solution
(l’indice désignant l’ordre de l’approximation). Procédons de même avec x1 , on obtient la
seconde approximation, et ainsi de suite :

f (x1 ) f (x2 ) f (xn )

x2 = x1 − , x 3 = x2 − , ··· , xn+1 = xn −
f 0 (x1 ) f 0 (x2 ) f 0 (xn )
En négligeant les termes d’ordre supérieur, on a remplacé la courbe représentative de
f (x) par une tangente en x = x0 , la méthode consiste géométriquement à tracer de proche
en proche les tangentes à la courbe et à rechercher leur point d’intersection avec l’axe des
abscisses. En effet, on a l’équation

f (x̄) = y = f (x0 ) + hf 0 (x0 ) = f (x0 ) + (x̄ − x0 )f 0 (x0 )

qui est l’équation de la tangente au point (x0 , f (x0 )).
Cette interprétation graphique fait que la méthode de Newton est souvent appelée mé-
thode de la tangente.
b) Remarques :
? Lorsque f (x0 )f 00 (x0 ) > 0, la suite xn+1 = xn − ff0(x n)
(xn )
est monotone, bornée et converge vers
x̄. Ce résultat donne une condition suffisante pour que la méthode de Newton converge
vers x̄. Mais la condition n’est pas nécessaire. Dans le cas où il n y a pas convergence,
on change la valeur de x0 .
? La méthode de Newton est une méthode par itération (ou méthode des approximations
successives). Ces méthodes se traduisent par répétition systématique d’une opération
déterminée, ce qui finit par donner de proche en proche la solution exacte.
? La formule de Newton est l’une des méthodes les plus utilisées parce qu’elle est relativement
simple et d’une convergence très suffisante dans la plupart des cas. Il faut toutefois que
x0 soit assez proche de la racine pour que la dérivée f 0 ne change pas de signe entre x0
et x. Si en effet f 0 s’annule sans que f s’annule aussi, la tangente à la courbe va couper
l’axe des x à l’infini. On risque alors d’osciller indéfiniment autour de la solution sans
jamais l’atteindre. Il est possible de modifier la formule de Newton et d’écrire

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 40

Fig. 5.1 – organigramme de la méthode de Newton-Raphson

f (xn )
xn+1 = xn − p
f 0 (xn )
où p est un nombre ou une fonction. On a aussi

g f (xn )
xn+1 = xn −
g f 0 (x
n ) + h f (xn )

c) Condition d’arrêt
La méthode de Newton converge vers la solution exacte x̄, mais on ne connaît pas explici-
tement à quel moment une certaine précision est atteinte. Pour cela, on utilise généralement
le résultat suivant :
Soit x̄ une racine de l’équation f (x) = 0 et soit xn une valeur approchée de x̄. Si sur [a, b]
contenant x̄ et xn on a
|f 0 (x)| ≥ m > 0 alors |xn − x̄| ≤ |f (x
m
n )|
.
Dans la pratique, on pourra prendre m comme m = min (|f 0 (a)| , |f 0 (b)|) et si est E la
précision exigée, s’arrêter lorsque |f (x
m
n )|
< E.
d) organigramme de la méthode (voire figure ci-dessus)
e) Remarques sur les formules du troisième et du premier ordre
En développant la série de Taylor jusqu’à l’ordre 2, on obtient

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 41

h2 00
f (x0 + h) = f (x0 ) + hf 0 (x0 ) +
f (x0 ) + O(h3 ).
2
En remplaçant h par h = x − x0 = − ff0 (d’après Newton), on obtient
µ ¶2
f (xn ) 1 f 00 (xn ) f (xn )
xn+1 = xn − 0 − .
f (xn ) 2 f 0 (xn ) f 0 (xn )
Soit l’équation f (x) = 0. De cette équation, on peut tirer x en fonction de x en écrivant
x = g(x), où g est défini à partir de f (x). Dans ce cas on peut écrire xn+1 = g(xn ). C’est la
formule d’itération du premier ordre.
e) Formule de Newton pour les racines complexes
Soit l’équation f (z) = 0 avec z = x + iy ;
on a
f (z0 )
f (z) ≈ f (z0 ) + (z − z0 )f 0 (z0 ) ⇒ z − z0 = −
f 0 (z0 )
Pour le calcul numérique, il faut séparer les parties réelles et imaginaires. On aura donc
df
f (z0 ) = G(x0 , y0 ) + iH(x0 , y0 ), = f 0 (z0 ) = G0 (x0 , y0 ) + iH 0 (x0 , y0 )
dz0
d’où
( GG0 +HH 0 ¯
¯
xn+1 ≈ xn − G02 +H 02 (xn,yn )
GH 0 −HG0 ¯
¯
yn+1 ≈ yn − G02 +H 02 (xn,yn )

5.2.2 Méthode de Dichotomie ou Méthode de Bissection

Supposons que f (a) < 0 < f (b) et f est strictement croissante sur [a, b], la solution x0
de l’équation f (x) = 0 à l’intervalle [a, b]. Le milieu de l’intervalle est c = (a + b)/2. Si f (a)
et f (c) ont des signes opposés, alors x0 ∈ [a, b]. Si par contre f (b) et f (c) sont de signes
contraires, alors x0 ∈ [b, c]. Dans l’un ou l’autre cas, la taille de l’intervalle est divisée par 2.
Une nouvelle bissection peut être faite sur le nouvel intervalle contenant la solution et ainsi
de suite. Ainsi, après n étapes, la taille de l’intervalle devient (a − b)/2n .

EXERCICE 2 : Ecrire l’organigramme de cette méthode.

5.3 CALCUL DE TOUTES LES RACINES D’UNE EQUATION POLYNOMIALE NON

LINEAIRE

5.3.1 Opération sur les polynômes

a) Schéma de calcul numérique d’un polynôme
Pour déterminer les valeurs numériques du polynome
Pn (x) = a0 xn + a1 xn−1 + ... + an−1 x + an ,
on peut calculer chaque terme et faire la somme. Mais, ce procédé n’est pas en général facile.
On peut utiliser le schéma dit de Horner. On calcule les nombres bj définis par :

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 42



 b0 = a0

 b1 = a1 + b0 x



 b2 = a2 + b1 x

 ..
.



 bp = ap + bp−1 x

 ..

 .


bn = an + bn−1 x
bn est alors la valeur numérique cherchée.

Algorithme :

Début
x←
b0 ← a0
Pour k allant de 1 à n faire
bk ← ak + bk−1 x
Fin pour
Ecrire bn = Pn (x)
Fin
La première méthode de calcul exige (2n − 1) multiplications alors que le schéma de
Horner n’en demande que n multiplications et le nombre d’addition est le même dans les
deux cas. De plus le schéma de Horner a l’avantage de n’introduire qu’un type d’opération :
bk = ak + bk−1 x, ce qui facilite la programmation.
Pour un polynôme à coefficients complexes,

Pn (z) = c0 z n + c1 z n−1 + · · · + cn−1 z + cn

où z = x + iy, et cp = ap + ibp ,
on a Pn (z) = Gn (x, y) + iH (x, y).
Le schéma de Horner peut être également être utilisé. On remplace tout simplement les
quantités bp par les quantités Pp = Gp + iHp avec

G0 = a0 , H0 = b0
Gp = ap + xGp−1 − yHp−1 , Hp = bp + yGp−1 + xHp−1
b) Division d’un polynôme par x − r
En divisant le polynôme
Pn (x) = a0 xn + a1 xn−1 + ... + an−1 x + an par C(x) = x − r on obtient un quotient Q(x)
et un reste R liés entre eux par l’identité P (x) = C(x)Q(x) + R.
Q(x) est un polynôme de degré n − 1 et R est une constante. Ils sont de la forme

Q(x) = b0 xn−1 + b1 xn−2 + ... + bn−2 x + bn−1 ; R = bn

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 43



 b0 = a0

 b1 = a1 + b0 r



 ..

 .
avec bp = ap + rbp−1

 ..

 .



 b = an−1 + rbn−2

 n−1
bn = an + rbn−1
Si bn = 0 alors le polynôme est divisible par x − r c’est-à-dire que r est une racine de
P (x).

Remarques :

? Ces formules sont semblables pour des grandeurs complexes. La seule différence apparaît
dans le calcul effectif qui exige un dédoublement de colonne comme ci-dessus.
? Quand les coefficients ap sont réels et que r = α + iβ on peut diviser P (x) par le facteur
quadratique (x − r)(x − r̄) où r̄ est le complexe conjugué de r.
c) Division d’un polynôme par un facteur quadratique
Soit le polynôme Pn (x) que l’on désire diviser par le facteur quadratique Q2 (x) = x2 −
Sx + P . On a Pn (x) = G2 (x)Qn−2 (x) + R1 (x). Qn−2 (x) et R1 (x) sont respectivement les
polynômes de degré n − 2 et 1 que l’on peut écrire sous la forme

Qn−2 (x) = b0 xn−2 + b1 xn−3 + ... + bn−3 x + bn−2 et R1 (x) = bn−1 (x − S) + bn

avec


 b0 = a0

 b1 = a1 + Sb0



 b2 = a2 + Sb1 − P b0



 ..

 .
bq = aq + Sbq−1 − P bq−2
 .


 .
.



 bn−2 = an−2 + Sbn−3 − P bn−4



 b = an−1 + Sbn−2 − P bn−3

 n−1
bn = an + Sbn−1 − P bn−2
Si S et P sont respectivement la somme et le produit de deux racines du polynôme Pn (x)
alors le reste R1 est nul, c’est-à-dire bn−1 = 0 et bn = 0.

5.3.2 Calcul d’une racine réelle, passage à l’équation de degré n − 1

L’idée la plus simple est de calculer la suite des valeurs bp qui conduisent à la valeur du
an
polynôme et comme on devrait avoir bn = 0, on obtiendrait x = − bn−1 .
On recommence le calcul des quantités bp avec cette nouvelle valeur de x jusqu’à obtenir
la précision voulue. A ce moment les coefficients b0 , b1 ,..., bn−1 donnent le polynôme restant
de degré n − 1. On peut alors continuer le processus de recherche des racines sur le polynôme
de degré n − 1

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 44

EXERCICE 3 : Etablir l’algorithme de ce schéma.

L’inconvénient de ce processus est qu’il n’est pas toujours convergent. Cette méthode
appelée méthode de Lin fut développé en 1941. Elle peut également s’appliquer sur le plan
complexe. Pour cela, on part d’une valeur arbitraire z0 = x0 +iy0 et on calcule les coefficients
de Horner complexes jusqu’à Gn−1 et Hn−1 . On obtient alors des nouvelles valeurs z1 =
x1 + iy1 en écrivant que Gn = 0 et Hn = 0. On a en effet
½ (
an + x1 Gn−1 − y1 Hn−1 = 0 x1 = a n G
Gn−1 +bn Hn−1
2 2
n−1 +Hn−1
⇒ an Hn−1 −bn Gn−1
bn + x1 Hn−1 + y1 Gn−1 = 0 y1 = G2 +H 2
n−1 n−1

On recommence le processus jusqu’à ce que

¯ (k+1) ¯
¯z − z (k) ¯ < ε.

5.3.3 Calcul de deux racines et passage à l’équation de degré n − 2

Dans le cas de la division d’un polynôme par le facteur quadratique G2 (x) = x2 − Sx + P ,
si S et P sont la somme et le produit de deux racines du polynôme Pn (x), alors on aura

bn = 0 et bn−1 = 0
an
⇒ P = bn−2 et S = P bn−3 −an−1
bn−2
.
Le processus peut alors continuer jusqu’à bn = bn−1 = 0 à la précision choisie. On obtient
finalement un nouveau polynôme de degré n − 2 pour lequel on peut reprendre le calcul pour
retrouver les deux racines suivantes et ainsi de suite, on parvient à résoudre l’équation.

5.3.4 Formule d’itération de Newton sur l’axe réel

Soit l’équation Pn (x) = 0, la dérivée de ce polynôme en x est Pn0 (x). Les valeurs numériques
Pn (x) et de Pn0 (x) se calculent par le schéma de Horner suivant :
 

 b 0 = a 0 
 c0 = b0

 b = a + xb 


 1

1 0  c1 = b1 + xc0


 
 .b2 = a2 + xb1
  c. 2 = b2 + xc1

.. et ..

 


 bp = ap + xbp−1 
 cp = bp + xcp−1

 . 
 ..

 .. 
 .

 

bn = an + xbn−1 cn = bn + xcn−1
on a Pn (x) = bn et Pn0 (x) = cn−1 . La formule d’itération de Newton devient
bn
xk+1 = f (xk ) où f (x) = x − cn−1 .
On peut alors poursuivre le processus jusqu’à obtenir |xk+1 − xk | < ε.
Après avoir trouvé une racine x∗ de l’équation, on divise le polynôme par x − x∗ et on
continue le processus jusqu’à l’obtention de toutes les racines.

5.3.5 Méthode de Bairstow

Soit à diviser le polynôme Pn (x) = a0 xn +a1 xn−1 +· · ·+an−1 x+an par G2 (x) = x2 −Sx+P .

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 45

On obtient Pn (x) = G2 (x)Qn−2 (x) + R1 (x) où

Qn−2 (x) = b0 xn−2 + b1 xn−3 + · · · + bn−3 x + bn−2 et R1 (x) = bn−1 (x − S) + bn ,



 b0 = a0

 b1 = a1 + Sb0



 b2 = a2 + Sb1 − P b0



 ..

 .
avec bq = aq + Sbq−1 − P bq−2

 ..

 .



 bn−2 = an−2 + Sbn−3 − P bn−4



 bn−1 = an−1 + Sbn−2 − P bn−3


bn = an + Sbn−1 − P bn−2
Si l’on a bn = 0 et bn−1 = 0, alors S et P sont la somme et le produit de deux racines du
polynôme et les grandeurs bq sont les coefficients du polynôme de degré n − 2. La méthode
de Bairstow consiste à partir des valeurs arbitraires de S et P et de calculer la suite des
coefficients bq .
Soient
F (S, P ) = bn−1 et G(S, P ) = bn − Sbn−1 ,
Si F et G ne sont pas nuls à la précision choisie, alors on améliore S0 et P0 (S0 et P0 étant
les valeurs initiales de S et de P ) par la formule d’itération de Newton à deux variables. Les
deux fonctions de S et de P sont F et G et les formules d’itération se mettent sous la forme
α β
Sk+1 = Sk + et Pk+1 = Pk +
∆ ∆
∂G ∂F ∂F ∂G ∂G ∂F ∂F ∂G
avec α = F −G , β=G −F et ∆ = − .
∂P ∂P ∂S ∂S ∂S ∂P ∂S ∂P
Pour obtenir les dérivées partielles par rapport à S, il suffit de dériver les relations bq par
rapport à S. Pour cela nous posons
∂bq ∂bq
= cq−1 et − bn−1 = cn−1
∂S ∂S
On obtient alors la suite suivante


 c0 = b0

 c1 = a1 + Sc0




 .c2 = a2 + Sc1 − P c0

..

 cq = aq + Scq−1 − P cq−2



 .

 ..


cn−2 = Scn−3 − P cn−4
∂bq
Posons ∂P
= −cq−2 ,

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 46

On constate que l’on obtient une suite identique à celle définie ci-dessus. Il vient finale-
ment : 
 α = bn cn−3 − bn−1 cn−2
β = bn cn−3 − bn−1 cn−1

∆ = cn−2 − cn−1 cn−3
Ces formules ont été établies par Bairstow en 1914 pour des besoins d’aéronautique. Avec
ces nouvelles valeurs de S et de P , on calcule à nouveau les coefficients bq et cq ainsi de suite
jusqu’à ce que bn−1 et bn soient nuls à la précision choisie. Une fois que l’on a terminé le
calcul de S et de P , on détermine les deux racines correspondantes par :
1h √ i 1h √ i
x= S ± S 2 − 4P ou x + iy = S ± i −S 2 + 4P .
2 2
Les coefficients restant bq sont ceux de l’équation de degré n − 2. On peut alors recom-
mencer le processus et calculer par paire toutes les racines. Toutefois, si n est impaire, il
restera à la fin une équation de la forme :
b0
b0 + b1 x = 0 ⇒ x = − .
b1

EXERCICE 4 :

1. Etablir l’algorithme et le programme en Pascal de la méthode de Newton sur l’axe réel

(qui permet de trouver toutes les racines réelles) d’une équation polynomiale.
2. Etablir l’algorithme de la méthode de Bairstow et traduire cet algorithme en Fortran.

5.3.6 Calcul simultané de toutes les racines

Pour calculer simultanément toutes les racines x1 , x2 ,..., xn d’un polynôme de degré n,
il faut disposer de n équations entre ces racines. Malheureusement, c’est souvent difficile
d’obtenir ce système d’équations entre les racines. Mais pour les équations de degré inférieur
à 5, on peut trouver des relations : relations de Viete
a) Cas où toutes les racines sont réelles
Si nous considérons l’équation du 3ieme degré x3 + a1 x2 + a2 x + a3 = 0, on a les relations
de Viete suivantes entre les racines :
 
 x1 + x2 + x3 = −a1  x1 = −a1 − x2 − x3
x1 x2 + x1 x3 + x2 x3 = a2 ⇔ x2 = (a2 − x1 x3 − x2 x3 )/x1
  x = − a3
x1 x2 x3 = −a3 3 x1 x2

On obtient ainsi un schéma itératif du 1er ordre qui permet par approximations successives
d’avoir les trois racines de l’équation.
 (k+1) (k) (k)

 x1 =³−a1 − x2 − x3 ´.
 (k+1) (k) (k) (k) (k) (k)
x2 = a2 − x1 x3 − x2 x3 x1

 .³ ´
 x(k+1) = −a3 x(k) x(k)
3 1 2

b) Cas où il y a des racines complexes

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 47

Si les coefficients de l’équation sont réels, les racines complexes se présentent par paires
conjuguées l’une de l’autre. En associant les racines par couple et en faisant apparaître les
sommes et les produits, on opère uniquement sur les quantités réelles.
Pour une équation du 3ieme degré x3 + a1 x2 + a2 x + a3 = 0, si on pose x1 + x2 = S,
x1 x2 = P et x3 = r, les relations de Viete s’écrivent :
 
 S + r + a1 = 0  S = −r − a1
P + rS − a2 = 0 ⇔ P = a2 − rS
 rP + a = 0  r = −a /P
3 3

On peut alors partir des valeurs quelconques de S, P et r pour faire le processus itératif
et déterminer ces grandeurs à la précision choisie.
Pour une équation du 4ieme ordre du type x4 +a1 x3 +a2 x2 +a3 x+a4 = 0, on a les relations
suivantes :
 

 S1 + S 2 + a1 = 0 
 S = x1 + x2
  1
S1 S2 + P1 + P2 − a2 = 0 P1 = x1 x2
avec

 S1 P2 + S2 P1 + a3 = 0 
 S2 = x3 + x4
 P P −a =0  P =x x
1 2 4 2 3 4
er
On peut alors tirer le schéma itératif du 1 ordre suivant :

 S1 = −a1 − S2


P1 = a2 − P2 − S1 S2

 S2 = −(a3 + S1 P2 )/P1
 P = a /P
2 4 1

Et pour une équation du 5ieme degré on a :



 S1 = − (a1 + S2 + r)


 P1 = a2 − P2 − S1 S2 − r (S1 + S2 )
S2 = −[a3 + S1 P2 + r (P1 + P2 + S1 S2 )]/P1

 P2 = [a4 − r (S1 P2 + S2 P1 )]/P1


 r = −a /P P
5 1 2

où r = x5 .

Remarques

? Relations de Viete
Soit le polynôme de degré n en z défini par : Pn (z) = a0 z n + a1 z n−1 + · · · + an−1 z + an
où les aj et z sont des réels ou complexes.
Soit zj les racines de ce polynôme, on a alors Pn (z) = a0 (z − z1 )(z − z2 ) · · · (z − zn ).
En développant cette expression et après identification des coefficients on obtient les
relations suivantes :

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 48

 P

 a1 = −a0 zi

 Pi

 a = a zi zj

 2 0
 i≺j P
a3 = −a0 zi zj zk

 i≺j≺k



 ..

 .

an = (−1)n a0 z1 z2 ...zn
? Localisation des racines dans le plan
– Règle des signes de Descartes
Le nombre de racines réelles positives ne peut pas dépasser le nombre de changements
de signe que l’on observe en parcourant la suite des coefficients du polynôme. Il ne peut
différer de ce nombre que par un nombre pair. On a une règle analogue pour les racines
négatives. Il suffit de changer x en −x.
Exemple : L’équation x4 −8x3 +24x2 −32x+15 = 0 n’a pas de racines réelles négatives.
De même, l’équation x4 − 2x3 − 3x2 + 8x − 4 = 0 n’a pas plus de trois racines réelles
positives.
– Règle de Gua
Si le carré d’un coefficient intermédiaire est inférieur ou égal au produit des coefficients
voisins, alors il y a des racines imaginaires.
Exemple : l’équation x3 + 7x2 + x + 7 = 0 possède une racine imaginaire 1 × 1 < 7 × 7
– Règle des lacunes
S’il manque un terme entre deux termes de même signe ou plusieurs termes entre deux
termes de signe quelconque, alors il y a des racines imaginaires.
Exemple : l’équation x3 + 6x + 20 = 0 a le terme en x2 absent et les coefficients des
termes x3 et x sont positifs. On peut donc affirmer qu’il y a des racines complexes.
– Règle de Sturn
Soit a un nombre réel positif, si l’équation (x − a)Pn (x) = 0 présente (2k + 1) variations
de signe de plus que l’équation Pn (x) = 0 alors il y a au moins 2k racines complexes
Exemple : Soit l’équation x3 + 7x2 + x + 7 = 0 et l’équation (x − 1)Pn (x) = x4 + 6x3 −
6x2 + 6x + 7 = 0 présente trois changements de signe, il y a donc 2 racines complexes
pour l’équation Pn (x) = 0.

5.4 SYSTEMES D’EQUATIONS NON LINEAIRES

Soit un système d’équation non linéaires

½
i = 1, ..., n
fj (xi ) = 0 avec
j = 1, ..., n
On peut aussi élaborer des processus itératifs du 1er, 2ième et 3ième ordre pour ce type
de système d’équations. En particulier, pour le système de deux équations non linéaires
½
f1 (x, y) = 0
f2 (x, y) = 0
on peut trouver x à partir de f1 (x, y) et y à partir de f1 (x, y) et utiliser le processus
itératif du 1er ordre suivant :

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

SYSTEMES D’EQUATIONS ALGEBRIQUES LINEAIRES ET NON LINEAIRES 49

½
xk+1 = f1 (xk , yk )
yk+1 = f2 (xk , yk )
On peut également généraliser les formules itératives de Jacobi et de Gauss-Seidel. De
même, à partir de la méthode de Newton, le système de deux équations non linéaires à 2
inconnues donne les formules itératives suivantes :
 ¯
 xk+1 = xk + 1 (f2 f1,y − f1 f2,y ) ¯¯ ½
∆ ¯ (xk ,yk ) ∆ = f1,x f2,y − f1,y f2,x
¯ avec
 yk+1 = yk + (f1 f2,x − f2 f1,x ) ¯
1 fi,x = ∂f
∂x
i
∆ (xk ,yk )

Dans le cas d’un système de degré supérieur à 2, on élabore le schéma de Newton de

(0)
la manière suivante : supposons que la solution d’ordre zéro soit le vecteur xi et que la
(1)
solution d’ordre 1 est xi , On a
³ ´ ³ ´ ³ ´
(1) (0) (1) (0) ∂fj
f j xi = f j xi + xi − xi | (0)
∂xi xi
³ ´
(1) (1)
Or xi étant solution de l’équation, on a fj xi = 0 ; d’où le système matriciel
³ ´ ∂f ³ ´
(1) (0) j (0)
xi − x i |x(0) = −fj xi
∂xi i

(1)
En évaluant ce système matriciel, on obtient le vecteur xi . On reprend le processus et
ainsi de suite. A chaque étape on résout le système matriciel :
³ ´ ∂f ³ ´
(k+1) (k) j (k)
xi − xi |x(k) = −fj xi
∂xi i

(k)
Dès que les quantités xi sont trouvées à la précision choisie, on a alors l’ensemble des
racines du système algébrique non linéaire.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 Chapitre VI

EQUATIONS AUX DERIVEES PARTIELLES ET

PROBLEMES DE STABILITE

6.1 NOTIONS GENERALES

6.1.1 Définitions
Considérons une équation aux dérivées partielles d’ordre 2 de la forme :

∂ 2U ∂ 2U ∂ 2U ∂U ∂U
A 2
+ 2B + C 2
+ F ( x, y, U, , )=0
∂x ∂x∂y ∂y ∂x ∂y
où A, B et C sont des fonctions de x et y et sont de classe C 2 .
¨ Si B 2 − AC = 0 alors l’équation est de type parabolique ;
¨ Si B 2 − AC < 0, l’équation est de type elliptique ;
¨ Si B 2 − AC > 0, l’équation est de type hyperbolique.
Pour résoudre une équation aux dérivées partielles, il faut connaître les conditions aux
limites ou aux frontières ainsi que les conditions initiales. La nature de ces conditions permet
de définir les types de problèmes suivants :
- Problème de DIRICHLET pur : les valeurs de U sur les frontières du domaine sont connues.
- Problème de NEUMANN pur : les valeurs des dérivées de U sur les frontières sont connues.
- Problème mixte DIRICHLET-NEUMANN : les valeurs de U sont connues sur une partie
des frontières et celle des dérivées de U sur la partie restante.
- Problème de CAUCHY : les valeurs de U et de ses dérivées sont connues sur les frontières
du domaine.

6.1.2 Quelques équations aux dérivées partielles en Physique et en Mécanique

a) Equation d’ondes
2
Elle a la forme : ∂∂tU2 = a2 ∆U + f ( ∂U
∂t
, x, y, z, t )
où a est la vitesse de l’onde et f une perturbation quelconque.
En dimension 1, cette équation décrit :
(ı)− Les vibrations transversales d’une corde. L’équation générale dans ce cas est :

∂ 2U ∂U ∂2U
ρ(x) + k = T + f (x, t)
∂t2 ∂t ∂x2

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABILITE 51

où ρ est la masse volumique, k le coefficient

p de dissipation, T la tension dans la
corde . Ici la vitesse des ondes est a = T /ρ.
p
(ıı)− Les vibrations longitudinales dans une barre avec a = E/ρ, E étant le module
d’élasticité.
p E
(ııı)− Les vibrations de torsion d’une barre avec a = µ/ρ, µ = 2(1+σ) où σ est le
coefficient de Poisson.
(ıv)− La propagation d’ondes électromagnétiques.
(v)− La propagation d’ondes électriques (équation des télégraphistes).
En dimension 2, elle peut décrire les vibrations d’une membrane plane. En dimension
3, elle décrit les ondes sonores en hydrodynamique, les ondes électromagnétiques (cas des
cavités résonantes et du rayonnement des antennes).
b) Equation de Poisson et de Laplace
L’équation de Poisson est une équation de la forme ∆U = f (x, y, z). Lorsque le second
membre de cette équation devient nul, elle se ramème à l’équation de Laplace : ∆U = 0.
Ce sont les prototypes des équations elliptiques pouvant décrire les problèmes station-
naires suivants :
(ı)− La répartition des charges électriques sur la surface d’un conducteur ;
(ıı)− L’équilibre thermique dans un corps homogène ;
(ııı)− Le potentiel des vitesses dans un fluide non visqueux, incompressible et en écou-
lement irrotationnel.
c) Equation de la chaleur
−−→
Sa forme générale est cρ ∂U∂t
= div ( k gradU ) + F (x, y, z, t ) où U est la température
en un point P (x, y, z) au temps t, ρ la masse volumique, c la chaleur spécifique, k la
conductivité interne et F l’apport des sources ou des puits de chaleur. Si c, ρ et k sont
constantes, cette équation se réduit à ∂U ∂t
= a2 ∆U + f ( x, y, z, t ).
Dans le cas où ∂U∂t
= 0, on est en présence d’un problème de transfert de chaleur en
régime permanent (conduction et convection).
Dans le cas général où ∂U ∂t
6= 0, on est en présence d’un problème de diffusion de la
chaleur.

6.2 DISCRETISATION

6.2.1 L’équation de Poisson dans le plan

Soit l’équation

∂ 2U ∂ 2U
+ = f (x, y ) (P )
∂x2 ∂y 2
définie sur un domaine D rectangulaire : D = {(x, y) ∈ R2 \ a < x < b et c < y < d}.
Soit (S) la frontière de D, la condition suivante : U (x, y) = g(x, y) est imposée sur (S). Nous
considérons que les fonctions f et g sont continues sur le domaine D. Soient h et k les pas
de calcul suivant x et y respectivement. Soient n et m les nombres entiers tels que

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABILITE 52

h = (b − a)/n et (d − c)/m. Un point du domaine D est défini par les cordonnées sui-
vantes :
½
xi = a + ih avec i = 1, 2..., n
yj = c + jk avec j = 1, 2, ..., m
D’après les définitions des dérivées, on a :

∂ 2U Ui+1, j − 2Ui, j + Ui−1, j h ∂ 4U

= − (ξij yj ) (P 1)
∂x2 h2 12 ∂x4
où Uij = U (xi , yj ) ; ξij ∈ ]xi−1, xi+1 [
2
En écrivant une équation similaire pour ∂∂yU2 , l’équation (P ) devient avec une erreur de
troncature d’ordre o (h2 + k 2 ), le système suivant :

" µ ¶ # µ ¶2
2
h h
2 + 1 Ui, j − (Ui+1, j + Ui−1, j ) − (Ui, j+1 + Ui, j−1 ) = −h2 fi,j (P 2)
k k

où fi,j = f (xi , yj ) ; i = 1, 2, . . . , n − 1 ; j = 1, 2, · · · , m − 1.
La discrétisation des conditions aux frontières est la suivante :

 U0, j = g(x0 , yj )
j = 0, 1, · · · , m

Un, j = g(xn , yj ) (P 3)

 Ui, 0 = g(xi , y0 )
i = 0, 1, · · · , n

Ui, m = g(xi , ym )
En associant les équations (P 2) et (P 3), on obtient le système d’équations algébrique
linéaire en Uij que l’on peut résoudre par la méthode d’élimination de Gauss.

EXERCICE 1 : Ecrire le schéma d’itération de Gauss-Seidel et de Jacobi pour résoudre

l’équation de poisson.

Remarques :

R1 : Lorsque les conditions aux limites sont de type Neumann, le schéma de discrétisation
doit se faire avec beaucoup de précautions. Généralement, il est conseillé de combiner
cette condition avec l’équation aux dérivées partielles de
¯ départ.
¯
Par exemple, si la condition de Neumann s’écrit, ∂U ∂y ¯
= g(x, y), alors on utilise le
S
schéma suivant. On développe U (xi , y1 ) autour du point U (xi , y0 ) ; c’est à dire que l’on
écrit :

∂U (xi , y0 ) k 2 ∂ 2 U (xi , y0 )
U (xi , y1 ) = U (xi , y0 ) + + 2
+ O(k 3 )
∂y 2 ∂y
2
La quantité ∂∂yU2 en (xi , y0 ) peut être remplacée par son équivalent à partir de l’équation
aux dérivées partielles de départ.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABILITE 53

∂2U 2
Par exemple, dans le cas de l’équation de Laplace, on a : ∂y 2
= − ∂∂xU2 . Il vient alors :
2 2
U (xi , y1 ) = U (xi , y0 ) + k ∂U (x∂yi ,y0 ) − k2 ∂ U∂x
(xi ,y0 )
2
∂U (xi ,y0 ) k2
= U (xi , y0 ) + k ∂y − 2h2 (Ui+1,0 − 2Ui,0 + Ui−1,0 )
Finalement la condition de Neumann donne :
" µ ¶2 # µ ¶2
k 1 k
Ui,1 − 1 + Ui,0 + (Ui+1,0 + Ui−1,0 ) = kg(xi , y0 ) (P 4)
h 2 h

On peut alors associer à (P 4) l’équation (P 2) pour avoir le système discret global à

résoudre.
R2 : Théorème du maximum et du minimum discret. Dans le cas de l’équation de Laplace,
on a

k2 h2
Ui,j = (U i+1,j + U i−1,j ) + (Ui,j+1 + Ui,j−1 )
2(k 2 + h2 ) 2(k 2 + h2 )
Donc Ui,j est le barycentre des nombres Ui+1,j , Ui−1,j , Ui,j+1 et Ui,j−1 affectés des poids
2 2
respectifs : α = 2(k2k+h2 ) , α, β = 2(k2h+h2 ) , et β.
Puisque ces poids sont tous positifs, Ui,j est donc compris entre la plus grande et la plus
petite des quantités Ui+1,j , Ui−1,j , Ui,j+1 et Ui,j−1 . D’où le théorème suivant :

Théorème : Il n’y a pas de maximum et de minimum strict de U à l’intérieur du

domaine D. Le maximum et le minimum sont atteints sur les frontières de D.
R3 : D’autres schémas de différences finies peuvent être établis à partir des formules de
2
discrétisation des dérivées. Par exemple, on peut définir ∂∂xU2 par des schémas présentant
des erreurs d’ordre O(h4 + k 4 ).

6.2.2 Discrétisation de l’équation de la chaleur

Nous nous limiterons au cas unidimensionnel. On a l’équation :
2
∂U 2∂ U
=a (C1)
∂t ∂x2
Dans le domaine D = T × L avec T = R+ et L = ]0, l[.
Les conditions aux limites sont de la forme U (0, t) = 0, U (l, t) = 0 et la condition initiale
est U (x, 0) = f (x).
Soit h le pas spacial et k le pas temporel et soit m le nombre entier tel que mh = l

(i) Méthode des différences progressives ou schéma explicite

U −U 2 U −2U +U
on a ∂U
∂t
= i, j+1k i, j et ∂∂xU2 = i+1, j hi,2 j i−1, j .
Dans le domaine D, on obtient le système discret :
³ 2
´ 2
Ui,j+1 = 1 − 2ah2k Ui,j + ah2k (Ui+1,j + Ui−1,j )
i = 1, 2, . . . , m − 1; j = 1, 2, . . . , ∞ (C2)

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABILITE 54

Le schéma (C2) présente une erreur de troncature d’ordre o(k + h2 ). C’est un schéma
explicite car pour calculer Ui,j+1 , il suffit de remplacer les quantités figurant au second
membre par leurs valeurs connues et stockées en mémoire. Les conditions initiales et
aux limites se discrétisent de la manière suivante :

Ui,0 = f (xi ), i = 0, 1, . . . , m
U0,j = 0 et Um,j = 0, j = 0, 1, 2, . . . , ∞ (C3)
La première condition (C3) peut être substituée dans (C3) pour déterminer Ui,1 pour i
allant de 1 jusqu’à m − 1. La condition U0,j = 0, et Um,j = 0 ⇒ U0,1 = Um,1 = 0. On a
toutes les valeurs de Ui,1 . En procédant de la même manière, on obtient les valeurs de
Ui,2 , Ui,3 ainsi de suite jusqu’à Ui,m−1 . L’inconvénient de cette méthode est qu’elle est
conditionnellement stable. Pour le démontrer il faut déterminer la condition de stabilité.

Condition de stabilité :

Considérons le schéma aux différences finies progressives de l’équation de la chaleur,

c’est-à-dire :

Ui, j+1 − Ui, j (Ui+1, j − 2Ui, j + Ui−1, j )

= a2 ,
k h2
Cette méthode consiste à introduire une perturbation εi,j de la solution numérique Ui,j .
0 0
La solution numérique devient alors : Ui,j = Ui,j + εi,j . En écrivant que Ui,j vérifie le
schéma discret, on obtient que la perturbation εi,j vérifie aussi l’équation discrète :

εi,j+1 − εi,j εi+1,j − 2εi,j + εi−1,j

= a2
k h2
Le schéma discret est stable si εi,j reste borné ou tend vers 0 au fur et à mesure
que le temps avance. Pour avoir la condition de stabilité, nous supposons que εi, j =
sin(iph) e−rjk . En introduisant cette expression dans l’équation vérifiée par εi,j , on ob-
2
tient : erk = 1 − 4ka
h2
sin2 ph
2
; or εi, j+1 = εi, j e−rk donc e−rk est le facteur d’amplification.
Le schéma est donc stable si r > 0, i.e. que le facteur d’amplification est inférieur à 1.
dans le cas particulier analysé ici, cette condition implique que :
¯ ¯ ¯ ¯
¯ 4ka 2
ph ¯ ¯ 4ka 2¯ 2
¯1 − sin 2 ¯ < 1 ⇒ ¯1 − ¯ < 1 ⇒ ka < 1
¯ h2 2¯ ¯ h2 ¯ h2 2

(ii) Différences finies régressives ou méthode implicite :

Afin d’éviter le problème de stabilité, on peut utiliser les différences finies régressives
pour la dérivée temporelle, i.e.

∂U Ui,j − Ui,j−1
= (C4)
∂t k
L’équation d’onde (C1) donne alors le système discret suivant :
(1 + λ)Ui,j − Ui+1,j − Ui−1,j = Ui,j−1
1 ≤ i ≤ n − 1; 1≤j≤∞ (C5)

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABILITE 55

où λ = a2 k/h2 . Compte tenu des conditions initiales, et des conditions aux limites, on
obtient un système matriciel tri diagonal que l’on peut résoudre par les méthodes itéra-
tives ou la méthode d’élimination de Gauss. Pour déterminer la condition de stabilité,
on procède de la même façon que précédemment. On trouve pour ce schéma implicite,
que le facteur d’amplification est
Áµ ¶
rk 4ka2 2 ph
e =1 1 + 2 sin
h 2
Il est toujours inférieur à 1, donc le schéma est inconditionnellement stable.

(iii) Autres méthodes :

? Méthode de Richardson
Les méthodes ci-dessus ont une erreur de troncature d’ordre o(k) en k. Si on veut
obtenir o(k 2 ), on utilise la formule suivante :

∂U Ui, j+1 − Ui, j−1

= (C6)
∂t 2k
Ce schéma n’est pas très utilisé parce qu’il est également instable.

? Méthode de Cranck-Nicolson
Cette méthode fait la somme des différences finies progressives et regressives. Après
avoir écrit la formule des différences regressives à l’instant j + 1, on obtient :

Ui,j+1 − Ui,j a2
− 2 (Ui+1,j − 2Ui,j + Ui+1,j+1 − 2Ui,j+1 + Ui−1,j+1 ) = 0
k 2h
C’est une méthode plus précise qui présente une erreur d’ordre o(h2 + k 2 ).

? Méthode de Dufort-Frankel
Pour résoudre le problème d’instabilité de la méthode de Richardson, la quantité
Ui,j est remplacée par (Ui,j+1 + Ui,j−1 )/2 et on obtient

Ui,j+1 − Ui,j−1
− a2 (Ui+1,j − Ui,j+1 − Ui,j−1 + Ui−1,j ) /h2 = 0
2k
Ce schéma est inconditionnellement stable.

EXERCICE 2 : Etablir l’algorithme pour le schéma explicite de même que les algorithmes
pour les méthodes des différences régressives et de Cranck-Nicolson.

6.2.3 Discrétisation de l’équation des ondes

Nous considérons le problème suivant :

∂2U 2
2∂ U
− a =0 (O1)
∂t2 ∂x2

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABILITE 56

dans le domaine D = T × L et soumis aux conditions suivantes :

∂U (x, 0)
U (0, t) = U (l, t) = 0; U (x, 0) = f (x); = g(x)
∂t
1◦ - Discrétisation directe
La discrétisation de (O1) par les différences finies donne :

Ui,j+1 = 2(1 − λ2 )Ui,j + λ2 (Ui+1,j + Ui−1,j ) − Ui,j−1 (O2)

l
où λ = ak/h ; i = 1, 2, . . . , m − 1 ; j = 1, 2, . . . , ∞ et m = h
.
Les conditions aux limites donnent :

U0,j = Um,j = 0 (O3)

et la première condition initiale donne

Ui,0 = f (xi ) (O4)

La deuxième condition initiale permet d’écrire :
Ui,1 = Ui,0 + kg(xi ) (O5).

Mais (O5) présente une erreur de troncature d’ordre o(k). Pour avoir une erreur d’ordre
o(k 2 ) comme l’équation (O2), on procède de la manière suivante :

∂U (xi , 0) U (xi , t1 ) − U (xi , 0) k ∂ 2 U (xi , 0)

= − 2
+ o(k 2 )
∂t k 2 ∂t
or
∂ 2 U (xi , 0) 2
2 ∂ U (xi , 0)
2
2∂ f fi+1 − 2fi + fi−1
2
= a 2
= a 2
= a2
∂t ∂x ∂x h2
En utilisant la condition initiale, on obtient finalement :

Ui,1 = (1 − λ2 )fi + λ2 (fi+1 + fi−1 )/2 + kgi (O6)

avec i = 1, 2, . . . , m − 1 ;
Il s’agit donc de résoudre le système discret formé des équations (O2) et (O6). La re-
cherche des valeurs Ui,j+1 nécessite la connaissance des valeurs Ui,j et Ui,j−1 . Au début, il
faut connaître Ui,0 et Ui,1 . Les valeurs de Ui,0 découlent de (O4) et celles de Ui,1 découlent
de (O6). On peut ainsi lancer le processus itératif pour calculer toutes les valeurs de Ui,j . Il
faut noter que ce schéma n’est stable que si λ < 1 ⇒ ak/h < 1.
2◦ - Décomposition en un système d’équations du premier ordre
Dans certains cas, il est avantageux de décomposer l’équation d’onde en deux équations
du premier ordre avant de discrétiser. Par exemple en posant V = ∂U
∂x
et W = a1 ∂U
∂t
, l’équation
d’onde se transforme en un système d’équations de la forme
∂V ∂W ∂W ∂V
=a ; =a
∂t ∂x ∂t ∂x

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABILITE 57

a)- Discrétisation explicite

On peut utiliser la discrétisation par les différences finies progressives dans le temps et
par les différences centrées dans l’espace. On obtient alors
½
Vi,j+1 = Vi,j + ak (Wi+1,j − Wi−1,j )/2h
Wi,j+1 = Wi,j + ak (Vi+1,j − Vi−1,j )/2h
Pour établir la condition de stabilité, posons
½
Vi,j = Aj eIpih
(I 2 = −1)
Wi,j = Bj eIpih
En substituant dans le système discret, on obtient
 h i
 Vi,j+1 = Aj + Bj Iak sin ph eIpih
jh
h i
 Wi,j+1 = Bj + Aj Iak sin ph eIpih
jh

or Vi,j+1 = Aj+1 eIpih et Wi,j+1 = Bj+1 eIpih

En comparant ces expressions, on obtient
µ ¶ µ ¶µ ¶
Ai,j+1 1 Ic Ai,j
=
Bi,j+1 Ic 1 Bi,j
où c = ak
h
sin ph
La stabilité du schéma exige que les valeurs propres de la µmatrice d’amplification
¶ A
1 Ic
aient des modules inférieurs à 1. Dans le cas présent, A = et les modules
Ic 1
des valeurs propres de A sont tous supérieurs à 1, donc le schéma est inconditionnelle-
ment instable.

b)- Amélioration du schéma explicite

A partir du schéma ci-dessus, on peut trouver tous les Vi,j+1 de la première équation
et les utiliser ensuite dans la deuxième équation. On obtient alors le schéma discret
suivant :
½
Vi,j+1 = Vi,j + ak(Wi+1,j − Wi−1,j )/2h
Wi,j+1 = Wi,j + ak (Vi+1,j+1 − Vi−1,j+1 )/2h
µ ¶
1 Ic
Pour un tel schéma, la matrice d’amplification est A = .
Ic 1 − c2
On trouve que le schéma est stable ssi ak h
< 2.

c)- Schéma implicite

On discrétise les dérivées spatiales à l’instant j + 1 et les dérivées temporelles par les
différences progressives. On obtient alors
½
Vi,j+1 = Vi,j + ak (Wi+1,j+1 − Wi−1,j+1 )/2h
Wi,j+1 = Wi,j + ak (Vi+1,j+1 − Vi−1,j+1 )/2h

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EQUATIONS AUX DERIVEES PARTIELLES ET PROBLEMES DE STABILITE 58

Les valeurs propres de la matrice d’amplification sont dans ce cas 1/(1+Ic) et 1/(1−Ic).
Elles ont toutes un module <1. Donc ce schéma implicite est inconditionnellement stable.

d)- Schéma de Lax

Si dans le schéma explicite on remplace Vi,j et Wi,j par leur valeurs moyennes respectives,
entre, i − 1 et i + 1, alors on obtient le schéma de Lax qui s’écrit
½
Vi,j+1 = (Vi+1,j + Vi−1,j )/2 + ak (Wi+1,j − Wi−1,j )/2h
Wi,j+1 = (Wi+1,j + Wi−1,j )/2 + ak (Vi+1,j+1 − Vi−1,j+1 )/2h
µ ¶
cos ph Ic
Dans ce cas la matrice d’amplification est A =
Ic cos ph
et la condition de stabilité est ak/h < 1.

e)- Schéma de Lax-Wendroff

Pour ce schéma, on développe Vi,j+1 et Wi,j+1 , jusqu’à l’ordre 2. Par exemple

∂Vi,j k 2 ∂ 2 Vi,j ∂Wi,j k 2 a2 ∂ 2 Vi,j

Vi,j+1 = Vi,j + k
+ = V i,j + ka +
∂t 2 ∂t2 ∂x 2 ∂x2
En faisant la même chose sur Wi,j+1 , on obtient le schéma discret de Lax-Wendroff
suivant :

½
Vi,j+1 = Vi,j + ak (Wi+1,j+1 − Wi−1,j )/2h + a2 k 2 (Vi+1,j − 2Vi,j + Vi−1,j )/2h
Wi,j+1 = Wi,j + ak (Vi+1,j+1 − Vi−1,j+1 )/2h + a2 k 2 (Wi+1,j − 2Wi,j + Wi−1,j )/2h

La matrice d’amplification est

µ 2 ¶
(h + a2 k 2 )(cos ph − 1)/h2 (Iak sin ph)/h
A=
(Iak sin ph)/h (h2 + a2 k 2 )(cos ph − 1)/h2

Condition de stabilité : akh

< 1.
Il faut noter cependant que lorsqu’une équation hyperbolique présente des conditions
initiales discontinues, la méthode des différences finies n’est plus adaptée ; on doit faire
appel à d’autres méthodes telle que la méthode des caractéristiques.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 Chapitre VII

RESOLUTION NUMERIQUE DES EQUATIONS DE

NAVIER STOKES ET DE BURGERS

7.1 EQUATION DE NAVIER STOKES

7.1.1 Généralités
D’après le principe fondamental de la dynamique, on a pour un écoulement l’équation

ργi = fi + σi,j
où :
fi = résultante des forces autre que les contraintes (Exemplr : forces de gravité)
σi,j = tenseur de contraintes
γi = accélération
ρ = masse volumique du fluide
dui ∂ui ∂~u
γi = = + ui,j uj ⇒ ~γ = = grad (~u.~u)
dt ∂t ∂t
Finalement,
∂~u 1 ¡ ¢
+ rot~u ∧ ~u + grad u2
~γ =
∂t 2
Dans le cas d’un fluide incompressible, on a

σij = −P δij + 2µεij

où
εij = tenseur de déformations ;
P = Pression.

σij,j = (−P δij ),j + µ(ui,j + uj,i ),j

= −P,j δij + µui,jj
= −P,i + µ∆ui
car div (~v ) = 0 ⇒ (uj,j ),i = 0.
Ainsi, pour un fluide incompressible, l’équation de Navier Stokes se met sous la forme :
· ¸
∂~u 1 ¡ 2¢
ρ + rot (~u) ∧ ~u + grad u = f~ − grad (P ) + µ∆~u (1.a)
∂t 2
div (~v ) = 0 (1.b)

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DE NAVIER STOKES ET DE BURGERS 60

a) Equation du tourbillon et de la fonction de courant ψ

~ = 1 rot (~u).
Définissons le tourbillon par : Ω 2
En appliquant l’opérateur 12 rot à l’équation (1a), on obtient l’équation suivante :
~
∂Ω ³ ´
~ ∧ ~u = υ∆Ω
+ rot Ω ~ (2)
∂t
où υ = µ/ρ.
Dans le cas d’un³ écoulement
´ plan, le vecteur tourbillon n’a qu’une seule composante
1 ∂u2 ∂u1
définie par : Ω = 2 ∂x − ∂y . L’équation du tourbillon prend alors la forme
∂Ω ∂(u1 Ω) ∂(u2 Ω)
+ + = υ∆Ω
∂t ∂x ∂y
Or le fluide est incompressible donc
∂u2 ∂u1
div (~v ) = 0 ⇒ + =0
∂x ∂y
on obtient alors
∂Ω ∂Ω ∂Ω
+ u1 + u2 = υ∆Ω (3)
∂t ∂x ∂y
Soit ψ(x, y, t) la fonction de courant définie par :
∂ψ ∂ψ
u1 = et u2 = −
∂y ∂x
On a alors

∆ψ = −Ω (4)
Les Equations (3) et (4) portant sur Ω et ψ sont couplées. L’équation (4) est linéaire,
tandis que (3) est non linéaire à cause des termes croisés u1 ∂Ω
∂x
et u2 ∂Ω
∂y
. Connaissant ψ, on
peut déduire u1 et u2 .
b) Equation de la pression
La connaissance du vecteur vitesse ~u peut permettre de déterminer la pression. A partir
de l’équation (1.a), On obtient
" µ 2 ¶2 #
∂2ψ ∂2ψ ∂ ψ
∆P = 2ρ 2
+ 2 −
∂x ∂y ∂x∂y
Ainsi connaissant ψ on peut aussi trouver la pression en intégrant l’équation de Poisson
ci-dessus.
c) Conditions initiales et aux limites
Pour les conditions initiales, il s’agit de donner par exemple dans le cas d’un problème plan
les valeurs de u1 et u2 ou de ψ à l’instant initial. Pour ce qui est des conditions aux limites,
désignons par (S) la frontière du domaine D dans lequel les équations de Navier Stokes sont
vérifiées. La vitesse sur la frontière (S) est définie par : ~u = ~us (t) avec la condition de flux

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DE NAVIER STOKES ET DE BURGERS 61

Z Z
~us (t).~nds = div~us dD = 0
s D

7.1.2 Discrétisation de l’équation d’advection-diffusion

Il s’agit de l’équation
∂Ω ∂Ω ∂Ω
+ u1 + u2 = υ∆Ω
∂t ∂x ∂y
où u1 (x, y, t) et u2 (x, y, t) sont supposées connues. Cette équation décrit par exemple la
distribution des particules où des gaz émis par une source plane en présence du vent. Dans
ce cas Ω est la concentration des particules, u1 et u2 sont les composantes de la vitesse du
vent. C’est un problème intéressant car il permet d’expliquer et de mesurer la pollution de
l’air par les déchets industriels.
– Si υ = 0 l’équation se réduit à l’équation de transport.
– Si u1 = u2 = 0 on a l’équation de la chaleur
a) Schéma de discrétisation centré en espace
Elle se fait de la manière suivante
(n+1) (n) (n) (n)
∂Ω Ωij − Ωij ∂Ω Ωi+1,j − Ωi−1,j
= et =
∂t k ∂x 2h
En discrétisant les autres termes, l’équation d’advection-diffusion se met sous la forme
discrète suivante :
³ ´ ³ ´
(n+1) (n) ku1 (n) (n) ku2 (n) (n)
Ωi,j = Ωi,j − 2h Ωi+1,j − Ωi−1,j − 2h Ωi,j+1 − Ωi,j−1
h i
(n) (n) (n) (n) (n)
+ υk
h2 Ω i+1,j + Ω i−1 + Ω i,j+1 + Ω i,j−1 − 4Ω 2
i,j + o(k + h )

υk
On peut établir que la condition de stabilité d’un tel schéma est h2
≤ 14 .

EXERCICE 1 : Etablir cette condition de stabilité.

b) Schéma explicité décentré en espace
Ce schéma s’écrit de la manière suivante :
(n+1) (n) (n) (n)
Ωi,j = Ωi,j h− kuh1 Pi,j − kuh2 Qi,j i
(n) (n) (n) (n) (n)
+ υk
h2
Ωi+1,j + Ω i−1 + Ω i,j+1 + Ω i,j−1 − 4Ω i,j

avec
( (
(n) (n) (n) (n)
(n) Ωi,j − Ωi−1,j si u1 > 0 (n) Ωi,j − Ωi−1,j si u2 > 0
Pi,j = (n) (n) Qi,j = (n) (n)
Ωi+1,j − Ωi,j si u1 < 0 Ωi,j − Ωi−1,j si u2 < 0
Ce schéma présente des conditions de stabilité similaires au schéma précédent.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DE NAVIER STOKES ET DE BURGERS 62

7.2 DISCRETISATION DE L’EQUATION DE NAVIER-STOCKES DANS LE PLAN

Elles peuvent se mettre sous la forme suivante :

( ³ 2 ´
∂Ω ∂Ω ∂Ω ∂ Ω ∂2Ω
∂t
+ u 1 ∂x + u 2 ∂y = υ ∂x2
+ ∂y 2
∂2ψ ∂2ψ
∂x2
+ ∂x2
= −Ω
où u1 = ∂ψ∂y
et u2 = − ∂ψ ∂x
.
Ici on ne connaît pas u1 et u2 comme dans le cas de l’équation d’advection-diffusion. On
peut alors considérer le problème comme un système quasi-linéaire et le mettre sous la forme
 ∂Ω(t) ³ 2 ´
(t−k) ∂Ω(t) (t−k) ∂Ω(t) ∂ Ω(t) ∂ 2 Ω(t)

 ∂t + u 1 ∂x
+ u 2 ∂y
= υ ∂x 2 + ∂y 2

∂ 2 ψ(t) ∂ 2 ψ(t)
 + = −Ω(t − k)
 ∂x2 ∂ψ ∂x2
u1 = ∂y , u2 = − ∂ψ
∂x
.
La discrétisation de ce nouveau système différentiel peut alors se faire comme dans le cas
de l’équation d’ advection-diffusion. Par exemple, en utilisant la discrétisation centrée en
espace, on obtient le schéma discret suivant :

(n+1)
Ωi,j
(n)
−Ωi,j
(n+1) (n) (n+1) (n) h i
(n−1) Ωi+1,j −Ωi−1,j (n−1) Ωi,j+1 −Ωi,j−1 υk (n) (n) (n) (n) (n)
k
+ u1 2h
+ u2 2h
= h2
Ωi+1,j + Ωi−1 + Ωi,j+1 + Ωi,j−1 − 4Ωi,j
(n) (n) (n) (n) (n)
ψi+1,j + ψi−1,j + ψi,j+1 + ψi,j−1 − 4ψi,j = −h2 Ωi,j

Pour ce qui est de la discrétisation des conditions aux limites sur la frontière (S), on a :

~u = ~us = us1~i + us2~j

Supposons que la frontière (S) est parallèle à l’axe des x. On a alors sur (S) les conditions
suivantes sur la fonction de courant ψ
(
∂ψ(x,0)
∂x
= −us2 (x) ∂ 2 ψ(x, 0) dus2 (x)
∂ψ(x,0) ⇒ =−
∂y
= us1 (x) ∂x 2 dx
∂ 2 ψ(x,0)
Pour déterminer les valeurs aux limites de Ω sur (S), il nous faut connaître ∂y 2
. Or
nous avons
2 2
ψ(x, h) = ψ(x, 0) + h ∂ψ(x,0) + h2 ∂ ψ(x,0)
h ∂y ∂y 2 i
∂ 2 ψ(x,0) h2 ∂ψ(x,0)
⇒ ∂y2 = 2 ψ(x, h) − ψ(x, 0) − h ∂y
∂ 2 ψ(x,0) ∂ 2 ψ(x,0)
Ainsi connaissant ∂x2
et ∂y 2
, on peut en déduire la valeur de Ω sur (S).

7.3 EQUATIONS DE BURGERS

7.3.1 Présentation de l’équation

Introduite en 1948, par Burgers afin de généraliser les équations non linéaires fréquemment
rencontrées en écoulement des fluides, elle est décrite par

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DE NAVIER STOKES ET DE BURGERS 63

∂u ∂u ∂ 2u
+u = µ 2,
∂t ∂x ∂x
C’est une équation parabolique. On la rencontre dans les problèmes de couches limites où
les équations de Navier Stokes se réduisent à

∂u ∂u ∂ 2u
+u =µ 2
∂x ∂y ∂y
Lorsque le terme de viscosité est absent, l’équation de Burgers se réduit à la forme hy-
perbolique suivante
∂u ∂u
+u =0
∂t ∂x
Cette dernière forme est l’équivalent de l’équation d’Euler de l’écoulement de fluide non
visqueux. De plus, c’est un modèle rencontré en dynamique des gaz dans les tubes ou les
tuyères à section variable. Elle décrit également le son violent que traîne tout avion super-
sonique et qui provoque une gène insupportable pour l’environnement.

7.3.2 Schémas de discrétisation de l’équation de Burgers sans terme de viscosité

Elle peut se mettre sous la forme :
∂u ∂F
+ =0 (1)
∂t ∂x
u2
où F = 2
. Elle peut également sous la forme
∂u ∂u
+A =0 (2)
∂t ∂x
∂F
avec A = ∂x
.
a) Schéma de Lax (1954)
La discrétisation de la forme (1) de l’équation de Burgers avec les différences centrées dans
l’espace donne :
k
ui,j+1 = ui,j −
(Fi+1,j − Fi−1,j ).
2h
Afin d’améliorer la précision de ce schéma explicite, Lax a proposé de remplacer ui,j par
(ui+1,j + ui−1,j )/2. Le schéma de Lax est alors
ui+1,j + ui−1,j k
ui,j+1 = − (Fi+1,j − Fi−1,j )
2 2h
b) Schéma du 2nd ordre en temps
il consiste à discrétiser la partie temporelle par les différences centrées avec une erreur
d’ordre o(k 2 ).
k
ui,j+1 = ui,j−1 − (Fi+1,j − Fi−1,j )
2h
c) Schéma de Lax-Wendroff (1960)

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

RESOLUTION NUMERIQUE DES EQUATIONS DE NAVIER STOKES ET DE BURGERS 64

Le schéma de Lax ci-dessus est un schéma du 1er ordre. Une amélioration a été faite en 1960
par Lax et Wendroff par l’établissement d’un schéma du 2nd ordre. La procédure utilisée est
la suivante :

∂u(x, t) k 2 ∂ 2 u(x, t)
u(x, t + k) = u(x, t) + k +
∂t 2 ∂t2
Et comme µ ¶
∂u ∂F ∂ 2u ∂ ∂F
=− ⇒ 2 =− .
∂t ∂x ∂t ∂x ∂t
De même
½ ½ ∂F
∂u
∂t
= −A ∂u
∂x ⇒ ∂t
= −A¡∂F
∂x ¢
∂F ∂2F
∂t
= A ∂u
∂t ∂t2
∂
= ∂x A ∂F
∂x
Il vient alors :
µ ¶
∂F k2 ∂ ∂F
u(x, t + k) = u(x, t) + k + A
∂x 2 ∂x ∂x
Et finalement

µ ¶2
k 1 k £ ¤
ui,j+1 = ui,j − (Fi+1,j − Fi−1,j ) + Ai+1/2,j (Fi+1,j − Fi,j ) − Ai−1/2,j (Fi,j − Fi−1,j )
2h 2 h
avec Ai+1/2,j = (ui,j + ui+1,j )/2 et Ai−1/2,j = (ui,j + ui−1,j )/2.
d) Autres Schémas
Il existe d’autres schémas de discrétisation de l’équation de Burgers tels que ceux de Mac
Cormack, Rusanov, Burstein-Mirin, Warming-Kulter, etc. Parmi ces schémas, le plus utilisé
est celui de Mac Cormack qui est de type prédicteur correcteur. Sous sa forme simple, il est
décrit par le schéma suivant :
Prédicteur :
ui,j+1 = ui,j − hk (Fi+1,j − Fi,j )
Correcteur £ : ¡ ¢¤
ui,j+1 = 12 ui,j + ui,j+1 − hk F i,j+1 − F i−1,j+1 .
Les barres indiquent les valeurs prédites.

7.3.3 Schéma de discrétisation de l’équation de Burgers visqueux

a) Schéma de Dufort-Frankel
Ici la discrétisation se met sous la forme :
ui,j+1 − ui,j−1 ui+1,j − ui−1,j ui+1,j − ui,j+1 − ui,j−1 + ui−1,j
+ Ai,j =µ
2k 2h h2
Ici on a remplacé ui,j de la discrétisation par le terme (ui,j+1 + ui,j−1 )/2.
b) Autres Schémas
Comme dans le cas de l’équation sans terme visqueux, de nombreux autres schémas ont
été élaborés pour l’équation de Burgers visqueux (il existe également un schéma de type
prédicteur correcteur de Mac Cormack).

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

 EXERCICES

EXERCICE 1 : La répartition de l’énergie rayonnée par une étoile en fonction de la longueur

d’onde est donnée par la loi de Planck
µ µ ¶ ¶
5 hc
I(λ) = (8πhc) / λ exp −1
λkT
où h, k et c sont respectivement les constantes de Planck, de Boltzmann et la célérité de
la lumière. On pose x = hc/λkT . D’après les résultats expérimentaux, I est maximale pour
xm .
1. Etablir l’équation vérifiée par xm .
2. On veut trouver xm par la méthode de Newton-Raphson. Etablir l’algorithme permet-
tant de calculer la valeur numérique de xm .

EXERCICE 2 : L’équation de Van der Waals pour une mole de gaz de masse M est
a
(P + 2 )(V − b) = RM T
V
a, b et R sont des constantes. P , T et V sont respectivement la pression, la température
et le volume. Etablir un algorithme qui permet de calculer le volume V chaque fois que l’on
connaît la température et la pression.

EXERCICE 3 : En 1976, Desantis a établi que le facteur de compressibilité b des gaz réels
vérifie la relation
z = (1 + y + y 2 − y 3 )/(1 − y 3 )
où y = b/4v , v est le volume molaire et z est une constante. On veut trouver b par la
méthode de Newton-Raphson.
1. Expliquer le principe de la méthode.
2. Etablir alors un algorithme permettant de résoudre le problème.
3. Traduire cet algorithme en Fortran.

EXERCICE 4 : Une grandeur y vérifie l’équation suivante :

dy
= −ay − by 4 + c avec y (0) = y0
dx
a , b et c sont des constantes . On se propose de résoudre cette équation numériquement par
la méthode de prédicteur-correcteur (Adam-BasForth du second ordre associé à la méthode
de Runge-Kutta d’ordre 2).

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EXERCICES 66

1. Etablir le schéma itératif correspondant.

2. Etablir l’algorithme de résolution de l’équation.
3. Traduire cet algorithme en Pascal et en Fortran.

EXERCICE 5 : Dans le domaine x ∈ ] 0 , 1 [, une grandeur vérifie l’équation différentielle

y 00 + xa y 0 + b ecy = 0
y 0 (0) = 0; y (1) = 1
a , b et c sont des constantes. On veut résoudre cette équation par la méthode de tir
(associé à la méthode de Runge-Kutta d’ordre 4)
1. En utilisant la méthode de Newton (pour les équations algébriques non linéaires), donner
le schéma itératif qui permet de calculer les valeurs successives de la condition de tir en
x = 0.
2. Etablir la méthode itérative de Runge - Kutta d’ordre 4
3. Etablir l’algorithme complet (incluant le schéma des conditions initiales) permettant de
résoudre le problème.

EXERCICE 6 : Dans l’espace KA entre une cathode et une anode A parallèle à K, le

potentiel à la distance x de K satisfait l’équation différentielle
V 00 = aV −1/2 (1)
où a est une constante. On se propose de trouver numériquement les valeurs de V .
1. Dans cette première question, on suppose que le potentiel et le champ électrique sont
nuls sur la cathode(problème à valeurs initiales). On résout alors (1) par la méthode de
Runge-Kutta d’ordre 2.
(a) Sachant que RK2 dérive de la formule intégrale des trapèzes, donner l’expression
générale des itérations RK2 pour une équation différentielle du premier ordre.
(b) L’équation (1) peut se décomposer en deux équations différentielles du premier
ordre. Etablir l’algorithme de résolution numérique de (1) par RK2.
(c) Traduire cet algorithme en Pascal et en Fortran.
2. Dans cette deuxième question, on suppose connues les conditions aux limites VK = 0
et VA = V 0. la distance KA = l. Le problème peut alors se résoudre par la méthode de
tir.
(a) Expliquer le principe de la méthode de tir.
(b) Ecrire alors l’algorithme de résolution de (1) par la méthode de tir.
(c) Traduire cet algorithme en Pascal et en Fortran.

EXERCICE 7 : Soit une équation de la forme :

y 000 = f (x, y, y 0 , y 00 ) ; y (a) = y1 , y (b) = y2 .
On veut résoudre cette équation par la méthode de tir.
1. Etablir le problème aux valeurs initiales résultant.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EXERCICES 67

2. Le problème aux valeurs initiales résultant doit être résolu par la méthode prédicteur-
correcteur.
(a) Expliquer le principe de la méthode prédicteur - correcteur sur une équation diffé-
rentielle du premier ordre.
(b) En déduire l’algorithme de résolution du problème différentiel du troisième ordre
ci-dessus.

EXERCICE 8 :
1. On considère une fonction y de classe C 4 sur un intervalle [a, b].
(a) En utilisant la formule d’interpolation polynomiale de dégré 2, établir que y 0 =
(y(x + h) − y(x − h))/2h où h est le pas de dérivation.
(b) Déduire de (a) une expression pour y 00
2. Soit le problème aux valeurs aux limites y 00 = p(x)y 0 + q(x)y + g(x) défini sur [a, b] avec
y(a) = c et y(b) = d. On veut résoudre ce problème par la méthode des différences
finies. Etablir à partir de question 1. le système algébrique découlant du problème.
3. On veut résoudre le système algébrique ci-dessus par la méthode itérative de Jacobi.
(a) Expliquer le principe de la méthode itérative de Jacobi
(b) Ecrire alors l’algorithme permettant de résoudre le problème.
(c) Traduire l’algorithme en Pascal et en Fortran.

EXERCICE 9 : Une particule de masse M se déplace dans un potentiel de la forme :

1 1
u (y) = ky 2 + k1 y 4 .
2 4
Le milieu présente une viscosité linéaire de coefficient α et la particule est soumise à une
force extérieure f (t) = F cos ωt.
1. Etablir l’équation du mouvement de la particule et adimensionner cette équation. (rendre
l’équation sans dimension).
2. On donne y (0) = y 0 (0) = 0 et on veut utiliser la méthode de Runge-Kutta d’ordre 4
pour trouver les valeurs de y et ẏ en fonction du temps. Etablir l’algorithme de résolution
numérique de l’équation adimensionnée trouvée à la question 1.
3. A partir de la résolution numérique, on obtient un ensemble de valeurs discrètes yi =
y (ti ). L’expression analytique y(t) étant inconnue, on propose la formule approchée
suivante :
N
X
y (t) = (an cos n ω t + bn sin n ω t)
n=1

(a) Expliquer comment on peut utiliser la méthode des moindres carrés pour déterminer
les valeurs des coefficients an et an .
(b) Etablir les algorithmes permettant de déterminer ces coefficients.
(c) Traduire l’algorithme en Pascal et en Fortran.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EXERCICES 68

EXERCICE 10 : On considère un système algébrique linéaire n × n de la forme

AX = b
où A = (aij ), b = (bj ) avec i, j = 1, 2, . . . , n. On veut résoudre ce système par la méthode
d’élimination de Gauss ou du pivot .
1. Indiquer le principe de la méthode.
2. Etablir l’algorithme correspondant.
3. De l’algorithme ci-dessus, on obtient un système triangulaire. Etablir l’algorithme de
résolution du système triangulaire obtenu.
4. Traduire l’algorithme en Pascal et en Fortran.

EXERCICE 11 : Le mouvement d’un pendule est décrit par l’équation différentielle

d2 Y
+ ω0 sin Y = 0
dt2
où ω0 est la pulsation propre du pendule et Y le déplacement angulaire.
1. On veut résoudre numériquement cette équation par la méthode de Runge-Kutta 2.
(a) Etablir l’algorithme de résolution de cette équation.
(b) Traduire cet algorithme en Fortran.
2. Lorsque le pendule oscille entre Y0 et −Y0 , sa période est définie par la relation
Zπ/2
2T0 dx
T (Y0 ) = p (2)
π 1 − A(sin x)2
0

où T0 = 2π/ω0 et A = sin2 Y0 /2.

(a) Après avoir établi la formule composite des trapèzes, donner l’algorithme de calcul
de l’intégrale (2)
(b) raduire cet algorithme en Pascal et en Fortran.
3. De l’intégration numérique de l’intégrale (2), on obtient le tableau suivant :
Y0 π/18 π/9 π/6 2π/9
T (Y0 )/T0 1,002 1,008 1,020 1,030
On exprime alors T (Y0 ) sous la forme T (Y0 ) = T0 (a + bY02 ). En utilisant la méthode des
moindres carrés, déterminer les valeurs numériques de a et de b

EXERCICE 13 : On mesure la capacité calorifique Cv d’un corps à volume constant et on

obtient les résultats suivants
T(K) 0.1 0.3 0.5
Cv × 104 cal/[Link] 0.211 0.694 1.361
1. On se propose de déterminer la chaleur Q absorbée par une mole du corps lorsque sa
température passe de 0.1K à 0.5K ; Calculer Q par la formule de Simpson.
2. En réalité, le tableau ci-dessus est une partie d’un tableau plus général contenant de
nombreuses valeurs expérimentales. On propose la formule Cv = a + bT . On veut alors
déterminer les coefficients a et b par la méthode des moindres carrés.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EXERCICES 69

(a) Définir le système d’équations vérifié par a et b.

(b) Etablir des algorithmes permettant d’évaluer numériquement les coefficients du
système d’équations et de calculer les valeurs de a et de b.
(c) Traduire ces algorithmes en Pascal.
(d) On veut calculer Q par la formule composite de Simpson à partir des valeurs expé-
rimentales de Cv . Etablir l’algorithme de calcul.

EXERCICE 14 : On considère l’équation intégrale :

Z b
φ (x) = f (x) + k (x, t) φ (t) dt (1)
a

où k (x, t) est une fonction connue. On subdivise le domaine [ a , b ] en n points équidistants

ti = ih.
1. L’équation (1) peut être approximée par la formule composite des trapèzes. Elle prend
alors la forme :
n
X
φ (x) = f (x) + (b − a) ci k (x, ti ) φi (2)
i=1

où φi = φ (ti ). Donner les expressions des ci

2. Sachant que l’équation (2) est vérifiée en tout point xj = tj , montrer que l’équation
intégrale (1) se réduit au système algébrique :
n
X
φj = fj + (b − a) ci k (tj , ti ) φi (3)
i=1

où les φi et φj sont les inconnues.

3. Le système algébrique (3) peut se mettre sous la forme :

AX = b (4)

(4) où A est une matrice n × n, X et b sont les matrices colonnes de n éléments. Définir
ces matrices.
4. On veut résoudre le système (4) par la méthode de Gauss.
(a) Donner l’algorithme de triangularisation du système.
(b) Ecrire l’algorithme de résolution du système triangularisé.
(c) Traduire ce dernier algorithme en Pascal et Fortran.
5. φ (x) possède un zéro dans l’intervalle [a, b] et on se propose de trouver ce zéro en
utilisant la méthode de bissection. Décrire cette méthode et donner son organigramme.

EXERCICE 15 : A l’instant initialt = 0, la température d’un mur d’épaisseur L vaut T0 .

Sur la surface x = 0, la température varie suivant la loi :
π
(x = 0, t) = T0 + T1 sin t; t>0
2

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EXERCICES 70

Par un mécanisme approprié, la face x = L est maintenue à la température constante TL .

L’équation de la chaleur dans le milieu est alors :
∂T ∂2T
=α 2
∂t ∂x
où α est la diffusivité.
1. On se propose de résoudre numériquement le problème ainsi posé par la méthode des
différences finies progressives.
(a) Donner la forme discrète du problème. Quel est l’ordre de l’erreur de discrétisation ?
(b) Ecrire l’algorithme de résolution du système discret obtenu.
(c) Traduire cet algorithme en Pascal et en Fortran.
2. On veut maintenant utiliser la méthode de Richardson pour résoudre le même problème.
(a) Donner le nouveau système discret et indiquer l’ordre et l’erreur de discrétisation.
(b) En admettant que l’on utilise le schéma discret de la question (1) pour approximer
les températuresTi,1 = T (x = i∆x, ∆t), établir l’algorithme permettant de trouver
la température aux instants ultérieurs.
(c) Traduire cet algorithme en Fortran.

EXERCICE 16 :
2
1. On considère l’équation de la chaleur ∂u ∂t
= a2 ∂∂xu2 définie sur D = R+ × ]0, l[ avec
u(0, t) = u(l, t) = 0 et u(x, 0) = f (x)
u −u
On veut résoudre cette équation par le schéma explicite avec ∂u ∂t
= i,j+1k i,j .
2
Etablir que la condition de stabilité de ce schéma est ah2k = 12 où k et h sont les pas
temporel et spatial.
2. Que devient cette condition de stabilité dans le cas d’un problème plan : ∂u
∂t
= a2 ∆u ?
On prendra h1 et h2 les pas spatiaux.
3. On considère l’équation d’advection-diffusion dans le plan :

∂Ω ∂Ω ∂Ω
+ u1 + u2 = υ∆Ω
∂t ∂x ∂y
On la discrétise par le schéma centré en espace. Etablir que la condition de stabilité est
υk
h2
≤ 14 où k est le pas temporel et h est le pas spatial suivant x et sur y. Que devient
cette condition dans le cas unidimensionnel ?

EXERCICE 17 : Soit l’équation d’onde suivante :

∂2u ∂u 2
2∂ u
+ a = c + du3
∂t2 ∂t ∂x2
∂u
où (t, x) ∈ R+ ∪ ]0, l[ ; u (0, t) = u (l, t) = 0 ; u (x, 0) = f (x) ; ∂t
(x, 0) = g (x).
a , c et d sont des constantes.
1. Définir chaque terme de cette équation.
2. On veut résoudre numériquement cette équation par les méthodes de différences finies
centrées. Donner la forme discrète du problème ainsi posé.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EXERCICES 71

3. Exprimer ui,2 = u (x = ih, t = 2k) en fonction des valeurs des fonctions f et g aux points
i − 2, i − 1, i + 1, et i + 2 et des pas h et k.
4. Etablir un algorithme de résolution du système discret obtenu à la question (2).
5. Traduire l’algorithme en Pascal et en Fortran.

EXERCICE 18 : On considère l’équation hyperbolique

∂ 2u ∂ 2u ∂ 2u
a1 + a2 + a3 = a4
∂t2 ∂x∂t ∂x2
définie sur D = R+ × ]0, l[ avec u(x, 0) = f (x) ; ∂u(x,0)
∂t
= g(x) ; u(0, t) = u(l, t) = 0 (t > 0)
Les coefficients ai sont des fonctions de x, t, u, ∂u
∂t
et ∂u
∂x
.
On veut résoudre cette équation par la méthode des caractéristiques. Etablir la méthode
et l’algorithme correspondant.

EXERCICE 19 : Le dispositif représente la coupe d’un tube rectangulaire dans lequel

s’écoule dans sa partie centrale un liquide avec une pression ps . La pression sur la sur-
face extérieure est nulle. A chaque point de la région hachurée, la distribution des pressions
vérifie l’équation :

∂ 2p ∂ 2p
+ =0
∂x2 ∂y 2
1. Faites un adimensionnement de cette équation et précisez les nouvelles conditions sur
la frontière. La pression adimensionnée est u.
2. (a) En utilisant la méthode des différences finies, établir l’équation discrète vérifiée par
ui,j en tout point (xi = i h, yj = j k).
(b) Discrétiser les conditions aux limites.
3. On se propose de trouver les ui,j par la méthode itérative de Gauss-Seidel
(a) Décrire cette méthode itérative.
(b) Etablir l’algorithme de calcul des ui,j .
4. On suppose que le tube a maintenant une forme cylindrique de rayon intérieur R1 =
400mm et de rayon extérieur R2 = 800mm.
(a) Quelle est alors l’équation vérifiée par la pression p ou la grandeur adimensionnée
u.
(b) Proposer une méthode de discrétisation de cette nouvelle équation.
On donne : L = 1000mm, L1 = 500mm, l = 800mm, l1 = 400mm

EXERCICE 20 : La répartition d’une grandeur physique u sur la plaque représentée vérifie

l’équation de poisson

∂ 2u ∂ 2u
+ = f (x, y)
∂u2 ∂y 2
avec u = g (x, y) sur la frontière.
1. En utilisant la méthode des différences finies, établir l’équation discrète vérifiée par ui,j .

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EXERCICES 72

Fig. 7.1 – exercice 19

Fig. 7.2 – Exercice 20

2. On se propose de trouver les ui,j par la méthode itérative de Jacobi.

(a) Décrire cette méthode.
(b) Etablir l’algorithme de calcul de ui,j .

EXERCICE 21 : La représentation d’une grandeur physique u sur la plaque triangulaire

ci-dessus vérifie l’équation de Laplace :

∂2u ∂ 2u
+ =0
∂x2 ∂y 2
∂u
avec u = 0 sur l’hypoténuse et ∂n
= a sur les autres côtés ( n est la normale par rapport
aux côtés considérés).
1. Discrétiser ce problème par la méthode des différences finies.
2. Etablir l’algorithme de calcul des valeurs discrètes ui,j sur toute la plaque.

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I

EXERCICES 73

Fig. 7.3 – Exerce 21

UY1/FS/DPH/LaMSEBP Cours de Méthodes Numériques, Master I