RESOLUTION DES EQUATIONS DIFFERENTIELLES

ORDINAIRES

V.1 Problème de Cauchy

On limite notre cours au cas d’équations différentielles ordinaires de la forme :
𝐓𝐫𝐨𝐮𝐯𝐞𝐫 𝐥𝐚 𝐟𝐨𝐧𝐜𝐭𝐢𝐨𝐧 𝒚 𝒕 𝐭𝐞𝐥𝐥𝐞 𝐪𝐮𝐞:
𝒚′ 𝒕 = 𝒇 𝒕, 𝒚 𝒕 ∀ 𝒕 ∈ 𝒂, 𝒃 (5.1)
𝒚 𝒕𝟎 = 𝒚 𝟎
Ce type de problème est appelé problème de Cauchy. Si la fonction 𝑓 𝑡, 𝑦 est continue et si
elle vérifie la condition de Lipschitz suivante :
∃L > 0 tel que ∀t ∈ [a, b] et ∀𝑦1 , 𝑦2 on ait : 𝑓 𝑡, 𝑦1 − 𝑓(𝑡, 𝑦2 ) ≤ 𝐿 𝑦1 − 𝑦2 (5.2)
alors le problème admet une solution unique pour toute valeur initiale.

V.2 Principe général des méthodes numériques

La solution exacte d’une équation différentielle du type (1) est une fonction continue. Les
ordinateurs ne peuvent fournir qu’un nombre fini de résultats numériques. Tout commence
donc par un choix préalable d’un nombre fini de points 𝑡𝑖 sur [a, b]. Ceci s’appelle une
discrétisation ou un maillage du domaine géométrique (ici le segment [a, b]). On limitera le
calcul au calcul approché de la solution en ces points.
Le choix des points 𝑡𝑖 est évidemment crucial pour la qualité de la solution numérique
obtenue. Le maillage doit permettre de représenter de façon précise la solution. Comme cette
solution est inconnue, on procède par des techniques d’adaptation de maillage à posteriori. On
calcule une première solution sur un premier maillage. On déduit de ce premier calcul, les
zones de fortes variations de la solution. On raffine par la suite le maillage dans ces zones. Par
souci de simplicité, on présente ici les méthodes numériques dans le cas de maillage à pas
uniformes. Si on partage l’intervalle 𝑎, 𝑏 en 𝑚 segments, le pas de discrétisation ℎ est égal :
𝑏−𝑎
ℎ= (5.3)
𝑚

On pose :
𝒕𝒏+𝟏 = 𝒕𝒏 + 𝒉 (𝑛 = 0,1, … . , 𝑚) avec 𝑡𝑜 = 𝑎
On peut construire des schémas d’intégration d’équations différentielles de diverses
manières :
– Soit en utilisant des formules d’intégration numérique approchée, du chapitre précédent.
– Soit en utilisant le développement de Taylor.

Deux critères principaux gouvernent le choix d’une méthode :

– l’ordre, qui mesure la précision du schéma.
– la stabilité, qui concerne le comportement de la solution approchée. Le schéma est stable si
la solution discrète reste bornée quel que soit la taille du pas de discrétisation.

V.3 Méthodes à un pas

Dans ces méthodes la valeur approchée 𝑦𝑛+1 de la fonction inconnue y pour 𝑡𝑛+1 est calculée
en fonction des seules valeurs de l’abscisse précédente 𝑡𝑛 , de l’approximation 𝑦𝑛 et du pas de
discrétisation h.
Si 𝑦𝑛 +1 s’obtient par une formule explicite de la forme :
𝑦𝑛+1 = 𝑦𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 , ℎ (5.4)
on dit que la méthode est explicite.
Si par contre 𝑦𝑛+1 est donnée par une relation de la forme :
𝑦𝑛+1 = 𝑦𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 , 𝑦𝑛 +1 , ℎ (5.5)
on dit que la méthode est implicite.

Ces schémas sont obtenus, par exemple, en intégrant l’équation différentielle et en utilisant
des formules d’intégration numérique pour le second membre. L’ordre du schéma sera égal au
degré du polynôme pour lequel l’intégration est exacte +1.
A titre d’exemple, on obtient les schémas suivants:
𝑡 𝑛 +1 𝑡 𝑛 +1
′
𝑦 𝑡 𝑑𝑡 = 𝑦 𝑡𝑛+1 − 𝑦 𝑡𝑛 = 𝑓 𝑡, 𝑦(𝑡) 𝑑𝑡
𝑡𝑛 𝑡𝑛
𝑡 𝑛 +1
Pour calculer 𝑡𝑛
𝑓 𝑡, 𝑦(𝑡) 𝑑𝑡, on utilise les méthodes d’intégration décrites dans le chapitre

précédent.
En choisissant la méthode des rectangles (point à gauche), on obtient :
𝑦𝑛+1 − 𝑦𝑛 = ℎ𝑓 𝑡𝑛 , 𝑦𝑛
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉𝒇 𝒕𝒏 , 𝒚𝒏 (5.6)
On obtient le schéma d’Euler explicite, ce schéma est d’ordre 1, car on a remplacé, la
fonction 𝑓 par une Constante (polynôme de degré 0).
En utilisant toujours la méthode des rectangles mais en choisissant le point à droite on
obtient:
𝑦𝑛+1 − 𝑦𝑛 = ℎ𝑓 𝑡𝑛+1 , 𝑦𝑛+1
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉𝒇 𝒕𝒏+𝟏 , 𝒚𝒏+𝟏 (5.7)
On obtient le schéma d’Euler implicite, ce schéma est d’ordre 1.

En utilisant toujours la méthode des rectangles mais en choisissant le point milieu

𝑡 𝑛 +𝑡 𝑛 +1 ℎ
= 𝑡𝑛 + 2 on obtient:
2

ℎ ℎ
𝑦𝑛 +1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 + , 𝑦 𝑡𝑛 +
2 2
ℎ
𝑦 𝑡𝑛 + 2 est calculé à l’aide du schéma d’Euler explicite :
ℎ ℎ 𝑘
𝑦 𝑡𝑛 + = 𝑦 𝑡𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 = 𝑦 𝑡𝑛 +
2 2 2
Ce qui donne :
𝒌 = 𝒉𝒇 𝒕𝒏 , 𝒚𝒏
𝒉 𝒌 (5.8)
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉𝒇 𝒕𝒏 + 𝟐 , 𝒚𝒏 + 𝟐

On obtient le "Schéma de Runge-Kutta d’ordre 2", appelé aussi "Schéma du point

milieu". C’est un schéma explicite d’ordre 2.

En utilisant la méthode des trapèzes, on obtient :

ℎ
𝑦𝑛+1 − 𝑦𝑛 = 𝑓 𝑡𝑛 , 𝑦𝑛 + 𝑓 𝑡𝑛+1 , 𝑦𝑛+1
2
𝒉
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝟐 𝒇 𝒕𝒏 , 𝒚𝒏 + 𝒇 𝒕𝒏+𝟏 , 𝒚𝒏+𝟏 (5.9)

On obtient le schéma de Crank-Nicolson, c’est un schéma d’ordre 2.

Le schéma de Crank-Nicholson est implicite car le calcul de la nouvelle valeur 𝑦𝑛+1 nécessite
la résolution d’une équation algébrique. Si l’on veut obtenir une méthode explicite du même
ordre, on peut procéder de la manière suivante :
𝒚∗𝒏+𝟏 = 𝒚𝒏 + 𝒉𝒇 𝒕𝒏 , 𝒚𝒏
𝒉 (5.10)
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝟐 𝒇 𝒕𝒏 , 𝒚𝒏 + 𝒇 𝒕𝒏+𝟏 , 𝒚∗𝒏+𝟏

On obtient le "Schéma de Runge-Kutta d’ordre 2", appelé aussi "Schéma d’Euler

modifié". C’est un schéma explicite d’ordre 2. Ce schéma est différent de celui du point
milieu mais il est très proche.
De même, l’utilisation de la formule d’intégration de Simpson est à la base de la formule de
Runge et Kutta d’ordre 4.
En général, un schéma explicite est facile à utiliser mais impose, pour sa stabilité, des valeurs
très petites du pas ℎ.
On peut obtenir les schémas précédents en utilisant le développement de Taylor d’ordre 1 :
𝑦 𝑡𝑛+1 = 𝑦 𝑡𝑛 + ℎ = 𝑦 𝑡𝑛 + ℎ𝑦 ′ 𝜏 𝜏 ∈ [𝑡𝑛 , 𝑡𝑛+1 ]
Si on prend 𝜏 = 𝑡𝑛 , on obtient :
𝑦 𝑡𝑛+1 = 𝑦 𝑡𝑛 + ℎ𝑦 ′ 𝑡𝑛
Sachant que 𝑦 ′ 𝑡𝑛 = 𝑓 𝑡𝑛 , 𝑦𝑛 , on obtient le schéma d’Euler Explicite.
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉𝒇 𝒕𝒏 , 𝒚𝒏
Géométriquement, le schéma d’Euler Explicite revient à remplacer localement en chaque
point 𝑡𝑛 , la courbe solution par sa tangente (Figure 22). On voit donc qu’au cours du
processus numérique, on va passer à chaque pas, d’une courbe solution à une courbe voisine
correspondant à une condition initiale légèrement différente. Si le problème est stable, on
pourra obtenir la convergence.

Figure 22 : Interprétation de la méthode d’Euler Explicite

Si on prend 𝜏 = 𝑡𝑛 +1 , on obtient :
𝑦 𝑡𝑛+1 = 𝑦 𝑡𝑛 + ℎ𝑦 ′ 𝑡𝑛+1

Sachant que 𝑦 ′ 𝑡𝑛+1 = 𝑓 𝑡𝑛+1 , 𝑦𝑛+1 , on obtient le schéma d’Euler Implicite.

𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉𝒇 𝒕𝒏+𝟏 , 𝒚𝒏+𝟏
On peut construire un schéma intermédiaire entre le schéma d’Euler Explicite et d’Euler
Implicite de la façon suivante :

𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉 𝟏 − 𝜶 𝒇 𝒕𝒏 , 𝒚𝒏 + 𝜶𝒇 𝒕𝒏+𝟏 , 𝒚𝒏+𝟏 (5.11)

Si 𝛼 = 0 on obtient le schéma d’Euler Explicite

Si 𝛼 = 1 on obtient le schéma d’Euler Implicite
Si 𝛼 = 0.5 on obtient le schéma implicite de Crank-Nicholson

𝑡 𝑛 +𝑡 𝑛 +1 ℎ
Si on choisit 𝜏 = = 𝑡𝑛 + 2 , on obtient le Schéma du Point Milieu :
2

ℎ
𝑦 𝑡𝑛+1 = 𝑦 𝑡𝑛 + ℎ𝑦 ′ 𝑡𝑛 +
2
ℎ ℎ ℎ
𝑦 ′ 𝑡𝑛 + = 𝑓 𝑡𝑛 + , 𝑦 𝑡𝑛 +
2 2 2
ℎ
𝑦 𝑡𝑛 + 2 est calculé à l’aide du schéma d’Euler explicite :
ℎ ℎ 𝑘
𝑦 𝑡𝑛 + = 𝑦 𝑡𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 = 𝑦 𝑡𝑛 +
2 2 2
avec
𝒌 = 𝒉𝒇 𝒕𝒏 , 𝒚𝒏
𝒉 𝒌
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉𝒇 𝒕𝒏 + , 𝒚𝒏 +
𝟐 𝟐

En utilisant le développement de Taylor à l’ordre 2, on obtient:

ℎ2 ′′ ′
𝑦 𝑡𝑛+1 = 𝑦 𝑡𝑛 + ℎ𝑦 𝑡𝑛 + 𝑦 𝜏 𝜏 ∈ [𝑡𝑛 , 𝑡𝑛 +1 ]
2
𝑑𝑦′(𝑡) 𝑑𝑓(𝑡, 𝑦) 𝜕𝑓(𝑡, 𝑦) 𝜕𝑓 𝑡, 𝑦
𝑦 ′′ 𝑡 = = = + 𝑓(𝑡, 𝑦)
𝑑𝑡 𝑑𝑡 𝜕𝑡 𝜕𝑦
En choisissant 𝜏 = 𝑡𝑛 , on obtient le schéma explicite de Taylor d’ordre 2 :
𝒉𝟐 𝝏𝒇 𝒕𝒏 , 𝒚𝒏 𝝏𝒇 𝒕𝒏 , 𝒚𝒏
𝒚𝒏+𝟏 = 𝒚𝒏 + 𝒉𝒇 𝒕𝒏 , 𝒚𝒏 + + 𝒇 𝒕𝒏 , 𝒚 𝒏
𝟐 𝝏𝒕 𝝏𝒚

L’inconvénient de ce schéma est qu’il exige la dérivabilité de la fonction 𝑓(𝑡, 𝑦).

V.4 Algorithmes des différents schémas

V.4.1 Schémas Explicites

a- Euler Explicite (d’ordre 1) :
𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛+1 = 𝑡𝑛 + ℎ (5.12)
𝑦𝑛 +1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 , 𝑦𝑛
b- Runge-Kutta (d’ordre 2) : Il existe deux versions pour ce schéma.
Schéma de Lax-Wendroff ou Schéma du point milieu
𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑘 = ℎ𝑓 𝑡𝑛 , 𝑦𝑛 (5.13)
ℎ 𝑘
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 + 2 , 𝑦𝑛 + 2

Schéma d’Euler modifié

𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛+1 = 𝑡𝑛 + ℎ
∗
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 , 𝑦𝑛
ℎ ∗
𝑦𝑛+1 = 𝑦𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 + 𝑓 𝑡𝑛+1 , 𝑦𝑛+1
2
Schéma de Taylor d’ordre 2
𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛+1 = 𝑡𝑛 + ℎ
ℎ2 𝜕𝑓 𝑡𝑛 , 𝑦𝑛 𝜕𝑓 𝑡𝑛 , 𝑦𝑛
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 , 𝑦𝑛 + + 𝑓 𝑡𝑛 , 𝑦𝑛
2 𝜕𝑡 𝜕𝑦

c- Runge-Kutta (d’ordre 4) :
Ce schéma est très précis, il est d’ordre 4.
𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛 +1 = 𝑡𝑛 + ℎ
𝑘1 = ℎ𝑓 𝑡𝑛 , 𝑦𝑛
ℎ 𝑘1
𝑘2 = ℎ𝑓 𝑡𝑛 + 2 , 𝑦𝑛 + 2 (5.14)
ℎ 𝑘2
𝑘3 = ℎ𝑓 𝑡𝑛 + 2 , 𝑦𝑛 + 2
𝑘4 = ℎ𝑓 𝑡𝑛 + ℎ, 𝑦𝑛 + 𝑘3
1
𝑦𝑛 +1 = 𝑦𝑛 + 6 𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4

V.4.2. Schémas Implicites

a- Euler Implicite (d’ordre 1) :
 𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛+1 = 𝑡𝑛 + ℎ (5.15)
𝑦𝑛 +1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛+1 , 𝑦𝑛+1

b- Crank-Nichoson (d’ordre 2)
𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛+1 = 𝑡𝑛 + ℎ (5.16)
ℎ
𝑦𝑛 +1 = 𝑦𝑛 + 2 𝑓 𝑡𝑛 , 𝑦𝑛 + 𝑓 𝑡𝑛+1 , 𝑦𝑛+1

V.5 Exemples et programmes

V.5.1 Exemple 1
Soit à résoudre l’équation différentielle :
𝑦 0 =1
′
𝑦 𝑡 = −𝑦 𝑡 + 𝑡 + 1 𝑡 ∈ [0,1]
1- Calculer la solution analytique.
2- Calculer une valeur approchée en utilisant les schémas décrits précédemment (prendre
h=0.2).

Solution
1- Solution exacte
L’équation à résoudre est une équation différentielle linéaire de 1 er ordre, sa solution
analytique est :
𝑦 𝑡 = 𝑒 −𝑡 + 𝑡
Les valeurs de y pour 𝑡 ∈ [0,1] avec ℎ = 0.2 sont données par le tableau suivant :

𝑡 0 0.2 0.4 0.6 0.8 1

𝑦 1.0000 1.0187 1.0703 1.1488 1.2493 1.3679

2- Solution approchée
On a 𝑡0 = 0, 𝑦0 = 1, ℎ = 0.2 et 𝑓 𝑡, 𝑦 = −𝑦 + 𝑡 + 1
V.5.1.1 Schémas Explicites
Euler explicite
𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑦𝑛 +1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 , 𝑦𝑛 = 𝑦𝑛 + ℎ −𝑦𝑛 + 𝑡𝑛 + 1
𝑦1 = 𝑦0 + ℎ −𝑦0 + 𝑡0 + 1 = 1 + 0.2 −1 + 0 + 1 = 1

Attention : 𝑦1 ≠ 𝑦(1). 𝑦1 = 𝑦 𝑡1 = 𝑦 𝑡0 + ℎ = 𝑦(0.2)

Les valeurs de y pour 𝑡 ∈ [0,1] avec ℎ = 0.2 sont données par le tableau suivant :
𝑡 0 0.2 0.4 0.6 0.8 1
𝑦 1.0000 1.0000 1.0400 1.1120 1.2096 1.3277

Point Milieu
𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑘 = ℎ𝑓 𝑡𝑛 , 𝑦𝑛 = ℎ −𝑦𝑛 + 𝑡𝑛 + 1
ℎ 𝑘 𝑘 ℎ
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 + , 𝑦𝑛 + = 𝑦𝑛 + ℎ − 𝑦𝑛 + + 𝑡𝑛 + +1
2 2 2 2
Détail de calcul pour 𝑦1
𝑘 = ℎ −𝑦0 + 𝑡0 + 1 = 0.1 −1 + 0 + 1 = 0
𝑘 ℎ
𝑦1 = 𝑦0 + ℎ − 𝑦0 + + 𝑡0 + + 1 = 1.0200
2 2
Les valeurs de y pour 𝑡 ∈ [0,2] avec ℎ = 0.2 sont données par le tableau suivant :
𝑡 0 0.2 0.4 0.6 0.8 1
𝑦 1.0000 1.0200 1.0724 1.1514 1.2521 1.3707

Euler modifié
𝑦0 = 1
𝑡𝑛 +1 = 𝑡𝑛 + ℎ
∗
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 , 𝑦𝑛 = 𝑦𝑛 + ℎ −𝑦𝑛 + 𝑡𝑛 + 1
ℎ ∗
ℎ ∗
𝑦𝑛 +1 = 𝑦𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 + 𝑓 𝑡𝑛+1 , 𝑦𝑛+1 = 𝑦𝑛 + −𝑦𝑛 + 𝑡𝑛 + 1 + −𝑦𝑛+1 + 𝑡𝑛 +1 + 1
2 2
Détail de calcul pour 𝑦1
𝑦1∗ = 𝑦0 + ℎ −𝑦0 + 𝑡0 + 1 = 1
ℎ
𝑦1 = 𝑦0 + −𝑦0 + 𝑡0 + 1 + −𝑦1∗ + 𝑡1 + 1 = 1.0200
2
Les valeurs de y pour 𝑡 ∈ [0,2] avec ℎ = 0.2 sont données par le tableau suivant :
 𝑡 0 0.2 0.4 0.6 0.8 1
𝑦 1.0000 1.0200 1.0724 1.1514 1.2521 1.3707

Runge-Kutta d’ordre 4
𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑘1 = ℎ𝑓 𝑡𝑛 , 𝑦𝑛 = ℎ −𝑦𝑛 + 𝑡𝑛 + 1
ℎ 𝑘1 𝑘1 ℎ
𝑘2 = ℎ𝑓 𝑡𝑛 + , 𝑦𝑛 + = ℎ − 𝑦𝑛 + + 𝑡𝑛 + +1
2 2 2 2
ℎ 𝑘2 𝑘2 ℎ
𝑘3 = ℎ𝑓 𝑡𝑛 + , 𝑦𝑛 + = ℎ − 𝑦𝑛 + + 𝑡𝑛 + +1
2 2 2 2
𝑘4 = ℎ𝑓 𝑡𝑛 + ℎ, 𝑦𝑛 + 𝑘3 = ℎ − 𝑦𝑛 + 𝑘3 + 𝑡𝑛 + ℎ + 1
1
𝑦𝑛 +1 = 𝑦𝑛 + 𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4
6
Détail de calcul pour 𝑦1
𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ = 0.2
𝑘1 = ℎ −𝑦0 + 𝑡0 + 1
𝑘1 ℎ
𝑘2 = ℎ − 𝑦0 + + 𝑡0 + + 1 = 0.02
2 2
𝑘2 ℎ
𝑘3 = ℎ − 𝑦0 + + 𝑡0 + + 1 = 0.018
2 2
𝑘4 = ℎ − 𝑦𝑛 + 𝑘3 + 𝑡𝑛 + ℎ + 1 = 0.0364
1
𝑦1 = 𝑦0 + 𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 = 1.0187
6
Les valeurs de y pour 𝑡 ∈ [0,2] avec ℎ = 0.2 sont données par le tableau suivant :
𝑡 0 0.2 0.4 0.6 0.8 1
𝑦 1.0000 1.0187 1.0703 1.1488 1.2493 1.3679

V.5.1.2 Schémas Implicites

Euler implicite
𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛+1 , 𝑦𝑛+1 = 𝑦𝑛 + ℎ −𝑦𝑛+1 + 𝑡𝑛+1 + 1
𝑦1 = 𝑦0 + ℎ −𝑦0 + 𝑡0 + 1 = 1 + 0.1 −1 + 0 + 1 = 1
Pour un schéma implicite, le calcul de 𝑦𝑛+1 nécessite la résolution de l’équation :
𝑦𝑛+1 = 𝑦𝑛 + ℎ −𝑦𝑛+1 + 𝑡𝑛+1 + 1
Ce qui donne :
 𝑦𝑛 + ℎ 𝑡𝑛 + ℎ + 1
𝑦𝑛 +1 =
1+ℎ
𝑦0 + ℎ 𝑡0 + ℎ + 1 1 + 0.1 0 + 0.1 + 1
𝑦1 = = = 1.0333
1+ℎ 1 + 0.1
Les valeurs de y pour 𝑡 ∈ [0,2] avec ℎ = 0.2 sont données par le tableau suivant :

𝑡 0 0.2 0.4 0.6 0.8 1

𝑦 1.0000 1.0333 1.0944 1.1787 1.2823 1.4019

Crank-Nicholson
𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
ℎ ℎ
𝑦𝑛+1 = 𝑦𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 + 𝑓 𝑡𝑛+1 , 𝑦𝑛+1 = 𝑦𝑛 + −𝑦𝑛 + 𝑡𝑛 + 1 + −𝑦𝑛+1 + 𝑡𝑛 +1 + 1
2 2
ℎ
𝑦𝑛 +1 = 𝑦𝑛 + −𝑦𝑛 + 𝑡𝑛 + 1 + −𝑦𝑛+1 + 𝑡𝑛+1 + 1
2
ℎ ℎ
1 − 2 𝑦𝑛 + 2 2𝑡𝑛 + ℎ + 2
𝑦𝑛+1 =
ℎ
1+2
ℎ 0.1
1 − 2 𝑦0 + ℎ 2𝑡0 + ℎ + 2 1 − 2 1 + 0.1 0 + 0.1 + 2
𝑦1 = = = 1.0182
ℎ 0.1
1+2 1+ 2

Les valeurs de y pour 𝑡 ∈ [0,2] avec ℎ = 0.2 sont données par le tableau suivant :
𝑡 0 0.2 0.4 0.6 0.8 1
𝑦 1.0000 1.0182 1.0694 1.1477 1.2481 1.3666

Exemple 2
𝜋
Résoudre dans l’intervalle [0, 𝜋] l’équation différentielle (prendre ℎ = 10) :

𝑦 0 =1
′ 2
𝑦 𝑡 = 𝑡 sin⁡
(𝑦) 𝑡 ∈ [0, 𝜋]

Euler explicite
𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 , 𝑦𝑛 = 𝑦𝑛 + ℎ𝑡𝑛2 sin⁡
(𝑦𝑛 )
 𝜋
Les valeurs de y pour 𝑡 ∈ [0, 𝜋] avec ℎ = 10 sont données par le tableau suivant :

t 0 𝜋 2𝜋 3𝜋 4𝜋 5𝜋 6𝜋 7𝜋 8𝜋 9𝜋 𝜋
10 10 10 10 10 10 10 10 10
y 1. 1. 1.0261 1.1322 1.3848 1.8724 2.6125 3.1759 3.1238 3.1591 3.1151

Lax-Wendroff ou Point Milieu

𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑘 = ℎ𝑓 𝑡𝑛 , 𝑦𝑛 = 𝑦𝑛 + ℎ𝑡𝑛2 sin⁡(𝑦𝑛 )
2
ℎ 𝑘 ℎ 𝑘
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 + , 𝑦𝑛 + = 𝑦𝑛 + ℎ 𝑡𝑛 + sin⁡
(𝑦𝑛 + )
2 2 2 2
𝜋
Les valeurs de y pour 𝑡 ∈ [0, 𝜋] avec ℎ = 10 sont données par le tableau suivant :

t 0 𝜋 2𝜋 3𝜋 4𝜋 5𝜋 6𝜋 7𝜋 8𝜋 9𝜋 𝜋
10 10 10 10 10 10 10 10 10
y 1. 1.0065 1.066 1.2404 1.6128 2.2145 2.7572 2.9854 3.0517 3.0536 2.9910

Euler modifié

𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
∗
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛 , 𝑦𝑛 = 𝑦𝑛 + ℎ𝑡𝑛2 sin⁡
(𝑦𝑛 )
ℎ ℎ
𝑦𝑛+1 = 𝑦𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 + 𝑓 𝑡𝑛+1 , 𝑦𝑛∗ +1 = 𝑦𝑛 + 𝑡𝑛2 sin⁡(𝑦𝑛 ) + 𝑡𝑛 + ℎ 2 sin⁡ ∗
(𝑦𝑛+1 )
2 2
𝜋
Les valeurs de y pour 𝑡 ∈ [0, 𝜋] avec ℎ = 10 sont données par le tableau suivant :

t 0 𝜋 2𝜋 3𝜋 4𝜋 5𝜋 6𝜋 7𝜋 8𝜋 9𝜋 𝜋
10 10 10 10 10 10 10 10 10
y 1. 1.0065 1.066 1.2404 1.6128 2.2145 2.7572 2.9854 3.0517 3.0536 2.9910

Runge-Kutta d’ordre 4
 𝑦0 𝑑𝑜𝑛𝑛é
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑘1 = ℎ𝑓 𝑡𝑛 , 𝑦𝑛 = ℎ 𝑡𝑛2 sin⁡
(𝑦𝑛 )
ℎ 𝑘1 ℎ 2 𝑘1
𝑘2 = ℎ𝑓 𝑡𝑛 + , 𝑦𝑛 + = ℎ 𝑡𝑛 + sin⁡ 𝑦𝑛 +
2 2 2 2
2
ℎ 𝑘2 ℎ 𝑘2
𝑘3 = ℎ𝑓 𝑡𝑛 + , 𝑦𝑛 + = ℎ 𝑡𝑛 + sin⁡ 𝑦𝑛 +
2 2 2 2
2
𝑘4 = ℎ𝑓 𝑡𝑛 + ℎ, 𝑦𝑛 + 𝑘3 = ℎ 𝑡𝑛 + ℎ sin⁡𝑦𝑛 + 𝑘3
1
𝑦𝑛+1 = 𝑦𝑛 + 𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4
6

𝜋
Les valeurs de y pour 𝑡 ∈ [0, 𝜋] avec ℎ = 10 sont données par le tableau suivant :

t 0 𝜋 2𝜋 3𝜋 4𝜋 5𝜋 6𝜋 7𝜋 8𝜋 9𝜋 𝜋
10 10 10 10 10 10 10 10 10
y 1. 1.0087 1.0711 1.2508 1.6274 2.2090 2.7515 3.0267 3.1103 3.128 3.128

L’avantage des schémas explicites est le fait qu’ils soient très faciles à utiliser, la valeur de
𝑦𝑛+1 dépend uniquement des valeurs au pas précédent. Pour calculer les valeurs de y, on a
utilisé le même programme que précédemment, on a juste changé, les données et le fichier
suivant :
% fichier functest.m
function y=functest(a,b)
y=b^2*sin(a) ;

Comme la solution exacte n’existe pas, dans la figure 23 on a tracé les solutions données par
les différents schémas explicites
 Figure 23 : Résultats des schémas explicites

Euler implicite

𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑓 𝑡𝑛+1 , 𝑦𝑛+1 = 𝑦𝑛 + ℎ 𝑡𝑛 + ℎ 2 sin⁡𝑦𝑛+1
Pour calculer 𝑦𝑛+1 , il faut résoudre une équation non linéaire :
𝑦𝑛+1 = 𝑦𝑛 + ℎ 𝑡𝑛 + ℎ 2 sin⁡𝑦𝑛+1
𝑦𝑛 , 𝑡𝑛 𝑒𝑡 ℎ sont connus, si on pose 𝑥 = 𝑦𝑛+1 , l’équation précédente devient :
𝐹 𝑥 = 𝐴 sin 𝑥 − 𝑥 + 𝐵 = 0 𝐴 = ℎ 𝑡𝑛 + ℎ 2 , 𝐵 = 𝑦𝑛
Qu’on résolve par la méthode de Newton-Raphson :
𝑥0 = 𝑦𝑛
𝑓(𝑥𝑛 )
𝑥𝑛+1 = 𝑥𝑛 −
𝑓′(𝑥𝑛 )
𝑐𝑟𝑖𝑡è𝑟𝑒 𝑑′𝑎𝑟𝑟ê𝑡
Pour chaque valeur de 𝑦𝑛 , il faut résoudre une équation non linéaire. C’est l’inconvénient
majeur des schémas implicites.
𝜋
Les valeurs de y pour 𝑡 ∈ [0, 𝜋] avec ℎ = 10 sont données par le tableau suivant :
t 0 𝜋 2𝜋 3𝜋 4𝜋 5𝜋 6𝜋 7𝜋 8𝜋 9𝜋 𝜋
10 10 10 10 10 10 10 10 10
y 1. 1.0265 1.1392 1.4148 1.8864 2.4064 2.7904 3.0019 3.0948 3.1283 3.1383

Crank-Nicholson
𝑦0 = 1
𝑡𝑛+1 = 𝑡𝑛 + ℎ
ℎ ℎ 2
𝑦𝑛+1 = 𝑦𝑛 + 𝑓 𝑡𝑛 , 𝑦𝑛 + 𝑓 𝑡𝑛+1 , 𝑦𝑛+1 = 𝑦𝑛 + 𝑡 sin⁡𝑦𝑛 + 𝑡𝑛 + ℎ 2 sin⁡𝑦𝑛+1
2 2 𝑛

Pour calculer 𝑦𝑛+1 , il faut résoudre une équation non linéaire :

ℎ ℎ
𝐹 𝑥 = 𝐴 sin 𝑥 − 𝑥 + 𝐵 = 0; 𝑥 = 𝑦𝑛+1 , 𝐴 = 𝑡𝑛 + ℎ 2 , 𝐵 = 𝑦𝑛 + 𝑡𝑛 2 sin⁡𝑦𝑛
2 2

𝜋
Les valeurs de y pour 𝑡 ∈ [0, 𝜋] avec ℎ = 10 sont données par le tableau suivant :

t 0 𝜋 2𝜋 3𝜋 4𝜋 5𝜋 6𝜋 7𝜋 8𝜋 9𝜋 𝜋
10 10 10 10 10 10 10 10 10
y 1. 1.0132 1.0810 1.2690 1.6495 2.2082 2.7386 3.0369 3.1289 3.1415 3.1416

Sur la figure 24, on a tracé les résultats des deux schémas implicites
Figure 24 : Résultats des schémas explicites