代码实现矩阵求逆的三种方式(超详细、已实现)

一、高斯消元法

         对于任意一个矩阵Anxn,其满足。基于此,高斯消元法具体步骤是先构造一个增广矩阵W=[E],则W为一个n x 2n的矩阵。我们需要对矩阵W进行矩阵行之间的变换,将其变为 [E|B] 的形势,如果能够成功变换,则B就为A矩阵的逆矩阵。

具体操作过程如下:

         (1)、将初始矩阵A右半部分进行扩增,得到矩阵W= [A | E],W为 nx2n。

         (2)、将首行作为基准,从上往下做行变换,将W前半部分转化为一个上三角矩阵。

         (3)、将W前半部分由上三角矩阵转化为对角矩阵。

         (4)、对W前半部分的对角矩阵乘以一个系数将其转化为单位矩阵。

         (5)、提出W矩阵的后半部分,就得到了我们想要的A矩阵的逆矩阵。

         转化过程中,如果发现W矩阵前半部分不能单位化,则判定A矩阵无可逆矩阵。

  具体代码实现如下:

//使用高斯消元法对矩阵进行求逆
void Gaussian_elimination(int arr[N][N])
{
	int i, j, k;
	float W[N][2*N], result[N][N];
	float tem_1, tem_2, tem_3;

	// 对矩阵右半部分进行扩增
	for(i = 0;i < N; i++){
		for(j = 0;j < 2 * N; j++){
			if(j<N){
				W[i][j] = (float) arr[i][j];
			}
			else{
				W[i][j] = (float) (j-N == i ? 1:0);
			}
		}
	}

	for(i=0;i<N;i++)
	{
		// 判断矩阵第一行第一列的元素是否为0,若为0,继续判断第二行第一列元素,直到不为0,将其加到第一行
		if( ((int) W[i][i]) == 0)
		{ 
			for(j=i+1;j<N;j++)
			{
				if( ((int) W[j][i]) != 0 ) break;
			}
			if(j == N)
			{
				printf("这个矩阵不能求逆");
				break;
			}
			//将前面为0的行加上后面某一行
			for(k=0;k<2*N;k++)
			{
				W[i][k] += W[j][k];
			}
		}

		//将前面行首位元素置1
		tem_1 = W[i][i];
		for(j=0;j<2*N;j++)
		{
			W[i][j] = W[i][j] / tem_1;
		}

		//将后面所有行首位元素置为0
		for(j=i+1;j<N;j++)
		{
			tem_2 = W[j][i];
			for(k=i;k<2*N;k++)
			{
				W[j][k] = W[j][k] - tem_2 * W[i][k];
			}
		}
	}

	// 将矩阵前半部分标准化
	for(i=N-1;i>=0;i--)
	{
		for(j=i-1;j>=0;j--)
		{
			tem_3 = W[j][i];
			for(k=i;k<2*N;k++)
			{
				W[j][k] = W[j][k] - tem_3*W[i][k];
			}
		}
	}

	//得出逆矩阵
	for(i=0;i<N;i++)
	{
		for(j=N;j<2*N;j++)
		{
			result[i][j-N] = W[i][j];
		}
	}

	printf("使用高斯消元法矩阵求逆的结果为:\n");
	for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			printf("%7.2f", result[i][j]);
		}
		printf("\n");
	}

	
}

二、伴随矩阵法求逆

        一个方阵A如果满足 |A| != 0,则可以认为矩阵A可逆,其逆矩阵为:

         使用伴随矩阵求逆法最关键的一步是如何求矩阵A的伴随矩阵A*,A*求解如下图:

 具体代码实现过程如下:

// 采用伴随矩阵的方法对矩阵进行求逆
void adjoint_matrix(int arr[N][N])
{
	int i, j, k, d, m, n, x, y;
	float W[N][N], result[N][N], tem_[N-1][N-1];
	float det_value, tmp;

	// 将整数类型强转为浮点类型方便计算
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		{
			W[i][j] = (float) arr[i][j];
		}
	}

	// 计算矩阵的行列式 |A|
	det_value = calculate_metrix(W);
	if((int) det_value == 0) printf("这个矩阵没有逆矩阵!");

	// 求矩阵的余子式矩阵
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		{
			//计算每个元素的余子式
			k = 0;
			for(m=0;m<N-1;m++)
			{
				if(k == i) k++;
				d = 0;
				for(n=0;n<N-1;n++)
				{
					if(d==j) d++;
					tem_[m][n] = W[k][d];
					d++;
				}
				k++;
			}

			//开始计算余子式矩阵
			tmp = calculate_metrix_1(tem_);

			//计算完后进行转置  
			result[j][i] = pow(-1, i+j) * tmp;
		}
	}
	
	printf("使用伴随矩阵求逆的结果为:\n");
	for(x=0;x<N;x++)
	{
		for(y=0;y<N;y++)
		{
			printf("%7.2f", result[x][y]/det_value);
		}
		printf("\n");
	}
}

三、LU分解法求逆

        LU分解法:

  

 算法的具体步骤:

        (1)、将 矩阵分解为 下三角矩阵和 上三角矩阵。

step1. L对角线填充1

step2. 

step3.  

step4. U是按行迭代计算,L是按列迭代计算,UL交错计算,且U先L一步。

        (2)、分别对 L U 求逆矩阵,得 L_inv 和 U_inv

step1. 列顺序行顺序, for  j = 0  to  m-1, i = j  to  m-1

 step2. 列顺序行倒序, for j = 0 to m-1, i = j to 0

 

        (3)、A_inv = U_inv * L_inv

具体代码实现:

// 采用LU分解法来对矩阵求逆
void LU_decomposition(int arr[N][N])
{
	float W[N][N], W_n[N][N], L[N][N], U[N][N], L_n[N][N], U_n[N][N];
	int i, j, k, d;
	float s;

	// 赋初值
	for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			W[i][j] = (float) arr[i][j];
			L[i][j] = 0;
			U[i][j] = 0;
			L_n[i][j] = 0;
			U_n[i][j] = 0;
			W_n[i][j] = 0;
		}
	}

	for(i=0;i<N;i++)  // L对角置1
	{
		L[i][i] = 1.0;
	}

	for(j=0;j<N;j++)  
	{
		U[0][j] = W[0][j];
	}

	for(i=1;i<N;i++)
	{
		L[i][0] = W[i][0] / U[0][0];
	}

	for(i=1;i<N;i++)
	{
		for(j=i;j<N;j++) // 求U
		{
			s = 0;
			for(k=0;k<i;k++)
			{
				s += L[i][k] * U[k][j];
			}
			U[i][j] = W[i][j] - s;
		}

		for(d=i;d<N;d++) // 求L
		{
			s = 0;
			for(k=0;k<i;k++)
			{
				s += L[d][k] * U[k][i];
			}
			L[d][i] = (W[d][i] - s) / U[i][i];
		}
	}

	for(j=0;j<N;j++)  //求L的逆
	{
		for(i=j;i<N;i++)
		{
			if(i==j) 
				L_n[i][j] = 1 / L[i][j];
			else if(i<j) 
				L_n[i][j] = 0;
			else
			{
				s = 0.;
				for(k=j;k<i;k++)
				{
					s += L[i][k] * L_n[k][j];
				}
				L_n[i][j] = -L_n[j][j] * s;
			}
		}
	}

	for(i=0;i<N;i++)  //求U的逆
	{
		for(j=i;j>=0;j--)
		{
			if(i==j)
				U_n[j][i] = 1 / U[j][i];
			else if(j>i) 
				U_n[j][i] = 0;
			else
			{
				s = 0.;
				for(k=j+1;k<=i;k++)
				{
					s += U[j][k] * U_n[k][i];
				}
				U_n[j][i] = -1 / U[j][j] * s;
			}
		}
	}


	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		{
			for(k=0;k<N;k++)
			{
				W_n[i][j] += U_n[i][k] * L_n[k][j];
			}
		}
	}

	printf("使用LU分解法求逆的结果为:\n");
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		{
			printf("%7.2f", W_n[i][j]);
		}
		printf("\n");
	}
}

四、三种方法求逆的结果:

参考链接:矩阵求逆-高斯消元法介绍及其实现_Your Blog-CSDN博客_高斯消元法求逆矩阵

                    求逆矩阵 —— LU分解法_Lancetop的博客-CSDN博客_lu分解求逆矩阵

  • 4
    点赞
  • 0
    评论
  • 20
    收藏
  • 打赏
    打赏
  • 扫一扫,分享海报

参与评论 您还未登录,请先 登录 后发表或查看评论
©️2022 CSDN 皮肤主题:深蓝海洋 设计师:CSDN官方博客 返回首页

打赏作者

爱学习的杨子

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值