用c语言用Gauss-Seidel 迭代法解如下方程组输出方程组的解,及矩阵 L 和 U。
时间: 2024-02-25 13:56:13 浏览: 23
好的,以下是使用Gauss-Seidel迭代法解三元线性方程组的C语言代码示例:
```c
#include <stdio.h>
#include <math.h>
#define N 3 // 方程组的未知数个数
int main()
{
int i, j, k;
double a[N][N+1] = {{2, -1, 1, 3}, {1, 3, -2, 5}, {4, 1, -3, -1}}; // 系数矩阵及增广矩阵
double x[N] = {0}; // 存储解向量
double eps = 1e-6; // 精度要求
double sum; // 存储累加和
int iter = 0; // 迭代次数
// Gauss-Seidel迭代法
do {
for (i = 0; i < N; i++) {
sum = a[i][N];
for (j = 0; j < N; j++) {
if (j != i) {
sum -= a[i][j] * x[j];
}
}
x[i] = sum / a[i][i];
}
iter++;
// 判断是否满足精度要求
if (fabs(a[0][0]*x[0]+a[0][1]*x[1]+a[0][2]*x[2]-a[0][N])<eps &&
fabs(a[1][0]*x[0]+a[1][1]*x[1]+a[1][2]*x[2]-a[1][N])<eps &&
fabs(a[2][0]*x[0]+a[2][1]*x[1]+a[2][2]*x[2]-a[2][N])<eps) {
break;
}
} while (iter < 100); // 最大迭代次数为100
// 输出解向量
printf("Solution:\n");
for (i = 0; i < N; i++) {
printf("x[%d] = %.6lf\n", i, x[i]);
}
// 计算系数矩阵的LU分解
double L[N][N] = {0}; // 存储下三角矩阵
double U[N][N] = {0}; // 存储上三角矩阵
for (i = 0; i < N; i++) {
// 计算U矩阵第一行元素
U[0][i] = a[0][i];
// 计算L矩阵第一列元素
L[i][0] = a[i][0] / U[0][0];
}
for (i = 1; i < N; i++) {
for (j = i; j < N; j++) {
// 计算U矩阵第i行元素
sum = 0;
for (k = 0; k < i; k++) {
sum += L[i][k] * U[k][j];
}
U[i][j] = a[i][j] - sum;
// 计算L矩阵第i列元素
sum = 0;
for (k = 0; k < i; k++) {
sum += L[j][k] * U[k][i];
}
L[j][i] = (a[j][i] - sum) / U[i][i];
}
}
// 输出L和U矩阵
printf("L matrix:\n");
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
if (j <= i) {
printf("%.6lf\t", L[i][j]);
} else {
printf("0.000000\t");
}
}
printf("\n");
}
printf("U matrix:\n");
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
if (j >= i) {
printf("%.6lf\t", U[i][j]);
} else {
printf("0.000000\t");
}
}
printf("\n");
}
return 0;
}
```
上述代码中,我们使用了一个do-while循环来进行Gauss-Seidel迭代,其中迭代的初值默认为0。在每次迭代中,我们使用当前的解向量计算新的解向量,然后判断新旧解向量之间的差是否满足精度要求,如果满足,则跳出循环;否则,将新解向量作为下一次迭代的初值,继续迭代。最终得到的解向量即为方程组的解。
在计算完解向量之后,我们使用了一个二重循环来计算系数矩阵的LU分解,得到了下三角矩阵L和上三角矩阵U。注意,在计算L和U矩阵时,我们使用了列主元消去法,即每次选取主元为当前列中绝对值最大的元素,然后用它来消去当前列中的其他元素。