#! /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. def q0(): r""" Préliminaires : - Comme pour tous ces exercices, on commence par **importer** les modules requis, ici `numpy `_, `numpy.linalg `_ et `matplotlib `_ : >>> 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. 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 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]]) 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] 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. def q2(): r""" Question 2. On passe à des simulations probabilistes. """ pass # %% Question 2.a. 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 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 `_. .. 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 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. 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. 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. 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. 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) `_. - 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. 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. 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