CPGE PTSI/PT - Sciences Industrielles de l'Ingénieur PT

Ingénierie Numérique Cours

v1.0
Lycée Jean Zay  21 rue Jean Zay  63300 Thiers  Académie de Clermont-Ferrand

Compétences visées :
A3-07 Analyser un algorithme.
Choisir une démarche de résolution d'un problème d'ingénierie numérique ou d'intelligence
C1-03 articielle.
C3-02 Résoudre numériquement une équation ou un système d'équations.
C3-03 Résoudre un problème en utilisant une solution d'intelligence articielle.
D3-03 Eectuer des traitements à partir de données.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 1 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Table des matières

1 Introduction 3
1.1 Objectif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Enjeux de la simulation numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Mise en ÷uvre pratique de la simulation numérique . . . . . . . . . . . . . . . . . . . . 3

2 Intégration et dérivation numériques 4

2.1 Intégration numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Dérivation numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Inuence du bruit de mesure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

3 Résolution d'équations diérentielles 10

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2 Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3 Résolution numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.4 Schémas numériques classiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.5 Cas des équations diérentielles d'ordre supérieur à 1 . . . . . . . . . . . . . . . . . . . 16
3.6 Bibliothèques de calcul numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4 Résolution d'équations de type f (x) = 0 17

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.2 Méthode de la dichotomie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.3 Méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.4 Réexion sur la convergence de ces 2 algorithmes . . . . . . . . . . . . . . . . . . . . . 20
4.5 Bibliothèques de calcul numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5 Résolution de systèmes linéaires 21

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2 Méthode du pivot de Gauss avec pivot partiel . . . . . . . . . . . . . . . . . . . . . . . 21
5.3 Bibliothèques de calcul numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

6 Traitement de chiers de mesures 23

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
6.2 Lecture de chiers et mise en forme des données . . . . . . . . . . . . . . . . . . . . . . 24
6.3 Filtrage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 2 / 28

CPGE PT - S2I Ingénierie Numérique Cours

1 Introduction
1.1 Objectif
Ce cours a pour objectif de lister les diérentes méthodes d'analyse numérique évoquées dans le pro-
gramme de Sciences de l'Ingénieur en PTSI-PT. Elles formeront un premier étage dans la construction
des connaissances en ce domaine, connaissances importantes pour la plupart des ingénieurs qui tra-
vaillent dans le secteur de la simulation numérique notamment. L'ensemble des algorithmes présentés
seront écrit en Python, conformément aux recommandations du programme.
Ce cours n'a en revanche pas vocation à mener une étude approfondie des concepts mathématiques
sous-jacents et des limites de chacune des méthodes abordées. Il pourra donc être utile de se référer à
la littérature spécialisée pour les étudier plus nement.

Attention
Pour faire le lien avec le cours d'Informatique du Tronc Commun, la plupart des fonctions Python
proposées dans ce cours seront codées en Python pur. Seuls quelques exemples feront appels à des
bibliothèques particulières telles que numpy (notamment au Ÿ3).
Cependant, en règle générale, il est bien plus simple et plus ecace d'utiliser les bibliothèques spé-
cialisées (numpy, scipy, etc).

1.2 Enjeux de la simulation numérique

Lors de la phase de conception d'un système, il est long, complexe et coûteux de réaliser un
prototype an de le soumettre à des tests pour valider les choix de conception.
La simulation numérique, qui est une représentation de phénomènes physiques complexes à l'aide
de modèles mathématiques simulés par des opérations faites sur un ordinateur, est une solution plus
économique. En eet, le prototype est virtuel et facilement modiable. La procédure de simulation est :
• rapide à mettre en place (quelques heures) ;
• peut être faite n'importe où et n'importe quand ;
• permet d'avoir une plus-value en orant une meilleure compréhension et interprétation du com-
portement du produit dans des situations parfois diciles voire impossibles à réaliser expérimen-
talement (par exemple la connaissance de la température à l'intérieur d'un matériau).

1.3 Mise en ÷uvre pratique de la simulation numérique

La plupart des modèles mathématiques font intervenir des équations diérentielles linéaires ou non
linéaires : il s'agit donc d'être capable de déterminer l'évolution d'une grandeur physique au cours du
temps mais également en fonction de l'espace.
On procède alors à une discrétisation du problème qui consiste à découper les signaux et le temps
en petits intervalles de manière à simplier les études et on notera pour la suite de ce cours :
• t = [ti ], i ∈ [1, n] la suite des valeurs discrètes du temps en partant de zéro (t0 = 0) ;
• y = [yi = y(ti )] la suite des évaluations à chaque instant discret ti d'une grandeur physique y(t).

À partir de ces dénitions, il sera possible de résoudre des équations diérentielles linéaires selon
diérentes méthodes sachant que les opérations numériques élémentaires sont la dérivation et l'inté-
gration.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 3 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Remarque f (x)
9
Dans la suite de ce cours, nous utiliserons dans
les exemples la fonction polynômiale suivante
7
(voir Figure 1 ci-contre) :
1 3 5
f (x) = x3 − x2 + 4x + 4
8 2
3
Dans les diérents exemples de code Python, on
considèrera comme prédénie la fonction f dénie 1
comme suit : x
0 1 2 3 4 5 6 7 8 9
1 def f(x):
2 return 1/8*x**3 - 3/2*x**2 + 4*x + 4 Figure 1  Fonction polynômiale exemple

2 Intégration et dérivation numériques

2.1 Intégration numérique
2.1.1 Principe
L'intégration numérique consiste à intégrer (de façon approchée) une fonction f (x) sur un intervalle
borné [a, b], c'est-à-dire à calculer l'aire I sous la courbe représentant la fonction, à partir d'un calcul
ou d'une mesure en un nombre ni de points.
Z b n−1
X Z xi+1
I= f (x)dx = f (x)dx
a i=0 xi

L'approche numérique pour l'intégration est nécessaire pour le traitement de données expérimen-
tales ou de résultats de simulation complexes, ou encore pour intégrer des fonctions n'ayant pas de
primitive connue analytiquement. Elle ne permet néanmoins pas d'obtenir des résultats corrects pour
n'importe quelle fonction : en particulier, la fonction ne doit pas présenter de branche innie dans
l'intervalle d'intégration.
La répartition des points en abscisse est généralement uniforme (pas d'échantillonnage constant h),
mais il existe des méthodes à pas variable, ou encore à pas adaptatif.
La précision de l'intégration numérique peut s'améliorer en augmentant le nombre de points n (en
diminuant le pas d'échantillonnage h) ou en augmentant le degré de l'interpolation polynomiale (sous
réserve de bonnes propriétés de continuité de la courbe).

2.1.2 Méthode des rectangles

Dans cette méthode, on calcule l'intégrale numérique en réalisant une somme de surfaces de rec-
tangles. Le domaine d'intégration est découpé en intervalles et on fait comme si la fonction restait
constante sur chaque intervalle.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 4 / 28

CPGE PT - S2I Ingénierie Numérique Cours

On a alors 2 approximations possibles :

n−1 n
X
X
f (x) I=h yi f (x) I=h yi
i=0 • i=1 •

• • • •
• •
yi • • yi • •
• •
• • • •
x x
a xi b a xi b
h h
(a) Méthode des rectangles à gauche (b) Méthode des rectangles à droite

Figure 2  Méthodes des rectangles

Cette méthode est très utilisée pour les intégrations en temps réel. L'intégrale calculée à l'instant t
est en réalité légèrement retardée de h2 (en moyenne), ce qui n'est pas gênant lorsque h est très faible.
La programmation classique de l'intégration numérique d'une courbe dénie par deux tableaux
d'abscisses x et d'ordonnées y s'écrit alors (pour un échantillonnage constant ou non) :
1 # Création des listes, avec a et b les bornes d'intégration et n le nombre de points calculés
2 # On peut aussi obtenir x et y en lisant un fichier de mesures
3 x = np.linspace(a,b,n) # Création d'un vecteur abcisses
4 y = f(x) # Calcul de l'image de x par la fonction f
1 def integration_rectangles_gauche(x,y): 1 def integration_rectangles_droite(x,y):
2 """ 2 """
3 Paramètres : 3 Paramètres :
4 - x : liste des abscisses 4 - x : liste des abscisses
5 - y : liste des ordonnées 5 - y : liste des ordonnées
6 """ 6 """
7 I = 0 7 I = 0
8 for i in range(len(x)-1): 8 for i in range(1,len(x)):
9 I += y[i]*(x[i+1]-x[i]) 9 I += y[i]*(x[i]-x[i-1])
10 return I 10 return I

2.1.3 Méthode du point milieu f (x)

•
La méthode du point milieu (hors programme)
consiste à considérer la fonction constante sur
chaque intervalle de largeur 2h et égale à la valeur • •
•
prise par le point au milieu de l'intervalle (Figure
3). yi • •
La valeur approchée de l'intégrale s'écrit alors •
(pour n impair) : • •
x
a xi b
n−1
2
h
X
I = 2h y2i−1 Figure 3  Méthode du point milieu
i=1

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 5 / 28

CPGE PT - S2I Ingénierie Numérique Cours

La programmation classique de la méthode du point milieu pour une courbe dénie par deux
tableaux d'abscisses x et d'ordonnées y s'écrit alors (pour un échantillonnage constant ou non, avec un
nombre impair de points) :
1 def integration_point_milieu(x,y):
2 """
3 Paramètres :
4 - x : liste des abscisses
5 - y : liste des ordonnées
6 """
7 I = 0
8 for i in range(1,len(x)-1,2):
9 I += y[i]*(x[i+1]-x[i-1])
10 return I

2.1.4 Méthode des trapèzes

Comme son nom l'indique, cette méthode d'in- f (x)
tégration utilise une somme de surfaces de trapèzes, •
de base h.
La valeur approchée de l'intégrale s'écrit alors : • •
yi •
n−1
X yi+1 + yi yi+1 • •
I=h
2 •
i=0
• •
x
La programmation classique de la méthode des a xi xi+1
h b
trapèzes pour une courbe dénie par deux tableaux
d'abscisses x et d'ordonnées y s'écrit alors (pour un Figure 4  Méthode des trapèzes
échantillonnage constant ou non) :
1 def integration_trapezes(x,y):
2 """
3 Paramètres :
4 - x : liste des abscisses
5 - y : liste des ordonnées
6 """
7 I = 0
8 for i in range(len(x)-1):
9 I += (y[i]+y[i+1])*(x[i+1]-x[i])/2
10 return I

2.1.5 Comparaison des diérentes méthodes

La fonction utilisée comme exemple pour ce cours étant un polynôme, il est aisé de calculer la valeur
exacte de son intégrale entre 2 bornes dénies. Dès lors, on peut tracer le graphique de la Figure 5,
représentant l'écart obtenu (en %) en fonction du nombre de points pour chaque méthode numérique.
On constate alors que le choix de la valeur du pas est important et peut avoir de sérieuses réper-
cussions sur la précision du calcul numérique :
• si le pas est trop grand, la Figure 5 indique que l'erreur croît rapidement à mesure que le pas
augmente ;

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 6 / 28

CPGE PT - S2I Ingénierie Numérique Cours

ε%
R. à gauche
16 R. à droite
14 Point milieu
Trapèzes
12
10
8
6
4
2
n
0 5 20 40 60 80 100

Figure 5  Comparaison des méthodes de calcul d'intégrale

• si le pas est trop petit, on risque d'augmenter les erreurs numériques liées à la représentation des
ottants en Python.
Le choix du pas de calcul se fera donc à l'aune de ces deux remarques et sera le plus souvent
guidé par l'expérience, tant qu'il n'est pas contraint par la nature des données recueillies lors d'une
acquisition par exemple (ou le choix sera limité par le nombre de points mesurés).

2.2 Dérivation numérique

2.2.1 Principe
La dérivation numérique permet de poser les bases utiles pour le paragraphe suivant sur l'inté-
gration des équations diérentielles. Elle s'avère par ailleurs très utile pour le traitement de données
expérimentales.
La dérivation numérique consiste à dériver (de façon approchée) une fonction sur un intervalle borné
[a, b], c'est-à-dire à déterminer la pente de la courbe représentant la fonction, à partir d'un calcul ou
d'une mesure en un nombre ni de points.
La répartition des points en abscisse est généralement uniforme (pas d'échantillonnage constant h),
mais il existe des méthodes à pas variable, ou encore à pas adaptatif.

2.2.2 Méthode à 1 pas

L'estimation de la dérivée la plus simple consiste à calculer la pente à partir du point courant et
du point précédent ou suivant (voir Figure 6a). L'estimation de la dérivée au point i peut alors se
calculer par :
• diérence avant : Di = h1 (yi+1 − yi ) (pente de la droite ∆+ ) ;
• diérence arrière : Di = h1 (yi − yi−1 ) (pente de la droite ∆− ).
Bien entendu, lorsqu'il s'agit de dériver une fonction  en temps réel , le point (xi+1 , yi+1 ) n'est
pas connu. La seule méthode possible est donc celle de la diérence arrière.
Notons aussi que le calcul de la dérivée conduit à un tableau de valeurs de dimension n-1.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 7 / 28

CPGE PT - S2I Ingénierie Numérique Cours

f (x) f (x)

yi+1 yi+1
∆
∆+
yi−1 ∆− yi−1
yi yi

x x
xi−1 xi xi+1 xi−1 xi xi+1

(a) Méthodes à 1 pas (b) Méthode à 2 pas

Figure 6  Dérivation numérique

La programmation classique des deux méthodes de dérivation à 1 pas d'une courbe dénie par deux
tableaux d'abscisses x et d'ordonnées y s'écrit alors (pour un échantillonnage constant ou non) :
1 def derivation_1_pas_arriere(x,y): 1 def derivation_1_pas_avant(x,y):
2 """ 2 """
3 Paramètres : 3 Paramètres :
4 - x : liste des abscisses 4 - x : liste des abscisses
5 - y : liste des ordonnées 5 - y : liste des ordonnées
6 """ 6 """
7 deriv = [] 7 deriv = []
8 for i in range(1,len(x)): 8 for i in range(len(x)-1):
9 deriv.append((y[i]-y[i-1])/(x[i 9 deriv.append((y[i+1]-y[i])/(x[i
]-x[i-1])) +1]-x[i]))
10 return deriv 10 return deriv

2.2.3 Méthode à 2 pas

On peut aussi approximer la pente du segment reliant le point précédent au point suivant (voir
Figure 6b). L'estimation de la dérivée au point i peut alors se calculer par : Di = h1 (yi+1 − yi−1 )
(pente de la droite ∆).
Cette méthode faisant appel au point (xi+1 , yi+1 ) , elle ne peut être utilisée pour un traitement de
données en temps réel.
La programmation classique de la méthode de dérivation à 2 pas d'une courbe dénie par deux
tableaux d'abscisses x et d'ordonnées y s'écrit alors (pour un échantillonnage constant ou non) :
1 def derivation_2_pas(x,y):
2 """
3 Paramètres :
4 - x : liste des abscisses
5 - y : liste des ordonnées
6 """
7 deriv = []
8 for i in range(1,len(x)-1):
9 deriv.append((y[i+1]-y[i-1])/(x[i+1]-x[i-1]))
10 return deriv

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 8 / 28

CPGE PT - S2I Ingénierie Numérique Cours

La Figure 7 montre les erreurs constatées entre la valeur exacte de la dérivée et les valeurs obtenues
par les 3 méthodes en chaque point de l'intervalle [1, 9], en fonction du nombre de points de calcul.

f ′ (x) analytique f (x) f ′ (x) 1 pas avant f ′ (x) 1 pas arrière f ′ (x) 2 pas

f ′ (x) f ′ (x)

x x

(a) Points calculés : n = 10 (b) Points calculés : n = 50

Figure 7  Inuence du pas de calcul sur la précision de la dérivée numérique

On constate immédiatement sur la Figure 7a que la précision du calcul numérique de la dérivée

en un point est directement liée aux nombres de points choisis ou disponibles dans le cas d'une mesure.
Ces résultats peuvent donc être très éloignés de la réalité quand les données sont trop peu nombreuses.

2.3 Inuence du bruit de mesure

Lorsque la courbe est issue d'une mesure, elle est généralement entachée d'un léger bruit, qui peut
devenir catastrophique pour l'évaluation de la dérivée (voir Figure 8).

Données brutes Données ltrées

θ(t) θL (t)

t t
θ̇(t) θ̇L (t)

t t

θ̇L (t) θ̇LL (t)

t t

(a) Scénario 1 (b) Scénario 2

Figure 8  Mesure d'un angle θ(t) au cours du temps par un capteur potentiomé-
trique (un lissage donne θL (t)). Calcul numérique de la dérivée (θ̇(t) et θ̇L (t)).

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 9 / 28

CPGE PT - S2I Ingénierie Numérique Cours

En eet, si les points de mesure restent  en moyenne  au voisinage de la valeur mesurée, il existe
des uctuations rapides, à la fréquence d'échantillonnage, entre les points successifs.
Le calcul de la dérivée conduit à déterminer la pente entre deux points successifs, ce qui per-
turbe fortement le signal dérivé et cache les évolutions  lentes  du signal (lentes devant la période
d'échantillonnage).
Deux solutions sont possibles :
• calculer la dérivée sur un temps plus long que le temps d'échantillonnage, par exemple avec une
méthode à 2 pas. Solution adoptée sur la Figure 8a ;
• ltrer (ou lisser) le signal d'origine pour supprimer l'essentiel du bruit (voir Ÿ6.3), puis dériver
(et ltrer/lisser une nouvelle fois le résultat), comme montré sur la Figure 8b.
Dans les deux cas, le signal dérivé sera entaché d'un retard sur le signal d'origine. Pour les traite-
ments diérés, le retard ne pose pas de problème mais, pour les systèmes de commande en temps réel,
il est nécessaire de trouver un compromis entre la qualité du signal dérivé et le retard.

3 Résolution d'équations diérentielles

3.1 Introduction
Une grande part des problèmes scientiques se modélisent par une équation diérentielle dont on
cherche la solution pour dimensionner ou comprendre le phénomène.
Quand les équations sont simples (linéaires d'ordre 1 ou 2) la résolution analytique est aisée, mais
nombre de modélisations conduisent à des équations non linéaires. Ce cours a pour objectif de proposer
des méthodes pour déterminer une solution approchée.

Remarque
Dans ce chapitre, les courbes et tableaux de données correspondent à un problème de mécanique
classique : la chute libre d'une bille dans l'huile. En appliquant le PFD à la bille, soumise à la pesanteur
et à du frottement visqueux, on aboutit à l'équation diérentielle suivante, qui régit l'évolution de la
vitesse de la bille :
dv(t)
m = −f v(t) + mg
dt
Sous forme adimensionnée, cette équation peut s'écrire :
dy(t)
= −y(t)
dt
Cette équation présente une solution analytique y(t) = y0 e−t où y0 est la condition initiale y(0) = y0 .

3.2 Problème de Cauchy

Le problème de Cauchy consiste à trouver les fonctions Y de [0, T ] → RN , telles que :

 dY
= F (t, Y )
dt où t0 ∈ [0, T ] et Y0 ∈ RN sont des données.
 Y (t ) = Y
0 0

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 10 / 28

CPGE PT - S2I Ingénierie Numérique Cours

La plupart des systèmes d'équations diérentielles de tout ordre peuvent se mettre sous cette
forme de système d'équations diérentielles du premier ordre (un exemple sera traité plus loin). Or, le
théorème de Cauchy-Lipschitz montre que pour les équations classiquement abordées en SI, le problème
de Cauchy admet une solution unique dénie sur [0, T ].

3.3 Résolution numérique

La résolution numérique consiste à retrouver la courbe d'évolution de la solution Y (t), en calculant
un certain nombre de points sur un intervalle de temps donné. Il y a donc échantillonnage du temps.
La forme de l'équation diérentielle traduit la loi d'évolution de la grandeur Y : à l'instant t , dY
dt
est la pente de Y (t), donc la valeur de F (t, Y ) permet d'estimer Y (t + ∆t), pour ∆t petit.
D'un point de vue plus mathématique, l'intégration de l'équation diérentielle entre deux pas de
temps de l'échantillonnage ti et ti+1 conduit à :
Z ti+1 Z ti+1 Z ti+1
dY
dt = F (t, Y )dt on en déduit que : Yi+1 = Yi + F (t, Y )dt
ti dt ti ti
Z ti+1
1
L'intégrale sur un pas de temps F (t, Y )dt correspond au calcul de la pente moyenne de
∆t ti
Y (t) sur le petit intervalle ∆t. Le principe des divers schémas numériques consiste à évaluer le plus
précisément possible cette valeur de la pente moyenne pour estimer la valeur de Yi+1 lorsque Yi est
connu.
Connaissant les conditions initiales Y0 , l'intégration numérique devient le calcul d'une suite, où, à
chaque temps ti , un calcul approché permet de calculer une nouvelle valeur de Y en ti+1 .

3.4 Schémas numériques classiques

3.4.1 Le schéma d'Euler explicite
Le schéma d'Euler explicite utilise comme valeur approchée de la pente moyenne sur [ti , ti+1 ], la
valeur de F en ti (en notant ∆t = ti+1 − ti ). Ainsi :
dY Yi+1 − Yi
= F (t, Y ) est approché par = F (ti , Yi ) ⇒ Yi+1 = Yi + ∆tF (ti , Yi )
dt ∆t

La résolution numérique consiste à calculer la y(t)

suite : 1
(
Y0 (condition initiale donnée)
Yi+1 = Yi + ∆tF (ti , Yi )
Yi
La Figure 9 montre la solution exacte ainsi que Yi+1 t
la solution approchée par le schéma d'Euler explicite. ti ti+1
La fonction F (t, Y ), traduisant la pente imposée à
la solution, est représentée par des segments incli- Figure 9  Schéma d'Euler explicite
nés selon une certaine pente. Ainsi, avec le schéma
d'Euler explicite, la pente du segment [ti , ti+1 ] est imposée par la valeur de F au point approché (ti , Yi ).

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 11 / 28

CPGE PT - S2I Ingénierie Numérique Cours

La programmation du schéma d'Euler conduit à dénir la fonction F(t,y) renvoyant la valeur de

F (t, y) et la fonction euler_explicite(F,y0,t) renvoyant le tableau des y(t) (y) pour chaque valeur
du tableau de temps t, avec pour condition initiale y0. Le code suivant est l'application pour le cas
d'étude adimensionnel (voir Ÿ3.1) :
1 import numpy as np
2
3 # Définition de la fonction F
4 def F(t,y):
5 return -y
6
7 # Schéma d'Euler explicite
8 def euler_explicite(F,y0,t):
9 y = np.zeros([t.shape[0],y0.shape[0]])
10 y[0,:]= y0
11 for i in range(t.shape[0]-1):
12 y[i+1,:]= y[i,:]+(t[i+1]-t[i])*F(t[i],y[i,:])
13 return y
14
15 # Résolution
16 t = np.linspace(0,6,100) # Temps
17 y0 = np.array([1]) # Conditions initiales
18 y = euler_explicite(F,y0,t)

La fonction euler_explicite est traitée vectoriellement (grâce à la bibliothèque numpy) de façon

à pouvoir l'utiliser aussi bien sur une équation simple que sur un système de plusieurs équations
diérentielles (voir Ÿ3.5).

Analyse de la convergence du schéma L'évolution de la solution exacte et des points calculés

par le schéma d'Euler explicite pour diérents pas de temps sont fournis Figure 10.

1.5 1
exacte
2.1
1
0.8 2.0
1.5
0.5 1.0
0.6
0.5
0 0.1
0.4
-0.5

-1 0.2

-1.5 0
2 4 6 8 t 1 2 t

Figure 10  Solution exacte et numérique pour le schéma d'Euler explicite

Les courbes montrent que le résultat est d'autant plus précis que le pas de temps est petit et que
le schéma diverge au-delà d'un pas de 2 s.
L'équation de la chute libre dans l'huile est bien connue en physique et en sciences de l'ingénieur : elle
présente une constante de temps de 1 s. Il est donc raisonnable de penser qu'une intégration numérique
ne peut traduire correctement le comportement de l'équation diérentielle que pour des pas de temps
petits devant cette constante de temps, donc pour δt ≪ 1 s.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 12 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Si 1 s < δt < 2 s, le schéma converge mais conduit à des oscillations qui n'ont pas lieu d'être. Pour
δt ≫ 2 s, le schéma ne converge pas.
Cette observation se généralise pour les équations diérentielles quelconques : le schéma d'Euler
explicite présente des conditions de convergence liées aux constantes de temps du modèle.

Performances du schéma Le tableau de la Figure 11 donne l'erreur comme étant le maximum

de l'écart entre la solution exacte et la solution approchée ainsi que le temps de calcul pour diérentes
valeurs du pas de temps.

h 10−1 10−2 10−3 10−4 10−5 10−6

N pas de temps 100 1000 10 000 100 000 1 000 000 10 000 000
Erreur 4,1 · 10−1 4,9 · 10−2 5,0 · 10−3 5,0 · 10−4 5,0 · 10−5 5,0 · 10−6
Temps (s) 5,2 · 10−3 1,6 · 10−2 7,8 · 10−2 5,7 · 10−1 5,0 50

Figure 11  Performances du schéma d'Euler explicite

L'erreur évolue linéairement en fonction du pas de temps : le schéma est dit  d'ordre 1 . Il faut
descendre à un pas de temps très faible avant d'obtenir un niveau d'erreur satisfaisant.
Le temps de calcul évolue lui aussi linéairement en fonction du pas de temps (complexité en O(N ),
car une seule boucle).
Pour améliorer la précision d'un facteur 10000, il faut diminuer d'un facteur 10000 le pas de temps
et donc augmenter d'un facteur 10000 le temps de calcul.

3.4.2 Le schéma d'Euler implicite

Le schéma d'Euler implicite utilise comme valeur approchée de la pente moyenne sur [ti , ti+1 ], la
valeur de F en ti+1 (en notant ∆t = ti+1 − ti ). Ainsi :
dY Yi+1 − Yi
= F (t, Y ) est approché par = F (ti+1 , Yi+1 ) ⇒ Yi+1 = Yi + ∆tF (ti+1 , Yi+1 )
dt ∆t

La résolution numérique consiste à calculer la y(t)

suite : 1
(
Y0 (condition initiale donnée)
Yi
Yi+1 = Yi + ∆tF (ti+1 , Yi+1 )
Yi+1
La Figure 12 montre la solution exacte ainsi que
t
la solution approchée par le schéma d'Euler impli- ti ti+1
cite. La fonction F (t, Y ), traduisant la pente impo-
sée à la solution, est représentée par des segments in- Figure 12  Schéma d'Euler implicite
clinés selon une certaine pente. Ainsi, avec le schéma
d'Euler explicite, la pente du segment [ti , ti+1 ] est imposée par la valeur de F au point approché
(ti+1 , Yi+1 ).
À la diérence du schéma d'Euler explicite, l'expression de Yi+1 devient implicite, c'est-à-dire qu'il
se calcule comme la solution d'une équation. Il convient alors de distinguer trois cas :

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 13 / 28

CPGE PT - S2I Ingénierie Numérique Cours

1. Pour une équation simple, il est parfois possible d'exprimer Yi+1 en fonction des données. Le
calcul est alors aussi simple que pour le schéma d'Euler explicite.
2. Pour un système d'équations linéaires de grande dimension de type F (ti+1 , Y i + 1) = AYi+1 + B ,
où A et B sont une matrice et un vecteur de grande dimension, le calcul de Yi+1 est explicite
mais nécessite de résoudre un système à chaque pas de temps.
3. Pour une équation quelconque, non linéaire, il est nécessaire de faire appel à un algorithme de
recherche de zéro à chaque pas de temps (voir Ÿ4). Le programme Python proposé ci-après couvre
ce cas.
Dans la plupart des cas, une méthode implicite est plus coûteuse en temps de calcul qu'une méthode
explicite, à pas de temps égal. Cependant, lorsque les conditions de convergence du schéma explicite
obligent à un pas très petit, il peut être avantageux d'employer une méthode implicite avec un pas de
temps raisonnable.
La programmation du schéma d'Euler dans le troisième cas (le plus général) conduit à dé-
nir la fonction equation(t,y) renvoyant la valeur de l'équation Yi+1 − Yi − ∆tF (ti+1 , Yi+1 ) à ré-
soudre, la fonction resoudre(equation,y0,eps) renvoyant la solution de l'équation, et la fonction
euler_implicite(F,y0,t) renvoyant le tableau des y(t) pour chaque valeur du tableau de temps t,
avec pour condition initiale y0.
1 import numpy as np
2
3 # Définition de la fonction F
4 def F(t,y):
5 return -y
6
7 # Equation implicite à résoudre pour déterminer Yi +1
8 def equation(F,t,y,dt,yi):
9 return y-yi-dt*F(t,y)
10
11 # Fonction de recherche de zéro par la méthode de la sécante
12 def resoudre(F,t,dt,yi,eps):
13 y = [yi] # La recherche débute à yi
14 y.append(yi+10*eps) # Second point proche de yi
15 i=0
16 while (abs(y[i+1]-y[i])>eps):
17 f1 = equation(F,t,y[i],dt,yi)
18 f2 = equation(F,t,y[i+1],dt,yi)
19 a = (f2-f1)/(y[i+1]-y[i])
20 b = f1-a*y[i]
21 y.append(-b/a)
22 i=i+1
23 return y[-1]
24
25 # Schéma d'Euler implicite
26 def euler_implicite(F,y0,t):
27 y = np.zeros([t.shape[0],y0.shape[0]])
28 y[0,:]= y0
29 for i in range(t.shape[0]-1):
30 y[i+1,:] = resoudre(F,t[i+1],(t[i+1]-t[i]),y[i,:],1e -4)
31 return y
32
33 # Résolution
34 t = np.linspace(0,6,100) # Temps
35 y0 = np.array([1]) # Conditions initiales
36 y = euler_implicite(F,y0,t)

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 14 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Attention
Contrairement à l'intégration par le schéma d'Euler explicite, cet algorithme n'est pas utilisable pour
un système d'équations diérentielles, car la recherche de zéro dans ce cas n'est valable que pour une
équation scalaire.

Analyse de la convergence du schéma L'évolution de la solution exacte et des points calculés

par le schéma d'Euler implicite pour diérents pas de temps est fournie Figure 13.

1 exacte
2.1
2.0
0.8 1.5
1.0
0.5
0.6
0.1

0.4

0.2

0
2 4 6 8 t

Figure 13  Solution exacte et numérique pour le schéma d'Euler implicite

Les courbes montrent que le résultat est d'autant plus précis que le pas de temps est petit et que
le schéma converge pour tous les pas de temps.
Cette observation se généralise pour les équations diérentielles quelconques : le schéma d'Euler
implicite ne présente pas de condition de convergence.

Performances du schéma Le tableau de la Figure 14 donne l'erreur comme étant le maximum

de l'écart entre la solution exacte et la solution approchée ainsi que le temps de calcul pour diérentes
valeurs du pas de temps.

h 10−1 10−2 10−3 10−4 10−5 10−6

N pas de temps 100 1000 10 000 100 000 1 000 000 10 000 000
Erreur 5,9 · 10−1 5,0 · 10−2 5,0 · 10−3 5,0 · 10−4 5,0 · 10−5 5,0 · 10−6
Temps (s) 5,2 · 10−3 3,1 · 10−2 3,7 · 10−1 3,7 37 370

Figure 14  Performances du schéma d'Euler implicite

L'erreur évolue linéairement en fonction du pas de temps : le schéma est dit  d'ordre 1 . Il faut
descendre à un pas de temps très faible avant d'obtenir un niveau d'erreur satisfaisant.
Le temps de calcul évolue lui aussi linéairement en fonction du pas de temps (complexité en O(N ),
car une seule boucle) et s'avère plus long que dans le cas explicite du fait de la résolution de l'équation
implicite.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 15 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Pour améliorer la précision d'un facteur 10000, il faut diminuer d'un facteur 10000 le pas de temps
et donc augmenter d'un facteur 10000 le temps de calcul.

3.5 Cas des équations diérentielles d'ordre supérieur à 1

Les équations diérentielles mises en ÷uvre en physique ou en ingénierie sont généralement d'ordre
supérieur à 1 (souvent d'ordre 2). Il est cependant possible d'exprimer une équation diérentielle
scalaire d'ordre n sous la forme d'un système de n équations diérentielles d'ordre 1.
Soit l'équation diérentielle :
dn y dn−1 y
 
dy
=f t, y, , . . . , n−1
dtn dt dt

En posant un vecteur Y (appelé vecteur d'état) contenant y1 = y et toutes les dérivées contenues
2 n−1
dans f : y2 = dydt1 , y3 = ddt2y et ainsi de suite jusqu'à yn = ddtn−1y , l'équation diérentielle s'exprime sous
la forme d'un système diérentiel :

dy1


 dt = y2
 dy2
 = y3 dY
..
dt
soit = F (t, Y )
 .

 dt
dyn

= f (t, y1 , y2 , . . . , yn )

dt

Avec :    
y1 y2
 y2 y3
   
 . ..
  
 .. et F (t, Y ) = 
  
Y = 
  . 

   
 yn−1   yn 
yn f (t, y1 , y2 , . . . , yn )

Autre solution Une autre solution classique permettant de résoudre cette équation diérentielle
avec la méthode d'Euler explicite est d'appliquer deux fois la dénition.
Yi+1 − Yi
On note Y p la dérivée première. On a Y pi =
h
Y pi+1 − Y pi Yi+2 − 2Yi+1 + Yi
On note Y pp la dérivée seconde. On a Y ppi = = .
h h2
On obtient alors une relation de récurrence directe qui peut être programmée. Ce schéma est un
schéma d'intégration à deux pas qui pourrait s'inscrire dans un cadre plus général des méthodes à
plusieurs pas.

Exemples de mise en équation La Figure 15 présente 3 exemples de mise en équation de pro-

blèmes classiques :
1. charge d'un condensateur : évolution de la tension u(t) dans un circuit RC ;
2. chute libre d'une bille dans l'huile : évolution de l'altitude z(t) ;
3. pendule simple sans frottement : évolution de l'angle θ(t) ;
4. système masse-ressort-amortisseur : évolution de l'abscisse x(t).

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 16 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Ex. Equa. di. Y (t) Y ′ (t) F

1 1
1 u̇(t) = (E − u(t)) - - f (u, t) = (E − u(t))
RC ! ! RC
f z(t) ż(t) f
2 z̈(t) = g − ż(t) Z(t) = Z ′ (t) = F (Z, t) = g − ż(t)
m ż(t) ! F (Z, t) ! m
g θ(t) θ̇(t) g
3 θ̈(t) = − sin θ(t) Θ(t) = Θ′ (t) = F (Θ, t) = − sin θ(t)
l θ̇(t) ! F (Θ, t) ! l
x(t) ẋ(t) 1
4 mẍ(t) + cẋ(t) + kx(t) = f (t) X(t) = X ′ (t) = F (X, t) = (f (t) − cẋ(t) − kx(t))
ẋ(t) F (X, t) m

Figure 15  Quelques exemples de mise en équation

3.6 Bibliothèques de calcul numérique

Python dispose de bibliothèques de calcul numérique optimisées : il n'est donc pas nécessaire de
refaire le programme d'intégration numérique à chaque problème. Sous Python, la bibliothèque à
charger est scipy.integrate. Il est nécessaire de dénir une fonction F admettant comme arguments
le temps t et le vecteur d'état Y, ainsi que les temps de début et de n et un vecteur de conditions
initiales. On peut ajouter explicitement le vecteur de temps pour lequel on souhaite récupérer les
résultats.
Les 2 exemple ci-dessous traitent les cas 1 et 4 de la Figure 15. Ainsi, vous pourrez voir décrit le
code pour des équations diérentielles du premier et du deuxième ordre.
1 """ 1 """
2 Charge d'un condensateur 2 Masse-ressort-amortisseur
3 """ 3 """
4 import scipy.integrate as integrate 4 import scipy.integrate as integrate
5 import numpy as np 5 import numpy as np
6 6
7 # Définition de la fonction f 7 # Définition de la fonction F
8 def f(t,u): 8 def F(t,Y):
9 R = 10 # Résistance 9 m = 1 # Masse
10 C = 0.01 # Capacité 10 c = 0.4 # Amortissement
11 E = 5 # Echelon de tension 11 k = 1 # Raideur
12 return (E-u)/(R*C) 12 F0 = 10 # Echelon de force
13 13 f1 = Y[1]
14 # Résolution 14 f2 = 1/m*(F0-c*Y[1]-k*Y[0])
15 T = 0.8 # Durée de l'étude 15 return array([f1,f2])
16 u0 = np.array([0]) # Conditions initiales 16
17 t = np.linspace(0,T,100) # Temps 17 # Résolution
18 u = integrate.solve_ivp(f,[0,T],u0) 18 T = 20 # Durée de l'étude
19 Y0 = np.array([0,0]) # Conditions initiales
20 t = np.linspace(0,T,100) # Temps
21 z = integrate.solve_ivp(F,[0,T],Y0)

4 Résolution d'équations de type f (x) = 0

4.1 Introduction
Un grand nombre de problèmes de physique ou d'ingénierie débouchent sur des équations non
linéaires qu'il est impossible de résoudre analytiquement. L'outil numérique est alors indispensable

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 17 / 28

CPGE PT - S2I Ingénierie Numérique Cours

pour trouver des solutions approchées de ces équations. Ce chapitre est dédié à la présentation de
deux méthodes classiques pour résoudre des équations (linéaires ou non) à une inconnue, c'est à dire
déterminer x solution d'une équation f (x) = 0, où f est une fonction quelconque.

Remarque
Pour ce chapitre 4, nous utiliserons dans les exemples la fonction polynômiale du début de ce cours,
et on cherchera à résoudre l'équation suivante :
1 3 3 2 1 3 3 2
f (x) = 3 ⇔ x − x + 4x + 4 = 3 ⇔ x − x + 4x + 1 = 0
8 2 8 2
On cherchera une solution sur l'intervalle [2, 6].

On considère donc pour la suite de ce chapitre 4 que la fonction g est dénie :

1 # Fonction g
2 def g(x):
3 return 1/8*x**3 - 3/2*x**2 + 4*x + 1

4.2 Méthode de la dichotomie

La méthode de la dichotomie est la plus simple et convient poiur les fonctions f continues pour
lesquelles est a priori connu un intervalle [a, b] tel que les signes de f (a) et f (b) soient opposés. Par
conséquent, il existe au moins une solution de l'équation f (x) = 0 sur l'intervalle [a, b]. On suppose
par la suite qu'il n'y a qu'une seule solution sur cet intervalle.

g(x) + + 
5 [ • ] Intervalle 1
+  
4 • [ • ] Intervalle 2
+ 
3 [ ] Intervalle 3
2

1 •
0 x
1 2 3 4 5 6 7
-1 •

-2 •

Figure 16  Recherche du zéro par dichotomie

Le principe de la méthode est simple (voir Figure 16 : l'image par f du milieu c de l'intervalle
est calculée, puis le signe est comparé aux signes de f (a) et f (b) pour déterminer si le zéro est sur
l'intervalle [a, c] ou [c, b]. Seul l'intervalle contenant la racine est conservé et le processus est répété
jusqu'à obtenir une précision susante sur la valeur approchée de la solution.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 18 / 28

CPGE PT - S2I Ingénierie Numérique Cours

1 def dichotomie(f,a,b,epsilon):
2 """
3 Paramètres :
4 - f : fonction
5 - a,b : bornes de l'intervalle
6 - epsilon : précision recherchée
7 """
8 while (abs(b-a)) > epsilon:
9 c = (a+b)/2
10 if f(a)*f(c) <= 0:
11 b = c # [a,c] sera le prochain intervalle
12 else :
13 a = c # [c,b] sera le prochain intervalle
14 return (a+b)/2

4.3 Méthode de Newton

La méthode de Newton s'appuie sur la dérivée de la fonction f pour s'approcher progressivement
de la solution. Elle convient pour les fonctions f continues et dérivables. À partir d'une estimation xi
de la solution, la courbe est assimilée à sa tangente en xi pour déterminer une estimation de la solution
xi+1 et ainsi de suite (voir Figure 17).
g(x)
5

4
•
3

1 •
x3 x1
0 1 • 7x
x0 x2
-1
•
-2

Figure 17  Recherche du zéro par la méthode de Newton

L'expression de la tangente en x0 s'écrit : y = f ′ (x)(x − xi ) + f (xi ). Cette droite coupe l'axe des
abscisses en xi+1 = xi − ff′(x
(xi ) , ce qui donne une nouvelle estimation de la solution.
i)

1 def newton(f,df,x0,epsilon):
2 """
3 Paramètres :
4 - f : fonction
5 - df : dérivée de f
6 - x0 : solution initiale
7 - epsilon : précision recherchée
8 """
9 xi = x0
10 xi_plus_1 = x0-f(x0)/df(x0)
11 while abs(xi_plus_1-xi) > epsilon:
12 xi = xi_plus_1
13 xi_plus_1 = xi-f(xi)/df(xi) # Calcul de la nouvelle estimation de la solution
14 return xi_plus_1

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 19 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Cette méthode converge plus rapidement que la méthode de la dichotomie, mais elle nécessite
l'expression de la dérivée de f , qui n'est pas toujours disponible. Dans ce cas, on pourra utiliser la
méthode de la sécante, qui est très similaire à la méthode de Newton, mais qui calcule la pente de la
tangente par un calcul numérique.

4.4 Réexion sur la convergence de ces 2 algorithmes

Ces 2 algorithmes sont toujours opérants sur des fonctions continues et monotones. Cependant, il
existe des contre-exemples. Par exemple, la dichotomie fonctionne sur les fonctions non monotones, à
condition qu'elle ne passe qu'une seule fois par zéro. La méthode de Newton nécessite des critères plus
exigeants :
• la dérivée en un point ne peut être nulle (sinon, il y aura une division par zéro lors du calcul
de la nouvelle solution approchée). Ce cas de gure est néanmoins plutôt rare du fait du calcul
numérique de la dérivée, qui donne rarement 0 ;
• si la valeur de départ est trop éloignée du vrai zéro, la méthode de Newton peut entrer en boucle
innie sans produire d'approximation améliorée ;
• si la pente de la courbe est très forte, du fait de la présence d'un critère d'arrêt sur deux valeurs
successives de la solution, il y aura arrêt sur une solution éloignée de la vraie solution.
La méthode de la dichotomie semble donc plus robuste, mais plus lente à converger. On peut alors
choisir une solution hybride en commençant la recherche par une méthode de la dichotomie et en
anant le résultat par la méthode de Newton.
A titre indicatif, la Figure 18 ci-dessous montre l'évolution de l'erreur en fonction du nombre
d'itérations (en échelle log). On constate l'ecacité supérieure de la méthode de Newton. Les courbes
ne peuvent descendre sous 10−16 du fait de la précision du codage des ottants dans les calculs.
ε
100 Dichotomie
10−2 Newton
10−4
10−6
10−8
10−10
10−12
10−14
10−16 n
6 101 52 102

Figure 18  Comparaison de l'évolution de l'erreur en fonction du nombre d'itérations

4.5 Bibliothèques de calcul numérique

Les 2 codes suivants permettent de mettre en ÷uvre les méthodes de la dichotomie et de Newton
en utilisant la bibliothèque scipy.
1 """ 1 """
2 Méthode de la dichotomie 2 Méthode de Newton
3 """ 3 """
4 from scipy.optimize import bisect 4 from scipy.optimize import newton
5 a,b = 1,6 # Intervalle de recherche 5 x0 = 3 # Point de départ
6 resultat = bisect(f,a,b) 6 resultat = newton(f,x0,fprime=None)

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 20 / 28

CPGE PT - S2I Ingénierie Numérique Cours

5 Résolution de systèmes linéaires

5.1 Introduction
Beaucoup de problèmes de simulation numérique conduisent in ne à la résolution d'un système
linéaire. Un exemple typique est celui d'une équation aux dérivées partielles (décrivant un système
physique, par exemple), dont la transposition numérique peut se faire en remplaçant les dérivées par
des diérences nies (voir Ÿ3.4), ce qui conduit alors généralement à un système d'équations couplées.
Si l'équation de départ est linéaire, le système obtenu l'est également et, si le nombre d'inconnues
n'est pas trop grand (au maximum quelques centaines de variables), ce système linéaire peut être résolu
numériquement de manière ecace par la méthode du pivot de Gauss.

5.2 Méthode du pivot de Gauss avec pivot partiel

Considérons un système linéaire Ax = b, où A est une matrice carrée d'ordre n et b un vecteur
colonne de taille n. Si A est une matrice inversible (ce que nous supposerons dans toute la suite),
il s'agit d'un système de Cramer qui admet donc une unique solution x donnée par x = A−1 b. La
résolution numérique d'un tel système ne se fait quasiment jamais en calculant d'abord l'inverse de A
(problème qui n'est pas plus simple que la résolution du système Ax = b).
La méthode de Gauss repose sur deux remarques simples :
1. si la matrice A est triangulaire supérieure, alors la résolution du système Ax = b est très facile ;
2. On ne change pas la solution d'un système linéaire en eectuant les mêmes opérations élémentaires
sur les lignes de A et de b. Par opération élémentaire, nous entendons ici la permutation de deux
lignes d'indices i et j (notées Li et Lj ), ou le remplacement de la ligne Li par la ligne Li + αLj
(opération terme à terme que l'on notera Li ←− Li + αLj , α étant un réel quelconque).
La méthode de Gauss consiste donc à transformer le système initial Ax = b en un système trian-
gulaire, en eectuant des opérations élémentaires adéquates sur les lignes de A et b , puis à résoudre
simplement le système obtenu.

Transformation du système (élimination de Gauss) Supposons qu'à l'étape j la matrice A

(transformée par les étapes antérieures) ait la structure suivante :
 
a11 · · · · · · · · · · · · a1n
 
 0 a22 · · · · · · · · · a2n 
 .. ... ... .. 

 . . 

A= .
 ..
 
0 ajj · · · ajn 

 . .. .. .. 
 .
 . . . . 

0 · · · 0 anj · · · ann

Une façon simple d'obtenir une matrice de même structure pour l'étape j + 1 consiste à eectuer
les opérations élémentaires :
aij
Li ←− Li − Lj , j+1≪i≪n
ajj
de façon à faire apparaître naturellement des zéros sur la colonne j pour les indices i ∈ {j + 1, . . . , n}.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 21 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Cette méthode n'est pas toujours satisfaisante, car le pivot (le coecient ajj qui apparaît au
dénominateur de l'équation précédente) peut être très petit (ce qui rend alors les calculs très sensibles
aux erreurs d'arrondis), voire même nul (auquel cas la division ne peut bien sûr pas être eectuée). Dans
la méthode de Gauss avec pivot partiel, on cherche donc dans un premier temps l'indice k ∈ {j, . . . , n}
pour lequel |akj | est maximal, puis (si k ̸= j ) les lignes j et k sont permutées avant d'eectuer les
opérations élémentaires énoncées plus haut. La stabilité numérique de l'algorithme est ainsi assurée.
L'inversibilité de la matrice A garantit l'inégalité maxk∈{j,...,n} |akj | > 0.

Remarque Méthode du pivot total

Signalons au passage l'existence de la méthode du pivot total, qui repose sur la sélection du plus
grand coecient (en valeur absolue) parmi les akl , j ≤ k ≤ n, j ≤ l ≤ n et nécessite en outre une
permutation de colonnes ; cette méthode est un peu plus robuste que celle du pivot partiel, mais
le gain est généralement considéré comme faible par rapport à la complexité supplémentaire qu'elle
introduit dans l'algorithme.

La transformation de la matrice A (et du vecteur b, soumis aux mêmes opérations élémentaires)

peut s'implémenter sous la forme suivante :
1 import numpy as np
2 def transforme_gauss(A0,b0):
3 """
4 Paramètres :
5 - A0 : matrice carrée (array numpy)
6 - b0 : vecteur (array numpy)
7 """
8 A,b = A0.copy(),b0.copy()
9 n = np.size(A,0)
10 # Boucle principale
11 for j in range(n-1):
12 # Recherche du pivot partiel (position k)
13 # avec correction pour repasser dans le référentiel de A
14 k = j + argmax (abs(A[j:n,j]))
15 # échange des lignes j et k de b et de A
16 b[[k,j]] = b[[j,k]]
17 A[[k,j],:] = A[[j,k],:]
18 # Transformation élémentaire des lignes j+1 à n
19 for i in range (j+1,n):
20 alpha = A[i,j]/A[j,j]
21 b[i] -= alpha*b[j]
22 A[i ,:] -= alpha*A[j,:]
23 return A,b

Résolution du système transformé Une fois le système linéaire initial transformé en un système
linéaire triangulaire (A étant alors une matrice triangulaire supérieure), la résolution devient immédiate
puisque xn est obtenu directement par :
bn
xn =
ann
et les valeurs de xn−1 , xn−2 , . . ., x1 successivement au moyen de la récurrence (descendante) :
n
!
1 X
x1 = bi − aik xk
aii
k=i+1

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 22 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Cette récurrence s'implémente simplement sous la forme suivante :

1 def solution_triangulaire(A,b):
2 """
3 Paramètres :
4 - A : matrice carrée (array numpy)
5 - b : vecteur (array numpy)
6 """
7 n = np.size(A,0)
8 x = np.zeros(n)
9 for i in range(n-1,-1,-1):
10 s = b[i]
11 for k in range(i+1,n):
12 s -= A[i,k]*x[k]
13 x[i] = s/A[i,i]
14 return x

Résolution du système linéaire initial La combinaison des deux étapes ci-dessus permet donc
de résoudre un système linéaire par la méthode de Gauss avec pivot partiel :
1 def gauss_pivot_partiel(A,b):
2 """
3 Paramètres :
4 - A : matrice carrée (array numpy)
5 - b : vecteur (array numpy)
6 """
7 AA,bb = transforme_gauss(A,b)
8 x = solution_triangulaire(AA ,bb)
9 return x

5.3 Bibliothèques de calcul numérique

Pour résoudre un système linéaire avec Python, on peut utiliser la fonction linalg.solve de numpy.
Il est ainsi possible de vérier le fonctionnement de la fonction gauss_pivot_partiel en comparant
sur l'exemple simple le résultat obtenu à celui donné par la commande linalg.solve(A,b) de numpy
(avec A un array correspondant à une matrice et b un array vecteur), qui renvoie également la solution
du système Ax = b.
>>> A = np.array([[1,2,3],[1,2,1],[3,1,1]],float)
>>> b = np.array([4,5,6],float)
>>> np.linalg.solve(A,b)
array([1.5,2.,-0.5])
>>> gauss_pivot_partiel(A,b)
array([1.5,2.,-0.5])

6 Traitement de chiers de mesures

6.1 Introduction
Très souvent, lors d'acquisition de mesures, on peut récupérer les données sous la forme d'un chier
csv (données séparées par les virgules). De plus, comme nous avons pu le voir dans le cours sur
les capteurs ou au Ÿ2.3, les signaux recueillis sont très souvent bruités, et empêchent un traitement
numérique ecace.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 23 / 28

CPGE PT - S2I Ingénierie Numérique Cours

L'objectif de ce chapitre est donc double : reprendre les méthodes permettant de lire les données
dans un chier, et présenter quelques méthodes de ltrage numérique.

6.2 Lecture de chiers et mise en forme des données

L'objectif n'est pas ici de faire un cours sur la lecture/écriture de chiers en Python. Néanmoins,
nous allons proposer un code simple et transposable qui permet de gérer le traitement d'un chier de
mesures.
Considérons le chier de mesures capsuleuse.csv ci-dessous :
1 =====================================
2 Acquisition capsuleuse
3 =====================================
4 t;omega_e;omega_s;C_e;C_s
5 0,0;3,083881578947369;52,68297697368422;0,0;0,0
6 0,01;2,569901315789474;46,51521381578948;0,0;0,0
7 0,02;2,8268914473684212;46,258223684210535;0,0;0,0
8 0,03;2,3129111842105265;47,28618421052632;0,0;0,0
9 0,04;2,055921052631579;46,51521381578948;0,0;0,0

Le code ci-dessous permet alors d'extraire de ce chier les données sous forme de trois listes de
même dimensions : une pour le temps, et une pour psipoint (ωe ) et thetapoint (ωs ).
1 mon_fichier=open("capsuleuse.csv","r") # Ouverture du fichier
2 for i in range(4): # On élimine les 4 lignes d'en-tête
3 poubelle=mon_fichier.readline()
4 temps=[] # initialisation des listes
5 psiPoint=[]
6 thetaPoint=[]
7 # Pour chaque ligne du fichier
8 for ligne in mon_fichier:
9 # On élimine le retour chariot et les séparateurs de colonne (point virgule)
10 # On remplace aussi les virgules par des points (séparateur décimal en Python)
11 data = ligne.replace(",",".").rstrip("\n\r").split(";")
12 # On récupère la valeur des colonnes temps, omega_e (thetaPoint) et omega_s (psiPoint)
13 temps.append(float(data[0]))
14 psiPoint.append(float(data[1]))
15 thetaPoint.append(float(data[2]))
16 mon_fichier.close() # Fermeture du fichier

6.3 Filtrage
Les capteurs renvoient un signal qui peut être bruité, du fait de leur mode de fonctionnement, de la
quantication, ou d'éléments perturbateurs pendant la transmission du signal. Que ce soit dans un dans
le cadre d'un asservissement numérique (en temps réel) ou dans l'analyse de mesures pré-enregistrées
(post-traitement), le bruit dans le signal fourni par le capteur est toujours néfaste, car il peut dégrader
les performances d'un système ou amener à des erreurs parfois spectaculaires de résolution numérique.
Très souvent, le signal brut issu du capteur est d'abord ltré de manière analogique, avec une
fréquence de coupure liée à la fréquence d'échantillonnage du convertisseur analogique-numérique (on
parle alors de ltre anti-repliement) puis par un ltre numérique sur le signal échantillonné dans le
microcontrôleur.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 24 / 28

CPGE PT - S2I Ingénierie Numérique Cours

6.3.1 Filtre à moyenne glissante

Le ltre à moyenne glissante consiste à calculer pour chaque yi la moyenne des n termes centrés
sur yi . On pourra donc utiliser la méthode suivante :
n−1
2
1 X
yi = yi+k
n
k=− n−1
2

En python, on peut coder cette fonction de la manière suivante (si on doit utiliser ce ltre pour
un traitement de données en temps réel, la moyenne sera calculée sur les n valeurs précédant la valeur
acquise en dernier) :
1 def moyenne_glissante(x,y,n):
2 """
3 Paramètres :
4 - x : abscisses
5 - y : ordonnees
6 - n : nombre de points
7 """
8 if n%2 == 0: # On souhaite un n impair
9 n -= 1
10 a = int((n-1)/2)
11 y_filtres = []
12 for i in range(a,len(y)-a):
13 y_filtres.append(sum(y[i-a:i+a+1])/n) # Calcul de la moyenne
14 return x[a:-a],y_filtres

ω(t) Données brutes

ω(t) n=8

t t

ω(t) n=4
ω(t) n = 32

t t

Figure 19  Inuence du nombre de points utilisés pour la moyenne glissante

6.3.2 Filtre passe-bas du premier ordre

Ce ltre a une fonction de transfert F (p) = YYf(p)
(p)
avec Y le signal d'entrée et Yf le signal ltré.
Pour obtenir l'équation de récurrence, il sut de discrétiser l'équation diérentielle correspondant à
ce ltre en partant de l'analyse de la relation entrée-sortie écrite sous la forme Yf (p) = F (p)Y (p).

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 25 / 28

CPGE PT - S2I Ingénierie Numérique Cours

Pour un ltre du premier ordre, la fonction de transfert est F (p) = 1+τ p ,

1
d'où :

L−1
Yf (p)(1 + τ p) = Y (p) −−→ τ ẏf (t) + yf (t) = y(t)

L'équation diérentielle est discrétisée en évaluant les diérents termes pour ti = iTe et la dérivée
yf −yf
est approchée par ẏf (ti ) = ẏfi = i Te i−1 (voir Ÿ2.2). On obtient alors :

yfi − yfi−1 Te τ
τ + y fi = y i ⇒ y fi = yi + yf
Te τ + Te τ + Te i−1

Pour respecter le théorème de Shannon, il faut prendre fe ≥ 2

τ soit Te ≤ τ2 , on obtient alors :

yfi = 0, 66yi + 0, 33yfi−1

L'implémentation de ce type de ltre est alors très simple, et peut s'écrire comme suit en Python :
1 def passe_bas_1er_ordre(x,y,tau):
2 """
3 Paramètres :
4 - x : abscisses
5 - y : ordonnees
6 - tau : constante de temps (idéalement, Te/2)
7 """
8 y_filtres = [y[0]]
9 for i in range(1,len(y)):
10 Te = x[i]-x[i-1]
11 y_filtres.append(Te/(tau+Te)*y[i]+tau/(tau+Te)*y_filtres[i-1])
12 return x[1:],y_filtres[1:]

La courbe suivante montre l'inuence du choix de la constante de temps du ltre sur la réponse
ltrée obtenue : plus la constante de temps est grande plus la courbe ltrée s'écarte de la courbe
d'origine non ltrée essentiellement en régime transitoire.

ω(t) Données brutes

ω(t) τ = 4Te

t t

ω(t) τ = 2Te ω(t) τ = 8Te

t t

Figure 20  Inuence de la constante de temps d'un ltre passe-bas numérique

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 26 / 28

CPGE PT - S2I Ingénierie Numérique Cours

6.3.3 Filtre passe-bas du deuxième ordre

Le ltre passe-bas du deuxième ordre est plus sélectif que le ltre passe-bas du premier ordre. La
1
fonction de transfert du deuxième ordre est F (p) = 2ξ 1
, d'où :
1+ ω0 p + ω0 2
p2
 
2ξ 1 L−1 1 2ξ
Yf (p) 1 + p + 2 p2 = Y (p) −−→ 2
ÿf (t) + ẏf (t) + yf (t) = y(t)
ω0 ω0 ω0 ω0

On approche au temps iTe , la dérivée première de la même façon que pour le premier ordre (ẏf (ti ) =
yf −yf ẏf −ẏf yf −2yfi−1 +yfi−2
ẏfi = i Te i−1 ) et la dérivée seconde par ÿf (ti ) = ÿfi = i Te i−1 = i Te 2
(voir Ÿ2.2).
L'équation discrétisée est alors :
1 yfi − 2yfi−1 + yfi−2 2ξ yfi − yfi−1
2 + + y fi = y i
ω0 2 Te ω 0 Te

Ce qui donne l'équation de récurrence suivante :

yfi = a0 yi + b0 yfi−1 + b1 yfi−2

avec :
ω0 2 Te 2 1 + ξω0 Te 1
a0 = b0 = 2 b1 = −
1 + 2ξω0 Te + ω0 2 Te 2 1 + 2ξω0 Te + ω0 2 Te 2 1 + 2ξω0 Te + ω0 2 Te 2

Il faut alors choisir la période Te telle que T1e ≥ 2ω0 (théorème de Shannon) : on choisit généralement
un facteur d'amortissement ξ = 1 pour assurer un aaiblissement optimal de −6 dB pour la pulsation
ω0 .
En prenant ω0 = 1
2Te (valeur limite), on obtient :
yfi = 0, 444yi + 0, 667yfi−1 − 0, 111yfi−2

L'implémentation de ce ltre en Python est alors relativement aisée :

1 def passe_bas_2eme_ordre(x,y,xi,omega0):
2 """
3 Paramètres :
4 - x : abscisses
5 - y : ordonnees
6 - xi : coefficient d'amortissement (idéalement, 1)
7 - omega0 : pulsation propre (idéalement, 1/(2Te))
8 """
9 y_filtres = [y[0],y[1]]
10 for i in range(2,len(y)):
11 Te = (x[i]-x[i-2])/2
12 a0 = omega0**2*Te**2/(1+(omega0*Te)**2+2*xi*omega0*Te)
13 b0 = 2*(1+xi*omega0*Te)/(1+(omega0*Te)**2+2*xi*omega0*Te)
14 b1 = -1/(1+(omega0*Te)**2+2*xi*omega0*Te)
15 y_filtres.append(a0*y[i]+b0*y_filtres[i-1]+b1*y_filtres[i-2])
16 return x[2:],y_filtres[2:]

La courbe suivante montre l'inuence du choix de la fréquence de coupure du ltre : on constante

qu'on observe peu de diérence quant à la valeur de ω0 , même si les résultats semblent faussés quand
ω0 est trop petit.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 27 / 28

CPGE PT - S2I Ingénierie Numérique Cours

ω(t) Données brutes

ω(t) ω0 = 1
4Te

t t

ω(t) ω0 = 1 ω(t) ω0 = 1
2Te 8Te

t t

Figure 21  Inuence de la pulsation de coupure d'un ltre passe-bas numérique

Références
[1] A. Caignot, V. Crespel et D. Violeau : Sciences Industrielles de l'Ingénieur - Tout en un -
MP/MP*, PSI/PSI*, PT/PT*. Vuibert, 2022.
[2] A. Caignot, M. Dérumaux, J. Labasque et L. Moisan : Informatique pour Tous - CPGE
scientiques - 1ere et 2eme année. Vuibert, 2015.
[3] E. Billette : Cours d'informatique de tronc commun, 2021. PTSI - Lycée Jean Zay - Thiers.
[4] A. Caignot et M. Dérumeaux : Cours d'ingénierie numérique, 2014. UPSTI.
[5] D. Defauchy : Cours d'ingénierie numérique, 2022. UPSTI.

Emmanuel BIGEARD - https://s2i.bigeard.me

Lycée Jean Zay - Thiers Page 28 / 28