C++实现复数矩阵求逆 matlab inv

最新推荐文章于 2024-03-12 15:11:21 发布
最新推荐文章于 2024-03-12 15:11:21 发布
阅读量6.8k 收藏 66
点赞数 10
分类专栏： c/c++ 文章标签： matlab c++ 复数矩阵 复数矩阵求逆 矩阵求逆
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/srg1011/article/details/106411604
版权

C++实现复数矩阵求逆 matlab inv

一、引言

之前偶尔一次有用到将matlab转为C++语言的需求，其中matlab有一个inv函数可以非常方便的求矩阵的逆，甚至是复数矩阵。而C++中没有类似的函数。在csdn上有一个matlab2c的库的博客（github地址第80行开始），但是只有实数矩阵求逆的代码，而我又在百度上搜到一篇文献中写到将复数矩阵转为两个实数矩阵然后在进行求解。然后我就将两个代码一结合，就实现了复数矩阵求逆。经验证与matlab的inv函数结果一致。

二、原理

2.1 实数矩阵求逆

这里直接给出matlab2c中求实数矩阵的代码,其中函数inv的两个参数：a代表实数矩阵，用double二维数组存储、num代表这个矩阵的行数/列数。

void swap(double* a, double* b);  //声明子程序
double** inv(double** a, int num)
{
	int* is, * js, i, j, k;
	int n = num;
	double temp, fmax;
	double** tp = new double* [num];
	for (int i = 0; i < num; i++) tp[i] = new double[num];
	for (i = 0; i < n; i++) {
		for (j = 0; j < n; j++) {
			tp[i][j] = a[i][j];
		}
	}
	is = new int[n];
	js = new int[n];
	for (k = 0; k < n; k++)
	{
		fmax = 0.0;
		for (i = k; i < n; i++) {
			for (j = k; j < n; j++)
			{
				temp = fabs(tp[i][j]);//找最大值
				if (temp > fmax)
				{
					fmax = temp;
					is[k] = i; js[k] = j;
				}
			}
		}

		if ((fmax + 1.0) == 1.0)
		{
			delete[] is;
			delete[] js;
			return NULL;
		}
		if ((i = is[k]) != k)
			for (j = 0; j < n; j++)
				swap(&tp[k][j], &tp[i][j]);//交换指针
		if ((j = js[k]) != k)
			for (i = 0; i < n; i++)
				swap(&tp[i][k], &tp[i][j]);  //交换指针
		tp[k][k] = 1.0 / tp[k][k];
		for (j = 0; j < n; j++)
			if (j != k)
				tp[k][j] *= tp[k][k];
		for (i = 0; i < n; i++)
			if (i != k)
				for (j = 0; j < n; j++)
					if (j != k)
						tp[i][j] = tp[i][j] - tp[i][k] * tp[k][j];
		for (i = 0; i < n; i++)
			if (i != k)
				tp[i][k] *= -tp[k][k];
	}
	for (k = n - 1; k >= 0; k--)
	{
		if ((j = js[k]) != k)
			for (i = 0; i < n; i++)
				swap(&tp[j][i], &tp[k][i]);
		if ((i = is[k]) != k)
			for (j = 0; j < n; j++)
				swap(&tp[j][i], &tp[j][k]);
	}
	delete[] is;
	delete[] js;
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
		}
	}
	return tp;
}
void swap(double* a, double* b)
{
	double c;
	c = *a;
	*a = *b;
	*b = c;
}

2.2 复数矩阵求逆

在那篇文献中，我们可以看到将复数矩阵C拆分为实部矩阵A+虚部矩阵Bi，则
C − 1 = ( A + i B ) − 1 = ( A + B A − 1 B ) − 1 − i A − 1 B ( A + B A − 1 B ) − 1 C^{-1}=(A+iB)^{-1}=(A+BA^{-1}B)^{-1}-iA^{-1}B(A+BA^{-1}B)^{-1} C−1=(A+iB)−1=(A+BA−1B)−1−iA−1B(A+BA−1B)−1
这样就将复数矩阵求逆的过程转为求实数矩阵的逆、加、乘运算了。
代码：

complex<double>** GetMatrixInverse(complex<double>** src, int n) {
	double** A = new double* [n];
	double** B = new double* [n];
	for (int i = 0; i < n; i++) {
		A[i] = new double[n];
		B[i] = new double[n];
		for (int j = 0; j < n; j++) {
			A[i][j] = src[i][j].real();
			B[i][j] = src[i][j].imag();
		}
	}
	double** A1 = inv(A, n);
	double** A1B = Mul(A1, B, n, n, n);
	double** BA1B = Mul(B, A1B, n, n, n);
	double** AjBA1B = Add(A, BA1B, n, n);
	double** AjBA1B_1 = inv(AjBA1B, n);
	double** A1B_AjBA1B_1 = Mul(A1B, AjBA1B_1, n, n, n);

	complex<double>** res = new complex<double> * [n];
	for (int i = 0; i < n; i++) {
		res[i] = new complex<double>[n];
		for (int j = 0; j < n; j++) {
			res[i][j].real(AjBA1B_1[i][j]);
			res[i][j].imag(-1.0 * A1B_AjBA1B_1[i][j]);
		}
	}
	return res;
}

三、代码

完整代码如下：（包含main函数中测试代码）

#include <bits/stdc++.h>
using namespace std;

// 矩阵加法 a+b,其中a、b均为n*m型矩阵
double** Add(double** a, double** b, int n, int m) {
	double** res = new double* [n];
	for (int i = 0; i < n; i++) res[i] = new double[m];
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < m; j++) {
			res[i][j] = a[i][j] + b[i][j];
		}
	}
	return res;
}

// 矩阵乘法，a*b,a为n*m型矩阵，b为m*o型矩阵
double** Mul(double** a, double** b, int n, int m, int o) {
	double** res = new double* [n];
	double temp = 0.0;
	for (int i = 0; i < n; i++) res[i] = new double[o];
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < o; j++) {
			temp = 0.0;
			for (int k = 0; k < m; k++) {
				temp += a[i][k] * b[k][j];
			}
			res[i][j] = temp;
		}
	}
	return res;
}

void swap(double* a, double* b);  //声明子程序
// 实数矩阵求逆，返回a的逆，其中a为num型方阵
double** inv(double** a, int num)
{
	int* is, * js, i, j, k;
	int n = num;
	double temp, fmax;
	double** tp = new double* [num];
	for (int i = 0; i < num; i++) tp[i] = new double[num];
	for (i = 0; i < n; i++) {
		for (j = 0; j < n; j++) {
			tp[i][j] = a[i][j];
		}
	}
	is = new int[n];
	js = new int[n];
	for (k = 0; k < n; k++)
	{
		fmax = 0.0;
		for (i = k; i < n; i++) {
			for (j = k; j < n; j++)
			{
				temp = fabs(tp[i][j]);//找最大值
				if (temp > fmax)
				{
					fmax = temp;
					is[k] = i; js[k] = j;
				}
			}
		}

		if ((fmax + 1.0) == 1.0)
		{
			delete[] is;
			delete[] js;
			return NULL;
		}
		if ((i = is[k]) != k)
			for (j = 0; j < n; j++)
				swap(&tp[k][j], &tp[i][j]);//交换指针
		if ((j = js[k]) != k)
			for (i = 0; i < n; i++)
				swap(&tp[i][k], &tp[i][j]);  //交换指针
		tp[k][k] = 1.0 / tp[k][k];
		for (j = 0; j < n; j++)
			if (j != k)
				tp[k][j] *= tp[k][k];
		for (i = 0; i < n; i++)
			if (i != k)
				for (j = 0; j < n; j++)
					if (j != k)
						tp[i][j] = tp[i][j] - tp[i][k] * tp[k][j];
		for (i = 0; i < n; i++)
			if (i != k)
				tp[i][k] *= -tp[k][k];
	}
	for (k = n - 1; k >= 0; k--)
	{
		if ((j = js[k]) != k)
			for (i = 0; i < n; i++)
				swap(&tp[j][i], &tp[k][i]);
		if ((i = is[k]) != k)
			for (j = 0; j < n; j++)
				swap(&tp[j][i], &tp[j][k]);
	}
	delete[] is;
	delete[] js;
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
		}
	}
	return tp;
}
void swap(double* a, double* b)
{
	double c;
	c = *a;
	*a = *b;
	*b = c;
}

complex<double>** GetMatrixInverse(complex<double>** src, int n) {
	double** A = new double* [n];
	double** B = new double* [n];
	for (int i = 0; i < n; i++) {
		A[i] = new double[n];
		B[i] = new double[n];
		for (int j = 0; j < n; j++) {
			A[i][j] = src[i][j].real();
			B[i][j] = src[i][j].imag();
		}
	}
	double** A1 = inv(A, n);
	double** A1B = Mul(A1, B, n, n, n);
	double** BA1B = Mul(B, A1B, n, n, n);
	double** AjBA1B = Add(A, BA1B, n, n);
	double** AjBA1B_1 = inv(AjBA1B, n);
	double** A1B_AjBA1B_1 = Mul(A1B, AjBA1B_1, n, n, n);

	complex<double>** res = new complex<double> * [n];
	for (int i = 0; i < n; i++) {
		res[i] = new complex<double>[n];
		for (int j = 0; j < n; j++) {
			res[i][j].real(AjBA1B_1[i][j]);
			res[i][j].imag(-1.0 * A1B_AjBA1B_1[i][j]);
		}
	}
	return res;
}

int main() {
	complex<double>** temp = new complex<double> * [3];
	for (int i = 0; i < 3; i++) {
		temp[i] = new complex<double>[3];
	}
	complex<double> t1(1, 0);
	temp[0][0] = t1;
	temp[2][1] = t1;
	temp[2][2] = t1;
	complex<double> t2(1, 1);
	temp[0][2] = t2;
	temp[1][1] = t2;
	complex<double> t3(0, 1);
	temp[1][0] = t3;
	complex<double> t4(2, -1);
	temp[0][1] = t4;
	complex<double> t5(1, 2);
	temp[1][2] = t5;
	complex<double> t6(-1, 1);
	temp[2][0] = t6;

	cout << "原方程：" << endl;
	for (int i = 0; i < 3; i++) {
		for (int j = 0; j < 3; j++) {
			cout << temp[i][j]<<" ";
		}
		cout << endl;
	}
	cout << endl <<"求逆："<<endl;
	complex<double>** res = GetMatrixInverse(temp, 3);
	for (int i = 0; i < 3; i++) {
		for (int j = 0; j < 3; j++) {
			cout << res[i][j] << " ";
		}
		cout << endl;
	}
	system("pause");
	return 0;
}

四、测试

构造复数矩阵：
[[1, 2-i, 1+i],
[i, 1+i, 1+2i],
[-1+i, 1, 1]]
使用本代码测试结果截图：
在这里插入图片描述
使用matlab运行结果截图：
在这里插入图片描述
可以看出运算结果是对的。

跑在路上的蜗牛
关注
  • 10
    点赞
  • 踩
  • 66
    收藏
    觉得还不错? 一键收藏
  • 7
    评论
专栏目录
复数矩阵求逆
09-16
复数矩阵求逆
复数矩阵求逆的 C 语言程序
wangxunweihua的博客
05-20 7015
关于复数矩阵求逆，如果使用 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+...
求逆矩阵（C++）
m0_50793313的博客
02-08 1337
使用C++实现两种求逆矩阵的方法
求逆矩阵程序(c语言)
qq_31605553的博客
12-01 2328
逆矩阵求解
verilog实现四阶复数矩阵求逆
最新发布
Karoski的博客
03-12 844
国科大高等数字集成电路大作业，利用Verilog实现四阶复数矩阵求逆，采用分块矩阵法对四阶常数矩阵求逆，再利用公式完成对复数矩阵的求逆，里面包括二阶矩阵求逆，求积，三级矩阵求行列式，四阶矩阵求积，求逆，求行列式和四阶复矩阵求逆的代码。代码附在txt文档中，请自行复制到程序中
43-左旋转字符串
菜鸟学编程
04-23 231
C++代码： class Solution { public: //实现子字符串的逆序 void InverseString(string &str, int begin, int end){//形参接收字符串&str while(begin < end){ swap(str[begin++],str[end--])...
用初等变换求复数矩阵的逆矩阵
05-22
复数矩阵求逆，用了STL中的vector表示矩阵，用了complex类表示复数，求逆矩阵使用了初等变换求逆矩阵的方法
复数矩阵求逆的C语言程序
m0_37043496的博客
10-16 1044
https://blog.csdn.net/wangxunweihua/article/details/90369914
C语言实现复数矩阵基本操作+LU分解求逆
qq_43756950的博客
03-13 3552
C语言实现（复数）矩阵运算 代码并不是很严谨，没有考虑内存释放之类的问题，初次尝试还有待完善，希望能够帮助到有需要的同学。 1.首先创建复数类型结构体 typedef struct { float real; float imag; }Complex_F; 结构体Complex_F中定义了复数的实部与虚部，这里我使用Complex_F *A或者Complex_F A[row*col]来表示一个矩阵，由于是行排列，所以没有使用二维数组，只需知道A[i][j]=A[i*col+j]即可。 下面
MATLAB中求矩阵的逆矩阵方法（2种）
热门推荐
weixin_40857506的博客
06-25 12万+
方法一：使用inv()函数求矩阵的逆 第一步：打开matlab之后，在命令行窗口中输入a=[1 2 3;4 5 6; 7 8 9]，新建一个a方矩阵，如下图所示： 第二步：在命令行窗口中输入inv(a)，按回车键，可以看到得到了矩阵的逆，如下图所示： 注意：a矩阵可逆的条件是非奇异 方法二：使用a^-1格式求矩阵的逆 第一步：在命令行窗口中输入a^-1，按回车键，可以得到矩阵的逆，如下图所示： 其实，还可以给-1加括号“()”，a^(-1)，如下图： 注：a必须是方针，即行数和列数相等。 ...
Matlab求矩阵的逆(三种方法)
m0_63823993的博客
10-19 5万+
Matlab求矩阵的逆(三种方法)
复数矩阵求逆c/c++实现
06-03
支持任意阶的复数矩阵求逆，a是复数矩阵的实属部分，b是虚数部分，c、d分别是输出实属和虚数部分
c++矩阵求逆的lu分解实现.cpp
01-08
1.求矩阵的crout（LU）分解，其中L为下三角矩阵，U为上三角矩阵 2.求L，U矩阵的伴随阵，参考文献：三角形矩阵求伴随矩阵的一种方法（曾月新） 3.求L，U矩阵的逆（伴随阵A* /det（A）） 4.inv_A = inv_U * inv_L
C++矩阵求逆
05-18
该资源是对于多入多出矩阵求逆的C++编程语言实现，用于多入多出的矩阵求逆算法，本人是用机器人系统设计控制算法，亲自验证代码的可行性
C++实现矩阵的求逆
07-09
利用C++实现矩阵的求逆运算,结果验证跟matlab计算结果差异很小
使用c++写的android矩阵求逆程序
07-27
使用c++写的android矩阵求逆程序，使用c++写的android矩阵求逆程序，使用c++写的android矩阵求逆程序使用c++写的android矩阵求逆程序使用c++写的android矩阵求逆程序使用c++写的android矩阵求逆程序使用c++写的...
北京理工大学复试上机题汇总
srg1011的博客
03-16 8189
秉承前人精神，祝学弟学妹们顺利上岸！！！！！ （00-14和16-17是我转载的别的博主的，15、18和19是我的原创。题意大都描述正确，代码只具有一定参考性，希望大家都自己实现一下） 历年北京理工大学复试上机题（00-14和16/17年） 2015年北理复试上机题 2018年北理复试上机题 2019年北理复试上机题 ...
2019年北理复试上机题
srg1011的博客
03-16 4454
1、碎片字符串 形如aabbaaacaa的字符串，可分为五个相同连续字母组成的碎片：'aa','bb','aaa','c','aa'，其中每个碎片只出现一次，即该字符串包含'aa','bb','aaa','c'四个碎片，且输出时按字典序排序。 样例： 输入：a 输出：a 输入：aabbaaacaa 输出：aa ...
c++矩阵类利用增广矩阵实现矩阵求逆
06-03
好的，您的问题是如何使用增广矩阵来实现矩阵求逆。那么首先需要明确的是，只有方阵才具有逆矩阵，也就是说只有方阵才能使用矩阵求逆的方法。其次，使用增广矩阵来实现矩阵求逆是一种常见的方法，但是并不是唯一的方法。 具体实现的步骤如下： 1. 将原始矩阵 A 与单位矩阵 I 进行水平拼接，得到增广矩阵 B。 2. 对 B 进行初等行变换，将 B 转化为一个上三角矩阵。 3. 对 B 进行进一步的初等行变换，将 B 转化为一个对角线为 1 的矩阵。 4. 对 B 进行进一步的初等行变换，将 B 转化为一个下三角矩阵。 5. 将得到的下三角矩阵的右半部分截取下来，即为所求的矩阵 A 的逆矩阵。 需要注意的是，这种方法需要保证矩阵 A 是可逆的，即其行列式不为 0。否则，该方法将无法求得矩阵 A 的逆矩阵。另外，增广矩阵求逆的时间复杂度为 O(n^3)，其中 n 表示矩阵的阶数。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 7
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值