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];//储存D逆
double E[MaxSize][MaxSize];//储存-(D-1)*(L+U)
double F[MaxSize];//(D-1)*B
double X[MaxSize];
double X1[MaxSize];
double Y[MaxSize];
#define e 1e-6//10的负6次方
int n;//矩阵尺寸大小
//所有迭代格式X(k+1)=B*X(k)+F,B范数<1,必收敛,B范数>1,未必发散
//雅克比迭代只适用于A[i][i]非0
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;//(D-1)*(L+U)
}
//求E的F范数,判断是否收敛
//E的F范数>1,不一定发散;但是E的F范数<1,一定收敛
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;//(D-1)*B
}
//X迭代
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;//更新X[i]
}
for(j=0;j<n;j++)
if(fabs(X[j]-X1[j])<e)//abs整型求绝对值,fabs浮点型求绝对值
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