计算方法之用雅克比法求线性方程组

/*************************************
* 用雅克比法求线性方程组
*
*  5*x1 + 2*x2 + 1*x3 = -12
*{-1*x1 + 4*x2 + 2*x3 = 20
*  2*x1 - 3*x2 +10*x3 = 3
*
**************************************/
#include<stdio.h>
#include<math.h>
#include<conio.h>

#define N 3
#define kmax 100
#define eps 1e-5

static double aa

= { { 5, 2, 1 }, { -1, 4, 2 }, { 2, -3, 10 } };
static double bb
= { -12, 20, 3 };

int main() {
int k, i, j;
double d, sum, s, norm;
double a[N + 1][N + 1], b[N + 1], x[N + 1], y[N + 1];

for (i = 1; i <= N; i++) {
for (j = 1; j <= N; j++)
a[i][j] = aa[i - 1][j - 1];
b[i] = bb[i - 1];
}

for (i = 1; i <= N; i++) {
x[i] = 0;
}
k = 0;

do {
k++;
if (k > kmax) {
printf("\nThe iterate failed !\n");
break;
}
norm = 0.0;
for (i = 1; i <= N; i++) {
sum = 0.0;
for (j = 1; j <= N; j++) {
if (j != i)
sum += a[i][j] * x[j];
}

y[i] = (b[i] - sum) / a[i][i];
d = fabs(y[i] - x[i]);
if (norm < d)
norm = d;
}
for (i = 1; i <= N; i++)
x[i] = y[i];
} while (norm >= eps);
if (norm < eps) {
printf("\nThe result is :\n");
printf("\nk = %d \n", k);
for (i = 1; i <= N; i++)
printf("x[%d] = %f\n", i, x[i]);
}
return 0;
}