Pôle Universitaire Léonard de Vinci

Méthodes numériques

TD No. 1 et 2 (durée 3h)

Méthodes de Newton et du point fixe. Ordre de convergence.

Exercice I (Méthode du point fixe)

On considère l’équation ex + 4x2 = 4 dans l’intervalle [0, 1].
(1) Montrer que l’équation admet une solution unique (notée x∗ ) dans [0, 1].
(2) Déterminer une fonction g afin d’approcher x∗ par la méthode de point fixe, i.e. trouver
une fonction itérative g telle que la suite générée par xn+1 = g(xn ) converge vers x∗ .
Justifier votre choix.
(3) On choisit x0 = 0. Donner une estimation théorique du nombre d’itérations suffisantes
afin d’approcher x∗ à 10−5 près.

Exercice II (Méthode de Newton vue comme une méthode du point fixe)

Soit f : R → R la fonction définie par f (x) = x2 − cos x (en unité

radian par la suite).
On considère l’équation f (x) = 0 sur l’intervalle [0, 1]. Soit xn n∈N la suite générée par
l’algorithme de Newton appliqué à cette équation en choisisant x0 = π/4.
(1) Montrer que l’équation admet une seule racine dans l’intervalle [0, 1], notée x∗ par la
suite.
f (x)
(2) On introduit la fonction g : [π/4, 1] → R définie par g(x) = x − 0 . Montrer que
f (x)
g(x) ∈ [π/4, 1] pour tout x ∈ [π/4, 1].
(3) Montrer que la fonction g est strictement contractante sur l’intervalle [π/4, 1].
(4) En déduire que la suite xn n∈N converge vers la racine x∗ . (Indice : on peut considérer

l’algorithme de Newton comme une méthode particulière du point fixe.)

En déduire une estimation d’erreur de |xn − x∗ |. (Remarque : on n’a pas besoin
d’optimiser cette estimation.)

Exercice III (Ordre de convergence)

Déterminer l’ordre de convergence des suites suivantes :
1 n xn 1
(1) xn = sin ; (2) xn = 10−2 ; (3) xn+1 = + , avec x0 > 0.
n 2 xn
(4) Soit(xn )n∈N une suite convergente vers x∗ à l’ordre α.
On suppose que |xn+1 − x∗ | = C|xn − x∗ ||xn−1 − x∗ | pour n >> 1, C 6= 0.
Déterminer la valeur de α.

1
Exercice IV (Prise en main de Python)

(1) Tester les instructions suivantes et comprendre comment on manipule des polynômes :
import numpy as np
A=[ 1, 2, 3]
P1=np.poly1d(A, variable=’x’)
P2=np.poly1d(A, variable=’s’, r=True)
print(P1)
print(P2)
Tester les instructions P1 + P2 , P1 ∗ P2 , P 2/P 1 et P 1 ∗ ∗2 et essayer de comprendre les
opérations.
(2) Calculer (à l’aide de Python) le résultat suivant :
4x3 + 2 + (x − 2)(x + 4)(x3 + 6x + 1).

(3) Tester les fonctions suivantes et comprendre comment on utilise des fonctions :
def saisie_coef(nc):
coef=[]
for i in range(nc):
c=float( input(’coef de degre ’+str(nc-i-1)+’ = ’) )
[Link](c)
return coef

def saisie_polynome():
print(’Saisir un nouveau polynome par coef.’)
deg=int( input(’degre = ? ’) )
Nb_coef=deg+1
coef_P=saisie_coef( Nb_coef )
P=np.poly1d(coef_P)
return P

P=saisie_polynome()
print(P)
Créer le polynôme 4x3 + 6x2 + 1 en utilisant ces fonctions.
(4) On veut générer un polynôme P défini par le résultat de calcul
P (x) = (x − x1 )(x − x2 ) · · · (x − xn ) avec n et x1 , x2 , · · · , xn donnés.
Ecrire un programme similaire à la question (3) afin de réaliser cet objectif.
Tester votre programme avec des cas concrets.

(5) Enfin, on voudrait tracer une courbe sur un intervalle [a, b] donné d’un polynôme choisi.
Pour cet objectif, étudier le programme suivant

2
 import [Link] as plt
import numpy as np

def tracer_courbe(a,b,P):
xmin=min(a,b)
xmax=max(a,b)
Y=[]
Interval=[Link](xmin, xmax, 0.01)
for x in Interval:
[Link]( P(x) )

[Link](Interval, Y, ’b’)
[Link](True)
[Link]()
Tester ce programme au cas où a = −5, b = 5 et P = x2 − 2x + 1.

(6) De façon similaire, tracer la courbe de la fonction sin(x2 + 3x) sur l’intervalle [0, 10].
(Indice : Définir d’abord la nouvelle fonction sin(x2 + 3x).)

3
Pôle Universitaire Léonard de Vinci
Méthodes numériques

TD No.3 et 4 (durée 3h)

Interpolation polynômiale. Méthode de moindres carrés.

Exercice I (Application d’interpolation)

Soit f la fonction définie par f (x) = sin(πx) + cos(2πx). On considère 3 points x0 = 0,
x1 = 1/4 et x2 = 1/2.
n o
(1) Déterminer le polynôme P (x) d’interpolation lié aux 3 points (xi , f (xi ) .
i=0,1,2

(2) En déduire une valeur approchée de sin(π/5) + cos(2π/5).

(3) Donner une estimation des erreurs de la valeur approchée obtenue en question (2) par
rapport à sa valeur exacte, sans calculer cette dernière.
(4) Que ce que l’on peut faire si l’on veut obtenir une valeur approchée plus précise ?

Exercice II (Interpolation de Newton)

Etant donné 10 points suivant
(x0 , y0 ) = (0, 2), (x1 , y1 ) = (1, 6), (x2 , y2 ) = (2, 12), (x3 , y3 ) = (3, 20), (x4 , y4 ) = (4, 30),
(x5 , y5 ) = (5, 42), (x6 , y6 ) = (6, 56), (x7 , y7 ) = (7, 72), (x8 , y8 ) = (8, 90), (x9 , y9 ) = (9, 110).

(1) Déterminer le polynôme d’interpolation lié aux points (noté P (x)).

(2) Supposons qu’il y a une perturbation sur le point (x5 , y5 ) qui devient (5, 42+ε). On note
Pε (x) le polynôme d’interpolation associé au nuage des points perturbés. Déterminer
la valeur de Pε (10).

Exercice III        
−1 0 1 2
Etant donné quatre points : , , , , on cherche des polynômes des
−1 0 1 8
moindres carrés.
(1) Déterminer la parabole (i.e. un polynôme de degré 2) des moindres carrés.
(2) Déterminer le polynôme de degré 3 des moindres carrés.
(3) Calculer le polynôme d’intrepolation de ces 4 points. Comparer à la question (2) : que
l’on peut conclure?

Exercice IV (Application de la méthode de moindres carrés)

Dix adolescents s’exercent à lancer le poids, du bras droit puis du bras gauche. Les résultats
(distances en mètres) obtenus sont les suivants :

4
 Construire un modèle qui permet d’établir un lien linéaire entre les distances de lance par le
bras droit et par le bras gauche. Utiliser le modèle afin de répondre les questions suivantes.
(1) Un adolescent qui lance à 6.5m du bras droit peut espérer lancer à combien de mètres
avec le bras gauche ?
(2) Un adolescent lance le poids à 4.2m du bras gauche. Quelle sera sa performance avec
le bras droit ?

Exercice V (Interpolation en Python)

Soit (xi , yi )i=0,··· ,n un nuage des n + 1 points donnés. Le but est de
— construire le polynôme d’interpolation (de Langrange) associé à ces points ;
— appliquer votre programme dans différents cas,
— et savoir analyser les résultats obtenus.

(1) Construction du programme

• Analyser le programme suivant
import numpy as np

def Li(i,X):
x=np.poly1d([1, 0])
L=1
n=len(X)
for j in range(n):
if j != i:
L=L*(x-X[j])/(float(X[i]-X[j]))
return L

def Interpol_Lagrange(X,Y):
s=0
n=len(X)
for i in range(n):
s=s+Li(i,X)*Y[i]
return s

Les X et Y représentent quels donnés dans la formule mathématique ? Chaque boucle

dans ce programme correspond quelle opération de la formule mathématique ?

5
 • Créer le programme principal qui
— initialise le nuage des points et les stocke dans les tableaux X et Y ;
— détermine le polynôme d’interpolation en appelant la fonction
Interpol Lagrange(X, Y );
— trace sur le même graphe les points saisis et la courbe du polynôme d’interpo-
lation dans un intervalle qui contient tous points saisis.

(2) Application et analyse du cas 1

On considère le nuage des points suivant :
(0, 1), (−1, −1), (1, 5), (3, 19), (4, 29).
• Déterminer le polynôme d’interpolation (noté P (x)) et tracer sa courbe.
• Changer légèrement la coordonnée d’un point (remplacer (0, 1) par (0, 1 + ε), par
exemple), puis déterminer le polynôme d’interpolation (noté Pε (x)). Tracer sa courbe
sur le même graphe. Observer comment la courbe de Pε (x) change en fonction de ε.

(3) Application et analyse du cas 2

Soit (xi , yi )i=0,··· ,n un nuage des n + 1 points donnés par xi = i et yi = i3 .
• Déterminer son polynôme d’interpolation et tracer sa courbe en utilisant le programme
réalisé en (1) et en supposant n = 5.
• Si on choisit n = 20 ? n = 50 ? Comment expliquer les résultats obtenus ?

(4) Pour ceux qui veulent aller plus loin ...

On a vu en cours que l’interpolation polynômiale peut se servir de l’approximation d’une fonc-
tion f . On va montrer ici que la stratégie du choix des points à interpoler a une conséquence
très importante sur la qualité de cette approximation.
On considère la fonction f définie par
1
f (x) = ,
1 + x2
et on voudrait approcher cette f sur l’intervalle [−5, 5] par des polynômes d’interpolation.
(4.A) D’abord, on travaille avec N abscisses d’interpolation équidistantes.

(a) On prend N = 3, i.e. on choisit 3 points équidistants x0 = −5, x1 = 0 et x2 = 5.

Utiliser votre programme afin de
— donner le polynôme d’interpolation dans ce cas ;
— tracer, dans la même figure, la courbe de la fonction f en couleur rouge et la courbe
du polynôme d’interpolation en bleu, tous les deux sur l’intervalle [−5, 5].
(b) Re-faire la question (a) en prenant N = 6 puis N = 11, où N désigne le nombre des
abscisses d’interpolation équidistantes.
(c) En comparant les trois cas étudiés, quelle est votre conclusion ?

6
(4.B) Maintenant, on travaille avec des points de Tchebycheff {xi }i=1,··· ,N définis par
 2i − 1 
xi = 5 cos π , i = 1, · · · , N.
2N
(d) Re-faire la question (a) en prenant respectivement N = 3, N = 6 et N = 11 avec ce
choix de xi .
(e) Quelle amélioration est apportée par le choix des points de Tchebycheff ?

7
Pôle Universitaire Léonard de Vinci
Méthodes numériques

TD No.5 et 6 (durée 3h)

Dérivation numérique

Exercice I (Approximation de dérivée seconde)

Soient f ∈ C ∞ (R), x0 ∈ R et h << 1.
(1) Proposer une formule d’approximation de f 00 (x0 ) en faisant intervenir les 3 points
suivants (x0 , f (x0 )), (x0 − h, f (x0 − h)) et (x0 + h, f (x0 + h)).
(2) Pour cette formule obtenue, établir une estimation de l’erreur.
De même manière, vous pouvez proposer une formule d’approximation améliorée de f 00 (x0 )
en faisant intervenir les 5 points x0 − 2h, x0 − h, x0 , x0 + h et x0 + 2h.

Exercice II (Influence des erreurs d’arrondi)

On veut utiliser la formule centrée
f (x + h) − f (x − h)
f 0 (x) = (1)
2h
afin d’approcher la dérivée d’une fonction f définie sur l’intervalle [a, b] vérifiant
max |f 000 (x)| ≤ M.
x∈[a,b]

On suppose que l’utilisation d’un ordinateur produit une erreur e(x) dans l’évaluation de
f (x), i.e. on a
f (x) = f ∗ (x) + e(x),
où f ∗ (x) représente la valeur effectivement calculée en pratique.
Par conséquent, l’erreur totale due à l’utilisation de la formule centrée (1) (avec f ∗ au lieu
de f ) sera
f ∗ (x + h) − f ∗ (x − h) e(x + h) − e(x − h) h2 000
f 0 (x) − = − f (θ), où θ ∈ [a, b].
2h 2h 6

(1) En supposant que |e(x)| < ε, montrer que la valeur absolue de l’erreur totale commise
est bornée par
ε h2
+ M.
h 6

(2) On considère une fonction dont certaines valeurs sont données ci-dessous.
En vous servant de la formule (1), calculer deux approximations
de f 0 (0.9) en prenant d’abord h = 0.002 et ensuite h = 0.005.
Sachant que la valeur exacte de f 0 (0.9) = 0.62161, calculer les
erreurs commises et expliquer les résultats obtenus.

8
 (3) On suppose que f (x) = sin(x) et que tous les chiffres des approximations de f (x)
du tableau sont significatifs. Déterminer analytiquement la valeur de h qui donne la
meilleure approximation de f 0 (0.9) en utilisant la formule (1). Tester-le numériquement.
(Indice : utiliser la question (1).)

Exercice III
On considère la formule aux différences
f (x + 2h) − 2f (x + h) + 2f (x − h) − f (x − 2h)
A(h) = ' f (3) (x),
2h3
une approximation de la dérivée troisième f (3) (x).
(1) On dispose des valeurs suivantes de la fonction f (x).
En vous servant de la formule aux différences
A(h), calculer deux approximations de f (3) (1).
Estimer numériquement l’ordre de précision
de cette formule aux différences sachant que
f (3) (1) = −1.103 435 06.

(2) En vous servant des développements de Taylor appropriés, montrer que

h2 (5)
A(h) = f (3) (x) + f (x) + O(h4 ),
4
et en déduire l’ordre de précision de l’approximation A(h).
(3) Si pour une certaine fonction f (x) dont on connaı̂t la valeur en certains points, on
n’obtenait pas l’ordre d’approximation obtenu en (2), quelles pourraient en être les
causes ?

Exercice IV (Schéma non centré)

Le tableau suivant donne l’altitude d’un objet en chute mesurée à intervalles de 5 secondes

Trouver une approximation d’ordre 1 de l’accélération en t = 25 secondes.

Exercice V (Expérience numérique en Python)

L’objectif de cet exercice est de réaliser un programme qui permet de calculer la valeur
approchée de la dérivée d’une fonction donnée (par exemple f (x) = cos x) sur un point
donné.
(1) Demander à l’utilisateur d’entrer un point x sur lequel on veut calculer la valeur ap-
prochée de la dérivée de f .

9
(2) Calculer la valeur approchée de f 0 (x) en utilisant le schéma centré à 3 points avec
respectivement
h = 10−i , i = 1, 2, · · · , N.
Tracer la courbe qui représente les valeurs approchées de f 0 (x) en fonction de i.
(Pourquoi pas en fonction de h ?)
Que peut-on observer, quand i = 1, 2, · · · , 6 ? quand i = 12, · · · , 20 ?
D’après votre observation numérique, l’approximation s’améliore quand la valeur de h
devient de plus en plus petit ? Pourquoi ?

(3) Analyse d’erreur : On note ε(h) = f 0 (x) − f (x+h)−f

2h
(x−h)
l’erreur d’approximation en
fonction de valeur de h.
Comment représenter une courbe de ε(h) en fonction de h afin d’apporter une meilleur
visibilité ?
Comment on peut reconnaitre numériquement que le schéma centré à 3 points est un
schéma d’ordre 2 ?
(4) Faire la même chose qu’en (2) et (3) en utilisant le schéma centré à 5 points.

10
Pôle Universitaire Léonard de Vinci
Méthodes numériques

TD No.7 et 8 (durée 3h)

Intégration numérique

Exercice I (Formule de Gauss)

Soit f une fonction régulière. On veut utiliser une formule de quadrature suivante afin
d’approcher l’intégrale Z 1  
f (x) dx ≈ k f (a) + f (b) .
−1
(1) Déterminer les trois constantes a, b et k de façon que la formule de quadrature soit
exacte sur l’espace vectoriel des polynômes de degré le plus élevé possible.
(2) Préciser le degré de précision de cette formule.

Exercice II (Construction de nouvelle formule quadrature)

On propose une formule de quadrature
Z 1
√ 1 2
xf (x) dx ≈ Af + Bf ,
0 3 3
Z 1
√  1  2 
et on note R(f ) l’erreur d’approximation xf (x) dx − Af + Bf .
0 3 3
(1) Déterminer A et B de telle sorte que cette formule de quadrature soit exacte (i.e.
R(f ) = 0) sur un espace vectoriel des polynômes de degré le plus élevé possible.
(2) Préciser le degré de précision de cette formule.

Exercice III (Application concrète)

La navigation inertielle est une technique utilisant des capteurs d’accélération et de rotation
afin de déterminer le mouvement absolu d’un objet. Cet exercice permet d’expliquer le
principe de la navigation inertielle.
Afin de simplifier le calcul, on suppose que le mouvement de l’objet est limité en une ligne
droite, et un référentiel est choisi tel que la position (notée x) et la vitesse (notée v) de l’objet
sont égales à 0 à l’instant t = 0. L’accélération (notée a, unité en m/s2 ) a été observée à
différents moments t (en s) :
t 0 1 2 3 4 5
a 1 2 0 -2 0 3
On suppose que max (|a0 (t)|, |a00 (t)|) ≤ M (M étant une constante connue).
t∈[0,5]

(1) Déterminer des approximations de la vitesse v(ti ) au moment ti = 1, 2, 3, 4, 5.

Déterminer une estimation d’erreur de ces approximations en fonction de M .

11
 (2) Déterminer une approximation de la distance parcourue entre les moments t = 0 et
t = 5.
Déterminer une estimation d’erreur de cette approximation en fonction de M .

Exercice IV (Méthode composite)

Z 2
1
On considère l’intégrale I = dx.
1 x

(1) Evaluer numériquement cette intégrale par la méthode composite des trapèzes pour le
pas h = 13 (on note Iapproche le résultat obtenu).
Et donner une estimation d’erreur de cette approximation, i.e. déterminer M tel que
|I − Iapproche | ≤ M.
(2) Calculer la valeur exacte de I.
(3) Comparer vos résultats des questions (1) et (2) :
• L’erreur “réelle” est inférieure à M ?
• Pourquoi Iapproche obtenue à la question (1) est supérieure à I ?
• Est-ce vrai quelque soit h ?
• Proposer une autre fonction dont la valeur de l’intégrale évaluée par la méthode
des trapèzes est toujours supérieure à la valeur exacte de l’intégrale.
(4) Si on souhaite évaluer I avec la méthode composite de Simpson, quel pas faut-il choisir
pour avoir une erreur inférieure à 10−4 ?

Exercice V (Expérience numérique en Python)

L’objectif de cet exercice est de réaliser un programme qui permet de calculer la valeur
approchée d’un intégrale d’une fonction donnée (par exemple f (x) = cos x) sur un intervalle
donné en utilisant des méthodes composites.
(1) Ecrire une fonction Trapeze(f, a, b, N) où f est la fonction à intégrer, [a, b] l’intervalle
de l’intégrale Ret N le nombre de points d’intégration. Elle doit retourner une valeur
b
approchée de a f (x) dx en utilisant la méthode composite des trapèzes.

(2) Tester votre fonction sur des exemples d’intégrales dont vous pouvez calculer la valeur.

(3) Pour les valeurs N = 21 , 22 , · · · , 210 , déterminer l’erreur (notée Erreur(N )) entre les
valeurs exactes de l’intégrale et celles d’approximation. Tracer Erreur(N ) en fonction
de N à double échelle logarithmique.
Que peut-on observer ? Quelle est votre conclusion au niveau de l’ordre de précision ?

(4) Faire la même chose qu’en (2) et (3) en utilisant la méthode composite de Simpson.

12
Pôle Universitaire Léonard de Vinci
Méthodes numériques

TD No.9 et 10 (durée 3h)

Schémas numériques à 1 pas pour les EDO

Exercice I
y 0 = −y + t + 1,

t ∈ [0, 1],
On considère le problème de Cauchy
y(0) = 1.

(1) Montrer que la fonction f (t, y) (à préciser) est lipschitzienne en y.

En déduire que ce problème admet une solution unique.

(2) Calculer analytiquement cette solution, notée y(t).

(3) Calculer la solution approchée à l’aide du schéma d’Euler avec le pas h = 31 .

(4) Calculer la solution approchée à l’aide du schéma de Taylor d’ordre 2 avec le pas h = 13 .
(Question calculatoire à laisser aux élèves, ne pas traiter en cours.)

(5) Comparer les deux solutions approchées à la solution exacte.

Qu’est-ce que vous observez ? Comment expliquez-vous le phénomène ?

Exercice II (Construction et application du schéma de Crank-Nicolson)

 0
y = f (t, y), t ∈ [a, b],
On considère le problème de Cauchy
y(a) = α.
L’objectif de cet exercice consiste à construire et utiliser (manuellement) un schéma implicite
de Crank-Nicolson.

On suppose que l’on discrétise [a, b] avec des points équidistants, notés ti , et on note h le
pas entre deux points voisins. On voudrait alors approcher les valeurs y(ti ) de la solution.

(1) Montrer que l’on a Z ti+1

y(ti+1 ) − y(ti ) = f (s, y(s)) ds.
ti
En approchant l’intégrale à l’aide du schéma du trapèze, on déduit le schéma de Crank-
Nicolson h 
yi+1 = yi + f (ti , yi ) + f (ti+1 , yi+1 ) .
2
(2) On considère dans cette question une équation différentielle linéaire :
 0
y = −y + t + 1, t ∈ [0, 5],
(A)
y(0) = 1.
On prend h = 1 et on calcule manuellement la solution approchée en utilisant le schéma
de Crank-Nicolson ci-dessus.

13
(3)∗ On considère maintenant une équation différentielle non-linéaire :
 0 2
y = et − ey + 2t, t ∈ [0, 1],
(B)
y(0) = 0.
Quelle est la difficulté rencontrée au cours de la mise en place du schéma ?
Expliquer comment la méthode de point fixe permet de surmonter cette difficulté.
La méthode de point fixe converge-t-elle dans notre cas? sous quelle condition?

Exercice III
On s’intéresse aux approximations de l’EDO (d’ordre 2) suivante
(
x00 (t) + 3x0 (t) + 2x(t) = t2 , t > 0,
(2)
x(0) = 1, x0 (0) = 0.

(1) Ecrire le système (2) comme un système du premier ordre et déterminer la méthode
d’Euler explicite pour calculer les approximations de x(tn+1 ) et x0 (tn+1 ) en fonction de
x(tn ) et x0 (tn ).

(2) En éliminant y, montrer que le système suivant

(
x0 (t) = y(t) − 2x(t), t > 0,
(3)
y 0 (t) = t2 − y(t),

possède la même solution x(t) que (2) si x(0) = 1 et y(0) est bien choisi.
Préciser le bon choix de valeur de y(0).

(3) Appliquer la méthode d’Euler explicite à (3) et donner les formules pour calculer les
approximations de x(tn+1) et y(tn+1 ) en fonction des approximations de x(tn ) et y(tn ).

(4) En utilisant le même paramètre de discrétisation h, calculer les approximations des

2 premiers termes x(t1 ) et x(t2 ) par les 2 formules obtenues en questions (1) et (3).
Montrer que les 2 formules donnent les mêmes résultats d’approximation pour x(t1 ) et
x(t2 ).

(5) Généralisation de la question (4). Montrer que les 2 formules obtenues en questions
(1) et (3) donnent toujours les mêmes approximations de x(tn ) (pour tout n possible),
si le même paramètre de discrétisation h est utilisé.

Exercice IV (Expérience numérique en Python)

L’objectif de cet exercice consiste à
– programmer en Python différents schémas (d’Euler, de Heun, de R.K.4) ;
– les appliquer dans certains exemples d’équation différentielle
– et obtenir une première expérience numérique en faisant d’analyse des résultats.

• Le schéma d’Euler est implémenté en Python par le programme suivant :

14
 import numpy as np

def Euler(a,b, N, alpha, f):

h = (b-a)/float(N)
t = [Link](N+1)
y = [Link](N+1)
t[0] = a
y[0] = alpha
i = 1
for i in range(N):
y[i+1]= y[i]+h*f(t[i], y[i] )
t[i+1]= t[i]+h
return t, y
où on suppose que :
— les données d’entrée de la fonction (intervalle [a, b], nombre de points N , condition
initiale α et la fonction f ) sont connues ;
— la sortie de la fonction est deux tableaux [t, y], où t est le tableau qui sauvegarde
les points discrétisés ti et y est le tableau qui sauvegarde la solution approchée
sur les points ti .
(1) Prendre un exemple de problème de Cauchy (par exemple le pb. (A) ou le pb. (B)).
— Implémenter la fonction f pour cet exemple concret.
— Implémenter la solution exacte (nommée SoluExacte) du problème correspondant
dans une troisième fonction.
(2) Ecrire un programme principal dans lequel
– on initialise des paramètres nécessaires ;
– on calcule la solution approchée en utilisant la fonction “Euler” et la solution
exacte sur les points ti en utilisant la fonction “SoluExacte” ;
– on trace les solutions exacte et approchée en fonction de ti dans la même figure ;
– on trace les erreurs entre les solutions exacte et approchée en fonction de ti dans
une autre figure.
(3) Lancer ce programme principal avec différentes valeurs de N (par exemple 10, 100, 1000).
Observer une corrélation entre les erreurs d’approximation et les valeurs de N .
(Cette observation nous permet de voir que la méthode d’Euler est d’ordre 1.)
(4) Refaire la question (3) avec le schéma de Heun.
Observer une corrélation entre les erreurs d’approximation et les valeurs de N .
(Cette observation nous permet de voir que la méthode de Heun est d’ordre 2.)

Exercice V (Système d’équations différentielles)

On considère le problème de Cauchy
 0 
x (t) = y(t), x(0) = 1,
t ∈ [0, 100], avec une condition initiale
y 0 (t) = −x(t), y(0) = 0.

15
(1) Déterminer la solution analytique du problème.
n o


En déduire la courbe paramétrée définie par x(t), y(t) t ∈ [0, 10] .

(2) On voudrait approcher la solution analytique par la méthode d’Euler.

— Expliciter la méthode d’Euler dans ce cas.
— Appliquer la méthode d’Euler au cas actuel avec h = 1.


On note xh (ti ), yh (ti ) les solutions approchées obtenues.


(3) Tracer la courbe reliant les points xh (ti ), yh (ti ) (pour une valeur h = 1 fixée).
Comparer cette courbe à la courbe analytique obtenue à la question (1).
Expliquer cette différence.

(4) On veut tester l’approximation par la méthode d’Euler pour différentes (petites) valeurs
h. Sachant qu’il faut approcher 2 fonctions inconnues (x(t) et (y(t)), l’implémentation
de la méthode d’Euler doit être adaptée comme suit :
import numpy as np

def Euler_Sys(a,b, N, alpha, F):

h = (b-a)/float(N)
t = [Link](N+1)
y = [Link]((N+1, 2))
t[0] = a
y[0,:] = alpha
i = 1
for i in range(N):
y[i+1,:]= y[i,:]+h*F(t[i], y[i,:] )
t[i+1]= t[i]+h

return y, t

def F(t, y):

val_f=[Link]((1,2))
val_f[:,0] = y[1]
val_f[:,1] = -y[0]
return val_f

Quelles modifications ont faites sur la fonction Euler de l’[Link] ? Pourquoi on modifie
de cette manière ?


Appliquer ces programmes adaptés avec différentes valeurs de h. On note xh (ti ), yh (ti )
les solutions approchées obtenues.


Tracer la courbe reliant les points xh (ti ), yh (ti ) (pour une valeur h fixée).
Que-ce que vous pouvez observer quand on choisit une valeur de h de plus en plus
petite ?

16
Pôle Universitaire Léonard de Vinci
Méthodes numériques

TD No.11 et 12 (durée 3h)

Schémas numériques à pas multiple

Exercice I (Application d’un schéma à multi-pas)

 0
y = −y + t + 1, t ∈ [0, 2],
On considère le problème de Cauchy Le schéma de Nyström
y(0) = 1.
à 2 pas est donné comme suit :
yi+2 = yi + 2hf (ti+1 , yi+1 ).
Calculer la solution approchée à l’aide de ce schéma avec le pas h = 21 .

Exercice II (Construction d’une méthode à 3 pas)

 0
y = f (t, y), t ∈ [a, b],
On considère le problème de Cauchy On discrétise [a, b] avec
y(a) = α.
N +1 points équidistants, notés ti , et on note h le pas entre deux points voisins. On voudrait
construire une méthode explicite à 3 pas de la forme suivante afin d’approcher les valeurs
y(ti ) :
yi+1 = yi + h(af (ti , yi ) + bf (ti−1 , yi−1 ) + cf (ti−2 , yi−2 )), i ≥ 2.
Déterminer les coefficients a, b et c tels que la méthode construite soit d’ordre 3.

Exercice III (Méthodes de prédiction-correction)

 0
y = f (t, y), t ∈ [a, b],
On considère le problème de Cauchy On discrétise [a, b] avec
y(a) = α.
N +1 points équidistants, notés ti , et on note h le pas entre deux points voisins. On voudrait
alors approcher les valeurs y(ti ) de la solution. L’objectif de cet exercice consiste à construire
et étudier une méthode de prédiction-correction.
(1) Sur l’intervalle [ti , ti+1 ], on intégre l’équation différentielle et on obtient
Z ti+1 Z ti+1
0
y(ti+1 ) − y(ti ) = y (s) ds = f (s, y(s)) ds.
ti ti

Montrer que, en approchant f (s, y(s)) par son polynôme d’interpolation correspondant
aux 3 points ti−1 , ti et ti+1 , on a la méthode d’Adams-Moulton à 2 pas comme suit :
(
y0 = α, y1 =h β ; i
h
yi+1 = yi + 12 5f (ti+1 , yi+1 ) + 8f (ti , yi ) − f (ti−1 , yi−1 ) ,
pour i = 1, 2, ..., N − 1, où on suppose que yi est la valeur approchée de y(ti ).

(2) Montrer que l’erreur de troncature de cette méthode est de l’ordre O(h3 ). (Ceci permet
de déduire que la méthode d’Adams-Moulton à 2 pas est une méthode d’ordre 3.)

17
 (•) Pour l’implémenter, il y a deux difficultés dont nous allons discuter par la suite :
— sur le choix de la valeur de β, inconnue jusque là ;
— sur le schéma : il est implicite, car yi+1 se trouve des deux côtés du schéma.

(3) Du fait que le schéma d’Adams-Moulton à 2 pas est implicite, on peut lui associer
un autre schéma explicite afin de prédire la valeur de yi+1 avant de l’utiliser dans le
calcul de f (ti+1 , yi+1 ). Pour ceci, on choisit (pour le moment) le schéma de Heun et on
peut construire une méthode de prédiction-correction avec
— le schéma de Heun pour le prédicteur ;
— le schéma d’Adams-Moulton à 2 pas pour le correcteur.
Expliciter le schéma composé de prédicteur-correcteur ainsi construit, en initialisant
la valeur de β à l’aide du schéma R.K.4.
 0
y = −y + t + 1, t ∈ [0, 5],
(4) On reprend le problème de Cauchy
y(0) = 1.
Appliquer le schéma développé en question (2) au problème ci-dessus avec h = 1.

(•) Remarque : On peut montrer que l’influence du prédicteur est nettement moindre que
celle du correcteur de point de vue de l’erreur de troncature locale si la fonction f (t, y)
est lipschitzienne en y. Ainsi, le correcteur étant d’ordre p, il convient de choisir un
prédicteur d’ordre p − 1 afin que l’ordre global de la méthode de prédiction-correction
soit égal à p.

Exercice IV (Méthodes de prédiction-correction : un autre exemple, à vous de pratiquer!)

 0
y = f (t, y), t ∈ [a, b],
On considère le problème de Cauchy
y(a) = α.
On discrétise [a, b] avec N + 1 points équidistants, notés ti , et on note h le pas entre deux
points voisins. On voudrait alors approcher les valeurs y(ti ) de la solution. L’objectif de cet
exercice consiste à construire une méthode de prédiction-correction et l’appliquer.

(1) Sur l’intervalle [ti , ti+2 ], on intégre l’équation différentielle et on obtient

Z ti+2 Z ti+2
0
y(ti+2 ) − y(ti ) = y (s) ds = f (s, y(s)) ds.
ti ti

Montrer que, en approchant f (s, y(s)) par son polynôme d’interpolation correspondant
aux 3 points ti , ti+1 et ti+2 , on a la méthode suivante :
(
y0 = α, y1 = hβ; i
h
yi+2 = yi + 3 f (ti+2 , yi+2 ) + 4f (ti+1 , yi+1 ) + f (ti , yi ) , i = 0, · · · , N − 2.

où on suppose que yi est la valeur approchée de y(ti ).

(2) Du fait que le schéma obtenu en (1) est implicite, on peut lui associer un autre
schéma explicite afin de prédire la valeur de yi+2 avant de l’utiliser dans le calcul
de f (ti+2 , yi+2 ). Pour ceci, on choisit le schéma d’Euler et on peut construire une
méthode de prédiction-correction :

18
 — pour le prédicteur, on utilise le schéma d’Euler ;
— pour le correcteur, on utilise le schéma obtenu en (1).
Expliciter l’algorithme de prédiction-correction ainsi construit.

(3) Le schéma obtenu en (1) est une méthode d’ordre 4. A partir de cette information,
expliquer pourquoi le choix du schéma d’Euler pour le prédicteur en question (2) n’est
pas un bon choix. (Mais on continue à l’utiliser en question (4) afin de simplifier le
calcul.)

(4) Appliquer la méthode de prédiction-correction obtenue en question (2) au problème

suivant  0
y = −y + t + 1, t ∈ [0, 2],
y(0) = 1,
avec h = 1/2 (et vous pouvez initialiser y1 à l’aide de la méthode d’Euler).

Exercice V (Expérience numérique en Python)

On continue notre expérience avec la méthode de prédiction-correction dans cet exercice.
Pour un exemple, on reprend la méthode construite dans l’[Link].
(1) Programmer (en Python) le schéma de prédicteur-correcteur en une fonction dont
l’appelle doit être comme suit
[t, y] = P redicCorrection(a, b, N, α, β, f )
(2) Ecrire un programme principal dans lequel :
– on initialise les paramètres nécessaires (surtout valeur de β à l’aide du schéma
R.K.4) ;
– on calcule la solution approchée en utilisant la fonction “PredicCorrection” ;
– on trace les solutions exacte et approchée en fonction des ti sur la même figure ;
– on trace les erreurs pour chaque ti entre les solutions exacte et approchée sur une
autre figure.
(3) On reprend le problème dans l’[Link] (Q4).
— Appliquer les programmes développés au problème choisi.
— Lancer le programme avec différentes valeurs de N (nombre de points discrétisés).
— Observer une corrélation entre les erreurs d’approximation et les valeurs de N
afin de confirmer numériquement que l’ordre de la méthode est égal à 3.
(4) Le choix de prédicteur n’est pas unique. On peut alors proposer, par exemple, d’utiliser
le schéma d’Euler pour le prédicteur à la place du schéma de Heun.
Re-faire les questions (1)–(3) avec ce nouveau choix de prédicteur.
Confirmer numériquement que l’ordre du nouveau schéma est égal à 2. Que peut-on
conclure ?

19
Pôle Universitaire Léonard de Vinci
Méthodes numériques

TD No.13 et 14 (durée 3h)

Analyse des erreurs des schémas numériques

Exercice I (Analyse d’erreurs)

Soit f une fonction de classe C 2 , lipschitzienne en y. On considère le problème de Cauchy
 0
y = f (t, y), t ∈ [a, b],
y(a) = α.
On suppose que, comme toujours, on discrétise [a, b] avec des points équidistants, notés ti , et
on note h le pas entre deux points voisins. On propose le schéma du point milieu suivant
afin d’approcher numériquement la solution
 h h 
yi+1 = yi + hf ti + , yi + f (ti , yi ) .
2 2
(1) Montrer que le schéma du point milieu est consistant.

(2) Montrer que le schéma du point milieu est stable.

En déduire que le schéma du point milieu est convergent.

Exercice II
Montrer le théorème suivant

Soit f : [a, b]×R → R une fonction continue et lipschitzienne en y (avec la constante

L). On note y(t) la solution du problème de Cauchy
 0
y (t) = f (t, y(t) ) (t ∈ [a, b])
,
y(a) = α
et yi les valeurs approchées de y(ti ) en utilisant la méthode d’Euler. On suppose
qu’il existe M telle que
∀t ∈ [a, b], |y 00 (t)| ≤ M.
Alors on a une estimation d’erreur suivante :
hM L(ti −a)
|y(ti ) − yi | ≤ [e − 1].
2L

Exercice III
Soit f : [a, b] × R → R une fonction suffisamment régulière et uniformément Lipschitzienne
en y (seconde variable), c’est à dire
∃L > 0, ∀t ∈ [a, b], ∀y1 , y2 ∈ R, on a |f (t, y1 ) − f (t, y2 )| ≤ L|y1 − y2 |.

20
Soit le problème de Cauchy
y 0 (t) = f (t, y(t) ) (t ∈ [a, b])

.
y(a) = α
On pose h = (b − a)/N , où N ∈ N∗ , et ti = a + ih pour 0 ≤ i ≤ N . On travaille avec le
schéma implicite suivant pour approcher la solution y(t) :
h 
yi+1 = yi + f (ti , yi ) + f (ti+1 , yi+1 ) .
2
(1) Montrer que le schéma est bien défini si h < 2/L.

(2) Déterminer l’ordre de ce schéma.

(3) Montrer que, pour tout h < h0 < 2/L assez petit, on a
 hL  ch3
y(ti+1 ) − yi+1 ≤ 1 + y(ti ) − yi + ,
1 − hL/2 1 − hL/2
où c une constante positive.

(4) En déduire que

max y(ti ) − yi ≤ Ch2 .
0≤i≤N

Exercice IV
y 0 = f (t, y),

t ∈ [a, b].
On considère le problème de Cauchy On discrétise [a, b] avec
y(a) = α.
N + 1 points équidistants, notés ti , et on note h le pas entre deux points voisins. Le schéma
d’Adams-Bashforth à 2 pas (noté AB2 ) est donné comme suit
(
y0 = α, y1 = h β; i
yi+1 = yi + h 32 f (ti , yi ) − 12 f (ti−1 , yi−1 ) , 1 ≤ i ≤ N − 1.

De façon similaire pour la stabilité (théorique) des schémas à 1 pas, on peut définir la stabilité
théorique du schéma AB2 comme suit :
Définition 1 (Stabilité théorique ( du schéma AB2 )
z0 = α̃, z1 =h β̃ ;
On considère son problème perturbé i
zi+1 = zi + h 32 f (ti , zi ) − 12 f (ti−1 , zi−1 ) + εi , 1 ≤ i ≤ N − 1.

On dit que le schéma AB2 est théoriquement stable ssi il existe deux constantes M1 et
M2 (indépendantes de h) telles que

max |zi − yi | ≤ M1 max(|α̃ − α|, |β̃ − β|) + M2 max |εi |.

i=0,··· ,N i=1,··· ,N −1

(Remarque : en fonction de référence, cette définition peut être donnée en différentes formes,
qui restent équivalentes entre elles.)
Montrer que le schéma AB2 est théoriquement stable si f est uniformément Lipschitzienne
en y.

21
Exercice V
y 0 = f (t, y),

t ∈ [a, b].
On considère le problème de Cauchy On discrétise [a, b] avec
y(a) = α.
N + 1 points équidistants, notés ti , et on note h le pas entre deux points voisins. Le schéma
d’Adams-Moulton à 2 pas (noté AM2 , vu en TD précédent) est donné comme suit :
(
y0 = α, y1 =h β ; i
h
yi+1 = yi + 12
5f (ti+1 , yi+1 ) + 8f (ti , yi ) − f (ti−1 , yi−1 ) ,

Etudier la stabilité théorique du schéma AM2 en basant sur l’[Link].

(1) Rédiger explicitement la définition de la stabilité théorique du schéma AM2 .

(2) Montrer que le schéma AM2 est théoriquement stable si f est uniformément Lipschitzi-
enne en y.

22
Pôle Universitaire Léonard de Vinci
Méthodes numériques

TD No.15 et 16 (durée 3h)

Méthode de différences finies

Exercice I (Différences finies en 1D)

 00
 u = 16u − 32,
 x ∈ [0, 1],
On considère le problème aux conditions aux limites u(0) = 1, Pb.(A)

 u(1) = −2.

On voudrait approcher la solution u(x) par la méthode des différences finies.

Pour cela, on discrétise l’intervalle [0, 1] avec un pas h = 1/4, et on note xi = i ∗ h avec
i = 0, 1, 2, 3, 4 et ui la valeur approchée de u(xi ).
(1) A partir du schéma centré à 3 points de dérivation numérique de u00 (xi ), proposer un
schéma afin d’approcher l’équation différentielle du problème (A) sur le point xi .

 un système linéaire Au = F , avec A une matrice 3 × 3, F un vecteur et

(2) En déduire

u1
u =  u2 . Déterminer explicitement la matrice A et le vecteur F .
u3

(3) Déterminer les valeurs approchées de u(1/4), u(1/2) et u(3/4).

(4) On passe au cas général où h = N1+1 , avec N un paramètre saisi par l’utilisateur. Re-
faire les questions (1) et (2) en précisant la matrice A et le vecteur F (qui dépendent
de N ).

(5) Une implémentation en Python de la question (4) est donnée ci-dessous. Analyser
cette implémentation et compléter ce qu’il faut afin de

• déterminer numériquement la solution (discrétisée) approchée u(xi ) avec différents

choix du paramètre N ;
• tracer la fonction “solution” et observer le comportement de cette solution en
fonction de N .

Une implémentation en Python :

import numpy as np
import [Link] as plt

### Ex1. (pas pour cas général)

def Diff_Finie_1D(N) :
h = 1.0/(N+1)
u0 = 1

23
 u1 = -2

# construire la matrice A
A = [[0 for _ in range(N)] for _ in range(N)]
for i in range(N) :
A[i][i]=-2/(h*h)-16
if (i<N-1) :
A[i+1][i]=1/(h*h)
A[i][i+1]=1/(h*h)
## print("Vérification : A=", A)

# construire le vecteur F
F = [0 for _ in range(N)]
for i in range(N) :
if(i == 0) :
F[i]= -u0/(h*h)-32
elif(i == N-1) :
F[i]= -u1/(h*h)-32
else :
F[i] = -32
## print("Vérification : F=", F)

# résoudre Au=F, u solution approchée

u = [Link](A, F)
u = [Link](u, 0, u0)
u = [Link](u, N+1, u1)
## print("Vérification : F=", F)
x = [Link](0, 1, N+2)
return x,u

Exercice II (Différences finies en 2D)

On considère le problème aux conditions aux limites
(
∆u(x, y) = x(1 − x) + y(1 − y) dans Ω =]0, 1[×]0, 1[;
u(x, y) = 0 sur ∂Ω,
∂ 2 ∂ 2
où ∆ = ∂x2 + ∂y 2 (l’opérateur laplacien) et où ∂Ω est la frontière du domaine Ω. On peut

montrer que ce problème admet une solution unique, qui est

1
u(x, y) = − x(1 − x)y(1 − y).
2
On s’intéresse à l’approximation numérique du problème dans cet exercice.

24
Pour cela, on discrétise le domaine Ω avec une grille montrée sur la
y
figure ci-contre, i.e. on a N points internes en suivant les directions y4
x et y, avec un pas h = 1/(N + 1). Tous points internes de la grille
y3
possédent donc les coordonnées
(xi , yj ) où xi = i ∗ h, yj = j ∗ h, y2

avec i, j = 1, 2, · · · , N . y1

On montre ci-après comment approcher la solution u(x, y) sur ces y0 x

x0 x1 x2 x3 x4
points internes (xi , yj ).
∂ 2u ∂ 2u
(1) Donner les schémas centrés de dérivation numérique de (xi , yj ) et de (xi , yj )
∂x2 ∂y 2
en utilisant 5 points qui sont (xi−1 , yj ), (xi , yj ), (xi+1 , yj ), (xi , yj−1 ) et (xi , yj+1 ).

(2) On considère l’équation originale sur chaque point (xi , yj ). Montrer que, en remplacant
les termes dérivées par les formules de dérivation numérique correspondante, on obteint
un système linéaire comme suit :
uh (x1 , y1 )
 
 uh (x1 , y2 ) 
 .. 

 . 

 u (x , y ) 
 h 1 N 
 u (x , y ) 
 h 2 1 
 uh (x2 , y2 ) 
 
 .. 
 . 
Aω = F, avec ω =  u (x , y )  .
 
 h 2 N 
 .. 
 . 
..
 
.
 
 
 uh (xN , y1 ) 
 
 uh (xN , y2 ) 
 
 .. 
 . 
uh (xN , yN )
Déterminer explicitement la matrice A et le vecteur F . Ainsi, la solution uh (xi , yj ) est
une approximation de la solution u(xi , yj ).

Avant de passer au cas général, penser à commencer par un cas concret où N = 2 afin
d’obtenir un premier expérience.

(3) (On traite cette question s’il reste de temps.) Implémenter ce résultat en Python et
déterminer numériquement la solution approchée uh (xi , yj ).
Représenter graphiquement la solution approchée (en utilisant le python module “mat-
[Link]” et la fonction “plot surface(x, y, z)”).
Comparer la solution approchée à la solution analytique (en représentant graphique-
ment l’écart u(xi , yj ) − uh (xi , yj )).

25