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("第几次计算\t\t uk\t\t\t\t yk\t\t βk\n"); 33 for (int i = 0; i < N; i++) 34 { 35 printf("%d\t", 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]=%lf\n",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("第几次计算\t\t uk\t\t\t\t yk\t\t βk\n"); 33 for (int i = 0; i < N; i++) 34 { 35 printf("%d\t", 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]=%lf\n",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 }
精确结果是 绝对值最大的特征值为45，对应的特征向量为（0，-0.5，-1）
