c語言實現求一個矩陣特徵值和特徵向量
阿新 • • 發佈：2019-02-06
前言
求矩陣的特徵值，主要是用的QR分解，在我的有一次部落格裡，我已經詳細地給出了計算的過程，大家有興趣可以去看下，經過幾天的鑽研，終於完成了整個的eig演算法。下面我將把我的整個程式碼附上，有不懂的可以問我，歡迎一起討論學習！
這是對上一次的修改版，上一次寫的程式是在C++編譯環境下編譯的，所以放在c裡面編譯有些會出錯。
最後，如果有不對的地方希望大家不吝賜教，謝謝！
#include<stdio.h> #include<stdlib.h> #include <math.h> #include <stdbool.h> //定義一個結構體，用來表示一個二維的矩陣 typedef struct { int row; int column; double *data;//用來存放矩陣的元素 }Matrix; /************************************************************************ 函式功能：初始化一個矩陣 輸入：要初始化的矩陣matrix、矩陣的行row、矩陣的列column 輸出：初始化成功：true；初始化失敗：false ************************************************************************/ bool InitMatrix(Matrix *matrix, int row, int column) { int matrix_size = row*column*sizeof(double); if (matrix_size <= 0) return false; matrix->data = (double*)malloc(matrix_size);//給矩陣分配空間 if (matrix->data) { matrix->row = row; matrix->column = column; return true; } else { matrix->row = 0; matrix->column = 0; return false; } } /************************************************************************ 函式功能：打印出一個矩陣 輸入：一個矩陣matrix 輸出：無 ************************************************************************/ void PrintMatrix(Matrix *matrix) { int matrix_num = matrix->row*matrix->column; for (int i = 0; i < matrix_num; i++) { printf("%12.4g ", matrix->data[i]); if ((i + 1) % (matrix->column) == 0) printf("\n"); } printf("\n"); } /************************************************************************ 函式功能：獲取一個矩陣的大小 輸入：一個矩陣matrix 輸出：矩陣的大小size ************************************************************************/ int GetMatrixSize(Matrix *matrix) { return matrix->row*matrix->column; } /************************************************************************ 函式功能：清零，使矩陣每個元素為0 輸入：需要清零的矩陣matrix 輸出：無 ************************************************************************/ void SetMatrixZeros(Matrix *matrix) { int matrix_num = GetMatrixSize(matrix); for (int i = 0; i < matrix_num; i++) matrix->data[i] = 0; } /************************************************************************ 函式功能：判斷一個矩陣是否為空 輸入：一個矩陣matrix 輸出：為空則true，否則為false ************************************************************************/ bool IsNullMatrix(Matrix *matrix) { int matrix_num =GetMatrixSize(matrix); if ((matrix_num <= 0) || (matrix->data == NULL)) return true; else return false; } /************************************************************************ 函式功能：釋放掉一個矩陣 輸入：一個矩陣matrix 輸出：無 ************************************************************************/ void DestroyMatrix(Matrix *matrix) { if (!IsNullMatrix(matrix)) { matrix->data = NULL; matrix->row = 0; matrix->column = 0; free(matrix->data); } } /************************************************************************ 函式功能：計算一個矩陣的2範數，即求模 輸入：一個矩陣matrix 輸出：所求的範數結果norm2_ans ************************************************************************/ double MatrixNorm2(Matrix *matrix) { double norm2_ans = 0.0; int matrix_num = GetMatrixSize(matrix); for (int i = 0; i < matrix_num; i++) norm2_ans+=(matrix->data[i]) * (matrix->data[i]); norm2_ans = (double)sqrt(norm2_ans); return norm2_ans; } /************************************************************************ 函式功能：把一個矩陣複製 輸入：需要進行復制的矩陣matrix_A，複製得到的一個矩陣matrix_B 輸出：無 ************************************************************************/ void CopyMatrix(Matrix *matrix_A, Matrix *matrix_B) { if (matrix_B->row != matrix_A->row) matrix_B->row = matrix_A->row; if (matrix_B->column != matrix_A->column) matrix_B->column = matrix_A->column; int size_A = GetMatrixSize(matrix_A); for (int i = 0; i < size_A; i++) matrix_B->data[i] = matrix_A->data[i]; } /************************************************************************ 函式功能：對一個方陣A進行QR分解 輸入：需要分解的矩陣A、分解後的正交矩陣Q和上三角矩陣R 輸出：無 ************************************************************************/ void QR(Matrix *A, Matrix *Q, Matrix *R) { Matrix col_A, col_Q; InitMatrix(&col_A, A->row, 1); SetMatrixZeros(&col_A); //用來存A的每一列 InitMatrix(&col_Q, A->row, 1); SetMatrixZeros(&col_Q); //用來存Q的每一列 if (A->row != A->column) printf("A is not a square matrix!"); int A_size = GetMatrixSize(A); int Q_size = GetMatrixSize(Q); int R_size = GetMatrixSize(R); if (Q_size != A_size) { DestroyMatrix(Q); InitMatrix(Q, A->row, A->column); SetMatrixZeros(Q); } else { Q->row = A->row; Q->column = A->column; SetMatrixZeros(Q); } if (R_size != A_size) { DestroyMatrix(R); InitMatrix(R, A->row, A->column); SetMatrixZeros(R); } else { R->row = A->row; R->column = R->column; SetMatrixZeros(R); } //施密特正交化 for (int j = 0; j < A->column; j++) { for (int i = 0; i < A->column; i++)//把A的第j列存入col_A中 { col_A.data[i] = A->data[i * A->column + j]; col_Q.data[i] = A->data[i * A->column + j]; } for (int k = 0; k < j; k++)//計算第j列以前 { R->data[k * R->column + j] = 0; for (int i1 = 0; i1 < col_A.row; i1++) {//R=Q'A(Q'即Q的轉置) 即Q的第k列和A的第j列做內積 R->data[k * R->column + j] += col_A.data[i1] * Q->data[i1 * Q->column + k];//Q的第k列 } for (int i2 = 0; i2 < A->column; i2++) { col_Q.data[i2] -= R->data[k * R->column + j] * Q->data[i2 * Q->column + k]; } } double temp = MatrixNorm2(&col_Q); R->data[j * R->column + j] = temp; for (int i3 = 0; i3 < Q->column; i3++) { //單位化Q Q->data[i3 * Q->column + j] = col_Q.data[i3] / temp; } } DestroyMatrix(&col_A); DestroyMatrix(&col_Q); } /************************************************************************ 函式功能：給特徵值排序，當flag=1時，則升序，當flag=0，則降序 輸入：需要排序的序列eValue，升序還是降序的選擇flag 輸出：排序成功則返回true，否則返回false ************************************************************************/ bool SortEigenValues(Matrix *eValue, int flag) { int size = GetMatrixSize(eValue); for (int i = 0; i < size - 1; i++) { int k = i; for (int j = i + 1; j < size; j++) { if (flag == 1) { if (eValue->data[k] > eValue->data[j]) { k = j; } } else if (flag == 0) { if (eValue->data[k] < eValue->data[j]) { k = j; } } else return false; } if (k != i) { double temp; temp = eValue->data[i]; eValue->data[i] = eValue->data[k]; eValue->data[k] = temp; } } return true; } /************************************************************************ 函式功能：計算兩個矩陣相乘C=A*B 輸入：用來存計算結果的矩陣C、需要進行乘法計算的兩個矩陣A和B 輸出：計算成功則輸出true，失敗則false ************************************************************************/ bool MatrixMulMatrix(Matrix *C, Matrix *A, Matrix *B) { if ((IsNullMatrix(A)) || (IsNullMatrix(B))) return false; int A_col = A->column; int B_row = B->row; InitMatrix(C, A->row, B->column); SetMatrixZeros(C); if (A_col != B_row) { printf("A_col!=B_row!"); return false; } for (int i = 0; i < A->row; i++) { for (int j = 0; j < B->column; j++) { for (int k = 0; k < A->column; k++) C->data[i*C->row + j] += A->data[i*A->row + k] * B->data[k*B->column + j]; } } return true; } /************************************************************************ 函式功能：已知一個矩陣的特徵值求對應的特徵向量 輸入：一個矩陣A、用來存結果的特徵向量eigenVector、已知的特徵值eigenValue 輸出：無 ************************************************************************/ void Eig(Matrix *A, Matrix *eigenVector, Matrix *eigenValue) { int num = A->column; double eValue; Matrix temp; InitMatrix(&temp, A->row, A->column); //CopyMatrix(A, &temp); for (int count = 0; count < num; count++) { eValue = eigenValue->data[count];//當前的特徵值 CopyMatrix(A, &temp);//這個每次都要重新複製，因為後面會破壞原矩陣(剛開始沒注意到這個找bug找了好久。。) for (int i = 0; i < temp.row; i++) { temp.data[i * temp.column + i] -= eValue; } //將temp化為階梯型矩陣(歸一性)對角線值為一 for (int i = 0; i < temp.row - 1; i++) { double coe = temp.data[i * temp.column + i]; for (int j = i; j < temp.column; j++) { temp.data[i * temp.column + j] /= coe;//讓對角線值為一 } for (int i1 = i + 1; i1 < temp.row; i1++) { coe = temp.data[i1 * temp.column + i]; for (int j1 = i; j1 < temp.column; j1++) { temp.data[i1 * temp.column + j1] -= coe * temp.data[i * temp.column + j1]; } } } //讓最後一行為1 double sum1 = eigenVector->data[(eigenVector->row - 1) * eigenVector->column + count] = 1; for (int i2 = temp.row - 2; i2 >= 0; i2--) { double sum2 = 0; for (int j2 = i2 + 1; j2 < temp.column; j2++) { sum2 += temp.data[i2 * temp.column + j2] * eigenVector->data[j2 * eigenVector->column + count]; } sum2 = -sum2 / temp.data[i2 * temp.column + i2]; sum1 += sum2 * sum2; eigenVector->data[i2 * eigenVector->column + count] = sum2; } sum1 = sqrt(sum1);//當前列的模 for (int i = 0; i < eigenVector->row; i++) { //單位化 eigenVector->data[i * eigenVector->column + count] /= sum1; } } DestroyMatrix(&temp); } int main() { const unsigned NUM = 50; //最大迭代次數，讓資料更準確 Matrix mymatrix,temp,temp_Q,temp_R, eValue; int row,col; while (1)//VS裡面執行寫scanf會報錯，改成scanf_s就可以了 { printf("請輸入要進行計算的矩陣行與列(以逗號隔開)："); scanf_s("%d,%d", &row, &col); InitMatrix(&mymatrix, row, col); InitMatrix(&temp, row, col); SetMatrixZeros(&temp); SetMatrixZeros(&mymatrix); int num = row*col; printf("按照以行的輸入順序依次輸入矩陣內的元素，一共輸入%d個元素：", num); int data; for (int i = 0; i < num; i++) { scanf_s("%d", &data); mymatrix.data[i] = data; } printf("輸入矩陣如下：\n"); PrintMatrix(&mymatrix); CopyMatrix(&mymatrix, &temp); InitMatrix(&temp_Q, mymatrix.row, mymatrix.column); InitMatrix(&temp_R, mymatrix.row, mymatrix.column); InitMatrix(&eValue, mymatrix.row, 1); //使用QR分解求矩陣特徵值 for (int k = 0; k < NUM; ++k) { QR(&temp, &temp_Q, &temp_R); MatrixMulMatrix(&temp, &temp_R, &temp_Q);//R*Q } /*printf("Q&R:\n"); PrintMatrix(&temp_Q); PrintMatrix(&temp_R);*/ //獲取特徵值，將之儲存於eValue for (int k = 0; k < temp.column; ++k) { eValue.data[k] = temp.data[k * temp.column + k]; } SortEigenValues(&eValue, 1);//給特徵值排序，1為升序，0為降序 printf("特徵值：\n"); PrintMatrix(&eValue); Eig(&mymatrix, &temp_Q, &eValue); printf("特徵向量：\n"); PrintMatrix(&temp_Q); DestroyMatrix(&eValue); DestroyMatrix(&mymatrix); DestroyMatrix(&temp); DestroyMatrix(&temp_Q); DestroyMatrix(&temp_R); } return 0; }
下面是執行結果：
matlab的結果：
