计算方法
内容见《数值计算方法》，过程有点长懒得照了，本来打算用Python写的，结果基本写完了发现Python的保留有效数字输出不自动补零，懒得重新写函数，一看又转回了c++,实在是无语
有时间再把Python的保留有效数子输出的函数补全吧
总体思想就是把矩阵拆成下三角和上三角，利用过程矩阵更快的求解
可以看看航电大兄弟的解释：=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 mathdef 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()