C语言求逆矩阵案例详解

时间:2021-09-22 06:00:55

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

  1. 对可逆矩阵A进行QR分解:A=QR
  2. 求上三角矩阵R的逆矩阵
  3. 求出A的逆矩阵:A^(-1)=R^(-1)Q^(H)

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

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
#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;
}

另附上矩阵乘法的子函数

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/*/
函数名: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<high1;++i)
    {
        for(j=0;j<weight2;j++)
        {
            for(k=0;k<weight;++k)
            {
                mulMatrixOut[i][j]+=matrix1[i][k]*matrix2[k][j];
            }
        }
    }
}

到此这篇关于C语言求逆矩阵案例详解的文章就介绍到这了,更多相关C语言求逆矩阵内容请搜索服务器之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持服务器之家!

原文链接:https://blog.csdn.net/dongdong_csdn/article/details/75045003