一、题目分析
给定501*501的方阵，主对角线上元素值为ai = (1.64 - 0.024i)sin(0.2i) - e(0.1/i), (i = 1, 2, ..., 501), b = 0.16, c = -0.064。首先由于题目是带状矩阵，下带宽和上带宽均为2，且是对称矩阵，因而对于后续分析处理cond2值非常方便。 求λ和λ時应当採用幂法来解决，题目的其他的λ应当採用反幂法来解决。为了求解特征值、cond2值和反幂法过程中的方程的值，可以採用带状矩阵的LU分解来解决。由于不需要存储大量的0元素，因而可以採用压缩矩阵的存储方式（有区别于CSR等压缩方法）。题目如下所示：
二、函数设计
注意到幂法和反幂法当中均存在的一些步骤：向量的点积和单位化，因而设计函数分别处理这个部分。分别设计幂法和反幂法的迭代函数。幂法函数的迭代注意，元素下标越界，所以在本函数当中使用了分类讨论的情况。反幂法首先应当做好LU分解，带状矩阵进行行压缩后的LU分解公式可以参考教材内容。之后是利用LU分解的结果求解方程组。这里，由于正常数学上的操作是行列从1开始，而计算机中的矩阵行列都是从0开始，这里处理的時候要特别注意，下标值非常容易出错。下面针对代码进行说明：
1 void uy_initialization(double u[], double y[], double *beta_k, double *beta_km, int size) {
2 memset(u, 1, size * sizeof(double));
3 memset(y, 0, size * sizeof(double));
4 *beta_k = 65535;
5 *beta_km = -1;
6 }
View Code
这段代码是对于迭代过程中产生的u向量和y向量的初始化过程。初始化可以初始化为任何值，这里使用memset函数对u赋值(1, 1, ..., 1)T。同時操作β的迭代初值。
1 double norm_2(double x[], double y[], int n) {
2 double eta = 0.0;
3 for (int i = 0; i < n; i++) {
4 eta += x[i] * y[i];
5 }
6 return eta;
7 }
View Code
这段代码是求解向量的2-范数。
1 void unitization(double x[], double y[], int n) {
2 double eta = 0.0;
3 for (int i = 0; i < n; i++) {
4 eta += x[i] * x[i];
5 y[i] = x[i];
6 }
7 eta = sqrt(eta);
8 for (int i = 0; i < n; i++) {
9 y[i] /= eta;
10 }
11 }
View Code
向量的单位化。
1 void iteration(double Matrix[][501], double b, double c, double u[], double y[]) {
2 for (int i = 0; i < 501; i++) {
3 switch (i)
4 {
5 case 0:
6 u[0] = Matrix[2][0] * y[0] + Matrix[1][1] * y[1] + Matrix[0][2] * y[2];
7 break;
8 case 1:
9 u[1] = Matrix[3][0] * y[0] + Matrix[2][1] * y[1]
10 + Matrix[1][2] * y[2] + Matrix[0][3] * y[3];
11 break;
12 case 499:
13 u[499] = Matrix[4][497] * y[497] + Matrix[3][498] * y[498]
14 + Matrix[2][499] * y[499] + Matrix[1][500] * y[500];
15 break;
16 case 500:
17 u[500] = Matrix[4][498] * y[498] + Matrix[3][499] * y[499] + Matrix[2][500] * y[500];
18 break;
19 default:
20 u[i] = Matrix[4][i - 2] * y[i - 2] + Matrix[3][i - 1] * y[i - 1] + Matrix[2][i] * y[i]
21 + Matrix[1][i + 1] * y[i + 1] + Matrix[0][i + 2] * y[i + 2];
22 break;
23 }
24 }
25
26 }
View Code
这一段是幂法的迭代计算过程。出于方便起见，没有使用标准化的处理方式，因为输入的矩阵情况已知。因为最前面和最後面的计算可能出现矩阵越界的情况，因此手动写出了相乘的方式。
以上部分属于幂法部分的计算重点，求解2-范数和单位化的函数在反幂法当中还会用到，不再重复。
1 void LU_decomposition(double Matrix[][501]) {
2 for (int k = 0; k < 501; k++) {
3 for (int j = k; j <= ((k + 2) < 500? (k + 2): 500); j++) {
4 double temp = 0;
5 for (int t = max(max(0, k - 2), j - 2); t < k; t++) {
6 temp += Matrix[k - t + 2][t] * Matrix[t - j + 2][j];
7 }
8 Matrix[k - j + 2][j] = Matrix[k - j + 2][j] - temp;
9 }
10
11 for (int i = k + 1; i <= ((k + 2) < 500? (k + 2): 500); i++) {
12 double temp = 0;
13 for(int t = max(max(0, i - 2), k - 2); t < k; t++) {
14 temp += Matrix[i - t + 2][t] * Matrix[t - k + 2][k];
15 }
16 Matrix[i - k + 2][k] = (Matrix[i - k + 2][k] - temp) / Matrix[2][k];
17 }
18 }
19 }
View Code
这里是LU分解的内容，直接在原矩阵上迭代。第0至2行属于U矩阵，第3至4行属于L矩阵（如果是分离的，则L矩阵的主对角线元素全部为1）。则由矩阵性质易知，第二行的元素是U是特征值。特別说明，这里也是採用了简化的写法，传入了二维数组而非指针。同時上下的带宽也已知是2，不用特别复杂的下标计算，直接使用即可。然而这里仍然非常容易出错，一定要多调试几回。这里可以把原矩阵的行列式也求解了。
1 void inverse_iteration(double Matrix[][501], double u[], double y[], int size) {
2 double *mid_vector;
3 mid_vector = (double *) malloc(size * sizeof(double));
4 for (int i = 0; i < 501; i++) {
5 mid_vector[i] = y[i];
6 for (int t = max(0, i - 2); t < i; t++) {
7 mid_vector[i] -= Matrix[i - t + 2][t] * mid_vector[t];
8 }
9 }
10 u[500] = mid_vector[500] / Matrix[2][500];
11 for (int i = 499; i >= 0; i--) {
12 double temp = 0;
13 for (int t = i + 1; t <= ((i + 2) < 500? (i + 2): 500); t++) {
14 temp += Matrix[i - t + 2][t] * u[t];
15 }
16 u[i] = (mid_vector[i] - temp) / Matrix[2][i];
17 }
18 }
View Code
求解完LU分解，利用其结果解Ls=y，Uu=s的线性方程组。注意这里的LU矩阵是在同一矩阵下的，是前一个函数中的在原矩阵上迭代得到的LU分解矩阵。同样注意下标。
至此，本实验当中的重要函数已经说明。
三、实验结果与总结分析
採用全1的u0作为初始迭代向量，可迭代以下结果：
| λ1 | -10.700113615150 | λ501 | 9.724634098140 | λs | -0.005557910794 | detA | 2.77279e+118 |
| cond2(A) | 1925.204273931527 | λi1 | -10.182934033146 | λi2 | -9.585707425068 | λi3 | -9.172672423928 |
| λi4 | -8.652284007898 | λi5 | -8.093483808675 | λi6 | -7.659405407692 | λi7 | -7.119684648691 |
| λi8 | -6.611764339397 | λi9 | -6.066103226595 | λi10 | -5.585101052628 | λi11 | -5.114083529812 |
| λi12 | -4.578872176865 | λi13 | -4.097829307903 | λi14 | -3.554211215751 | λi15 | -3.041090018133 |
| λi16 | -2.533970311130 | λi17 | -2.003230769564 | λi18 | -1.503557611227 | λi19 | -0.993558606008 |
| λi20 | -0.487042673885 | λi21 | 0.022317362496 | λi22 | 0.532417474207 | λi23 | 1.052898962693 |
| λi24 | 1.589445881881 | λi25 | 2.059487502799 | λi26 | 2.558075597073 | λi27 | 3.080240509307 |
| λi28 | 3.613620867692 | λi29 | 4.091378510451 | λi30 | 4.603035378279 | λi31 | 5.132924283898 |
| λi32 | 5.594906348083 | λi33 | 6.080933857027 | λi34 | 6.680354092112 | λi35 | 7.293877448126 |
| λi36 | 7.717111714236 | λi37 | 8.225220014050 | λi38 | 8.648666065194 | λi39 | 9.254200344575 |
使用不同的初始向量可以得到略有区别的迭代结果和不同的迭代收敛速度。
写函数時，可以注意函数的通用性，更多地使用指针进行内存分配和数组操作，从而方便其他类似的程序调用。