[数值算法]线性方程组求解算法---基于LU分解法的追赶法

[数值算法]线性方程组求解算法---基于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