计算方法 - 矩阵的三角分解 - 解对称正定矩阵方程组的平方根法（线性方程组的解法）

最新推荐文章于 2023-11-24 21:09:13 发布

Spikeeee-

最新推荐文章于 2023-11-24 21:09:13 发布

阅读量554

点赞数 2

分类专栏：计算方法

本文链接：https://blog.csdn.net/weixin_45919985/article/details/116056521

版权

计算方法专栏收录该内容

13 篇文章 1 订阅

订阅专栏

详细过程见《数值计算方法》
题目：
在这里插入图片描述
C++题解：

#include <iostream>
#include <cstdio>
#include <queue>
#include <cstring>
#include <cmath>
using namespace std;
#define ll long long
const int maxn = 205;
const int INF = 0x3f3f3f;
double temp[105 * 105] = {0.0};
double x[105] = {0.0}, y[105] = {0.0};
double L[105][105] = {0.0};
int n = 0, n2 = 0;

inline void IN() { //输入
    int t = 0;
    while(scanf("%lf",&temp[++t]) != EOF) {}
    while(n*(n+1) != (t-1)) {n++;}
    n2 = n * n;
}

inline double *A(int i, int j) {
    return &temp[j + (i-1) * n];
}

inline double *b(int j) {
    return &temp[j + n2];
}

inline void OUT() {
    printf("L:\n");
    for(int i = 1; i <= n; ++i) {
        for(int j = 1; j <= n; ++j) {
            printf("%.5lf ", L[i][j]);
        }
        printf("\n");
    }
    printf("y:\n");
    for(int i = 1; i <= n; ++i) {
        printf("%.5lf\n", y[i]);
    }
    printf("x:\n");
    for(int i = 1; i <= n; ++i) {
        printf("%.5lf\n", x[i]);
    }
}

inline void get_L(){
    L[1][1] = sqrt(*A(1,1));
    for(int i = 2; i <= n; ++i) {
        L[i][1] = *A(i,1)/L[1][1];
    }
    for(int j = 2; j <= n; ++j) {
        L[j][j] = *A(j,j);
        for(int k = 1; k <= j-1; ++k) 
            L[j][j] -= L[j][k] * L[j][k];
        L[j][j] = sqrt(L[j][j]);
        for(int i = j+1; i <= n; ++i) if(j != n) {
            L[i][j] = *A(i,j);
            for(int k = 1; k <= j-1; ++k) 
                L[i][j] -= L[i][k] * L[j][k];
            L[i][j] /= L[j][j];
        }
    }
}

inline void get_xy() {
    y[1] = *b(1) / L[1][1];
    for(int i = 2; i <= n; ++i) {
        y[i] = *b(i);
        for(int k = 1; k <= i-1; ++k) 
            y[i] -= L[i][k] * y[k];
        y[i] /= L[i][i];
    }
    x[n] = y[n] / L[n][n];
    for(int i = n-1; i >= 1; --i) {
        x[i] = y[i];
        for(int k = i+1; k <= n; ++k) 
            x[i] -= L[k][i] * x[k];
        x[i] /= L[i][i];
    }
}

int main() {
    IN();
    get_L();
    get_xy();
    OUT();
    return 0;
}