Le TP12

MPSI 831, PCSI 833
Lycée Masséna
TP 12 : Intégration et Analyse Numérique.
1
Méthodes des rectangles à gauche et du point milieu
Il faut savoir écrire ces deux méthodes presque sans réfléchir.
Question 1. Écrire deux fonctions int_rec_g(f,a,b,n) et int_rec_m(f,a,b,n) prenant toutes deux en entrée une
Rb
fonction f : [a, b] → R, les bornes de l’intervalle, et un entier n > 0, et retournant l’approximation de a f (t) dt obtenue
avec la méthode des rectangles à gauche ou la méthode du point milieu, avec l’intervalle [a, b] découpé en n morceaux.
Question 2. Vérifiez que la méthode du point milieu est exacte pour les fonction affines t 7→ at + b, contrairement à
la méthode des rectangles à gauche. Vérifiez que la méthode du point milieu est inexacte pour t 7→ t2 .
>>> print(int_rec_g(lambda
1.0
>>> print(int_rec_m(lambda
1.0
>>> print(int_rec_g(lambda
0.45
>>> print(int_rec_m(lambda
0.5
>>> print(int_rec_m(lambda
0.33499999999999996
2
x:1,0,1,10))
x:1,0,1,10))
x:x,0,1,10))
x:x,0,1,10))
x:x*x,0,1,10))
Intégrale d’une fonction issue d’une série de mesures
On rappelle ici les principales fonctions pour lire un fichier.
fonction
f=open(nom_du_fichier,'r')
f.read()
f.readlines()
f.readline()
desciption
Ouvre le fichier nom_de_fichier (donné sous la forme d’une chaîne de caractères indiquant son emplacement) en lecture (r comme read). Le fichier doit
exister et seule la lecture est autorisée.
Lit tout le fichier d’un coup et le renvoie sous forme de chaîne de caractères (à
ne réserver qu’aux fichiers de taille raisonnable).
Pareil que précédemment, mais le résultat est un tableau de chaîne de caractères, avec un élément du tableau pour une ligne. Attention, le saut de ligne \n
est présent à la fin de la ligne.
Lit une unique ligne du fichier et la renvoie sous forme de chaîne (avec \n
au bout). Le curseur de lecture (virtuel !) est positionné au début de la ligne
suivante.
On se donne un fichier représentant une série de mesures (à récupérer sur le site web). Les lignes du fichier sont
de la forme x;y où x et y sont des flottants. On suppose que ces lignes sont le résultat de mesures, c’est à dire qu’il
existe une fonction f telle que yi = f (xi ) pour chaque ligne xi ; yi . Les abscisses xi sont supposées ordonnées de façon
croissante, par contre elles ne sont pas supposées ordonnées régulièrement : la distance xi+1
R x − xi n’est pas constante.
En supposant que le fichier a ` lignes (x0 ; y0 jusqu’à x`−1 ; y`−1 ), on souhaite approcher x0`−1 f (t) dt, en utilisant la
méthode des trapèzes, c’est à dire :
Svartz
Page 1/4
2014/2015
MPSI 831, PCSI 833
Lycée Massena
Z
x`−1
f (t) dt '
x0
`−2
X
(xk+1 − xk )
k=0
f (xk ) + f (xk+1 )
2
Question 3. Écrire une fonction integrale(fichier) permettant d’automatiser le procédé, la fonction integrale
prenant en entrée un fichier ouvert en lecture. On pourra par exemple utiliser readlines, split pour scinder une
chaîne de caractères autour d’un élément (par exemple, chaine.split(';') fournit un tableau correspondant à la
séparation de la chaîne de caractères chaine autour de ;), et on n’oubliera pas float permettant de convertir une
chaîne de caractères en flottant.
J’obtiens ceci :
>>> f1=open("fichier1.txt",'r') ; integrale(f1)
464.5123984199802
>>> f2=open("fichier2.txt",'r') ; integrale(f2)
0.20764416913028605
Attention : pour ouvrir le fichier fichier.txt il faut préciser le chemin absolu au lycée, qui commence par quelque
chose comme "U:\Documents\" et se termine donc par "fichier.txt".
3
Coefficients dans la méthode de Newton-Cotes
La méthode de Newton-Cotes à N points d’interpolation sur l’intervalle [−1, 1] consiste à faire l’approximation
suivante :
Z
1
f (t) dt ' 2
−1
N
−1
X
ωi f (ξi )
avec
ξi = −1 +
i=0
2i
.
N −1
où les (ωi ) sont choisis de sorte que cette approximation soit exacte pour les fonctions polynomiales de degré au plus
R1
PN −1
N − 1. Par linéarité de l’application f 7→ −1 f (t) dt − 2 i=0 ωi f (ξi ), il suffit d’imposer l’égalité pour les fonctions
polynomiales t 7→ tk pour 0 ≤ k ≤ N − 1. On se ramène donc au système linéaire à N équations et N inconnues
suivant :
Z
N
−1
X
2i
1 1 k
ωi ξik
avec ξi = −1 +
.
t dt =
∀ 0 ≤ k ≤ N − 1,
2 −1
N
−1
i=0
Matriciellement, on a donc à résoudre M X = Y , avec :




M =


1
ξ0
ξ02
..
.
1
ξ1
ξ12
..
.
1
ξ2
ξ22
..
.
ξ0N −1
ξ1N −1
ξ2N −1
···
···
···
···
···
1


ξN −1
2
ξN
−1
..
.









X=


N −1
ξN
−1
ω0
ω1
ω2
..
.
ωN −1










1
Y = 
2


et
R1
dt
R 1−1
t
dt
R 1−1 2
t dt
−1
..
.
R 1 N −1
t
dt
−1








On va résoudre ce système par pivot de Gauss. Il se trouve que la matrice M est mal conditionnée (c’est-à-dire que
les erreurs d’arrondis dues à l’utilisation de flottants nous éloigne énormément de la solution du système), mais fort
heureusement les coefficients de M et du second membre Y sont des rationnels, on va donc pouvoir résoudre notre
système de manière exacte sans crainte des erreurs d’arrondis. On rappelle que pour travailler avec des rationnels, on
procède de la façon suivante :
import fractions as frac
frac.Fraction(1,2) #le rationnel 1/2.
Avant de commencer, récupérez sur le site web l’algorithme complet du pivot de Gauss si vous ne l’avez pas dans
vos fichiers personnels.
Svartz
Page 2/4
2014/2015
MPSI 831, PCSI 833
Lycée Massena
Question 4. Écrire une fonction systeme(N) prenant en entrée un entier N > 1 et retournant la matrice M et le
second membre Y du système ci-dessus. On fera attention à ce que les coefficients soient tous rationnels.
Question 5. Résoudre le système pour N = 2 et N = 3 à l’aide du pivot de Gauss. Vérifiez que,
• pour N = 2, la solution est :
1
ω0 = ω1 =
2
C’est la méthode des trapèzes ;
• pour N = 3, la solution est :
1
2
ω0 = ω2 =
et
ω1 =
6
3
C’est la méthode de Simpson.
Question 6. Vérifiez expérimentalement que les coefficients ωi sont symétriques, c’est à dire que ωN −1−i = ωi . Vérifiez
aussi que pour N ≥ 9, des coefficients ωi négatifs apparaissent.
4
Le phénomène de Runge
On travaille toujours sur l’intervalle [−1, 1]. Pour f une fonction [−1, 1] −→ R, on peut se poser la question de la
R1
convergence de la méthode de Newton-Cotes vers −1 f (t) dt lorsque N tend vers +∞. On va voir qu’il n’en est rien :
il vaut mieux découper l’intervalle en morceaux plus petits plutôt que de chercher à augmenter l’ordre de la méthode
sans diminuer la taille de l’intervalle. On considère pour cela la fonction de Runge
R:
[−1, 1] −→
x
7−→
R
1
1 + 25x2
Question 7. On rappelle que pour tracer le graphe d’une fonction, il suffit de se donner deux tableaux X et Y de même
taille, les éléments de X étant des abscisses, et les éléments de Y les valeurs de la fonction aux points de X, et d’utiliser
le module matplotlib.pyplot :
import matplotlib.pyplot as plt
plt.plot(X,Y)
plt.show() #pour afficher.
On rappelle que la fonction linspace du module numpy produit des flottants régulièrement espacés. Ainsi
import numpy as np
X=np.linspace(-1,1,100)
produit un tableau numpy comprenant 100 points de l’intervalle [−1, 1]. Tracer le graphe de la fonction de Runge.
Question
8. Écrire une fonction int_newton(N) prenant en entrée un entier N > 1, et calculant l’approximation de
R1
R(t)
dt
obtenue avec la méthode de Newton-Cotes à N points. La calculer pour quelques valeurs de N . Attention, la
−1
résolution du système linéaire prend du temps pour N à partir de 50, c’est la difficulté de résoudre exactement avec des
rationnels : comme les coefficients sont des rapports d’entiers qui peuvent grandir assez vite, les opérations élémentaires
(addition, multiplication...) ne se font plus en temps constant contrairement aux opérations sur des flottants.
Question 9. Comparer les résultats obtenus pour plusieurs N (jusqu’à 50) avec les résultats donnés par la méthode
des rectangles (avec l’intervalle [−1, 1] découpé en N morceaux), et la valeur donnée par scipy.integrate.quad, qui
s’utilise comme suit :
import scipy.integrate as sci
sci.quad(R,-1,1) #avec R la fonction ci-dessus
Svartz
Page 3/4
2014/2015
MPSI 831, PCSI 833
Lycée Massena
La méthode de Newton-Cotes à RN points est exacte
polynomiales de degré au plus N − 1. Ainsi,
R 1 sur les fonctions
PN −1
1
elle consiste à remplacer l’intégrale −1 f (t) dt par −1 Pf,N (t) dt = 2 k=0 ωk Pf,N (ξk ), où ξk = −1 + N2k
−1 et :
Pf,N :
[−1, 1] −→
t
7−→
R
N
−1
X
i=0
f (ξi )
N
−1
Y
j=0
t − ξj
ξi − ξj
j6=i
En effet, la fonction Pf,N est une fonction polynomiale 1 de degré au plus N − 1, dont les valeurs en les ξi sont les
f (ξi ).
Question 10. Définir une fonction interp(N) prenant en entrée un entier N > 1 et retournant la fonction PR,N ,
c’est à dire le polynôme interpolateur à N points obtenu pour la fonction R en les ξi = −1 + N2i−1 . Tracer sur un
même graphique la fonction R et quelques PR,N (voir figure suivante). Pour tracer un graphe et lui donner un nom
(label), on ajoute simplement label="nom" à l’usage de plt.plot. Pour positionner la légende comme sur la figure
avant le plt.show(), on utilise plt.legend(loc="lower center"). Vérifiez expérimentalement que bien que les PR,N
coïncident avec la fonction R sur N points, PR,N s’éloigne drastiquement de R en d’autres points.
Figure 1 – Le phénomène de Runge
1. C’est le polynôme interpolateur de Lagrange, que vous verrez sûrement en cours ou en exercice de mathématiques.
Svartz
Page 4/4
2014/2015

Download Report