【数值计算】幂法与反幂法

2019-04-14 12:01发布

站内文章 / 模拟电子
16689 0 1375

幂法

求矩阵模最大的特征值及其对应特征向量
注：需要模最大特征值唯一，矩阵各列线性无关 // 幂法求特征值 // 需要保证各列线性无关 #include #include #include #include #include using namespace std; #define N 5 class MT { public: int x, y; double mt[N][N]; MT() {} MT(int x, int y) { this->x = x; this->y = y; for (int i = 0; i < this->x; i++) for (int j = 0; j < this->y; j++) mt[i][j] = 0.0; } }; MT operator * (MT a, MT b) { MT c(a.x, b.y); for (int i = 0; i < c.x; i++) for (int j = 0; j < c.y; j++) { c.mt[i][j] = 0; for (int k = 0; k < a.y; k++) c.mt[i][j] += a.mt[i][k] * b.mt[k][j]; } return c; } MT A(3, 3); MT xk(3, 1), yk(3, 1); double lastm, m; int n; double eps; void init() { n = 3; A.mt[0][0] = 2; A.mt[0][1] = 4; A.mt[0][2] = 6; A.mt[1][0] = 3; A.mt[1][1] = 9; A.mt[1][2] = 15; A.mt[2][0] = 4; A.mt[2][1] = 16; A.mt[2][2] = 36; xk.mt[0][0] = xk.mt[0][1] = xk.mt[0][2] = 1.0; m = 1.0; lastm = -9999.0; eps = 0.0001;//精度设置 } int main() { init(); while (fabs(lastm - m) > eps) { lastm = m; yk = A * xk; m = yk.mt[0][0]; for (int i = 1; i < n; i++) m = max(m, yk.mt[i][0]); for (int i = 0; i < n; i++) { xk.mt[i][0] = (1.0 / m) * yk.mt[i][0]; } for (int i = 0; i < n; i++) { cout << right << setw(6) << yk.mt[i][0] << " "; } cout << right << setw(6) << m << " "; for (int i = 0; i < n; i++) { cout << right << setw(6) << xk.mt[i][0] << " "; } cout << endl; } return 0; }

反幂法

求矩阵模最小的特征值及其对应特征向量
注：与幂法几乎完全相同 // 反幂法求特征值 // 需要保证各列线性无关 #include #include #include #include #include using namespace std; #define N 5 class MT { public: int x, y; double mt[N][N]; MT() {} MT(int x, int y) { this->x = x; this->y = y; for (int i = 0; i < this->x; i++) for (int j = 0; j < this->y; j++) mt[i][j] = 0.0; } }; MT operator * (MT a, MT b) { MT c(a.x, b.y); for (int i = 0; i < c.x; i++) for (int j = 0; j < c.y; j++) { c.mt[i][j] = 0; for (int k = 0; k < a.y; k++) c.mt[i][j] += a.mt[i][k] * b.mt[k][j]; } return c; } MT A(3, 3); MT xk(3, 1), yk(3, 1); double lastm, m; int n; double eps; void init() { n = 3; A.mt[0][0] = 2; A.mt[0][1] = 0; A.mt[0][2] = 0; A.mt[1][0] = 2; A.mt[1][1] = 3; A.mt[1][2] = 2; A.mt[2][0] = 1; A.mt[2][1] = 2; A.mt[2][2] = 3; xk.mt[0][0] = 1.0; m = 1.0; lastm = 0.0001; eps = 0.0001;//精度设置 } MT gauss(MT a, MT y) { for (int i = 0; i < n; i++) { for (int j = i + 1; j < n; j++) { double t = a.mt[j][i] / a.mt[i][i]; for (int k = i; k < n; k++) { a.mt[j][k] -= t * a.mt[i][k]; } y.mt[j][0] -= t * y.mt[i][0]; } } MT rt(3, 1); for (int i = n - 1; i >= 0; i--) { double tmp = y.mt[i][0]; for (int j = n - 1; j > i; j--) tmp -= rt.mt[j][0] * a.mt[i][j]; rt.mt[i][0] = tmp / a.mt[i][i]; } return rt; } int main() { init(); while (fabs(1.0 / lastm - 1.0 / m) > eps) { lastm = m; yk = gauss(A, xk); m = yk.mt[0][0]; for (int i = 1; i < n; i++) m = max(m, yk.mt[i][0]); for (int i = 0; i < n; i++) { xk.mt[i][0] = (1.0 / m) * yk.mt[i][0]; } for (int i = 0; i < n; i++) { cout << right << setw(10) << yk.mt[i][0] << " "; } cout << right << setw(10) << m << " "; for (int i = 0; i < n; i++) { cout << right << setw(10) << xk.mt[i][0] << " "; } cout << endl; } return 0; }

热门文章