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

78 阅读 0 评论 52 点赞

我是靠谱客的博主含蓄小猫咪，最近开发中收集的这篇文章主要介绍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//最大迭代次数

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//最大迭代次数

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 }