[功能]
用全选主元高斯-约当(Gauss-Jordan)消去法求n阶实(或复)矩阵A的逆矩阵。
【方法说明】
高斯-约当法求逆过程可以用以下计算步骤表示。
对于k=1,2,…,n做如下运算。
(1)归一化计算。
(2)消元计算。
矩阵求逆的计算工作量(乘除法次数)为O()。
为了数值计算的稳定性,整个求逆过程需要全选主元。全选主元的过程如下:对于矩阵求逆过程中的第k步,首先,在右下方(包括
)的n-k+1阶子阵中选取绝对值最大的元素,并将该元素所在的行号记录在IS(k)中,列号记录在JS(k)中。然后,通过行交换与列交换将该绝对值最大的元素交换到主对角线
的位置上,即做以下两步:
经过全选主元后, 矩阵求逆的计算过程是稳定的。但最后需要对结果进行恢复。恢复 的过程如下:对于 k 从 n 到 1 , 分别做以下两步。
//矩阵求逆.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;
}