求逆矩阵的方法「建议收藏」

一般求逆矩阵的方法有两种，伴随阵法和初等变换法。但是这两种方法都不太适合编程。伴随阵法的计算量大，初等变换法又难以编程实现。
适合编程的求逆矩阵的方法如下：
1、对可逆矩阵A进行QR分解：A=QR
2、求上三角矩阵R的逆矩阵
3、求出A的逆矩阵：A^(-1)=R^(-1)Q^(H)
以上三步都有具体的公式与之对应，适合编程实现。
C语言实现代码：

#include 

#include 



#define SIZE 8



double b[SIZE][SIZE]={

  

  0};//应该读作“贝尔塔”，注释中用B表示

double t[SIZE][SIZE]={

  

  0};//求和的那项

double Q[SIZE][SIZE]={

  

  0};//正交矩阵

double QH[SIZE][SIZE]={

  

  0};//正交矩阵的转置共轭

double R[SIZE][SIZE]={

  

  0};//

double invR[SIZE][SIZE]={

  

  0};//R的逆矩阵

double invA[SIZE][SIZE]={

  

  0};//A的逆矩阵，最终的结果

//={0};//

double matrixR1[SIZE][SIZE]={

  

  0};

double matrixR2[SIZE][SIZE]={

  

  0};



//double init[3][3]={3,14,9,6,43,3,6,22,15};

double init[8][8]={  

    0.0938  ,  0.5201 ,   0.4424  ,  0.0196  ,  0.3912  ,  0.9493 ,   0.9899  ,  0.8256,

    0.5254  ,  0.3477 ,   0.6878  ,  0.3309 ,   0.7691  ,  0.3276 ,   0.5144  ,  0.7900,

    0.5303  ,  0.1500 ,   0.3592  ,  0.4243 ,   0.3968  ,  0.6713 ,   0.8843  ,  0.3185,

    0.8611  ,  0.5861 ,   0.7363  ,  0.2703 ,   0.8085  ,  0.4386 ,   0.5880  ,  0.5341,

    0.4849  ,  0.2621 ,   0.3947  ,  0.1971 ,   0.7551  ,  0.8335 ,   0.1548  ,  0.0900,

    0.3935  ,  0.0445 ,   0.6834  ,  0.8217 ,   0.3774  ,  0.7689 ,   0.1999  ,  0.1117,

    0.6714  ,  0.7549 ,   0.7040  ,  0.4299 ,   0.2160  ,  0.1673 ,   0.4070  ,  0.1363,

    0.7413  ,  0.2428 ,   0.4423  ,  0.8878 ,   0.7904  ,  0.8620 ,   0.7487  ,  0.6787

};

/*/ 函数名:int main() 输入: 输出: 功能:求矩阵的逆 pure C language 首先对矩阵进行QR分解之后求上三角矩阵R的逆阵最后A-1=QH*R-1,得到A的逆阵。 作者：HLdongdong *//////////////////////////////////////////////////////////////////////

int main()

{

    int i;//数组 行

    int j;//数组 列

    int k;//代表B的角标

    int l;//数组 列

    double dev;

    double numb;//计算的中间变量

    double numerator,denominator;

    double ratio;

    /////////////////求B/////////////////

    for(i=0;i
    {

        for(j=0;j
        {

            b[j][i]=init[j][i];

        }

        for(k=0;k
        {

            if(i)

            {

                numerator=0.0;

                denominator=0.0;

                for(l=0;l
                {

                    numerator+=init[l][i]*b[l][k];

                    denominator+=b[l][k]*b[l][k];

                }

                dev=numerator/denominator;

                t[k][i]=dev;

                for(j=0;j
                {

                    b[j][i]-=t[k][i]*b[j][k];//t init =0 !!!

                }

            }

        }

    }

    ///////////////////对B单位化,得到正交矩阵Q矩阵////////////////////

    for(i=0;i
    {

        numb=0.0;

        for(j=0;j
        {

            numb+=(b[j][i]*b[j][i]);

        }

        dev=sqrt(numb);

        for(j=0;j
        {

            Q[j][i]=b[j][i]/dev;

        }

        matrixR1[i][i]=dev;

    }

    /////////////////////求上三角R阵///////////////////////

    for(i=0;i
    {

        for(j=0;j
        {

            if(j
            {

                matrixR2[j][i]=t[j][i];

            }

            else if(j==i)   

            {

                matrixR2[j][i]=1;

            }

            else

            {

                matrixR2[j][i]=0;

            }

        }

    }

    mulMatrix(matrixR1,matrixR2,SIZE,SIZE,SIZE,R);

///////////////////////QR分解完毕//////////////////////////

    printf("QR分解：\n");

    printf("Q=\n");

    for(i=0;i
    {

        for(j=0;j
        {

            printf("%2.4f ",Q[i][j]);

        // 

        }

        printf("\n");

    }

    printf("R=\n");

    for(i=0;i
    {

        for(j=0;j
        {

            printf("%2.4f ",R[i][j]);

        // 

        }

        printf("\n");

    }

/////////////////////求R的逆阵//////////////////////////

    for(i=SIZE-1;i>=0;--i)

    {

        invR[i][i]=1/R[i][i];

        //R[i][i]=invR[i][i];

        if(i!=(SIZE-1))//向右

        {

            for(j=i+1;j
            {

                invR[i][j]=invR[i][j]*invR[i][i];

                R[i][j]=R[i][j]*invR[i][i];

            }

        }

        if(i)//向上

        {

            for(j=i-1;j>=0;--j)

            {

                ratio=R[j][i];

                for(k=i;k
                {

                    invR[j][k]-=ratio*invR[i][k];

                    R[j][k]-=ratio*R[i][k];

                }

            }   

        }

    }



///////////////////////////////////////////////////////



    printf("inv（R）=\n");

    for(i=0;i
    {

        for(j=0;j
        {

            printf(" %2.4f ",invR[i][j]);

        // 

        }

        printf("\n");

    }

////////////////////结果和MATLAB差一个负号，神马鬼？？？？？？？？/////////////////////

/////////////////////求QH//////////////////////////

    for(i=0;i
    {

        for(j=0;j
        {

            QH[i][j]=Q[j][i];

        }

    }

///////////////////////求A的逆阵invA/////////////////////////////



    mulMatrix(invR,QH,SIZE,SIZE,SIZE,invA);



    printf("inv（A）=\n");

    for(i=0;i
    {

        for(j=0;j
        {

            printf(" %2.4f ",invA[i][j]);

        // 

        }

        printf("\n");

    }



///////////////////////结果与MATLAB的结果在千分位后有出入，但是负号都是对的^v^///////////////////////////

    return 0;

}

另附上矩阵乘法的子函数

/*/

函数名:void mulMatrix(double matrix1[SIZE][SIZE],double matrix2[SIZE][SIZE],int high1,int weight,int weight2,double mulMatrixOut[SIZE][SIZE])

输入:依次是 左矩阵，右矩阵，左矩阵高度，左矩阵宽度，右矩阵宽度，输出矩阵

输出:

功能:矩阵乘法

作者：HLdongdong

*//

void mulMatrix(double matrix1[SIZE][SIZE],double matrix2[SIZE][SIZE],int high1,int weight,int weight2,double mulMatrixOut[SIZE][SIZE])

{

 int i,j,k;

 for(i=0;i
 { 

   

 for(j=0;j
 { 

   

 for(k=0;k
 { 

   

 mulMatrixOut[i][j]+=matrix1[i][k]*matrix2[k][j];

 }

 }

 }

}

求逆矩阵的方法「建议收藏」

相关推荐