1. No category

MPSI2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Table des mati`
eres
1 Sch´
ema d’Euler
1.1 Position du probl`eme et exemples . .
1.1.1 Radioactivit´e . . . . . . . . .
1.1.2 Croissance logistique . . . . .
1.1.3 Chute libre . . . . . . . . . .
1.2 Principe de la m´ethode d’Euler . . .
´
1.3 Ecriture
du programme Python . . .
1.4 Test du programme sur les exemples
1.4.1 Radioactivit´e . . . . . . . . .
1.4.2 Croissance logistique . . . . .
1.4.3 Chute libre . . . . . . . . . .
1.5 Deux exemples probl´ematiques . . .
1.6 Qualit´e d’un sch´ema num´erique . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
1.6.1 Consistance . . . . . . . . . . . . .
1.6.2 Stabilit´e . . . . . . . . . . . . . . .
1
1.6.3 Convergence . . . . . . . . . . . .
1
1.6.4 Erreurs d’arrondi . . . . . . . . . .
2
2
2 Equations d’ordre 2 ou plus
2
2.1 Vectorisation du probl`eme . . . . . . . . .
4
2.2 Adaptation du programme . . . . . . . . .
5
2.3 Exemple de l’´equation du pendule simple
5
2.4 Syst`eme proies pr´edateurs . . . . . . . . .
6
7 3 Autres sch´
emas num´
eriques
7
3.1 M´ethode d’Euler implicite . . . . . . . . .
3.2 Sch´ema de Runge-Kutta d’ordre 4 : RK4 .
9
.
.
.
.
.
.
.
.
.
.
.
.
9
9
9
9
.
.
.
.
.
.
.
.
.
.
.
.
10
10
10
11
13
15
. . . 15
. . . 16
Hormis quelques cas d’´ecole simples, on ne sait pas d´eterminer d’expressions analytiques pour les solutions d’´equations
diff´erentielles.
• L’´etude du mouvement d’un pendule simple sans frottement nous am`ene `a l’´equation :
θ¨ + ω 2 sin(θ) = 0
• Les probl`emes de cin´
etique chimique conduisent `a des syst`emes diff´erentiels non lin´eaires d´ecrivant les ´evolutions des
diff´erents r´eactifs au cours du temps.
• Les ph´enom`enes que l’on observe en m´
ecanique des fluides sont en partie d´ecrits par les ´equations non lin´eaires aux
d´eriv´ees partielles de Navier-Stokes.
Le but de ce chapitre est de pr´esenter des m´ethodes qui permettent d’avoir des solutions approch´ees d’´equations
diff´erentielles.
1
Sch´
ema d’Euler
1.1
Position du probl`
eme et exemples
La mod´elisation d’un grand nombre de probl`emes ayant leur origine en g´eom´etrie, m´ecanique, physique, sciences de
l’ing´enieur, chimie, biologie, ´economie conduit `
a chercher les solutions de syst`emes de la forme :
0
y (t) = F (t, y(t))
(C) :
y(t0 ) = y0
La fonction y est celle que l’on cherche, elle est de classe C 1 sur un intervalle [a, b]. La condition y(t0 ) = y0 est appel´ee
condition initiale, on a t0 ∈ [a, b]. Un tel probl`eme est appel´e probl`
eme de Cauchy.
Le th´eor`eme de Cauchy-Lipschitz garantit l’existence et l’unicit´e d’une solution au probl`eme de Cauchy si l’on ajoute
des conditions sur la fonction F .
Afin d’appr´ehender la notation y 0 (t) = F (t, y(t)), voyons quelques exemples.
1.1.1
Radioactivit´
e
On observe que si N (t) d´esigne le nombre d’atomes radioactifs pr´esents dans un ´echantillon de mati`ere donn´e, le nombre
des d´esint´egrations d’atomes de ce type est proportionnel `a N (t) et `a la dur´ee de l’observation. C’est-`a-dire, avec λ une
constante r´eelle :
N (t1 ) − N (t0 ) = −λN (t0 )(t1 − t0 )
On passe `
a la limite quand t1 tend vers t0 apr`es avoir divis´e par t1 − t0 pour obtenir :
N 0 (t) = −λN (t)
Cette ´equation rentre dans le cadre de la forme annonc´ee avec F (t, N ) = −λN . On remarque d’ailleurs que la fonction F
ne d´epend pas de t, ce sera souvent le cas dans les exemples que nous allons prendre. On dit que l’´equation diff´erentielle est
autonome.
1
Chapitre 11
R´
esolution num´
erique d’´
equations diff´
erentielles
MPSI2
1.1.2
Informatique
Croissance logistique
En dynamique des populations, un mod`ele de croissance d´ecrit la fa¸con dont une population ´evolue dans le temps. Le
mod`
ele de croissance logistique prend en compte l’environnement et l’auto-r´egulation des populations dans l’expression
du taux de croissance, ce qui donne :
P (t) P (t)
P 0 (t) = r 1 −
K
• P (t) est le nombre d’individus en fonction du temps.
• r est le taux d’accroissement intrins`eque de la population.
• K est la capacit´e d’accueil du milieu.
On remarque que P 0 (t) < 0 lorsque P (t) > K et P 0 (t) > 0 lorsque P (t) < K. L`a encore, cela rentre dans le cadre
P
th´eorique pr´esent´e avec F (t, P ) = r 1 −
P.
K
1.1.3
Chute libre
Un objet en chute libre dans l’atmosph`ere a sa vitesse qui v´erifie l’´equation diff´erentielle :
(
α
v(t)2 − g
m
= v0
v 0 (t)
=
v(t0 )
Dans cette ´equation, on a tenu compte de la r´esistance de l’air que nous supposons proportionnelle au carr´e de la vitesse. Le
coefficient α est proportionnel `
a la densit´e de l’air ρ, `a la surface de la section du corps en chute S et `a Cx : le coefficient de
1
r´esistance a´erodynamique qui d´epend de la forme. Ce qui donne α = SCx ρ.
2
α
Ici, on a : F (t, v) = v 2 − g.
m
1.2
Principe de la m´
ethode d’Euler
On reprend le probl`eme de Cauchy :
(C) :
y 0 (t) = F (t, y(t))
y(t0 ) = y0
i) On ne connaˆıt pas la solution exacte, y, cependant on dispose de deux informations : y(t0 ) = y0 et y 0 (t0 ) = F (t0 , y(t0 )).
ii) On se donne un pas h (qui a vocation `
a ˆetre proche de 0) et on pose t1 = t0 + h. On esp`ere qu’une valeur approch´ee
de y(t1 ) va ˆetre :
y1 = y0 + hy 0 (t0 ) = y0 + hF (t0 , y(t0 ))
C’est-`
a-dire l’on utilise l’approximation :
y(t0 + h) − y(t0 )
≈ y 0 (t0 )
h
Graphiquement, cela revient `
a approcher la courbe repr´esentative de y par sa tangente en t0 , cette approximation sera
d’autant plus pertinente que le pas sera proche de 0.
iii) On poursuit le proc´ed´e en partant de y1 . On pose t2 = t1 + h et on donne une valeur approch´ee de y(t2 ) de la
mˆeme fa¸con :
y2 = y1 + hy 0 (t1 ) = y1 + hF (t1 , y1 )
Plus g´en´eralement, on se donne la subdivision de [t0 , t0 + T ] d´efinie par :
∀n ∈ J0, N K, tn = t0 + n
Le pas h vaut
T
N
T
. Pour tout n ∈ J0, N − 1K :
N

 y0
tn+1

yn+1
=
=
=
y(t0 )
tn + h
yn + hF (tn , yn )
2
Chapitre 11
MPSI2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
La ligne bris´ee reliant les points {(tn , yn ), n ∈ J0, N K} donnera une solution approch´ee de notre ´equation diff´erentielle
avec condition initiale.
Remarque : Cette m´ethode revient `
a consid´erer que si h est petit :
y(tn + h) = y(tn ) + hy 0 (tn )
Exercice 1 : On s’int´eresse au probl`eme de Cauchy suivant dont on connaˆıt bien entendu la solution exacte :
0
y (t) = y(t)
y(0) = 1
On cherche `
a trouver une solution approch´ee sur l’intervalle [0, 1]. On subdivise l’intervalle [0, 1] avec un pas de h =
N ∈ N∗ .
´
1. Ecrire
la suite (yn ) donn´ee par la m´ethode d’Euler.
2. Exprimer yn en fonction de n et du pas h.
3. En examinant la valeur obtenue en t = 1 dire si l’approximation obtenue de la fonction semble satisfaisante ?
1
o`
u
N
Voici les graphiques obtenus avec N = 4 et N = 10 :
3
Chapitre 11
MPSI2
1.3
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
´
Ecriture
du programme Python
La fonction que l’on va ´ecrire prend en entr´ee 5 param`etres :
• La fonction F , qui d´epend des variables t et y, mˆeme si dans les exemples que l’on a vu la fonction F ne d´epend pas
de t.
• t0 est la borne inf´erieure de l’intervalle d’´etude.
• y0 est la valeur initiale, on a y(t0 ) = y0 .
• T est la longueur de notre intervalle d’´etude : [t0 , t0 + T ].
T
• N est le nombre de subdivisions de notre intervalle, ainsi le pas h vaut .
N
Remarques : i) Le fait de renvoyer X et Y nous permettra ensuite de faire rapidement un trac´e de la ligne polygonale
bris´ee obtenue.
ii) Il faudra d´efinir la fonction F qui correspond au probl`eme ´etudi´e au pr´ealable.
4
Chapitre 11
MPSI2
1.4
1.4.1
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Test du programme sur les exemples
Radioactivit´
e
On reprend l’´equation d´efinie dans la partie 1.1.1 avec λ = 0.5, c’est-`a-dire que l’on consid`ere :
0
N (t) = −0.5N (t)
N (0) = 1
De plus, on choisit T = 7, on travaillera donc sur l’intervalle [0, 7]. Voici les trac´es obtenus pour N = 5, N = 10, N = 20 et
N = 100. On y a fait ´egalement figurer la solution exacte : N (t) = e−0.5t que l’on connaˆıt dans ce cas-l`a.
Pour N = 100 les courbes semblent confondues. Les r´esultats obtenus sont instantan´es. Pour obtenir ces courbes, on a
d´efini la fonction F ainsi :
Puis, on fait appel `
a notre fonction Euler :
5
Chapitre 11
MPSI2
1.4.2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Croissance logistique
On reprend l’´equation logistique :
P (t) P 0 (t) = r 1 −
P (t)
K
avec r = 0.1 et K = 1000.
Cette fois-ci, on fait varier la condition initiale P (0) (qui correspond `a la population `a l’instant t = 0) de 200 `
a 2000. On
choisit l’intervalle [0, 12] et N = 1000. Voici le script utilis´e et les graphiques obtenus :
On constate que si la population initiale est au-dessus de la capacit´e d’accueil, la population d´ecroˆıt et semble tendre
vers K = 1000. Inversement, si la population initiale est en-dessous de la capacit´e d’accueil, la population semble tendre vers
1000 en croissant.
6
Chapitre 11
MPSI2
1.4.3
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Chute libre
On reprend l’´equation diff´erentielle :
(
v 0 (t)
=
v(t0 )
=
α
v(t)2 − g
m
v0
Choisissons α = 0.29, m = 100, g = 9.81, v0 = 0, T = 60 et N = 1000, voici la courbe obtenue :
r
gm
). Notre graphique, nous indique que cette vitesse
α
limite vaut environ 58m.s−1 , c’est-`
a-dire 216km.h−1 , cela correspond `a ce que l’on peut mesurer en pratique.
Comme attendu, il y a une vitesse limite (elle vaut th´eoriquement
1.5
Deux exemples probl´
ematiques
On consid`ere le probl`eme de Cauchy :
(
y(1) = 1
y 0 (t) = 3
y(t)
5
− 3
t
t
1
On peut v´erifier que l’unique solution est y : t 7→ 2 . Pourtant, on constate que la solution approch´ee obtenue `
a l’aide de
t
la m´ethode d’Euler s’´ecarte assez rapidement de la solution exacte.
7
Chapitre 11
MPSI2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
1
et la condition y(1) = 1 impose λ = 0. D`es que l’on s’´ecarte de la
t2
solution, on suit les tangentes des solutions qui comportent un terme en λt3 . Ce n’est pas la m´ethode d’Euler qui est en
cause, c’est juste que le probl`eme est mal pos´e.
La forme g´en´erale des solutions est t 7→ λt3 +
Exercice 2 : On consid`ere le probl`eme de Cauchy :
y(1) = 1
y 0 (t) = 100(sin(t) − y(t))
On applique le sch´ema num´erique d’Euler et l’on suppose que l’on a fait une erreur εn sur yn .
1. Quelle sera l’erreur sur yn+1 ?
2. Et sur les valeurs suivantes ?
3. Expliquer le graphique suivant qui est obtenu avec h = 0.02.
On dit que le probl`eme est mal conditionn´e.
8
Chapitre 11
MPSI2
1.6
1.6.1
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Qualit´
e d’un sch´
ema num´
erique
Consistance
L’erreur de consistance donne l’ordre de grandeur de l’erreur effectu´ee `a chaque pas. Dans le cas de la m´ethode
d’Euler, l’erreur de consistance mesure l’erreur qu’entraˆıne le fait d’approcher le nombre d´eriv´e par un taux d’accroissement.
La somme des erreurs de consistance `
a chaque pas donnera l’ordre de grandeur de l’erreur globale. On peut d´emontrer en
faisant certaines hypoth`eses sur la fonction F que si le pas h tend vers 0, l’erreur globale tend vers 0. D’un point de vue
th´eorique, cela signifie que plus le pas est petit, meilleure sera l’approximation.
Un sch´ema num´erique est dit consistant si la somme des erreurs de consistance tend vers 0 quand le pas tend vers 0.
C’est le cas de la m´ethode d’Euler.
1.6.2
Stabilit´
e
La stabilit´
e contrˆ
ole la diff´erence entre deux solutions approch´ees correspondant `a des conditions initiales voisines : un
sch´ema est dit stable si un petit ´ecart entre les conditions initiales et de petites erreurs d’arrondi dans le calcul des (yn )
provoquent une erreur finale control´ee lin´eairement par la perturbation initiale. La m´ethode d’Euler est une m´ethode stable.
1.6.3
Convergence
Un sch´ema num´erique est dit convergent lorsque l’erreur globale qui est le maximum des ´ecarts entre la solution exacte
et la solution approch´ee tend vers 0. Pour que le sch´ema soit convergent, il suffit que la m´ethode soit consistante et stable.
1.6.4
Erreurs d’arrondi
Pour la mise en application du sch´ema, il faut aussi prendre en compte les erreurs d’arrondi. Pour minimiser l’erreur
globale, on pourrait ˆetre tent´e d’appliquer la m´ethode d’Euler avec un pas tr`es petit, par exemple de l’ordre de 10−16 . Mais
en faisant ceci, le temps de calcul deviendrait irr´ealiste et les erreurs d’arrondi feraient diverger la solution approch´ee tr`es
rapidement puisque pour les flottants de l’ordre de 10−16 , les calculs ne sont plus du tout exacts.
En pratique, il faut prendre le pas h assez petit pour que la solution approch´ee soit proche de la solution exacte mais pas
trop petit pour ´eviter des erreurs d’arrondi et avoir un temps de calcul raisonnable.
9
Chapitre 11
MPSI2
2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Equations d’ordre 2 ou plus
2.1
Vectorisation du probl`
eme
L’´equation d’un pendule simple sans frottement est de la forme :
θ00 (t) + ω 2 sin(θ(t)) = 0
avec ω =
g
l
o`
u g est l’acc´el´eration de la pesanteur et l la longueur de la tige. Une telle ´equation ne peut pas s’´ecrire directement sous la
forme y 0 (t) = F (t, y(t)) puisqu’une d´eriv´ee seconde intervient dans l’expression.
L’id´ee est alors de consid´erer le vecteur Y = (θ, θ0 ), on a Y 0 = (θ0 , θ00 ). Ce qui donne la relation Y 0 = G(Y ) avec
G : (α, β) 7→ (β, −ω 2 sin(α))
Par ce proc´ed´e de vectorisation, on ram`ene les ´equations diff´erentielles d’ordre p `a des ´equations d’ordre 1 mais avec
des fonctions `
a valeurs dans Rp .
Exercice 3 : Lorsque l’´equation diff´erentielle est lin´eaire, il est ´egalement possible d’utiliser une ´ecriture matricielle. Par
exemple, on consid`ere l’´equation diff´erentielle :
y (3) + ty 00 (t) − t2 y 0 (t) + 2y(t) = cos(t)
´
l’´equation diff´erentielle sous la forme Y 0 = AY + B o`
u A est une matrice 3 × 3.
On pose le vecteur Y = (y, y 0 , y 00 ). Ecrire
2.2
Adaptation du programme
La m´ethode d’Euler s’adapte pour des fonctions a` valeurs vectorielles. Voici cette adaptation en dimension 2, c’est-`
a-dire
que la fonction F est a` valeurs dans R2 et les valeurs approch´ees de la solution yn poss`edent deux coordonn´ees donc vont
ˆetre rang´ees dans une matrice `
a deux lignes. On utilise les matrices du module numpy (type array).
Prenons l’exemple de l’´equation diff´erentielle :
 00
 y (t) + y(t) = 0
y(0) = 0
 0
y (0) = 1
Si l’on pose Y = (y, y 0 ), ce syst`eme peut se r´e´ecrire :
0
Y (t) = F (t, Y (t))
Y (0) = (0, 1)
o`
u F (t, (α, β)) = (β, −α)
On commence par cr´eer la fonction F :
Puis, on fait appel `
a la fonction Euler2d :
Ce qui nous int´eresse est la premi`ere ligne de la matrice Y ainsi renvoy´ee puisqu’elle contient des valeurs approch´ees de
la fonction y. On peut faire le trac´e :
10
Chapitre 11
MPSI2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Comme attendu, on trouve une approximation de la fonction sinus qui est bien l’unique solution du probl`eme de Cauchy
pr´ec´edent.
2.3
Exemple de l’´
equation du pendule simple
On reprend l’´equation du pendule simple sans frottement
θ00 (t) + ω 2 sin(θ(t)) = 0
avec ω =
g
l
On a vu que le syst`eme va s’´ecrire sous la forme Y 0 = G(Y ) avec
G : (α, β) 7→ (β, −ω 2 sin(α))
Prenons comme param`etres g = 9.81 m.s−2 , l = 1 m, on d´efinit la fonction G :
hπ i
On utilise ensuite notre fonction Euler2d avec comme param`etre initial y0 =
, 0 , ce qui correspond `a un angle initial
4
π
de et une vitesse initiale nulle. La fonction nous renvoie le vecteur X qui correspond `a la discr´etisation des temps ainsi que
4
le vecteur Y qui comporte une approximation de θ et θ0 . Il est int´eressant de tracer le graphique dans l’espace des phases,
c’est-`
a-dire θ en abscisses et θ0 en ordonn´ees.
11
Chapitre 11
MPSI2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Le mauvais raccord, qui ne devrait pas avoir lieu puisqu’il n’y a pas de frottement, est dˆ
u `a l’approximation dans la
m´ethode d’Euler. Si l’on augmente le nombre d’it´erations, c’est-`a-dire si l’on diminue le pas de la m´ethode, cette imperfection
disparaˆıt :
Le graphique obtenu est conforme `
a l’intuition, quand l’angle est maximal ou minimal, la vitesse s’annule.
12
Chapitre 11
MPSI2
2.4
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Syst`
eme proies pr´
edateurs
Nous venons de voir qu’il est possible de ramener une ´equation diff´erentielle d’ordre sup´erieur `a 1 `a une ´equation d’ordre
1 vectorielle et que l’adaptation de la m´ethode d’Euler permet d’obtenir ´egalement une solution approch´ee. Il est possible
de se servir de notre fonction Euler2d afin de r´esoudre un syst`eme d’´equations diff´erentielles, comme le montre l’exemple
suivant, tr`es classique.
Lotka et Volterra ont propos´e au d´ebut du 20i`eme si`ecle un mod`ele simple pour d´ecrire les relations entre proies et
pr´edateurs. Consid´erons deux populations animales de proies (effectif y1 (t)) et de pr´edateurs (effectif y2 (t)). On choisit la
mod´elisation suivante :
• en l’absence de pr´edateurs, l’´evolution des proies est r´egie par l’´equation
y10 (t) = ry1 (t)
o`
u r est le taux d’accroissement intris`eque des proies.
• en pr´esence de pr´edateurs, la population des proies subit un pr´el`evement proportionnel au produit des deux populations :
y10 (t) = ry1 (t) − py1 (t)y2 (t)
• en l’absence de proies, la population des pr´edateurs p´ericlite selon l’´equation :
y20 (t) = −dy2 (t)
• en pr´esence de proies, on ajoute un terme d’accroissement proportionnel au produit des deux populations :
y20 (t) = −dy2 (t) + qy1 (t)y2 (t)
Finalement, lorsque les deux populations sont pr´esentes sur un mˆeme territoire, leur ´evolution conjointe est donc r´egie
par le syst`eme :
 0
y (t) = ry1 (t) − py1 (t)y2 (t)


 10
y2 (t) = −dy2 (t) + qy1 (t)y2 (t)
y1 (t0 ) = n1



y2 (t0 ) = n2
Comme dans les exemples pr´ec´edents, on vectorise la situation en posant Y (t) = (y1 (t), y2 (t)) et Y 0 (t) = (y10 (t), y20 (t)).
Le syst`eme pr´ec´edent se r´esume sous la forme suivante :
0
Y (t) = F (t, Y (t))
Y (t0 ) = (n1 , n2 )
avec F (t, y1 , y2 ) = (ry1 − py1 y2 , −dy2 + qy1 y2 )
Afin d’illustrer la situation, on choisit les param`etres r = 2.5, p = 0.05, d = 0.8 et q = 0.05. On peut d´efinir la fonction
F :
Voici les r´esultats obtenus avec n1 = 100, n2 = 100 et T = 30 :
13
Chapitre 11
MPSI2
R´
esolution num´
erique d’´
equations diff´
erentielles
Evolution des proies en fonction du temps :
Informatique
Evolution des pr´edateurs en fonction du temps :
Si on les repr´esente sur un mˆeme graphique, on obtient :
Ce mod`ele peut ˆetre v´erifi´e par des observations, comme le montrent ces mesures effectu´ees au Canada entre 1845 et
1955. On constate que les courbes sont globalement p´eriodiques, l’une ´etant retard´ee par rapport `a l’autre.
14
Chapitre 11
MPSI2
3
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Autres sch´
emas num´
eriques
On reprend la forme g´en´erale de notre probl`eme de Cauchy :
0
y (t) = F (t, y(t))
(C) :
y(t0 ) = y0
On rappelle que la m´ethode d’Euler consiste `
a choisir un pas h et `a construire deux suites (tn )0≤n≤N et (yn )0≤n≤N de la
fa¸con suivante :

= y(t0 )
 y0
tn+1 = tn + h

yn+1 = yn + hF (tn , yn )
3.1
M´
ethode d’Euler implicite
On garde les mˆemes notations mais on modifie l´eg`erement le sch´ema num´erique :

= y(t0 )
 y0
tn+1 = tn + h

yn+1 = yn + hF (tn+1 , yn+1 )
La derni`ere ligne est une ´equation en yn+1 , d’o`
u le nom de sch´
ema d’Euler implicite. On peut imaginer la r´esoudre
explicitement ou appliquer par exemple la m´ethode de Newton pour trouver yn+1 .
Sous certaines hypoth`eses le sch´ema d’Euler implicite est plus stable, c’est-`a-dire que les erreurs num´eriques se propagent
moins vite que le sch´ema d’Euler classique (dit sch´ema d’Euler explicite).
Cette m´ethode se g´en´eralise ´egalement aux dimensions sup´erieures mais elle peut n´ecessiter par exemple d’inverser une
matrice si l’´equation est lin´eaire.
15
Chapitre 11
MPSI2
3.2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Sch´
ema de Runge-Kutta d’ordre 4 : RK4
Karl Runge
Martin Wilhem Kutta
La m´ethode de Runge-Kutta d’ordre 4 est celle qui pr´esente le meilleur rapport pr´ecision/complexit´e, elle est tr`es utilis´ee
en pratique. Avec les notations usuelles, on a :

k1 = hF (tn , yn )









k1 h


y
=
y(t
)
0


 k2 = hF tn + 2 , yn + 2
 0
tn+1 = tn + h
u
o`



 yn+1 = yn + 1 k1 + 2k2 + 2k3 + k4
k2 h


,
y
+
k
=
hF
t
+

n
3
n
6

2
2





 k = hF (t + h, y + k )
4
n
n
3
Les coefficients qui apparaissent dans ces formules sont judicieusement ajust´es pour obtenir une m´ethode d’ordre 4, c’esta`-dire que l’erreur de consistance va ˆetre major´ee par une quantit´e de l’ordre de grandeur de h4 . En bref, quand le pas tend
vers 0, la solution approch´ee converge rapidement vers la solution exacte. C’est `a comparer avec les m´ethodes d’Euler qui
sont d’ordre 1.
On reprend le probl`eme de Cauchy dont on sait que la fonction exponentielle est solution :
0
y (t) = y(t)
y(0) = 1
Voici les solutions approch´ees obtenues par la m´ethode d’Euler et la m´ethode RK4 pour diff´erentes valeurs de N .
16
Chapitre 11
MPSI2
R´
esolution num´
erique d’´
equations diff´
erentielles
Informatique
Cette m´ethode se g´en´eralise ´egalement pour des ´equations diff´erentielles d’ordre sup´erieur et pour des syst`emes d’´equations
diff´erentielles.
17
Chapitre 11