关于复数矩阵求逆,如果使用 MATLAB,就非常简单。我们先用一个 MATLAB 的例子来说明,等会要将 C 语言的程序和 MATLAB 的程序进行对比。
close all;
clear all;
clc;
%定义矩阵a为复数矩阵
a = [[4+2*i,3+1*i,4+3*i,5+5*i];
[1+7*i,8+2*i,2+2*i,9+3*i];
[4+4*i,5+6*i,1+7*i,7+2*i];
[6+1*i,7+8*i,1+4*i,2+5*i]];
a_real = real(a); %求矩阵实部
a_imag = imag(a); %求矩阵虚部
a_inv = inv(a); %求矩阵的逆
可以看到 a_inv 是这样的:
下面列出 C 语言的矩阵求逆的代码。由于VS2012中不能使用 complex.h,我们在这里将复数矩阵的实部和虚部分开定义,数值和上面定义的 a 一样。
a_real[0][0] = 4; a_real[0][1] = 3; a_real[0][2] = 4; a_real[0][3] = 5;
a_real[1][0] = 1; a_real[1][1] = 8; a_real[1][2] = 2; a_real[1][3] = 9;
a_real[2][0] = 4; a_real[2][1] = 5; a_real[2][2] = 1; a_real[2][3] = 7;
a_real[3][0] = 6; a_real[3][1] = 7; a_real[3][2] = 1; a_real[3][3] = 2;
a_imag[0][0] = 2; a_imag[0][1] = 1; a_imag[0][2] = 3; a_imag[0][3] = 5;
a_imag[1][0] = 7; a_imag[1][1] = 2; a_imag[1][2] = 2; a_imag[1][3] = 3;
a_imag[2][0] = 4; a_imag[2][1] = 6; a_imag[2][2] = 7; a_imag[2][3] = 2;
a_imag[3][0] = 1; a_imag[3][1] = 8; a_imag[3][2] = 4; a_imag[3][3] = 5;
完整代码及其注释如下:
#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#define N 4
#define DEBUG 1 //debug label,0即不打印相关结果,非0打印相关输出结果
//计算实数矩阵的逆
void matrix_inverse_LU(double a[][N],double a_inverse[N][N])
{
double l[N][N], u[N][N];
double l_inverse[N][N], u_inverse[N][N];
//double a_inverse[N][N];
int i, j, k;
double s, t;
memset(l, 0, sizeof(l));
memset(u, 0, sizeof(u));
memset(l_inverse, 0, sizeof(l_inverse));
memset(u_inverse, 0, sizeof(u_inverse));
memset(a_inverse, 0, sizeof(u_inverse));
for (i = 0; i < N;i++) //计算l矩阵对角线
{
l[i][i] = 1;
}
for (i = 0;i < N;i++)
{
for (j = i;j < N;j++)
{
s = 0;
for (k = 0;k < i;k++)
{
s += l[i][k] * u[k][j];
}
u[i][j] = a[i][j] - s; //按行计算u值
}
for (j = i + 1;j < N;j++)
{
s = 0;
for (k = 0; k < i; k++)
{
s += l[j][k] * u[k][i];
}
l[j][i] = (a[j][i] - s) / u[i][i]; //按列计算l值
}
}
for (i = 0;i < N;i++) //按行序,行内从高到低,计算l的逆矩阵
{
l_inverse[i][i] = 1;
}
for (i= 1;i < N;i++)
{
for (j = 0;j < i;j++)
{
s = 0;
for (k = 0;k < i;k++)
{
s += l[i][k] * l_inverse[k][j];
}
l_inverse[i][j] = -s;
}
}
#if DEBUG
printf("test l_inverse:\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
s = 0;
for (k = 0; k < N; k++)
{
s += l[i][k] * l_inverse[k][j];
}
printf("%f ", s);
}
putchar('\n');
}
#endif
for (i = 0;i < N;i++) //按列序,列内按照从下到上,计算u的逆矩阵
{
u_inverse[i][i] = 1 / u[i][i];
}
for (i = 1;i < N;i++)
{
for (j = i - 1;j >=0;j--)
{
s = 0;
for (k = j + 1;k <= i;k++)
{
s += u[j][k] * u_inverse[k][i];
}
u_inverse[j][i] = -s / u[j][j];
}
}
#if DEBUG
printf("test u_inverse:\n");
for (i = 0;i < N;i++)
{
for (j = 0;j < N;j++)
{
s = 0;
for (k = 0;k < N;k++)
{
s += u[i][k] * u_inverse[k][j];
}
printf("%f ",s);
}
putchar('\n');
}
#endif
for (i = 0;i < N;i++) //计算矩阵a的逆矩阵
{
for (j = 0;j < N;j++)
{
for (k = 0;k < N;k++)
{
a_inverse[i][j] += u_inverse[i][k] * l_inverse[k][j];
}
}
}
#if DEBUG
printf("test a:\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
s = 0;
for (k = 0; k < N; k++)
{
s += a[i][k] * a_inverse[k][j];
}
printf("%f ", s);
}
putchar('\n');
}
#endif
}
//矩阵乘法,由于这里计算的是长度N的方阵,所以实际上ROW,MID,COL的值都是N
void MulMatrix(const double (*arr1)[N], const double (*arr2)[N], int ROW, int MID, int COL, double (*arr3)[N])
{
double sum=0.0;
int i,j,m;
for (i = 0; i<ROW; i++)
{
for (j = 0; j<COL; j++)
{
sum = 0.0;
for (m = 0; m<MID; m++)
{
sum = sum + arr1[i][m] * arr2[m][j];
}
arr3[i][j] = sum;
}
}
}
//矩阵加法
void PlusMatrix(const double (*arr1)[N], const double (*arr2)[N], int ROW, int COL, double (*arr3)[N])
{
int i,j;
for(i=0;i<N;i++)//控制行
{
for(j=0;j<N;j++)
{
arr3[i][j]=arr1[i][j]+arr2[i][j];
}
}
}
//矩阵减法
void MinusMatrix(const double (*arr1)[N], const double (*arr2)[N], int ROW, int COL, double (*arr3)[N])
{
int i,j;
for(i=0;i<N;i++)//控制行
{
for(j=0;j<N;j++)
{
arr3[i][j]=arr1[i][j]-arr2[i][j];
}
}
}
void main()
{
int i, j, k;
double a[N][N]; //N表示矩阵维度,为4
double a_real[N][N];
double a_imag[N][N];
double Ainv[N][N],Binv[N][N],BAinv[N][N],BAinvB[N][N],A_P_BAinvB[N][N],A_P_BAinvB_inv[N][N],AinvB[N][N],AinvB_A_P_BAinvB[N][N],AinvB_A_P_BAinvB_inv[N][N];
//将a矩阵的实部和虚部分开定义
a_real[0][0] = 4;a_real[0][1] = 3;a_real[0][2] = 4;a_real[0][3] = 5;
a_real[1][0] = 1;a_real[1][1] = 8;a_real[1][2] = 2;a_real[1][3] = 9;
a_real[2][0] = 4;a_real[2][1] = 5;a_real[2][2] = 1;a_real[2][3] = 7;
a_real[3][0] = 6;a_real[3][1] = 7;a_real[3][2] = 1;a_real[3][3] = 2;
a_imag[0][0] = 2;a_imag[0][1] = 1;a_imag[0][2] = 3;a_imag[0][3] = 5;
a_imag[1][0] = 7;a_imag[1][1] = 2;a_imag[1][2] = 2;a_imag[1][3] = 3;
a_imag[2][0] = 4;a_imag[2][1] = 6;a_imag[2][2] = 7;a_imag[2][3] = 2;
a_imag[3][0] = 1;a_imag[3][1] = 8;a_imag[3][2] = 4;a_imag[3][3] = 5;
//这些计算公式来源于这里:https://wenku.baidu.com/view/2de4c1bc284ac850ad024244.html
matrix_inverse_LU(a_real,Ainv);
matrix_inverse_LU(a_imag,Binv);
MulMatrix(a_imag,Ainv,N,N,N,BAinv);
MulMatrix(BAinv,a_imag,N,N,N,BAinvB);
PlusMatrix(a_real,BAinvB,N,N,A_P_BAinvB);
matrix_inverse_LU(A_P_BAinvB,A_P_BAinvB_inv);
MulMatrix(Ainv,a_imag,N,N,N,AinvB);
MulMatrix(AinvB,A_P_BAinvB_inv,N,N,N,AinvB_A_P_BAinvB_inv);
//最后一步要把AinvB_A_P_BAinvB_inv每个数都求相反数,也是上面链接的文献里说的
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
AinvB_A_P_BAinvB_inv[i][j] = -AinvB_A_P_BAinvB_inv[i][j];
}
}
//输出a的逆矩阵的实部矩阵
printf("test a_inv_real:\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
printf("%f ", A_P_BAinvB_inv[i][j]);
}
printf("\n");
}
//输出a的逆矩阵的虚部矩阵
printf("test a_inv_imag:\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
{
printf("%f ", AinvB_A_P_BAinvB_inv[i][j]);
}
printf("\n");
}
//使得窗口不要关闭
getchar();
}
最后输出结果如下:
通过对比,可以发现该结果中的实部和虚部与 MATLAB 输出的结果是一致的。这证明我们的程序是正确的。