最近在上计算方法这门课，要求是用MATLAB做练习题，但是我觉得C语言也很棒棒啊~
题目：
高斯消元法实际上可以看成是将系数矩阵A分解为一个单位上三角矩阵L和一个下三角矩阵U的乘积，当高斯消元无法使用的时候当然LU分解也不能用，只要A各阶顺序主子式不为0就行了。
使用VS，代码如下：
//使用LU分解法求解线性方程组#include "stdafx.h"#include<stdlib.h>#include "math.h"double **A, *b, *x, *y, **L, **U;unsigned int RANK = 4;unsigned int makematrix(){unsigned int r, c;printf("请输入矩阵行列数，用空格隔开：");scanf_s("%d %d", &r, &c);A = (double**)malloc(sizeof(double*)*r);//创建一个指针数组，把指针数组的地址赋值给a ,*r是乘以r的意思for (int i = 0; i < r; i++)A[i] = (double*)malloc(sizeof(double)*c);//给第二维分配空间for (int i = 0; i < r; i++) {for (int j = 0; j < c; j++)A[i][j] = 0.0;}b = (double*)malloc(sizeof(double)*r);for (int i = 0; i < r; i++){b[i] = 0.0;}x = (double*)malloc(sizeof(double)*c);for (int i = 0; i < c; i++){x[i] = 0.0;}L = (double**)malloc(sizeof(double*)*r);//创建一个指针数组，把指针数组的地址赋值给a ,*r是乘以r的意思for (int i = 0; i < r; i++)L[i] = (double*)malloc(sizeof(double)*c);//给第二维分配空间for (int i = 0; i < r; i++) {for (int j = 0; j < c; j++)L[i][j] = 0.0;}U = (double**)malloc(sizeof(double*)*r);//创建一个指针数组，把指针数组的地址赋值给a ,*r是乘以r的意思for (int i = 0; i < r; i++)U[i] = (double*)malloc(sizeof(double)*c);//给第二维分配空间for (int i = 0; i < r; i++) {for (int j = 0; j < c; j++)U[i][j] = 0.0;}y = (double*)malloc(sizeof(double)*c);for (int i = 0; i < c; i++){y[i] = 0.0;}return r;}void getmatrix(void)//输入矩阵并呈现{printf("请按行从左到右依次输入系数矩阵A，不同元素用空格隔开\n");for (int i = 0; i < RANK; i++){for (int j = 0; j<RANK; j++){scanf_s("%lf", &A[i][j]);}}printf("系数矩阵如下\n");for (int i = 0; i < RANK; i++){for (int j = 0; j<RANK; j++){printf("%g\t", A[i][j]);}printf("\n");}printf("请按从上到下依次输入常数列b，不同元素用空格隔开\n");for (int i = 0; i<RANK; i++){scanf_s("%lf", &b[i]);}printf("常数列如下\n");for (int i = 0; i<RANK; i++){printf("%g\t", b[i]);}printf("\n");}void LUrush_calculation(void)//LU法解线性方程组{double get_add=0.0;printf("利用以上A与b组成的增广阵进行LU分解法计算方程组\n");for (int i = 0; i < RANK; i++){U[0][i] = A[0][i];L[i][i] = 1;}for (int i = 1; i < RANK; i++){L[i][0] = A[i][0]/U[0][0];}for (int i = 1; i < RANK; i++){for (int j = 0; j<RANK; j++){if (i <= j){get_add = 0.0;for (int k = 0; k <= i - 1; k++) get_add = get_add + L[i][k] * U[k][j];U[i][j] = A[i][j] - get_add;}if (i < j){get_add = 0.0;for (int k = 0; k <= i - 1; k++) get_add = get_add + L[j][k] * U[k][i];L[j][i] = (A[j][i] - get_add) / U[i][i];}}}printf("由系数矩阵A得到的L U矩阵如下\n");for (int i = 0; i < RANK; i++){for (int j = 0; j<RANK; j++){printf("%g\t", L[i][j]);}printf("\t\t");for (int j = 0; j<RANK; j++){printf("%g\t", U[i][j]);}printf("\n");}printf("利用 Ly = d求解y，解得：\n");y[0] = b[0];for (int i = 1; i<RANK; i++){get_add = 0.0;for (int k = 0; k <= i - 1; k++) get_add = get_add + L[i][k] *y[k];y[i] = b[i] - get_add;}for (int i = 0; i<RANK; i++){printf("y%d = %g\n", i + 1, y[i]);}printf("利用 Ux = y求解x，解得：\n");x[RANK - 1] = y[RANK - 1] / U[RANK - 1][RANK - 1];for (int i = RANK - 2; i >= 0; i--){get_add = 0.0;for (int k = i+1; k <= RANK - 1; k++) get_add = get_add + U[i][k] * x[k];x[i] = (y[i] - get_add) / U[i][i];}for (int i = 0; i<RANK; i++){printf("x%d = %g\n", i + 1, x[i]);}printf("计算完成，按回车退出程序或按1重新输入矩阵\n");}int main(){_again:RANK = makematrix();getmatrix();LUrush_calculation();getchar();if ('1' == getchar())goto _again;return 0;}
按设计的提示老老实实 输入题目的系数矩阵和常数向量后，得到运行结果：
可以看出来，X4其实就是0，但是显示不是0，这并不是程序有问题，而是这就是直接法求线性方程组会失真的体现。
