[数值算法]线性方程组的求解---迭代法小结.
by EmilMatthew 05/9/08
当迭代法做多了之后，便感受到这种思想的好用了,今天再次实现用迭代法求解问题，感觉轻松了不少。
当某个关于x的迭代式xk=fi(xk_1)保持收敛时，则可利迭代来求出最后的结果，直到满足精度要求为止.
对于线性方程组的求解中迭代法的应用亦是如此,当矩阵的规模大于100*100时，迭代法的优势就突现出来了.
1) Jacobi
对于线性方程组，可以很方便的对第i行中的元素xi用该行中的其余x来表示，
即:
Xi_k+1=1/aii*(bi-sum(aij*Xj_k&&j1=i))
在数学上可以证明它满足迭代条件,所以，这便有了第一个求线性方程组的迭代算法:
Jacobi算法
/*
JacobiMethod, coded by EmilMathew 05/9/8, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.
*/
/*main aid function,to determine the end of the iterator*/
int preciseReached(Type* xList1,Type* xList2,int len,Type precise)
{
int i=0;
int flag;
assertF(xList1!=NULL,"in preciseReached,xList1 is NULL\n");
assertF(xList2!=NULL,"in preciseReached,xList2 is NULL\n");
flag=1;
for(i=0;i<len;i++)
if(fabs(xList1[i]-xList2[i])>precise)
{
flag=0;
break;
}
return flag;
}
void JacobiMethod(Type** inMatrixArr,Type* bList,Type* xAnsList,Type precise,int size,
FILE* outputFile)
{
int iteratorNum;
int i=0,j=0;
Type sum=0;
Type* tmpXList;
/*pointer data assertion*/
assertF(inMatrixArr!=NULL,"in JacobiMethod,matrixArr is NULL\n");
assertF(bList!=NULL,"in JacobiMethod,bList is NULL\n");
assertF(xAnsList!=NULL,"in JacobiMethod,xAnsList is NULL\n");
tmpXList=(Type*)malloc(sizeof(Type)*size);
assertF(tmpXList!=NULL,"in JacobiMethod,tmpXList is NULL\n");
/*init data*/
iteratorNum=0;
for(i=0;i<size;i++)
xAnsList[i]=0;
fprintf(outputFile,"k\t\t");
for(i=1;i<=size;i++)
fprintf(outputFile,"x%d\t\t",i);
fprintf(outputFile,"\r\n");
fprintf(outputFile,"%-15d",iteratorNum);
for(i=0;i<size;i++)
fprintf(outputFile,"%-15f",xAnsList[i]);
fprintf(outputFile,"\r\n");
iteratorNum++;
do
{
listArrCopy(xAnsList,tmpXList,size);
for(i=0;i<size;i++)
{
sum=0;
for(j=0;j<size;j++)
if(i!=j)
sum+=inMatrixArr[i][j]*tmpXList[j];
xAnsList[i]=(float)1/inMatrixArr[i][i]*(bList[i]-sum);
}
fprintf(outputFile,"%-15d",iteratorNum);
for(i=0;i<size;i++)
fprintf(outputFile,"%-15f",xAnsList[i]);
fprintf(outputFile,"\r\n");
iteratorNum++;
}
while(!preciseReached(xAnsList,tmpXList,size,precise));
free(tmpXList);
}
2 Gauss_Seidel
接下来是一个改进版本，主要提高了迭代速度.注意到在逐个求x_k=1的分量时，当计算到分量xi_k+1时，分量x1_k+1~xi-1_k+1都已得到，而用新分量求下一步分量会比原来更准确，速度上也会有提高,而这个动作在编程上更易实现，原来在迭代的过程中还要用到一个tmpXList数组，此时已经不需要在主迭代过程中使用了.
这便有了下面的Gauss_Seidel方法:
/*
Gauss_Seidel, coded by EmilMathew 05/9/8, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.
*/
void Gauss_SeidelMethod(Type** inMatrixArr,Type* bList,Type* xAnsList,Type precise,int size,FILE* outputFile)
{
int iteratorNum;
int i=0,j=0;
Type sum=0;
Type* tmpXList;
/*pointer data assertion*/
assertF(inMatrixArr!=NULL,"in JacobiMethod,matrixArr is NULL\n");
assertF(bList!=NULL,"in JacobiMethod,bList is NULL\n");
assertF(xAnsList!=NULL,"in JacobiMethod,xAnsList is NULL\n");
tmpXList=(Type*)malloc(sizeof(Type)*size);
assertF(tmpXList!=NULL,"in JacobiMethod,tmpXList is NULL\n");
/*init data*/
iteratorNum=0;
for(i=0;i<size;i++)
xAnsList[i]=0;
fprintf(outputFile,"k\t\t");
for(i=1;i<=size;i++)
fprintf(outputFile,"x%d\t\t",i);
fprintf(outputFile,"\r\n");
fprintf(outputFile,"%-15d",iteratorNum);
for(i=0;i<size;i++)
fprintf(outputFile,"%-15f",xAnsList[i]);
fprintf(outputFile,"\r\n");
iteratorNum++;
do
{
listArrCopy(xAnsList,tmpXList,size);
for(i=0;i<size;i++)
{
sum=0;
for(j=0;j<size;j++)
if(i!=j)
sum+=inMatrixArr[i][j]*xAnsList[j];
xAnsList[i]=(float)1/inMatrixArr[i][i]*(bList[i]-sum);
}
fprintf(outputFile,"%-15d",iteratorNum);
for(i=0;i<size;i++)
fprintf(outputFile,"%-15f",xAnsList[i]);
fprintf(outputFile,"\r\n");
iteratorNum++;
}
while(!preciseReached(xAnsList,tmpXList,size,precise));
free(tmpXList);
}
3逐次超松弛迭代法:
逐次超松弛迭代法主要是在前面的Gauss_Seidel基础上加以改进，考虑引入适当的松驰因子及做适当的线性变换可得.
/*
SORMethod, coded by EmilMathew 05/9/8, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.
*/
void SORMethod(Type** inMatrixArr,Type* bList,Type* xAnsList,Type w,Type precise,int size,FILE* outputFile)
{
int iteratorNum;
int i=0,j=0;
Type sum=0;
Type* tmpXList;
/*pointer data assertion*/
assertF(inMatrixArr!=NULL,"in JacobiMethod,matrixArr is NULL\n");
assertF(bList!=NULL,"in JacobiMethod,bList is NULL\n");
assertF(xAnsList!=NULL,"in JacobiMethod,xAnsList is NULL\n");
tmpXList=(Type*)malloc(sizeof(Type)*size);
assertF(tmpXList!=NULL,"in JacobiMethod,tmpXList is NULL\n");
/*init data*/
iteratorNum=0;
for(i=0;i<size;i++)
xAnsList[i]=0;
fprintf(outputFile,"k\t\t");
for(i=1;i<=size;i++)
fprintf(outputFile,"x%d\t\t",i);
fprintf(outputFile,"\r\n");
fprintf(outputFile,"%-15d",iteratorNum);
for(i=0;i<size;i++)
fprintf(outputFile,"%-15f",xAnsList[i]);
fprintf(outputFile,"\r\n");
iteratorNum++;
do
{
listArrCopy(xAnsList,tmpXList,size);
for(i=0;i<size;i++)
{
sum=0;
for(j=0;j<size;j++)
if(i!=j)
sum+=inMatrixArr[i][j]*xAnsList[j];
xAnsList[i]=xAnsList[i]+w/inMatrixArr[i][i]*(bList[i]-sum);
}
fprintf(outputFile,"%-15d",iteratorNum);
for(i=0;i<size;i++)
fprintf(outputFile,"%-15f",xAnsList[i]);
fprintf(outputFile,"\r\n");
iteratorNum++;
}
while(!preciseReached(xAnsList,tmpXList,size,precise));
free(tmpXList);
}
测试：
用以上三种方法对同一个线性方程组进行求解，精度为0.0001,得到的迭代结果如下:
testProgram:
inputData:
3,0.0001;
10,-1,-2,7.2;
-1,10,-2,8.3;
-1,-1,5,4.2;
Jacobi:
k x1 x2 x3
0 0.000000 0.000000 0.000000
1 0.720000 0.830000 0.840000
2 0.971000 1.070000 1.150000
3 1.057000 1.157100 1.248200
4 1.085350 1.185340 1.282820
5 1.095098 1.195099 1.294138
6 1.098337 1.198337 1.298039
7 1.099442 1.199442 1.299335
8 1.099811 1.199811 1.299777
9 1.099936 1.199936 1.299924
10 1.099978 1.199979 1.299975
Guass_SeidelMethod
k x1 x2 x3
0 0.000000 0.000000 0.000000
1 0.720000 0.902000 1.164400
2 1.043080 1.167188 1.282054
3 1.093130 1.195724 1.297771
4 1.099126 1.199467 1.299719
5 1.099890 1.199933 1.299965
6 1.099986 1.199992 1.299996
SORMethod 松驰系数w为1.05
k x1 x2 x3
0 0.000000 0.000000 0.000000
1 0.756000 0.950880 1.240445
2 1.078536 1.197696 1.297986
3 1.100408 1.199735 1.300131
4 1.099979 1.200039 1.299997
5 1.100004 1.199998 1.300001
可以看出，各种方法在迭代的速度上是逐个增加的,与理论是相一致的.
