Table des matières

Table des matières I

I Cours et Exercices 1

1 Approximation et Interpolation 3
1.1 Interpolation linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Interpolation de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Interpolation de Neville-Aitken . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Approximation par les moindres carrés . . . . . . . . . . . . . . . . . . . . 5
1.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Différentiation et Intégration 9
2.1 Différentiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 Dérivées d’ordre 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.2 Dérivées d’ordre 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Intégration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Méthodes simples . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1.1 Méthodes des rectangles . . . . . . . . . . . . . . . . . . 10
2.2.1.2 Méthodes des trapèzes . . . . . . . . . . . . . . . . . . . 10
2.2.1.3 Méthodes de Simpson . . . . . . . . . . . . . . . . . . . 11
2.2.2 Méthodes composées ou composites . . . . . . . . . . . . . . . . . 11
2.2.2.1 Méthodes des rectangles . . . . . . . . . . . . . . . . . . 11
2.2.2.2 Méthodes des trapèzes . . . . . . . . . . . . . . . . . . . 11
2.2.2.3 Méthodes de Simpson . . . . . . . . . . . . . . . . . . . 12
2.3 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Racines d’équations 15
3.1 Méthode de la bissection . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Méthode de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . 16
3.3 Méthode de la sécante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.4 Racines d’un polynôme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
I
II Table des matières

4 Équations différentielles 21
4.1 Problèmes avec conditions initiales . . . . . . . . . . . . . . . . . . . . . . 21
4.1.1 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1.2 Méthode de Runge-Kutta d’ordre 2 . . . . . . . . . . . . . . . . . 22
4.1.3 Méthode de Runge d’ordre 4 . . . . . . . . . . . . . . . . . . . . . 23
4.1.4 Méthode de Runge-Kutta d’ordre 4 . . . . . . . . . . . . . . . . . 23
4.1.5 Méthode d’Adams-Bashforth . . . . . . . . . . . . . . . . . . . . 24
4.1.5.1 Méthodes d’Adams-Bashforth à deux pas . . . . . . . . . 24
4.1.5.2 Méthode d’Adams-Bashforth à trois pas . . . . . . . . . 25
4.1.5.3 Méthode d’Adams-Bashforth à quatre pas . . . . . . . . 25
4.1.6 Méthode prédicteur-correcteur . . . . . . . . . . . . . . . . . . . . 26
4.2 Problèmes avec conditions aux limites . . . . . . . . . . . . . . . . . . . . 26
4.2.1 Méthodes des différences finies . . . . . . . . . . . . . . . . . . . 26
4.2.2 Méthodes des fonctions complémentaires . . . . . . . . . . . . . . 27
4.2.3 Méthode de tir . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.3 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

5 Résolution des systèmes algébrique 33

5.1 Méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.1.1 Méthode de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.1.2 Pivot partiel, pivot total . . . . . . . . . . . . . . . . . . . . . . . . 34
5.1.3 Décomposition LU - Méthode de Crout . . . . . . . . . . . . . . . 37
5.1.4 Cas des systèmes tridiagonales . . . . . . . . . . . . . . . . . . . . 37
5.1.5 Méthode de Cholesky . . . . . . . . . . . . . . . . . . . . . . . . 38
5.2 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.2.1 La méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.2.2 La méthode de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . 40
5.2.3 Convergence des méthodes de Jacobi et de Gauss-Seidel . . . . . . 41
5.3 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

6 Équations aux dérivées partielles 45

6.1 Équation de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.2 Équation de la chaleur . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.2.1 Méthode de différences progressives ou schéma explicite . . . . . . 48
6.2.2 Méthode des différences finies régressives ou méthodes implicites . 50
6.2.3 Méthode de Cranck-Nichelson . . . . . . . . . . . . . . . . . . . . 50
6.2.4 Méthode de Richardson . . . . . . . . . . . . . . . . . . . . . . . 51
6.2.5 Méthode de Dufort-Frankel . . . . . . . . . . . . . . . . . . . . . 51
6.3 L’équation d’onde . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.4 Équation d’advection-diffusion . . . . . . . . . . . . . . . . . . . . . . . . 52
6.4.1 Schéma de discrétisation centré en espace . . . . . . . . . . . . . . 53
6.4.2 Schéma explicite décentré en espace . . . . . . . . . . . . . . . . . 53
Table des matières III

6.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

II Correction des Exercices 57

Exercices du chapitre 1 59

Exercices du chapitre 2 61

Exercices du chapitre 3 63

Exercices du chapitre 4 67
 I
Cours et Exercices
Chapitre 1

Approximation et Interpolation

L’interpolation est une opération mathématique qui consiste à établir l’expression mathé-
matique d’une fonction à partir des valeurs de cette fonction en des points bien précis. Ces
points sont souvent disposes sur un réseau régulier. Si le point à évaluer se trouve à l’in-
térieur d’un intervalle élémentaire constitue de deux points consécutifs du réseau, on parle
d’interpolation. Quand le point à évaluer se trouve en dehors des bornes du réseau, on parle
d’extrapolation.
L’approximation d’une fonction consiste à chercher la fonction la plus proche possible,
selon certains critères, d’une fonction donnée. Dans ce cas, il n’est en général plus imposé
de passer exactement par les points donnes initialement.

1.1 Interpolation linéaire

Soit f une fonction continue sur l’intervalle ra, bs, xi et xi`1 deux points. Dans l’interpo-
lation linéaire, on considéré que les points consécutifs du réseau sont reliés par des segments
de droite de la forme :
Ppxq “ ax ` b où a et b sont à determiner (1.1)
On obtient : "
axi ` b “ Ppxi q “ f pxi q
axi`1 ` b “ Ppxi`1 q “ f pxi`1 q
On résout facilement ce système et on obtient :
f pxi`1 q ´ f pxi q xi`1 f pxi q ´ xi f pxi`1 q
a“ et b“
xi`1 ´ xi xi`1 ´ xi
On a donc : „  „ 
f pxi`1 q ´ f pxi q xi`1 f pxi q ´ xi f pxi`1 q
Ppxq “ x` (1.2)
xi`1 ´ xi xi`1 ´ xi

1.2 Interpolation de Lagrange

L’interpolation de Lagrange est une généralisation de l’interpolation linéaire à un po-
lynôme d’ordre n. Ce qui nécessite l’utilisation de (n+1) points consécutifs. Le polynôme
3
4 Approximation et Interpolation

d’interpolation de Lagrange, dont la valeur coïncide avec f aux points xi , c’est-a-dire véri-
fiant Pn pxi q “ f pxi q est donnée par la formule :
n ź n
ÿ px ´ x j q
Pn pxq “ f pxi q (1.3)
i“0 j“0
px i ´ x j q
j‰i

L’algorithme de calcul des polynômes de Lagrange est donnée par :

Algorithme lagrange
Début
SÐ0
Pour i allant de 0 à n Faire
PÐ1
Pour j allant de 0 à n Faire
Si (j‰i) Alors
px ´ x j q
PÐP
pxi ´ x j q
FinSi
FinPour
S Ð S ` P f pxi q
FinPour
Fin

1.3 Interpolation de Neville-Aitken

L’algorithme de Neville-Aitken permet de calculer le polynôme d’interpolation de La-
grange sur n points à partir d’une expression portant (n-1) points.
Soient f0 , f1 , f2 , . . . fn , les valeurs connues de f aux points x0 , x1 , x2 , . . . xn . On a le schéma
suivant :
f0 “ P0,0 P1,0 P2,0 ¨ ¨ ¨ Pi,0 ¨ ¨ ¨ Pn,0
f1 “ P0,1 P1,1 P2,1 ¨ ¨ ¨
f2 “ P0,2 P1,2 ¨ ¨ ¨
.. .. .. pour i “ 0, ¨ ¨ ¨ , n (1.4)
. . .
.. ..
. . P1,n´1
fn “ Pn,0
Pn,0 est la valeur recherchée. Ce schéma s’obtient par la formule suivante :
"
pxi` j`1 ´ xq px j ´ xq i “ 0, ¨ ¨ ¨ , n
Pi`1, j pxq “ Pi, j ´ Pi, j`1 avec (1.5)
pxi` j`1 ´ x j q pxi` j`1 ´ x j q j “ 0, ¨ ¨ ¨ , n
1.4 Approximation par les moindres carrés 5

Algorithme Neville-Aitken
Début
Pour i allant de 0 à n Faire
P0,i Ð fi
FinPour
Pour i allant de 0 à n-1 Faire
Pour j allant de 0 à n-i-1 Faire
pxi` j`1 ´ xq px j ´ xq
Pi`1, j Ð Pi, j ´ Pi, j`1
pxi` j`1 ´ x j q pxi` j`1 ´ x j q
FinPour
FinPour
Fin

Remarque 1.1
Notons que l’algorithme de Neville-Aitken est mieux adapté pour connaitre la valeur du
polynôme d’interpolation de Lagrange en un point donné.

1.4 Approximation par les moindres carrés

Soit f pxq une fonction dont on connait les valeurs aux points x1 , x2 , . . . , xm . La méthode
des moindres carrés permet d’approximer la fonction f pxq par par une autre fonction Ppxq.
En posant f pxi q “ fi , le principe consiste donc à minimiser le reste :
m
ÿ
L“ r fi ´ Pn pxi qs2 (1.6)
i“1

Ce qui peut se faire en posant :

BL
“0 avec k “ 0, . . . , n (1.7)
Bak
où les ak sont les paramètres de la fonction d’approximation Pn pxq.
Posons par exemple :
n
ÿ
Pn pxq “ a0 g0 pxq ` ¨ ¨ ¨ ` an gn pxq “ ak gk pxq avec n!m
k“0

Les gk peuvent être de nature différente. On a donc :

« ff2
ÿm n
ÿ
L“ fi ´ ak gk pxi q
i“0 k“0
m n m
BL ÿ ÿ ÿ
“ 0 ðñ g j pxi q ak gk pxi q “ g j pxi q fi avec j “ 0, . . . , n
Bak ,k“0,...,n i“1 k“0 i“1
6 Approximation et Interpolation

Ce qui donne le système matriciel suivant :

» fi » fi » řm fi
A00 A01 ¨ ¨ ¨ ¨ ¨ ¨ A0n a0 i“1 g0 pxi q f i
— A10 A11 ¨ ¨ ¨ ¨ ¨ ¨ A1n ffi — řm
a1 ffi — i“1 g1 pxi q fi
ffi — ffi
— .. .. .. ffi — .. ..
— ffi — ffi
— . . . ffi — . . (1.8)
ffi — ffi
ffi “ — ffi
— . . .. ffi .. ..
– .. ..
— ffi — ffi
. fl – . fl –
řm .
fl
An0 An1 ¨ ¨ ¨ ¨ ¨ ¨ Ann an i“1 gn pxi q f i

où les éléments de la matrice du système sont donnés par :

m "
ÿ j “ 0, . . . , n
A jk “ g j pxi qgk pxi q avec (1.9)
k “ 0, . . . , n
i“1

L’algorithme de calcul des A jk est donné par :

Début
Pour j allant de 0 à n Faire
Pour k allant de 0 à n Faire
A jk Ð 0
Pour i allant de 1 à m Faire
A jk Ð A jk ` g j pxi qgk pxi q
FinPour
FinPour
FinPour
Fin

Le système matriciel peut être résolu par la méthode de Jacobi ou de Gauss-Seidel.

1.5 Exercices
Exercice 1.1
Soit la table de la fonction erreur donnée ci-dessus :
x 0,0 0,4 0,8 1,2 1,6
f 0,0 0,428392 0,742101 0,910314 0,970348
1. Écrire un programme fortran qui utilise l’interpolation de Lagrange pour évaluer la
fonction erreur aux points : x=0,3 et x=0,5
2. Comparer vos résultats avec les valeurs théoriques :
f p0, 3q “ 0, 328627 et f p0, 5q “ 0, 520500

Exercice 1.2
Les résultats d’une expérience sont donnés dans le tableau suivant :
1.5 Exercices 7

x 1 2 3 4 5
f 1,6 2 2,3 2,4 2,5
On se propose de déterminer l’expression y “ f pxq. Pour cela on choisit les expressions
suivantes :
b b
y “ a` y “ a ` ` cex
x x
En utilisant la méthode des moindres carrés, calculer les coefficients a, b et c pour chaque
approximation.
Chapitre 2

Différentiation et Intégration

2.1 Différentiation

Soit f pxq une fonction de classe Cn`1 sur un intervalle I, on a le développement de Taylor
autour du point x0 suivant :

1 px ´ x0 q2 2 px ´ x0 qn pnq
f pxq “ f px0 q ` px ´ x0 q f px0 q ` f px0 q ` ¨ ¨ ¨ ` f (2.1)
2 n
La fonction f pxq peut être estimée à l’aide d’un polynôme d’interpolation Pn pxq de degré n.
On a :
f pxq » Pn pxq

2.1.1 Dérivées d’ordre 1

On admet les formules de dérivation suivante :

f px ` hq ´ f pxq
f 1 pxq “ ` 0phq difference avant ou progressive d’odre 1
h

f pxq ´ f px ´ hq
f 1 pxq “ ` 0phq difference arriere ou regressive d’ordre 1
h
´ f px ` 2hq ` 4 f px ` hq ´ 3 f pxq
f 1 pxq “ ` 0ph2 q différence progressive d’ordre 2
2h
f px ` hq ´ f px ´ hq
f 1 pxq “ ` 0ph2 q différence centrée d’ordre 2
2h
3 f pxq ´ 4 f px ´ hq ` f px ´ 2hq
f 1 pxq “ ` 0ph2 q différence regressive d’ordre 2
2h
f px ´ 2hq ´ 8 f px ´ hq ` 8 f px ` hq ´ f px ` 2hq
f 1 pxq “ ` 0ph4 q diff centrée d’ordre 4
12h
9
10 Différentiation et Intégration

2.1.2 Dérivées d’ordre 2

On admet les formules de dérivation suivante :

f px ´ 2hq ´ 2 f px ´ hq ` f pxq
f 2 pxq “ ` 0phq difference regressive d’ordre 1
h2
f px ` 2hq ´ 2 f px ` hq ` f pxq
f 2 pxq “ ` 0phq difference progressive d’ordre 1
h2
f px ` hq ´ 2 f pxq ` f px ´ hq
f 2 pxq “ ` 0ph2 q difference centrée d’ordre 2
h2
´ f px ` 2hq ` 16 f px ` hq ´ 30 f pxq ` 16 f px ´ hq ´ f px ´ 2hq
f 2 pxq “ ` 0ph4 q
12h2

2.2 Intégration

2.2.1 Méthodes simples

Soit f une fonction continue sur l’intervalle ra, bs. On cherche a évaluer l’intégrale sui-
vante :
żb
f pxqdx
a

2.2.1.1 Méthodes des rectangles

C’est la méthode la plus simple qui consiste à interpoler la fonction f à intégrer par une
fonction constante. ż b
– si f pxq » f paq, on a : f pxqdx » pb ´ aq f paq
żab
– si f pxq » f pbq, on a : f pxqdx » pb ´ aq f pbq
a

2.2.1.2 Méthodes des trapèzes

En interpolant f par un polynôme de degré 1, les deux points d’interpolation pa, f paqq et
pb, f pbqq suffisent à tracer un segment dont l’intégrale correspond à l’aire d’un trapèze. On
a:
żb
pb ´ aq
f pxqdx » p f paq ` f pbqq (2.2)
a 2
2.2 Intégration 11

2.2.1.3 Méthodes de Simpson

En interpolant f par un polynôme de degré 3, trois points sont nécessaires pour le carac-
a`b
tériser : les valeurs aux extrémités a, b et celle choisie en leur milieu . On a :
2
żb ˆ ˆ ˙ ˙
1 a`b
f pxqdx » pb ´ aq f paq ` 4 f ` f pbq (2.3)
a 6 2

2.2.2 Méthodes composées ou composites

Les méthodes simples proposées ci-dessus sont réputées être très peu précises. Pour amé-
liorer leur précision, on effectue une subdivision régulière de l’intervalle ra, bs, comme suit :

a “ x0 ă x1 ă ¨ ¨ ¨ ă xn´1 ă xn “ b

On a donc :
żb n´1
ÿ ż xi`1
f pxqdx » f pxqdx (2.4)
a i“0 xi

2.2.2.1 Méthodes des rectangles

En posant f pxq » f pxi q pour x P rxi , xi`1 s, on a :

żb n´1
ÿ
f pxqdx » pxi`1 ´ xi q f pxi q (2.5)
a i“0

ou bien en prenant f pxq » f pxi`1 q pour x P rxi , xi`1 s, on a :

żb n´1
ÿ
f pxqdx » pxi`1 ´ xi q f pxi`1 q (2.6)
a i“0

2.2.2.2 Méthodes des trapèzes

Ici on remplace la fonction f par une droite passant par les points pxi , f pxi qq et pxi`1 , f pxi`1 qq.
On a donc :
żb n´1
ÿ f pxi q ` f pxi`1 q
f pxqdx » pxi`1 ´ xi q (2.7)
a i“0
2
12 Différentiation et Intégration

2.2.2.3 Méthodes de Simpson

Dans la méthodes de Simpson, la fonction f est remplacée par un polynôme du second

degré passant par les points : pxi , f pxi qq, pxi`1 , f pxi`1 qq, pxi`2 , f pxi`2 qq. On a :
żb n´1
ÿ pxi`1 ´ xi q ´ ´x ` x ¯ ¯
i`1 i
f pxqdx » f pxi`1 q ` 4 f ` f pxi q (2.8)
a i“0
6 2

D’autres méthodes beaucoup plus précises que celle utilisées ci-dessus ont été élaboré. On
par exemple citer :
– La méthode de Romberg
– Les méthodes de Gauss basées sur les polynômes orthogonaux et qui sont certainement
les plus utilisées et les plus précises.

2.3 Exercices
Exercice 2.1
Établissez les dérivées d’ordre 1 suivantes :
1. différence progressive d’ordre 2
2. différence régressive d’ordre 2
3. différence centrée d’ordre 2
Établissez les dérivées d’ordre 2 suivantes :
1. différence régressive d’ordre 1
2. différence progressive d’ordre 1
3. différence centrée d’ordre 2
Exercice 2.2
Établissez les formules composites des trapèzes et de Simpson.
Exercice 2.3
On mesure la capacité calorifique Cv d’un corps à volume constant et on obtient les résul-
tats suivants :
TpKq 0.1 0.3 0.5
4
Cv ˆ 10 cal{K.mole 0.211 0.694 1.361
1. On se propose de déterminer la chaleur Q absorbée par une mole du corps lorsque la
température passe de 0.1 K à 0.5 K ; 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.
2.3 Exercices 13

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

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

Racines d’équations

Que ce soit en science expérimentale ou en ingénierie, on est souvent amené à résoudre

les équations de la forme :

f pxq “ 0 où f pxq est une fonction quelconque

Nous allons présenter dans ce chapitre, un certain nombre de méthodes qui permettent de
résoudre ce type de problème.

3.1 Méthode de la bissection

Soit f une fonction numérique strictement monotone sur un intervalle ra, bs. On suppose
que l’équation f pxq “ 0 n’a qu’une et une seule solution dans cet intervalle. Les étapes de
cette méthode sont les suivantes :
– poser c “ a`b 2
– si f pcq ˆ f pbq ă 0, alors a “ c sinon b “ c
– on répète ces opérations jusqu’à ce que |a ´ b| ă ε
L’algorithme de cette méthode est le suivant :

Algorithme bissection
Début
TantQue (|a ´ b| ą ε) Faire
c Ð a`b
2
Si ( f pcq ˆ f pbq ă 0) Alors
aÐc
Sinon
bÐc
FinSi
FinTantQue
Fin

15
16 Racines d’équations

3.2 Méthode de Newton-Raphson

La méthode de Newton encore appelée méthode des tangentes est une méthode par ap-
proximation successives. Son schéma itératif est donnée par :
f pxn q
xn`1 “ xn ´ (3.1)
f 1 pxn q
On arrête l’itération lorsque la différence entre deux pas consécutifs est inférieur à la préci-
sion souhaitée |xn`1 ´ xn | ă ε. Si la racine recherchée est double, la méthode converge très
lentement. Pour remédier à ce problème, on pose :
f pxq
f pxq “
f 1 pxq
Remarque 3.1
Le choix d’une valeur très éloignée de la racine recherchée peut amener la méthode a ne
pas converger.
L’algorithme de la méthode de Newton est le suivant :
Algorithme Newton-Raphson
Début
{ on choisit une valeur initiale }
SÐ?
Répéter
PÐS
{ d f est la dérivée de f }
S Ð P ´ p f ppq{d f ppqq
Jusque (|S ´ P| ă ε)
Fin

3.3 Méthode de la sécante

Dans certaines situations, la dérivée de f est très compliquée ou même impossible à ex-
pliciter. On remplace alors f 1 pxq par le taux d’accroissement suivant :
f pxn q ´ f pxn´1 q
f 1 pxq »
xn ´ xn´1
Son schéma itératif est le suivant :
xn´1 f pxn q ´ xn f pxn´1 q
xn`1 “ (3.2)
f pxn q ´ f pxn´1 q
Il faut donc deux valeurs pour initialiser ce schéma.
3.4 Racines d’un polynôme 17

3.4 Racines d’un polynôme

La méthode de Newton-Raphson permet aussi de déterminer les racines d’un polynôme.
Elle est souvent appelée méthode de Newton-Hörner. Écrivons notre polynôme sous la forme
suivante :
Ppxq “ a0 xn ` a1 xn´1 ` a2 xn´2 ` ¨ ¨ ¨ ` an´1 x ` an (3.3)
La division de Ppxq par le monôme x ´ xk où xk est une valeur arbitraire donne :
Ppxq “ px ´ xk qpb0 xn´1 ` b1 xn´2 ` b2 xn´3 ` ¨ ¨ ¨ ` bn´2 x ` bn´1 q ` R
où R est une constante de valeur R “ Ppxk q.
En divisant de nouveau le polynôme de degré n-1 par x ´ xk , on obtient :
b0 xn´1 ` b1 xn´2 ` b2 xn´3 ` ¨ ¨ ¨ ` bn´2 x ` bn´1 “
“ px ´ xk qpc0 xn´2 ` c1 xn´3 ` c2 xn´4 ` ¨ ¨ ¨ ` cn´3 x ` cn´2 q ` S
Donc :
Ppxq “ px ´ xk q2 pc0 xn´2 ` c1 xn´3 ` c2 xn´4 ` ¨ ¨ ¨ ` cn´3 x ` cn´2 q ` px ´ xk qS ` R (3.4)
La valeur de la constante S est donnée par : S “ P1 pxk q. L’application de la formule de
Newton peut donc se mettre sous la forme : xk`1 “ xk ´ RS .
Les valeurs de R et S s’obtiennent donc au moyen des relations de récurrence qu’on établit
en égalant les puissances successives de x dans les diverses expressions du polynôme :
$
’
’ b0 “ a0
’
’
’
’ b1 “ a1 ` xk b0
’
& b
2 “ a2 ` xk b1
.
.. (3.5)
’
’
’
’
’
’ bn´1 “ an´1 ` xk bn´2
’
% R “ a `x b
n k n

De même, les coefficients ci s’obtiennent au moyen des formules analogues :

$
’
’ c 0 “ b0
’
’
’
’ c1 “ b1 ` xk c0
’
& c2 “ b2 ` xk c1
’
’
c3 “ b3 ` xk c2 (3.6)
’ ..
’
’
’ .
’
’
’
’
’ cn´2 “ bn´2 ` xk cn´3
% S “ b `x c
n´1 k n´2

Remarque 3.2
– Le processus ci-dessus est répété pour les polynômes d’ordre inférieur jusqu’à l’obten-
tion de toutes les racines recherchées.
18 Racines d’équations

– D’autres méthodes beaucoup plus efficaces existent. On peut par exemple citer : la
méthode de Bairstow, la méthode de Frobenius.
L’algorithme de la méthode de Newton-Hörner est le suivant :
Algorithme Newton-Horner
Début
Pour k allant de n à 1 Faire
b0 Ð a0
c0 Ð b0
sÐ?
Répéter
pÐs
Pour i allant de 1 à k Faire
bi Ð ai ` pbi´1
ci Ð bi ` pci´1
FinPour
s Ð p ´ c bk
k´1
Jusque |s ´ p| ă ε
Pour i allant de 1 à k-1 Faire
ai Ð ai ` sbi´1
FinPour
FinPour
Fin

3.5 Exercices
Exercice 3.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
8πhc
Ipλ q “ ` ` ˘ ˘
λ 5 exp λhc
kT ´ 1
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. Établir l’équation vérifiée par xm
2. On veut trouver xm par la méthode de Newton-Raphson. Établir l’algorithme permettant
de calculer la valeur numérique de xm

Exercice 3.2
L’équation de Van der Waals pour une mole de gaz est
a
pP ` 2 qpV ´ bq “ RMT
V
3.5 Exercices 19

a, b, et R sont des constantes. P, T et V sont respectivement la pression, la température et

le volume. Établir un algorithme qui permet de calculer le volume V chaque fois que l’on
connait la température et la pression.
Exercice 3.3
En 1976, Desantis a établi que le facteur de compressibilité b des gaz réels vérifie la
relation
z “ p1 ` y ` y2 ´ y3 q{p1 ´ y3 q
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 de Newton
2. Établir alors un algorithme permettant de résoudre le problème
3. Traduire cet algorithme en Fortran
Exercice 3.4
Une poutre de longueur L encastrée aux deux extrémités subit une déformation au temps
t “ 0 et se met par la suite à vibrer. La déformation upx,tq de la poutre à la position x et au
temps t est solution de :
B2u 4
4B u
`c 4 “ 0
Bt 2 Bx
qui est une équation aux dérivées partielles d’ordre 4. Des calculs permettent de montrer que
l’équation caractéristique des fréquences de la poutre est donnée par :
1 ´ cospKLq coshpKLq “ 0
1. Écrire un programme fortran qui permet de calculer les trois premiers modes de vibra-
tions de la poutre
2. comparer avec les résultats présent dans la littérature
4, 730 7, 853 10, 996
K1 “ K2 “ K3 “
L L L
Chapitre 4

Équations différentielles

La résolution numérique des équations différentielles est probablement le domaine de

l’analyse numérique où les applications sont les plus nombreuses. Que ce soit en méca-
nique des fluides, en transfert de chaleur ou en analyse de structures, on aboutit souvent à
la résolution d’équations différentielles. Dans ce chapitre, nous allons exposé les différentes
méthodes permettant de résoudre les problèmes avec conditions initiales et les problèmes
avec conditions aux limites.

4.1 Problèmes avec conditions initiales

On cherche à résoudre le problème suivant, appelé problème de Cauchy :
" 1
y pxq “ f px, ypxqq
(4.1)
yp0q “ y0

4.1.1 Méthode d’Euler

Dans la méthode, on remplace la dérivée première y1 pxq par la formule des différences
finies progressives. On a donc :
ypx ` hq ´ ypxq
où h est le pas de discretisation
h
Ce qui donne :
ypx ` hq “ ypxq ` h f px, ypxqq (4.2)
Sous sa forme discrète, on obtient :
yi`1 “ yi ` h f pxi , yi q formule d’euler progressive (4.3)
Une autre technique pour obtenir la formule d’Euler serait d’intégrer l’équation différentielle
sur les intervalles rxi , xi`1 s :
ż xi`1 ż xi`1 ż xi`1
1
y pxqdx “ f px, yqdx ùñ yi`1 “ yi ` f px, yqdx
xi xi xi

21
22 Équations différentielles

or d’après la formule des rectangles, on a :

ż xi`1
f px, yqdx “ pxi`1 ´ xi q f pxi , yi q “ h f pxi , yi q
xi

Ce qui donne finalement :

yi`1 “ yi ` h f pxi , yi q (4.4)

Remarque 4.1
– La méthode d’Euler est très peu utilisée du fait de sa très faible précision
– D’autres variantes de la méthode d’Euler existent, telles que : la méthode d’Euler ré-
gressive, décentrée, . . .

4.1.2 Méthode de Runge-Kutta d’ordre 2

Dans la méthode de Runge-Kutta d’ordre 2, on utilise les formules d’intégration des tra-
pèzes. On a :
ż xi`1
yi`1 “ yi ` f px, yqdx (4.5)
xi

Or ż xi`1
h
f px, yqdx “ r f pxi , yi q ` f pxi`1 , yi`1 qs
xi 2
Ce qui nous donne donc :
h
yi`1 “ yi ` r f pxi , yi q ` f pxi`1 , yi`1 qs
2 (4.6)
h
“ yi ` r f pxi , yi q ` f pxi`1 , yi ` h f pxi , yi qqs
2
Son algorithme est le suivant :
Algorithme RK2
Début
{ on donne la valeur initiale }
y0 “?
Pour i allant de 0 à n-1 Faire
K1 Ð h f pxi , yi q
K2 Ð h f pxi ` h, yi ` K1 q
yi`1 Ð yi ` 21 pK1 ` K2 q
xi`1 Ð xi ` h
FinPour
Fin
4.1 Problèmes avec conditions initiales 23

4.1.3 Méthode de Runge d’ordre 4

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

h” ´x ` x
i`1 i
¯ ı
yi`1 “ yi ` f pxi , yi q ` 4 f , yi` 1 ` f pxi`1 , yi`1 q (4.7)
6 2 2

Or

xi`1 “ xi ` h
h
yi` 1 “ yi ` f pxi , yi q
2 2
yi`1 “ yi ` h f pxi ` h, yi ` h f pxi , yi qq

On obtient donc :
„ 
h h h
yi`1 “ yi ` f pxi , yi q ` 4 f pxi ` , yi ` f pxi , yi qq ` f pxi ` h, yi ` h f pxi , yi qq (4.8)
6 2 2
En posant :
L1 “ h f pxi , yi q
L2 “ h f pxi ` h, yi ` L1 q
L3 “ h f pxi ` h, yi ` L2 q
L4 “ h f pxi ` 2h , yi ` L21 q
On obtient :
1
yi`1 “ yi ` pL1 ` 4L4 ` L3 q (4.9)
6

4.1.4 Méthode de Runge-Kutta d’ordre 4

C’est une amélioration de la méthode de Runge du 4ème ordre. Elle fait usage du dévelop-
pement en série d’ordre 4 de ypx ` hq. Son schéma itératif est donné par :
1
yi`1 “ yi ` pL1 ` 2L2 ` 2L3 ` L4 q (4.10)
6
avec
L1 “ h f pxi , yi q
L2 “ h f pxi ` 2h , yi ` L21 q
L3 “ h f pxi ` 2h , yi ` L22 q
L4 “ h f pxi ` h, yi ` L3 q
Son algorithme est donnée par :
24 Équations différentielles

Algorithme RK4
Début
x0 Ð ? y0 Ð ? h Ð ?
Pour i allant de 0 à n-1 Faire
L1 Ð h f pxi , yi q
L2 Ð h f pxi ` h2 , yi ` L21 q
L3 Ð h f pxi ` h2 , yi ` L22 q
L4 Ð h f pxi ` h, yi ` L3 q
yi`1 Ð yi ` 61 pL1 ` 2L2 ` 2L3 ` L4 q
xi`1 Ð xi ` h
FinPour
Fin

4.1.5 Méthode d’Adams-Bashforth

Ces méthodes sont aussi appelées méthodes à pas multiples par opposition aux méthodes
précédentes qui sont des méthodes à un pas. Les méthodes d’Adams-Bashforth sont bien
adaptées aux systèmes à faibles non linéarités mais dont les paramètres ne varient pas beau-
coup dans le temps.

4.1.5.1 Méthodes d’Adams-Bashforth à deux pas

On veut connaitre yk`1 connaissant yk et yk´1 . Pour cela on fait les développements sui-
vants :
1 h2 2
yk`1 » yk ` hyk ` yk (4.11)
2
y1k ´ y1k´1
y1k´1 » y1k ´ hy2k ùñ y2k » (4.12)
h
En remplaçant p4.12q dans p4.11q, on obtient :
h
yk`1 “ yk ` r3y1k ´ y1k´1 s (4.13)
2
Or y1k “ f pxk , yk q et y1k´1 “ f pxk´1 , yk´1 q
h
ùñ yk`1 “ yk ` r3 f pxk , yk q ´ f pxk´1 , yk´1 qs (4.14)
2
Remarque 4.2
Sa mise en oeuvre nécessite la connaissance de yk et yk´1 .
4.1 Problèmes avec conditions initiales 25

4.1.5.2 Méthode d’Adams-Bashforth à trois pas

Pour obtenir cette méthode, on utilise les développements suivants :

h2 2 h3 3
yk`1 » yk ` hy1k `y ` y (4.15)
2 k 6 k
1 1 2 h2 3
yk´1 » yk ´ hyk ` yk (4.16)
2
yk´2 » yk ´ 2hyk ` 2h2 y3
1 1 2
k (4.17)
p4.16q et p4.17q permettent d’obtenir le système suivant :
h2 3 2 1 1
2 yk ´ hyk “ yk´1 ´ yk
2h2 y3 2 1 1
k ´ 2hyk “ yk´2 ´ yk
La résolution de ce système nous donne :
3{2y1k ` 1{2y1k´2 ´ 2y1k y1k ´ 2y1k´1 ` y1k´2
“ y2k 3
et yk “
h h
Le remplacement de ces relations dans p4.15q donne :
h “ 1 ‰
yk`1 “ yk ` 23yk ´ 16y1k´1 ` 5y1k´2 (4.18)
12
h
“ yk ` r23 f pxk , yk q ´ 16 f pxk´1 , yk´1 q ` 5 f pxk´2 , yk´2 qs (4.19)
12
Remarque 4.3
Ce schéma itératif nécessite le calcul des valeurs yk , yk´1 , yk´2

4.1.5.3 Méthode d’Adams-Bashforth à quatre pas

On utilise les développements suivants :

h2 2 h3 3 h4 4
yk`1 » yk ` hy1k ` y ` y ` y (4.20)
2 k 6 k 24 k
h2 h3 4
y1k´1 » y1k ´ hy2k ` y3 ´ y (4.21)
2 k 6 k
4h3 4
y1k´2 » y1k ´ 2hy2k ` 2h2 y3 k ´ y (4.22)
3 k
1 1 2 9h2 3 9h3 4
yk´3 » yk ´ 3hyk ` y ´ y (4.23)
2 k 2 k
Un raisonnement identique à la méthode précédente permet d’obtenir :
h
yk`1 “ yk ` r55 f pxk , yk q ´ 56 f pxk´1 , yk´1 q ` 37 f pxk´2 , yk´2 q ´ 9 f pxk´3 , yk´3 qs (4.24)
24
Remarque 4.4
La mise en oeuvre de ce schéma nécessite le calcul de yk , yk´1 , yk´2 et yk´3
26 Équations différentielles

4.1.6 Méthode prédicteur-correcteur

Les formules de Runge, Runge-Kutta sont dites implicites en ce sens que les relations qui
permettent d’évaluer yk`1 dépendent de yk`1 lui-même. Pour contourner cette difficulté, on
combine les formules de Runge-Kutta et les formules d’Adams-Bashforth en des schémas
dits prédicteur-correcteurs. Il s’agit simplement d’utiliser les schémas d’Adams-Bashforth
pour obtenir une première approximation de yk`1 , qui est l’étape de prédiction. On fait appel
ensuite aux formules de Runge-Kutta pour corriger et éventuellement améliorer cette ap-
proximation. Par exemple, une méthode de prédicteur-correcteur peut utiliser les formules
suivantes :
h
yk`1 “ yk ` r3 f pxk , yk q ´ f pxk´1 , yk´1 qs predicteur(AB à deux pas) (4.25)
2
h
yk`1 “ yk ` r f pxk , yk q ` f pxk`1 , yk`1 qs correcteur(RK d’ordre 2) (4.26)
2
Son algorithme est le suivant :
Algorithme AB2-RK2
Début
x0 Ð ? h Ð ? y0 Ð ? y1 Ð y0 ` h f px0 , y0 q
Pour k allant de 1 à n-1 Faire
yk`1 Ð yk ` 2h r3 f pxk , yk q ´ f pxk´1 , yk´1 qs
xk`1 Ð xk ` h
yk`1 Ð yk ` 2h r f pxk , yk q ` f pxk`1 , yk`1 qs
FinPour
Fin

4.2 Problèmes avec conditions aux limites

Dans cette section, nous faisons l’étude des équations différentielles linéaires d’ordre 2
avec conditions aux limites de la forme :

y2 pxq “ppxqy1 pxq ` qpxqypxq ` rpxq (4.27)

ypaq “ ya et ypbq “ yb (4.28)

4.2.1 Méthodes des différences finies

L’objectif est de trouver une approximation yi de ypxi q en certains points xi de l’intervalle

ra, bs. On divise d’abord cet intervalle en n sous-intervalles de longueur h “ pb ´ aq{n et on
note x0 “ a, xi`1 “ xi `h et enfin xn “ b. Les conditions aux limites imposent immédiatement
que y0 “ ya et yn “ yb . Il reste donc à déterminer les pn ´ 1q inconnues y1 , y2 , . . . , yn´1 .
4.2 Problèmes avec conditions aux limites 27

Pour ce faire, on construit un système linéaire de dimension pn ´ 1q. Au point xi , l’équation

différentielle s’écrit :
y2i ´ pi y1i ´ qi yi “ ri (4.29)
En posant
yi`1 ´ yi yi`1 ´ 2yi ` yi´1
y1i “ et y2i “ (4.30)
h h2
On obtient l’équation discrète suivante :
p1 ´ hpi qyi`1 ´ p2 ´ hpi ´ h2 qi qyi ` yi´1 “ h2 ri (4.31)
avec i “ 1, . . . , n ´ 1
Ce schéma peut se mettre sous la forme algébrique suivante :
» fi »
α1 β1 0 ¨ ¨ ¨ ¨ ¨ ¨ 0
fi » fi
y1 h2 r1 ´ y0
. .. ffi — y ffi —
— 1 α2 β2 . . . ffi — 2 ffi — 0
— ffi
. . . . . ffi — .. ffi — ..
ffi
— 0 .. .. .. .. .. ffi — . ffi —
—
.
ffi
(4.32)
ffi — . ffi “ — ffi
— .. . . . . . . . . . . . . ..
—
ffi — .. ffi — ffi
— . 0 ffi — . ffi
— .. . . . . ffi
. . . . . . βn´2 fl – .. fl –
ffi —
0
ffi
– . fl
2
h rn´1 ´ βn´1 yn
0 ¨ ¨ ¨ ¨ ¨ ¨ 0 1 αn´1 yn´1

où αi “ ´p2 ´ hpi ´ h2 qi q et βi “ p1 ´ hpi q

Ce système algébrique peut se résoudre par les méthodes itératives de Jacobi ou de Gauss-
Seidel.

4.2.2 Méthodes des fonctions complémentaires

4.2.3 Méthode de tir

Elle est utilisée dans le cas des problèmes aux limites non linéaire. 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 l’équation différentielle suivante :
y2 “ f px, y, y1 q aăxăb (4.33)
avec ypaq “ ya et ypbq “ yb . La méthode de tir consiste à choisir un nombre réel γ tel que :
y1 pαq “ γ et résoudre le problème aux valeurs initiales suivant :
y2 “ f px, y, y1 q ypaq “ ya y1 paq “ γ (4.34)
Les méthodes de résolution numérique des problèmes aux valeurs initiales peuvent donc
s’appliquer ici pour trouver les valeurs de y en tout point du domaine ra, bs. Le choix de γ
est bon lorsque en x “ b, on a :
ypbq » yb (4.35)
28 Équations différentielles

Une procédure à utiliser pour le choix de γ est la suivante : on sait que la valeur ypbq dépend
de γ, c’est-à-dire ypbq “ ypb, γq. Pour un bon choix de γ, on doit avoir
ypb, γq ´ yb » 0 ðñ f pγq “ 0 où f pγq “ ypb, γq ´ yb
or d’après la méthode de Newton, le schéma itératif donnant la solution de f pγq “ 0 est :
f pγn q
γn`1 “ γn ´ (4.36)
f 1 pγn q

ypb, γn q ´ yb
ùñ γn`1 “ γn ´
y1 pb, γn q
pγn ´ γn´1 qpypb, γn q ´ yb q
“ γn ´
ypb, γn q ´ ypb, γn´1 q
Ainsi pratiquement on choisit une première valeur γ0 de manière arbitraire, une deuxième
valeur de γ1 aussi arbitraire, mais assez proche de γ0 et les autres valeurs de γ sont ensuite
évaluées à partir de la méthode de Newton-Raphson. Le code numérique s’arrête lorsque
γ “ γn tel que :
|ypb, γn q ´ yb | ă ε (4.37)
où ε est une précision que nous avons imposée

4.3 Exercices
Exercice 4.1
Le mouvement d’un pendule est décrit par l’équation différentielle
d 2Y
` ω0 sinY “ 0
dt 2
où ω0 est la pulsation propre et Y le déplacement angulaire.
1. On veut résoudre numériquement cette équation par la méthode de Runge-Kutta2
(a) Établir 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
2T0 π{2
ż
dx
T pY0 q “ a
π 0 1 ´ Apsin xq2
où T0 “ 2π{ω0 et A “ sin2 pY0{2q
(a) Après avoir établi la formule composite des trapèzes, donner l’algorithme de calcul
de l’intégrale ci-dessus
(b) Traduire cet algorithme en fortran
4.3 Exercices 29

3. De l’intégration numérique de l’intégrale ci-dessus, on obtient le tableau suivant :

Y0 π{18 π{9 π{6 2π{9
T pY0 q{T0 1,002 1,008 1,020 1,030
On exprime alors T pY0 q sous la forme T pY0 q “ T0 pa ` bY02 q. En utilisant la méthode des
moindres carrés, déterminer les valeurs numériques de a et de b.
Exercice 4.2
Une grandeur y vérifie l’équation suivante :
dy
“ ´ay ´ by4 ` c avec yp0q “ 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 (Adams-Bashforth du second degré associé à la méthode
de Runge-Kutta d’ordre 2)
1. Établir le schéma itératif correspondant
2. Établir l’algorithme de résolution de l’équation
3. Traduire cet algorithme en fortran
Exercice 4.3
Dans le domaine x Ps0, 1r, une grandeur vérifie l’équation différentielle
a
y2 ` y1 ` becy “ 0
x
1
y p0q “ 0 yp1q “ 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 algebriques non linéaires), don-
ner le schéma itératif qui permet de calculer les valeurs successives de la condition de
tir en x “ 0
2. Établir la méthode itérative de Runge-Kutta d’ordre 4
3. Établir l’algorithme complet (incluant le schéma des conditions initiales) permettant de
résoudre le problème.
Exercice 4.4
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.
´1{2
V 2 “ aV
où a est une constante. On se propose de trouver numériquement les valeurs de V .
30 Équations différentielles

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 l’équation différen-
tielle par la méthode de Runge-Kutta d’ordre 2.
(a) Sachant que RK2 dérive de la formule des trapèzes, donner l’expression générale
des itérations RK2 pour une équation différentielle du premier ordre.
(b) L’equation différentielle peut se décomposer en deux équations différentielles du
premier ordre. Établir son algorithme de résolution numérique par RK2.
(c) Traduire cet algorithme en fortran.
2. Dans cette deuxième question, on suppose connues les conditions aux limites VK “ 0 et
VA “ V0 , la distance de 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) Écrire alors l’algorithme de résolution de l’équation différentielle par la méthode de
tir
(c) traduire cet algorithme en fortran.
Exercice 4.5
Soit une équation de la forme
y3 “ f px, y, y1 , y2 q ypaq “ y1 ypbq “ y2
On veut résoudre cette équation par la méthode de tir.
1. Établir le problème aux valeurs initiales résultant
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éren-
tielle du premier ordre.
(b) En déduire l’algorithme de résolution du problème différentiel du troisième ordre
ci-dessus
Exercice 4.6
Une particule de masse M se déplace dans un potentiel de la forme :
1 1
upyq “ ky2 ` k1 y4
2 4
Le milieu présente une viscosité linéaire de coefficient λ et la particule est soumise à une
force extérieure f ptq “ F cos ωt.
1. Établir l’équation du mouvement de la particule et adimensionner cette équation (rendre
l’équation sans dimension).
4.3 Exercices 31

2. On donne yp0q “ y1 p0q “ 0 et on veut utiliser la méthode de Runge-Kutta d’ordre 4 pour

trouver les valeurs de y et y9 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 “
ypti q. L’expression analytique yptq étant inconnue, on propose la formule approchée sui-
vante :
ÿN
yptq “ pan cos nωt ` bn sin nωtq
n“1
(a) Expliquer comment on peut utiliser la méthode des moindres carrés pour déterminer
les valeurs des coefficients an et bn
(b) Établir les algorithmes permettant de déterminer ces coefficients
(c) Traduire l’algorithme en fortran
Chapitre 5

Résolution des systèmes algébrique

5.1 Méthodes directes

5.1.1 Méthode de Gauss

R R
Soit A P nˆn et b P n . On suppose que la matrice A est inversible c’est-à-dire que aii ‰ 0
pour tout i. On veut résoudre le système Ax “ b
¨ ˛¨
a11 a12 ¨ ¨ ¨ a1n´1 a1n
˛ ¨ ˛
x1 b1
˚ a21 a22 ¨ ¨ ¨ a2n´1 a2n ‹ ˚ x2
‹ ‹ ˚ b2 ‹
.. .. .. ‹ ˚
˚
˚
. . . ‹˚ .. ‹ ˚
‹“˚ .. ‹
(5.1)
˚ . .
˚ ‹
.. .. .. ‚˝ x
˚ ‹ ‹ ˚ ‹
˝ . . . n´1
‚ ˝ bn´1 ‚
an1 an2 ¨ ¨ ¨ ¨ ¨ ¨ ann xn bn

par la méthode de Gauss. Le procédé d’élimination de Gauss pour résoudre un système li-
néaire n ˆ n consiste à utiliser la première équation pour exprimer la première inconnue
x1 en fonction des autres x2 , . . . , xn puis à reporter la valeur ainsi trouvée dans les équations
suivantes. On obtient alors un système n ´ 1 ˆ n ´ 1 en les inconnues x2 , . . . , xn auquel on ap-
plique la même méthode. Ce procédé permet d’obtenir, après n ´ 1 telles étapes, un nouveau
système qui est triangulaire et équivalent au premier.

¨ ˛¨ ˛ ¨ ˛
a11 a12 ¨¨¨ a1n´1 a1n x1 b1
˚ 0 a22 ¨¨¨ a2n´1 a2n ‹˚ x2 ‹ ˚ b2 ‹
˚ .. .. .. .. ‹˚ .. ‹ ˚
‹“˚ .. ‹
(5.2)
˚
˚ . . . . ‹˚
‹˚ . ‹ ˚ . ‹
‹
˝ 0 0 ¨ ¨ ¨ an´1n´1 an´1n ‚˝ x
n´1
‚ ˝ bn´1 ‚
0 0 ¨¨¨ 0 ann xn bn

L’algorithme de l’élimination de Gauss est donnée par :

33
34 Résolution des systèmes algébrique

Algorithme elemination de Gauss

Début
Pour k allant de 1 à n-1 Faire
Pour i allant de k+1 à n Faire
bi Ð bi ´ aaik bk
kk
Pour j allant de k+1 à n Faire
ai j Ð ai j ´ aaik ai j
kk
FinPour
FinPour
FinPour
Fin

Les inconnues x1 , x2 , . . . , xn sont déterminées à partir du système triangulaire p5.2q. On cal-

cule successivement xn à partir de la dernière équation, puis xn´1 à partir de l’avant dernière
et ainsi de suite. Ce qui donne
$
’
’ xn “ bn {ann
xn´1 “ pbn´1 ´ an´1,n xn q{an´1,n
&
(5.3)
’
’ ¨¨¨
x1 “ pb1 ´ a12 x2 ´ ¨ ¨ ¨ ´ a1n xn q{a11
%

Cette méthode est appelée méthode de la remontée. Son algorithme est le suivant
Algorithme methode de la remontée
Début
xn Ð bn {ann
Pour i allant de n-1 1 parPasDe -1 Faire
SÐ0
Pour j allant de i+1 à n Faire
S Ð S ` ai j x j
FinPour
xi Ð pbi ´ Sq{aii
FinPour
Fin

Remarque 5.1
L’algorithme complet de résolution est obtenu en associant l’algorithme d’élimination de
Gauss et l’algorithme de la méthode de la remontée.

5.1.2 Pivot partiel, pivot total

La première étape de la méthode d’élimination de Gauss fait jouer un rôle particulier

au coefficient a11 de la matrice A. On l’appelle pivot de la méthode et, pour cette raison, on
parle souvent de la méthode du pivot de Gauss. Le choix de a11 comme pivot n’est pas le seul
5.1 Méthodes directes 35

possible : on obtient système équivalent en permutant entre elles les différentes équations ou
bien en prenant les inconnues dans un ordre différent et, à ces nouveaux systèmes, on peut
aussi appliquer la méthode du pivot de Gauss. Ces différentes stratégies ne sont pas sans
influence du point de vue du calcul des erreurs.
La méthode du pivot partiel consiste à choisir pour pivot le coefficient de plus grand
module de la première colonne de A

|ai1 | “ max |ak1 | (5.4)

1ďkďn

Notons que ai1 ‰ 0 puisque l’on suppose que A est inversible. On pourra donc diviser par
ai1 . On obtient un nouveau système en permutant dans A les lignes 1 et i et en laissant les
autres inchangées. L’algorithme de Gauss complet avec pivot partiel est donné par :

Algorithme Gauss-pivot-partiel
Début
Pour k allant de 1 à n-1 Faire
{ recherche du pivot max }
amax Ð akk i0 Ð k
Pour i allant de k+1 à n Faire
Si (|aik |>|amax|) Alors
amax Ð aik
i0 Ð i
FinSi
FinPour
Pour j allant de k à n Faire
a1 Ð ai0, j
ai0, j Ð ak j
ak j Ð a1
FinPour
b1 Ð bi0 bi0 Ð bk bk Ð b1
{ triangularisation du systeme }
Pour i allant de k+1 à n Faire
rik Ð aik {akk
bi Ð bi ´ rik {bk
Pour j allant de k+1 à n Faire
ai j Ð ai j ´ rik ak j
FinPour
FinPour
FinPour
36 Résolution des systèmes algébrique

{ recherche des solutions }

xn Ð bn {ann
Pour i allant de n-1 1 parPasDe -1 Faire
SÐ0
Pour j allant de i+1 à n Faire
S Ð S ` ai j x j
FinPour
xi Ð pbi ´ Sq{aii
FinPour
Fin

La methode du pivot total prend pour pivot ai j tel que

|ai j | “ max |akl | (5.5)

1ďk,lďn

On doit donc permuter les lignes 1 et i ainsi que les colonnes 1 et j. L’algorithme complet
de Gauss avec pivot total est donné par :

Algorithme Gauss-pivot-total
Début
Pour k allant de 1 à n-1 Faire
{ recherche du pivot }
amax Ð akk ; i0 Ð k ; j0 Ð i
Pour i allant de k à n Faire
Pour j allant de k à n Faire
Si (|ai j |>|amax|) Alors
amax Ð ai j ; j0 Ð j ; i0 Ð i
FinSi
FinPour
FinPour
{ permutation des lignes i0 et k dans A et b }
Pour j allant de k à n Faire
a1 Ð ai0, j ; ai0, j Ð ak j ; ak j Ð a1
FinPour
b1 Ð bi0 ; bi0 Ð bk ; bk Ð b1
{ permutation des colonnes j0 et k dans A et inc }
Pour i allant de k à n Faire
a1 Ð ai, j0 ; ai, j0 Ð aik ; aik Ð a1
FinPour
c1 Ð inc j0 ; inc j0 Ð inck ; inck Ð c1
{ triangularisation du systeme }
{ identique au cas du pivot partiel }
FinPour
5.1 Méthodes directes 37

{ recherche des solutions }

xincn Ð bn {ann
Pour i allant de n-1 1 parPasDe -1 Faire
SÐ0
Pour j allant de i+1 à n Faire
S Ð S ` ai j xinc j
FinPour
xinci Ð pbi ´ Sq{aii
FinPour
Fin

5.1.3 Décomposition LU - Méthode de Crout

Pour résoudre Ax “ b, on cherche à écrire A “ LU où

– L est une matrice triangulaire inférieure avec des 1 sur la diagonale
– U est une matrice triangulaire supérieure
La résolution de Ax “ b est alors ramenée aux résolutions successives des systèmes éche-
lonnes Ly “ b et Ux “ y. On calcule les matrices L et U avec les formules suivantes :
i´1
ÿ
ui j “ai j ´ lik uk j j “ i, . . . , n (5.6)
k“1
j´1
1 ÿ
li j “ pai j ´ lik uk j q i “ j ` 1, . . . , n (5.7)
ujj k“1

Remarque 5.2
Une autre méthode serait plutôt de choisir la matrice U avec des 1 sur la diagonale. Cette
méthode est appelée méthode de Doolittle. Le calcul des matrices L et U se fait avec les
formules suivantes :
i´1
1 ÿ
ui j “ pai j ´ lik uk j q j “ i, . . . , n (5.8)
li j k“1
j´1
ÿ
li j “ai j ´ lik uk j i “ j ` 1, . . . , n (5.9)
k“1

5.1.4 Cas des systèmes tridiagonales

Les méthodes de discrétisation des équations aux dérivées partielles conduisent souvent
à la résolution de systèmes dont les matrices ont des structures bandes, creuses ou par blocs.
Exploiter la structure de la matrice permet une réduction considérable du cout des algo-
rithmes de factorisation et de substitution. En algèbre linéaire numérique, l’algorithme des
38 Résolution des systèmes algébrique

matrices tridiagonales, aussi connu sous le nom d’algorithme de Thomas, est une forme
simplifiée de l’élimination de Gauss utilisée pour résoudre les systèmes tridiagonals d’équa-
tions. Un système tridiagonal pour n ´ 1 inconnues peut s’écrire sous la forme

αi xi´1 ` βi xi ` γi xi`1 “ bi i “ 1, . . . , n (5.10)

où α1 “ 0 et γn´1 “ 0. Sous la forme matricielle, le système s’écrit :

» fi »
β1 γ1 0 ¨ ¨ ¨ ¨ ¨ ¨ 0
fi » fi
x1 b1
... .. ffi —
. ffi — x2 ffi — b2
—
— α2 β2 γ2 ffi — ffi
.. ffi
ffi
— 0 ... ... ... ...
— — x3 ffi — b3
. ffi
ffi
ffi — .. ffi “ — .. (5.11)
ffi — ffi
— .. . . . . . . . . . . . .
—
ffi — . ffi — .
ffi
— . 0 ffi — — . ffi — .
ffi
— .. ... ... ... . .
ffi
γn´2 fl – . fl – .
ffi
– . fl
0 ¨ ¨ ¨ ¨ ¨ ¨ 0 αn´1 βn´1 xn´1 bn´1

La première étape de résolution consiste à modifier les éléments non nuls de la matrice
tridiagonale.
$ γ
1
’
& ; i“1
1 β1
γi “ γi (5.12)
’
% 1 α ; i “ 2, 3, . . . , n ´ 2
βi ´ γi´1 i

et $
b1
’
’
& ; i“1
1 β1
bi “ bi ´ b1i´1 αi (5.13)
’
’
% 1 α ; i “ 2, 3, . . . , n ´ 1
βi ´ γi´1 i

Les solutions du système sont obtenues par les relations suivantes :

xn´1 “b1n´1
(5.14)
xi “b1i ´ γi1 xi`1 ; i “ n ´ 2, n ´ 3, . . . , 1

Remarque 5.3
– L’algorithme de Thomas n’utilise que Opnq opérations contrairement à l’élimination de
Gauss qui nécessite Opn3 q opérations.
– Le système p4.32q du chapitre 4 peut être avantageusement résolu par l’algorithme de
Thomas.

5.1.5 Méthode de Cholesky

La méthode de Cholesky s’applique aux matrices symétriques et définies positives. Elle

consiste à mettre A sous la forme A “ Lt L. En appliquant l’algorithme à la i-ème étape, on
5.2 Méthodes itératives 39

calcule la matrice L par :

j´1
ÿ
li j “pai j ´ lik lk j q{l j j j “ 1, . . . , i ´ 1
k“1
(5.15)
i´1
ÿ
lik2 q {2
1
lii “pai j ´
k“1

On résout ensuite le système Ly “ b puis Lt x “ y par un double balayage

i´1
ÿ
yi “pbi ´ lik yk q{lii i “ 1, . . . , n
k“1
n
(5.16)
ÿ
xi “py j ´ lki xk q{lii
k“i`1

5.2 Méthodes itératives

L’idée des méthodes itératives est de construire une suite de vecteur xpkq qui converge
vers le vecteur x, solution du système Ax “ b. L’intérêt des méthodes itératives, comparées
aux méthodes directes, est d’être simples à programmer et de nécessiter moins de place en
mémoire. En revanche le temps de calcul est souvent plus long

5.2.1 La méthode de Jacobi

On remarque que si les éléments diagonaux de A sont non nuls, le système linéaire Ax “ b
est équivaut à : ¨ ˛
n
1 ÿ
xi “ ˝bi ´ ai j x j ‚ i “ 1, . . . , n (5.17)
aii j“1, j‰i

Pour une donnée initiale xp0q choisie, on calcule xpk`1q par

¨ ˛
n
pk`1q 1 ÿ pkq
xi “ ˝bi ´ ai j x j ‚ i “ 1, . . . , n (5.18)
aii j“1, j‰i

Dans la pratique on arrêtera les itérations lorsque les changements dans le vecteur solution,
c’est-à-dire l’erreur devient suffisamment petite
řn k`1
i“1 |xi ´ xik |
řn k`1
ăε (5.19)
i“1 |xi |
avec ε une tolérance donnée
40 Résolution des systèmes algébrique

Algorithme Jacobi
Début
{ on initialise les valeurs de x0 }
Pour i allant de 1 à n Faire
x0i Ð 0
FinPour
Répéter
s1 Ð 0 ; s2 Ð 0
Pour i allant de 1 à n Faire
sÐ0
Pour j allant de 1 à n Faire
Si (j‰i) Alors
s Ð s+ai j x0 j
FinSi
FinPour
bi ´ s
xi Ð
aii
s1 Ð s1+|xi ´ x0i |
s2 Ð s2+|xi |
x0i Ð xi
FinPour
s1
dÐ
s2
Jusque (d<ε)
Fin

5.2.2 La méthode de Gauss-Seidel

Pour cette méthode, on a :

¨ ˛
i´1 n
pk`1q 1 ÿ pk`1q
ÿ pkq
xi “ ˝bi ´ ai j x j ´ ai j x j ‚ i “ 1, . . . , n (5.20)
aii j“1 j“i`1

L’algorithme de Gauss-Seidel ne nécessite qu’un vecteur de stockage, xpkq étant remplace

par xpk`1q au cours de l’itération. Il est en général plus rapide que l’algorithme de Jacobi,
donc préférable. Dans la pratique, l’itération s’arrête pour :
řn k`1
i“1 |xi ´ xik |
řn k`1
ăε
i“1 i |x |
5.3 Exercices 41

Algorithme Gauss-Seidel
Début
Répéter
s2 Ð 0 ; s3 Ð 0
Pour i allant de 1 à n Faire
sÐ0
Pour j allant de 1 à i-1 Faire
s Ð s+ai j x j
FinPour
s1 Ð 0
Pour j allant de i+1 à n Faire
s1 Ð s1+ai j x0 j
FinPour
bi ´ s ´ s1
xi Ð
aii
s2 Ð s2+|xi ´ x0i |
s2 Ð s3+|xi |
FinPour
s2
dÐ
s3
Pour i allant de 1 à n Faire
x0i Ð xi
FinPour
Jusque (d<ε)
Fin

5.2.3 Convergence des méthodes de Jacobi et de Gauss-Seidel

Les méthodes de Jacobi et de Gauss-Seidel convergent pour les cas suivants :

– A est une matrice positive
– A est une matrice diagonalement dominante, c’est-à-dire
ÿ
|aii | ą |ai j |
j‰i

– A est une matrice tridiagonale par blocs

5.3 Exercices
Exercice 5.1
1. On considère une fonction y de classe C4 sur un intervalle ra, bs.
(a) En utilisant la formule d’interpolation polynomiale de degré 2, établir que y1 “
pypx ` hq ´ ypx ´ hqq{2h où h est le pas de dérivation.
42 Résolution des systèmes algébrique

(b) Déduire de paq une expression pour y2

2. Soit le problème aux valeurs aux limites y2 “ ppxqy1 ` qpxqy ` gpxq défini sur ra, bs avec
ypaq “ c et ypbq “ d. On veut résoudre ce problème par la méthode des différences
finies. Établir à partir de la 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) Écrire alors l’algorithme permettant de résoudre le problème
(c) Traduire l’algorithme en Fortran
Exercice 5.2
On considère un système algébrique linéaire n ˆ n de la forme
Ax “ b
où A “ pai j q, b “ pb j q 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. Établir l’algorithme correspondant
3. De l’algorithme ci-dessus, on obtient un système triangulaire. Établir l’algorithme de
résolution du système triangulaire obtenu.
4. Traduire l’algorithme en Fortran
Exercice 5.3
On considère l’équation intégrale :
żb
φ pxq “ f pxq ` kpx,tqφ ptqdt p1q
a
où kpx,tq est une fonction connue. On subdivise le domaine ra, bs en n points équidistants
ti “ ih
1. L’equation p1q peut être approximée par la formule composite des trapèzes. Elle prend
alors la forme
n
ÿ
φ pxq “ f pxq ` pb ´ aq ci kpx,ti qφi p2q
i“1
où φi “ φ pti q. Donner les expressions des ci
2. Sachant que l’équation p2q est vérifiée en tout point x j “ t j , montrer que l’équation
intégrale p1q se réduit au système algébrique :
n
ÿ
φ j “ f j ` pb ´ aq ci kpt j ,ti qφi p3q
i“1
où les φi et φ j sont les inconnues.
5.3 Exercices 43

3. Le système p3q peut se mettre sous la forme :

Ax “ b p4q
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 p4q par la méthode de Gauss
(a) Donner l’algorithme de triangularisation du système
(b) Écrire l’algorithme de résolution du système triangularisé
(c) Traduire ce dernier algorithme en Fortran
5. φ pxq possède un zéro dans l’intervalle ra, bs et on propose de trouver ce zéro en utilisant
la méthode de bissection. Décrire cette méthode et donner son organigramme
Chapitre 6

Équations aux dérivées partielles

Les équations aux dérivées partielles sont aux fonctions de plusieurs variables ce que les
équations différentielles sont aux fonctions d’une variable réelle. L’objectif de ce chapitre est
d’introduire à l’étude et à la résolution approchée sur ordinateur des équations aux dérivées
partielles les plus simples qu’on rencontre dans les applications. Les équations aux dérivées
partielles sont classées en trois grandes classes fondamentales d’équations : les équations
elliptiques, les équations paraboliques et les équations hyperboliques.
La résolution des équations aux dérivées partielles nécessitent la connaissance des condi-
tions aux limites et des conditions initiales. La nature des conditions permet de définir les
types de problèmes suivants :
– Problème de Dirichlet : les valeurs des points sur les frontières sont connues
– Problème de Newmann : les valeurs des dérivées des points sur les frontières sont
connues
– Problème mixte Dirichlet-Newmann : les valeurs des points sont connues sur une partie
des frontières, et celle des dérivées sur la partie restante
– Problème de Cauchy : les valeurs des points et de ses dérivées sont connues sur les
frontières du domaine

6.1 Équation de Poisson

Considérons l’équation de Poisson suivante en dimension 2 :
B2u B2u
` “ f px, yq (6.1)
Bx2 By2
défini sur le domaine D suivant :
D “ tpx, yq P R2{a ă x ă b, c ă y ă du
sur la frontière S du domaine D, on donne la condition suivante :
upx, yq “ gpx, yq (6.2)
pour discrétiser l’équation p6.1q, en prenant xi “ a ` ih et y j “ c ` jk, on pose
B 2 u ui`1, j ´ 2ui j ` ui´1, j B 2 u ui, j`1 ´ 2ui j ` ui, j´1
“ ` Oph2 q et “ ` Opk2 q
Bx2 h2 By2 k2
45
46 Équations aux dérivées partielles

avec h “ pb ´ aq{n et k “ pd ´ cq{m. Ce qui nous donne

«ˆ ˙ ff
2 ˆ ˙2
h h `
ui, j`1 ` ui, j´1 “ ´h2 fi j
` ˘ ˘
2 ` 1 ui j ´ ui`1, j ` ui´1, j ´ (6.3)
k k
avec "
i “ 1, . . . , n ´ 1
j “ 1, . . . , m ´ 1
L’erreur de troncature pour cette équation est h2 ` k2 . La discrétisation des conditions aux
frontières donne :
" "
u0 j “ gpx0 , y j q ui0 “ gpxi , y0 q
j “ 0, . . . , m i “ 0, . . . , n (6.4)
un j “ gpxn , y j q uim “ gpxi , ym q
En posant n “ m, les relations p6.3q et p6.4q permettent d’obtenir le système algébrique
suivant :
AU “ F ` G (6.5)
”` ˘ ı
2 ` ˘ 2
En posant α “ 2 hk ` 1 et β “ hk , on a :
» fi » fi
D ´I 0 ¨ ¨ ¨ ¨ ¨ ¨ 0 α ´β 0 ¨¨¨ ¨¨¨ 0
. . ... ..
— ´I D ´I . . .. ffi .
— ffi — ffi
— ´β α ´β ffi
. . .
— 0 . . . . . . . . .. ffi
— . . ffi
— 0 ... ...
— ... ... .. ffi
. ffi
A“— . .
—
. . .
ffi avec D “ —
— .. ... ... ... ...
ffi
— .. . . . . . . . . 0 ffi ffi
— . 0
ffi
ffi
— .. ... ... ... ffi — .. ... ... ... ffi
– . ´I fl – . ´β fl
0 ¨ ¨ ¨ ¨ ¨ ¨ 0 ´I D 0 ¨¨¨ ¨¨¨ 0 ´β α
et
» fi
0 0 0 ´h2 f11
» fi » fi
´1 ¨¨¨ ¨¨¨ u1n
... .. .. ..
0 ´1 0 .
— ffi
— ffi — . ffi — . ffi
. . . ...
— ffi — ffi
—
0 ... ... ... ffi — u ffi — 2
´h f1,n´1 ffi
— ffi — 1,n´1
U “— F “—
ffi — ffi
´I “ — .. ... ... ... ...
ffi .. ffi .. ffi
—
— . 0
ffi
ffi — . ffi — . ffi
— .. ... ... ... ffi — .. ffi — .. ffi
– . 0 fl – . fl – . fl
un´1,n´1 2
´h fn´1,n´1
0 ¨¨¨ ¨¨¨ 0 0 ´1
G est un vecteur creux comportant les valeurs de u aux frontières du domaine. Dans la
pratique, on peut utiliser la méthode de Jacobi ou de Gauss-Seidel pour résoudre ce système.
La matrice A étant tridiagonale par bloc, profitons de cette structure particulière et déduisons
de p6.3q la relation suivante :
1
p´h2 fi j ` pui`1, j ` ui´1, j q ´ β pui, j`1 ` ui, j´1 qq
ui j “ (6.6)
α
Le schéma discret de Jacobi pour l’équation de poisson est donnée par :
1
uk`1
ij “ p´h2 fi j ` puki`1, j ` uki´1, j q ´ β puki, j`1 ` uki, j´1 qq (6.7)
α
6.1 Équation de Poisson 47

Le schéma discret de Gauss-Seidel pour l’équation de poisson est donnée par :

1
uk`1
ij “ p´h2 fi j ` puki`1, j ` uk`1 k k`1
i´1, j q ´ β pui, j`1 ` ui, j´1 qq (6.8)
α

Algorithme Poisson-Gauss-Seidel
Début
{ on definit les conditions aux limites }
Pour i allant de 0 à n Faire
ui0 Ð gpxi , y0 q
uim Ð gpxi , ym q
FinPour
Pour j allant de 0 à m Faire
u0 j Ð gpx0 , x j q
un j Ð gpxn , y j q
FinPour
{ on initialise ui j }
Pour i allant de 1 à n Faire
Pour j allant de 1 à m Faire
ui j Ð 0
FinPour
FinPour
Répéter
S1 Ð 0 ; S2 Ð 0
Pour i allant de 1 à n-1 Faire
Pour j allant de 1 à m-1 Faire
v Ð α1 p´h2 fi j ` ui`1, j ` ui´1, j ´ β pui, j`1 ` ui, j´1 qq
S1 Ð S1+|v ´ ui j |
S2 Ð S2+|v|
ui j Ð v
FinPour
FinPour
d Ð S1S2
Jusque (d ă ε)
Fin

Remarque 6.1
Lorsque les conditions aux limites sont de types Neumann, le schéma de discrétisation
doit être fait avec beaucoup de précaution. Généralement il est conseillé d’associer cette
condition de Neumann avec l’équation aux dérivées partielles de départ. Par exemple si la
condition de Neumann s’écrit ˆ ˙
Bu
“ gpx, yq (6.9)
By S
48 Équations aux dérivées partielles

alors on utilise le schéma suivant : on développe upxi , y1 q autour de upxi , y0 q pour obtenir :
Bu k2 B 2 u
upxi , y1 q “ upxi , y0 q ` k pxi , y0 q ` 2
pxi , y0 q ` Opk3 q (6.10)
By 2 By
B2u
La quantité 2 pxi , y0 q peut être remplacée par son équivalent à partir de l’EDP de départ.
By
Par exemple, dans le cas de l’équation de Laplace, on a :
B2u B2u
“´ 2
By2 Bx
Le développement précédent donne alors :
Bu k2 B 2 u
upxi , y1 q “ upxi , y0 q ` k pxi , y0 q ´ pxi , y0 q
By 2 Bx2
k2 ui`1,0 ´ 2ui,0 ` ui´1,0
ˆ ˙
“ ui,0 ` kgpxi , y0 q ´
2 h2
Finalement la discrétisation de la condition de Neumann permet d’obtenir l’équation discrète
suivante :
k2 k2
ˆ ˙
ui,1 ´ 1 ` 2 ui,0 ` 2 pui`1,0 ` ui´1,0 q “ kgpxi , y0 q (6.11)
h 2h
On peut alors associer p6.9q à l’équation p6.3q pour avoir le système discret global à ré-
soudre.

6.2 Équation de la chaleur

On cherche à résoudre numériquement l’équation de la chaleur donnée par :

R` ,
2
Bu 2B u
“ c 2, t P x Ps0, lr (6.12)
Bt Bx
les conditions aux limites sont données par :
up0,tq “ upl,tq “ 0 (6.13)
et la condition initiale est donnée par :
upx, 0q “ f pxq (6.14)

6.2.1 Méthode de différences progressives ou schéma explicite

Pour cette méthode, on utilise les discrétisations suivantes :

Bu ui, j`1 ´ ui, j B 2 u ui`1, j ´ 2ui, j ` ui´1, j
“ ` Opkq 2
“ 2
` Oph2 q
Bt k Bx h
6.2 Équation de la chaleur 49

En remplaçant ces relations dans p6.10q, on obtient :

ui, j`1 ´ ui, j ui`1, j ´ 2ui, j ` ui´1
“ c2 (6.15)
k h2
Ce qui donne
2c2 k c2 k
ˆ ˙
ui, j`1 “ 1 ´ 2 ui, j ` 2 pui`1, j ` ui´1, j q (6.16)
h h
avec "
i “ 1, . . . , n ´ 1
j “ 1, . . . , 8
Le schéma p6.14q présente une erreur de troncature d’ordre 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
leur valeur connue et stockée en mémoire de la machine.
Les conditions initiales et aux limites se discrétisent de la manière suivante :
ui,0 “ f pxi q i “ 0, 1, . . . , n
"
u0, j “ 0 (6.17)
j “ 0, . . . , 8
un, j “ 0
La première condition de p6.15q est substituée dans p6.14q pour déterminer les valeurs
ui,1 pour i “ 1, . . . , n ´ 1. La condition aux limites u0, j “ 0 et um, j “ 0 permet d’obtenir
u0,1 “ um,1 “ 0. On a donc toutes les valeurs ui,1 . En procédant de la même manière on ob-
tient toutes les valeurs ui,2 , ui,3 ¨ ¨ ¨
L’inconvénient de cette méthode explicite est qu’elle est conditionnellement stable, c’est-à-
dire qu’une condition doit exister sur h et k enfin que les résultats obtenus soit bon. Pour
établir cette condition de stabilité, nous réécrivons le schéma discrétisé sous la forme sui-
vante :
ui, j`1 ´ ui, j ui`1, j ´ 2ui, j ` ui´1, j
“ c2
k h2
Supposons qu’a cause des erreurs de troncature et d’arrondi, la valeur de ui, j subit une per-
turbation εi, j . La nouvelle solution numérique devient alors u1i, j “ ui, j ` εi, j . Elle satisfait
aussi à l’équation discrète. En écrivant que u1i, j vérifie le schéma discret, on obtient que εi, j
vérifie l’équation discrète suivante :
εi, j`1 ´ εi, j εi`1, j ´ 2εi, j ` εi´1, j
“ c2 (6.18)
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 (lorsque j croit). Supposons que εi, j est de la forme suivante :
εi, j “ sinpiphqe´r jk (6.19)
où p et r peuvent être des nombres complexes ou réels.
En introduisant cette expression dans l’équation vérifiant εi, j , on obtient

´rk 4kc2 2 ph
e “ 1 ´ 2 sin (6.20)
h 2
50 Équations aux dérivées partielles

Donc e´rk est le facteur d’amplification. Le schéma est stable si le facteur d’amplification
est inférieur à 1. Dans le cas particulier analysé ici, cette condition implique que :
2
ˇ ˇ
ˇ
ˇ1 ´ 4kc 2 ph ˇ
2
sin ˇă1
ˇ h 2ˇ
2ˇ kc2 1
ˇ ˇ
ˇ 4kc
ùñ ˇˇ1 ´ 2 ˇˇ ă 1 ùñ 2 ă
h h 2

6.2.2 Méthode des différences finies régressives ou méthodes implicites

Enfin d’éviter les problèmes de stabilité que l’on rencontre dans la méthode explicite, on
peut utiliser les différences finies régressives pour la dérivée temporelle, c’est-à-dire
Bu ui, j ´ ui, j´1
“ (6.21)
Bt k
Dans ce cas on obtient le schéma discret suivant :
c2 k
p1 ` λ qui, j ´ λ ui`1, j ´ λ ui´1, j “ ui, j´1 avec λ “ 2 (6.22)
h
En prenant en compte les conditions initiales et aux limites, on obtient un système matriciel
tridiagonal qu’on peut résoudre par les méthodes itératives ou la méthode d’élimination de
Gauss. Le facteur d’amplification de ce schéma est donné par :
1
e´rk “ (6.23)
4kc2
1` h2
sin2 ph
2

Il apparait que le facteur d’amplification dans ce schéma est toujours inférieur à 1 ; c’est-à-
dire que le schéma est inconditionnellement stable

6.2.3 Méthode de Cranck-Nichelson

Cette méthode fait une sommation des méthodes des différences finies progressives et
régressives après avoir écrit la formule de différences régressives à l’instant j ` 1. On obtient
alors le schéma suivant :
ui, j`1 ´ ui, j c2
´ 2 pui`1, j ´ 2ui, j ` ui´1, j ` ui`1, j`1 ´ 2ui, j`1 ` ui´1, j`1 q “ 0 (6.24)
k 2h
Le facteur d’amplification de ce schéma est donné par :
2
1 ´ 2ch2k sin2 ph
2
e´rk “ 2 (6.25)
1 ` 2ch2k sin2 ph
2

On voit bien qu’on a toujours |e´rk | ă 1, donc la méthode de Cranck-Nichelson est incondi-
tionnellement stable.
6.3 L’équation d’onde 51

6.2.4 Méthode de Richardson

Pour la méthode de Richardson, on fait en sorte que l’erreur de troncature pour la dérivée
temporelle soit d’ordre Opk2 q. Pour cela on utilise le schéma discret suivant :
Bu ui, j`1 ´ ui, j´1
“ ` Opk2 q (6.26)
Bt 2k
Le schéma de Richardson prend alors la forme suivante :
ui, j`1 ´ ui, j´1 ui`1, j ´ 2ui, j ` ui´1, j
´ c2 “0 (6.27)
2k h2
C’est un schéma conditionnellement stable.

6.2.5 Méthode de Dufort-Frankel

Pour résoudre le problème de d’instabilité de la méthode de Richardson, Dufort et Frankel

ont proposé de remplacer la quantité ui, j par
ui, j`1 ` ui, j´1
(6.28)
2
Le schéma discret se met alors sous la forme suivante :
ui, j`1 ´ ui, j´1 ui`1, j ´ ui, j`1 ´ ui, j´1 ` ui´1, j
´ c2 (6.29)
2k h2

6.3 L’équation d’onde

Nous considérons le problème suivant
B2u 2
2B u
´ c “0 sur D “ T ˆ L (6.30)
Bt 2 Bx2
et soumis aux conditions suivantes
up0,tq “ upl,tq “ 0
Bu (6.31)
upx, 0q “ f pxq px, 0q “ gpxq
Bt
L’utilisation de la méthode de différences finies permet de mettre l’EDP sous la forme dis-
crète suivante :
ui, j`1 “ 2p1 ´ λ 2 qui, j ` λ 2 pui`1, j ` ui´1, j q ´ ui, j´1 (6.32)
avec "
i “ 1, . . . , n ´ 1 ck l
; λ“ et n “
j “ 1, . . . , 8 h h
La discrétisation des conditions aux limites donne :
u0, j “ un, j “ 0 (6.33)
52 Équations aux dérivées partielles

La première condition initial donne :

ui,0 “ f pxi q (6.34)
La deuxième condition initiale permet d’écrire :
ui,1 “ ui,0 ` kgpxi q (6.35)
Le schéma discret p6.33q présente un erreur de troncature d’ordre Opkq. Pour avoir une erreur
de troncature d’ordre Opk2 q comme dans p6.30q, on procède de la manière suivante :
Bupxi , 0q upxi ,t1 q ´ upxi , 0q k B 2 u
“ ´ 2
pxi , 0q ` Opk2 q
Bt k 2 Bt
Or
B2u 2
2B u 2 p f i`1 ´ 2 f i ` f i´1 q
px i , 0q “ c px i , 0q “ c
Bt 2 Bx2 h2
on obtient donc que la deuxième condition initial donne le schéma discret suivant :
λ2
ui,1 “ p1 ´ λ 2 q fi ` p fi`1 ` fi´1 q ` kgi (6.36)
2
La résolution numérique de l’équation d’onde revient à résoudre numériquement le système
discret formé des équations p6.30q et p6.34q. La recherche des valeurs ui, j`1 nécessite la
connaissance des valeurs ui, j et ui, j`1 . Au début, il faut connaitre les ui,0 et ui,1 . Les valeurs
de ui,0 découlent de p6.32q et celle de ui,1 découlent de p6.34q. On peut ainsi lancer le pro-
cessus itératif pour calculer toutes les valeurs de ui, j . Il faut cependant noter que ce schéma
ck
de discrétisation a une condition de stabilité qui est ă1
h
Remarque 6.2
On peut établir d’autres schémas discrets en décomposant l’équation d’onde en 2 équa-
tions du 1er ordre avant de leur discrétiser. Dans cette perspective, on a plusieurs types de
schémas : schéma explicite, schéma implicite, schéma de Lax, et le schéma de Lax-Wendroff

6.4 Équation d’advection-diffusion

L’equation d’advection-diffusion est de la forme :
BΩ BΩ BΩ
` u1 ` u2 “ v∆Ω (6.37)
Bt Bx By
où u1 et u2 sont des fonctions connues qui dépendent de x, y et t.
Cette équation décrit par exemple la distribution des particules ou des gaz émis par une
source plane en présence du vent. Dans ce cas Ω est la concentration des particules ou
des gaz. u1 et u2 sont les composantes de la vitesse du vent. Il s’agit donc d’une équation
intéressante qui permet d’expliquer et de déterminer la concentration des polluants dans
l’air ; lesquels sont émis par des industries. Si v “ 0,l’équation se réduit à celle de transport.
Si u1 “ u2 “ 0, on a l’équation de la chaleur. Pour discrétiser cette équation nous prenons
∆t “ k, ∆x “ ∆y “ h
6.5 Exercices 53

6.4.1 Schéma de discrétisation centré en espace

Pour ce schéma, les dérivées spatiales s’expriment par les relations

n n n n
BΩ Ωi`1, j ´ Ωi´1, j BΩ Ωi, j`1 ´ Ωi, j´1
“ “
Bx 2h By 2h
La discrétisation dans le temps se fait par la relation :
n`1 n
BΩ Ωi, j ´ Ωi, j
“
Bt k
Définitivement, l’équation d’advection-diffusion se met sous la forme discrète suivante :
ku1 ` n ˘ ku2 ` n ˘
Ωn`1 n
i, j “ Ωi, j ´ Ωi`1, j ´ Ωni´1, j ´ Ωi, j`1 ´ Ωni, j´1
2h 2h (6.38)
vk ` ˘
` 2 Ωni`1, j ` Ωni´1, j ` Ωni, j`1 ` Ωni, j´1 ´ 4Ωni, j
h
On établit que la condition de stabilité du schéma est :
vk 1
ď (6.39)
h2 4

6.4.2 Schéma explicite décentré en espace

Enfin d’obtenir des résultats plus précis, l’expérience des simulations numériques a per-
mis d’établir que les signes des grandeurs u1 et u2 peuvent avoir une influence sur la qualité
des résultats. C’est ainsi qu’est né le schéma explicite décentré en espace qui a la forme
suivante :
ku1 n ku2 n
Ωn`1 n
i, j “ Ωi, j ´ P ´ Q
h i, j h i, j (6.40)
vk “ ‰
` 2 Ωni`1, j ` Ωni´1, j ` Ωni, j ` Ωni, j´1 ´ 4Ωni, j
h
où
Ωni, j ´ Ωni´1, j si
"
u1 ą 0
Pi,nj “ (6.41)
Ωni`1, j ´ Ωni, j si u1 ă 0

Ωni, j ´ Ωni, j`1 si

"
u2 ą 0
Qni, j “ (6.42)
Ωni, j`1 ´ Ωni, j si u2 ă 0

6.5 Exercices
Exercice 6.1
54 Équations aux dérivées partielles

A l’instant t “ 0, la température d’un mur d’épaisseur L vaut T0 . Sur la surface x “ 0, la

température varie suivant la loi :
π
T px “ 0,tq “ T0 ` T1 sin t t ą 0
2
Par un mécanisme approprié, la face x “ L est maintenue à la température constante TL .
L’equation de la chaleur dans le milieu est alors :
BT B2T
“α 2
Bt Bx
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) Écrire l’algorithme de résolution du système discret obtenu
(c) Traduire cet algorithme 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ératures Ti, j “ T px “ i∆x, ∆tq, établit l’algorithme permettant de trouver la
température aux instants ultérieurs.
(c) Traduire cet algorithme en Fortran
Exercice 6.2
1. On considère l’équation de la chaleur
2
Bu 2B u
“a 2
Bt Bx
R
définie sur D “ ` ˆs0, lr avec up0,tq “ upl,tq “ 0 et upx, 0q “ f pxq. On veut résoudre
cette équation par le schéma explicite avec
Bu ui, j`1 ´ ui, j
“
Bt k
Établir que la condition de stabilité de ce schéma est
a2 k 1
ď où k et h sont les pas temporel et spatial
h2 2
2. Que devient cette condition de stabilité dans le cas d’un problème plan :
Bu
“ a2 ∆u
Bt
On prendra h1 et h2 les pas spatiaux.
6.5 Exercices 55

3. On considère l’équation d’advection-diffusion dans le plan :

BΩ BΩ BΩ
` u1 ` u2 “ v∆Ω
Bt Bx By
On la discrétise par le schéma centré en espace. Établir que la condition de stabilité est
vk
h2
ď 14 où k est le pas temporel et h est le pas spatial suivant x et y. Que devient cette
condition dans le cas unidimensionnel ?

Exercice 6.3
Soit l’équation d’onde suivante :
B2u Bu 2
2B u
2
` a “ c 2 ` du3
Bt Bt Bx
R`Ys0, lr ; up0,tq “ upl,tq “ 0 ; upx, 0q “ f pxq ; BuBt px, 0q “ gpxq. a, c et d sont des
où pt, xq P
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é.
3. Exprimer ui,2 “ upx “ ih,t “ 2kq 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. Établir un algorithme de résolution du système discret obtenu à la question (2).
5. Traduire cet algorithme en Fortran

Exercice 6.4
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 surface extérieure est
nulle. A chaque point de la région hachurée, la distribution des pressions vérifie l’équation :
B2 p B2 p
` “0
Bx2 By2
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 pxi “ ih, y j “ jkq
(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) Établir l’algorithme de calcul des ui, j
56 Équations aux dérivées partielles

4. On suppose que le tube a maintenant une forme cylindrique de rayon intérieur R1 “

400 mm et de rayon extérieur R2 “ 800 mm
(a) Quelle est alors l’équation vérifiée par la pression p ou la grandeur grandeur adi-
mensionnée u.
(b) Proposer une méthode de discrétisation de cette nouvelle équation
On donne : L “ 1000 mm, L1 “ 500 mm, l “ 800 mm, l1 “ 400 mm
Exercice 6.5
La répartition d’une grandeur physique u sur la plaque représentée vérifie l’équation de
poisson
B2u B2u
` “ f px, yq
Bx2 By2
avec u “ gpx, yq sur la frontière.
1. En utilisant la méthode des différences finies, établit l’équation discrète vérifiée par ui, j
2. On se propose de trouver les ui, j par la méthode itérative de Jacobi.
(a) Décrire cette méthode
(b) Établir l’algorithme de calcul de ui, j
Exercice 6.6
La répartition d’une grandeur physique u sur la plaque triangulaire ci-dessus vérifie l’équa-
tion de Laplace :
B2u B2u
` “0
Bx2 By2
Bu
avec u “ 0 sur l’hypoténuse et Bn “ 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. Établir l’algorithme de calcul des valeurs discrètes ui, j sur toute la plaque
 II
Correction des Exercices
Exercices du chapitre 1

Exercice 1.1

59
Exercices du chapitre 2

61
Exercices du chapitre 3

Exercice 3.1
La répartition de l’énergie rayonnée par un étoile en fonction de la longueur d’onde est
donnée par la loi de Planck :
8πhc 1
Ipλ q “ 5 ˆ
λ expp λhc
kT q
où h, k et c sont respectivement les constantes de Planck, de Boltzmann et la célérité de la
lumière.
Posons x “ λhc kT , on a : ˆ ˙4
kT x5
Ipxq “ 8πkT ˆ
hc exppxq ´ 1
ˆ ˙4 „ 4 x
5x pe ´ 1q ´ x5 ex

dIpxq kT
ùñ “ 8πkT
dx hc pex ´ 1q2
1. Établissons l’équation vérifiée par xm
dIpxq
“ 0 ðñ xm exm ´ 5exm ` 5 “ 0
dx
d’où l’équation vérifiée par xm
2. Établissons l’algorithme permettant de calculer la valeur numérique de xm
Il faut résoudre l’équation suivante :
f pxq “ 0 avec f pxq “ xex ´ 5ex ` 5

Algorithme Planck
{ définit la fonction f }
{ définit la derivée d f }
Début
eps Ð ? ; s Ð ?
Répéter
pÐs
s Ð p- dffppq
ppq
Jusque (|s-p|<eps)
Fin

63
64 Exercices du chapitre 3

Exercice 3.2

L’algorithme de Van der Waals pour une mole de gaz de masse M est :
´ a¯
p ` 2 pv ´ bq “ RMT
v
Etablissons un algorithme qui permet de calculer le volume v chaque fois que l’on connait
la temperature et la pression.
On a :
´ a¯
f pvq “ p ` 2 pv ´ bq ´ RMT “ 0
v
On utilise pour cela l’algorithme de Newton-Raphson

Algorithme Van-der-Waals
{ définit la fonction f }
{ définit la dérivée d f }
Début
eps Ð ? ; v Ð ?
Répéter
sÐv
v Ð s- dffpsq
psq
Jusque (|v-s|<eps)
Fin

Exercice 3.3

En 1976, Desantis a établi que le facteur de compressibilité b des gaz réels vérifie la
relation :
z “ p1 ` y ` y2 ´ y3 q{p1 ´ y3 q

où y “ b{4v, v est le volume molaire et z est une constante

1. Voir cours

2. Établissons un algorithme permettant de résoudre le problème.

On a :
f pyq “ pz ´ 1qy3 ` y2 ` y ` 1 ´ z “ 0
Exercices du chapitre 3 65

Algorithme Desantis
{ definit la fonction f }
{ definit la derivee d f }
Début
eps Ð ? ; y Ð ?
Répéter
sÐy
y Ð s- dffpsq
psq
Jusque (|y-s|<eps)
b Ð 4vy
Fin

3. Traduisons cet algorithme en Fortran

program desantis
implicit none
real :: y ,b ,s , eps
eps =? ; y =?
!
do
s=y
y =s -( f ( s )/ df ( s ))
if ( abs (y - s ) < eps ) exit
end do
b = y *4* v
contains
real function f ( y )
implicit none
real :: y , z
f =(1 - z )* y **3+ y **2+ y +1 - z
end function
real function df ( y )
implicit none
real :: y , z
df =3*(1 - z )* y **2+2* y +1
end function
!
end program

Exercice 3.4
1. On est amené à rechercher les racines de la fonction :
f pxq “ 1 ´ cos x cosh x
66 Exercices du chapitre 3

Cette fonction varie fortement, car cosh x prend des valeurs très grandes tandis que cos x
oscille du positif au négatif. Pour simplifier le traitement numérique, un bonne stratégie
consiste à considérer la fonction :
1 ´ cos x cosh x 1
f1 pxq “ “ ´ cos x
cosh x cosh x
qui possède les mêmes racines que f pxq. On constate à l’aide d’une étude graphique
qu’il y’a changement de signe dans les intervalles r3, 5s, r6, 8s et r10, 12s. La méthode
de bissection appliquée à chacun de ces trois intervalles donne x1 “ 4, 730, x2 “ 7, 853
et x3 “ 10, 996
program bissection
implicit none
integer :: k
real :: eps ,a ,b , c
real , dimension (1:6):: x
data x /3.0 ,5.0 ,6.0 ,8.0 ,10.0 ,12.0/
eps =1.0 e -6
do k =1 ,6 ,2
a = x ( k ) ; b = x ( k +1)
do while ( abs (a - b ) > eps )
c =( a + b )/2.0
if ( f ( c )* f ( b ) <0.0) then
a=c
else
b=c
end if
end do
print * , c
end do
contains
real function f ( x )
implicit none
real :: x
f =(1 - cos ( x )* cosh ( x ))/ cosh ( x )
end function
end program

2. On obtient donc finalement :

4, 730 7, 853 10, 996
K1 “ K2 “ K3 “
L L L
Ce qui correspond donc aux valeurs de la littérature
Exercices du chapitre 4

Exercice 4.1
Le mouvement d’un pendule est décrit par l’équation différentielle
d 2Y
` ω0 sinY “ 0
dt 2
où ω0 est la pulsation propre et Y le déplacement angulaire
1. (a) On a :
d 2Y
“ ´ω0 sinY “ f pY q
dt 2
Posons
dY
“v
dt
On a donc : " dY
dt “v
dv
dt “ ´ω0 sinY “ f pY q
L’algorithme de résolution est le suivant :
Algorithme RK2
{ On definit la fonction f }
Début
y0 Ð ? ; v0 Ð ? ; h Ð ?
Pour k allant de 0 à n-1 Faire
l1 Ð hf(yk )
c1 Ð hvk
l2 Ð hf(yk +c1)
c2 Ð h(vk +l1)
vk`1 Ð vk + 21 (l1+l2)
yk`1 Ð yk + 21 (c1+c2)
FinPour
Fin

(b) Traduisons cet algorithme en Fortran

program rk2
implicit none
67
68 Exercices du chapitre 4

integer , parameter :: n =?
integer :: k
real :: h , l1 , l2 , c1 , c2
real , dimension (0: n ):: y , v
y (0)=? ; v (0)=? ; h =?
do k =0 ,n -1
l1 = h * f ( y ( k ))
c1 = h * v ( k )
l2 = h * f ( y ( k )+ c1 )
c2 = h *( v ( k )+ l1 )
v ( k +1)= v ( k )+( l1 + l2 )/2.0
y ( k +1)= y ( k )+( c1 + c2 )/2.0
end do
contains
real function f ( y )
implicit none
real :: y
f = - omega * sin ( y )
end function
end program

2. Lorsque le pendule oscille entre Y0 et ´Y0 , sa période est définie par la relation
2T0 π{2
ż
dx
T pY0 q “ a
π 0 1 ´ Apsin xq2
où T0 “ 2π{ω0 et A “ sin2 pY0{2q
(a) Pour l’établissement de la formule composite des trapèzes, voir exercice 2.2 du
chapitre 2.
On a : ż π{2
2T0 1
T pY0 q “ f pxqdx avec f pxq “ a
0 π 1 ´ Apsin xq2
On a donc :
ż π{2 n´1
ÿh π{2
f pxqdx » p f pxi q ` f pxi`1 qq avec h “
0 i“0
2 n

L’algorithme est le suivant :