[Calcul scientifique / technique par Python] Solution numérique du problème des valeurs propres de la matrice par multiplication de puissance, algèbre linéaire numérique

introduction

En science et en ingénierie, nous résolvons souvent des problèmes de valeurs propres en mécanique quantique et en systèmes oscillateurs.

Méthode d'alimentation est un problème de valeur propre standard.

A\mathbf{u} = \lambda \mathbf{u}

C'est l'un des algorithmes pour trouver la valeur propre maximale (valeur propre dominante) de. ** ** Lorsque vous ouvrez le livre d'introduction autour de cette zone, vous le verrez souvent en premier.

** C'est un contenu très pédagogique pour apprendre la solution numérique rudimentaire du problème des valeurs propres. Dans cet article, nous examinerons cette méthode. ** **

Cette méthode est très simple.

Procédure de calcul et points

  1. Définissez un vecteur propre estimé initial approprié $ u_0 $.
  2. Multipliez $ u_0 $ par A ($ AAAAAA ... \ mathbf {u_0} = A ^ k \ mathbf {u_0} \ (k-> ∞) $.
  3. Ensuite, le vecteur de direction converge vers un vecteur constant $ \ mathbf {u '} $.
  4. ** Dans de nombreux cas, $ \ mathbf {u '} $ est le vecteur propre correspondant à la valeur propre avec la plus grande valeur absolue! ** **
  5. u'Si est un vecteur propre, la valeur propre correspondante\lambdaEst\mathbf{u^T}A\mathbf{u}/(|u|^2)Peut être donné.

Donc, l'essentiel est de multiplier le ** A pour vérifier la convergence du vecteur direction. ** ** Ensuite, l'origine du nom "power method" est due à la méthode d'évaluation de la puissance de la matrice de coefficients A, comme on peut le voir en 2 ci-dessus.

Si la base mathématique de 4 ci-dessus vous préoccupe, veuillez vous référer à l'addendum de cet article

De plus, la méthode de puissance est une méthode permettant de trouver la valeur propre maximale de la valeur absolue de A. ** Avec un peu d'ingéniosité, il est possible de trouver la valeur propre avec la plus petite valeur absolue, la valeur propre la plus proche du nombre complexe donné $ z $, la deuxième valeur propre, etc. ** Je prévois de publier un bon article.


Contenu

Résoudre un problème de valeurs propres simple qui applique la méthode de multiplication de puissance, Cependant, assurez-vous que la valeur propre avec la valeur absolue maximale et le vecteur propre correspondant sont obtenus.

Comme $ A $ de $ A \ mathbf {u} = \ lambda \ mathbf {u} $

A=
\left[
\begin{matrix}
1 & 2 \\
3 & 4 
\end{matrix}
\right]

Penser à.

La solution exacte pour la plus grande valeur propre absolue est $ \ frac {5+ \ sqrt {33}} {2} = 5.372281323 ... $.


code

En tant que condition de jugement de convergence Valeur absolue du changement de $ \ lambda $ en répétant les étapes $ k $ et $ k + 1 $ \frac{|\lambda_{k+1}-\lambda_k|}{| \lambda_k|} \lt \epsilon=0.0001Répétez le calcul jusqu'à.

"""
Problème de valeur propre de la matrice:Méthode d'alimentation:
"""

import numpy as np

A=np.array([[1,2],[3,4]])

x0 = np.array([1,0]); x1 = np.array([0,1])
u = 1.0*x0+2.0*x1 #Le vecteur propre initial. Approprié.

rel_eps = 0.0001 #Condition de convergence de la valeur propre
#Génération de colonnes Krylov

rel_delta_u=100.0
while rel_delta_u >= rel_eps :  #Boucle principale
    uu = u/np.linalg.norm(u) #Normalisation(Réglez la norme sur 1)
    print("u=",uu)
    
    u = np.dot(A,uu.T)

    eigen_value=np.dot(uu,u)/(np.dot(uu,uu.T))

    print("eigen_value=",eigen_value)
    
    delta_u_vec = uu-u/np.linalg.norm(u)
    abs_delta_u_value= np.linalg.norm(delta_u_vec)
    rel_delta_u=abs_delta_u_value/np.abs(eigen_value) #Changement relatif de la valeur propre par rapport aux étapes répétées
    print("rel_delta_u_vec = ",rel_delta_u)


résultat



u= [ 0.41612395  0.9093079 ]
eigen_value= 5.37244655582
rel_delta_u_vec =  3.29180183204e-05

Vous pouvez voir que cela correspond très bien à la solution exacte 5.372281323 ...


Les références

Les livres suivants ont été utiles dans la rédaction de cet article. [1] est une description simple et facile à comprendre. [2] est un résumé de la façon de résoudre les problèmes de valeurs propres en utilisant numpy et scipy.

[1] Gilbert Strang, ["World Standard MIT Textbook Strang: Linear Algebra Introduction"](https://www.amazon.co.jp/%E4%B8%96%E7%95%8C%E6%A8%99 % E6% BA% 96MIT% E6% 95% 99% E7% A7% 91% E6% 9B% B8-% E3% 82% B9% E3% 83% 88% E3% 83% A9% E3% 83% B3% E3% 82% B0-% E7% B7% 9A% E5% BD% A2% E4% BB% A3% E6% 95% B0% E3% 82% A4% E3% 83% B3% E3% 83% 88% E3 % 83% AD% E3% 83% 80% E3% 82% AF% E3% 82% B7% E3% 83% A7% E3% 83% B3-% E3% 82% AE% E3% 83% AB% E3% 83% 90% E3% 83% BC% E3% 83% 88 / dp / 4764904055 / ref = pd_lpo_sbs_14_t_0? _ Encoding = UTF8 & psc = 1 & refRID = 9817PCQXDR5497M5GPS2), Modern Science, 2015.

[2] [[Calcul scientifique / technique par Python] Résolution d'un problème de valeur propre (généralisé) en utilisant numpy / scipy, en utilisant la bibliothèque] (http://qiita.com/sci_Haru/items/034c6f74d415c1c10d0b)


Addenda

La multiplication de puissance est un problème de valeur propre standard

A\mathbf{u} = \lambda \mathbf{u}

Dans

À partir du vecteur initial $ u_0 $

u_i = Au_{i-1}, (i=1,2,...)

Il s'agit d'une méthode pour trouver la valeur propre / le vecteur propre avec la valeur absolue maximale. (Cette $ u_0, Au_0, A ^ 2u_0, A ^ 3u_0 $, ... est [colonne Kurilov](https://ja.wikipedia.org/wiki/%E3%82%AF%E3%83%AA% Il s'appelle E3% 83% AD% E3% 83% 95% E9% 83% A8% E5% 88% 86% E7% A9% BA% E9% 96% 93))

Maintenant, l'équation des valeurs propres $Au_i = \lambda_i u_i$

En **, il n'y a pas de réduction de la valeur propre **, c'est-à-dire $\lambda_i \neq \lambda_j \ (i \neq j)$

Supposer. Aussi, concernant l'ordre de grandeur des valeurs propres

|\lambda_1| > |\lambda_2|>|\lambda_3|...

Supposer.|\lambda_1|Est la valeur propre maximale de la valeur absolue, et nous voulons la trouver par la méthode de multiplication de puissance.

Considérons maintenant un vecteur approprié $ u_0 $. Supposons que le coefficient $ c_i $ soit fixé de sorte que $ u_0 $ puisse être développé comme suit.

u_0 = c_1 \mathbf{u_1} + c_2 \mathbf{u_2}+c3 \mathbf{u_3}+...

Quand $ A ^ k $, qui devrait être $ A $, est appliqué à ceci,

A^k u_0 = c_1 A^k \mathbf{u_1} + c_2 A^k \mathbf{u_2} + c_3 A^k \mathbf{u_3}+...
= c_1 \lambda_1^k \mathbf{u_1} + c_2 \lambda_2^k \mathbf{u_2} + c_3 \lambda_3^k\mathbf{u_3} +...
=\lambda_1^k (c_1 \mathbf{u_1} + c_2 \lambda_2^k/\lambda_1^k \mathbf{u_2} + c_3 \lambda_3^k/\lambda_1^k\mathbf{u_3}) +...

Ce sera.

|\lambda_1|Est la plus grande valeur propre, donc pour un grand nombre de puissance k\lambda_2^k/\lambda_1^kOu\lambda_3^k/\lambda_1^k ...Devrait converger vers zéroest. En d'autres termes

(k-> ∞ )A^k u_0 => \lambda_1^k (c_1 \mathbf{u_1} )

Vous pouvez vous attendre à ce que ce soit le cas.

En appliquant à plusieurs reprises $ A $ au ** vecteur d'essai initial $ u_0 $ de cette manière, un vecteur parallèle au vecteur propre correspondant à la valeur propre avec la valeur absolue maximale est obtenu. ** **

Puisque le vecteur propre a toujours une indéfinité de multiple constant (multiple numérique complexe arbitraire) (cela signifie que la valeur propre ne change pas même si les deux côtés de l'équation propre sont multipliés par un complexe), celui qui convient à l'objectif est sélectionné. Dans la plupart des cas, un avec une norme de 1 est choisi.

La valeur propre maximale absolue $ \ lambda $ provient de l'équation des valeurs propres pour un $ k $ suffisamment grand.

\lambda = \frac{u_k^T A u_k}{u_k^T u_k} =\frac{u_k^T u_{k+1}}{u_k^T u_k}

Peut être calculé comme. Il s'agit du [quotient de Rayleigh](https://ja.wikipedia.org/wiki/%E3%83%AC%E3%82%A4%E3%83%AA%E3%83%BC%E5%95 Il s'appelle% 86).

Comme décrit ci-dessus, il est possible de calculer la valeur absolue maximale du problème des valeurs propres et le vecteur propre correspondant.

Mise en garde

  1. Lors de l'utilisation de la méthode de puissance, la condition qu'il n'y a pas de duplication dans la valeur unique de $ A $ doit être remplie. ** Si les valeurs propres sont rétractées ou proches, la solution ne convergera pas ou nécessitera trop d'étapes itératives pour converger **. 2.La matrice réelle est une valeur propre complexezSi, son nombre complexe conjuguéz*Est également une valeur unique. Ce temps|z|=|z*|Par conséquent, les valeurs uniques seront les mêmes. S'ils ont la valeur absolue maximale, la situation de 1 ci-dessus se produira et la solution ne convergera pas.

Recommended Posts

[Calcul scientifique / technique par Python] Solution numérique du problème des valeurs propres de la matrice par multiplication de puissance, algèbre linéaire numérique
[Calcul scientifique / technique par Python] Solution numérique d'un problème d'oscillateur harmonique unidimensionnel par vitesse Méthode de Berle
[Calcul scientifique / technique par Python] Liste des matrices qui apparaissent dans Hinpan en algèbre linéaire numérique
[Calcul scientifique / technique par Python] Calcul de somme, calcul numérique
[Calcul scientifique / technique par Python] Résolution du problème de la valeur aux limites des équations différentielles ordinaires au format matriciel, calcul numérique
[Calcul scientifique / technique par Python] Solution numérique de l'équation de Laplace-Poisson bidimensionnelle pour la position électrostatique par la méthode Jacobi, équation aux dérivées partielles elliptiques, problème des valeurs aux limites
[Calcul scientifique / technique par Python] Résolution d'équations linéaires simultanées, calcul numérique, numpy
[Calcul scientifique / technique par Python] Marche aléatoire 2D (problème de marche ivre), calcul numérique
[Calcul scientifique / technique par Python] Solution numérique d'équations différentielles ordinaires du premier ordre, problème de valeur initiale, calcul numérique
[Calcul scientifique / technique par Python] Solution numérique d'équations d'ondes unidimensionnelles et bidimensionnelles par méthode FTCS (méthode explicite), équations aux dérivées partielles bi-courbes
[Calcul scientifique / technique par Python] Calcul de matrice inverse, numpy
[Calcul scientifique / technique par Python] Interpolation de Lagrange, calcul numérique
[Calcul scientifique / technique par Python] Calcul du produit de la matrice par l'opérateur @, python3.5 ou supérieur, numpy
[Calcul scientifique / technique par Python] Résolution de problèmes de valeurs propres (généralisés) en utilisant numpy / scipy, en utilisant des bibliothèques
[Calcul scientifique / technique par Python] Résolution de l'équation différentielle ordinaire du second ordre par la méthode Numerov, calcul numérique
[Calcul scientifique / technique par Python] Calcul numérique pour trouver la valeur de la dérivée (différentielle)
[Calcul scientifique / technique par Python] Solution analytique sympa pour résoudre des équations
[Calcul scientifique / technique par Python] Solution numérique d'une équation de conduction thermique non stationnaire unidimensionnelle par méthode Crank-Nicholson (méthode implicite) et méthode FTCS (méthode de solution positive)
[Calcul scientifique / technique par Python] Fonctionnement de base du tableau, numpy
[Calcul scientifique / technique par Python] Intégration Monte Carlo, calcul numérique, numpy
[Calcul scientifique / technique par Python] Intégration numérique, loi trapézoïdale / Simpson, calcul numérique, scipy
[Calcul scientifique / technique par Python] Simulation de Monte Carlo par la méthode metropolis de la thermodynamique du système de spin ascendant 2D
[Calcul scientifique / technique par Python] Ajustement par fonction non linéaire, équation d'état, scipy
[Calcul scientifique / technique par Python] histogramme, visualisation, matplotlib
[Calcul scientifique / technique par Python] Dessin d'animation de mouvement parabolique avec locus, matplotlib
[Calcul scientifique / technique par Python] Dessin de surface courbe 3D, surface, fil de fer, visualisation, matplotlib
[Calcul scientifique / technique par Python] Résolution de l'équation de Newton unidimensionnelle par la méthode Runge-Kutta du 4ème ordre
Calcul numérique du fluide compressible par la méthode des volumes finis
[Calcul scientifique / technique par Python] Graphique logistique, visualisation, matplotlib
Calcul de la matrice d'homographie par la méthode des moindres carrés (méthode DLT)
[Calcul scientifique / technique par Python] Graphique de coordonnées polaires, visualisation, matplotlib
[Calcul scientifique / technique par Python] Tracé, visualisation, matplotlib de données 2D lues à partir d'un fichier
[Calcul scientifique / technique par Python] Dessin, visualisation, matplotlib de lignes de contour 2D (couleur), etc.
[Python] Calcul numérique d'une équation de diffusion thermique bidimensionnelle par méthode ADI (méthode implicite de direction alternative)
[Calcul scientifique / technique par Python] Résolution de l'équation de Schledinger à l'état stationnaire dans le potentiel d'oscillateur isotrope tridimensionnel par la méthode matricielle, problème des valeurs aux limites, mécanique quantique
[Calcul scientifique / technique par Python] Interpolation spline de troisième ordre, scipy
[Méthode de calcul numérique, python] Résolution d'équations différentielles ordinaires par la méthode Eular
[Calcul scientifique / technique par Python] Liste des utilisations des fonctions (spéciales) utilisées en physique en utilisant scipy
[Calcul scientifique et technique par Python] Dessin de figures fractales [Triangle de Shelpinsky, fougère de Bernsley, arbre fractal]
Calcul scientifique / technique avec Python] Dessin et visualisation d'isoplans 3D et de leurs vues en coupe à l'aide de mayavi
[Calcul scientifique / technique par Python] Exemple de visualisation de champ vectoriel, champ magnétique électrostatique, matplotlib
[Calcul scientifique / technique par Python] Transformation de Fourier à grande vitesse discrète en 3D unidimensionnelle, scipy
Histoire d'approximation de puissance par Python
[Calcul scientifique / technique par Python] Comparaison des vitesses de convergence de la méthode SOR, de la méthode Gauss-Seidel et de la méthode Jacobi pour l'équation de Laplace, équation différentielle partielle, problème des valeurs aux limites
[Calcul scientifique / technique par Python] Dérivation de solutions analytiques pour équations quadratiques et cubiques, formules, sympy
[Calcul scientifique / technique par Python] Résolution d'équations différentielles ordinaires, formules mathématiques, sympy
[Calcul scientifique / technique par Python] Résolution de l'équation de Schrödinger unidimensionnelle à l'état stationnaire par méthode de tir (1), potentiel de type puits, mécanique quantique
[Calcul scientifique / technique par Python] Tracer, visualiser, matplotlib des données 2D avec barre d'erreur
[Calcul scientifique / technique par Python] Résolution d'une équation de Schrödinger unidimensionnelle en régime permanent par méthode de tir (2), potentiel d'oscillateur harmonique, mécanique quantique
Calcul des indicateurs techniques par TA-Lib et pandas
Matrice unitaire et matrice inverse: Algèbre linéaire en Python <4>
Calcul matriciel et équations linéaires: Algèbre linéaire en Python <3>