用c语言用Jacobi 迭代法解三阶方程组输出方程组的解,及矩阵 L 和 U。
时间: 2024-02-25 09:56:07 浏览: 101
好的,以下是使用Jacobi迭代法解三元线性方程组的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 x0[N] = {0}; // 存储迭代初值
double eps = 1e-6; // 精度要求
double sum; // 存储累加和
int iter = 0; // 迭代次数
// Jacobi迭代法
do {
for (i = 0; i < N; i++) {
sum = a[i][N];
for (j = 0; j < N; j++) {
if (j != i) {
sum -= a[i][j] * x0[j];
}
}
x[i] = sum / a[i][i];
}
iter++;
// 判断是否满足精度要求
if (fabs(x[0]-x0[0])<eps && fabs(x[1]-x0[1])<eps && fabs(x[2]-x0[2])<eps) {
break;
}
// 更新迭代初值
for (k = 0; k < N; k++) {
x0[k] = x[k];
}
} 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循环来进行Jacobi迭代,其中迭代的初值默认为0。在每次迭代中,我们使用当前的迭代初值计算新的解向量,然后判断新旧解向量之间的差是否满足精度要求,如果满足,则跳出循环;否则,将新解向量作为下一次迭代的初值,继续迭代。最终得到的解向量即为方程组的解。
在计算完解向量之后,我们使用了一个二重循环来计算系数矩阵的LU分解,得到了下三角矩阵L和上三角矩阵U。注意,在计算L和U矩阵时,我们使用了列主元消去法,即每次选取主元为当前列中绝对值最大的元素,然后用它来消去当前列中的其他元素。
阅读全文