C语言实现雅克比迭代法求根
雅克比迭代法求根
C语言实现雅克比迭代法求根问题描述算法思想C语言程序实验结果
问题描述
设方程组
A
x
=
b
Ax = b
Ax=b的系数矩阵
A
A
A非奇异 ，且
a
i
i
≠
0
{a_{ii}} \ne 0
aii=0将
A
A
A分裂为：
A
=
D
+
L
+
U
A = D + L + U
A=D+L+U
D
=
d
i
a
g
(
a
11
,
a
22
,
⋯
a
n
n
)
D = diag({a_{11}},{a_{22}}, \cdots {a_{nn}})
D=diag(a11,a22,⋯ann),
L
L
L和
U
U
U分别是将
A
A
A对角元素置为
0
0
0的下三角和上三角矩阵，进行一系列转化形成迭代格式可求解方程组。
算法思想
将方程组
∑
i
=
1
i
=
n
a
i
j
x
i
=
b
i
,
i
=
1
,
2...
,
n
\sum_{i=1}^{i=n}{{a_{ij}}{x_i}} = {b_i}, {i = 1,2...,n}
i=1∑i=naijxi=bi,i=1,2...,n 乘以
1
a
i
i
\frac{1}{{{a_{ii}}}}
aii1，得到等价方程组
x
i
=
1
a
i
i
(
b
i
−
∑
i
=
1
i
=
n
b
i
−
a
i
j
x
i
)
{x_i} = \frac{1}{{{a_{ii}}}}({b_i} - \sum_{i=1}^{i=n}{b_i}- {{a_{ij}}{x_i}} )
xi=aii1(bi−i=1∑i=nbi−aijxi) 简记为
x
=
B
x
+
f
x=Bx+f
x=Bx+f. 其中，
B
=
I
−
D
−
1
=
−
D
−
1
(
L
+
U
)
,
f
=
D
−
1
b
{B = I - {D^{-1}}}=-D^{-1}(L+U),f=D^{-1}b
B=I−D−1=−D−1(L+U),f=D−1b 我们称
φ
(
x
)
=
B
x
+
f
\varphi (x) = Bx + f
φ(x)=Bx+f为迭代函数，任取初始向量
x
=
x
(
0
)
x = {x^{(0)}}
x=x(0)，按照
x
(
k
+
1
)
=
B
x
(
k
)
+
f
{x^{(k + 1)}} = B{x^{(k)}} + f
x(k+1)=Bx(k)+f形成的迭代格式，称这种迭代格式为雅克比迭代法。
C语言程序
#include <stdio.h>
#include <stdlib.h>
#include<math.h>
#define MaxSize 100
double A
[MaxSize
][MaxSize
];
double B
[MaxSize
];
double C
[MaxSize
][MaxSize
];
double D
[MaxSize
][MaxSize
];
double E
[MaxSize
][MaxSize
];
double F
[MaxSize
];
double X
[MaxSize
];
double X1
[MaxSize
];
double Y
[MaxSize
];
#define e 1e-6
int n
;
void InitMatrix()
{
int i
,j
;
for(i
=0;i
<n
;i
++)
for(j
=0;j
<n
;j
++)
{
if(i
==j
)
{
D
[i
][j
]=1/A
[i
][i
];
E
[i
][j
]=0;
}
if(i
<j
)
E
[i
][j
]=A
[i
][j
];
if(i
>j
)
E
[i
][j
]=A
[i
][j
];
}
}
void Jacobi()
{
int i
,j
,k
,r
;
double sum1
,sum2
,sum4
=0;
for(i
=0;i
<n
;i
++)
for(j
=0;j
<n
;j
++)
{
sum1
=0;
for(k
=0;k
<n
;k
++)
{
sum1
=sum1
+D
[i
][k
]*E
[k
][j
];
}
E
[i
][j
]=-sum1
;
}
for(i
=0;i
<n
;i
++)
for(j
=0;j
<n
;j
++)
sum4
=sum4
+pow(E
[i
][j
],2);
if(sqrt(sum4
)<1)
printf("/*****B矩阵F范数<1，必有数值解！*****/\n");
else
printf("/*****B矩阵F范数>1，不一定有数值解！*****/\n");
for(i
=0;i
<n
;i
++)
{
sum2
=0;
for(k
=0;k
<n
;k
++)
{
sum2
=sum2
+D
[i
][k
]*B
[k
];
}
F
[i
]=sum2
;
}
for(r
=1;r
<100;r
++)
{
int flag
=0;
double sum3
;
for(i
=0;i
<n
;i
++)
X1
[i
]=X
[i
];
for(i
=0;i
<n
;i
++)
{
sum3
=0;
for(k
=0;k
<n
;k
++)
{
sum3
=sum3
+E
[i
][k
]*X1
[k
];
}
X
[i
]=F
[i
]+sum3
;
}
for(j
=0;j
<n
;j
++)
if(fabs(X
[j
]-X1
[j
])<e
)
flag
++;
if(flag
==n
)
{
printf("迭代第%d次满足精度！\n",r
);
break;
}
printf("第%d次迭代向量X：\n",r
);
for(i
=0;i
<n
;i
++)
printf("%lf\n",X
[i
]);
}
}
void input()
{
int i
,j
;
printf("请输入系数矩阵A：\n");
for(i
=0;i
<n
;i
++)
for(j
=0;j
<n
;j
++)
scanf("%lf",&A
[i
][j
]);
printf("请输入向量B：\n");
for(i
=0;i
<n
;i
++)
scanf("%lf",&B
[i
]);
printf("请输入初始向量X：\n");
for(i
=0;i
<n
;i
++)
scanf("%lf",&X
[i
]);
}
void print()
{
int i
,j
;
printf("方程组的近似解：\n");
for(i
=0;i
<n
;i
++)
printf("%lf\n",X
[i
]);
}
int main()
{
void InitMatrix();
void Jacobi();
void input();
void print();
printf("请输入矩阵行数：\n");
scanf("%d",&n
);
printf("\n");
input();
InitMatrix();
Jacobi();
print();
return 0;
}
实验结果
{
x
1
+
2
x
2
−
2
x
3
=
1
x
1
+
x
2
+
x
3
=
3
2
x
1
+
2
x
2
+
x
3
=
5
\begin{cases} {{x_1} + 2{x_2} - 2{x_3} = 1}\\ {{x_1} + {x_2} + {x_3} = 3}\\ {2{x_1} + 2{x_2} + {x_3} = 5} \end{cases}
⎩⎪⎨⎪⎧x1+2x2−2x3=1x1+x2+x3=32x1+2x2+x3=5 初始向量为：
x
(
0
)
=
(
0
,
0
,
0
)
T
{x^{(0)}} = {(0,0,0)^T}
x(0)=(0,0,0)T
