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

幂法

求矩阵模最大的特征值及其对应特征向量
注：需要模最大特征值唯一，矩阵各列线性无关

// 幂法求特征值
// 需要保证各列线性无关
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <iomanip>

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 <iostream>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <iomanip>

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;
}