PSI_Mat2_2015_25

  • Thème : algèbre linéaire, construction de matrices, calculs de valeurs propres, polynômes et diagonalisations.

Documentation

Correction d’un exercice de maths avec Python, issu des annales pour les oraux du concours Centrale-Supélec.

PSI_Mat2_2015_25.X = Polynomial([ 0., 1.], [-1, 1], [-1, 1])

On définit le monôme X, on pourra ensuite travailler avec plus facilement. Il faut utiliser la classe 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 \(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])
PSI_Mat2_2015_25.q0()[source]

Préliminaires :

  • Comme pour tous ces exercices, on commence par importer les modules requis, ici numpy et numpy.linalg :
>>> import numpy as np
>>> import numpy.linalg as LA  # Convention recommandée
PSI_Mat2_2015_25.q1()[source]

Question 1.

Rien de difficile, on doit écrire une fonction A_na(n, a) qui renvoie cette matrice \(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 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. ]]
PSI_Mat2_2015_25.A_na(n, a)[source]

Implémente la fonction décrite en Q1.

PSI_Mat2_2015_25.q2()[source]

Question 2.

On utilise numpy.linalg.eigvals() pour calculer de façon approchée les valeurs propres de \(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)))  
- 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 \(A_{n,a}\) ne semble pas dépendre de la valeur de \(a\) ! On le démontrera à la question Q3.c.

On peut faire sans itertools.product(), avec deux boucles for imbriquées :

>>> 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)))
- 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]
...
PSI_Mat2_2015_25.q3()[source]

Question 3.

On va devoir utiliser ici le module numpy.polynomial.

On n’a jamais travaillé avec durant l’année, je recommande donc la lecture de ce petit tutoriel (en anglais).

Cette fiche de révision (en français) peut aider aussi.

Je recommande d’importer la classe numpy.polynomial.Polynomial du module 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 \(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 à \(1 + 17 X^3\) = 1 + 17 * X ** 3.

On pourrait modifier la méthode 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 implémente et démontre cette idée.

PSI_Mat2_2015_25.q3_a()[source]

Question 3.a.

On commence par définir manuellement les premiers polynômes \(P_{1}\) et \(P_{2}\) :

>>> P1 = X
>>> P2 = X**2 - 1

puis une fonction qui construit \(P_{n+2}\) en fonction de \(P_{n+1}\) et \(P_{n}\) :

>>> def Pnextnext(Pn, Pnext): return X * Pnext - Pn

On peut dès lors afficher les coefficients des polynômes \(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.])
PSI_Mat2_2015_25.Pnextnext(Pn, Pnext)[source]

Calcule et renvoie \(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 !

PSI_Mat2_2015_25.q3_b()[source]

Question 3.b.

Et calculer leurs racines est très facile, avec la méthode roots (numpy.polynomial.Polynomial.roots()) :

>>> P1 = X; P2 = X**2 - 1
>>> Pn, Pnext = P1, P2  # n = 1
>>> for n in range(3, 8 + 1):  
...     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 \(0\) est racine pour les \(n\) impairs, mais pas pour les \(n\) pairs, ce qu’on démontre à la question Q4.

PSI_Mat2_2015_25.q3_c()[source]

Question 3.c.

Le lien naturel entre \(P_n\) et \(A_{n,a}\) est bien-sûr que : \(\chi_{A_{n,a}}(X) = P_n(X)\).

On a observé en Q2 que le spectre de \(A_{n, a}\) ne dépend pas de la valeur numérique de \(a \neq 0\). (Par contre, les vecteurs propres de \(A_{n, a}\) pour chacune de ses valeurs propres dépendront bien-sûr de \(a\).)

On peut procéder par récurrence, pour montrer que :

\[\forall a \neq 0, \forall n \geq 1, \chi_{A_{n, a}}(X) = P_n(X).\]

  • Pour \(n = 1\), \(A_{1, a} = [[0]]\) qui a bien \(X = P_1(X)\) comme polynôme caractéristique (peu importe \(a\)).

  • Pour \(n = 2\), \(A_{2, a} = [[0, 1/a], [a, 0]]\) qui a bien \(X^2 - (1/a) a = X^2 - 1 = P_2(X)\) comme polynôme caractéristique (peu importe \(a\)).

  • Soit \(a \neq 0\). Maintenant, si c’est vrai pour \(A_{n, a}\) et \(A_{n+1, a}\), montrons que c’est encore vrai pour \(A_{n+2, a}\). On peut développer \(\chi_{A_{n+2, a}}(X)\) selon sa première ligne, pour obtenir (en écriture par bloc pour le bloc de fin \(A_{n, a}\)) :

    \[\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 \((-a)\) grâce au \((-1/a)\); et par hypothèse de récurrence on a remplacé \(\chi_{A_{n+1, a}}(X)\) par \(P_{n+1}(X)\) et on remplace \(\chi_{A_{n, a}}(X)\) par \(P_{n}(X)\) :

    \[\begin{split}\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)\end{split}\]

Ainsi, on a bien montré le résultat annoncé. (Sans même avoir besoin de montrer d’abord l’indépendance en \(a\)).

PSI_Mat2_2015_25.q4()[source]
  • Question 4.a) : \(A_{n,a}\) est-elle inversible ?

D’après le lien entre \(P_n\) et \(A_{n,a}\) démontré en Q3.c, on sait que \(A_{n,a}\) est inversible ssi \(0 \not\in \mathrm{Sp}(A_{n,a})\) ssi \(0 \not\in \mathrm{Roots}(P_n)\) ssi \(P_n(0) \neq 0\).

Or, par construction, on a \(P_{n+2}(0) = - P_{n}(0)\) donc par récurrence on obtient que \(P_n(0) \neq 0\) ssi \(P_1(0) \neq 0\) si \(n\) est impair et \(P_2(0) \neq 0\) si \(n\) est pair.

Or \(P_1, P_2\) sont connus : \(P_1(0) = 0\) et \(P_2(0) = -1\).

Ainsi, \(P_n(0) \neq 0\) ssi \(n\) est pair.

Conclusion : \(A_{n,a}\) est inversible ssi \(n\) est pair.

  • Question 4.b) : \(A_{n,a}\) est-elle diagonalisable ?

D’après le lien entre \(P_n\) et \(A_{n,a}\) démontré en Q3.c, on sait que \(A_{n,a}\) est diagonalisable ssi \(P_n\) est scindé à racines simples \(\mathrm{Roots}(P_n)\) est de cardinal exactement \(n\).

Nos observations en Q3.b nous laissent suggérer que c’est toujours le cas.

À vous de faire ce dernier morceau : montrez que \(P_n\) est toujours scindé à racines simples.

Indice

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 \(\alpha\) est (au moins) double ssi elle annule aussi \(P_{n}'\).

On a \(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 \(P_{n+2}\)), donc si \(\alpha\) annule \(P_{n+2}(X)\), on a déjà \(\alpha P_{n+1}(\alpha) = P_{n}(\alpha)\), et donc si elle est (au moins) double, on a aussi \(P_{n+1}(\alpha) = P_{n}'(\alpha) - P_{n+1}'(\alpha)\). En combinant les deux, on obtient que : \(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.

Conclusion : \(A_{n,a}\) est diagonalisable pour tout \(n\).

PSI_Mat2_2015_25.q5()[source]

Question 5.

Soit \(a \in \mathbb{R}^{*}\) et \(n \geq 1\). On va d’abord utiliser la borne de Cauchy pour la localisation des racines du polynôme \(P_n(X)\) (qui sont exactement les valeurs propres de \(A_{n, a}\), d’après Q3.c).

Note

Plus de détails sur ce résultat se trouvent par exemple sur cette page Wikipédia (comme application du théorème de Rouché).

La borne de Cauchy, pour un polynôme unitaire (entier, réel, ou complexe) \(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 \(1 + R_n\), défini par :

\[R_n := \sum_{i=0}^{n-1} |a_{i}|.\]

Ici, \(P_n\) est bien-sûr unitaire, et on doit donc borner cette somme \(\sum_{i=0}^{n-1} |a_{i}^{(n)}|\) (où l’on note \(a_{i}^{(n)}\) le ième coefficient du polynôme \(P_n\) de degré \(n\)). Par récurrence, on a \(a_{i}^{(n+2)} = a_{i-1}^{(n+1)} - a_{i}^{(n)}\), donc \(|a_{i}^{(n+2)}| \leq |a_{i-1}^{(n+1)}| + |a_{i}^{(n)}|\) par une inégalité triangulaire.

Et comme \(a_{n+1}^{(n)} = 0\), et en posant \(a_{-1}^{(n+1)} := 0\) (pour simplifier les notations), on peut directement obtenir que :

\[\begin{split}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)\end{split}\]

On a \(R_1 = 1\) et \(R_2 = 2\), donc on peut montrer par récurrence immédiate que \(R_n \leq 2^n\) d’après l’équation précédente \((1)\).

On doit encore montrer que les racines de \(P_n(X)\) sont réelles, mais on le sait déjà puisqu’il est scindé (sur \(\mathbb{R}\)) à racines simples réelles (d’après Q4).

Conclusion : pour tout \(a \neq 0\), et \(n \in \mathbb{N}\), les valeurs propres de \(A_{n, a}\) sont contenues dans l’intervalle \([-1 - 2^n, 1 + 2^n]\).

Note

C’est un résultat très large : on a observé en Q3.b que les racines restent petites. Par exemple, pour \(n = 80\), \(\max_{i} |\lambda_i(A_{n, a})|\) valait \(\simeq 2.368\), beaucoup (beaucoup) plus petit que notre borne \(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 \(\mathbf{x}\) ou d’une matrice \(\mathbf{A}\) est définie par :

\[\begin{split}\| \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}).\end{split}\]

Soit \(n \geq 1\), et \(a := 1\). Prenons \(\lambda\) une valeur propre de \(A := A_{n, a}\), et \(X\) un vecteur propre (donc non nul) associé à \(A\) et \(\lambda\).

On a donc \(A X = \lambda X\). Dès lors, \(\| A X \|_{\infty} = |\lambda| \| X \|_{\infty}\), donc si on prend un \(X\) de norme infinie \(= 1\), on obtient que \(|\lambda| = \| A X \|_{\infty}\). Mais \(A\) a une forme bien spéciale : elle n’a qu’un ou deux coefficients par ligne, qui sont tous égaux à \(a = 1\). Donc \(\| A X \|_{\infty} \leq 2 \| X \|_{\infty} = 2\) par une simple inégalité triangulaire. Ainsi \(|\lambda| \leq 2\). Et voilà. BOUM !

Donc on vient de montrer que le segment \([-2, 2]\) contient le spectre de \(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.

Sortie du script

$ python PSI_Mat2_2015_25.py
Test automatique de toutes les doctests ecrites dans la documentation (docstring) de chaque fonction :
Trying:
    import numpy as np
Expecting nothing
ok
Trying:
    import numpy.linalg as LA  # Convention recommandée
Expecting nothing
ok
Trying:
    a = 5   # Exemple
Expecting nothing
ok
Trying:
    n = 4  # Exemple
Expecting nothing
ok
Trying:
    diagsup = np.diag([1/a]*(n-1), k=1)  # n-1 fois 1/a sur la diag sup
Expecting nothing
ok
Trying:
    diaginf = np.diag([a]*(n-1), k=-1)   # n-1 fois a sur la diag inf
Expecting nothing
ok
Trying:
    A_na = diagsup + diaginf  # Somme des deux matrices (n,n)
Expecting nothing
ok
Trying:
    print(A_na)  # Elle a la bonne forme, OK
Expecting:
    [[ 0.   0.2  0.   0. ]
     [ 5.   0.   0.2  0. ]
     [ 0.   5.   0.   0.2]
     [ 0.   0.   5.   0. ]]
ok
Trying:
    from itertools import product
Expecting nothing
ok
Trying:
    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
Expecting:
    - 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]
    ...
ok
Trying:
    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)))
Expecting:
    - 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]
    ...
ok
Trying:
    from numpy.polynomial import Polynomial as P  # Convention recommandée
Expecting nothing
ok
Trying:
    X = P([0, 1])  # On définit le monôme X
Expecting nothing
ok
Trying:
    Q1 = 1 + 2*X + 17*X**3
Expecting nothing
ok
Trying:
    Q2 = Q1 - 2*X
Expecting nothing
ok
Trying:
    print(Q2)
Expecting:
    poly([  1.   0.   0.  17.])
ok
Trying:
    P1 = X
Expecting nothing
ok
Trying:
    P2 = X**2 - 1
Expecting nothing
ok
Trying:
    def Pnextnext(Pn, Pnext): return X * Pnext - Pn
Expecting nothing
ok
Trying:
    Pn, Pnext = P1, P2  # n = 1
Expecting nothing
ok
Trying:
    for n in range(3, 8+1):
        Pn, Pnext = Pnext, Pnextnext(Pn, Pnext)  # One more
        print("P{} =".format(n), Pnext)
Expecting:
    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.])
ok
Trying:
    P1 = X; P2 = X**2 - 1
Expecting nothing
ok
Trying:
    Pn, Pnext = P1, P2  # n = 1
Expecting nothing
ok
Trying:
    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()))
Expecting:
    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...]
ok
6 items had no tests:
    __main__
    __main__.A_na
    __main__.Pnextnext
    __main__.q3_c
    __main__.q4
    __main__.q5
6 items passed all tests:
   2 tests in __main__.q0
   6 tests in __main__.q1
   3 tests in __main__.q2
   5 tests in __main__.q3
   5 tests in __main__.q3_a
   3 tests in __main__.q3_b
24 tests in 12 items.
24 passed and 0 failed.
Test passed.

Plus de details sur ces doctests peut etre trouve dans la documentation de Python:
https://docs.python.org/3/library/doctest.html (en anglais)

Le fichier Python se trouve ici : PSI_Mat2_2015_25.py.