【数值分析期末大作业/第四次大作业】

 

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

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

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

 

 

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

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

//矩阵乘法
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;
}

//矩阵乘法
vecRow multi(vec A, vecRow 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{ 1 };
    vecRow res(row1);
    double temp{};
    assert(col1 == row2 && "第一个矩阵的列和第二个矩阵的行不相等");
    for (int i = 0;i != row1;++i) {
        temp = 0;
        for (int j = 0;j != row2;++j) {
            temp += A[i][j] * B[j];
        }
        res[i] = temp;
    }
    return res;
}


//高斯法求解方程
/*
    A：左式
    B：右式
*/
vecRow equationSolver(vec A, vecRow B) {   //如果这里为vec&A vecRow&B ，则A和B的值会改变。如果没有&，则不变，改变的为A和B的副本
    int row1{ static_cast<int> (A.size()) };
    int col1{ static_cast<int>(A[0].size()) };
    int row2{ static_cast<int>(B.size()) };
    assert((row1 == row2) || (row1 > 0) && "请检查输入");
    int point;
    double temp;
    double m;
    vecRow x(row1);
    if (row1 == 1) {
        x[row1] = B[row1] / A[row1][row1];
        return x;
    }
    //一列一列来
    for (int k = 0;k != col1 - 1;++k) {
        //选主元
        point = k;
        //寻找k列中绝对值最大的数的行point
        for (int r = k;r != row1;++r) {
            if (fabs(A[point][k]) < fabs(A[r][k])) {
                point = r;
            }
        }
        //判断A是否为奇异矩阵
        assert(fabs(A[point][k]) > 1e-6 && "A为奇异矩阵");
        //换行
        temp = B[k];
        B[k] = B[point];
        B[point] = temp;
        for (int i = 0;i != col1;++i) {
            temp = A[point][i];
            A[point][i] = A[k][i];
            A[k][i] = temp;
        }
        //消去
        m = 0.0;
        for (int i = k + 1;i != col1;++i) {
            m = A[i][k] / A[k][k];
            B[i] -= m * B[k];
            for (int j = 0;j != col1;++j) {
                A[i][j] -= m * A[k][j];
            }
        }
    }
    //解方程
    double sum{ 0.0 };
    x[row1 - 1] = B[row1 - 1] / A[row1 - 1][row1 - 1];
    for (int i = row1 - 2;i >= 0;--i) {
        sum = 0.0;
        for (int j = i + 1;j != col1;++j) {
            sum += A[i][j] * x[j] / A[i][i];
        }
        x[i] = B[i] / A[i][i] - sum;
    }
    /*for (double i : x) {
        std::cout << i << "  ";
    }*/
    return x;
}


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

//计算多项式拟合的结果
/*
    c：系数矩阵
    y：数组y，需要求解得出该数组中每一个y相对应的x
*/
vecRow cal_x(vecRow c, vecRow y) {
    int row_c{ static_cast<int> (c.size()) };
    int row_y{ static_cast<int> (y.size()) };
    double temp{};
    vecRow result(row_y);
    for (int i = 0;i != row_y;++i) {
        temp = 0.0;
        for (int j = 0;j != row_c;++j) {
            temp += pow(y[i], j)*c[j];
        }
        result[i] = temp;
    }
    return result;
}

//计算x=φ(y)
double cal_x(vecRow c, double y) {
    int k{ static_cast<int> (c.size()) };
    double x;
    x = 0.0;
    for (int j = 0;j < k;j++) {
        x += pow(y, j) * c[j];
    }
    return x;
}

//计算误差
/*
    A和B分别时表格中给定的x值和计算拟合出的x值
*/
double epsilon(vecRow A, vecRow B) {
    int row_A{ static_cast<int> (A.size()) };
    int row_B{ static_cast<int> (B.size()) };
    assert(row_A == row_B, "please check the input");
    double mid{};
    for (int i = 0;i != row_A;++i) {
        mid += pow(A[i] - B[i], 2);
    }
    return mid;
}

//返回方程求解的结果vecRow c，然后具体的k可以通过计算c的大小得出
//传入x=φ(y)的x和y
/*
    x,y由题目给出
*/
vecRow qxnh(vecRow x,vecRow y) {
    int k{ 2 };
    int row_x{ static_cast<int> (x.size()) };
    int row_y{ static_cast<int> (y.size()) };
    assert(row_x == row_y, "please check the input");
    double ep{};
    vecRow result(row_x);
    while (1) {
        //构建A
        vec A(row_x, vecRow(k));
        for (int i = 0;i!= row_y;++i) {
            for (int j = 0;j != k;++j) {
                A[i][j] = pow(y[i], j);
            }
        }
        vec AT{ transpose(A) };
        //ATA
        vec ATA{ multi(AT,A) };
        //ATy,该例中，y为x。函数重载
        vecRow ATy{ multi(AT, x) };
        vecRow c{ equationSolver(ATA,ATy) };
        result = cal_x(c, y);
        ep = epsilon(result, x);
        if (ep <= 1e-6) {
            std::cout << "停止循环，此时k = " << k << ";多项式的阶数为：" << k - 1 << '\n';
            return c;
        }
        else {
            ++k;
        }
    }
}

/*
    打印出矩阵（调试用）
*/
void show(vec A) {
    for (vecRow i : A) {
        for (double j : i) {
            std::cout << j << " ";
        }
        std::cout << '\n';
    }
}

/*
    打印出矩阵（调试用）
*/
void show(vecRow A) {
    for (double i : A) {
        std::cout << i << " ";
        std::cout << '\n';
    }
}

//牛顿法解方程
/*
    已知t，y，x，求解y'
*/
double Newton(double t, double y, double x) {
    double epsilon = 1e-14;
    int M = 10000;
    double F;
    double dF;
    double ddy;
    int count;
    double dy;

    double dy_ = 0.0;
    int flag = 0;
START:
    dy = dy_;
    count = 0;
    while (1) {
        F = pow(dy, 3) - 6 * pow(t, 3) * x * pow(y, 2) * dy + pow(t, 3) * pow(x, 3) + 8 * pow(t, 6) * pow(y, 6);
        dF = 3 * pow(dy, 2) - 6 * pow(t, 3) * x * pow(y, 2);
        ddy = F / dF;
        if (fabs(ddy / dy) <= epsilon) {
            return dy;
        }
        else {
            dy -= ddy;
            count++;
        }
        if (count < M) {
            continue;
        }
        else {
            flag++;
            if (dy_ == 0) {
                dy_ -= 2;
            }
            else if (dy_ < 0) {
                dy_ = -dy_;
            }
            else {
                dy_ += 2;
                dy_ = -dy_;
            }
            if (flag >= 1000) {
                std::cout << "牛顿法失败" << '\n';
                return 0;
            }

        }
        goto START;

    }
    return 0;
}



//R-K法求解方程：已知t0和y0，多项式曲线拟合得出x0，牛顿法得出y'，可求得对应的y(t)
/*
    c:系数矩阵
    t0：t的初值
    y0：y的初值
    h：计算步长
    t_end：t的终值
    除c外，其余参数可使用默认值
*/
vecRow RK(vecRow c, double t0=0.5,double y0=8.0,double h=0.001,double t_end = 5.5) {
    int k{ static_cast<int> (c.size()) };
    double x0{ cal_x(c, y0) };
    double k1, k2, k3, k4;
    double t1, t2, t3;
    double y1, y2, y3;
    double x1, x2, x3;
    int count{ 0 };
    int flag{ 0 };
    int num{ static_cast<int> (((t_end - t0) / 0.1)+1) };
    vecRow t(num);
    vecRow y(num);
    while (1) {
        if (count % 100 == 0) {
            flag = count / 100;
            //std::cout << "flag: " << flag << '\n';
            t[flag] = t0;
            y[flag] = y0;
        }
        x0 = cal_x(c, y0);
        k1 = Newton(t0, y0, x0);
        t1 = t0 + 0.5 * h;
        y1 = y0 + 0.5 * h * k1;
        x1 = cal_x(c, y1);
        k2 = Newton(t1, y1, x1);
        y2 = y0 + 0.5 * h * k2;
        x2 = cal_x(c, y2);
        t2 = t1;
        k3 = Newton(t2, y2, x2);
        t3 = t0 + h;
        y3 = y0 + h * k3;
        x3 = cal_x(c, y3);
        k4 = Newton(t3, y3, x3);
        y0 += (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
        t0 += h;
        if (k1 == 0 || k2 == 0 || k3 == 0 || k4 == 0) {
            continue;
        }
        count++;
        if (t0 > t_end+0.001) {

            break;
        }
        else {
            continue;
        }
    }
    //std::cout << "x0: " << x0 << '\n';
    return y;
}

//复化simpson求解积分
/*
    y：y(t)
    m、h的定义见simpson公式
    n:n个等距节点

*/
double simpson(vecRow y, double h=0.1,int n=51) {
    double result;
    double ya{ y[0] };
    double yb{ y[n-1] };
    int m{ (n - 1) / 2 };
    double sum1, sum2;
    sum1 = 0.0;
    for (int i = 1;i < m + 1;i++) {
        sum1 += y[2 * i - 1];
    }
    sum2 = 0.0;
    for (int i = 1;i < m;i++) {
        sum2 += y[2 * i];
    }
    result = (h / 3) * (ya + yb + 4 * sum1 + 2 * sum2);
    return result;
}


//创建t矩阵
/*
    t0:t的初值
    n:t的个数
*/
vecRow create_t(double t0 = 0.5,int n = 51) {
    vecRow t(n);
    for (int i = 0;i != n;++i) {
        t[i] = t0 + 0.5 * i;
    }
    return t;
}


int main()
{
    vecRow x{ 0.000000,0.3149803,0.9999999,1.9655560,3.1748020,4.6050390,6.2402510,
                     8.068264,10.07937,12.26556,14.62009,17.13712,19.81156,22.63891,
                     25.61514,28.73660,31.99999,35.40231,38.94074,42.61273,46.41589 };
    //int k ;
    int row_x{ static_cast<int> (x.size() )};
    vecRow y(row_x);
    for (int i = 0;i != row_x;++i) {
        y[i] = 0.5 * i;
    }
    vecRow c{ qxnh(x,y) };
    vecRow y2{ RK(c) };
    vecRow t{ create_t()};
    int n{ 51 };
    vecRow x2(n);
    for (int i = 0;i!=51;i++) {
        x2[i] = cal_x(c,y2[i]);
        printf("t[%d]=%.10f    y[%d]=%.10f   x[%d]=%.10f\n", i, t[i], i, y2[i], i, x2[i]);
    }
    double S{ simpson(y2) };
    printf("积分值S = %.10f\n", S);
}