这篇文章主要讲解了“C语言求逆矩阵的实现方法”,文中的讲解内容简单清晰,易于学习与理解,下面请大家跟着小编的思路慢慢深入,一起来研究和学习“C语言求逆矩阵的实现方法”吧!
一般求逆矩阵的方法有两种,伴随阵法和初等变换法。但是这两种方法都不太适合编程。伴随阵法的计算量大,初等变换法又难以编程实现。
适合编程的求逆矩阵的方法如下:
对可逆矩阵A进行QR分解:A=QR
求上三角矩阵R的逆矩阵
求出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语言求逆矩阵的实现方法这一问题有了更深刻的体会,具体使用情况还需要大家实践验证。这里是亿速云,小编将为大家推送更多相关知识点的文章,欢迎关注!
免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。