以下算法的适用条件：A的各阶主子式不为零

另外还可以采用

直接法：

消元法：Gauss-Jordan消元法，

分解法：Dolittle分解 （我用的是Courant分解法），追赶法，对称正定矩阵的LDL‘分解

----------

迭代法：

Jacobi迭代

Gauss-Seidel迭代

松弛迭代

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

你上网可以搜索一下，或者看看数值计算方面的书

OK, 你看看这个, 另外还加了注释 ：

Courant分解算法：

aX = b, 作 A=LU, L是下三角矩阵， U是上三角矩阵

即L =

| L11

| L21 L22

| L31 L32 L33

| ..............

| Ln1 Ln2 ........Lnn

U =

| 1 U12 ..... U1n

| 空格 1 ..... U2n

| 空格 空格 ........

| 空格 空格 空格 空格 空格1

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

aX = b -----> LUX = b

记 UX = y,

由Ly = b得到

因为无法输出数学符号，以下采用[i, j]Ai 表示对Ai从i到j求和

yi = (bi - [j=1, i-1]Lij yj) / Lii i = 1, 2, ..., n

由UX = y得到

xi = yi - [j=i+1, n]uij xj j = n, n-1, ..., 2, 1

你在纸上验证一下就明白了

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

以下采用Courant分解 解 aX = b, 经检查，程序运行正确

这是运行结果：

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

Input n value(dim of Ax=b): 3

Now input the matrix a(i, j), i, j = 0, ..., 2:

1 2 1 -2 -1 -5 0 -1 6

Now input the matrix b(i), i = 0, ..., 2:

24 -63 50

Solve...x_i =

7.000000

4.000000

9.000000

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

#include "stdafx.h"

#include

#include

#define MAX_N 20

int main(int argc, char* argv[])

{

int n; // 未知数个数

int i, j, k;

static double a[MAX_N][MAX_N], b[MAX_N], x[MAX_N], y[MAX_N];

static double l[MAX_N][MAX_N], u[MAX_N][MAX_N];

printf("\nInput n value(dim of Ax=b): ");

scanf("%d", &n);

if(n >MAX_N)

{

printf("The input n is larger than MAX_N, please redefine the MAX_N.\n");

return 1;

}

if(n <= 0)

{

printf("Please input a number between 1 and %d.\n", MAX_N);

return 1;

}

// {{ 程序输入

printf("Now input the matrix a(i, j), i, j = 0, ..., %d:\n", n-1);

for (i=0; i

for (j=0; j

scanf("%lf", &a[i][j]);

printf("Now input the matrix b(i), i = 0, ..., %d:\n", n-1);

for(i=0; i

scanf("%lf", &b[i]);

// }} 程序输入

for(i=0; i

u[i][i] = 1; //

for(k=0; k

{

for(i=k; i

{

l[i][k] = a[i][k];

for(j=0; j<=k-1; j++)

l[i][k] -= (l[i][j]*u[j][k]);

}

for(j=k+1; j

{

u[k][j] = a[k][j];

for(i=0; i<=k-1; i++)

u[k][j] -= (l[k][i]*u[i][j]);

u[k][j] /= l[k][k];

}

}

for(i=0; i

{

y[i] = b[i];

for(j=0; j<=i-1; j++)

y[i] -= (l[i][j]*y[j]);

y[i] /= l[i][i];

}

for(i=n-1; i>=0; i--) // 解UX = Y

{

x[i]=y[i];

for(j=i+1; j

x[i] -= (u[i][j]*x[j]);

}

printf("Solve...x_i = \n"); // 输出结果

for(i=0; i

printf("%f\n", x[i]);

return 0;

}