#! /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