Mathématiques générales ISA4

Mourad ABOUZAID

ISABTP 4e année
2019-2020
version enseignant

Version du 28 août 2019

2
Table des matières

1 Systèmes différentiels 5
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.1 Définitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.2 Courbes intégrales et portrait de phase . . . . . . . . . . . . . . . . 8
1.1.3 Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2 Systèmes différentiels linéaires . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.1 Définitions et représentation matricielle . . . . . . . . . . . . . . . . 9
1.2.2 Structure de l’ensemble des solutions . . . . . . . . . . . . . . . . . 10
1.2.3 Matrice résolvante . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.4 Système linéaires à coefficients constants . . . . . . . . . . . . . . . 13
1.2.5 Classification des portraits de phases (systèmes linéaires 2 × 2) . . . 15
1.3 Étude qualitative des systèmes autonomes . . . . . . . . . . . . . . . . . . 17
1.3.1 Systèmes autonomes et champs de vecteurs . . . . . . . . . . . . . . 17
1.3.2 Points d’équilibre est stabilité . . . . . . . . . . . . . . . . . . . . . 18
1.3.3 Linéarisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.3.4 Cycles limites et attracteurs étranges . . . . . . . . . . . . . . . . . 20

2 Résolution numériques d’équations différentielles 21

2.1 Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2 Différences finies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Application à la résolution des EDO . . . . . . . . . . . . . . . . . . . . . 25
2.3.1 Schémas d’Euler (explicite) . . . . . . . . . . . . . . . . . . . . . . 25
2.3.2 Schémas numériques . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.3.3 Consistance, stabilité, convergence . . . . . . . . . . . . . . . . . . 28

3 Résolution numériques d’équations aux dérivées partielles 33

3.1 Notions de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.1.1 Définitions et notations . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.1.2 Exemples d’EDP classiques . . . . . . . . . . . . . . . . . . . . . . 34
3.1.3 Classification des EDP linéaires d’ordre 2 . . . . . . . . . . . . . . . 35
3.1.4 Conditions initiales et conditions limites . . . . . . . . . . . . . . . 36
3.2 Différences finie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3 Illustrations de la méthode des différences finies . . . . . . . . . . . . . . . 41
3.3.1 Équation de la chaleur 2D, stationnaire . . . . . . . . . . . . . . . . 41
3.3.2 Équation de la chaleur 1D instationnaire . . . . . . . . . . . . . . . 52

3
4 TABLE DES MATIÈRES

3.4 Résolution exacte par séparation des variables . . . . . . . . . . . . . . . . 58

Chapitre 1

Systèmes différentiels

Introduction
Les équations différentielles et les équations aux dérivées partielles sont au cœur de la
modélisation des phénomènes physiques.

L’étude théorique et pratique de ces équations et des systèmes associés a pour cela
rapidement constitué une part importantes des mathématiques modernes. Ces études on
en particulier produit de nombreuses méthodes d’études qualitatives et quantitatives pour
ce type de problèmes, ainsi que des méthodes de résolution approchées, spécifiquement
développées pour pouvoir être confiées à un ordinateur.

Dans les années précédentes, nous avons étudié différentes méthodes d’études pour les
équations différentielles (méthodes avant tout basées sur une classification rigoureuse de
ce type d’équations).

Le but de ce chapitre est de présenter les résultats spécifiques aux systèmes différen-
tiel, i.e. aux systèmes de plusieurs équations différentielles, portant sur un ensemble de
fonctions d’une même variable.

Nous verrons en particulier que de nombreux résultats théoriques vus pour les équa-
tions seules restent valables pour les systèmes différentiels.

Précisément, après quelques généralités, nous présenterons des méthodes de résolution

algébrique valable pour un certain type de systèmes : les systèmes différentiels linéaires à
coefficients constants.

Dans un second temps, nous présenterons des méthodes d’étude qualitative de sys-
tèmes plus généraux : les systèmes autonomes.

Les méthodes numériques de résolution approchées évoquées ci-dessus seront abordées

dans les chapitres suivants de cette année. Les méthodes d’étude qui suivent nous permet-
trons alors de comparer aux résultats attendus les résultats approchés obtenus à l’aide de
méthodes algorithmiques.

5
6 CHAPITRE 1. SYSTÈMES DIFFÉRENTIELS

1.1 Généralités
1.1.1 Définitions
Formellement, une équations différentielle est une équation liant une fonction d’une
variable réelle à ses dérivées. D’autre part, on appelle système différentiel la donnée d’un
ensemble d’équations différentielles, liant entre elles différentes fonctions d’une même va-
riable réelle et leurs dérivées. Précisément,

Définition (équation différentielle - système différentiel)

Soit n ∈ N∗ .
1. On appelle équation différentielle d’ordre n, portant sur une fonction y d’une va-
riable t, toute équations de la forme

(E) : y (n) = ϕ(y, y 0 , . . . , y (n−1) , t)

où ϕ est une fonction connue.

2. On appelle système différentiel d’ordre n, portant sur un ensemble y1 , . . . , yp de

fonctions d’une même variable t, un ensemble d’équations de la forme
 (n) (n−1) (n−1)
 y1 = ϕ1 (y1 , . . . , y1
 , . . . , y p , . . . , yp , t)
 y (n) = ϕ (y , . . . , y (n−1) , . . . , y , . . . , y (n−1) , t)

2 2 1 1 p p
(S) : .

 .
.

 (n) (n−1) (n−1)
yp = ϕm (y1 , . . . , y1 , . . . , yp , . . . , y p , t)

Les fonctions inconnues peuvent représenter différentes quantités évoluant conjoin-

tement au sein d’un même système dynamiques (coordonnées, température, pression,
position, etc) en fonctions, par exemple du temps. Résoudre un tel système consiste à dé-
terminer explicitement chacune des fonctions inconnues. Cela permet alors de connaître
l’état du système étudié à chaque instant.

L’objet de ce chapitre est donc de présenter des méthodes mathématiques permettant

de résoudre ces systèmes différentiels.

Dans le cas d’un système linéaire, il s’agit précisément de déterminer l’ensemble des
p-uplets de fonctions (y1 , . . . , yp ) qui vérifient l’ensemble des équations différentielles du
système.

Comme dans le cas d’une équations seule, la première étape permettant la mise en
place de méthodes d’étude consiste en une classification des différents types de systèmes
différentiels.
1.1. GÉNÉRALITÉS 7

Ainsi, comme dans le cas d’une seule équation, on peut distinguer les systèmes diffé-
rentiels en fonction de leur ordre, i.e. la plus haute dérivée intervenant dans ses équations.
Cependant, la proposition ci-dessous permet de transformer n’importe quel système diffé-
rentiel en un système d’ordre 1, quitte à augmenter le nombre d’équations et d’inconnues :

Proposition

Soit (E) une équation différentielle d’ordre n :

(E) : ϕ(y, y 0 , . . . , y (n) , t) = 0

En posant 

 y1 = y
 y2 = y 0

..

 .
 y = y (n−1)

n

toute solution de (E) peut être associée à une et une seule solution du système

0
 y1
 = y2
0
 y2 = y3



(S) : ..
.
y0

 = yn
 yn−1


0

n = ϕ(y1 , y2 , . . . , yn , t)

Ce résultat s’étendant aux équations d’ordre n > 1 issues d’un système, nous pourrons
nous contenter, dans la suite, d’étudier les systèmes différentiels d’ordre 1 :


 y10 = ϕ1 (y1 , . . . , yp , t)
 y 0 = ϕ2 (y1 , . . . , yp , t)

2
(S) : ..

 .
 y 0 = ϕ (y , . . . , y , t)

p p 1 p

Notons en particulier qu’un tel système est carré.

Par ailleurs, en note Y la fonction vectorielle

Y : R −→ Rp
t 7−→ (y1 (t), . . . , yp (t)
on peut construire une forme vectorielle pour le système (S) ci-dessus sous la forme

(S) ⇐⇒ Y 0 = F (Y, t)
8 CHAPITRE 1. SYSTÈMES DIFFÉRENTIELS

basée sur une fonction F : Rp+1 → Rp définie à partir des fonctions ϕ1 à ϕp .

1.1.2 Courbes intégrales et portrait de phase

En petites dimensions (i.e. si p = 2 ou p = 3), on peut définir une représentation
géométrique de l’ensemble des solutions d’un tel système. En effet, dans le cas p = 2, on
a
 0
x = ϕ(x, y, t)
(S) :
y 0 = ψ(x, y, t)

Toute solution de (S) est alors un couple de fonctions (x, y) dépendant d’une même
variable t et définie sur un intervalle I ⊂ R. En se plaçant dans le plan muni d’un repère
(orthonormé), on peut alors représenter une telle solution par l’ensemble

T(x,y) = {(x(t), y(t)), t ∈ I}

Les fonctions x et y étant continues, cet ensemble est une courbe du plan appelée
courbe intégrale du système (S) ou trajectoire de (S).

Enfin, dans ce cas, l’ensemble des solutions de (S) peut être représenté par une famille
de courbes formant le portrait de phase du système (S).

Note : dans le cas p = 3, chaque solution de (S) peut être associée à une courbe
de l’espace muni d’un repère. L’ensemble de ces courbes forme, là encore, le portrait de
phase de (S). Enfin, bien que l’on n’ait plus accès à l’outil géométrique, il est formelle-
ment possible de définir les courbes intégrales et le portrait de phase pour des systèmes
différentiels de toutes tailles.

1.1.3 Problème de Cauchy

Comme dans le cas p = 1, on peut montrer que, sauf cas particulier, tout système
différentiel de p équations d’ordre 1 admet une infinité de solutions. On peut alors isoler
une unique solution en ajoutant au système un ensemble de p conditions initiales

y1 (t0 ) = y10 , y2 (t0 ) = y20 , . . . , yp (t0 ) = yp0

en une valeur t0 fixée du paramètre.

Note : via l’écriture vectoriel du système (S), l’ensemble des conditions initiales
peuvent se donner sous la forme Y (t0 ) = Y0 ∈ Rp fixé.

On obtient ainsi un problème de Cauchy et l’on peut énoncer le théorème suivant :

1.2. SYSTÈMES DIFFÉRENTIELS LINÉAIRES 9

Théorème

Soit (P ) un problème de Cauchy de la forme

 0
Y = F (Y, t) (S)
(P ) :
Y (t0 ) = Y0 (C.I.)

Si la fonction F vérifie certaines hypothèses de continuité (que l’on supposera toujours

vraies), alors le problème (P ) admet une unique solution.

D’un point de vue géométrique, cela signifie qu’en chaque point Y0 ∈ Rp passe une
et une seule courbe intégrale. Autrement dit, dans le portrait de phase d’un système
différentiel, les courbes intégrales ne peuvent se couper.

1.2 Systèmes différentiels linéaires

1.2.1 Définitions et représentation matricielle

Définition

1. Un système différentiel est dit linéaire si chacune des fonctions ϕi qui le définie est
linéaire :

0
 y1 = a11 (t)y1 + a12 (t)y2 + . . . + a1n (t)yn + b1 (t)

(S) : .. .. ..
. . .
 y 0 = a (t)y + a (t)y + . . . + a (t)y + b (t)

p p1 1 p2 2 pp p p

2. Les fonctions aij sont appelées les coefficients de (S).

3. Les fonctions bi forment le second membre.

Si elles sont toutes nulles, le système est dit homogène.

Le caractère linéaire d’un tel système permet alors d’en donner une représentation
sous la forme d’une matricielle. Précisément, avec les notations ci-dessus, en notant

 
  a11 (t) a12 (t) · · · a1p (t)  
y1 (t) b 1 (t)

..  , A(t) =  a21 (t) a22 (t) · · · a2n (t)   . 
Y (t) =   , B(t) =  .. 
 
.   .. ..
 . . 
yp (t) bp (t)
ap1 (t) ap2 (t) · · · app (t)

on a
10 CHAPITRE 1. SYSTÈMES DIFFÉRENTIELS

(S) ⇐⇒ Y 0 = A(t)Y + B(t)

1.2.2 Structure de l’ensemble des solutions

Le caractère linéaire des systèmes différentiels linéaires permet d’étudier en détail la
structure de l’ensemble Σ des solutions d’un tel système. En effet, on peut montrer (notam-
ment via la notation matricielle) qu’un système différentiel de la forme Y 0 = A(t)Y + B(t)
vérifie le principe de superposition :
Théorème

Soient
(S) : Y 0 = A(t)Y + B(t) et (H) : Y 0 = A(t)Y
un système différentiel linéaire et son système homogène associé.

1. Pour toutes solutions Y1 et Y2 de (H), et pour tout λ ∈ R, la fonction vectorielle

Y : t 7−→ λ.Y1 (t) + Y2 (t)

est une solution de (S).

2. Si Yp est une solution de (S), pour toute solution Yh de (H), la fonction vectorielle

Y : t 7−→ Yp (t) + Yh (t)

est une solution de (H).

Exercice : à démontrer.

D’un point de vue pratique, cela signifie que :

1. L’ensemble Σh des solutions d’un système différentiel linéaire homogène est un es-
pace vectoriel.

2. L’ensemble Σ des solutions d’un système différentiel quelconque est un espace af-
fine, dirigé par l’ensemble Σh des solutions du système homogène associé.

On peut en outre montrer que ces espaces sont de dimension p, soit le nombre d’in-
connue du système. En pratique cela signifie que

1. L’ensemble Σh des solutions d’un système différentiel linéaire homogène à p incon-

nues est engendré par p solutions indépendantes Y1 , . . . , Yp :
1.2. SYSTÈMES DIFFÉRENTIELS LINÉAIRES 11

Σh = {Y : t 7−→ α1 .Y1 (t) + . . . + αp .Yp (t), αk ∈ R}

p
( )
X
= Y : t 7−→ αk .Yk (t), αk ∈ R
k=1

2. Si Yp est solution d’un système différentiel linéaire (H), l’ensemble Σ des solutions
de (S) est

Σ = {Y : t 7−→ Yp (t) + Yh (t), Yh ∈ Σh }

Autrement dit,
1. Pour résoudre un système différentiel linéaire et homogène à p inconnues, on doit
déterminer p solutions vectorielles Y1 , . . . , Yp indépendantes.

2. Pour résoudre un système différentiel linéaire non homogène, on doit résoudre le

système différentiel homogène associé et trouver une solution du système complet.

On va voir dans la suite comment, dans le cas d’un système linéaire homogène à co-
efficients constants Y 0 = AY (avec A ∈ Mn (R)), il peut être possible de construire p
solutions d’un système homogène à l’aide des éléments propres (i.e. valeurs propres et
vecteurs propres) de la matrice A.

Dans le cas d’un système non homogène, la recherche d’une solution particulière Yp
est, en général, un problème compliqué. Bien qu’il existe quelques méthodes spécifiques à
certains types de second membre, la détermination d’une solution particulière s’appuie le
plus souvent sur le problème physique dont est issue le système étudié.

1.2.3 Matrice résolvante

Soit (H) : Y 0 = A(t)Y un système différentiel linéaire homogène à p inconnues.
Là encore, le fait que l’ensemble des solutions Σh de (H) s’exprime sous la forme de
combinaisons linéaires de p solutions de base induit une représentation matriciel de cet
ensemble de solutions. Précisément, si Y1 , . . . , Yp sont p solutions indépendantes de (H),
en notant
 
 
 
 Y1 (t) Y2 (t) · · · Yp (t)
M (t) =  

 
12 CHAPITRE 1. SYSTÈMES DIFFÉRENTIELS

on a
   

 α1 

 ..  p
Σh = Y : t 7−→ M (t) × C, C =  .  ∈ R
 
 αp 

On note en particulier que cet ensemble de solution dépend de p constantes α1 , . . . , αp

que l’on peut choisir librement pour construire une solution de (H). Il est alors possible
d’isoler une unique solution en ajoutant à (H) un ensemble de conditions initiales de la
forme
 
y10
(C.I.) : Y (0) = Y0 =  ...  ∈ Rp
 
yp0

On peut montrer qu’il existe alors une unique matrice R(t) carrée à p lignes telle que
l’unique solution du problème de Cauchy (H) + (C.I.) soit

Y : t 7−→ R(t) × Y0

Cette matrice est appelées la résolvante du système (H).

En pratique, on peut montrer que si M (t) est une matrice de solutions de base, alors

R(t) = M (t) × M (0)−1

En effet, si Y (t) = M (t) × C, alors

Y (0) = Y0 ⇐⇒ M (0) × C = Y0
⇐⇒ C = M (0)−1 × Y0

On a alors

Y (t) = M (t) × M (0)−1 ×Y0

| {z }
R(t)
1.2. SYSTÈMES DIFFÉRENTIELS LINÉAIRES 13

1.2.4 Système linéaires à coefficients constants

Solutions et éléments propres

Un système différentiel linéaires homogène (H) : Y 0 = A(t)Y est dit à coefficients

constants si aucun de ses coefficients aij de la matrice A ne dépend pas de t.

La matrice A de ce système est alors une matrice à coefficients réels (ou complexes) et,
comme annoncé plus haut, il est possible de construire des solutions de (H) en s’appuyant
sur les éléments propres de A.

Rappels : soit A une matrice carrée à p lignes. On appelle vecteur propre de A tout
vecteur V ∈ Rp tel que

V 6= 0 et ∃λ ∈ R / A × V = λ.V

Le coefficient λ de la relation liant A × V et V est alors appelé valeur propre de A

associée au vecteur propre V .

Mais alors, si A est la matrice d’un système différentiel linéaire, homogène à coefficients
constants, pour tout vecteur propre V ∈ Rp de A, en notant λ la valeur propre associée à
V , la fonction vectorielle

Y : t 7−→ eλt .V

est une solution du système Y 0 = AY .

En effet, si Y (t) = eλt .V , on a

Y 0 (t) = λ.eλt .V

et

A × Y (t) = A × (eλt .V ) = eλt .A × V = eλt .λ.V = Y 0 (t)

Ainsi, si A admet p vecteurs propres V1 , . . . , Vp indépendants (autrement dit, si la

matrice A est diagonalisable), en notant λi la valeur propre associée à Vi , les p fonctions
vectorielles

Y1 (t) = eλ1 t V1 , . . . , Yp (t) = eλp t Vp

forment une base de l’ensemble des solutions du système Y 0 = AY et

14 CHAPITRE 1. SYSTÈMES DIFFÉRENTIELS

p
( )
X
Σh = Y : t 7−→ αk .eλk t Vk , αk ∈ R
k=1

Note : avec les notations du paragraphe précédent, la matrice M (t) est ici
 
 
 λt λp t

 e V1 · · · e Vp
M (t) =  1 

 

Valeurs propres complexes

Si la matrice A du système Y 0 = AY admet des valeurs propres complexes, la méthode

présentée ci-dessus reste valable, les exponentielles apparaissant dans les solutions de base
étant alors des exponentielles complexes. Cependant, si la matrice A est à coefficients réels,
les solutions du problème sont alors des fonctions vectorielles réelles et il est alors possible
d’exprimer ces solutions à l’aide des fonctions trigonométriques.
Précisément, si A est à coefficients réels, on peut montrer que si λ ∈ C est une valeur
propre de A, alors λ est aussi une valeur propre de A. D’autre part, si Z ∈ Cn est un
vecteur propre de A associé à λ, alors Z ∈ Cn est un vecteur propre de A associé à λ.
Ainsi, les fonctions

Y1 : t 7→ eλt Z et Y2 : t 7→ eλt Z

sont des solutions de base du système Y 0 = AY .

Cependant, on peut montrer que dans la base de solutions, on peut remplacer ces deux
solutions par

Ỹ1 : t 7→ Re(eλt Z) et Ỹ2 : t 7→ Im(eλt Z)

dont les coefficients sont réels et s’expriment à l’aide des fonctions trigonométriques
sinus et cosinus.

Exemple : déterminer les solutions (réelles) du système différentiel Y 0 = AY avec

 
0 −1
A=
1 0
1.2. SYSTÈMES DIFFÉRENTIELS LINÉAIRES 15

Matrice non diagonalisable

Si la matrice A du système Y 0 = AY n’est pas diagonalisable, pour chaque couple (λ, X)
de (vap,vep) de A, la fonction vectorielle

Y : t 7→ eλt X

est encore une solution du système, mais A ne possède pas assez de vecteurs propres in-
dépendants pour pouvoir extraire n fonctions vectorielles de base.

Deux méthodes différentes permettent alors de déterminer l’ensemble des solutions du

système Y 0 = AY :

• L’algèbre linéaire montre qu’il est possible d’associer à A une matrice triangulaire
supérieure dont les termes diagonaux sont les valeurs propres de A et dont seuls
quelques coefficient de la surdiagonale sont non nuls et peuvent être pris égaux à
1. Il est alors possible de calculer directement la résolvante R(t) = eAt du système
en généralisant la notion d’exponentielle aux matrices à l’aide de la formule
+∞
M
X 1 n
e = M
n=0
n!

• Si la matrice A n’est pas diagonalisable, on ne peut trouver suffisamment de vec-

teurs propres pour constituer une base de Rn . Formellement, cela se produit lorsque
la dimension d d’un sous espace propre Eλ est strictement inférieure à l’ordre α
de la valeur propre λ associée dans la factorisation du polynôme caractéristique de
A. Cependant, à l’image de ce qui se fait en dimension 1, on peut montrer que la
valeur propre λ produit r = α − d solutions différentes, de la forme

Yi (t) = ti eλt Vi , i ∈ {0, . . . , r − 1}

On peut en outre obtenir les Vi ∈ Cn donnant ces solutions en cherchant une

solution particulière au système sous la forme

Y (t) = eλt V0 + tV1 + . . . + tr−1 Vr−1

que l’on injecte dans le système Y 0 = AY .

1.2.5 Classification des portraits de phases (systèmes linéaires

2 × 2)
Dans le cas des systèmes linéaires, homogènes, à coefficients constants, on peut mon-
trer qu’il n’existe qu’un nombre fini de type de portraits de phases possibles. Précisé-
ment, on peut montrer que la forme générale du portrait de phase d’un système Y 0 = AY
où A ∈ M2 (R) dépend

• du signe des valeurs propres de A si celles-ci sont réelles,

16 CHAPITRE 1. SYSTÈMES DIFFÉRENTIELS

• du signe de la partie réelle des valeurs propres de A si celles-ci sont complexes

conjuguées.

Ainsi,

• Si λ1 et λ2 sont réelles et de même signe, le portrait de phase du système (S) est

un nœud.

λ1 , λ2 > 0 λ1 , λ2 < 0

y y

x x

nœud répulsif nœud attractif

• Si λ1 et λ2 sont réelles et de signes opposés, le portrait de phase du système (S)

est un col.

• Si λ1 = a + ib et λ2 = a − ib = λ1 , le portrait de phase de (S) fait apparaître des

1.3. ÉTUDE QUALITATIVE DES SYSTÈMES AUTONOMES 17

phénomènes tournants :

a>0 a=0 a<0

y y y

x x x

spirale divergente cercles spirale convergente

1.3 Étude qualitative des systèmes autonomes

Hormis les système différentiels linéaires, il n’existe que très peu de systèmes différen-
tiels que l’on sait résoudre de façon exacte. Cependant, il est possible, dans certains cas,
d’effectuer une étude qualitative des solutions.

1.3.1 Systèmes autonomes et champs de vecteurs

Définition (système différentiel autonome)

Un système différentiel (S) est dit autonome s’il est de la forme

(S) : Y 0 = F (Y )

Autrement dit, un système différentiel est dit autonome si la fonction F le définissant

(sous forme vectorielle) ne dépend que des inconnues y1 , . . . , yp (et non de t).

Exemple : les système différentiels linéaires, homogènes à coefficients constants sont

des systèmes autonomes. Si l’on ajoute un second membre constant à ce type de système,
on obtient encore un système autonome.

D’un point de vue théorique, la fonction F définissant un système différentiel autonome

est un champ de vecteurs de Rp :

F : Rp −→ Rp
X = (x1 , . . . , xp ) 7−→ F (X) = (f1 (X), . . . , fp (X))
18 CHAPITRE 1. SYSTÈMES DIFFÉRENTIELS

On peut alors montrer que, d’un point de vue géométrique, il existe un lien fort entre
les vecteurs du champ F et les courbes intégrales du portrait de phase de l’équation (S).
Précisément, en chaque point X = (x1 , . . . , xp ) ∈ Rp , le vecteur F (X) est tangent à la
courbe intégrale de (S) passant par X.

Ainsi, dans les cas p = 2 et p = 3, le tracé du champ de vecteur F permet une

première approche du portrait de phase du système (S). Dans le cas des systèmes linéaires
autonomes, on retrouve ainsi les différents portraits de phases possibles :
Enfin, on verra que, dans le cas général (i.e. pour les systèmes non linéaires), bien
qu’il n’existe plus de méthode systématique de résolution, l’étude du champ de vecteurs F
permet d’obtenir de l’information qualitative sur certaines des solutions de (S).

1.3.2 Points d’équilibre est stabilité

Définition (point fixe)

Soit (S) : Y 0 = F (Y ) un système différentiel autonome, défini par un champ de vecteurs

F : Rp −→ Rp
X 7−→ (f1 (X) . . . , fp (X))

On appelle point d’équilibre de (S), ou point fixe de (S) tout vecteur X ∗ ∈ Rp tel que

F (X ∗ ) = 0Rp ⇐⇒ ∀i ∈ {1, . . . , p}, fi (X ∗ ) = 0

D’un point de vue théorique, il est clair que les points d’équilibre d’un système auto-
nome (S) correspondent à ses solutions constantes.

Exemple : si A ∈ Mp (R) est une matrice inversible, le système différentiel Y 0 = AY

admet pour unique point d’équilibre le vecteur X ∗ = 0Rp .

On va voir comment l’étude du champ de vecteurs F au voisinage de ses points fixes

permet d’obtenir de l’information sur les solutions de (S) passant “proche” d’un point
d’équilibre.

Nous verrons en particulier qu’il est possible d’établir une classification des points
d’équilibre d’un système différentiel autonome en fonction du comportement asympto-
tique des solutions voisines de ces points fixes.
1.3. ÉTUDE QUALITATIVE DES SYSTÈMES AUTONOMES 19

Définition (stabilité des points d’équilibres)

Soit (S) : Y 0 = F (Y ) un système différentiel autonome associé à un champ de vecteurs F

de Rp et soit X ∗ ∈ R∗ un point d’équilibre de (S).

1. X ∗ est dit stable si toute solution de (S) passant “proche” de X ∗ reste proche de X ∗ .

2. X ∗ est dit asymptotiquement stable si toute solution de (S) passant “proche” de X ∗

converge vers X ∗ .

Notes :
• Un point d’équilibre X ∗ qui n’est pas stable est dit instable. Dans ce cas, les solu-
tions de (S) passant proche de X ∗ ont tendance à s’éloigner de X ∗ .

• Il existe une définition plus rigoureuse de la notion de stabilité, basée sur la notion
de voisinage dans Rp .

Exemples : si (S) : Y 0 = AY est un système différentiel linéaire associé à une

matrice A ∈ Mp (R) inversible, la nature du point fixe X ∗ = 0Rp dépend des valeurs
propres de A :

• X ∗ = 0Rp est stable si toute valeur propre λ de A vérifie Re(λ) 6 0.

• X ∗ = 0Rp est asymptotiquement stable si toute valeur propre λ de A véri-
fie Re(λ) < 0.

1.3.3 Linéarisation
Si (S) : Y 0 = F (Y ) n’est pas un système linéaire, il n’est en général pas possible de
déterminer de façon explicite les solutions de (S). Cependant, si l’on parvient à détermi-
ner les points d’équilibre de (S), il est possible (comme cela a été annoncé) d’effectuer
une étude qualitative des solutions de (S) en étudiant de prêt le champ de vecteurs F au
voisinage de ses points fixes.

Précisément, soit X ∗ un point fixe de (S). Si les fonctions coordonnées admettent

des dérivées partielles en X ∗ , le champ de vecteurs F admet un développement limité
d’ordre 1 en X ∗ :

F (X) = F (X ∗ ) +∇(F )(X ∗ )(X − X ∗ ) + o(||X − X ∗ ||)

| {z }
=0

où ∇(F )(X ∗ ) est la matrice jacobienne (ou le gradient) de F en X ∗ :

20 CHAPITRE 1. SYSTÈMES DIFFÉRENTIELS

 
∂fi ∗
∇(F )(X ∗ ) = 
 
(X ) 
 ∂xj 

Au voisinage de X ∗ (i.e. tant que Y (t) − X ∗ reste petit), les courbes intégrales de (S)
sont proches des courbes intégrales du système linéaire à coefficients constants

Y 0 = ∇(F )(X ∗ )(Y − X ∗ )

La nature du point X ∗ (stable/asymptotiquement stable/instable) et la forme de ces

courbes dépend alors des valeurs propres de la matrice ∇(F )(X ∗ ).

Exemple : soit
x0 = (1 − y)x

(S) :
y 0 = −(1 − x)y
Le champ de vecteurs F associé à (S) est

F : R2 −→ R2
(x, y) 7−→ ((1 − y)x, −(1 − x)y)

Le point X ∗ = (0, 0) est alors un point d’équilibre de (S) et

 
1 0
∇(F )(0, 0) =
0 −1

Le point X ∗ = (0, 0) est donc un col et il est instable.

1.3.4 Cycles limites et attracteurs étranges

Chapitre 2

Résolution numériques d’équations

différentielles

Introduction

En modélisation, les lois de la physique se traduisent en général sous la forme d’équa-

tions faisant le lien entre la valeur d’une ou plusieurs quantités caractéristiques du pro-
blème étudié et les variations de ces mêmes quantités.

En pratique, ces lois se traduisent donc sous la forme d’équations différentielles. La

résolution de ces équations permet alors de déterminer les grandeurs caractéristiques as-
sociées et comprendre et prédire les phénomènes étudiés.

Formellement, il n’existe que peu d’équations différentielles que l’on sait résoudre de
façon exacte. Elles sont en général issue de problèmes simplifiés, qui permettent une pre-
mière approche du problème étudié.

Pour une étude plus poussée, on doit se rabattre sur des méthodes de résolution ap-
prochée.

Les méthodes que l’on va exposer ici ont été étudiées et développées au XVIIIième
siècle (par Euler, notamment) mais ont connu un essor important depuis l’apparition des
calculateurs.

Elles consistent à construire, à l’aide de suites numériques, des solutions approchées

pour les équations différentielles que l’on ne sait pas résoudre de façon exacte.

L’étude théorique de ces méthodes consiste d’abord à déterminer les processus mathé-
matiques permettant d’aboutir au résultat voulu puis à déterminer des critères permettant
d’évaluer la qualité de ces méthodes.

21
22CHAPITRE 2. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS DIFFÉRENTIELLES

2.1 Position du problème

Dans un premier temps, on considère un problème de Cauchy, défini sur un intervalle
de la forme I = [0, Tmax ] et constitué d’une équation différentielle d’ordre 1 et d’une
condition initiale :  0
y = f (t, y) (E)
(P ) :
y(0) = y0 (C.I.)
La théorie nous assure que, sous certaines hypothèses de régularité portant sur la
fonction f définissant l’équation (E) (hypothèses que l’on supposera toujours vérifiées),
le problème (P ) admet une unique solution, encore notée y.

Précisément, il existe une unique fonction y telle que

∀t ∈ I, y 0 (t) = f (t, y(t)) et y(0) = y0

Les méthodes de résolution commencent par fixer une discrétisation de l’intervalle

d’étude I : pour un entier n ∈ N∗ fixé, on découpe l’intervalle I en n intervalles de même
longueur
Tmax
h=
n
Ces intervalles sont définis par n + 1 points

0 = t0 < t1 < . . . < tn = Tmax

Il est alors possible d’approcher la solution y cherchée, sous la forme d’un nuage de
points (ti , yi ) ∈ R2 tels que

∀i ∈ {0, . . . , n}, yi ≈ y(ti )

2.2 Différences finies

Les différentes méthodes d’approximation varient en fonction de la méthode de calcul
des valeurs yi évoquées plus haut. Elles sont cependant toutes basées sur une approxima-
tion de la notion de dérivée à l’aide de la formule de Taylor : soient

• y une fonction d’une variable, dérivable autant que nécessaire,

• un point t du domaine de y,
• un “décalage” h > 0.

On a alors

h2 00 hn
y(t + h) = y(t) + h.y 0 (t) + y (t) + ... + y (n) (t) + O(hn+1 )
2 n!

À l’ordre 1, la formule de Taylor permet une approximation de la dérivée y 0 (t) en un

point t donné :
2.2. DIFFÉRENCES FINIES 23

Taylor ordre 1 : y(t + h) = y(t) + h.y 0 (t) + O(h2 )

y(t + h) − f (t)
⇒ y 0 (t) = + O(h)
h

D’un point de vue géométrique, cette approximation revient à approcher la tangente

à la courbe Cy en t en une corde :

Cette formule est donc appelée différence finie d’ordre 1, décentrée à droite.

Note : la formule de Taylor permet en outre de mesurer l’erreur commise en fonction

du décalage h. Précisément, l’approximation

y(t + h) − y(t)
y 0 (t) ≈
h

introduit une erreur de l’ordre de grandeur de h.

En s’appuyant encore sur le développement de Taylor, il est possible d’obtenir d’autres

approximations pour la dérivée y 0 (t). Ces différentes formes sont à la base des différentes
méthodes connues.

Ainsi, en appliquant la même formule en −h, on a :

f (x − h) = f (x) − h.f 0 (x) + O(h2 )

f (x) − f (x − h)
⇒ f 0 (x) = + O(h)
h
24CHAPITRE 2. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS DIFFÉRENTIELLES

On peut là encore donner une interprétation géométrique de l’approximation

y(t) − y(t − h)
y 0 (t) ≈
h

Cette formule donne ainsi la différence finie d’ordre 1, décentrée à gauche.

À l’aide du développement de Taylor à l’ordre 2 appliqué en x + h et en x − h, on peut

obtenir une troisième approximation pour la dérivée f 0 (x) :

h2 00
(1) : f (x + h) = f (x) + h.f 0 (x) + .f (x) + O(h3 )
2

h2 00
(2) : f (x − h) = f (x) − h.f 0 (x) + .f (x) + O(h3 )
2
f (x + h) − f (x − h)
(1) − (2) ⇒ f 0 (x) = + O(h2 )
2h

On obtient ainsi une approximation plus précise de la dérivée f 0 (x), que l’on peut là
encore illustrer graphiquement :
2.3. APPLICATION À LA RÉSOLUTION DES EDO 25

Enfin, les formules (1) et (2) ci-dessus permettent également d’obtenir une approxi-
mation pour la dérivée seconde f 00 (x) :

f (x + h) − 2f (x) + f (x − h)
(1) + (2) ⇒ f 00 (x) = + O(h2 )
h2

Note : d’un point de vue théorique, il est bien entendu possible de généraliser ces
approximation à tous ordre. Cependant, il est rare que l’on ait besoin d’appliquer ces
méthodes aux ordres supérieurs à 2.

2.3 Application à la résolution des EDO

2.3.1 Schémas d’Euler (explicite)
Revenons maintenant au problème de Cauchy
 0
y = f (t, y)
(P ) :
y(0) = y0

et son unique solution y, définie sur un intervalle I = [0, Tmax ].

On note encore un découpage de I sous la forme de n intervalles de même longueur

h:

0 = t0 < t1 < . . . < tn = Tmax

La fonction y étant solution de (P ), on a, pour tout i ∈ {0, . . . , n − 1}

y 0 (ti ) = f (ti , y(ti ))

Par ailleurs, la formule de différence finie d’ordre 1 décentrée à droite, appliquée en ti

nous donne

y(ti + h) − y(ti ) y(ti+1 ) − y(ti )

y 0 (ti ) ≈ =
h h

Ainsi, pour tout i ∈ {0, . . . , n − 1}, on a

y(ti+1 ) − y(ti )
≈ f (ti , y(ti )) ⇐⇒ y(ti+1 ) ≈ y(ti ) + h.f (ti , y(ti ))
h
26CHAPITRE 2. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS DIFFÉRENTIELLES

Cette approximation est à l’origine de la formule de récurrence caractérisant la méthode

d’Euler explicite :

yi+1 = yi + h.f (ti , yi )

Cette formule permet, à partir de la valeur initiale y0 , de calculer, de proche en proche,

n valeurs numériques y1 , . . . , yn qui, a priori, sont des valeurs approchées des ordon-
nées y(ti ).

Exemples :

1. Soit
y 0 = 2y

(P1 ) :
y(0) = 1
dont l’unique solution est
y : t 7−→ e2t
La méthode d’Euler, appliqué sur l’intervalle [0, 2], découpé en n = 10 intervalles
réguliers donne le graphe ci-dessous.

Approximation - Methode d'Euler

n = 10, h = 0.2
Solution exacte
50 Solution approchee

40

30

20

10

0
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

2. Soit
y 0 = 2(1 − y)

(P2 ) :
y(0) = 0
dont l’unique solution est
y : t 7−→ 1 − e−2t
La méthode d’Euler, appliqué sur l’intervalle [0, 2], découpé en n = 10 intervalles
réguliers donne le graphe ci-dessous.
2.3. APPLICATION À LA RÉSOLUTION DES EDO 27

Approximation - Methode d'Euler

n = 10, h = 0.2
1.0

0.8

0.6

0.4

0.2
Solution exacte
0.0 Solution approchee
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

3. Soit
y 0 = −2ty

(P3 ) :
y(0) = 1
dont l’unique solution est
2
y : t 7−→ e−t

La méthode d’Euler, appliqué sur l’intervalle [0, 2], découpé en n = 10 intervalles

réguliers donne le graphe ci-dessous.

Approximation - Methode d'Euler

n = 10, h = 0.2
1.0 Solution exacte
Solution approchee

0.8

0.6

0.4

0.2

0.0
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

2.3.2 Schémas numériques

Le schéma d’Euler explicite présenté ci-dessus est un exemple de schéma numérique
permettant de résoudre de façon approchée un problème de Cauchy.

En exploitant les différentes formes de différences finies, il est possible de construire

d’autres types de schémas numériques.
28CHAPITRE 2. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS DIFFÉRENTIELLES

Précisément, on appelle schéma numérique toute méthode de construction permettant

de générer, pour un problème de Cauchy, une approximation de la solution y sous la forme
d’une suite de points {(ti , yi ), i ∈ {0, . . . , n}} telle que

• la suite

t0 = 0 < t1 < . . . < tn = Tmax

réalise un découpage régulier de l’intervalle d’étude [0, Tmax ] en n intervalles de

longueur h = Tmax
n
,

• la suite (y0 , . . . , yn ), construite sur une relation de récurrence de la forme

yi+1 = yi + h.ϕ(ti , yi , h)

vérifie

∀i ∈ {0, . . . , n}, yi ≈ y(ti )

Exemple : le schéma d’Euler explicite correspond à

ϕ : (t, y, h) 7−→ f (t, y)

2.3.3 Consistance, stabilité, convergence

Convergence

Sur les exemples ci-dessus, on constate que, bien que l’allure de la courbe soit glo-
balement respectée, certaines approximations sont meilleures que d’autres, pour le même
intervalle et le même pas de temps.

À l’aide des outils de l’analyse, il est alors possible d’étudier la qualité des approxi-
mations obtenues.
2.3. APPLICATION À LA RÉSOLUTION DES EDO 29

Définition (convergence)

Soient n ∈ N∗ et P = {(ti , yi ), i ∈ {0, . . . , n}} un nuage de points produit par un schéma

numérique et approchant la solution y d’un problème de Cauchy.

• On appelle erreurs locales de l’approximation les distances

(n)
εi = |yi − y(ti )|, i ∈ {0, . . . , n}

• On appelle erreur globale de l’approximation la valeur

(n)
ε(n) = max εi = max {|yi − y(ti )|}
i∈{0,...,n} i∈{0,...,n}

Le schéma numérique associé au nuage de points P est alors dit convergent si

lim ε(n) = 0
n → +∞
h→0

Exemple : soit
y 0 = 2y

(P1 ) :
y(0) = 1
Le schéma d’Euler explicite appliqué pour différentes valeurs de n produit les approxi-
mations suivantes :
Methode d'Euler à différents pas de temps
Solution exacte
25 n = 10, h = 0.2
n = 20, h = 0.1
n = 50, h = 0.04
20 n = 100, h = 0.02

15

10

0
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

Consistance et stabilité
D’un point de vue pratique, les erreurs locales apparaissant dans un schéma numérique
(i.e. la différence entre les valeurs calculées yi et les valeurs souhaitées y(ti )) sont issues
30CHAPITRE 2. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS DIFFÉRENTIELLES

de deux sources :

1. Les erreurs “théoriques”, appelées erreurs de consistance, issues du schéma choisi.

2. Les erreurs numériques, issues de l’incertitude inhérente à tout calcul effectué par
un ordinateur.

Ces deux notions sont formalisées dans les définitions ci-dessous.

Définition (Consistance)

Considérons un schéma numérique basé sur une relation de la forme

yi+1 = yi + h.ϕ(ti , yi , h)

1. L’erreur de consistance locale du schéma (au temps ti ), relativement à une fonction

y est la différence

i (y) = y(ti+1 ) − (y(ti ) + +h.ϕ(ti , y(ti ), h))

2. L’erreur de consistance du schéma, relativement à y est

n−1
X
(n)
 (y) = |i (y)|
i=0

3. Le schéma numérique est alors dit consistant si, pour toute solution y de l’équation
différentielle y 0 = f (t, y), on a

lim (n) (y) = 0

n→+∞

4. Le schéma numérique est dit d’ordre d ∈ N∗ si

 
(n) 1
 =O
nd

Par ailleurs, un schéma numérique est dit stable si une petite perturbation dans lors
des calculs successifs ne conduit pas à une divergence de la méthode, malgré son caractère
consistant. Précisément,
2.3. APPLICATION À LA RÉSOLUTION DES EDO 31

Définition (stabilité)

Un schéma numérique, associé à une fonction ϕ, est dit stable s’il existe une constante
S (appelée constante de stabilité) telle que pour toutes suites (yn ) et (ỹn ) vérifiant les
relations 
yi+1 = yi + h.ϕ(ti , yi , h)
ỹi+1 = ỹi + h.ϕ(ti , ỹi , h) + δi
on a !
Xn
max |yi − ỹi | 6 S |y0 − ỹ0 | + |δi |
i∈{0,...,n}
i=0

Conditions de convergence
Considérons un schéma numérique construit sur une relation de récurrence de la forme
yi+1 = yi + h.ϕ(ti , yi , h)
appliqué à une équation différentielle
(E) : y 0 = f (t, y)
On peut montrer qu’un tel schéma, appliqué à la résolution de (E) est consistant si et
seulement si
ϕ(t, y, 0) = f (t, y)
C’est en particulier le cas pour le schéma d’Euler explicite puisque dans ce schéma,
on a
ϕ(t, y, h) = f (t, y) ∀h > 0
On peut également montrer qu’un tel schéma, appliqué à (E) est stable sur l’inter-
valle [0, Tmax ] si ϕ est lipschitzienne en y, i.e. il existe une constante L > 0, dite constante
de Lipschitz telle que
∀(x, y), |ϕ(t, x, h) − ϕ(t, y, h| 6 L.|x − y|
Dans ce cas, la constante de stabilité est
S = [Link]
Note : le caractère lipschitzien induit, en pratique, une contrainte sur les variations
de la fonction ϕ, relativement à la variable y. Cette condition est en réalité la condition
principale à l’existence et l’unicité d’une solution au problème de Cauchy. Cette condition
n’apporte donc pas de contrainte supplémentaire dans le cadre qui nous occupe.

D’autre part on peut noter que si ϕ est continue par rapport à chacune de ses variables,
le fait de travailler sur un intervalle fermé borné de la forme [0, Tmax ] induit automatique-
ment le caractère lipschitzien de ϕ.
32CHAPITRE 2. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS DIFFÉRENTIELLES

Enfin, notons que toutes ces notions sont bien cohérentes avec la notion de convergence
(qui est en définitive la seule condition qui importe) via le théorème suivant :
Théorème (Convergence d’un schéma numérique)

Tout schéma numérique stable et consistant est convergent.

Chapitre 3

Résolution numériques d’équations aux

dérivées partielles

3.1 Notions de base

Comme cela a été évoqué en introduction de ce cours, les équations aux dérivées par-
tielles (EDP) sont au cœur de la modélisation des phénomènes physiques.

L’un des rôles des mathématiques dans cette modélisation consiste à donner des outils
de résolution pour ce type d’équations.

Cela passe, comme à chaque fois, par une étude théorique des EDP, permettant dans
un premier temps d’établir différentes classifications, correspondant aux différents types
d’équations que l’on peut obtenir grâce aux lois de la physique.

L’objectif de ces classifications est de pouvoir développer des méthodes de résolutions

spécifiques à chaque classe, en s’appuyant à chaque fois sur les caractéristiques communes
aux équations de chaque classe.

Dans le cas des équations aux dérivées partielles, il n’existe en réalité que très peu
d’équations que l’on sait résoudre de façon exacte. Cependant, les efforts de classification
ont tout de même permis de mettre en évidence certains liens théorique existant entre
différents phénomènes physiques à priori indépendants.

Par ailleurs, cette classification se retrouve également dans les méthodes de résolution
approchée, qui font l’objet principal de ce cours.

3.1.1 Définitions et notations

Une équation aux dérivées partielles est une équations dont l’inconnue est une fonction
u de plusieurs variables (x, y, z, t, . . .), liant entre elles les dérivées partielles de u.

Notation : soit u une fonction dépendant des variables (x, y, z, t, ...). Pour alléger les
notations, on notera

33
34CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

∂u ∂ 2u ∂ 2u ∂ 2u
ux = , utt = , uxy = = ,...
∂x ∂t2 ∂x∂y ∂y∂x

Parmi les caractéristiques principales d’une EDP, on appelle ordre d’une EDP la plus
haute dérivée partielle intervenant dans l’équation. Les lois de la physique produisent en
général des équations d’ordre 2 (via notamment le Laplacien).

Comme dans le cas des équations différentielles ordinaires (EDO), on distingue éga-
lement les EDP linéaires, i.e. linéaires en l’inconnue u et ses dérivées. Une EDP linéaire
sera de plus dite à coefficients constants si les coefficients de la combinaison linéaire (qui
en toute généralité peuvent dépendre des inconnues du problème) sont des réels fixés.

On a vu dans les années précédentes qu’il existait des méthodes de résolution spéci-
fiques aux EDP linéaires d’ordre 1 à coefficients constants. Il est parfois possible d’étendre
ces techniques à certaines EDP linéaires (toujours d’ordre 1) à coefficients non constants.

Dans la suite, on verra que le caractère linéaire des EDP auxquelles on appliquera les
méthodes de résolution approchée est fondamental pour la mise en œuvre de ces méthodes.

3.1.2 Exemples d’EDP classiques

• Équation de Laplace 2D :

uxx + uyy = 0 ⇐⇒ ∆u = 0

Elle apparaît dans tous les domaines de la physique. Ainsi, si u représente une
distribution de chaleur dans un domaine 2D, cette équation régit par exemple la
répartition de chaleur dans le domaine, une fois l’équilibre atteint.

• Équation de diffusion de la chaleur régissant l’évolution de la température dans un

milieu homogène :

1D : ut = α2 uxx

2D : ut = α2 (uxx + uyy )

3D : ut = α2 (uxx + uyy + uzz )

Note : la forme générale de cette équation est

ut = ∆u

• Équation de propagation des ondes 3D :

3.1. NOTIONS DE BASE 35

utt = c2 (uxx + uyy + uzz ) ⇐⇒ utt = c2 ∆u

• Équation de Poisson (électrostatique, gravitation, ets) :

∆u = f (x, y, z)

Notes :

• Toutes les équations ci-dessus sont linéaires, à coefficients constants.

• Hormis la dernières, toutes les équations ci-dessus sont homogènes. Le terme f (x, y, z)
de l’équation de Poisson est appelé second membre de l’équation. Dans la modé-
lisation, il correspond en général à un apport extérieur au système, au cours de
l’expérience.

• Bien que formellement, toutes les variables d’un problème soient équivalentes, on
peut, d’un point de vue physique, distinguer les variables d’espace (x, y, z) de la
variable de temps (t). Cette distinction permet en particulier de distinguer deux
types de problèmes :

— Les problèmes stationnaires ou problèmes d’équilibre, indépendants du temps :

l’équation de Laplace et l’équation de Poisson.

— Les problèmes dynamiques (ou transitoires) dépendant du temps : les quatre

autres exemples ci-dessus.

On verra dans la suite que l’on retrouve cette distinction dans les différentes mé-
thodes de résolution approchée. On verra également réapparaître cette distinction
lorsque l’on évoquera les notions de conditions limites et/ou conditions initiales.

3.1.3 Classification des EDP linéaires d’ordre 2

Sous sa forme générale, une EDP linéaire d’ordre 2 portant sur deux variables (x, y)
est de la forme

(E) : Auxx + Buxy + Cuyy + Dux + Euy + F u = G

Par analogie avec les polynômes de degré 2 en deux variable et la classification des
formes quadratiques, on peut classer ces équations aux dérivées partielles en fonction des
coefficients A, B et C. Précisément,

• (E) est dite elliptique si B 2 − 4AC < 0.

36CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

• (E) est dite parabolique si B 2 − 4AC = 0.

• (E) est dite hyperbolique si B 2 − 4AC > 0.

Outre qu’il existe différentes méthodes de résolution liées à cette classification, chaque
type d’équation peut également être associé à une classe de problèmes physiques ana-
logues :

• Les équations elliptiques apparaissent dans de nombreux problèmes stationnaires

(i.e. indépendants du temps).

• Les équations paraboliques apparaissent dans des problèmes de conduction et dif-

fusion (chaleur, réactions chimiques, etc).

• Les équations hyperboliques apparaissent dans des problèmes de propagation d’ondes

(vibration, vagues, etc).

3.1.4 Conditions initiales et conditions limites

Comme pour les EDO, chaque EDP admet une infinité de solutions. Du point de vue
de la modélisation, cela signifie que la réponse attendue dépend de conditions supplémen-
taires portant d’abord sur l’état du système étudié aux bords ∂D du domaine d’étude D :
les conditions limites.

Ces conditions limites peuvent être de trois types différents selon le problème étudié :

• Conditions limites type Dirichlet (CLD) : la valeur de la fonction u est imposée

sur l’ensemble de la frontière ∂D du domaine d’étude D.
3.1. NOTIONS DE BASE 37

Exemple : CLD pour un problème 1D sur un segment de longueur L ramené au

segment [0, L] :

 u(t, 0) = u0 (t) ∀t > 0

u(t, L) = uL (t) ∀t > 0



• Conditions limites type Neumann (CLN) : la valeur de la dérivées directionnelle

∂u
normale est imposée sur l’ensemble de la frontière ∂D.
∂n

Exemple : CLN pour un problème 1D sur le segment [0, L] :


∂u
 ∂x (t, 0) = v0 (t) ∀t > 0





 ∂u

 (t, L) = vL (t) ∀t > 0
∂x

Du point de vue du modèle, cette condition correspond à “la porosité” de la fron-

tière ∂D.

• Conditions limite mixtes : des conditions limites de la forme

∂u
∀X0 ∈ ∂D, ∀t > 0, a.u(t, X0 ) + b. (t, X0 ) = u0 (t)
∂n
38CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

sont imposées au bord de D.

Par ailleurs, dans le cas d’un problème dynamique (i.e. dépendant du temps t), on peut
également ajouter des conditions initiales, donnant l’état général du système au début de
l’expérience (i.e. à t = 0).

Dans ce cas, la nature de ces conditions initiales dépendent alors de l’ordre de la plus
haute dérivée temporelle (ut , utt , etc) intervenant dans l’équation :

• Si (E) est d’ordre 1 en temps, u(X, 0) est imposée pour tout X ∈ D.

∂u
• Si (E) est d’ordre 2 en temps, u(X, 0) et (X, 0) sont imposées pour tout X ∈ D.
∂t
Exemple : problème transitoire 1D sur le segment [0, L] :

3.2 Différences finie

Faute de méthode générale permettant de résoudre toutes les EDP, il a fallu dévelop-
per des méthodes de résolution approchée.

L’objectif est ici de présenter la méthode des différences finies, permettant de trans-
former une EDP donnée en un ensemble de systèmes d’équations algébriques.

Le cœur de la méthode est basé sur une approximation de la notion de dérivée à l’aide
de la formule de Taylor : soient

• f une fonction d’une variable, dérivable autant que nécessaire,

• un point x du domaine de f ,
3.2. DIFFÉRENCES FINIE 39

• un “décalage” h > 0.

On a alors

h2 00 hn
f (x + h) = f (x) + h.f 0 (x) + f (x) + ... + f (n) (x) + O(hn+1 )
2 n!

Ainsi, à l’ordre 1, la formule de Taylor permet une approximation de la dérivée f 0 (x) :

Taylor ordre 1 : f (x + h) = f (x) + h.f 0 (x) + O(h2 )

f (x + h) − f (x)
⇒ Dérivée première décentrée à droite : f 0 (x) = + O(h)
h

D’un point de vue géométrique, cette approximation revient à approcher la tangente

à Cf en x en une corde :

La formule de Taylor permet alors de mesurer l’erreur commise en fonction du déca-

lage h.

En appliquant la même formule en −h, on obtient une seconde approximation possible :

Taylor ordre 1 : f (x − h) = f (x) − h.f 0 (x) + O(h2 )

f (x) − f (x − h)
⇒ Dérivée première décentrée à gauche : f 0 (x) = + O(h)
h
40CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

À l’aide du développement de Taylor à l’ordre 2 appliqué en x + h et en x − h, on peut

obtenir une troisième approximation pour la dérivée f 0 (x) :

h2 00
(1) : f (x + h) = f (x) + h.f 0 (x) + .f (x) + O(h3 )
2

h2 00
(2) : f (x − h) = f (x) − h.f 0 (x) + .f (x) + O(h3 )
2
f (x + h) − f (x − h)
(1) − (2) ⇒ f 0 (x) = + O(h2 )
2h

On obtient ainsi une approximation plus précise de la dérivée f 0 (x), que l’on peut là
encore illustrer graphiquement :

Enfin, les formules (1) et (2) ci-dessus permettent également d’obtenir une approxi-
mation pour la dérivée seconde f 00 (x) :
3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 41

f (x + h) − 2f (x) + f (x − h)
(1) + (2) ⇒ f 00 (x) = + O(h2 )
h2

Il est bien entendu possible de généraliser ces approximation à tous ordre.

D’autre part, il est également possible d’adapter rapidement ces formules aux dérivées
partielles de fonctions de plusieurs variables. Ainsi, si u est une fonction de deux variables
(x, y), on a

∂u u(x + h, y) − u(x, y)
(x, y) = ux (x, y) ≈
∂x h
u(x, y + k) − u(x, y)
uy (x, y) ≈
k
u(x + h, y) − 2u(x, y) + u(x − h, y)
uxx (x, y) ≈
h2
u(x, y + k) − 2u(x, y) + u(x, y + k)
uyy (x, y) ≈
k2

etc...

3.3 Illustrations de la méthode des différences finies

3.3.1 Équation de la chaleur 2D, stationnaire

On souhaite dans cet exemple résoudre l’équation de la chaleur stationnaire sur un do-
maine carré D du plan, sous des conditions limites de Dirichlet. Précisément, considérons
le carré [0, 1]2 du plan muni d’un repère (xOy).
42CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

y
u(x, 1) = 0
1

u(1, y) = 0
u(0, y) = 0
∆u(x, y) = 0

0 x
0 u(x, 0) = sin(πx) 1
On souhaite connaître la température u(x, y) en tout point (x, y) ∈ D une fois l’équi-
libre atteint. Cette distribution de température est déterminée par

• l’équation de la chaleur à l’intérieur du domaine D :

∀(x, y) ∈]0, 1[2 , ∆u(x, y) = 0

• les conditions limites :

Nord : ∀x ∈ [0, 1], u(x, 1) = 0

Sud : ∀x ∈ [0, 1], u(x, 0) = sin(πx)

Est : ∀y ∈ [0, 1], u(1, y) = 0

Ouest : ∀x ∈ [0, 1], u(0, y) = 0

Discrétisation du problème
Faute de pouvoir obtenir explicitement la fonction u, on remplace le calcul exacte en
chaque point du domaine D par un calcul (approché) en un nombre fini de points (xj , yi ).

Pour cela, on commence par discrétiser le domaine D :

• On découpe l’intervalle d’abscisses [0, 1] ⊂ (0x) en p intervalles de même longueur

3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 43

1
h =
p
(p ∈ N∗ fixé).

On obtient ainsi p + 1 bornes x0 < . . . < xp vérifiant

j
∀j ∈ {0, . . . , p}, xj = j.h =
p

• On découpe l’intervalle d’ordonnées [0, 1] ⊂ (Oy) en n intervalles de même longueur

1
k =
n
(n ∈ N∗ fixé).

On obtient ainsi n + 1 bornes y0 < . . . < yn vérifiant

i
∀i ∈ {0, . . . , n}, yi = i.k =
n

On obtient ainsi un maillage de (n + 1).(p + 1) points du domaine D :

M = {(xj , yi ), 0 6 i 6 n, 0 6 j 6 p}
44CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

y
yn

(xj , yi )
yi

(x1 , y1 )
y1

y0 x
x0 x1 xj xp

On cherche alors à calculer uij = u(xj , yi ) en chacun de ces points :

y
un0

ui + 1, j

ui0 ui, j − 1 uij ui, j + 1

ui − 1, j

u11
u10

x
u00 u01 u0j u0p

Les conditions limites donnent les valeurs exactes de la température en chacun des 2n + 2p
points situés aux bords de D :
3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 45

Condition Ouest Conditions Est

∀i ∈ {0, . . . , n}, ui0 = 0 uin = 0

Condition Sud Conditions Nord

∀j ∈ {0, . . . , p}, u0j = sin(π.xj ) = sin(π.j.h) upj = 0

y
un1 unj unp
un0

ui0 uip

u10 u1p

x
u00 u01 u0j u0p
Par ailleurs, les approximations des différentes dérivées partielles de u à l’aide des
formules de différences divisées permettent de traduire l’équation ∆u = 0 en un système
d’équations linéaire portant sur les valeurs uij cherchées. En effet, en tout point (x, y) ∈ D,
on a
u(x + h, y) − 2u(x, y) + u(x − h, y)
uxx (x, y) ≈
h2
u(x, y + k) − 2u(x, y) + u(x, y − k)
uyy (x, y) ≈
k2
Aux points du maillage intérieurs au domaine D, on a donc, pour tout i ∈ {1, . . . , n − 1}
et pour tout j ∈ {1, . . . , p − 1} :

 uxx (xj , yi ) ≈ ui,j+1 − 2uij + ui,j−1

h2


 u − 2uij + ui−1,j
 uyy (xj , yi ) ≈ i+1,j


k2
46CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

et

(Eij ) : ∆u(xj , yi ) = 0 ⇐⇒ uxx (xj , yi ) + uyy (xj , yi ) = 0

ui,j+1 − 2uij + ui,j−1 ui+1,j − 2uij + ui−1,j

⇐⇒ + =0
h2 k2

⇐⇒ k 2 (ui,j+1 + ui,j−1 ) − 2(h2 + k 2 )uij + h2 (ui+1,j + ui−1,j ) = 0

On transforme ainsi l’équation de la chaleur en un système linéaire de (n − 1).(p − 1)

équations portant sur les (n − 1).(p − 1) valeurs uij encore inconnues.

Dans la suite, on prendra n = p (et donc h = k) afin de simplifier les formules. Les
équations du système linéaires sont donc

∀i, j ∈ {1, . . . , n − 1}, ui,j+1 + ui,j−1 − 4uij + ui+1,j + ui−1,j = 0

Note : dans le système ci-dessus, on note que chaque équation (Eij ) peut être associée
au point (xj , yi ) du maillage et fait intervenir la température uij en ce point ainsi que les
températures aux quatre points les plus proches. On peut alors noter que

ui,j+1 + ui,j−1 + ui+1,j + ui−1,j

(Eij ) ⇐⇒ uij =
4

Cela traduit le fait qu’à l’équilibre, la température en chaque point du maillage est
égale à la moyenne des températures aux quatre points du maillage les plus proches.
3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 47

y
un1 unj unp
un0

ui + 1, j

ui0 ui, j − 1 uij ui, j + 1 uip

ui − 1, j

u10 u1p

x
u00 u01 u0j u0p

Construction du système linéaire

Pour terminer la résolution approchée du problème étudié, il reste à déterminer les ma-
trices du système linéaire (S) : M U = B correspondant aux équations (Eij ) obtenues lors
de la discrétisation. Le système portant sur (n − 1)2 inconnues liées par (n − 1)2 équations,
2
la matrice M est carrée à (n−1)2 lignes et les vecteurs U et B sont des vecteurs de R(n−1) .

Exemple : pour n = 3, on a le maillage présenté sur la figure ci-dessous. Les conditions

limites imposent les valeurs de la température aux 16 points du bord et les inconnues du
système à résoudre correspondent aux températures des 4 points du maillage intérieur au
domaine, parcours ligne par ligne du bas vers le haut.

y
u31 u32 u33
u30

u21 u22
u20 u23

u11 u12
u10 u13

x
u00 u01 u02 u03

L’équation de Laplace en chacun des quatre points inconnus donne :

48CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

• En (x1 , y1 ) :

(E11 ) : u01 + u10 − 4u11 + u12 + u21 = 0

• En (x2 , y1 ) :

(E12 ) : u02 + u11 − 4u12 + u13 + u22 = 0

• En (x1 , y2 ) :

(E21 ) : u11 + u20 − 4u21 + u22 + u31 = 0

• En (x2 , y2 ) :

(E22 ) : u12 + u21 − 4u22 + u23 + u32 = 0

les températures encadrées correspondant aux points du bord (connues).

On peut alors trier ce système linéaire de quatre équations portant sur les quatre
inconnues (u11 , u12 , u21 , u22 ) :



 −4u11 + u12 + u21 = −u01 − u10




u11 − 4u12 + u22 = −u02 − u13



(S) ⇔
u11 − 4u21 + u22 = −u20 − u31








u12 + u21 − 4u22 = −u23 − u32


sin π3

   

−4 1 1 0 u11
 1 −4 2π
0 1  u12  = −  sin 3
  
⇔   
forme matricielle  1 0 −4 1  u21   0 
0 1 1 −4 u22 0
| {z }| {z } | {z }
M U B

En confiant la résolution de ce système à la machine, on obtient

 
√ 3
3 3 

U=
16  1 
1
que l’on peut représenter par une matrice de couleur :
3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 49

y
u31 u32 u33
u30

u21 u22
u20 u23

u11 u12
u10 u13

x
u00 u01 u02 u03

Dans le cas général (i.e. pour n ∈ N>2 quelconque), le même procédé produit un sys-
tème linéaire M U = B similaire, de taille (n − 1)2 où

• M est la matrice par blocs, carrée à (n − 1)2 lignes suivante :

 
A In−1 On−1 · · · On−1
... .. 
 In−1 A In−1 . 

 . . . On−1 

M =   On−1 In−1 A  ∈ M(n−1)2 (R)
 . ... ... ...
 ..

I 
On−1 · · · On−1 In−1 A

où A est la matrice carrée à n − 1 lignes définie par

 
−4 1 0 ··· 0
... .. 
 1 −4 1 . 

 . . .. 
A =   0 .. .. . 0  ∈ Mn−1 (R)
 . .
 .. . . . . . ... 
1 
0 ··· 0 1 −4

I est la matrice identité à n − 1 lignes et O est la matrice carrée nulle à n − 1 lignes.

• U est le vecteur inconnu

50CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

 
u11
 u12 
..
 
.
 
 
u1,n−1
 
 
 
 u21 
 
 u22 
 .. 
2
U = − .  ∈ R(n−1)
 
u2,n−1
 
 
 .. 

 . 


 un−1,1 


 un−1,2 

 .. 
 . 
un−1,n−1

• B est le second membre du système, construit à partir des conditions limites :

 
u10 + u01
 u02 
..
 
.
 
 
u0,n−2
 
 
 u0,n−1 + u1n 
 
 
 u20 
 
 0 

 .
..


2
B = −  ∈ R(n−1)
 
 0 
u
 
 2n 
 .. 

 . 

 u + u 
 n−1,0 n1 

 un2 


 .
..


 
 un,n−2 
un,n−1 + un−1,n

Notes :
• La matrice M du système à inverser est creuse (i.e. elle contient de très nombreux
zéros). L’utilisation d’algorithmes de résolution spécifique à ce type de matrice per-
met de réduire considérablement le nombre d’opérations nécessaires à son inversion
et permet ainsi de gagner en précision.

• Si l’on remplace les conditions limites de Dirichlet par des conditions de Neumann
3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 51

∂u
de la forme ∂n = Cste, on doit approcher ces dérivées normales par des diffé-
rences finies (décentrée) faisant là encore intervenir les températures aux points
du maillage. On ajoute ainsi au système autant d’équations que nécessaires ; ce
qui modifie quelque peut la matrice M , le second membre B ainsi que le vecteur
inconnue U .

• Il est également possible de tenir compte d’un éventuel second membre

∆u = f (x, y)

Les valeurs de f aux points du maillage (i.e. les valeurs fij = f (xj , yi ) connues)
apparaissent alors dans le second membre B du système ci-dessus.
52CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

3.3.2 Équation de la chaleur 1D instationnaire

Il est également possible d’appliquer la méthode des différences finie à la résolution

d’équations paraboliques, portant par exemple sur l’évolution au cours du temps de la tem-
pérature au sein d’une barre homogène de longueur 1. Précisément, on souhaite connaître
la température u(x, t) en tout point de l’intervalle x ∈ [0, 1] et à chaque instant t > 0,
connaissant la température en chaque point à t = 0 et la température aux extrémités
pour chaque t > 0. La fonction u de deux variables (x, t) est alors l’unique solution du
problème constitué des contraintes suivantes.

t
• L’équation de la chaleur 1D :

∀x ∈]0, 1[, ∀t > 0, ut (x, t) = uxx (x, t)

• Les conditions limites :

u(1, t) = 1
u(0, t) = 1
∀t > 0, u(0, t) = u(1, t) = 1 ut (x, t) = ∆u(x, y)
• Les conditions initiales :

∀x ∈]0, 1[, u(x, 0) = 0

0 x
0 u(x, 0) = 0 1

Discrétisation du domaine

Là encore, la résolution approchée à l’aide des différences finies commence par une
discrétisation du domaine de résolution. Pour cela, on commence par découper l’inter-
valle [0, 1] en n intervalles de même longueur h = n1 . On obtient ainsi n + 1 abscisses
distinctes

xj = j.h, j ∈ {0, . . . , n}

On choisit ensuite un pas de temps k = ∆t qui permettent une discrétisation de R+ :

ti = i.k, k∈N

On obtient ainsi un maillage du plan (xOt)

3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 53

ui0 uin
h

©
©
k
u10 u1n

x
u00 u01 u0j u0n

les conditions limites et initiales permettant de connaître les valeurs de la température

aux bords sud, est et ouest du domaine :

∀j ∈{1, . . . , n − 1}, u0j = 0

∀i ∈ N, ui0 = uin = 1

En exprimant les dérivées partielles de u à l’aide des différences finies, on peut alors
discrétiser l’équation de la chaleur, à l’image de ce que l’on a fait dans le cas stationnaire.
Cependant, dans le cas instationnaire, il existe différentes méthodes de discrétisation.

Discrétisation de l’équation (méthode explicite)

La méthode explicite est basée sur une approximation de la dérivée première ut par
une différence finie décentrée à droite (et la différence finie d’ordre 2 pour uxx ) :

u(x, t + k) − u(x, t)
∀(x, t), ut (x, t) ≈
k
u(x − h, t) − 2u(x, t) + u(x + h, t)
uxx (x, t) ≈
h2

En appliquant ces approximations aux point (xj , yi ) du maillage intérieurs au domaine,

on obtient, pour tout i > 0 et pour tout j ∈ {1, . . . , n − 1},
54CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

ui+1,j − uij ui,j−1 − 2uij + ui,j+1

ut (xj , ti ) = uxx (xj , yi ) ⇐⇒ =
k h2

k
⇐⇒ ui+1,j = uij + (ui,j−1 − 2uij + ui,j+1 )
h2

⇐⇒ ui+1,j = δui,j−1 + (1 − 2δ)uij + δui,j+1 (Eij )

k
où δ= .
h2
Cette mise en équation permet de montrer que la température en un point xj au
temps ti+1 = ti + ∆t dépend de la température aux points xj−1 , xj et xj+1 au temps ti :
t

ui + 1, j

ui0 ui, j − 1 uij ui, j + 1 uin

u10 u1n

x
u00 u01 u0j u0n
Les conditions initiales donnant la température en chaque point à l’instant t = 0, on
peut ainsi calculer, par récurrence, la température en tout point en tout temps.

Formellement, on peut alors représenter l’ensemble des équations (Eij ) sous la forme
d’équations vectorielles. Précisément, à i ∈ N fixé, on note
 
ui1
Ui =  ...  ∈ Rn−1
 
ui,n−1
3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 55

Les conditions initiales donnent alors la valeur du vecteur U0 et l’ensemble des équa-
tions {(Eij ), j ∈ {1, . . . , n − 1}} se traduit alors sous la forme d’une relation de récur-
rence matricielle, à l’aide de laquelle on peut calculer les vecteurs Ui de proche en proche :

Ui+1 = Mδ × Ui + Vi

où
 
1 − 2δ δ 0 ··· 0 
ui0

.. ..
δ 1 − 2δ δ . . 0
   
   

Mδ =  ... ... ... 
et Vi = δ. 
 .. 
 0 0 
 . 
.. .. .. ..
 

. . .
  0
. δ
  
0 ··· 0 δ 1 − 2δ uin

Discrétisation de l’équation (méthodes semi-implicites)

D’un point de vue qualitatif, la méthode ci-dessus étant issue des différences finies, on
peut l’erreur théorique εij à chaque pas de temps ti , au point xj vérifie

εij = O(k + h2 )

Autrement dit, la méthode ci-dessus est

• D’ordre 1 en temps.

• D’ordre 2 en espace.

Cependant, bien que naturelle et relativement simple à mettre en œuvre, elle peut être
particulièrement instable.

Précisément, on peut montrer que, pour éviter que les erreurs de calculs issus des
approximations numériques, le paramètre δ vérifie δ 6 12 . Ainsi, il faut que

1
Condition CFL : k 6 h2
2

Or la précision de la solution obtenue dépendant en particulier du nombre n + 1 de

points de discrétisation, une forte augmentation de n imposer restriction importante sur
la taille du pas de temps k et augmente d’autant le nombre de calculs nécessaires sur un
même intervalle de temps.
56CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

On peut palier à cette instabilité à l’aide de méthodes implicites ou semi-implicites.

Précisément, la méthode implicite consiste à approcher la dérivée première ut (x, t) par

la différence finie décentrée à gauche :

u(x, t) − u(x, t − k) u(x, t + k) − u(x, t)

ut (x, t) ≈ ⇔ ut (x, t + k) ≈
k k

En appliquant cette approximation aux points du maillage (xj , yi+1 ), 1 6 j 6 n − 1,

on remplace l’équation de la chaleur en un ensemble d’équations linéaires :

ui+1,j − uij ui+1,j+1 − 2ui+1,j + ui+1,j−1

ut (xj , yi+1 ) = uxx (xj , yi+1 ) ⇔ =
k h2

⇔ − δui+1,j+1 + (2δ + 1)ui+1,j − δui+1,j−1 = uij

 
  ui0
ui1 
 0 

En posant comme dans la méthode explicite Ui =  .. ..
 et Vi = δ.  , on
   
.  . 
u1,n−1  0 
uin
rassemble l’ensemble de ces équations en un système vectoriel

M−δ × Ui+1 = Ui + Vi+1

où
 
1 + 2δ −δ 0 ··· 0
. ..
−δ 1 + 2δ −δ . . .
 
 
M−δ

=  ... ... 
 0 −δ 0 

 .. .. .. .. 
 . . . . −δ 
0 ··· 0 −δ 1 + 2δ

On peut alors, comme dans le cas explicite, déterminer les températures au temps ti+1
à partir des températures au temps ti et des conditions limites, au prix ici de l’inversion
de la matrice M−δ . Du point de vue de l’efficacité, le coup de cette matrice est largement
compensé par la possibilité d’augmenter le pas de temps k (pour un même pas d’espace
h). On peut en particulier s’affranchir de la contrainte k 6 12 h2 .

Notons que là encore, la méthode est, d’un point de vue théorique, d’ordre 1 en temps
3.3. ILLUSTRATIONS DE LA MÉTHODE DES DIFFÉRENCES FINIES 57

et d’ordre 2 en espace.

Enfin, il est possible d’augmenter encore l’efficacité de la méthode en choisissant un

procédé hybride de discrétisation de l’équation de Laplace :

ui+1,j − uij ui,j−1 − 2uij + ui,j+1

Schéma explicite : =
k h2
ui+1,j − uij ui+1,j+1 − 2ui+1,j + ui+1,j−1
Schéma implicite : =
k h2
ui+1,j − uij ui+1,j+1 − 2ui+1,j + ui+1,j−1
θ-Schéma (θ ∈ [0, 1]) : = θ.
k h2
ui,j−1 − 2uij + ui,j+1
+ (1 − θ)
h2

Note : concernant le θ-schéma,

• On retrouve le schéma explicite pour θ = 0.

• On retrouve le schéma implicite pour θ = 1.

Du point de vue du problème modélisé, les θ-schémas permettent de prendre en compte,

en chaque point de la barre étudiée, la température aux points voisins au temps présent
en plus du temps précédent.

Par ailleurs, pour θ = 21 , on obtient le schéma de Crank-Nicholson. On peut alors mon-

trer que ce schéma est d’ordre 2 en espace et en temps et qu’il est inconditionnellement
stable, i.e. on peut s’affranchir de la contrainte CFL.

Dans le schéma de Crank-Nicholson se traduit alors, chaque équation s’écrit

δ δ δ δ
− .ui+1,j+1 + (1 + δ)ui+1,j − .ui+1,j−1 = .ui,j+1 + (1 − δ) .uij + .ui,j−1
2 2 2 2

On peut alors représenter l’ensemble de ces équations par le système vectoriel

1
M− δ × Ui+1 = M δ × Ui + (Vi + Vi+1 )
2 2 2
où les matrices Mδ et Vi ont été définies aux paragraphes précédents.

On obtient là encore la température en chaque point de la barre au temps ti+1 au prix

d’une inversion de matrice carrée à n − 1 lignes.
58CHAPITRE 3. RÉSOLUTION NUMÉRIQUES D’ÉQUATIONS AUX DÉRIVÉES PARTIELLES

3.4 Résolution exacte par séparation des variables

Bien qu’il ne soit pas possible de déterminer de façon analytique les solutions de la
plupart des équations aux dérivées partielles, il existe tout de même quelques méthodes
permettant, dans certains cas spécifiques, de déterminer la solution d’un problème consti-
tué d’une équation aux dérivées partielles complétée par des conditions limites et/ou des
conditions initiales.

L’une d’entre elle (dite par séparation des variable), basée sur la recherche de solutions
de la forme
u(x, t) = X(x).T (t)
permet par exemple de résoudre l’équation de la chaleur 1D instationnaire sous des condi-
tions limites linéaires et homogènes, i.e. de la forme
[Link] (0, t) + b.u(0, t) = 0 et [Link] (1, t) + d.u(1, t) = 0, a, b, c, d ∈ R
Ainsi, considérons le problème de diffusion de chaleur 1D suivant :

 (E) : ut = α2 uxx ∀x ∈]0, 1[, ∀t > 0
(C.L.) : u(0, t) = u(1, t) = 0 ∀t > 0
(C.I.) : u(x, 0) = Φ(x) ∀x ∈ [0, 1]


Φ étant une fonction supposée connue, donnant la distribution de chaleur le long de la

barre à t = 0.

1. Montrer que si u est une solution de (E) de la forme u : (x, t) 7→ X(x).T (t), alors
les fonctions X et T vérifient des équations différentielles (ordinaires !) de la forme
(EX ) : X 00 − kX = 0 et (ET ) : T 0 − k.α2 .T = 0
pour une certaine constante k.
2. On pose k = −λ2 < 0. Résoudre les équations (EX ) et (ET ).
3. En déduire une expression de u(x, t) en fonction de λ, α et deux constantes A, B ∈
R.
4. Montrer que pour tout n ∈ N, il existe une fonction un de la forme
2 π 2 α2 t
un (x, t) = e−n .An . sin(nπx)
vérifiant à la fois l’équation (E) et les conditions initiales (C.L.).
5. Montrer que la fonction
+∞
X
u : (x, t) 7−→ un (x, t)
n=1

est également une solution de (E) vérifiant les conditions initiales (C.L.).
6. En exploitant la condition limite (C.I.) établir un lien entre les coefficients An
obtenus à la question précédente et la fonction Φ.
7. Détermine les coefficients An dans le cas ou Φ(x) = sin(πx) pour tout x ∈ [0, 1].
8. Même chose pour Φ(x) = 1 pour tout x ∈ [0, 1].