计算方法–用雅克比迭代法和高斯赛德尔迭代法求解线性方程组
AX = b
雅克比迭代法的数学基础
x = Gx + b
x = (A ^ -1) b
A = L + D + U(L是一个下三角阵,他的对角线为1,D是一个aii的对角阵,U是一个上三角阵)
(L + D + U)X = b —-> DX = LX + UX + b——> DX^(k + 1) = LX^(k + 1) + UX^(k) + b —-> X^(k + 1) = (E - A)X^(k) + D^(-1)b
高斯–赛德尔的做法(同理)
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<math.h>
#define epx 0.5 * 1e-5
#define M 1000 //最大迭代次数
#define N 5 //方程组的阶数
int main() {
double a[N][N] = {{0,0,0,0,0}, {0,8,7,0,0}, {0,6, 12, 5, 0}, {0,0, 4, 9, 3}, {0, 0, 0, 1, 2} };
double b[N] = {0, 0, -2, 8,6};
double y[N]; //结果向量
double x[N];
double norm;
memset(y, 0 , sizeof(y));
int k = 1;
int i,j;
int flag ;
double result;
do {
flag = 0;
//保存老值--strcpy();
for (i = 0; i < N ; i++) {
x[i] = y[i];
}
for (i = 1; i < N ; i++) {
result = 0;
for (j = 1; j < N; j++) {
if (i != j) {
result += a[i][j] * y[j];
}
}
y[i] = (b[i] - result) / a[i][i]; //求出解空间
}
//求y -x 的范数
norm =0.0;
for (i = 1; i < N; i++) {
if (fabs(x[i] - y[i]) > norm) {
norm = fabs(x[i] - y[i]);
}
}
if (norm < epx) {
break;
}
}while(k++ < M );
printf("########answer###########\n");
//打印结果
for (i = 1; i < N; i++) {
printf("%-4lf\n",y[i]);
}
return 0;
}