Université Paris-Saclay - Année 2022-2023

Cours accéléré d'analyse numérique - M2 AMS

I - Schémas de diérences nies pour les équations diérentielles.

1 L'approximation numérique des équations diérentielles.

On s'intéresse à l'approximation numérique des équations diérentielles, ordinaires et aux dérivées
partielles. C'est surtout les deuxièmes qui nous intéresseront mais pour introduire la méthode des
diérences nies, on commencera par les premières.
Les équations diérentielles correspondent souvent à des modèles mathématiques décrivant des phéno-
mènes physiques. Dans beaucoup de cas, on ne sait pas calculer une solution explicite d'une équation
donnée, et on doit utiliser des techniques de résolution approchée. On s'intéresse alors à la discréti-
sation et à la résolution numérique de l'équation, autrement dit on remplace l'équation, qui est posée
sur un domaine continu, par un problème discret, en discrétisant le domaine. Si l'équation originale est
linéaire on obtient, en discrétisant, un système linéaire, si l'équation originale est non linéaire, on peut
obtenir une équation non linéaire à résoudre, par exemple par une méthode de résolution approchée
d'équations non linéaires comme la méthode de Newton. Les principales méthodes de discrétisation
pour une équation diérentielle sont les méthodes des diérences nies, la méthode des éléments nis et
la méthode des volumes nis. Il y a aussi d'autres classes de méthodes, comme les méthodes spectrales.
Au début de ce cours on s'intéresse à la méthode des diérences nies.

1.1 La méthode des diérences nies - approximation des dérivées.

On considère une équation diérentielle, dont l'inconnue est une fonction u, posée sur un domaine
physique Ω ⊆ Rd , d ≥ 1.
Le principe de la méthode des diérences nies est d'approcher la solution u de cette équation en un
ensemble discret (mais grand) {x1 , . . . , xN } de points du domaine Ω, en remplaçant les dérivés de u
aux points xi par des quotients de diérences faisant intervenir les points voisins de xi . Cette approche
sera valable si la solution du problème est régulière.
Par exemple, on a pour une fonction u d'une seule variable que sa dérivé en un point x de son domaine
vérie
u(x + h) − u(x) u(x) − u(x − h) u(x + h) − u(x − h)
u0 (x) = lim = lim = lim = ...,
h→0 h h→0 h h→0 2h
et on peut donc approcher u0 (x) par les quotients
u(x + h) − u(x) u(x) − u(x − h) u(x + h) − u(x − h)
, , ,...,
h h 2h
avec h petit.
Ces quocients de diérences nies font intervenir des valeurs de u en un nombre ni (dans chacun de
ces exemples, deux) de points du domaine proche de x. On va vérier un peu plus en bas que l'on peut
aussi approcher u0 (x) par exemple par la quantité
2u(x + h) + 3u(x) − 6u(x − h) + u(x − 2h)
,
6h
pour h petit.
On a aussi que si u est deux fois dérivable au point x, alors
u(x + h) − 2u(x) + u(x − h)
u00 (x) = lim
h→0 h2

1
et on peut donc aussi approcher u00 (x) par le quotient de diérences nies
u(x + h) − 2u(x) + u(x − h)
,
h2
avec h petit.
Une question que l'on se posera souvent dans le cadre de la méthode des diérences nies est celle de
l'ordre de l'approximation. On va expliciter cette notion plus loin, mais essayons de la comprendre dès
maintenant. Elle mesure l'erreur commise dans ces approximations.

Supposons a, b ∈ R et u : [a, b] −→ R une fonction de classe C 4 . Soit x ∈]a, b[.

Exercice 1.
1. Montrer, en utilisant des développements de Taylor de u autour de x, que, lorsque h → 0,
u(x + h) − u(x) u(x) − u(x − h)
u0 (x) − = O(h), u0 (x) − = O(h),
h h

u(x + h) − u(x − h)
u0 (x) − = O(h2 ).
2h
et que
2u(x + h) + 3u(x) − 6u(x − h) + u(x − 2h)
u0 (x) − = O(h3 ).
6h
2. Construction d'une approximation de diérences nies à un ordre précis.
Chercher une approximation par diérences nies de u0 (x) basée sur les points x, x − h et x − 2h
de la forme au(x) + bu(x − h) + cu(x − 2h) tel que

u0 (x) − au(x) + bu(x − h) + cu(x − 2h) = O(h2 ).

Les diérences ci-dessus s'appellent erreurs de troncature ou de consistance au point x. On dira

que cette erreur est d'ordre 1 si elle se comporte en O(h), d'ordre p > 0 si elle se comporte en O(hp ).
Plus p est grand, plus l'approximation choisie de la dérivée de u est précise.

L'exercice 1.2 nous fait penser que l'on peut approcher la dérivée d'une fonction u en un point x par
une formule de diérences nies à un ordre souhaité, en utilisant des valeurs de u en points autour de
x, en cherchant les bons poids pour chaque valeur. Mais on voit facilement que le plus élévé est l'ordre
souhaité, le plus de points il est nécéssaire d'utiliser.

Exercice 2. [Interpolation polynomial]

On peut approcher la dérivée d'une fonction u en un point x en approchant la fonction u par un
polynôme X 7→ P (X) et en calculant P 0 (x). Vérier que si l'on approche u par ses polynômes d'inter-
polation de dégré 1 aux points x et x − h, ou x et x + h, ou x − h et x + h, et que l'on calcule les dérivées
de ces pôlynômes au point x, on obtient les formules de diérences nies correspondantes que l'on a
vu précédamment. On peut faire le même exercice en approchant u par le polynôme d'interpolation
(de degré 2) aux points x, x − h et x − 2h et par le polynôme (de degré 3) qui interpole u aux points
x, x + h, x − h et x − 2h.

1.1.1 Approximation des dérivées dans un ensemble de points d'un intervale.

N = h. Considérons alors les N + 1 points xn =
Soit [a, b] ⊆ R et h > 0 tel qu'il existe N ∈ N tel que b−a
a + nh, n = 0, . . . , N . L'ensemble discret formé par les points x0 , . . . , xN s'appelle une discrétisation
régulière ou uniforme de [a, b] de pas h.

2
Soit u une fonction régulière dénie sur [a, b] et soit U ∈ RN +1 le vecteur (U0 , . . . , UN ) avec Un =
u(xn ), n = 0, . . . , N .
On souhaite ici approcher les dérivées de u en chaque point xn de la discrétisation de l'intervale [a, b]
au même ordre, en utilisant le même type de diérences nies en chaque point. La diculté en faire
ceci consiste dans l'approximation des dérivées aux points {a, b} du bord de l'intervale et en général
dans ces points on doit utiliser une approximation diérente de celle utilisée aux points de l'interieur
du domaine. Ce problème se pose lors de la discrétisation des conditions aux limites pour un problème
aux limites pour une équation aux dérivées partielles.

Exercice 3.
1. Déterminer en fonction de h et de U un vecteur V ∈ RN +1 vériant

max |Vn − u0 (xn )| = O(h).

n=0,...,N

2. Déterminer en fonction de h et de U un vecteur W ∈ RN +1 vériant

max |Wn − u0 (xn )| = O(h2 ).

n=0,...,N

3. A faire sur python.

Illustrer numériquement les deux approximations de la dérivée de u pour la fonction u(x) =
sin(x) dans l'intervalle [a, b] = [0, 2π]. Pour ce faire représenter dans la même gure la fonction
10 et puis avec h = 100 . Représenter
u0 , le vecteur V et le vecteur W , en fonction de x, avec h = 2π 2π

dans une autre gure les diérences |u − V | et |u − W |.

0 0

On peut construire une discrétisation non uniforme de l'intervale [a, b], c'est-à-dire telle que la diérence
entre chaque deux points consécutifs de la discrétisation ne soit pas toujours égale. Dans ce cas on
appelera pas de la discrétisation la plus grande distance entre chaque deux points consécutifs de la
discrétisation.

1.1.2 Analyse de l'ordre de précision.

La méthode des diérences nies pour approcher la solution u d'une équation diérentielle posée dans
un domaine Ω consiste à remplacer l'équation diérentielle par une équation aux diérences nies, en
remplaçant les dérivées de u par des quocients de diérences nies en des points obtenues en discrétisant
le domaine Ω avec un pas h > 0.
Nous avons vu qu'en utilisant des développements de Taylor relativement simples, on peut obtenir
l'ordre de l'erreur commmise entre la dérivée et son approximation par diérences nies. Dans les
exemples que l'on a vu, cette erreur se comporte en O(hp ), avec p ∈ N. En géneral on ne sait pas
exactement comment se comportent ses erreurs. Une manière d'avoir une idée de ce comportement
est de représenter graphiquement la diérence entre la fonction et ses dérivées, pour des valeurs de h
diérentes, et de chercher à comprendre comment se comporte cette erreur en fonction de h.
Plus précisément, considérons toujours le cas d'une fonction u dénie sur un intervalle [a, b] et, pour
h > 0, une discrétisation de l'intervalle [a, b] de pas h, donnée par N + 1 points x0 , . . . , xN (N dépend
de h).
Supposons de l'on approche la dérivée u0 (xn ) en chaque point xn de la discrétisation par un quocient de
diérences nies donnée par une formule Du(xn ). On note E(h) l'erreur globale commise (par exemple
en norme innie) avec le pas h : E(h) = maxn=0,...,N |u0 (xn ) − Du(xn )|. Supposons que l'on s'attend à
ce que E(h) se comporte comme hp pour un certain p ∈ N, c'est-à-dire qu'il existe C > 0 indépendant
de h tel que E(h) ∼ Chp . Alors on s'attend à ce que log E(h) se comporte comme

log E(h) ∼ log(C) + p log(h).

3
On peut alors représenter dans un graphique log(E(h) en fonction de log(h), pour des valeurs de h de
plus en plus petits. Si le comportement est l'attendu, on obtiendra que, dans cette échelle log-log, le
graphique de E(h) en fonction de h est une droite de pente p, et cette pente de la droite nous donne
l'ordre de la méthode.
Exercice 4. À faire en python.
Soit u(x) = sin(x).
1. Choisir un point quelconque x et représenter, pour plusieurs valeurs de h de plus en plus pe-
tites, la diérence en valeur absolue entre u0 (x) et, respectivement u(x+h)−u(x)
h , u(x+h)−u(x−h)
2h
et 2u(x+h)+3u(x)−6u(x−h)+u(x−2h)
6h , en fonction de h, en échèlle logarithmique. Représenter des
droites de pente 1, 2 et 3 sur le même graphique pour visualiser comment se comporte cette
diérence en fonction de h.
2. Faire le même exercice mais en considérant maintenant des discrétisations uniformes de pas h
de par exemple l'intervale [0, π], et calculer l'erreur globale des mêmes approximations sur cet
intervale. Pour approcher les points de bord dans la dernière formule on peut utiliser les valeurs
exactes ici.

2 Approximation des EDO.

On s'intéresse à des méthodes numériques pour approcher les équations diérentielles ordinaires.

On considère le problème de Cauchy pour une équation diérentielle ordinaire, ou pour un système
d'équations diérentielles ordinaires, de la forme
(
y 0 (t) = f (t, y(t)),
(1)
y(t0 ) = y0 ,

où f : I × Rn −→ Rn est une fonction continue et localement Lipschitzienne par rapport à sa deuxième

variable, avec I ⊆ R un intervalle ouvert de R et où (t0 , y0 ) ∈ I × Rn est donné. Le problème de Cauchy
(1) admet alors une unique solution maximale dénie dans un intervalle ouvert J ⊆ I : ce résultat est
une conséquence du théorème de Cauchy-Lipschitz.
On souhaite calculer une solution approchée (discrète) de ce problème dans un intervalle de la forme
[t0 , tf ] = [t0 , t0 + T ] ⊆ J , avec T > 0.
Pour cela on procède comme suit.
i) On discrétise (c'est-à-dire on remplace du continu par du discret) l'intervalle [t0 , t0 + T ]. Pour
ce faire on se donne N ∈ N et on dénit une sub-division de [t0 , t0 + T ] en N sous-intervalles
dénis par les N + 1 points, appelées points de la discrétisation,
t0 < t1 < · · · < tN = t0 + T.

On appelle aussi les points tn les instants temporels, car pour de nombreuses EDO modélisant
des phénomènes physiques, la variable t représente le temps.
Le point t0 est appelé l'instant initial et le point tN l'instant nal. On appelle aussi :
hn := tn+1 − tn le pas de temps ou le pas de discrétisation entre tn et tn+1 , pour
n = 0, . . . , N − 1 ;
h := max
n=0,...,N −1
hn le pas de temps ou pas de la discrétisation.
La discrétisation est dite uniforme si hn = h, pour tout n = 0, . . . , N − 1. Dans ce cas la pas
h de la discrétisation est déni par
T
h=
N

4
 et les N + 1 points de discrétisation sont dénis par
tn = t0 + nh, n = 0, . . . , N.
ii) On construit N + 1 valeurs y 0 , . . . , y N qui approchent la solution exacte y de (1), aux points
t0 , . . . , tN , autrement dit on construit des valeurs y 0 , . . . , y N tels que
y n ≈ y(tn ), pour n = 0, . . . , N.
Une méthode numérique ou un schéma numérique est donc la donnée d'une telle construction.
Plus précisément, c'est la donnée d'une suite (y 0 , . . . , y N )N ∈N , indexée par N ∈ N, de valeurs ap-
prochant les valeurs exactes (y(t0 ), . . . , y(tN )) de la solution y de (1) aux points de la discrétisation
(t0 , . . . , tN ). On appelle souvent la suite (y 0 , . . . , y N )N ∈N de solution approchée ou solution nu-
mérique et, pour N ∈ N donnée, le terme d'ordre N de cette suite solution approchée (ou solution
numérique) associée à une discrétisation de l'intervalle [t0 , t0 +T ] de N +1 points ou de pas h = max hn .
0,...,N
Comme la valeur de la solution exacte à l'instant initial y(t0 ) = y0 est connue, on prendra = y0 . y0
La dépendance de la suite (y 0 , . . . , y N ) par rapport à N n'est pas toujours spéciée mais il faut la gar-
der à l'esprit, ainsi que la dépendance de la suite (t0 , . . . , tN ) dénissant les points de la discrétisation
de l'intervalle [t0 , t0 + T ].

Lorsque N augment, le nombre de points de la discrétisation augmente et on s'attend à ce que les valeurs
approchées y 0 , . . . , y N approchent de manière de plus en plus précise les valeurs exactes y(t0 ), . . . , y(tN ).

Nous allons nous restreindre à des discrétisations uniformes de l'intervalle [t0 , t0 + T ] (le cas d'une
discrétisation à pas hn variable est souvent utilisé lorsque l'on veut adapter le pas de la discrétisation
pour améliorer la précision du schéma). Dans ce cas, comme h = NT , on peut de manière équivalente
voir la suite approchée (y 0 , . . . , y N ) comme une suite indexée par N ∈ N ou comme une suite indexée
par h ∈]0, 1] ; on a choisi l'intervalle ]0, 1] pour h mais on aurait pu choisir n'importe quel autre in-
tervalle de la forme ]0, hmax ]. Ce qui nous intéressera sera le comportement de la solution approchée
(y 0 , . . . , y N ) lorsque N tend vers l'inni (i.e., lorsque le nombre de points de la discrétisation tend vers
l'inni, ce qui correspond à ce que l'ensemble discret qui représente l'intervalle [t0 , t0 + T ] tend vers cet
intervalle), ou de manière équivalente, lorsque h tend vers 0 (i.e. lorsque la distance entre deux points
consécutifs de la discrétisation tend vers 0).

Dans cette feuille on va considérer les schémas numériques suivants an de calculer une approximation
y n de y(tn ) pour n = 0, . . . , N .
On pose y 0 = y0 , et, pour n = 0, . . . , N − 1,

Euler explicite : y n+1 = y n + hf (tn , y n )

Heun : y n+1 = y n + h n n + hf (tn , y n )



2 f (tn , y ) + f tn + h, y
Runge Kutta 4 : y n+1 = y n + h
6 (p1 + 2p2 + 2p3 + p4 ) où p1 = f (tn , y n ),
p2 = f (tn + h2 , y n + h2 p1 ),
p3 = f (tn + h2 , y n + h2 p2 ),
p4 = f (tn + h, y n + hp3 )

Ces trois méthodes sont des méthodes explicites (la valeur y n+1 se calcule explicitement à partir de
valeurs antérieures), à un pas (la valeur y n+1 se calcule uniquement à partir de la valeur précédente
y n ). Ces méthodes peuvent s'écrire sous une forme générique
(S) y n+1 = y n + hΦ(tn , y n , h),

5
pour n = 0, . . . , N − 1, avec y 0 donné (nous considérerons toujours y 0 = y0 ), où Φ : [t0 , t0 + T ] ×
Rn × [0, 1] −→ Rn est une fonction continue. Par exemple, dans le cas de la méthode d'Euler,
Φ(t, y, h) = f (t, y).

2.1 Questions théoriques - convergence d'une méthode numérique.

La convergence d'une méthode numérique pour une EDO est liée à deux notions : la consistance et
la stabilité. Nous allons retrouver ces notions dans le cadre des méthodes de diérences nies pour
les EDP.
Considérons une méthode numérique dénie par (S). Soit y la solution exacte du problème de Cauchy
(1). On appelle erreur locale à l'instant tn de la méthode, la quantité
en = y(tn ) − y n .

La méthode dénie par (S) est dite convergente sur l'intervalle [t0 , t0 +T ] si pour toute donnée initiale
y0 ,  
lim max ken k = 0,
h→0 n=0,...,N

autrement dit si pour toute donnée initiale y0 ,

 
lim max ky(tn ) − y n k = 0.
h→0 n=0,...,N

On remarque que la solution exacte y de (1) vérie

y(tn+1 ) = y(tn ) + hΦ(tn , y(tn ), h) + εn , n = 0, . . . , N − 1,

où la quantité
εn = y(tn+1 ) − y(tn ) − hΦ(tn , y(tn ), h)
s'appelle erreur de consistance de la méthode (S) à l'étape n, relative à la solution y . C'est l'erreur
commise par une itération du schéma sur la solution exacte. La méthode est dite consistante si la
somme des erreurs de consistance est petite :
N
X −1
lim kεn k = 0;
h→0
n=0

la méthode est dite consistante d'ordre au moins p s'il existe une constante C > 0, indépendante
de h, tel que
N
X −1
∀ N ∈ N, kεn k ≤ Chp .
n=0

La méthode est consistante si l'erreur commise par une étape du schéma est petite car la discrétisation
est cohérente avec l'EDO.

La méthode est dite stable s'il existe une constante M > 0, indépendante de h, tel que pour toutes
suites y n et yen dénies par
y n+1 = y n + hΦ(tn , y n , h), n = 0, . . . , N − 1, y 0 donné,
yen+1 = yen + hΦ(tn , yen , h) + ρn , n = 0, . . . , N − 1, ye0 donné,

avec ρn ∈ R, n = 0, . . . , N − 1, on ait
 N
X −1 
n n 0 0
max ky − ye k ≤ M ky − ye k + kρn k .
n=0,...,N
n=0

6
Une méthode stable est alors une méthode telle que, lorsque l'on introduit une perturbation ρn à chaque
étape du schéma, l'erreur totale de la méthode est contrôlée par le cumul des perturbations.

Exercice 5. Montrer que si la méthode est stable et consistante, alors elle est convergente. Montrer que
si la méthode est stable et consistante d'ordre au moins p alors il existe une constante C̃ indépendante
de h tel que
max ky(tn ) − y n k ≤ C̃hp .
n=0,...,N

On utilise souvent la formulation  méthode d'ordre p  au lieu de  méthode consistante d'ordre p pour
une méthode dont l'erreur de consistance est d'ordre p. Ce dernier exercice montre que si la méthode
est consistante d'ordre p, l'erreur globale du schéma se comporte aussi comme O(hp ).

Exercice 6. Supposons la solution y du problème de Cauchy (1) régulière. On peut montrer, en

utilisant des développements de Taylor appropriés, que la méthode d'Euler explicite est consistante
d'ordre 1 et que la méthode de Heun est consistante d'ordre 2 (celle de RK4 est consistante d'ordre 4
mais les calculs sont plus compliqués).

2.2 Implémentation d'une méthode numérique sous python.

Quelques conseils pour l'organisation de votre programme : vous allez programmer plusieurs
schémas numériques et éventuellement les tester sur plusieurs exemples d'EDO. Il est souhaitable de
bien organiser votre script pour ne pas se perdre. Par exemple vous pourrez dénir au tout début les
fonctions second membre des EDO que vous allez tester, en les appelant par exemple f 1, f 2, . . . . Vos
fonctions doivent avoir toujours deux variables, t et y , même si elles ne dépendent pas de t (comme
ça les fonctions dénissant les schémas seront adaptées à tous les cas).Vous pourrez ensuite dénir les
solutions exactes de ces équations, si elles sont connues, en les appelant par exemple yex1, yex2, . . . .
Également, vous pouvez dénir ensuite les fonctions que vous allez construire pour chacun des schémas.
Lorsque vous faites appelle à une fonction pour tester un certain schéma dans une situation concrète,
il faut dénir avant les paramètres du cas que vous allez tester. Par exemple, faire
t0=0
t f =1
y0=2
h=0.5
euler_exp ( f1 , t0 , t f , y0 , h )
au lieu de
euler_exp ( f1 , 0 , 1 , 2 , 0 . 5 )
Cela vous permet de facilement changer les paramètres s'il le faut, et la visibilité de votre programme
sera aussi facilitée.

Exercice 7.

1. Pour chacune des trois méthodes numériques données, écrire une fonction python de la forme
meth ( f c t , t0 , t f , y0 , h )
avec meth = euler_exp, Heun ou RK4, prenant en argument fct la fonction f (t, y) de (1), les
extrémités t0 et tf de l'intervalle de temps [t0 , tf ] = [t0 , t0 + T ], la donnée initiale y0 et le pas
de temps h. Cette fonction devra retourner deux tableaux :

7
  [t0 , t1 , ..., tN ], tableau numpy unidimensionnel de taille (N + 1) × 1 représentant la subdi-
vision de l'intervalle [t0 , tf ] de pas h considérée,
 [y 0 , y 1 , ..., y N ], tableau numpy de taille (N + 1) × n représentant la solution approchée aux
instants tn , n = 0, . . . , N , donnée par la méthode choisie.
Remarque 1 : votre fonction meth pourra plutôt dépendre de N, le nombre de pas de la subdivision
considérée, au lieu de dépendre du pas h. Dans les deux cas, faire très attention à ce que le
nombre de points de la discrétisation soit cohérent avec le pas de la discrétisation.
Remarque 2 : la structure de votre fonction pourra aussi être autre que celle proposée.
Remarque 3 : vous pouvez dans un premier temps traiter le cas n = 1, d'une EDO scalaire, puis
essayez d'étendre votre fonction au cas n > 1.
2. Testez vos trois fonctions sur le modèle logistique
(
y 0 (t) = cy(1 − yb ),
(P1 )
y(0) = a,

dont la solution exacte est

b
y(t) = b−a −ct
,
1+ a e
avec les données c = 1, b = 2, a = 0.1, dans l'intervalle [0, 15], avec un pas h = 0.2. Tracer sur
la même fenêtre la solution exacte et les solutions approchées, obtenue avec le pas h et dans
une autre fenêtre avec un pas égal à h2 .
3. Testez ensuite vos fonctions dans le cas vectoriel n > 1 (n = 2) sur le problème
(
y 00 (t) = −y(t) + cos(t)
(P2 )
y(0) = 5, y 0 (0) = 1,

dont la solution exacte est

1
y(t) = sin(t)t + 5 cos(t) + sin(t),
2
dans l'intervalle [0, 15], avec un pas h = 0.2.
Pour ce faire, il faut écrire l'équation d'ordre 2 de (P2 ) sous la forme d'un système de deux
équations d'ordre 1 dans les nouvelles variables u(t) = y(t) et v(t) = y 0 (t). On se ramènera alors
à la résolution d'une équation de la forme
X 0 = F (t, X), avec X = (u, v) = (y, y 0 )T .
Représenter à nouveau la solution exacte et les solutions approchées dans une même fenêtre
graphique.
La solution y de (P2 ) correspond à la première composante du vecteur X ci-dessus. Votre fonction
euler_exp retournera dans ce cas, si vous avez respecté la structure conseillée, un tableau de
taille (N + 1) × 2, N étant le nombre de points de la discrétisation. Ce tableau donne les valeurs
approchées de X au points de la discrétisation considérée, la première colonne correspondant à
la première composante de X , la seconde à la seconde composante de X . La solution approchée
de (P2 ) que l'on cherche correspond alors à la première colonne de ce tableau.

Exercice 8. [Étude de l'ordre de la méthode]

On se donne un problème de Cauchy de la forme (1) et une méthode numérique pour approcher la
solution de (1) dans un intervalle de la forme [t0 , t0 + T ]. Pour un certain pas de temps h xé (ou pour
un certain nombre de points de la discrétisation N xé), l'erreur globale entre la solution approchée
associée à une discrétisation de pas h de l'intervalle [t0 , t0 + T ] et la solution exacte est donnée par :
Eh = max (|y(tk ) − y k |).
k=0,··· ,N

8
Ci dessus, y(tk ) est la solution exacte à l'instant tk et y k la valeur approchée de y(tk ), donnée par le
schéma numérique.

Remarque : l'erreur globale E dépend du pas h, ou, de manière équivalente, du nombre de points N
de la discrétisation. Parfois dans la littérature on ne spécie pas la dépendance de E par rapport à h,
mais on doit toujours garder à l'esprit cette dépendance et que l'erreur est donc une fonction de h.

Considérons le problème (
y 0 (t) = cos(t)−y(t)
1+t ,
(P3 ) 1
y(0) = − 4 ,
dont la solution exacte est
sin(t) − 1/4
y(t) = .
1+t
1. Calculez les solutions approchées de (P3 ) obtenues avec le schéma d'Euler explicite, avec h =
1/2s pour s = 1, 2, ..., 8 ; Calculer, pour chaque valeur de h = 1/2s , s = 1, 2, ..., 8, l'erreur glo-
bale Eh correspondante. Représentez ensuite, en échelle logarithmique, l'erreur en fonction du
pas de temps h, autrement dit, représentez log(Eh ) en fonction de log(h). Vous devez obtenir
des points qui sont à peu près alignés sur une droite de pente 1. Vériez graphiquement que
c'est le cas, en estimant la pente de la droite passant au plus prêt des points (ou en représentant
une droite de pente 1 qui passe par un des points et en vériant que tous les points sont à peu
près sur cette droite).
Ceci signie que log(Eh ) ∼ C + log(h) et donc que Eh ∼ Ch e , pour certaines constantes C et C e.
On dit que la méthode d'Euler explicite est d'ordre 1 : c'est l'ordre de la puissance de h dans
cette relation. On a donc que l'erreur globale Eh tend vers 0 comme h. L'ordre d'une méthode
donne une indication sur sa vitesse de convergence. Une méthode d'ordre p est une méthode dont
l'erreur globale tend vers 0 comme hp . Donc plus l'ordre est élevé, plus la méthode converge plus
vite.
Remarque : pour étudier numériquement l'ordre d'une méthode, on utilise souvent l'échelle lo-
garithmique pour tracer l'erreur en fonction du pas de discrétisation h. La pente de la droite
obtenue donne l'ordre p de la méthode : si Eh ∼ Chp alors log(Eh ) ∼ log(C) + p log(h).
2. Vous pouvez refaire l'exercice pour la méthode de Heun et de RK4. Vous allez conclure que la
méthode de Heun est d'ordre 2 et la méthode de RK4 d'ordre 4.

2.2.1 Un exemple de méthode implicite : simulation du système de Van der Pol par la
méthode d'Euler implicite en utilisant la méthode de Newton.
On considère le problème de Cauchy pour le système de Van der Pol :

x(t)3
(  
x0 (t) = ε x(t) − 3 + y(t),
y 0 (t) = −x(t),
de donnée initiale
x(0) = x0 , y(0) = y0 .
Nous allons calculer une solution approchée de ce problème par la méthode d'Euler implicite, dont
l'itération est donnée par
Y n+1 = Y n + hf (tn+1 , Y n+1 ).
Pour ce faire, nous allons utiliser la méthode de Newton pour résoudre l'équation non linéaire dénissant
l'itération de la méthode d'Euler implicite.

9
Le système de Van der Pol est un exemple d'un problème dit  raide . Il s'agit de problèmes dont la
solution présente d'importantes variations dans un intervalle de temps très court. Nous verrons qu'au
contraire des méthodes explicites, les méthodes implicites se comportent en général bien vis à vis de
problèmes dits  raides .

Dénition des fonctions du système.

Exercice 9.
Écrire le système de Van der Pol sous la forme

Y 0 = F (Y ),

avec Y = (x, y)t et F : R2 −→ R2 . Dénir la fonction F , ainsi que sa matrice Jacobienne DF en

complétant les fonctions python suivantes (respecter la structure donnée : F retourne un array de taille
2 de composantes F1 et F2 et DF un tableau de taille 2 × 2 avec les dérivées partielles de F ) .
# ep= A c o m p l e t e r avec l a v a l e u r du parametre e p s i l o n

def F(Y) :
x=Y [ 0 ]
y=Y [ 1 ]
#F1 = . . . # A c o m p l e t e r : premiere composante de F
#F2 = . . . # A c o m p l e t e r : deuxieme composante de F
return np . a r r a y ( [ F1 , F2 ] )
def DF(Y) :
x=Y [ 0 ]
y=Y [ 1 ]
# DFL1=np . array ( [ . . . , . . . ] ) # A c o m p l e t e r : premiere l i g n e de l a ma tr ic e
# j a c o b i e n n e de F en Y=(x , y )
# DFL2 = . . . # A c o m p l e t e r : deuxieme l i g n e de l a m at r ic e j a c o b i e n n e de F
return np . a r r a y ( [ DFL1 , DFL2 ] )

La méthode de Newton.

Soit G : Rn −→ Rn . La méthode de Newton pour résoudre l'équation G(Y ) = 0 est la méthode itérative
dénie par la suite
Yn+1 = Yn − (DG(Yn ))−1 G(Yn ),
avec Y0 donné. Sous certaines conditions, si Y0 est proche d'un zéro de G, la suite {Yn }n converge vers
ce zéro.
On va dénir une fonction python qui calcule une solution approchée de l'équation G(Y ) = 0 par la
méthode de Newton. Une telle fonction doit dépendre de :
 une valeur eps dénissant la tolérance avec laquelle on va calculer une solution de G(Y ) = 0
(ici on va calculer une valeur approchée Y N de Y tel que G(Y ) = 0 vériant |G(Y N )| ≤ eps ) ;
 la fonction G et sa matrice jacobienne DG ;
 l'initialisation Y0 ;
 un nombre maximal nmax d'itérations à ne pas dépasser.

Exercice 10.
Compléter la fonction python suivante permettant de calculer une valeur approchée d'un zéro d'une
fonction G par la méthode de Newton :

10
def Newton ( eps ,G,DG, Y0 , nmax ) :
t o l=np . l i n a l g . norm (G(Y0 ) )
n=1 # n = nombre d ' i t e r a t i o n s q u i s e r o n t e f f e c t u e e s
Y=Y0
while ( t o l >eps and n<nmax ) :
# A completer . . . .
return n ,Y

Pour implémenter la méthode d'Euler implicite pour le système de Van der Pol, on remarque que
l'itération de la méthode s'écrit sous la forme

Y n+1 = Y n + hF (Y n+1 ).

On peut donc la re-écrire sous la forme

Gn (Y n+1 ) = 0,

où Gn (Y ) = Y − Y n − hF (Y ). Pour implémenter le méthode d'Euler implicite, on doit donc à chaque

itération trouver la solution Y de Gn (Y ) = 0. La fonction Gn dépend de Y n et elle change alors à
chaque itération. On peut la dénir à chaque itération à partir de la fonction F en utilisant la fonction
lambda de python.

Exercice 11. Construire une fonction donnant la solution approché Y de l'EDO Y 0 = F (Y ), de

condition initiale y(t0 ) = Y 0, associée à une discrétisation uniforme de pas h de l'intervalle [t0 , tf ], en
suivant le modèle suivant :
def EI (F , t0 , t f , Y0 , h ) :
#. . . A completer
return Y # v e c t e u r contenant l e s v a l e u r s a p p r o c h e e s
# o b t e n u e s par l a methode d ' Euler i m p l i c i t e

Exercice 12. On considère ε = 10, x0 = 2, y0 = 0. Compléter votre code pour obtenir la solution
approchée obtenue par la méthode d'Euler implicite dans l'intervalle [0, 30] et avec un pas h = 0.05.
Représenter la solution approchée x ainsi obtenue, en fonction du temps.
Dans un autre graphique, représenter dans le plan (x, y) les trajectoires des solutions approchées
obtenues.

Exercice 13. Refaire l'exercice en utilisant les méthodes programmées auparavant. On pourra aussi
utiliser un solveur d'équations diérentielles ordinaires fournis par python, la fonction odeint du
package [Link] de python. Comparer les résultats obtenus.

11