Code source de PSI_Mat2_2015_25

#! /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 2, http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/PSI-Mat2-2015-25.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 numpy.linalg as LA  # Convention recommandée
from numpy.polynomial import Polynomial as P  # Convention recommandée

#: On définit le monôme X, on pourra ensuite travailler avec plus facilement.
#: Il faut utiliser la classe :py:class:`numpy.polynomial.Polynomial` (importée comme ``P`` avec ``from numpy.polynomial import Polynomial as P``) :
#:
#: >>> X = P([0, 1])
#: >>> print(X)  # Affiche la liste des coefficients
#: poly([0, 1])
#:
#: Pour s'en servir, par exemple définir :math:`Q(X) = 1 + 2X + 17X^3`, on écrit les maths comme en Python :
#:
#: >>> Q = 1 + 2*X + 17*X**3
#: >>> print(Q)  # Affiche la liste des coefficients
#: poly([1, 2, 0, 17])
X = P([0, 1])


# %% 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/>`_ et `numpy.linalg <http://docs.scipy.org/doc/numpy/reference/routines.linalg.html>`_ : >>> import numpy as np >>> import numpy.linalg as LA # Convention recommandée """ pass
# %% Question 1.
[docs]def q1(): r""" Question 1. Rien de difficile, on doit écrire une fonction ``A_na(n, a)`` qui renvoie cette matrice :math:`A_{n,a}`. On peut le faire à la main, en partant d'une matrice ``np.zeros(n,n)`` plein de zéros, qu'on remplit avec des boucles ``for``. On peut gagner du temps en utilisant :py:func:`numpy.diag` et son argument optionnel ``k`` qui permet de dire sur quelle diagonale on place le vecteur (``k=0`` pour la diagonale principale, ``k=1`` pour la diagonale supérieure, et ``k=-1`` pour la diagonale inférieure) : >>> a = 5 # Exemple >>> n = 4 # Exemple >>> diagsup = np.diag([1/a]*(n-1), k=1) # n-1 fois 1/a sur la diag sup >>> diaginf = np.diag([a]*(n-1), k=-1) # n-1 fois a sur la diag inf >>> A_na = diagsup + diaginf # Somme des deux matrices (n,n) >>> print(A_na) # Elle a la bonne forme, OK [[ 0. 0.2 0. 0. ] [ 5. 0. 0.2 0. ] [ 0. 5. 0. 0.2] [ 0. 0. 5. 0. ]] """ pass
[docs]def A_na(n, a): r""" Implémente la fonction décrite en `Q1 <#PSI_Mat2_2015_25.q1>`_. """ return np.diag([1 / a] * (n - 1), k=1) + np.diag([a] * (n - 1), k=-1)
# %% Question 2.
[docs]def q2(): r""" Question 2. On utilise :py:func:`numpy.linalg.eigvals` pour calculer de façon approchée les valeurs propres de :math:`A_{n,a}`. >>> from itertools import product >>> for (n, a) in product(range(3, 8 + 1), [-2, -1, 1, 2, 3]): ... print("- Pour n =", n, "et a =", a, ": sp(A_na) =", LA.eigvals(A_na(n, a))) # doctest: +ELLIPSIS - Pour n = 3 et a = -2 : sp(A_na) = [ -1.41421356e+00 9.77950360e-17 1.41421356e+00] - Pour n = 3 et a = -1 : sp(A_na) = [ -1.41421356e+00 9.77950360e-17 1.41421356e+00] - Pour n = 3 et a = 1 : sp(A_na) = [ -1.41421356e+00 9.77950360e-17 1.41421356e+00] - Pour n = 3 et a = 2 : sp(A_na) = [ -1.41421356e+00 9.77950360e-17 1.41421356e+00] - Pour n = 3 et a = 3 : sp(A_na) = [ -1.41421356e+00 -1.66533454e-16 1.41421356e+00] ... On remarque que le spectre de :math:`A_{n,a}` ne semble pas dépendre de la valeur de :math:`a` ! On le démontrera à la question `Q3.c <#PSI_Mat2_2015_25.q3_c>`_. On peut faire sans :py:func:`itertools.product`, avec deux boucles for imbriquées : >>> for n in range(3, 8 + 1): # doctest: +ELLIPSIS ... for a in [-2, -1, 1, 2, 3]: ... print("- Pour n =", n, "et a =", a, ": sp(A_na) =", LA.eigvals(A_na(n, a))) - Pour n = 3 et a = -2 : sp(A_na) = [ -1.41421356e+00 9.77950360e-17 1.41421356e+00] - Pour n = 3 et a = -1 : sp(A_na) = [ -1.41421356e+00 9.77950360e-17 1.41421356e+00] ... """ for n in range(3, 8 + 1): for a in [-2, -1, 1, 2, 3]: print("- Pour n =", n, "et a =", a, ": sp(A_na) =", LA.eigvals(A_na(n, a)))
# %% Question 3.
[docs]def q3(): r""" Question 3. On va devoir utiliser ici le module `numpy.polynomial <https://docs.scipy.org/doc/numpy/reference/routines.polynomials.package.html>`_. On n'a jamais travaillé avec durant l'année, je recommande donc la lecture de `ce petit tutoriel (en anglais) <https://docs.scipy.org/doc/numpy/reference/routines.polynomials.classes.html>`_. `Cette fiche de révision (en français) <http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/Python-polynomes.pdf>`_ peut aider aussi. Je recommande d'importer la classe :py:class:`numpy.polynomial.Polynomial` du module :py:mod:`numpy.polynomial` comme suit : >>> from numpy.polynomial import Polynomial as P # Convention recommandée Dès lors, on peut définir un nouveau polynôme par ``P(list_coef)`` en donnant la liste de ses coefficients. On peut par exemple commencer par le monôme X : >>> X = P([0, 1]) # On définit le monôme X C'est très pratique de travailler avec ce X puisqu'on peut écrire des polynômes comme en maths : >>> Q1 = 1 + 2*X + 17*X**3 >>> Q2 = Q1 - 2*X >>> print(Q2) poly([ 1. 0. 0. 17.]) .. note:: Ce polynôme :math:`Q_2(X) = Q_1(X) - 2X = 1 + 17 X^3` n'est pas "joliment" affiché, mais on nous donne sa liste de coefficients, ``[1, 0, 0, 17]`` qui correspond bien à :math:`1 + 17 X^3` = ``1 + 17 * X ** 3``. On pourrait modifier la méthode :py:meth:`numpy.polynomial.Polynomial.__str__` pour afficher joliment mais ce n'est pas le but de ce sujet d'oral (et ce genre de chose est *clairement hors programme*). Si vous êtes curieux, `ce morceau de code <https://bitbucket.org/snippets/lbesson/j6dpz#file-nicedisplay_numpy_polynomial_Polynomial.py>`_ implémente et démontre cette idée. """ pass
# %% Question 3.a.
[docs]def q3_a(): r""" Question 3.a. On commence par définir manuellement les premiers polynômes :math:`P_{1}` et :math:`P_{2}` : >>> P1 = X >>> P2 = X**2 - 1 puis une fonction qui construit :math:`P_{n+2}` en fonction de :math:`P_{n+1}` et :math:`P_{n}` : >>> def Pnextnext(Pn, Pnext): return X * Pnext - Pn On peut dès lors afficher les coefficients des polynômes :math:`P_3, \dots, P_8` : >>> Pn, Pnext = P1, P2 # n = 1 >>> for n in range(3, 8+1): ... Pn, Pnext = Pnext, Pnextnext(Pn, Pnext) # One more ... print("P{} =".format(n), Pnext) P3 = poly([ 0. -2. 0. 1.]) P4 = poly([ 1. 0. -3. 0. 1.]) P5 = poly([ 0. 3. 0. -4. 0. 1.]) P6 = poly([-1. 0. 6. 0. -5. 0. 1.]) P7 = poly([ 0. -4. 0. 10. 0. -6. 0. 1.]) P8 = poly([ 1. 0. -10. 0. 15. 0. -7. 0. 1.]) """ pass
[docs]def Pnextnext(Pn, Pnext): r""" Calcule et renvoie :math:`P_{n+2} = X P_{n+1} - P_{n}` : :: return X * Pnext - Pn .. note:: Ce module permet de travailler avec des polynômes comme on le fait en maths, avec des opérations arithmétiques écrites naturellement. C'est très pratique ! """ return X * Pnext - Pn
# %% Question 3.b.
[docs]def q3_b(): r""" Question 3.b. Et calculer leurs racines est très facile, avec la méthode ``roots`` (:py:meth:`numpy.polynomial.Polynomial.roots`) : >>> P1 = X; P2 = X**2 - 1 >>> Pn, Pnext = P1, P2 # n = 1 >>> for n in range(3, 8 + 1): # doctest: +ELLIPSIS ... Pn, Pnext = Pnext, Pnextnext(Pn, Pnext) # One more ... print("roots P{} =".format(n), list(Pnext.roots())) roots P3 = [-1.4142..., 0.0, 1.4142...] roots P4 = [-1.6180..., -0.6180..., 0.6180..., 1.6180...] roots P5 = [-1.7320..., -0.9999..., 0.0, 0.9999..., 1.7320...] roots P6 = [-1.8019..., -1.2469..., -0.4450..., 0.4450..., 1.2469..., 1.8019...] roots P7 = [-1.8477..., -1.4142..., -0.7653..., 0.0, 0.7653..., 1.4142..., 1.8477...] roots P8 = [-1.8793..., -1.5320..., -1.0000..., -0.3472..., 0.3472..., 1.0000..., 1.5320..., 1.8793...] On constate que :math:`0` est racine pour les :math:`n` impairs, mais pas pour les :math:`n` pairs, ce qu'on démontre à la question `Q4 <#PSI_Mat2_2015_25.q4>`_. """ pass
# %% Question 3.c.
[docs]def q3_c(): r""" Question 3.c. Le lien naturel entre :math:`P_n` et :math:`A_{n,a}` est bien-sûr que : :math:`\chi_{A_{n,a}}(X) = P_n(X)`. On a observé en `Q2 <#PSI_Mat2_2015_25.q2>`_ que le spectre de :math:`A_{n, a}` ne dépend pas de la valeur numérique de :math:`a \neq 0`. (Par contre, les vecteurs propres de :math:`A_{n, a}` pour chacune de ses valeurs propres dépendront bien-sûr de :math:`a`.) On peut procéder par récurrence, pour montrer que : .. math:: \forall a \neq 0, \forall n \geq 1, \chi_{A_{n, a}}(X) = P_n(X). - Pour :math:`n = 1`, :math:`A_{1, a} = [[0]]` qui a bien :math:`X = P_1(X)` comme polynôme caractéristique (peu importe :math:`a`). - Pour :math:`n = 2`, :math:`A_{2, a} = [[0, 1/a], [a, 0]]` qui a bien :math:`X^2 - (1/a) a = X^2 - 1 = P_2(X)` comme polynôme caractéristique (peu importe :math:`a`). - Soit :math:`a \neq 0`. Maintenant, si c'est vrai pour :math:`A_{n, a}` et :math:`A_{n+1, a}`, montrons que c'est encore vrai pour :math:`A_{n+2, a}`. On peut développer :math:`\chi_{A_{n+2, a}}(X)` selon sa première ligne, pour obtenir (en écriture par bloc pour le bloc de fin :math:`A_{n, a}`) : .. math:: \chi_{A_{n+2, a}}(X) = X \times \chi_{A_{n+1, a}}(X) - (-a) \mathrm{det}([1/a, 0, \dots, 0], [-a , A_{n, a}]]) On développe une seconde fois la matrice de droite, qui annule le terme :math:`(-a)` grâce au :math:`(-1/a)`; et par hypothèse de récurrence on a remplacé :math:`\chi_{A_{n+1, a}}(X)` par :math:`P_{n+1}(X)` et on remplace :math:`\chi_{A_{n, a}}(X)` par :math:`P_{n}(X)` : .. math:: \chi_{A_{n+2, a}}(X) &= X P_{n+1}(X) - (-a) (-1/a) \chi_{A_{n, a}}(X) \\ &= X P_{n+1}(X) - P_{n}(X) = P_{n+2}(X) Ainsi, on a bien montré le résultat annoncé. (Sans même avoir besoin de montrer d'abord l'indépendance en :math:`a`). """ pass
# %% Question 4.
[docs]def q4(): r""" - **Question 4.a)** : :math:`A_{n,a}` est-elle inversible ? D'après le lien entre :math:`P_n` et :math:`A_{n,a}` démontré en `Q3.c <#PSI_Mat2_2015_25.q3_c>`_, on sait que :math:`A_{n,a}` est inversible *ssi* :math:`0 \not\in \mathrm{Sp}(A_{n,a})` *ssi* :math:`0 \not\in \mathrm{Roots}(P_n)` *ssi* :math:`P_n(0) \neq 0`. Or, par construction, on a :math:`P_{n+2}(0) = - P_{n}(0)` donc par récurrence on obtient que :math:`P_n(0) \neq 0` *ssi* :math:`P_1(0) \neq 0` si :math:`n` est impair et :math:`P_2(0) \neq 0` si :math:`n` est pair. Or :math:`P_1, P_2` sont connus : :math:`P_1(0) = 0` et :math:`P_2(0) = -1`. Ainsi, :math:`P_n(0) \neq 0` *ssi* :math:`n` est pair. **Conclusion :** :math:`A_{n,a}` est inversible *ssi* :math:`n` est pair. - **Question 4.b)** : :math:`A_{n,a}` est-elle diagonalisable ? D'après le lien entre :math:`P_n` et :math:`A_{n,a}` démontré en `Q3.c <#PSI_Mat2_2015_25.q3_c>`_, on sait que :math:`A_{n,a}` est diagonalisable *ssi* :math:`P_n` est scindé à racines simples :math:`\mathrm{Roots}(P_n)` est de cardinal exactement :math:`n`. Nos observations en `Q3.b <#PSI_Mat2_2015_25.q3_b>`_ nous laissent suggérer que c'est toujours le cas. À vous de faire ce dernier morceau : montrez que :math:`P_n` est *toujours* scindé à racines simples. .. hint:: Pour montrer qu'un polynôme réel n'a que des racines simples, on peut montrer qu'il n'a aucune (au moins) double. Et on sait qu'une racine :math:`\alpha` est (au moins) double *ssi* elle annule aussi :math:`P_{n}'`. On a :math:`P_{n+2}'(X) = P_{n+1}(X) + P_{n+1}'(X) - P_{n}'(X)` (d'après la définition par récurrence de :math:`P_{n+2}`), donc si :math:`\alpha` annule :math:`P_{n+2}(X)`, on a déjà :math:`\alpha P_{n+1}(\alpha) = P_{n}(\alpha)`, et donc si elle est (au moins) double, on a aussi :math:`P_{n+1}(\alpha) = P_{n}'(\alpha) - P_{n+1}'(\alpha)`. En combinant les deux, on obtient que : :math:`P_{n}(\alpha) = \alpha (P_{n}'(\alpha) - P_{n+1}'(\alpha))`. À vous de terminer ! **En fait, ça a l'air assez dur**. Nous n'avons pas réussi en TD avec toute la classe. Si vous avez une idée, `contactez-moi svp <http://perso.crans.org/besson/contact/>`_. **Conclusion :** :math:`A_{n,a}` est diagonalisable pour tout :math:`n`. """ pass
# %% Question 5.
[docs]def q5(): r""" Question 5. Soit :math:`a \in \mathbb{R}^{*}` et :math:`n \geq 1`. On va d'abord utiliser la **borne de Cauchy** pour la localisation des racines du polynôme :math:`P_n(X)` (qui sont exactement les valeurs propres de :math:`A_{n, a}`, d'après `Q3.c <#PSI_Mat2_2015_25.q3_c>`_). .. note:: Plus de détails sur ce résultat se trouvent par exemple `sur cette page Wikipédia <https://fr.wikipedia.org/wiki/Th%C3%A9or%C3%A8me_de_Rouch%C3%A9#Applications>`_ (comme application du théorème de Rouché). La borne de Cauchy, pour un polynôme unitaire (entier, réel, ou complexe) :math:`Q_n(X) = X^n + \sum_{i=0}^{n-1} a_{i} X^i`, énonce que toutes ses racines sont contenues dans le disque centré en zéro de rayon :math:`1 + R_n`, défini par : .. math:: R_n := \sum_{i=0}^{n-1} |a_{i}|. Ici, :math:`P_n` est bien-sûr unitaire, et on doit donc borner cette somme :math:`\sum_{i=0}^{n-1} |a_{i}^{(n)}|` (où l'on note :math:`a_{i}^{(n)}` le ième coefficient du polynôme :math:`P_n` de degré :math:`n`). Par récurrence, on a :math:`a_{i}^{(n+2)} = a_{i-1}^{(n+1)} - a_{i}^{(n)}`, donc :math:`|a_{i}^{(n+2)}| \leq |a_{i-1}^{(n+1)}| + |a_{i}^{(n)}|` par une inégalité triangulaire. Et comme :math:`a_{n+1}^{(n)} = 0`, et en posant :math:`a_{-1}^{(n+1)} := 0` (pour simplifier les notations), on peut directement obtenir que : .. math:: R_{n+2} &:= \sum_{i=0}^{n+2-1} |a_{i}^{(n+2)}| \\ &= |a_{n}^{(n+1)}| + \sum_{i=0}^{n} |a_{i}^{(n+2)}| \\ &\leq |a_{n}^{(n+1)}| + \sum_{i=0}^{n} |a_{i-1}^{(n+1)}| + |a_{i}^{(n)}| &= \sum_{i=0}^{n+1} |a_{i}^{(n+1)}| + \sum_{i=0}^{n} |a_{i}^{(n)}| \\ &= R_{n+1} + R_{n} \;\;\; (1) On a :math:`R_1 = 1` et :math:`R_2 = 2`, donc on peut montrer par récurrence immédiate que :math:`R_n \leq 2^n` d'après l'équation précédente :math:`(1)`. On doit encore montrer que les racines de :math:`P_n(X)` sont réelles, mais on le sait déjà puisqu'il est scindé (sur :math:`\mathbb{R}`) à racines simples réelles (d'après `Q4 <#PSI_Mat2_2015_25.q4>`_). **Conclusion :** pour tout :math:`a \neq 0`, et :math:`n \in \mathbb{N}`, les valeurs propres de :math:`A_{n, a}` sont contenues dans l'intervalle :math:`[-1 - 2^n, 1 + 2^n]`. .. note:: C'est un résultat très large : on a observé en `Q3.b <#PSI_Mat2_2015_25.q3_b>`_ que les racines restent petites. Par exemple, pour :math:`n = 80`, :math:`\max_{i} |\lambda_i(A_{n, a})|` valait :math:`\simeq 2.368`, beaucoup (beaucoup) plus petit que notre borne :math:`2^{80} \simeq 10^{24}` ! On va donc essayer d'obtenir une meilleure borne, et elle sera bien meilleure ! On rappelle que la norme infinie d'un vecteur :math:`\mathbf{x}` ou d'une matrice :math:`\mathbf{A}` est définie par : .. math:: \| \mathbf{x} \|_{\infty} &:= \max_{1 \leq i \leq n} |x_i|, \;\;\text{pour}\; \mathbf{x} = (x_i)_{1 \leq i \leq n} \in \mathbb{R}^n \\ \| \mathbf{A} \|_{\infty} &:= \max_{1 \leq i,j \leq n} |A_{i,j}|, \;\;\text{pour}\; \mathbf{A} = (A_{i,j})_{1 \leq i,j \leq n} \in \mathcal{M}_{n,n}(\mathbb{R}). Soit :math:`n \geq 1`, et :math:`a := 1`. Prenons :math:`\lambda` une valeur propre de :math:`A := A_{n, a}`, et :math:`X` un vecteur propre (donc non nul) associé à :math:`A` et :math:`\lambda`. On a donc :math:`A X = \lambda X`. Dès lors, :math:`\| A X \|_{\infty} = |\lambda| \| X \|_{\infty}`, donc si on prend un :math:`X` de norme infinie :math:`= 1`, on obtient que :math:`|\lambda| = \| A X \|_{\infty}`. Mais :math:`A` a une forme bien spéciale : elle n'a qu'un ou deux coefficients par ligne, qui sont tous égaux à :math:`a = 1`. Donc :math:`\| A X \|_{\infty} \leq 2 \| X \|_{\infty} = 2` par une simple inégalité triangulaire. Ainsi :math:`|\lambda| \leq 2`. Et voilà. BOUM ! Donc on vient de montrer que le segment :math:`[-2, 2]` contient le spectre de :math:`A = A_{n, a}`. .. note:: Pendant l'oral, ce n'est pas très grave de ne pas obtenir "le meilleur" résultat possible. Si l'examinateur attendait quelque chose de plus précis, il/elle vous le dira, et vous serez capable de travailler ensemble pour obtenir le résultat attendu. """ 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_25.py