【矩阵求逆】_逆矩阵数值稳定-CSDN博客

本文链接：https://blog.csdn.net/qq_45449625/article/details/126843026

[功能]
用全选主元高斯-约当(Gauss-Jordan)消去法求n阶实(或复)矩阵A的逆矩阵 $A^{-1}$ 。

【方法说明】
高斯-约当法求逆过程可以用以下计算步骤表示。

对于k=1,2,…,n做如下运算。

(1)归一化计算。

$1 / a_{k k} \Rightarrow a_{k k} \\ a_{k j} a_{k k} \Rightarrow a_{k j}, \quad j=1,2, \cdots, n ; j \neq k \\$
(2)消元计算。

$\begin{array}{ll}\\ a_{i j}-a_{i k} a_{k j} \Rightarrow a_{i j}, & i=1,2, \cdots, n ; i \neq k \\ & j=1,2, \cdots, n ; j \neq k \\ -a_{i k} a_{k k} \Rightarrow a_{i k}, & i=1,2, \cdots, n ; i \neq k\\ \end{array}$
矩阵求逆的计算工作量(乘除法次数)为O( $n^{3}$ )。
为了数值计算的稳定性，整个求逆过程需要全选主元。全选主元的过程如下:对于矩阵求逆过程中的第k步,首先,在 $a_{kk}$ 右下方(包括 $a_{kk}$ )的n-k+1阶子阵中选取绝对值最大的元素，并将该元素所在的行号记录在IS(k)中,列号记录在JS(k)中。然后,通过行交换与列交换将该绝对值最大的元素交换到主对角线 $a_{kk}$ 的位置上，即做以下两步:

$(1)\boldsymbol{A}(k, l) \Leftrightarrow \boldsymbol{A}(\operatorname{IS}(k), l), l=1,2, \cdots, n \\ (2)\boldsymbol{A}(l, k) \Leftrightarrow \boldsymbol{A}(l, J S(k)), l=1,2, \cdots, n$

经过全选主元后, 矩阵求逆的计算过程是稳定的。但最后需要对结果进行恢复。恢复的过程如下：对于 k 从 n 到 1 , 分别做以下两步。

$(1) \boldsymbol{A}(k, l) \Leftrightarrow \boldsymbol{A}(J \mathrm{~S}(k), l), l=1,2, \cdots n 。\\ (2) \boldsymbol{A}(l, k) \Leftrightarrow \boldsymbol{A}(l, I S(k)), l=1,2, \cdots, n 。$

//矩阵求逆.cpp
#include <cmath>
#include <iostream>
using namespace std; 
class complex
{
	private:
		double R;	//复数的实部 
		double I;	//复数的虚部
	
	public:
		
		complex(double real=0,double image=0)  //构造函数 
		{
			R=real;
			I=image;
		}
		
		void print()	//复数输出
		{
			cout<<"("<<R<<","<<I<<")" ;    //输出为（R,I）
			return ; 
		 } 
		 
		double cfabs()	//复数的模 
		{
			double y;
			y=sqrt(R*R+I*I);
			return 0;
		}
		
		double angle()	//复数幅角 
		{
			double y;
			y=atan2(I,R);
			return y;
		}
		
		complex operator + (complex &c2)	//复数加法;重载运算符+ 
		{
			complex c;
			c.R=R+c2.R;
			c.I=I+c2.I;
			return c;
		}
		
		complex operator - (complex &c2)	//复数减法;重载运算符-
		{
			complex c;
			c.R=R-c2.R;
			c.I=I-c2.I;
			return c;
		}
		
		complex operator * (complex &c2)	//复数乘法;重载运算符 *
		{
			complex c;
			double p,q,s;
			p=R*c2.R;
			q=I*c2.I;
			s=(R+I)*(c2.R+c2.I);
			c.R=p-q;
			c.I=s-p-q;
			return c;
		}
		
		complex operator / (complex &c2)	//复数除;重载运算符 /
		{
			complex c;
			double p,q,s,w;
			p=R*c2.R;
			q=-I*c2.I;
			s=(R+I)*(c2.R-c2.I);
			w=(c2.R)*(c2.R)+(c2.I)*(c2.I);
			if(w+1.0 !=1.0)
			{
				c.R=(p-q)/w;
				c.I=(s-p-q)/w;
			}
			else 
			{
				c.R=1e+300;
				c.I=1e+300;
			}
			return c;
		}
		
		complex Power(int n)	//复数乘幂
		{
			complex c;
			double r,q;
			q=atan2(I,R);
			r=sqrt(R*R+I*I );
			if(r+1.0 !=1.0)
			{
				r=n*log(r);
				r=exp(r);
			}
			c.R =r*cos(n*q);
			c.I =r*sin(n*q);
			return c;
		} 
		 
		void Root(int n,complex *p)	//复数的 n 次方根
		{
			complex c;
			int k;
			double r,q,t;
			if(n<1)	return ;
			q=atan2(I,R);
			r=sqrt(R*R+I*I);
			if(r+1.0 !=1.0)
			{
				r=(1.0/n)*log(r);
				r=exp(r);
			}
			for(k=0;k<n;k++)
			{
				t=(2.0*k*3.1415926987+q)/n;
				c.R=r*cos(t);
				c.I=r*sin(t);
				p[k]=c;
			}
			 	
		} 
		
		complex Exp()	//复数指数
		{
			complex c;
			double p;
			p=exp(R);
			c.R=p*cos(I);
			c.I=p*sin(I);
			return c;
		} 
		
		complex Log()	//复数对数
		{
			complex c;
			double p;
			p=R*R+I*I;
			p=log(sqrt(p));
			c.R =p;
			c.I =atan2(I,R);
			return c;	
		} 
		
		complex Sin()	//复数正弦
		{
			complex c;
			double p,q;
			p=exp(I);
			q=exp(-I);
			c.R =sin(R)*(p+q)/2;
			c.I =cos(R)*(p-q)/2;
			return c;	
		} 
		
		complex Cos()	//复数余弦
		{
			complex c;
			double p,q;
			p=exp(I);
			q=exp(-I);
			c.R =cos(R)*(p+q)/2;
			c.I =-sin(R)*(p-q)/2;
			return c;	
		} 		
 } ;
 
double ffabs (double p)
{
	double q; 
	q=fabs(p); 
	return(q);
}

double ffabs(complex p)	//计算复数的模	
{
	double q; 
	q=p.cfabs();
	return(q);
}
double ff(double p)	//计算1.0/p	
{
double q; 
q=1.0/p; return(q);
}

complex ff(complex p)	//计算(1+j0)/p	
{
	complex q;
	q=complex(1.0,0.0)/p; 
	return(q);
}
//a 原矩阵,返回逆矩阵
//n 矩阵阶数
template <class T>	//模板声明T为类型参数	
int inv(T a[], int n)	//若矩阵奇异,则返回标志值0,否则返回标志值非 0	
{
	int *is,*js,i,j,k,l,u,v; 
	double d, q; 
	T p;
	is=new int[n]; 
	js=new int[n];
	for (k=0; k<=n-1; k++)
  {
	d=0.0;
	for (i=k; i<=n-1; i++)	//选主元
	{
		d=0.0;
		for (i=k; i<=n-1; i++)
		for (j=k; j<=n-1; j++)
		{
		l=i*n+j;
		q=ffabs(a[l]);	//计算元素绝对值(模)	
		if(q>d)
			{
			d=q;
			is[k]=i;
			js[k]=j;
			}
		}
		if(d+1.0==1.0)	//矩阵奇异	
		{
		delete[] is; delete[] js; 
		cout<<"矩阵奇异!" <<endl;
		return(0);	//返回奇异标志值
		}	
		if (is[k]!=k)
		for (j=0;j<=n-1;j++) //行交换
		{
		u=k*n+j;v=is[k]*n+j;
	 	p=a[u];a[u]=a[v];
		a[v]=p;
		}
		if(js[k]!=k)
		for(i=0;i<=n-1;i++) //列交换
		{	
		 u=i*n+k;v=i*n+js[k];
		 p=a[u];a[u]=a[v]; 
		 a[v]=p;
		}
		l=k*n+k;
		a[l]=ff(a[l]);	//计算1/a[1]	
		for(j=0;j<=n-1;j++) //归一化 
		if(j!=k)
		{
		u=k*n+j;
		a[u]=a[u]*a[1];
		}
		for (i-0; i<=n-1;i++) //消元计算 
		if(i!=k)
		for (j=0;j<=n-1;j++) 
		if (j!=k)
		{
		u=i*n+j;
		a[u]=a[u]-a[i*n+k]*a[k*n+j];
		}
		for (i=0;i<=n-1;i++) 
		if (i!=k)
		{
		u=i*n+k; 
		a[u]=(a[u]-a[u]-a[u])*a[l];
		}
	}
 }
for (k=n-1; k>=0; k--)	//恢复行列交换	
	{
		if(js[k]!=k)
		for (j=0; j<=n-1; j++)
			{
			u=k*n+j; 
			v=js[k]*n+j; 
			p=a[u];a[u]=a[v];
			a[v]=p;
			}
		if (is[k]!=k)
		for (i=0; i<=n-1; i++)
			{	
			u=i*n+k;
			v=i*n+is[k]; 
			p=a[u];
			a[u]=a[v];
			a[v]=p;	
			}
	}
	delete[] is;
	delete[] js; 
	return(1);
}
int main()
{
	int i,j;
	double a[4][4]={{0.2368,0.2471,0.2568,1.2671},
					{1.1161,0.1254,0.1397,0.1490},
					{0.1582,1.1675,0.1768,0.1871},
					{0.1968,0.2071,1.2168,0.2271}};
	double b[4][4];
	for(i=0;i<=3;i++)
	for(j=0;j<=3;j++) b[i][j]=a[i][j]; 
	i=inv(&b[0][0],4); 
	if (i!=0)
	{
		cout<<"实矩阵 A:" <<endl; 
		for(i=0;i<=3;i++)
		{
			for (j=0; j<=3; j++) 
			cout <<a[i][j]<<"	"; 
			cout <<endl;
		}
			
		cout<<"逆矩阵A^(-1):"<<endl; 
		for(i=0;i<=3;i++)
		{
			for (j=0; j<=3; j++) cout <<b[i][j]<<"	";
			cout <<endl;
		}
	}
	return 0;
}