数值分析实验之矩阵特征值（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 import numpy as np
 3 
 4 
 5 def Solve(mat, max_itrs, min_delta):
 6     """
 7     mat 表示矩阵
 8     max_itrs 表示最大迭代次数
 9     min_delta 表示停止迭代阈值
10     """
11     itrs_num = 0
12     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 / m
24         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     int m;
 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 int main()
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 }

运行结果：

posted on 2020-05-27 10:27 ぺあ紫泪°冰封ヤ阅读(3160) 评论(0) 编辑收藏举报

公告