C语言求逆矩阵的实现方法

发布时间:2021-08-16 15:16:15 作者:chen
来源:亿速云 阅读:1265

这篇文章主要讲解了“C语言求逆矩阵的实现方法”,文中的讲解内容简单清晰,易于学习与理解,下面请大家跟着小编的思路慢慢深入,一起来研究和学习“C语言求逆矩阵的实现方法”吧!

一般求逆矩阵的方法有两种,伴随阵法和初等变换法。但是这两种方法都不太适合编程。伴随阵法的计算量大,初等变换法又难以编程实现。
适合编程的求逆矩阵的方法如下:

  1. 对可逆矩阵A进行QR分解:A=QR

  2. 求上三角矩阵R的逆矩阵

  3. 求出A的逆矩阵:A^(-1)=R^(-1)Q^(H)

以上三步都有具体的公式与之对应,适合编程实现。
C语言实现代码:

#include <stdio.h>
#include <math.h>

#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<SIZE;++i)
    {
        for(j=0;j<SIZE;++j)
        {
            b[j][i]=init[j][i];
        }
        for(k=0;k<i;++k)
        {
            if(i)
            {
                numerator=0.0;
                denominator=0.0;
                for(l=0;l<SIZE;++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<SIZE;++j)
                {
                    b[j][i]-=t[k][i]*b[j][k];//t  init  =0  !!!
                }
            }
        }
    }
    ///////////////////对B单位化,得到正交矩阵Q矩阵////////////////////
    for(i=0;i<SIZE;++i)
    {
        numb=0.0;
        for(j=0;j<SIZE;++j)
        {
            numb+=(b[j][i]*b[j][i]);
        }
        dev=sqrt(numb);
        for(j=0;j<SIZE;++j)
        {
            Q[j][i]=b[j][i]/dev;
        }
        matrixR1[i][i]=dev;
    }
    /////////////////////求上三角R阵///////////////////////
    for(i=0;i<SIZE;++i)
    {
        for(j=0;j<SIZE;++j)
        {
            if(j<i)
            {
                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<SIZE;++i)
    {
        for(j=0;j<SIZE;++j)
        {
            printf("%2.4f    ",Q[i][j]);
        //  
        }
        printf("\n");
    }
    printf("R=\n");
    for(i=0;i<SIZE;++i)
    {
        for(j=0;j<SIZE;++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<SIZE;++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<SIZE;++k)
                {
                    invR[j][k]-=ratio*invR[i][k];
                    R[j][k]-=ratio*R[i][k];
                }
            }   
        }
    }

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

    printf("inv(R)=\n");
    for(i=0;i<SIZE;++i)
    {
        for(j=0;j<SIZE;++j)
        {
            printf(" %2.4f  ",invR[i][j]);
        //  
        }
        printf("\n");
    }
////////////////////结果和MATLAB差一个负号,神马鬼????????/////////////////////
/////////////////////求QH//////////////////////////
    for(i=0;i<SIZE;++i)//实矩阵就是转置
    {
        for(j=0;j<SIZE;++j)
        {
            QH[i][j]=Q[j][i];
        }
    }
///////////////////////求A的逆阵invA/////////////////////////////

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

    printf("inv(A)=\n");
    for(i=0;i<SIZE;++i)
    {
        for(j=0;j<SIZE;++j)
        {
            printf(" %2.4f  ",invA[i][j]);
        //  
        }
        printf("\n");
    }

///////////////////////结果与MATLAB的结果在千分位后有出入,但是负号都是对的^v^///////////////////////////
    return 0;
}

另附上矩阵乘法的子函数

/*/
函数名:void mulMatrix(double matrix1[SIZE][SIZE],double matrix2[SIZE][SIZE],int high2,int weight,int weight2,double mulMatrixOut[SIZE][SIZE])
输入:依次是 左矩阵,右矩阵,左矩阵高度,左矩阵宽度,右矩阵宽度,输出矩阵
输出:
功能:矩阵乘法
作者:HLdongdong
*//
void mulMatrix(double matrix1[SIZE][SIZE],double matrix2[SIZE][SIZE],int high2,int weight,int weight2,double mulMatrixOut[SIZE][SIZE])
{
    int i,j,k;
    for(i=0;i<high2;++i)
    {
        for(j=0;j<weight2;j++)
        {
            for(k=0;k<weight;++k)
            {
                mulMatrixOut[i][j]+=matrix1[i][k]*matrix2[k][j];
            }
        }
    }
}

感谢各位的阅读,以上就是“C语言求逆矩阵的实现方法”的内容了,经过本文的学习后,相信大家对C语言求逆矩阵的实现方法这一问题有了更深刻的体会,具体使用情况还需要大家实践验证。这里是亿速云,小编将为大家推送更多相关知识点的文章,欢迎关注!

推荐阅读:
  1. C语言中怎么求大数的阶乘
  2. c语言求平方的函数怎么写?

免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。

c语言

上一篇:Java怎么实现Spring的IoC容器的操作

下一篇:Spring中如何实现给Bean属性赋值

相关阅读

您好,登录后才能下订单哦!

密码登录
登录注册
其他方式登录
点击 登录注册 即表示同意《亿速云用户服务条款》