L U分解在我之前写的文章里。
定义的变量有点多,但挺容易看的。
#include<stdio.h> #include<math.h> #define N 3 int main (void) { double
A[N][N] = {0}; double D[N][N] = {0}; double L[N][N] = {0}; double
U[N][N] = {0}; double C[N][N] = {0}; double H[N][N] = {0}; double
B[N][N] = {0}; double F[N][N] = {0}; double F1[N][N] = {0}; double
X[N][N] = {0}; double X1[N][N] = {0}; double X2[N][N] = {0}; double
X3[N][N] = {0}; double X4[N][N] = {0}; double sum1 = 0; int i,j,k;
printf("请输入你的矩阵A:\n"); for(i = 0;i<N;i++) {
printf("第%d行:",i+1); for(j=0;j<N;j++) {
scanf("%lf",&A[i][j]); } } printf("\n");
printf("请输入你的常数矩阵:\n"); for(i = 0;i<N;i++) {
printf("第%d行:",i+1); for(j=0;j<N;j++) {
scanf("%lf",&B[i][j]); } } printf("\n");
printf("请输入你的矩阵U:\n"); for(i = 0;i<N;i++) {
printf("第%d行:",i+1); for(j=0;j<N;j++) {
scanf("%lf",&U[i][j]); } } printf("\n");
printf("请输入你的矩阵L:\n"); for(i = 0;i<N;i++) {
printf("第%d行:",i+1); for(j=0;j<N;j++) {
scanf("%lf",&L[i][j]); } } printf("\n");
/*以迹组成的矩阵D的逆矩阵*/ for(i = 0;i<N;i++) { D[i][i] = 1/A[i][i];
} /*矩阵L+U = H*/ for(i = 0;i<N;i++) { for(j=0;j<N;j++)
{ H[i][j] = L[i][j] + U[i][j]; } }
/*以迹组成的矩阵D的逆矩阵与矩阵L+U = H的乘积矩阵C*/ for(i = 0;i<N;i++) { for(j =
0;j<N;j++) { for(k=0;k<N;k++) {
C[i][j] += D[i][k]*H[k][j]; } } }
/*以迹组成的矩阵D的逆矩阵与系数矩阵的乘积新的系数矩阵F*/ for(i = 0;i<N;i++) { for(j =
0;j<N;j++) { for(k=0;k<N;k++) {
F[i][j] += D[i][k]*B[k][j]; } } } for(i =
0;i<N;i++) { F1[i][0] = F[i][0]; } sum1 = 1; while(sum1
> 0.0001) { for(i = 0;i<N;i++) {
for(j=0;j<N;j++) { X1[i][j] = 0;
X2[i][j] = 0; X3[i][j] = 0; X4[i][j] = 0;
} } sum1 = 0; /*乘积矩阵C与初值矩阵X的新乘积矩阵X1*/
for(i = 0;i<N;i++) { for(j = 0;j<N;j++)
{ for(k=0;k<N;k++) {
X1[i][j] += C[i][k]*X[k][j]; } } }
/*新乘积矩阵X1的第一列矩阵X2*/ for(i = 0;i<N;i++) { X2[i][0]
= X1[i][0]; } /*核心雅可比迭代公式*/ for(i = 0;i<N;i++)
{ for(j = 0;j<N;j++) { X3[i][j] =
X2[i][j] + F1[i][j]; } } /*判断是否需要继续迭代*/
/*新初值矩阵X3与原始初值矩阵X的差矩阵X4,X3与X仅有第一列元素被赋值为非0*/ for(i = 0;i<N;i++)
{ for(j = 0;j<N;j++) { X4[i][j] =
X3[i][j] - X[i][j]; } } /*差矩阵X4的二范数*/ for(i
= 0;i<N;i++) { sum1 += X4[i][0]*X4[i][0]; }
sum1 = sqrt(sum1); /*把新初值矩阵赋给原始初值矩阵,代替其位置。*/ for(i =
0;i<N;i++) { for(j = 0;j<N;j++) {
X[i][j] = X3[i][j]; } } }
printf("最终利用雅可比迭代法求得的线性方程组的解是:\n"); for(i = 0;i<N;i++) {
printf("%lf\t",X3[i][0]); } printf("\n"); return 0; } That's all!?