【数值分析大作业1/第一次大作业】幂法、反幂法、矩阵行列式、谱范数（条件数）

---------------------------------------------------------------------------------------------------------------------

本博客提供的代码为在BUAA 2022秋季数值分析课提交的大作业基础上，再进行修改的产物。

10系人，提供的代码不一定完善/正确，欢迎提出意见/建议。

抄袭后果自负~课代表会查重滴~

---------------------------------------------------------------------------------------------------------------------

 

 【point 1】幂法主要用于计算矩阵的按模为最大的特征值和相应的特征向量。可通过构建带有位移的矩阵计算得出实际最大和最小的特征值。

　　　传送门：【c++/c】幂法求解绝对值最大的特征值、最大和最小的特征值 - 致命一姬 - 博客园 (cnblogs.com)

 【point 2】Doolittle分解和反幂法。Doolitle实现LU分解，以及方程的计算。反幂法用于计算矩阵的按模为最小的特征值和相应的特征向量。

　　　传送门：【c++/c】反幂法求解绝对值最小的特征值 - 致命一姬 - 博客园 (cnblogs.com)

【point 3】第二小问的求解。构造带有位移的不同A矩阵，再通过反幂法求解得到各个A矩阵按模最小的特征值。

【point 4】第三小问的求解：

　　1.谱范数：cond(A)2 = 绝对值最大的特征值/绝对值最小的特征值。利用前面计算所得的结果即可。

　　2.A行列式：

　　　　

 

 

#include <iostream>
#include <vector>
#include <math.h>
#include <cassert>
#include <iomanip>

using vec = std::vector<std::vector<double>>;
using vecRow = std::vector<double>;


/*
    定义一些全局变量
    g_r：矩阵A的下半带宽
    g_s：矩阵A的上半带宽
    g_time：迭代的最大次数，超出此数停止计算
    g_err：给定误差
*/
int g_r{ 2 };
int g_s{ 2 };
int g_time{ 5000 };
double g_err{ 10e-12 };

/*
    创建A矩阵
    b，c为题目给定量
    size为矩阵大小
*/
vec createA(double b = 0.16, double c = -0.064,int size = 501) {
    vec A(size, vecRow(size));
    double e{ exp(1) };
    for (int i = 0;i != size;++i) {
        for (int j = 0;j != size;++j) {
            if (i == j) {
                A[i][j] = (1.64 - 0.024 * (i+1)) * sin(0.2 * (i+1)) - 0.64 * exp(0.1/(i+1));
            }
            else if (fabs(i - j) == 1) {
                A[i][j] = b;
            }
            else if (fabs(i - j) == 2) {
                A[i][j] = c;
            }
            else {
                A[i][j] = 0;
            }
        }
    }
    return A;
}

/*
    给定的初始迭代矩阵u
*/
vecRow createU(int num = 1,int size=501) {
    vecRow u(size);
    for (int i = 0;i < size;++i) {
        u[i] = num;
    }
    return u;
}


/*
    计算A的转置矩阵AT
*/
vec T(vec A) {
    int row{ static_cast<int>(A.size()) };
    int col{ static_cast<int>(A[0].size()) };
    vec AT(col, vecRow(row));
    for (int i = 0;i != col;++i) {
        for (int j = 0;j != row;++j) {
            AT[i][j] = A[j][i];
        }
    }
    return AT;
}

/*
    矩阵乘法重载
*/
vec multi(vec A, vec B) {
    int row1{ static_cast<int> (A.size()) };
    int col1{ static_cast<int>(A[0].size()) };
    int row2{ static_cast<int>(B.size()) };
    int col2{ static_cast<int>(B[0].size()) };
    vec res(row1, vecRow(col2));
    double temp{};
    assert(col1 == row2 && "第一个矩阵的列和第二个矩阵的行不相等");
    for (int i = 0;i != row1;++i) {
        for (int j = 0;j != col2;++j) {
            temp = 0;
            for (int k = 0;k != col1;++k) {
                temp += A[i][k] * B[k][j];
            }
            res[i][j] = temp;
        }
    }
    return res;
}

double multi(vecRow A, vecRow B) {
    int row1{ static_cast<int> (A.size()) };
    int row2{ static_cast<int>(B.size()) };
    assert(row1 == row2 && "第一个矩阵的列和第二个矩阵的行不相等");
    double sum{0.0};
    for (int i = 0;i != row1;++i) {
        sum += A[i] * B[i];
    }
    return sum;
}

vecRow multi(vec A, vecRow y) {
    int row1{ static_cast<int> (A.size()) };
    int col1{ static_cast<int>(A[0].size()) };
    int row2{ static_cast<int>(y.size()) };
    vecRow result(row2);
    double sum{};
    assert(col1 == row2 && "第一个矩阵的列和第二个矩阵的行不相等");
    for (int i = 0;i != row1;++i) {
        sum = 0;
        for (int j = 0;j != col1;++j) {
            sum += A[i][j] * y[j];
        }
        result[i] = sum;
    }
    return result;
}

/*
    单位化向量
*/
vecRow normalize(vecRow u) {
    double mid{};
    mid = multi(u, u);
    double eta{ sqrt(mid) };
    int row{ static_cast<int> (u.size()) };
    vecRow y(row);
    for (int i = 0;i < row;++i) {
        y[i] = u[i] / eta;
    }
    return y;
}

/*
    将矩阵A变为零矩阵
*/
vec changeZero(vec A) {
    int row{ static_cast<int> (A.size()) };
    int col{ static_cast<int>(A[0].size()) };
    for (int i = 0;i!=row;++i) {
        for (int j = 0;j != col;++j) {
            A[i][j] = 0.0;
        }
    }
    return A;
}

/*
    将向量A变为零向量
*/
vecRow changeZero(vecRow A) {
    int row{ static_cast<int> (A.size()) };
    for (int i = 0;i < row;++i) {
        A[i] = 0;
    }
    return A;
}

/*
    为了节省存储量，A的带外元素不给存储，仅存A的带内元素。
    为此，设置一个二维数组C（m,n），用于存储A的带内元素。
    m = r+s+1
    数组C的第j列存放A的第j列带内元素，并使A的主对角线元素存放在C的第s+1行中。
*/
vec createC(vec A) {
    int row{ static_cast<int> (A.size()) };
    int col{ static_cast<int>(A[0].size()) };
    int m{ g_r + g_s + 1 }; //5
    int n{ col }; //501
    vec C(m, vecRow(n));
    C = changeZero(C);
    int temp{};
    int start{};
    int end{};
    int mid{ g_s + 1  };//3
    int i{};// A矩阵的行
    int k{};// C矩阵的行
    for (int j = 0;j < n;j++) {  //j为A\C矩阵的列
        if (j < mid - 1) {
            start = mid - j - 1;//实际矩阵存放位置需要减去1
            end = m - start;
i = 0;
k = start;
while ((k < m) && (i < end)) {
    C[k][j] = A[i][j];
    ++i;
    ++k;
}
        }
       else if (j >= mid - 1 && j < n - g_r) {
       start = j - g_s;
       end = j + g_r + 1;
       i = start;
       k = 0;
       while ((i < end) && (k < m)) {
           C[k][j] = A[i][j];
           ++i;
           ++k;
       }

        }
       else if (j >= n - g_r) {
       start = j - g_r;
       end = n - start + 1;
       i = start;
       k = 0;
       while ((i < n) && (k < end)) {
           C[k][j] = A[i][j];
           ++i;
           ++k;
       }
        }
       else {
       continue;
        }
    }
    return C;
}

/*
    求a，b中的最小值
*/
int min(int a, int b) {
    return (a < b) ? a : b;
}


/*
    求最大值的函数重载
*/
int max(int a, int b, int c) {
    int max1{};
    (a > b) ? max1 = a : max1 = b;
    (max1 > c) ? max1 = max1 : max1 = c;
    return max1;
}

/*
    本题需要对int ，double两种类型的数进行处理
    因此使用了功能模板
*/
template <typename T>
T max(T a, T b) {
    return (a > b) ? a : b;
}

//Doolittle 分解法 LU分解
vec dooLittleLU(vec C) {
    int row{ static_cast<int> (C.size()) };
    int col{ static_cast<int>(C[0].size()) };
    int m{ g_r + g_s + 1 };
    int n{ col };
    vecRow u(col);
    // LU分解，A->C
    //对于k=1,2，...，n计算
    int min1{};
    double sum{};
    int t0{};
    for (int k = 0;k < n;++k) {
        min1 = min((k + g_s), n - 1); //n需要-1，而k已经比原来小1了
        for (int j = k;j < min1 + 1;++j) {
            t0 = max(1 - 1, k - g_r, j - g_s);
            //t0 = (k > 2) * (k >= j) * (k - 2) + (j > 2) * (j > k) * (j - 2);
            for (int t = t0;t < k;++t) {
                //C[k - j + g_s + 1][j] = C[k - j + g_s + 1][j] - C[k - t + g_s + 1][t] * C[t - j + g_s + 1][j];
                //注意下标的转换，比公式上小1，因为从k=0开始计数
                C[k - j + g_s + 1 - 1][j] -= C[k - t + g_s + 1 - 1][t] * C[t - j + g_s + 1 - 1][j];
            }
        }
        if (k < n - 1) {
            min1 = min(k + g_r, n - 1);
            for (int i = k + 1;i < min1 + 1;++i) {
                //t0 = (k > 2) * (k >= i) * (k - 2) + (i > 2) * (i > k) * (i - 2);
                t0 = max(0, i - g_r, k - g_s);
                //sum = 0;
                for (int t = t0;t < k;++t) {
                    // std::cout << "mid: " << i - k + g_s + 1 - 1 << '\n';
                     // sum += C[i - t + g_s + 1][t] * C[t - k + g_s + 1][k];
                    C[i - k + g_s + 1 - 1][k] -= C[i - t + g_s + 1 - 1][t] * C[t - k + g_s + 1 - 1][k];
                }
                C[i - k + g_s + 1 - 1][k] = C[i - k + g_s + 1 - 1][k] / C[g_s + 1 - 1][k];
            }
        }

    }
    return C;
    //求解
}

//Doolittle 分解法 回代解方程
/*
    A(u_k) = y_k-1
*/
vecRow dooLittleSolve(vec C, vecRow uk) {
    int row{ static_cast<int> (uk.size()) };
    int n{ row };
    int t0{};
    double min1{};
    for (int i = 2 - 1;i < row + 1 - 1;++i) {
        t0 = static_cast<int>(max(1 - 1, i - g_r));
        for (int t = t0;t < i - 1 + 1;++t) {
            uk[i] -= C[i - t + g_s + 1 - 1][t] * uk[t];
        }
    }
    uk[n - 1] = uk[n - 1] / C[g_s + 1 - 1][n - 1];
    for (int i = n - 1 - 1;i >= 0;--i) {
        min1 = min(i + g_s, n-1);
        for (int t = i + 1;t < min1 + 1;++t) {
            uk[i] -= C[i - t + g_s + 1 - 1][t] * uk[t];
        }
        uk[i] /= C[g_s + 1 - 1][i];
    }
    return uk;
}

void print(vec C) {
    int row{ static_cast<int> (C.size()) };
    int col{ static_cast<int>(C[0].size()) };
    for (vecRow i : C) {
        for (double j : i) {
            std::cout << j << '\t';
        }
        std::cout << '\n';
    }
}

void print(vecRow C) {
    for (double i : C) {
        std::cout << i << '\t';
    }
    std::cout << '\n';
}


/*
    幂法求解绝对值最大的特征值
    A：要求的矩阵
    u0：初始迭代向量
*/
double mi(vec A, vecRow u0) {
    int k{ 1 };
    int row{ static_cast<int> (u0.size()) };
    vecRow y(row);
    vecRow uk(row);
    //double betaSet[static_cast<int>(g_time)]{};
    //vecRow betaSet(g_time);
   // betaSet[0]= 0;
    double beta0{};
    double betaStep{};
    double distance{};
    double beta{};
    while (k) {
        y = normalize(u0);
       // print(y);
        uk = multi(A, y);
        //betaSet[k] = multi(y, uk);
        //std::cout << betaSet[k];
        beta0 = betaStep;
        betaStep = multi(y, uk);
        distance = fabs((betaStep - beta0) / betaStep);
        if(distance<=g_err){
            beta = betaStep;
            //std::cout <<beta<< " is found\n";
            return beta;
        }
        if (k > g_time) {
            std::cout << "out of g_time" << '\n';
            return 0;
        }
        u0 = uk;
        //std::cout << "step " << k << '\n';
        ++k;
    }
    return 0;
}


/*
    幂法求解绝对值最小的特征值
    A：要求的矩阵
    u0：初始迭代向量
*/
double inverseMi(vec A, vecRow u0) {
    int row{ static_cast<int> (u0.size()) };
    vecRow y(row);
    vecRow uk(row);
    int m{ g_r + g_s + 1 };
    int n{ row };
    vec C(m, vecRow(n));
    double betak{};
    double beta0{};
    int t0{};
    double sum{};
    double min1{};
    C = createC(A);
    C = dooLittleLU(C); //LU分解
    //std::cout << "beta: " << betak << '\n';
    int k{ 1 };
    double distance{};
    while (k) {
        beta0 = betak;
        y = normalize(u0);
       // C = dooLittleLU(A); //LU分解
        uk = y;
        uk = dooLittleSolve(C, y);
        betak = 1 / multi(y, uk);
        distance = fabs((betak - beta0) / betak);
        if (distance <= g_err) {
            return betak;
        }
        if (k > g_time) {
            std::cout << "out of g_time" << '\n';
            return 0;
        }
        //std::cout << "step: " << k << '\n';
        u0 = uk;
        k++;
    }
    //return u;
    return betak;
 }

/*
    因为幂法求解的结果是带绝对值的。
    构建新的矩阵B，以求解最大/最小的特征值。
*/
vec createB(vec A, double beta) {
    int row{ static_cast<int> (A.size()) };
    int col{ static_cast<int>(A[0].size()) };
    for (int i = 0;i < row;++i) {
        for (int j = 0;j < col;++j) {
            if (i == j) {
                A[i][j] -= beta * 1;
            }
            else {
                continue;
            }
        }
    }
    return A;
}

/*
    1.幂法求解出绝对值最大的特征值
    2.如果特征值>0，则为最大的特征值，构建矩阵B来求解最小的特征值
    3.如果特征值<0，则为最小的特征值，构建矩阵B来求解最大的特征值
*/
vecRow chargeMi(vec A,double beta,vecRow u0) {
    vec B{ createB(A, beta) };
    double lambda{ mi(B, u0) + beta };
    vecRow lambdaMinMax(2);
    if (beta > 0) {
        lambdaMinMax[0] = lambda;
        lambdaMinMax[1] = beta;
    }
    else {
        lambdaMinMax[0] = beta;
        lambdaMinMax[1] = lambda;
    }
    return lambdaMinMax;
}


/*
    题目的第二小问
*/
vecRow createUK(double lambda1, double lambda501, int K) {
    vecRow miuk(K);
    miuk[0] = 0;
    for (int k = 1;k < K;++k) {
        miuk[k] = lambda1 + k *( (lambda501 - lambda1) / 40);
        //std::cout << "\nmiuk " << k << ": " << miuk[k] << '\n';
    }
    return miuk;
}


/*
    利用反幂法求解最接近的特征值
*/
void lambdaik(vecRow u0,vec A,vecRow miuk) {
    int row{ static_cast<int> (A.size()) };
    int col{ static_cast<int>(A[0].size()) };
    int K{ static_cast<int>(miuk.size()) };
    double result{};
    vec Ak(row, vecRow(col));
    for (int k = 1;k < K;++k) {
        Ak = A;
        for (int i = 0;i < row;++i) {
            Ak[i][i] -= miuk[k];
        }
        result = inverseMi(Ak, u0) + miuk[k];
        std::cout << "lambda_i" << k << ": " << result << '\n';
    }
}


/*
    求解A的谱范数（条件数）cond(A)2
    cond(A)2 = 绝对值最大的特征值/绝对值最小的特征值
*/
double cond2(vecRow lambdaMinMax,double lambdas) {
    double lambdaAbsMax{max(fabs(lambdaMinMax[0]),fabs(lambdaMinMax[1]))};
    double result{ abs(lambdaAbsMax / lambdas) };
    return result;
}

/*
    求解矩阵A的行列式detA
    1.对矩阵A进行doolittle的LU分解，A=LU；
    2.L为对角线上元素全为1的下三角矩阵，U为上三角矩阵。
    3.detL = 1，detU = 对角线元素的累乘
    4.det A = detL * detU
    */
double det(vec A) {
    int row{ static_cast<int> (A.size()) };
    int col{ static_cast<int>(A[0].size()) };
    int m{ g_r + g_s + 1 };
    int n{ row };
    vec C(m, vecRow(n));
    C = createC(A);
    C = dooLittleLU(C);
    double detA{1.0};
    for (int i = 0;i < row;++i) {
        detA *= C[g_s][i];
    }
    return detA;
}



int main()
{
    //创建初始迭代向量u0
    int size{501};
    double beta{};
    vecRow u{ createU(2,size) };
    //利用幂法计算绝对值最大的特征值
    vec A{ createA() };
    beta = mi(A, u);
    //计算最小和最大的特征值
    vecRow lambdaMinMax{ chargeMi(A,beta,u) };
    std::cout << '\n';
    std::cout << std::setprecision(12);
    std::cout << std::scientific;
    std::cout << "最小的特征值，也即lambda1为：" << lambdaMinMax[0] << '\n';
    std::cout << "最大的特征值，也即lambda501为：" << lambdaMinMax[1] << '\n';
    //计算绝对值最小的特征值
    std::cout << '\n';
    double beta2{ inverseMi(A, u)};
    std::cout << "\n绝对值最小的特征值，也即lambdas为：" << beta2 << '\n';
    //计算第二小问的内容
    std::cout << '\n';
    int K{ 40 };
    vecRow miuk(K);
    miuk = createUK(lambdaMinMax[0], lambdaMinMax[1], K);
    lambdaik(u, A, miuk);
    //计算A的谱范数
    std::cout << '\n';
    double cond2A{ cond2(lambdaMinMax,beta2) };
    std::cout << "A的谱范数：" << cond2A << '\n';
    //计算A的行列式
    std::cout << '\n';
    double detA{ det(A) };
    std::cout << "detA：" << detA<< '\n';
}