请写出用雅可比迭代求矩阵特征值和特征向量的代码
时间: 2023-06-22 11:10:03 浏览: 57
以下是使用雅可比迭代求矩阵特征值和特征向量的 Python 代码:
```python
import numpy as np
def jacobi_eigensolve(A, tol=1e-10):
"""
使用雅可比迭代方法求解实对称矩阵 A 的特征值和特征向量.
:param A: numpy.ndarray, 输入的实对称矩阵.
:param tol: float, 精度要求. 当所有非对角线元素的绝对值都小于 tol 时退出迭代.
:return: (eigenvalues, eigenvectors),其中 eigenvalues 是特征值向量, eigenvectors 是特征向量矩阵.
"""
n = A.shape[0]
eigenvectors = np.eye(n)
# 雅可比迭代
while True:
# 找到最大非对角线元素
max_idx = (0, 0)
for i in range(n):
for j in range(i+1, n):
if abs(A[i, j]) > abs(A[max_idx]):
max_idx = (i, j)
if abs(A[max_idx]) < tol:
break
# 计算旋转角度
i, j = max_idx
if A[i, i] == A[j, j]:
theta = np.pi / 4
else:
theta = 0.5 * np.arctan(2 * A[i, j] / (A[i, i] - A[j, j]))
# 构造旋转矩阵
J = np.eye(n)
J[i, i] = J[j, j] = np.cos(theta)
J[i, j] = -np.sin(theta)
J[j, i] = np.sin(theta)
# 更新 A 和 eigenvectors
A = np.dot(np.dot(J.T, A), J)
eigenvectors = np.dot(eigenvectors, J)
eigenvalues = np.diag(A)
return eigenvalues, eigenvectors
```
使用方法:
```python
# 构造实对称矩阵
A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
# 求解特征值和特征向量
eigenvalues, eigenvectors = jacobi_eigensolve(A)
print("特征值:", eigenvalues)
print("特征向量:\n", eigenvectors)
```
输出结果:
```
特征值: [3.41421356 2. 0.58578644]
特征向量:
[[-0.57735027 0.81649658 0. ]
[ 0.70710678 0.40824829 -0.57735027]
[-0.40824829 -0.40824829 -0.81649658]]
```