c语言计算矩阵特征值和特征向量-1(幂法）

77 阅读 0 评论 51 点赞

我是靠谱客的博主含蓄小猫咪，这篇文章主要介绍c语言计算矩阵特征值和特征向量-1(幂法），现在分享给大家，希望可以做个参考。

复制代码

1 #include <stdio.h>

2 #include <math.h>

3 #include <stdlib.h>

4 #define M 3
//方阵的行数 列数

5 #define ε0
0.00000001//ε0为要求的精度

6 #define N
100000//最大迭代次数

8 //函数预声明

9 void printMatrix(double a[][3], int m, int n);//矩阵的打印
 10 void printVector(double a[], int m);//向量的打印
 11 double dotVector(double a[], double b[], int m);//两个一维向量之积,结果为一个数
 12 void dotMatrVect(double a[][3], double yk0[], double uk1[], int m);//矩阵和向量点积u=a.*y，yk0对应于书上y（k-1）
 13 void unitVector(double a[], int η, int m);//向量的单位化
 14 double relaError(double lamada1, double lamada2);//计算相对误差
 15
 16 //主函数
 17 int main(void)
 18 {
 19
double a[M][M] = { { 6, -12, 6 }, { -21, -3, 24 }, { -12, -12, 51 } };//待求特征值和特征向量的矩阵
 20
double uk0[M] = { 1.0, 0.0, 0.0 };//迭代向量
 21
double uk1[M] = { 0.0, 0.0, 0.0 };//迭代向量
 22
double β0 = 0.0;//β（k-1）
 23
double β1 = 0.0;//βk
 24
double η0 = 0.0;//向量u（k-1）的二范数
 25
double ε = 0.0;//计算的精度
 26
printf("待求特征值和特征向量的矩阵A:n");
 27 
printMatrix(a, M, M);
 28
printf("n");
 29
printf("初始向量u0：n");
 30 
printVector(uk0, M);
 31
printf("n");
 32
printf("第几次计算tt uktttt yktt βkn");
 33
for (int i = 0; i < N; i++)
 34 
{
 35
printf("%dt", i);//***打印计算次数i
 36
printVector(uk0, M);//***打印uk
 37
printf("|");//***打印分隔
 38
η0 = sqrt(dotVector(uk0, uk0, M));//初始向量u0的2范数
 39
unitVector(uk0, η0, M);//将初始向量u0单位化作为y（k-1）也就是yk0
 40
printVector(uk0, M); //***打印单位化后的uk0，也就是y（k-1）
 41
dotMatrVect(a, uk0, uk1, M);//uk1 = A.*yk0;
 42
printf("|");//***打印分隔
 43
β1 = dotVector(uk0, uk1, M);//β1=y（k-1）.*uk1
 44
if (i>0)
 45 
{
 46
printf("%lf ", β1);//***打印βk
 47 
}
 48
printf("n");
 49
ε = relaError(β0, β1);
 50
//判断是否收敛
 51
if (ε < ε0) //若收敛
 52 
{
 53
printf("收敛n");
 54
break;
 55 
}
 56
else //若不收敛,则变量交换 uk0=uk1;
 57 
{
 58
//double tem = 0.0;
 59
for (int q = 0; q < M; q++)
 60 
{
 61
//uk0[q] = uk1[q];
 62
//tem = uk0[q];
 63
uk0[q] = uk1[q];
 64
uk1[q] = 0.0;//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！
 65 
}
 66
β0 = β1;
 67 
}
 68 
}
 69
 70
system("pause");
 71 }
 72
 73 //函数具体执行
 74
 75 //矩阵的打印
 76 void printMatrix(double a[][M], int m, int n)
 77 {
 78
for (int i = 0; i<m; i++)
 79 
{
 80
for (int j = 0; j<n; j++)
 81 
{
 82
printf("%lf
", a[i][j]);
 83 
}
 84
printf("n");
 85 
}
 86 }
 87 //向量的打印
 88 void printVector(double a[], int m)
 89 {
 90
for (int i = 0; i < m; i++)
 91 
{
 92
printf("%lf
", a[i]);
 93 
}
 94 }
 95 //两个一维向量之积
 96 double dotVector(double a[], double b[], int m)
 97 {
 98
double dotsum = 0.0;
 99
for (int i = 0; i < m; i++)
100 
{
101
dotsum = dotsum + a[i] * b[i];
102 
}
103
return(dotsum);
104 }
105 //矩阵和向量点积u=a.*y，yk0对应于书上y（k-1）
106 void dotMatrVect(double a[][M], double yk0[], double uk1[], int m)
107 {
108
double a1, b, c;
109
for (int i = 0; i < m; i++)
110 
{
111
uk1[i] = 0;//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！
112
for (int j = 0; j < m; j++)
113 
{
114
uk1[i] = uk1[i] + a[i][j] * yk0[j];//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！！！！
115
a1 = a[i][j];
116
b = yk0[j];
117
c = uk1[i];
118
//printf("a[%d][%d]=%lfn",i,j,a[i][j]);
119 
}
120
121 
}
122
//printVector(uk1, 3);
123 }
124 //向量的单位化
125 void unitVector(double a[], int η, int m)
126 {
127
for (int i = 0; i < m; i++)
128 
{
129
a[i] = a[i] / η;
130 
}
131 }
132 //计算误差
133 double relaError(double β1, double β2)
134 {
135
double ε;
136
ε = fabs(β2 - β1) / fabs(β2);
137
return ε;
138 }

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242

1 #include <stdio.h>

2 #include <math.h>

3 #include <stdlib.h>

4 #define M 3
//方阵的行数 列数

5 #define ε0
0.00000001//ε0为要求的精度

6 #define N
100000//最大迭代次数

7

8 //函数预声明

9 void printMatrix(double a[][3], int m, int n);//矩阵的打印
 10 void printVector(double a[], int m);//向量的打印
 11 double dotVector(double a[], double b[], int m);//两个一维向量之积,结果为一个数
 12 void dotMatrVect(double a[][3], double yk0[], double uk1[], int m);//矩阵和向量点积u=a.*y，yk0对应于书上y（k-1）
 13 void unitVector(double a[], int η, int m);//向量的单位化
 14 double relaError(double lamada1, double lamada2);//计算相对误差
 15
 16 //主函数
 17 int main(void)
 18 {
 19
double a[M][M] = { { 6, -12, 6 }, { -21, -3, 24 }, { -12, -12, 51 } };//待求特征值和特征向量的矩阵
 20
double uk0[M] = { 1.0, 0.0, 0.0 };//迭代向量
 21
double uk1[M] = { 0.0, 0.0, 0.0 };//迭代向量
 22
double β0 = 0.0;//β（k-1）
 23
double β1 = 0.0;//βk
 24
double η0 = 0.0;//向量u（k-1）的二范数
 25
double ε = 0.0;//计算的精度
 26
printf("待求特征值和特征向量的矩阵A:n");
 27 
printMatrix(a, M, M);
 28
printf("n");
 29
printf("初始向量u0：n");
 30 
printVector(uk0, M);
 31
printf("n");
 32
printf("第几次计算tt uktttt yktt βkn");
 33
for (int i = 0; i < N; i++)
 34 
{
 35
printf("%dt", i);//***打印计算次数i
 36
printVector(uk0, M);//***打印uk
 37
printf("|");//***打印分隔
 38
η0 = sqrt(dotVector(uk0, uk0, M));//初始向量u0的2范数
 39
unitVector(uk0, η0, M);//将初始向量u0单位化作为y（k-1）也就是yk0
 40
printVector(uk0, M); //***打印单位化后的uk0，也就是y（k-1）
 41
dotMatrVect(a, uk0, uk1, M);//uk1 = A.*yk0;
 42
printf("|");//***打印分隔
 43
β1 = dotVector(uk0, uk1, M);//β1=y（k-1）.*uk1
 44
if (i>0)
 45 
{
 46
printf("%lf ", β1);//***打印βk
 47 
}
 48
printf("n");
 49
ε = relaError(β0, β1);
 50
//判断是否收敛
 51
if (ε < ε0) //若收敛
 52 
{
 53
printf("收敛n");
 54
break;
 55 
}
 56
else //若不收敛,则变量交换 uk0=uk1;
 57 
{
 58
//double tem = 0.0;
 59
for (int q = 0; q < M; q++)
 60 
{
 61
//uk0[q] = uk1[q];
 62
//tem = uk0[q];
 63
uk0[q] = uk1[q];
 64
uk1[q] = 0.0;//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！
 65 
}
 66
β0 = β1;
 67 
}
 68 
}
 69
 70
system("pause");
 71 }
 72
 73 //函数具体执行
 74
 75 //矩阵的打印
 76 void printMatrix(double a[][M], int m, int n)
 77 {
 78
for (int i = 0; i<m; i++)
 79 
{
 80
for (int j = 0; j<n; j++)
 81 
{
 82
printf("%lf
", a[i][j]);
 83 
}
 84
printf("n");
 85 
}
 86 }
 87 //向量的打印
 88 void printVector(double a[], int m)
 89 {
 90
for (int i = 0; i < m; i++)
 91 
{
 92
printf("%lf
", a[i]);
 93 
}
 94 }
 95 //两个一维向量之积
 96 double dotVector(double a[], double b[], int m)
 97 {
 98
double dotsum = 0.0;
 99
for (int i = 0; i < m; i++)
100 
{
101
dotsum = dotsum + a[i] * b[i];
102 
}
103
return(dotsum);
104 }
105 //矩阵和向量点积u=a.*y，yk0对应于书上y（k-1）
106 void dotMatrVect(double a[][M], double yk0[], double uk1[], int m)
107 {
108
double a1, b, c;
109
for (int i = 0; i < m; i++)
110 
{
111
uk1[i] = 0;//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！
112
for (int j = 0; j < m; j++)
113 
{
114
uk1[i] = uk1[i] + a[i][j] * yk0[j];//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！！！！
115
a1 = a[i][j];
116
b = yk0[j];
117
c = uk1[i];
118
//printf("a[%d][%d]=%lfn",i,j,a[i][j]);
119 
}
120
121 
}
122
//printVector(uk1, 3);
123 }
124 //向量的单位化
125 void unitVector(double a[], int η, int m)
126 {
127
for (int i = 0; i < m; i++)
128 
{
129
a[i] = a[i] / η;
130 
}
131 }
132 //计算误差
133 double relaError(double β1, double β2)
134 {
135
double ε;
136
ε = fabs(β2 - β1) / fabs(β2);
137
return ε;
138 }

为啥上面的总是算的不是太精确呢？？

奥，因为二范数取的是int类型；

复制代码

1 #include <stdio.h>

2 #include <math.h>

3 #include <stdlib.h>

4 #define M 3
//方阵的行数 列数

5 #define ε0
0.00000001//ε0为要求的精度

6 #define N
100000//最大迭代次数

8 //函数预声明

9 void printMatrix(double a[][3], int m, int n);//矩阵的打印
 10 void printVector(double a[], int m);//向量的打印
 11 double dotVector(double a[], double b[], int m);//两个一维向量之积,结果为一个数
 12 void dotMatrVect(double a[][3], double yk0[], double uk1[], int m);//矩阵和向量点积u=a.*y，yk0对应于书上y（k-1）
 13 void unitVector(double a[], double η, int m);//向量的单位化
 14 double relaError(double lamada1, double lamada2);//计算相对误差
 15
 16 //主函数
 17 int main(void)
 18 {
 19
double a[M][M] = { { 6, -12, 6 }, { -21, -3, 24 }, { -12, -12, 51 } };//待求特征值和特征向量的矩阵
 20
double uk0[M] = { 2.0, 1.0, 6.0 };//迭代向量
 21
double uk1[M] = { 0.0, 0.0, 0.0 };//迭代向量
 22
double β0 = 0.0;//β（k-1）
 23
double β1 = 0.0;//βk
 24
double η0 = 0.0;//向量u（k-1）的二范数
 25
double ε = 0.0;//计算的精度
 26
printf("待求特征值和特征向量的矩阵A:n");
 27 
printMatrix(a, M, M);
 28
printf("n");
 29
printf("初始向量u0：n");
 30 
printVector(uk0, M);
 31
printf("n");
 32
printf("第几次计算tt uktttt yktt βkn");
 33
for (int i = 0; i < N; i++)
 34 
{
 35
printf("%dt", i);//***打印计算次数i
 36
printVector(uk0, M);//***打印uk
 37
printf("|");//***打印分隔
 38
η0 = sqrt(dotVector(uk0, uk0, M));//初始向量u0的2范数
 39
unitVector(uk0, η0, M);//将初始向量u0单位化作为y（k-1）也就是yk0
 40
printVector(uk0, M); //***打印单位化后的uk0，也就是y（k-1）
 41
dotMatrVect(a, uk0, uk1, M);//uk1 = A.*yk0;
 42
printf("|");//***打印分隔
 43
β1 = dotVector(uk0, uk1, M);//β1=y（k-1）.*uk1
 44
if (i>0)
 45 
{
 46
printf("%lf ", β1);//***打印βk
 47 
}
 48
printf("n");
 49
ε = relaError(β0, β1);
 50
//判断是否收敛
 51
if (ε < ε0) //若收敛
 52 
{
 53
printf("收敛n");
 54
break;
 55 
}
 56
else //若不收敛,则变量交换 uk0=uk1;
 57 
{
 58
//double tem = 0.0;
 59
for (int q = 0; q < M; q++)
 60 
{
 61
//uk0[q] = uk1[q];
 62
//tem = uk0[q];
 63
uk0[q] = uk1[q];
 64
uk1[q] = 0.0;//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！
 65 
}
 66
β0 = β1;
 67 
}
 68 
}
 69
 70
system("pause");
 71 }
 72
 73 //函数具体执行
 74
 75 //矩阵的打印
 76 void printMatrix(double a[][M], int m, int n)
 77 {
 78
for (int i = 0; i<m; i++)
 79 
{
 80
for (int j = 0; j<n; j++)
 81 
{
 82
printf("%lf
", a[i][j]);
 83 
}
 84
printf("n");
 85 
}
 86 }
 87 //向量的打印
 88 void printVector(double a[], int m)
 89 {
 90
for (int i = 0; i < m; i++)
 91 
{
 92
printf("%lf
", a[i]);
 93 
}
 94 }
 95 //两个一维向量之积
 96 double dotVector(double a[], double b[], int m)
 97 {
 98
double dotsum = 0.0;
 99
for (int i = 0; i < m; i++)
100 
{
101
dotsum = dotsum + a[i] * b[i];
102 
}
103
return(dotsum);
104 }
105 //矩阵和向量点积u=a.*y，yk0对应于书上y（k-1）
106 void dotMatrVect(double a[][M], double yk0[], double uk1[], int m)
107 {
108
double a1, b, c;
109
for (int i = 0; i < m; i++)
110 
{
111
uk1[i] = 0;//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！
112
for (int j = 0; j < m; j++)
113 
{
114
uk1[i] = uk1[i] + a[i][j] * yk0[j];//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！！！！
115
a1 = a[i][j];
116
b = yk0[j];
117
c = uk1[i];
118
//printf("a[%d][%d]=%lfn",i,j,a[i][j]);
119 
}
120
121 
}
122
//printVector(uk1, 3);
123 }
124 //向量的单位化
125 void unitVector(double a[], double η, int m)
126 {
127
for (int i = 0; i < m; i++)
128 
{
129
a[i] = a[i] / η;
130 
}
131 }
132 //计算误差
133 double relaError(double β1, double β2)
134 {
135
double ε;
136
ε = fabs(β2 - β1) / fabs(β2);
137
return ε;
138 }

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242

1 #include <stdio.h>

2 #include <math.h>

3 #include <stdlib.h>

4 #define M 3
//方阵的行数 列数

5 #define ε0
0.00000001//ε0为要求的精度

6 #define N
100000//最大迭代次数

7

8 //函数预声明

9 void printMatrix(double a[][3], int m, int n);//矩阵的打印
 10 void printVector(double a[], int m);//向量的打印
 11 double dotVector(double a[], double b[], int m);//两个一维向量之积,结果为一个数
 12 void dotMatrVect(double a[][3], double yk0[], double uk1[], int m);//矩阵和向量点积u=a.*y，yk0对应于书上y（k-1）
 13 void unitVector(double a[], double η, int m);//向量的单位化
 14 double relaError(double lamada1, double lamada2);//计算相对误差
 15
 16 //主函数
 17 int main(void)
 18 {
 19
double a[M][M] = { { 6, -12, 6 }, { -21, -3, 24 }, { -12, -12, 51 } };//待求特征值和特征向量的矩阵
 20
double uk0[M] = { 2.0, 1.0, 6.0 };//迭代向量
 21
double uk1[M] = { 0.0, 0.0, 0.0 };//迭代向量
 22
double β0 = 0.0;//β（k-1）
 23
double β1 = 0.0;//βk
 24
double η0 = 0.0;//向量u（k-1）的二范数
 25
double ε = 0.0;//计算的精度
 26
printf("待求特征值和特征向量的矩阵A:n");
 27 
printMatrix(a, M, M);
 28
printf("n");
 29
printf("初始向量u0：n");
 30 
printVector(uk0, M);
 31
printf("n");
 32
printf("第几次计算tt uktttt yktt βkn");
 33
for (int i = 0; i < N; i++)
 34 
{
 35
printf("%dt", i);//***打印计算次数i
 36
printVector(uk0, M);//***打印uk
 37
printf("|");//***打印分隔
 38
η0 = sqrt(dotVector(uk0, uk0, M));//初始向量u0的2范数
 39
unitVector(uk0, η0, M);//将初始向量u0单位化作为y（k-1）也就是yk0
 40
printVector(uk0, M); //***打印单位化后的uk0，也就是y（k-1）
 41
dotMatrVect(a, uk0, uk1, M);//uk1 = A.*yk0;
 42
printf("|");//***打印分隔
 43
β1 = dotVector(uk0, uk1, M);//β1=y（k-1）.*uk1
 44
if (i>0)
 45 
{
 46
printf("%lf ", β1);//***打印βk
 47 
}
 48
printf("n");
 49
ε = relaError(β0, β1);
 50
//判断是否收敛
 51
if (ε < ε0) //若收敛
 52 
{
 53
printf("收敛n");
 54
break;
 55 
}
 56
else //若不收敛,则变量交换 uk0=uk1;
 57 
{
 58
//double tem = 0.0;
 59
for (int q = 0; q < M; q++)
 60 
{
 61
//uk0[q] = uk1[q];
 62
//tem = uk0[q];
 63
uk0[q] = uk1[q];
 64
uk1[q] = 0.0;//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！
 65 
}
 66
β0 = β1;
 67 
}
 68 
}
 69
 70
system("pause");
 71 }
 72
 73 //函数具体执行
 74
 75 //矩阵的打印
 76 void printMatrix(double a[][M], int m, int n)
 77 {
 78
for (int i = 0; i<m; i++)
 79 
{
 80
for (int j = 0; j<n; j++)
 81 
{
 82
printf("%lf
", a[i][j]);
 83 
}
 84
printf("n");
 85 
}
 86 }
 87 //向量的打印
 88 void printVector(double a[], int m)
 89 {
 90
for (int i = 0; i < m; i++)
 91 
{
 92
printf("%lf
", a[i]);
 93 
}
 94 }
 95 //两个一维向量之积
 96 double dotVector(double a[], double b[], int m)
 97 {
 98
double dotsum = 0.0;
 99
for (int i = 0; i < m; i++)
100 
{
101
dotsum = dotsum + a[i] * b[i];
102 
}
103
return(dotsum);
104 }
105 //矩阵和向量点积u=a.*y，yk0对应于书上y（k-1）
106 void dotMatrVect(double a[][M], double yk0[], double uk1[], int m)
107 {
108
double a1, b, c;
109
for (int i = 0; i < m; i++)
110 
{
111
uk1[i] = 0;//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！
112
for (int j = 0; j < m; j++)
113 
{
114
uk1[i] = uk1[i] + a[i][j] * yk0[j];//在第二次使用前一定把uk1[i]的所有元素归零！！！！！！！！！
115
a1 = a[i][j];
116
b = yk0[j];
117
c = uk1[i];
118
//printf("a[%d][%d]=%lfn",i,j,a[i][j]);
119 
}
120
121 
}
122
//printVector(uk1, 3);
123 }
124 //向量的单位化
125 void unitVector(double a[], double η, int m)
126 {
127
for (int i = 0; i < m; i++)
128 
{
129
a[i] = a[i] / η;
130 
}
131 }
132 //计算误差
133 double relaError(double β1, double β2)
134 {
135
double ε;
136
ε = fabs(β2 - β1) / fabs(β2);
137
return ε;
138 }