Cours/TP0 bis: Importation de modules et Analyse

Download Report

Transcript Cours/TP0 bis: Importation de modules et Analyse

5/2
931,932,933,934
Lycée Masséna
Cours/TP0 bis: Importation de modules et Analyse numérique
1
Importer un module (et voir l’aide)
Prenons deux opérations mathématiques classiques : la racine carrée et l’exponentielle. En Python, elles sont
données par sqrt() et exp(). Seulement, ces fonctions nécessitent d’être importées de la bibliothèque standard, car
elles ne se trouvent pas dans le noyau Python. Sans surprise, on les trouve dans un module appelé math. Deux (quatre)
solutions pour les importer :
— importer le module math avec la commande import math. On accède ensuite aux fonctions avec math.sqrt()
et math.exp().
— une variante de la précédente : on peut utiliser un alias pour éviter de taper le nom du module en entier. Par
exemple, on importera le module sous le nom m, avec import math as m. On accède ensuite aux fonctions avec
m.sqrt() et m.exp(). C’est particulièrement utile lorsque les noms des modules sont longs.
— importer spécifiquement les fonctions sqrt et exp du module math, en tapant from math import sqrt,exp.
On accède aux fonctions simplement avec sqrt() et exp().
— une variante de la précédente : importer toutes les fonctions du module math : from math import *. Cette
solution est à éviter pour plusieurs modules : on risque d’avoir plusieurs fonctions avec le même nom. D’une
manière générale, ayez le réflexe de n’importer que ce dont vous avez besoin !
Il y a bien sûr d’autres fonctions usuelles dans le module math. Outre les fonctions trigonométriques et logarithmiques 1 qui ont les noms attendus, on trouve par exemple les définitions des constantes e et pi. Remarquez que
par défaut, Spyder importe les modules math, numpy, et matplotlib. Il est possible de désactiver cette importation
automatique, ce qu’on devrait faire bientôt pour vous forcer à importer les modules, au cas où à l’oral. Comme il existe
un grand nombre de modules, il est hors de question de les connaître par coeur. Python propose une documentation
très détaillée à l’adresse : http://docs.python.org/3/.
On va voir maintenant comment tracer des courbes avec matplotlib résoudre des équations numériques avec
scipy, des systèmes linéaires avec numpy et des équations différentielles avec scipy.
2
Tracé de courbes
Pour tracer des courbes, il faut importer le module matplotlib.pyplot. Voici une manière de tracer le graphe de
la fonction cosinus sur l’intervalle [−π, π] (on importe également le module numpy, bien qu’on n’en ait pas expressément
besoin).
from math import cos,pi
import matplotlib.pyplot as plt
import numpy as np
abscisses=np.linspace(-pi,pi,20) #20 points dans [-pi,pi], régulièrement espacés.
ordonnees=[cos(x) for x in abscisses] #on calcule les ordonnées associées.
plt.plot(abscisses,ordonnees) #on construit le graphe par couples (abscisse,ordonnée)
plt.show() #et on l’affiche !
abscisses et ordonnees sont des tableaux. La fonction linspace de l’environnement numpy permet d’obtenir
facilement des points régulièrement espacés dans un intervalle donnés, utilisez print pour vérifier. La construction
d’ordonnees se fait par compréhension, notion que l’on a vu dans le TP précédent.
1. Le graphe de la fonction cosinus obtenu est un peu anguleux aux alentours de π et −π, corrigez ce problème.
1. La liste des fonctions du module math : https://docs.python.org/3/library/math.html
Svartz
Page 1/4
2014/2015
5/2
931,932,933,934
Lycée Massena
2. Si l’on rajoute un autre graphe avant la commande plt.show(), les deux graphes seront superposés. Affichez
également le graphe de la fonction sinus.
On pourra se reporter au tutoriel officiel de matplotlib.pyplot 2 pour voir comment ajouter des options (grille,
couleurs des graphes, allure du tracé...).
3
Résolution d’équations numériques et module scipy
On a vu la dernière fois comment résoudre une équation de la forme f (x) = 0 avec la méthode de la dichotomie.
Une méthode plus rapide, lorsqu’elle s’applique, est la méthode de Newton. Elle consiste à effectuer un certain nombre
de fois l’itération :
f (x)
x ←− x − 0
f (x)
1. Si vous ne l’avez pas déja fait la dernière fois, codez la méthode de dichotomie pour résoudre une équation
de la forme f (x) = 0 sur un intervalle [a, b] tel que f (a)f (b) < 0. On prendra en paramètre un flottant tel
que l’itération s’arrête lorsque l’intervalle qui encadre le zéro est de longueur inférieure à . Votre fonction doit
renvoyer une approximation du résultat, mais vous pouvez utiliser print pour afficher les valeurs successives
des bornes de l’intervalle.g
def dicho(f,a,b,epsilon):
...
2. Codez la méthode de Newton. On prendra au choix un paramètre n indiquant le nombre d’itérations, ou alors un
flottant tel que l’itération s’arête lorsque la différence entre deux termes consécutifs est inférieure à . La fonction
f 0 sera explicitement fournie sous la forme d’une fonction g. Le préambule sera donc def newton(f,g,x0,n)
ou def newton(f,g,x0,epsilon).
Comparez la vitesse de convergence avec la méthode dichotmique, avec le
√
calcul de 2. On rappelle qu’on peut déclarer une fonction sans la nommer à l’aide du mot-clef lambda. Par
exemple lambda x:x*x-2.
3. Les méthodes de dichotomie et de Newton sont fournies par le module scipy.optimize, et s’utilisent comme
suit (exemple avec la fonction x 7→ x2 − 2).
import scipy.optimize as opt
opt.bisect(lambda: x*x-2,0,5)
opt.newton(lambda: x*x-2,10,lambda x:2*x)
Utilisez les méthodes de scipy avec la fonction exponentielle (qu’on importera du module math) pour trouver
une approximation de e. Remarquez qu’on peut ne pas spécifier la dérivée seconde dans la méthode de Newton
(c’est alors la méthode de la sécante qui est utilisée) ou au contraire préciser en plus la dérivée seconde (c’est
une autre méthode à la convergence encore plus rapide qui est utilisée).
Pour des problèmes à plusieurs dimensions, on peut utiliser la fonction fsolve, par exemple :
import math
def f(variables):
x,y = variables
return([x*y-2*y-2**x,math.log(x)-y-math.cos(x)])
opt.fsolve(f,(0.1,0))
4
Résolution d’équations linéaires
Le but de cette section est de recoder le pivot de Gauss 3 et de voir comment numpy peut le faire aussi bien que
nous. En Python, on représente les matrices par des tableaux de tableaux, les éléments du gros tableau étant les lignes
de la matrice. Par exemple :

1
M = 4
7
2
5
8

3
6 
9
M=[[1,2,3],
[4,5,6],
[7,8,9]]
2. http://matplotlib.org/users/pyplot_tutorial.html
3. C’est fort instructif !
Svartz
Page 2/4
2014/2015
5/2
931,932,933,934
Lycée Massena
On accède à l’élément (i, j) de la matrice par M[i][j] (attention, les éléments sont numérotés depuis 0). Si M est
bien une matrice (chaque ligne a le même nombre d’éléments), on obtient le nombre de lignes avec len(M) et le nombre
de colonnes avec len(M[0]).
On va écrire un algorithme de résolution de système par pivot.
1. Mise en garde : effectuez les opérations suivantes.
>>>
>>>
>>>
>>>
2.
3.
4.
5.
6.
L=[0,0]
M=L
M[0]=1
print(L)
et observez le résultat. Explication : L est ici un objet composé, et l’affectation M=L fait pointer l’alias M vers le
même objet en mémoire que L. Copier intégralement une liste constituée d’éléments de types simples est réalisable
par l’opération : M=L[:]. Si L est-elle même constituée de types composés (en particulier si L représente une
matrice), il faut copier également les éléments des éléments de L. Réalisez une fonction copy_matrice telle que
si M est une matrice, copy_matrice(M) retourne une matrice copiée 4 .
Écrire une fonction echange_lignes prenant en paramètre une matrice M, et deux indices i et j, tels que
echange_lignes(M,i,j) echange les deux lignes i et j de la matrice M.
Écrire une fonction transvection prenant en paramètre une matrice M, deux indices i et j, et un scalaire 5 mu,
tels que transvection(M,i,j,mu) réalise la transvection Li ← Li + µLj , où Li et Lj sont les lignes d’indices i
et j de M .
On va réaliser un algorithme de pivot sur les lignes. Pour plus de stabilité, au lieu de choisir le premier élément
non nul d’une colonne pour réaliser le pivot, on prend l’élément de la colonne ayant la plus grand valeur
absolue. Écrire une fonction choisir_pivot prenant en entrée une matrice M et un indice de colonne j tels que
choisir_pivot(M,j) renvoie l’indice i de la ligne de M tel que M[i][j] soit maximal en valeur absolue.
Écrire enfin la fonction principale : resolution, prenant en entrée une matrice M et une matrice colonne Y tels
que resolution(M,Y) renvoie une matrice colonne X, vérifiant M X = Y . On supposera que le système est de
Cramer. On effectuera dans l’ordre :
— une copie des matrices M et Y.
— un pivot sur les lignes de M (en faisant les mêmes opérations sur Y).
— et on effectuera la phase de remontée pour trouver les composantes de X
Évidemment, on ne recode par l’algorithme du pivot dés que l’on veut résoudre un système. On utilise plutôt le
module numpy.linalg. La syntaxe est la suivante :
import numpy as np
np.linalg.solve(M,Y)
Remarquez que numpy utilise une structure nommée array qui permet de faire de l’algèbre linéaire 6 :
>>> A=np.array([0,1]) ; B=np.array([1,2])
>>> print(3*A,A+B)
[0,3], [1,3]
Et on déclare de même des matrices (np.array(M), ou M est une liste de listes).
5
Résolution d’équations différentielles
Il existe un grand nombre de méthode d’approximations, on va refaire simplement la méthode d’Euler explicite,
avant d’utiliser le module scipy. On considère une équation différentielle de la forme
0
x (t) = f (x(t), t)
x(a) = x0
En supposant qu’il existe une solution sur l’intervalle [a, b], on essaie d’approcher la solution sur cet intervalle en
se donnant une subdivision (régulière) de l’intervalle [a, b] en un tableau T = [t0 = a < · · · tn−1 = b] et en calculant
des approximations de x(t1 ), x(t2 ).... en posant :
xi = xi−1 + f (xi−1 , ti−1 )(ti − ti−1 ) pour i ∈ {1, . . . , n − 1}
4. Là on le fait à la main, mais il faut savoir qu’il existe un module copy qui fait ça très bien.
5. Pourquoi pas lambda, comme scalaire ? :)
6. Qu’on reverra plus tard avec les 3/2, quand ils auront fait un peu d’algèbre.
Svartz
Page 3/4
2014/2015
5/2
931,932,933,934
Lycée Massena
1. Implémenter la méthode d’Euler sous la forme d’une fonction euler prenant en paramètre une fonction f, un
flottant x0 et un tableau T tels que euler(f,x0,T) retourne le tableau des xi , pour i ∈ {0, . . . , n − 1}.
Tracez sur un même graphique la fonction exponentielle sur l’intervalle [0, 2] et une approximation avec votre
méthode (par exemple pour un tableau T de taille 20 donné par la subdivision régulière de l’intervalle [0,2], on
pourra utiliser linspace du module numpy).
2. On va maintenant utiliser le module scipy.integrate, et plus particulièrement la fonction odeint 7 . La syntaxe
est la même que celle proposée pour la méthode d’Euler 8 .
from scipy.integrate import odeint
import numpy as np
def f(y,t): #la fonction qui mène à l’exponentielle
return y
t=np.linspace(0, 2, 20)
x=odeint(f,1,t) #le tableau des xi.
Bien sûr, il est possible de résoudre des équations différentielles où l’inconnue est une fonction vectorielle. Voici
un exemple de pendule amorti.
def pendule(x,t):
return (x[1],-3*x[0] - x[1])
t=np.linspace(0,10,100)
x=odeint(pendule,(1, 0),t) #(1,0): conditions initiales.
z=[u[0] for u in x] #on ne prend que x, pas x’
plt.plot(t,z)
plt.show()
3. Résoudre à la main l’équation différentielle suivante :
 00
 x (t) = x(t)0 + 6x(t)
x(0) = 2
 0
x (0) = −4
4. Résolvez là maintenant de manière approchée avec odeint, sur l’intervalle [0, 30], avec 100 points d’interpolation.
Expliquez le résultat.
7. ordinary differential equations integration !
8. C’est pour cette raison que je vous l’ai fait rédiger comme ça, en fait.
Svartz
Page 4/4
2014/2015