#! /usr/bin/env python3
# -*- coding: utf-8; mode: python -*-
r""" Correction d'un exercice de maths avec Python, issu des annales pour les oraux du concours Centrale-Supélec.
- *Sujet :* PSI 3, http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/PSI-Mat2-2015-26.pdf
- *Dates :* séances de révisions, vendredi 10 et samedi 11 juin (2016).
- *Auteur :* Lilian Besson, pour le cours d'informatique pour tous en prépa MP au Lycée Lakanal (http://perso.crans.org/besson/infoMP/),
- *Licence :* MIT Licence (http://lbesson.mit-license.org).
-----------------------------------------------------------------------------
"""
from __future__ import print_function, division # Python 2 compatibility
import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as LA # Convention recommandée
import random
# %% Question 0.
[docs]def q0():
r""" Préliminaires :
- Comme pour tous ces exercices, on commence par **importer** les modules requis, ici `numpy <https://docs.scipy.org/doc/numpy/>`_, `numpy.linalg <http://docs.scipy.org/doc/numpy/reference/routines.linalg.html>`_ et `matplotlib <http://matplotlib.org/users/beginner.html>`_ :
>>> import numpy as np
>>> import numpy.linalg as LA # Convention recommandée
>>> import matplotlib.pyplot as plt
- On aura aussi besoin de faire des tirages aléatoires, on peut utiliser au choix soit :py:mod:`numpy.random` soit :py:mod:`random`. Je préfère le plus simple, :py:mod:`random` :
>>> import random
"""
pass
# %% Question 1.
[docs]def q1():
r""" Question 1.
- On écrit d'abord une fonction qui créé cette matrice :math:`M(a, b)`,
- Puis une fonction qui calcule ses deux valeurs propres complexes :math:`\lambda_1, \lambda_2` (de façon approchée),
- Et enfin, on renvoit `ecart(a, b) = abs(lambda1 - lambda2)`.
"""
pass
[docs]def M(a, b):
r""" Renvoie la matrice :math:`M(a, b)` définie dans l'énoncé :
.. math::
M(a, b) := [[3a - 2b, - 6a + 6b + 3], [a - b, - 2a + 3b + 1]]
= \begin{bmatrix}
3a - 2b & - 6a + 6b + 3 \\
a - b & - 2a + 3b + 1
\end{bmatrix}
Exemple :
>>> print(M(1, 3))
[[-3 15]
[-2 8]]
"""
return np.array([[3 * a - 2 * b, - 6 * a + 6 * b + 3], [a - b, - 2 * a + 3 * b + 1]])
[docs]def lam12(a, b):
r""" Calcule et renvoie :math:`\lambda_1, \lambda_2 = \mathrm{Sp}(M(a, b))`.
Exemple :
>>> lam12(1, 3)
(2.0000000000000009, 2.9999999999999991)
>>> [round(l) for l in lam12(1, 3)]
[2.0, 3.0]
"""
Mab = M(a, b)
lambdas = LA.eigvals(Mab)
return lambdas[0], lambdas[1]
[docs]def ecart(a, b):
r""" Calcule et renvoie :math:`|\lambda_1 - \lambda_2|`, approchée à :math:`10^{-2}` près.
Exemple :
>>> ecart(1, 3)
1.0
.. note::
:math:`M(a, b)` est diagonalisable *ssi* ses deux valeurs propres sont distinctes.
Donc tester si ``ecart(a, b)`` est nul est un bon indicateur de la diagonalisabilité de :math:`M(a, b)` (inexact mais presque correct).
Cette idée relie les questions `Q2.a <#PSI_Mat2_2015_25.q2_a>`_, `Q2.c <#PSI_Mat2_2015_25.q2_c>`_, `Q3.c <#PSI_Mat2_2015_25.q3_c>`_.
"""
lambda_1, lambda_2 = lam12(a, b)
eps = abs(lambda_1 - lambda_2)
return round(eps, 3)
# %% Question 2.
[docs]def q2():
r""" Question 2.
On passe à des simulations probabilistes.
"""
pass
# %% Question 2.a.
[docs]def q2_a():
r""" Question 2.a.
- On va d'abord écrire une fonction :py:func:`geom` qui effectue un tirage selon une loi géométrique (de paramètre ``p``),
- Ensuite, la fonction :py:func:`hasard` est immédiate à écrire.
.. note:: ``hasard(p) / 500`` est une valeur approchée de la probabilité que ``ecart(A, B)`` soit non-nul, soit que :math:`M(A, B)` soit diagonalisable.
"""
pass
[docs]def geom(p):
r""" Effectue un tirage selon une loi géométrique de paramètre ``p``.
Par définition, en utilisant :py:func:`random.random`, on effectue ce tirage avec ce code : ::
n = 0
# Si 0 < p < 1, on est sûr que cette boucle termine (mais ça peut prendre du temps !)
while random.random() > p:
# L'expérience de probabilité p a échoué, on compte +1
n += 1
# On renvoie le nombre d'échec avant le succès
return n
Pour plus de détails, `cf. cette page Wikipédia si besoin <https://fr.wikipedia.org/wiki/Loi_g%C3%A9om%C3%A9trique>`_.
.. note:: Cette fonction est au programme, il faut la connaître !
.. seealso::
On aurait aussi pu utiliser :py:func:`numpy.random.geometric`, qui fonctionne directement :
>>> np.random.seed(0) # On enlève l'aléatoire, pour que l'exemple soit le même à chaque fois
>>> p = 0.3 # Valeur arbitraire de p = 0.3
>>> a = np.random.geometric(p)
>>> b = np.random.geometric(p)
>>> ecart(a, b)
0.0
"""
assert 0 < p < 1, "Erreur, non 0 < p < 1"
n = 0
# Si 0 < p < 1, on est sûr que cette boucle termine (mais ça peut prendre du temps !)
while random.random() > p:
# L'expérience de probabilité p a échoué, on compte +1
n += 1
# On renvoie le nombre d'échec avant le succès
return n
[docs]def hasard(p, N=500):
r""" Fonction assez élémentaire : pour ``N = 500`` fois, on tire ``a`` et ``b`` grâce à :py:func:`geom`, et on compte +1 si :py:func:`ecart` a donné ``ecart(a, b) >= 1e-1``.
.. seealso::
En utilisant :py:func:`numpy.random.geometric`, on peut générer ``N = 500`` tirages en une seule ligne :
>>> np.random.seed(0) # On enlève l'aléatoire, pour que l'exemple soit le même à chaque fois
>>> p = 0.3 # Valeur arbitraire de p = 0.3
>>> N = 500
>>> a = np.random.geometric(p, N)
>>> b = np.random.geometric(p, N)
>>> # Cette sum( ecart(a[i], b[i]) >= 1e-1 ...) compte le nombre de ecart(a, b) qui sont >= 1e-1
>>> nb = sum(ecart(a[i], b[i]) >= 1e-1 for i in range(N))
>>> nb # 434 sur 500
434
"""
nb = 0
epsilon = 1e-1 # 10e-1
for _ in range(N):
# On tire les deux valeurs a, b des v.a. A, B
a = geom(p)
b = geom(p)
# On calcule l'écart et on vérifie
e = ecart(a, b)
if e >= epsilon:
nb += 1 # Un de plus est >= epsilon !
return nb # Nb empirique de ecart(a, b) >= epsilon
# %% Question 2.b.
[docs]def q2_b(N=500):
r""" Question 2.b.
- On va d'abord créer la liste des valeurs de ``p`` :math:`= \frac{1}{100}, \dots, \frac{99}{100}`,
- Puis pour chaque ``p``, on calcule ``hasard(p) / 500`` grâce à :py:func:`hasard`, qu'on stocke dans une liste de valeurs,
- Enfin, on ouvre une fenêtre :py:mod:`matplotlib` et on affiche (avec :py:func:`matplotlib.pyplot.plot`) ces valeurs.
.. attention::
La fonction :py:func:`hasard` est, comme nom l'indique, aléatoire !
Donc chaque exécution de cette fonction sera différente, et donc chaque graphique sera différent !
Cela devrait donner une figure comme ça :
.. image:: PSI_Mat2_2015_26__figure1.png
"""
P = 100 # Nb-1 de valeurs de p à essayer
valeurs_p = [p / P for p in range(1, P)]
valeurs_hasard = [hasard(p, N=N) / N for p in valeurs_p]
plt.figure()
plt.grid(True)
plt.plot(valeurs_p, valeurs_hasard, 'b+-', label=r"$\frac{\mathrm{hasard}(p)}{%d}$" % N)
plt.title(r"$(p, \frac{\mathrm{hasard}(p)}{%d})$ pour $p = \frac{1}{100}, ..., \frac{99}{100}$" % N)
plt.savefig("PSI_Mat2_2015_26__figure1.png", dpi=180)
plt.show()
# %% Question 2.c.
[docs]def q2_c(N=500):
r""" Question 2.c.
On a juste à ajouter trois lignes dans la fonction précédente `Q2.b <#PSI_Mat2_2015_25.q2_b>`_, en calculant les valeurs de :math:`\frac{2 - 2 p + p^2}{2 - p}` pour ``p = np.linspace(1e-2, 1 - 1e-2, 400)`` (par exemple, pour ``400`` valeurs).
On est très content d'observer que cette formule colle très bien aux valeurs numériques simulées en `Q2.b <#PSI_Mat2_2015_25.q2_b>`_.
(et donc, ça donne une très bonne indication pour la dernière question `Q3.c <#PSI_Mat2_2015_25.q3_c>`_).
Cela devrait donner une figure comme ça :
.. image:: PSI_Mat2_2015_26__figure2.png
.. note:: J'ai fait ``N = 5000`` simulations plutôt que ``N = 500`` pour mieux voir que la formule colle aux valeurs aléatoires.
"""
P = 100 # Nb-1 de valeurs de p à essayer
valeurs_p = [p / P for p in range(1, P)]
valeurs_hasard = [hasard(p, N=N) / N for p in valeurs_p]
T = 500 # Nb de points à calculer pour la courbe "continue"
valeurs_t = np.linspace(1e-2, 1 - 1e-2, T)
valeurs_courbe = (2 - 2 * valeurs_t + valeurs_t**2) / (2 - valeurs_t)
plt.figure()
plt.hold(True)
plt.grid(True)
plt.plot(valeurs_p, valeurs_hasard, 'b+-', label=r"$\frac{\mathrm{hasard}(p)}{%d}$" % N)
plt.plot(valeurs_t, valeurs_courbe, 'r:', label=r"$\frac{2 - 2 p + p^2}{2 - p}$")
plt.title(r"$(p, \frac{\mathrm{hasard}(p)}{%d})$ pour $p = \frac{1}{100}, ..., \frac{99}{100}$ et courbe $\frac{2 - 2 p + p^2}{2 - p}$" % N)
plt.legend() # Utilise les label donnés dans les plt.plot(..., label="...")
plt.savefig("PSI_Mat2_2015_26__figure2.png", dpi=180)
plt.show()
# %% Question 3.
[docs]def q3():
r""" Question 3.
On définit la matrice :math:`M'(a, b)` suivante :
.. math::
M'(a, b) := [[a+1, 1], [0, b]]
= \begin{bmatrix}
a + 1 & 1 \\
0 & b
\end{bmatrix}
"""
pass
# %% Question 3.a.
[docs]def q3_a():
r""" Question 3.a.
C'est une pure question de maths, à faire sur papier.
Pour des matrices :math:`2 \times 2`, si elles sont diagonalisables,
il suffit de vérifier qu'elles ont même polynôme caractéristique (et qu'il soit scindé à racines simples), soit mêmes traces et mêmes déterminant :
- Trace :
.. math::
\mathrm{tr}(M(a, b)) &= (3a - 2b) + (-2a + 3b + 1) = a + b + 1. \\
\mathrm{tr}(M'(a, b)) &= (a + 1) + (b) = a + b + 1.
- Déterminant :
.. math::
\mathrm{det}(M(a, b)) &= (3a - 2b) \times (-2a + 3b + 1) - (a - b) \times (-6a + 6b + 3) = b (a + 1). \\
\mathrm{det}(M'(a, b)) &= (a + 1) \times b - 0 \times 1 = b (a + 1).
.. note:: Vous pouvez vérifier le calcul de :math:`\mathrm{det}(M(a, b))` vous-même, ou par exemple `sur ce site (SymPyGamma.com) <http://www.sympygamma.com/input/?i=simplify%28%283+*+a+-+2+*+b%29+*+%28-2+*+a+%2B+3+*+b+%2B+1%29+-+%28a+-+b%29+*+%28-6+*+a+%2B+6+*+b+%2B+3%29%29>`_.
- Si elles ne sont pas diagonalisables, en fait ça implique que :math:`a + 1 = b`, et donc en regardant la forme simplifiée de :math:`M(a, b)` dans ce cas, on arrive assez vite à montrer qu'elle est semblable à :math:`M'(a, b)`.
"""
pass
# %% Question 3.b.
[docs]def q3_b():
r""" Question 3.b.
Une matrice réelle :math:`2 \times 2` est diagonalisable *ssi* elle admet deux valeurs propres distinctes.
Donc la matrice :math:`M'(a, b)` est diagonalisable *ssi* :math:`a + 1 \neq b`.
"""
pass
# %% Question 3.c.
[docs]def q3_c():
r""" Question 3.c.
D'après `Q3.a <#PSI_Mat2_2015_25.q3_a>`_ et `Q3.b <#PSI_Mat2_2015_25.q3_b>`_,
:math:`M(a, b)` est diagonalisable
*ssi* :math:`M'(a, b)` l'est (elles sont semblables)
*ssi* :math:`a + 1 \neq b`.
Calculons donc la probabilité que :math:`A + 1 \neq B`
(on rappelle que les deux variables aléatoires :math:`A` et :math:`B` suivent :math:`\mathcal{G}(p)` une loi géométrique de paramètre :math:`p`).
.. math::
\mathbb{P}(A + 1 \neq B) &= \sum_{k=0}^{+\infty} \mathbb{P}(A + 1 \neq k, B = k) \\
&= \sum_{k=0}^{+\infty} \mathbb{P}(A \neq k - 1) \mathbb{P}(B = k) \;\;\text{car}\;A,B \;\text{indépendants}\\
&= \mathbb{P}(B = 0) + \sum_{k=1}^{+\infty} (1 - \mathbb{P}(A = k - 1)) ((1-p)^{k} p) \;\;\text{car}\;B \thicksim \mathcal{G}(p) \\
&= p + p \sum_{k=1}^{+\infty} (1 - ((1-p)^{k-1} p)) (1-p)^{k} \;\;\text{car}\;A \thicksim \mathcal{G}(p) \\
&= p + p \sum_{k=1}^{+\infty} (1-p)^{k} - p^2 \sum_{k=0}^{+\infty} (1-p)^{2k-1} \;\;\text{et}\;0 < (1-p) < 1 \\
&= p + p \frac{(1-p)}{1 - (1-p)} - p^2 \frac{1}{(1-p)} \frac{(1-p)^2}{1 - (1-p)^2} \\
&= p + (1-p) - \frac{p^2}{1-p} \frac{(1-p)^2}{2p-p^2} \\
&= 1 - \frac{p^2(1-p)^2}{p(2-p)(1-p)} \\
&= 1 - \frac{p(1-p)}{(2-p)} \\
&= \frac{(2-p) - p(1-p)}{2-p} \\
&= \frac{2-2p+p^2}{2-p}
On retrouve bien la formule donnée par l'énoncé en `Q3.c <#PSI_Mat2_2015_25.q3_c>`_. OUF !
"""
pass
# %% Fin de la correction.
if __name__ == '__main__':
from doctest import testmod
print("Test automatique de toutes les doctests ecrites dans la documentation (docstring) de chaque fonction :")
testmod(verbose=True)
# testmod()
print("\nPlus de details sur ces doctests peut etre trouve dans la documentation de Python:\nhttps://docs.python.org/3/library/doctest.html (en anglais)")
# Fin de PSI_Mat2_2015_26.py