import numpy as np
import matplotlib.pyplot as plt
On rappelle l'algorithme de recherche de la plus grande valeur propre de $A\in M_n(\mathbb{C})$. On pose $x_0\in\mathbb{C}^n$ puis $$x_{k+1}=\frac{Ax_k}{\Vert A x_{k}\Vert_2}.$$ Quand on arrête le processus à $k_0$, on pose alors $x_1=\frac{x_{k_0}}{\Vert x_{k_0}\Vert_2}$ et $$\lambda_1=\frac{\langle Ax_{k_0},x_{k_0}\rangle}{\Vert x_{k_0}\Vert_2^2}.$$
On considère la matrice $$A=\begin{pmatrix} 2&-1&0& 0&\ldots&0\\ -1&2&-1 &0&\ddots&\vdots\\ 0&\ddots&\ddots&\ddots&\ddots&0\\ 0&\ddots&\ddots&\ddots&\ddots&0\\ \vdots&\ddots&0&-1 &2&-1\\ 0&\dots&0&0&-1&2 \end{pmatrix}.$$
def create_matrix(n):
A = np.zeros((n, n))
for i in range(n):
A[i, i] = 2
if i > 0:
A[i, i - 1] = -1
A[i - 1, i] = -1
return A
n=5
A = create_matrix(n)
print(A)
[[ 2. -1. 0. 0. 0.] [-1. 2. -1. 0. 0.] [ 0. -1. 2. -1. 0.] [ 0. 0. -1. 2. -1.] [ 0. 0. 0. -1. 2.]]
eig
pour trouver une approximation des valeurs propres $\lambda_1>\dots>\lambda_n$ de la matrice $A$. Recopier $\lambda_1$ et $\lambda_n$.n = 10
A = create_matrix(n)
eigenvalues, eigenvectors = np.linalg.eig(A)
lambda_1 = max(eigenvalues)
lambda_n = min(eigenvalues)
print(lambda_1)
print(lambda_n)
print(eigenvalues)
3.9189859472289967 0.08101405277100575 [3.91898595 3.68250707 3.30972147 2.83083003 2.28462968 1.71537032 0.08101405 0.31749293 0.69027853 1.16916997]
def power_method(A, x0, tol=1e-6, max_iter=1000):
x = x0
for i in range(max_iter):
Ax = np.dot(A, x)
xnew = Ax / np.linalg.norm(Ax)
if np.linalg.norm(xnew - x) < tol:
break
x=xnew
lambda_max = np.dot(x, A @ x)
return lambda_max, i + 1
w = np.arange(1, n + 1)
x0 = w / np.linalg.norm(w)
lambda_max, iterations = power_method(A, x0)
print(f"\nMéthode de la puissance pour lambda_max : {lambda_max}")
print(f"Iterations nécessaires : {iterations}")
Méthode de la puissance pour lambda_max : 3.9189859471692072 Iterations nécessaires : 190
w_sin = np.sin(2 * np.pi * np.arange(1, n + 1) / (n + 1))
x0_sin = w_sin / np.linalg.norm(w_sin)
lambda_max_sin, iterations_sin = power_method(A, x0_sin)
print(f"\nMéthode de la puissance pour lambda_max avec w_sin : {lambda_max_sin}")
print(f"Iterations nécessaires avec w_sin : {iterations_sin}")
print(x0_sin)
Méthode de la puissance pour lambda_max avec w_sin : 0.31749293433763764 Iterations nécessaires avec w_sin : 1 [ 0.23053002 0.38786839 0.42206128 0.3222527 0.12013117 -0.12013117 -0.3222527 -0.42206128 -0.38786839 -0.23053002]
n_values = np.arange(5, 101, 5)
iterations_list = []
ratio_eigen=[]
for n in n_values:
A_n = create_matrix(n)
w_n = np.arange(1,n+1)
x0_n = w_n / np.linalg.norm(w_n)
_, iterations_n = power_method(A_n, x0_n,max_iter=10000)
iterations_list.append(iterations_n)
eigenvalues, eigenvectors = np.linalg.eig(A_n)
eigen=sorted(eigenvalues,reverse=True)
lambda_1=eigen[0]
lambda_2 = eigen[1]
ratio_eigen.append(lambda_1/lambda_2)
# Tracé du graphique
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('n')
ax1.set_ylabel('Iterations', color=color)
ax1.plot(n_values, iterations_list, color=color, marker='o')
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx()
color = 'tab:blue'
ax2.set_ylabel('Lambda_1 / Lambda_2', color=color)
ax2.plot(n_values, ratio_eigen, color=color, marker='x')
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout()
plt.title('Nombre d\'itérations et Ratio en fonction de n')
plt.show()
Pourquoi cet algorithme doit permettre de trouver $\lambda_n$ ? Reprendre la question 3. pour calculer $\lambda_n$ en implémentant l'algorithme précédent.
def inverse_power_method(A, x0, tol=1e-6, max_iter=1000):
x = x0
for i in range(max_iter):
Ainvx = np.linalg.solve(A, x)
x_new = Ainvx / np.linalg.norm(x)
if np.linalg.norm(x_new -x) < tol:
break
x=x_new
Ainvx= np.linalg.solve(A, x)
lambda_maxinv = np.dot(x,Ainvx)/np.linalg.norm(Ainvx)**2
lambda_min=lambda_maxinv
return lambda_min, i + 1
lambda_min, iterations_min = inverse_power_method(A, x0)
print(f"\nMéthode itérative pour lambda_min : {lambda_min}")
print(f"Iterations nécessaires : {iterations_min}")
eigenvalues,_ =np.linalg.eig(A)
print(f"Plus petite valeur propre : {min(eigenvalues)}")
Méthode itérative pour lambda_min : 0.08101405277100522 Iterations nécessaires : 13 Plus petite valeur propre : 0.08101405277100575
def custom_matrix(n):
matrix = np.zeros((n, n))
for i in range(n):
for j in range(n):
matrix[i, j] = (1/2) / (n - i - j + 3/2)
return matrix
8.Reprendre les questions 3, 5 et 6 avec cette matrice. Pourquoi le nombre d'itérations est-il plus élevé ? (Comparer par exemple les valeurs propres de cette matrice par rapport à celle de la question 3.).
## Question 3.
n=10
Atilde=custom_matrix(n)
w = np.ones(n)
x0 = w / np.linalg.norm(w)
lambda_max, iterations = power_method(Atilde, x0)
print(f"\nMéthode de la puissance pour lambda_max : {lambda_max}")
print(f"Iterations nécessaires : {iterations}")
#### question 5
n_values = np.arange(5, 101, 5)
iterations_list = []
ratio_eigen=[]
for n in n_values:
A_n = custom_matrix(n)
w_n = np.arange(1,n+1)
x0_n = w_n / np.linalg.norm(w_n)
_, iterations_n = power_method(A_n, x0_n,max_iter=10000)
iterations_list.append(iterations_n)
eigenvalues, eigenvectors = np.linalg.eig(A_n)
eigen=sorted(abs(eigenvalues),reverse=True)
lambda_1=eigen[0]
lambda_2 = eigen[1]
ratio_eigen.append(lambda_1/lambda_2)
# Tracé du graphique
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('n')
ax1.set_ylabel('Iterations', color=color)
ax1.plot(n_values, iterations_list, color=color, marker='o')
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx()
color = 'tab:blue'
ax2.set_ylabel('Lambda_1 / Lambda_2', color=color)
ax2.plot(n_values, ratio_eigen, color=color, marker='x')
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout()
plt.title('Nombre d\'itérations et Ratio en fonction de n')
plt.show()
#### question 6
lambda_min, iterations_min = inverse_power_method(Atilde, x0)
print(f"\nMéthode itérative pour lambda_min : {lambda_min}")
print(f"Iterations nécessaires : {iterations_min}")
eigenvalues,_ =np.linalg.eig(Atilde)
print(f"Plus petite valeur propre : {min(abs(eigenvalues))}")
Méthode de la puissance pour lambda_max : 1.0372923569170385 Iterations nécessaires : 1000
Méthode itérative pour lambda_min : 2.608679599748989e-05 Iterations nécessaires : 6 Plus petite valeur propre : 2.608679599760921e-05
On revient maintenant sur l matrice $$A=\begin{pmatrix} 2&-1&0& 0&\ldots&0\\ -1&2&-1 &0&\ddots&\vdots\\ 0&\ddots&\ddots&\ddots&\ddots&0\\ 0&\ddots&\ddots&\ddots&\ddots&0\\ \vdots&\ddots&0&-1 &2&-1\\ 0&\dots&0&0&-1&2 \end{pmatrix},$$ avec $n=10$. On souhaite estimer l'ensemble des valeurs propre de $A$. Pour cela on remarquera que comme $A$ est symétrique, ses vecteurs propres forment une base orthonormée. Ainsi, on peut facilement montrer que, en notant $v_{1} , \ldots , v_{n}$ des vecteurs propres unitaires asociés à $\lambda_{1}> ...> \lambda_{n}$, la matrice $$ A^{(1)} = A - \lambda_{1} v_{1}v_{1}^{T} $$ a $0 , \lambda_{2} , \ldots , \lambda_{n}$ comme valeurs propres, et $v_{1} , \ldots v_{n}$ comme vecteurs propres associés. De manière générale, la matrice $$ A^{(\ell)} = A - \sum_{i=1}^{\ell}\lambda_{i} v_{i}v_{i}^{T} = A^{(\ell-1)} - \lambda_{\ell}v_{\ell}v_{\ell}^{T} $$ a pour valeurs propres $0 , \ldots , 0 , \lambda_{\ell +1} , \ldots , \lambda_{n}$ et pour vecteurs propres $v_{1} , \ldots v_{n}$. Estimer toutes les valeurs propres. Pour assurer que l'intialisation vérifie les bonnes propriétés, on prendra $x_{0}= w/\|w\|$ avec $w$ tiré selon une loi normale grâce à la fonction np.random.normal. On pourra également utiliser la fonction np.outer pour le produit de vecteurs. Enfin, on pourra utiliser la fonction suivante:
def pow_method(A, x0, tol=1e-6, max_iter=1000):
x = x0
for i in range(max_iter):
Ax = np.dot(A, x)
xnew = Ax / np.linalg.norm(Ax)
if np.linalg.norm(xnew - x) < tol:
break
x=xnew
lambda_max = np.dot(x, A @ x)/np.linalg.norm(x)**2
return lambda_max, xnew/np.linalg.norm(xnew)
n=10
w = np.random.normal(0, 1,size=n)
x0=w/(np.linalg.norm(w))
x0@x0.T
A_n = create_matrix(n)
Ai=A
lambdalast=0
vlast=np.zeros(n)
listvect=[]
for i in range(n):
Ai=Ai-lambdalast*np.outer(vlast,vlast.T)
lambdalast, vlast = pow_method(Ai, x0)
listvect.append(lambdalast)
eigenvalues, eigenvectors = np.linalg.eig(A)
eigenvalues = sorted(eigenvalues, reverse=True)
error=np.linalg.norm(listvect-np.array(eigenvalues))
print(listvect)
print(eigenvalues)
print(error)
[3.9189859471709587, 3.682507065681182, 3.309721467902276, 2.8308300260103794, 2.2846296765490823, 1.7153703234537963, 1.1691699739969248, 0.6902785321104195, 0.31749293433772013, 0.08101405277102347] [3.9189859472289963, 3.6825070656623637, 3.309721467890574, 2.8308300260037744, 2.2846296765465772, 1.71537032345343, 1.1691699739962271, 0.6902785321094301, 0.3174929343376377, 0.08101405277100472] 6.253748155096576e-11