Table des matières

1 Interpolation - Approximation de fonctions 2

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1 Les 3 grandes classes de méthodes . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Les 3 grandes familles de fonctions approximantes . . . . . . . . . . . . . 3
1.3 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 Interpolation polynomiale : polynôme d’interpolation de Lagrange 3
1.3.2 Interpolation de Newton . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Approximation de fonctions . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.1 Approximation au sens des moindres carrés . . . . . . . . . . . . . 7
1.5 Travaux pratiques sur Matlab . . . . . . . . . . . . . . . . . . . . . . . . 8

2 Équations différentielles ordinaires 11

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1 Proposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Méthodes de résolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.1 Discrétisation du problème aux limites à 1D . . . . . . . . . . . . 12
2.2.2 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.3 Méthode de Runge-Kutta d’ordre 4 . . . . . . . . . . . . . . . . . 13
2.3 Autres schémas numériques . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Passage de l’équation différentielle ordinaire d’ordre 2 à une équation dif-
férentielle ordinaire d’ordre 1 . . . . . . . . . . . . . . . . . . . . . . . . . 14

1
 CHAPITRE 1

Interpolation - Approximation de
fonctions

Introduction
Les polynômes constituent une famille de fonctions tout à fait remarquable en Mathéma-
tiques. Ils sont aussi un outil essentiel du calcul et de l’analyse numérique, notamment
dans l’évaluation ou l’approximation des fonctions, dans les problèmes d’interpolation
et d’extrapolation, dans la résolution d’équations différentielles, etc. L’interpolation est
un processus qui génère de l’information là où elle n’est en principe pas disponible ( les
valeurs de la fonction ne sont connues qu’en certains points ).

L’interpolation doit être distinguée de l’approximation de fonction qui consiste à chercher

la fonction la plus proche possible, selon certains critères, d’une fonction donnée. Dans
le cas de l’approximation il n’est en général plus imposé de passer exactement par les
points donnés initialement. Ceci permet de mieux prendre en compte le cas des erreurs
de mesure, et c’est ainsi que l’exploitation de données expérimentales pour la recherche
de lois empiriques relève plus souvent de la régression linéaire, ou plus généralement de
la méthode des moindres carrés.

1.1 Les 3 grandes classes de méthodes

Il existe trois grandes classes de méthodes :
* l’interpolation
On approche la fonction f par un fonction g appartenant à une classe de fonctions
approximantes (souvent des polynômes) et qui coïncide avec f en n points directs.
* l’approximation aux moindres carrés
On minimise la somme des carrés des erreurs en n points x1 , x2 , ..., xn , c’est-à-dire

2
 1.2. Les 3 grandes familles de fonctions approximantes

n
X
qu’on cherche g qui minimise la quantité : ||f (xi ) − g(xi )||2
i=1

* l’approximation uniforme
On minimise le maximum de l’amplitude de l’erreur entre la fonction f et son
approximation. Les points de collocation s’expriment en fonction des racines des
polynômes de Tchebychev

1.2 Les 3 grandes familles de fonctions approximantes

Les trois grandes familles de fonctions approximantes sont :
* les polynômes (théorème de Stone-Weierstrass)
* les fractions rationnelles (approximants de Padé)
* les fonctions trigonométriques (séries de Fourrier)

1.3 Interpolation
1.3.1 Interpolation polynomiale : polynôme d’interpolation de
Lagrange
Le problème est alors de déterminer l’unique polynôme Pn de degré n, qui satisfait

Pn (xi ) = f (xi ), ∀ i = 0, ..., n. (8)

Le polynôme qui satisfait cette égalité, le polynôme d’interpolation de Lagrange

n
X
Pn (x) = lk (x)f (xk )
k=0

où les lk sont des polynômes de degré n qui forment une base de Pn

n
Y x − xi x − x0 x − xk−1 x − xk+1 x − xn
lk (x) = = ··· ··· (9)
i=0,i6=k
x k − xi xk − x0 xk − xk−1 xk − xk+1 x k − xn

Nous pouvons remarqué que pour chaque i = 1, ..., n, nous avons

n
(
Y x − xi 1 si x = xk ,
= (10)
i=0,i6=k
xk − xi 0 si x = xi

Exemple : Cherchons le polynôme qui passe par les points (1,-3), (3,3) et (4,9).

on a : p(x) = a(x − 3)(x − 4) + b(x − 1)(x − 4) + c(x − 1)(x − 3).

En substituant x = 1, x = 3, x = 4, on a :

3
 1.3. Interpolation

p(1) p(3) p(4)

a= (1−3)(1−4)
, b= (3−1)(3−4)
, c= (4−1)(4−3)

Donc

p(x) = p(1) (x−3)(x−4)

(1−3)(1−4)
+ p(3) (x−1)(x−4)
(3−1)(3−4)
+ p(4) (x−1)(x−3)
(4−1)(4−3)
.

Comme p(1) = −3, p(3) = 3 et p(4) = 9, un petit calcul donne p(x) = x2 − x − 3.

Exercice d’application

Étant donnés 3 points (0,1),(2,5),(4,17) nous allons déterminer le polynôme d’interpola-

tion de Lagrange de degré 2 passant par ces points.

Solution

Calculons l0 , l1 , l2 :

l0 (x) = (x−2)(x−4)
8
, l1 (x) = − x(x−4)
4
, l2 (x) = x(x−2)
8

Le polynôme d’interpolation de Lagrange est alors : P2 (x) = l0 (x) + 5l1 (x) + 17l2 (x) =
1 + x2

L’interpolation polynomiale est une technique d’un ensemble de données ou d’une

fonction par un polynôme. En d’autres termes, étant donné un ensemble de points (ob-
tenu, par exemple, à la suite d’une expérience), on cherche un polynôme qui passe par
tous ces points, et éventuellement vérifie d’autres conditions, de degré si possible le plus
bas.

L’interpolation de Lagrange, peut fort bien diverger même pour des fonctions très régu-
lières (phénomène de Runge)

1.3.2 Interpolation de Newton

Px (x) = a0
+ a1 (x − x0 )
+ a2 (x − x0 )(x − x1 )
+ a3 (x − x0 )(x − x1 )(x − x2 ) (11)
+ ...............................................
+ an−1 (x − x0 )(x − x1 )(x − x2 ) · · · (x − xn−2 )
+ an (x − x0 )(x − x1 )(x − x2 ) · · · (x − xn−2 )(x − xn−1 )

Théorème

4
 1.3. Interpolation

l’unique polynome de degré n passant par les (n+1) points de collocation (xi , f (xi )) pour
i = 0, 1, 2, ..., n peut s’écrire selon la formule de Newton (11) ou encore sous la forme
récursive :

n
X n−1
Y
Pn (x) = f [x0 ] + f [ x0 , ..., xk ]Ak (x) avec Ak (x) = (x − xk ) (12)
k=1 k=0

Les coefficients de ce polynôme sont les différences divisées :

ai = f [ x0 , x1 , x2 , ..., xi ] pour 0≤i≤n (13)

Remarques :

Une fois les coefficients ai connus, on peut évaluer le polynôme de Newton au moyen d’un
algorithme similaire au schéma de Horner. On écrit le polynôme (11) sous la forme :

Pn (x) = a0 +(x−x0 )(a1 +(x−x1 )(a2 +(x−x2 )(a3 +· · ·+(x−xn−2 )(an−1 +an (x−xn−1 )))))
(14)
De cette façon, on réduit le nombre d’opérations nécessaires à l’évaluation du polynôme.
En outre, cette forme est moins sensible aux effets des erreurs d’arrondis

Maintenant, il reste à calculer efficacement la valeur de ce polynôme. La manière la plus

simple consiste à construire une table dite de différences divisées de la façon suivante :

Principe de la méthode

Table des différences divisées

xi f (xi ) f [xi , xi+1 ] f [xi , xi+1 , xi+2 ] f [xi , xi+1 , xi+2 , xi+3 ]

x0 f (x0 )
f [x0 , x1 ]
x1 f (x1 ) f [x0 , x1 , x2 ]
f [x1 , x2 ] f [x0 , x1 , x2 , x3 ]
x2 f (x2 ) f [x1 , x2 , x3 ]
f [x2 , x3 ]
x3 f (x3 )

5
 1.3. Interpolation

La construction de cette table est simple ; on s’est arrêté aux 3ièmes différences divisées.
Les 1ères différences divisées découlent de la définition simple à savoir que :

f (x1 ) − f (x0 )
f [x0 , x1 ] = (15)
x1 − x0

Les deuxièmes différences divisées

f [x1 , x2 ] − f [x0 , x1 ]
f [x0 , x1 , x2 ] = (16)
x2 − x0

Les troisièmes différences divisées :

f [x1 , x2 , x3 ] − f [x0 , x1 , x2 ]
f [x0 , x1 , x2 , x3 ] = (17)
x3 − x0

La formule de Newton utilise la diagonale principale de cette table.

Exercice : Trouver le polynôme d’interpolation de Newton passant par les points (0, 1)
, (1, 2), (2, 9) et (3, 28).

Résolution :

P3 (x) = f (x0 ) + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 ) + f [x0 , x1 , x2 , x3 ](x −

x0 )(x − x1 )(x − x2 )

Établissons la table des différences divisées :

f [x0 , x1 ] = 2−1
1−0
= 1; 9−2
f [x1 , x2 ] = 2−1 = 7;
28−9 7−1
f [x2 , x3 ] = 3−2 = 19; f [x0 , x1 , x2 ] = 2−0 = 3;
f [x1 , x2 , x3 ] = 19−7
3−1
= 6; 6−3
f [x0 , x1 , x2 , x3 ] = 3−0 =1

xi f (xi ) f [xi , xi+1 ] f [xi , xi+1 , xi+2 ] f [xi , xi+1 , xi+2 , xi+3 ]

0 1
1
1 2 3
7 1
2 9 6
19
3 28

6
 1.4. Approximation de fonctions

De cette table, en appliquant :

Pn (x) = 1
+ 1(x − x0 )
+ 3(x − x0 )(x − x1 )
+ 1(x − x0 )(x − x1 )(x − x2 )

on obtient :

P3 (x) = 1 + 1(x − 0) + 3(x − 0)(x − 1) + 1(x − 0)(x − 1)(x − 2)

= x3 + 1

1.4 Approximation de fonctions

1.4.1 Approximation au sens des moindres carrés
On suppose de disposer de n + 1 points (x0 , y0 ), (x1 , y1 ), ......., (xn , yn ).Si le nombre de
données est grand, le polynôme interpolant peut présenter des oscillations importantes.
Pour avoir une meilleur représentation des données on peut chercher un polynôme de
degré m < n qui approche au mieux les données. La notion de meilleur proximité retenue
ici est le critère des moindres carrés.
Il consiste à rechercher le minimum de la fonction
n
X
[Pm (xi ) − yi ]2 (18)
i=0

Le minimum est réalisé pour des valeurs particulières des paramètres qui correspondent
donc à l’ajustement du modèle par apport aux points de données. On dit aussi que cal-
culer Pm c’est effectuer une régression ou un ajustement de nuages de points (linéaire
quand m = 1, quadratique quand m = 2, cubique quand m = 3 et quartique quand
m = 4).

Dans Matlab le polynôme aux moindres carrés peut être calculé toujours avec la com-
mande polyfit : polyfit(x, y, m) où x et y contiennent les valeurs approcher et m est le
degré du polynôme aux moindres carrés.

Exemple : considérons les points (1,1), (2,3), (3,4), (4,3), (5,4) et (6,2). Supposons
que nous souhaitons trouver une approximation des moindres carrés par un polynôme
de dégrée 1. Considérons alors les erreurs

i = |p(xi ) − yi | = |axi + b − yi |, où i = 1, 2, 3, 4, 5, 6,

et choisir a et b afin de minimiser

6
X 6
X
S(a, b) = 2
i = [axi + b − yi ]2 .
i=1 i=1

7
 1.5. Travaux pratiques sur Matlab

Pensons maintenant à S(a, b) comme une fonction de deux variables a et b. Nous

devons avoir alors :

∂S ∂S
∂a
= 0 et ∂b
=0

En substituant (xi , yi ) pour i = 1, 2, 3, 4, 5, 6, nous avons :

91a + 21b = 63

21a + 6b = 17,

enfin a = 1/5 et b = 32/15. par conséquent

p(x) = 15 x + 32
15

1.5 Travaux pratiques sur Matlab

Rappelons qu’un polynôme de degré n, noté Pn , est de la forme :

Pn (x) = a0 + a1 x + a2 x2 + a3 x3 + · · · + an xn (19)

où an 6= 0. Il a au plus n racines distinctes rk , réelles ou imaginaires, telles que Pn (rk ) = 0.

Ses dérivées sont aussi des polynômes, avec Pn (n) (x) = an n! et Pn (n+j) (x) = 0, j ≥ 1,pour
tout x.

Fonctions utilisées en Matlab

roots : calcul de racines à partir des coefficients

poly : calcul des coefficients à partir des racines
conv : produit de polynômes
deconv : division euclidienne des polynômes avec la syntaxe [q,r]=deconv(u,v)
polyval : évaluation du polynôme
polyfit : approximation par un polynôme

Opérations sur les polynômes dans Matlab

Les polynômes sont représentés par des vecteurs lignes dont les composantes sont don-
nées par ordre de puissances décroissantes. Un polynôme de degré n est représenté par
un vecteur de taille n + 1.

Soit p(x) = x3 − 6x2 + 11x − 6

»p = [1 − 6 11 − 6]
»r=roots(p) % calcul des racines
»p=poly(r) % calcul des coefficients du polynôme
Soit A la matrice donnée par A=[1 2 3 ; 4 5 6 ; 7 8 0]

8
 1.5. Travaux pratiques sur Matlab

»p=poly(A) calcul le polynôme caractéristique

»r=roots(p) % calcul des racines de p
»r=eig (A) % valeurs propres de A qui sont les racines du polynôme caractéristique

Produit et division
Soit f et g deux polynômes : Le produit de convolution h de f et g est donné par h =conv
(f, g)

La fonction deconv donne la division de f par g avec la syntaxe : [q, r] = deconv (f, g)
où q est le quotient et r le reste.

Dérivation d’un polynôme

»polyder(p) % calcul de la dérivée de p

Evaluation d’un polynôme

Pour évaluer un polynôme en un point, on utilise la fonction «polyval». Par exemple,

»v=[1 0 2 3] ;

polyval (v, 0 :5)% calcul des ordonnées correspondant aux abscisses 0,1,2,3,4,5

»x=roots(v)

»polyval (v, x(1)) % calcul de l’ordonnée d’une racine

Exercice d’application

Les masses volumiques d’un matériau pour différentes températures sont données ci-
dessous :

Température T en ˚C : 94 205 371

Masse volumique en Kg/m3 : 929 902 860

1. Ecrire la formule d’interpolation de Lagrange qui permet d’interpoler les différents
points des données
2. En utilisant l’interpolation de Lagrange, trouver les masses volumiques du matériau
pour les températures : 251 ˚C ; 351 ˚C et 800 ˚C
3. Utiliser l’interpolation de Newton pour répondre à la question 2.
Exemple
xm = 0 : 10;

p = [1 2 0] ;

9
 1.5. Travaux pratiques sur Matlab

ym = polyval(p, xm ) ;% évaluation du polynôme en xm

% figure (1)

%plot (xm , ym ,0 r0 ) ; % représentation graphique de ym

p1=polyfit(xm , ym , n) % approximation par un polynôme d’ordre n

y=polyval(p1,xm ) % calcul des ordonnées de ce nouveau polynôme en fonction des abs-

cisses

% figure (2)

plot (xm , ym ,0 r0 , xm , y,0 o0 ) % affichage et comparaison.

table= y − ym % dresse le tableau des différentes valeurs

Exercice : Utiliser la commande ’interp1’

x= 0 :10 ;

y= sin(x) ;

xi = 0 : .25 : 10 ;

yi =interp1 (x, y, xi );

plot(x, y,0 o0 , xi , yi ,0 r0 )

10
 CHAPITRE 2

Équations différentielles ordinaires

Introduction
Une équation différentielle est une équation impliquant une ou plusieurs dérivées d’une
fonction inconnue. Si toutes les dérivées sont prises par rapport à une seule variable,
on parle d’équation différentielle ordinaire. Une équation mettant en jeu des dérivées
partielles est appelée équation aux dérivées partielles.
On dit qu’une équation différentielle (ordinaire ou aux dérivées partielles) est d’ordre
p si elle implique des dérivées d’ordre au plus p. Dans le présent chapitre, nous considé-
rons des équations différentielles ordinaires d’ordre un.
On considère un problème dit de Cauchy de la forme suivante :
 0
 y (t) = f (t, y(t)) ∀ t∈I
Trouver y : I ⊂ R −→ R tel que (E) :
y(t0 ) = y0 t0 ∈ I

où f : I × R −→ R une fonction donnée.
Quelques méthodes de résolution de (E) seront présentées avec l’ordre de précision de
chaque méthode.

2.1 Proposition
On rappelle dans la proposition suivante un résultat classique d’analyse.
On suppose que la fonction f est
1. continue par rapport à ses deux variables ;
2. lipschitzienne par rapport à sa deuxième variable, c’est-à-dire qu’il existe une
constante positive L (appelée constante de lipschitz) telle que
|f (t, y1 ) − f (t, y2 )| ≤ L|y1 − y2 |, ∀ t ∈ I, ∀ y1 , y2 ∈ R (20)

11
 2.2. Méthodes de résolution

Alors la solution y = y(t) du problème de Cauchy existe, est unique et appartient à

C 1 (I).

2.2 Méthodes de résolution

2.2.1 Discrétisation du problème aux limites à 1D

b−a
h = tn+1 − tn : (pas constant) =
N
tn = n h; n = 1, ... N

2.2.2 Méthode d’Euler

Principe de la méthode
Trois expressions de récurrences peuvent être obtenues en suivant le méthode d’Euler :

¶ Méthode d’Euler progressif (ou explicite)

yi+1 = yi + hf (xi , yi ), i = 0, 1, 2, ..., n − 1 (21)

· Méthode d’Euler rétrograde (ou implicite) :

yi+1 = yi + hf (xi+1 , yi+1 ), i = 0, 1, 2, ..., n − 1 (22)

¸ Méthode d’Euler semi-implicite

f (xi , yi ) + f (xi+1 , yi+1 )

yi+1 = yi + h , i = 0, 1, 2, ..., n − 1 (23)
2

Exercice d’application
1. Montrer que y(t) = t + e−t est solution de l’équation différentielle :
 0
 y (t) = −y(t) + t + 1

y(0) = 1


2. En utilisant la méthode d’Euler, déterminer y(t)

3. Présenter un tableau de comparaison des 2 solutions (analytique et numérique)
4. Évaluer la précision des résultats.

12
 2.3. Autres schémas numériques

2.2.3 Méthode de Runge-Kutta d’ordre 4

Algorithme
t0 , y0 , h, N donnés.
Pour n ≤ N

k1 = hf (tn , yn )
h k1
k2 = hf (tn + , yn + )
2 2
h k2
k3 = hf (tn + , yn + )
2 2 (24)
k4 = hf (tn + h, yn + k3 )
1
yn+1 = yn + (k1 + 2k2 + 2k3 + k4 )
6
tn+1 = tn + h

Afficher tn+1 et yn+1

Arrêt
Exercice 1 : Résoudre par la méthode de Runge-Kutta d’ordre 4 l’équation différentielle
de l’ordinaire.  0
 y (t) = −y(t) + t + 1

y(0) = 1 (condition initiale)



2.3 Autres schémas numériques

¬ La formule d’Adams-Bashorth (AB3) est un exemple remarquable de méthode à
trois pas (p = 2), du troisième ordre (explicite)

h
un+1 = un + (23fn − 16fn−1 + 5fn−2 ) (25)
12
­ La formule implicite à trois pas et du quatrième ordre d’Adams-Moulton (AM4)

h
un+1 = un + (9fn+1 + 19fn − 5fn−1 + fn−2 ) (26)
24

® la méthode implicite à deux pas du second ordre dite de différentiation rétrograde

(BDF2 de l’anglais backward difference formula)

4 1 2h
un+1 = un − un−1 + fn+1 (27)
3 3 3
ou la méthode BDF3, implicite à trois pas du troisième ordre,

18 9 2 6h
un+1 = un − un−1 + un−2 + fn+1 (28)
11 11 11 11

13
 2.4. Passage de l’équation différentielle ordinaire d’ordre 2 à une équation différentielle
ordinaire d’ordre 1

2.4 Passage de l’équation différentielle ordinaire d’ordre

2 à une équation différentielle ordinaire d’ordre
1
Exercice 1
On considère l’équation de Matthieu amortie :

ÿ + bẏ + a(1 +  cos(t))y = 0

où a, b et  sont des paramètres. On prend comme conditions initiales y(0) = 10−3 et
ẏ = 0. En posant y1 = y et y2 = ẏ
on se ramène à la forme canonique :

 y˙1 = y2
(E) :
y˙2 = −by2 − a(1 +  cos(t))y1


Écrivons une fonction eqdif1.m pour résoudre (E).

function ypoint = eqdif1 (t, y, a, b, epsilon)

ypoint (1,1)= y(2)

ypoint (2,1)= -b*y(2) -a*(1+epsilon*cos(t))*y(1) ;

end

% Paramètres
a=1 ;
b=0.1 ;
epsilon=1 ;

% Temps final
tfinal= 10*pi ;

% Conditions initiales
y01=1e-3 ;
y02=0 ;

[t,y]=ode45(@(t,y))eqdif1(t, y, a, b, epsilon), [0 tfinal], [y01 y02] ;

% Résolution

y1=y( :,1) ; % Extraction de y1 et y2

y2=y( :,2) ;

subplot(221)

14
 2.4. Passage de l’équation différentielle ordinaire d’ordre 2 à une équation différentielle
ordinaire d’ordre 1

plot(t,y1) % y1 fonction du temps

% (représente y(t) pour l’EDO originale)
subplot(222)
plot(t,y2) % y2 fonction du temps
% (représente dy(t)/dt pour l’EDO originale)
subplot(212)
plot(y1,y2) % Plan de phase

Exercice 2
% méthode d’Euler
% Approximation de solution du problème de la valeur initiale
% dy/dt=y-t^2+1 ; 0<=t<=2 ; y(0)=0.5 ;

f=@(t,y) (y-t^2+1) ;
a=input (’Enter left end point, a : ’) ;
b=input (’Enter right end point, b : ’) ;
n=input (’Enter no. of subintervals, n : ’) ;
alpha=input (’Enter the initial condition, alpha : ’) ;

h=(b-a)/n ;
t= a ;
w= alpha ;
fprintf(’ t w\n ’) ;
fprintf(’%5.4f %11.8\n’,t,w) ;

for i= 1 :n
w=w+h*f(t,w) ;
t=a+i*h ;
fprintf(’%5.4f %11.8\n’,t,w) ;
plot(t,w,’r*’) ; grid on ;
xlabel(’t values’) ; vlabel(’w values’) ;
end

15