计算方法 - 矩阵的三角分解 - 杜利特尔分解法（线性方程组的解法）

内容见《数值计算方法》，过程有点长懒得照了，本来打算用Python写的，结果基本写完了发现Python的保留有效数字输出不自动补零，懒得重新写函数，一看又转回了c++,实在是无语
有时间再把Python的保留有效数子输出的函数补全吧
总体思想就是把矩阵拆成下三角和上三角，利用过程矩阵更快的求解
可以看看航电大兄弟的解释：https://blog.csdn.net/qq_37264323/article/details/104478778?utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7EOPENSEARCH%7Edefault-1.control&dist_request_id=&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7EOPENSEARCH%7Edefault-1.control
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;

void _pri(double x, long long n, bool f) { //预处理数字，保留位数，是否需要四舍五入
    double _ma = 1.0, _mi = 1.0, _t = 10.0;
    long long _n = n, i = 0, cnt = 0;
    while(_n) {
        if(_n % 2 == 1) _ma *= _t, _mi /= _t;
        _n >>= 1, _t *= _t;
    }
    if(x < 0) printf("-"), x = -x;
    if(x < 1) printf("0.");
    if(x < _mi) {
        for(i = 0; i < n; ++i) printf("0");
        return;
    }
    double t = 1.0, _x = x;
    if(x >= 1)
        while(_x >= 1)
            _x /= 10, cnt++;
    if(x < 1)
        while(_x < 1) {
            if(_x < 0.1) printf("0");
            _x *= 10;
        }
    while(_x < _ma) _x *= 10;
    long long a = (long long)_x;
    if(f == 1 && cnt <= n && a % 10 >= 5)
        if(a > 0) a += 10;
        else if(a < 0) a -= 10;
    a /= 10;
    char ans[105]; _n = 0;
    while(a) {
        ans[_n++] = a%10 + '0';
        a /= 10;
    }
    for(i = 0; i < _n; ++i) {
        if(i == cnt && cnt != 0) printf(".");
        printf("%c",ans[_n - i - 1]);
    }
}

double a[1005 * 1005] = {0.0}, b[1005] = {0.0};
double L[1005][10005] = {0.0}, U[1005][1005] = {0.0};
double x[1005] = {0.0}, y[1005] = {0.0};
int n = 0, n2 = 0;

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

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

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

void _swap(int a, int b) {
    for(int i = 1; i <= n; ++i) {
        swap(*A(a,i), *A(b,i));
    }
    swap(*B(a), *B(b));
}

inline void SIM() { //化简为三角行列式
    for(int k = 1; k <= n; ++k) {
        for(int i = k; i <= n; ++i) {
            if(i != k && k != n)
                L[i][k] = *A(i, k);
            U[k][i] = *A(k, i);
            for(int j = 1; j <= k-1; ++j) {
                if(i != k && k != n)
                    L[i][k] -= L[i][j] * U[j][k];
                U[k][i] -= L[k][j]*U[j][i];
            }
            if(i != k && k != n)
                L[i][k] /= U[k][k];
        }
    }
    y[1] = *B(1);
    for(int k = 2; k <= n; ++k) {
        y[k] = *B(k);
        for(int j = 1; j <= k-1; ++j) {
            y[k] -= L[k][j] * y[j];
        }
    }
    x[n] = y[n] / U[n][n];
    for(int k = n-1; k >= 1; --k) {
        x[k] = y[k];
        for(int j = k+1; j <= n; ++j) {
            x[k] -= U[k][j]*x[j];
        }
        x[k] /= U[k][k];
    }
}

inline void OUT() {
    printf("L:\n");
    for(int i = 1; i <= n; ++i) {
        for(int j = 1; j <= n; ++j) {
            if(i == j) L[i][j] = 1;
            if(L[i][j] == 0.0) printf("0 ");
            else _pri(L[i][j],(long long)5,1), printf(" ");
        }
        printf("\n");
    }
    printf("U:\n");
    for(int i = 1; i <= n; ++i) {
        for(int j = 1; j <= n; ++j) {
            if(U[i][j] == 0.0) printf("0 ");
            else _pri(U[i][j],(long long)5,1), printf(" ");
        }
        printf("\n");
    }
    printf("y:\n");
    for(int i = 1; i <= n; ++i) {
        if(y[i] == 0.0) printf("0\n");
        else _pri(y[i],(long long)5,1), printf("\n");
    }
    printf("x:\n");
    for(int i = 1; i <= n; ++i) {
        if(x[i] == 0.0) printf("0\n");
        else _pri(x[i],(long long)5,1), printf("\n");
    }
}

int main() {
    IN();
    SIM();
    //ANS();
    OUT();
    return 0;
}

Python不完全版本：

import numpy as np
import math

def main():
    A = np.array(input().split(" "))
    B = np.array(input().split(" "))
    n = int(math.sqrt(len(A)))
    A = A.reshape((n, n))
    L = np.identity(n)
    U = np.identity(n)
    B = B.astype(np.float64)
    A = A.astype(np.float64)
    for k in range(n):
        for i in range(k, n):
            if i != k and k != n-1:
                L[i][k] = A[i][k]
            U[k][i] = A[k][i]
            for j in range(0, k):
                if i != k and k != n-1:
                    L[i][k] -= L[i][j] * U[j][k]
                U[k][i] -= L[k][j]*U[j][i]
            if i != k and k != n-1:
                L[i][k] /= U[k][k]
    Y = np.zeros(n)
    X = np.zeros(n)
    Y[0] = B[0]
    for k in range(1, n):
        Y[k] = B[k]
        for j in range(0, k):
            Y[k] -= L[k][j] * Y[j]
    X[n-1] = Y[n-1] / U[n-1][n-1]
    for k in range(n-2, -1, -1):
        X[k] = Y[k]
        for j in range(k+1, n):
            X[k] -= U[k][j]*X[j]
        X[k] /= U[k][k]
    print("L:"), _Pri(L, n, n)
    print("U:\n", U)
    print("y:\n", Y)
    print("x:\n", X)

def _Pri(X, _i, _j):
    for i in range(_i):
        for j in range(_j):
            print('%05f' % X[i][j], end=' ')
        print('\n', end='')

if __name__ == '__main__':
    main()