python 矩阵特征值_数值分析实验之矩阵特征值（Python代码）

一、实验目的

1.求矩阵的部分特征值问题具有十分重要的理论意义和应用价值；

2.掌握幂法、反幂法求矩阵的特征值和特征向量以及相应的程序设计；

3.掌握矩阵QR分解

二、实验原理

幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法, 特别是用于大型稀疏矩阵。设实矩阵A=[aij]n×n有一个完全的特征向量组，其特征值为λ1 ，λ2 ,…,λn,相应的特征向量为x1 ,x2 ,…,xn.已知A的主特征值是实根，且满足条件

|λ1 |>|λ2 |≥|λ3 |≥…≥|λn |

现讨论求λ1 的方法。

幂法的基本思想是任取一个非零的初始向量ν0,由矩阵A构造一向量序列，称为迭代向量。由假设，ν0 可表示为

ν0 =α1 x1 +α2 x2 + … +αn xn (α≠0 ),

于是得到序列vk=Avk-1,序列νk /λ1k 越来越接近A的对应于λ1 的特征向量

三、实验内容

选取五级矩阵如下：

四、实验要求

利用幂法、反幂法求某个5阶矩阵的主特征值和特征向量，利用QR分解求一个5阶矩阵的所有特征值和特征向量

五、实验代码

幂法(Python)

1 #-*- coding:utf-8 -*-

2 importnumpy as np3

4

5 defSolve(mat, max_itrs, min_delta):6 """

7 mat 表示矩阵8 max_itrs 表示最大迭代次数9 min_delta 表示停止迭代阈值10 """

11 itrs_num =012 delta = float('inf')13 N =np.shape(mat)[0]14 #所有分量都为1的列向量

15 x = np.ones(shape=(N, 1))16 #x = np.array([[0],[0],[1]])

17 while itrs_num < max_itrs and delta >min_delta:18 itrs_num += 1

19 y =np.dot(mat, x)20 #print(y)

21 m =y.max()22 #print("m={0}".format(m))

23 x = y /m24 print("***********第{}次迭代*************".format(itrs_num))25 print("y =",y)26 print("m={0}".format(m))27 print("x^T为：",x)28 print(29 "——————————————分割线——————————————")30

31

32 IOS = np.array([[2, 3, 4, 5, 6],33 [4, 4, 5, 6, 7],34 [0, 3, 6, 7, 8],35 [0, 0, 2, 8, 9],36 [0, 0, 0, 1, 0]], dtype=float)37 Solve(IOS, 100, 1e-3)

运行结果：

反幂法(MATLAB)

A =[2,3,4,5,6;4,4,5,6,7;0,3,6,7,8;0,0,2,8,9;0,0,0,1,0];

I = eye(5,5);

p=13;

u0 = [1;1;1;1;1];

v = inv(A - p * I) * u0;

u = v / norm(v, inf);

i = 0;

while norm(u - u0, inf) > 1e-5

u0 = u;

v = inv(A - p * I) * u0;

u = v / norm(v, inf);

i=i+1;

end

x = p + 1 / norm(v, inf);

fprintf('迭代次数为%g\n',i);

fprintf('主特征值为%g\n',x);

u

运行结果：

QR分解(C)

1 void QR(Matrix* A, Matrix* Q, Matrix*R)2 {3 unsigned i, j, k, m;4 unsigned size;5 const unsigned N = A->row;6 matrixType temp;7

8 Matrix a, b;9

10 if (A->row != A->column || 1 != A->height)11 {12 printf("ERROE: QR() parameter A is not a square matrix!\n");13 return;14 }15

16 size =MatrixCapacity(A);17 if (MatrixCapacity(Q) !=size)18 {19 DestroyMatrix(Q);20 SetMatrixSize(Q, A->row, A->column, A->height);21 SetMatrixZero(Q);22 }23 else

24 {25 Q->row = A->row;26 Q->column = A->column;27 Q->height = A->height;28 }29

30 if (MatrixCapacity(R) !=size)31 {32 DestroyMatrix(R);33 SetMatrixSize(R, A->row, A->column, A->height);34 SetMatrixZero(R);35 }36 else

37 {38 R->row = A->row;39 R->column = A->column;40 R->height = A->height;41 }42

43 SetMatrixSize(&a, N, 1, 1);44 SetMatrixSize(&b, N, 1, 1);45

46 for (j = 0; j < N; ++j)47 {48 for (i = 0; i < N; ++i)49 {50 a.array[i] = b.array[i] = A->array[i * A->column +j];51 }52

53 for (k = 0; k < j; ++k)54 {55 R->array[k * R->column + j] = 0;56

57 for (m = 0; m < N; ++m)58 {59 R->array[k * R->column + j] += a.array[m] * Q->array[m * Q->column +k];60 }61

62 for (m = 0; m < N; ++m)63 {64 b.array[m] -= R->array[k * R->column + j] * Q->array[m * Q->column +k];65 }66 }67

68 temp = MatrixNorm2(&b);69 R->array[j * R->column + j] =temp;70

71 for (i = 0; i < N; ++i)72 {73 Q->array[i * Q->column + j] = b.array[i] /temp;74 }75 }76

77 DestroyMatrix(&a);78 DestroyMatrix(&b);79 }80

81 void Eigenvectors(Matrix* eigenVector, Matrix* A, Matrix*eigenValue)82 {83 unsigned i, j, q;84 unsigned count;85 intm;86 const unsigned NUM = A->column;87 matrixType eValue;88 matrixType sum, midSum, mid;89 Matrix temp;90

91 SetMatrixSize(&temp, A->row, A->column, A->height);92

93 for (count = 0; count < NUM; ++count)94 {95 //计算特征值为eValue，求解特征向量时的系数矩阵

96 eValue = eigenValue->array[count];97 CopyMatrix(&temp, A, 0);98 for (i = 0; i < temp.column; ++i)99 {100 temp.array[i * temp.column + i] -=eValue;101 }102

103 //将temp化为阶梯型矩阵

104 for (i = 0; i < temp.row - 1; ++i)105 {106 mid = temp.array[i * temp.column +i];107 for (j = i; j < temp.column; ++j)108 {109 temp.array[i * temp.column + j] /=mid;110 }111

112 for (j = i + 1; j < temp.row; ++j)113 {114 mid = temp.array[j * temp.column +i];115 for (q = i; q < temp.column; ++q)116 {117 temp.array[j * temp.column + q] -= mid * temp.array[i * temp.column +q];118 }119 }120 }121 midSum = eigenVector->array[(eigenVector->row - 1) * eigenVector->column + count] = 1;122 for (m = temp.row - 2; m >= 0; --m)123 {124 sum = 0;125 for (j = m + 1; j < temp.column; ++j)126 {127 sum += temp.array[m * temp.column + j] * eigenVector->array[j * eigenVector->column +count];128 }129 sum = -sum / temp.array[m * temp.column +m];130 midSum += sum *sum;131 eigenVector->array[m * eigenVector->column + count] =sum;132 }133

134 midSum =sqrt(midSum);135 for (i = 0; i < eigenVector->row; ++i)136 {137 eigenVector->array[i * eigenVector->column + count] /=midSum;138 }139 }140 DestroyMatrix(&temp);141 }142 intmain()143 {144 const unsigned NUM = 50; //最大迭代次数

145

146 unsigned N = 5;147 unsigned k;148

149 Matrix A, Q, R, temp;150 Matrix eValue;151

152 //分配内存

153 SetMatrixSize(&A, N, N, 1);154 SetMatrixSize(&Q, A.row, A.column, A.height);155 SetMatrixSize(&R, A.row, A.column, A.height);156 SetMatrixSize(&temp, A.row, A.column, A.height);157 SetMatrixSize(&eValue, A.row, 1, 1);158

159 A.array[0] = 2;A.array[1] = 3;A.array[2] = 4;A.array[3] = 5;A.array[4] = 6;160 A.array[5] = 4;A.array[6] = 4;A.array[7] = 5;A.array[8] = 6;A.array[9] = 7;161 A.array[10] = 0;A.array[11] = 3;A.array[12] = 6;A.array[13] = 7;A.array[14] = 8;162 A.array[15] = 0;A.array[16] = 0;A.array[17] = 2;A.array[18] = 8;A.array[19] = 9;163 A.array[20] = 0;A.array[21] = 0;A.array[22] = 0;A.array[23] = 1;A.array[24] = 0;164

165 //拷贝A矩阵元素至temp

166 CopyMatrix(&temp, &A, 0);167

168 //初始化Q、R所有元素为0

169 SetMatrixZero(&Q);170 SetMatrixZero(&R);171 //使用QR分解求矩阵特征值

172 for (k = 0; k < NUM; ++k)173 {174 QR(&temp, &Q, &R);175 MatrixMulMatrix(&temp, &R, &Q);176 }177 //获取特征值，将之存储于eValue

178 for (k = 0; k < temp.column; ++k)179 {180 eValue.array[k] = temp.array[k * temp.column +k];181 }182

183 //对特征值按照降序排序

184 SortVector(&eValue, 1);185

186 //根据特征值eValue，原始矩阵求解矩阵特征向量Q

187 Eigenvectors(&Q, &A, &eValue);188

189 //打印特征值

190 printf("特征值：");191 PrintMatrix(&eValue);192

193 //打印特征向量

194 printf("特征向量");195 PrintMatrix(&Q);196 DestroyMatrix(&A);197 DestroyMatrix(&R);198 DestroyMatrix(&Q);199 DestroyMatrix(&eValue);200 DestroyMatrix(&temp);201

202 return 0;203 }

运行结果：