[数值算法]线性方程组求解算法---基于LU分解法的追赶法
By EmilMatthew
05/9/08
基于LU分解法的追赶法主要用于求解一类形式上特殊的线性方程组。
Ax=b,其中,当|i-j|>1时，有aij=0.求解这类形式的线性方程组，在样条函数，常微分方程边值问题的数值解中都会有用.
原理上仍是考虑将这个矩阵分成L和U两个矩阵，不过，由于这个矩阵存大着大量的0值元素，所以，给最终的递推求解带来了相当大的便利.
更加详细的内容请去参考相应的数值算法书籍.
/*This function :chasingMethod is coded by EmilMatthew,you can use and modify it as you wish,
But I do not pomise it can fit all your need .*/
void chasingMethod(Type** inMatrixArr,Type* bList,Type* xAnsList,int size)
{
/*Core think of this algorithm,this algorithm is just a special version of the LUPationMethod:
matrix_L*yAnsList=bList;
matrix_U*xAnsList=yAnsList;
*/
Type* listA,* listB,* listC;/*to crack the in matrix in three one demision list*/
Type* listAlpha,* listBeta;/*The tween list for resolve the problem*/
Type* yAnsList;
int i,j;
/*pointer data assertion*/
assertF(inMatrixArr!=NULL,"in LUPationMethod,matrixArr is NULL\n");
assertF(bList!=NULL,"in LUPationMethod,bList is NULL\n");
assertF(xAnsList!=NULL,"in LUPationMethod,xAnsList is NULL\n");
/*The Assertion of the pass in matrix's format.*/
for(i=0;i<size;i++)
for(j=i;j<size;j++)
if(j-i>1)
assertF(inMatrixArr[i][j]==0,"In chasingMethod,the pass in matrix format is not correct\n");
/*Mem Apply*/
listArrMemApply(&listA,size);
listArrMemApply(&listB,size);
listArrMemApply(&listC,size);
listArrMemApply(&listAlpha,size);
listArrMemApply(&listBeta,size);
listArrMemApply(&yAnsList,size);
/*locate the inMatrixArr data into correct position*/
listB[0]=inMatrixArr[0][0];
listC[0]=inMatrixArr[0][1];
for(i=1;i<size-1;i++)
{
listA[i]=inMatrixArr[i][i-1];
listB[i]=inMatrixArr[i][i];
listC[i]=inMatrixArr[i][i+1];
}
listA[size-1]=inMatrixArr[size-1][size-2];
listB[size-1]=inMatrixArr[size-1][size-1];
/*main program of the chasingMethod*/
listAlpha[0]=listB[0];
listBeta[0]=listC[0]/listAlpha[0];
for(i=1;i<size-1;i++)
{
listAlpha[i]=listB[i]-listA[i]*listBeta[i-1];
listBeta[i]=listC[i]/listAlpha[i];
}
listAlpha[i]=listB[i]-listA[i]*listBeta[i-1];
/*To get the yAnsList*/
yAnsList[0]=bList[0]/listAlpha[0];
for(i=1;i<size;i++)
yAnsList[i]=(bList[i]-listA[i]*yAnsList[i-1])/listAlpha[i];
showArrListFloat(yAnsList,0,size);
/*To get the xAnsList*/
xAnsList[size-1]=yAnsList[size-1];
for(i=size-2;i>=0;i--)
xAnsList[i]=yAnsList[i]-listBeta[i]*xAnsList[i+1];
/*Mem Free*/
free(listA);
free(listB);
free(listC);
free(listAlpha);
free(listBeta);
free(yAnsList);
}
/*算法简单测试*/
//inputData
3;
2,1,0,3;
1,3,1,5;
0,1,2,3;
//ouputData:
after the chasingMethod:the ans x rows is:(from x0 to xn-1)
1.00000 1.00000 1.00000