工作记录——8

其他 2020-04-30 16:38:54 阅读次数: 0

Time : 2020.02.05-2020.02.10

记录c++矩阵求逆的三种方法。

方法1：

#include <iostream>
#include <stdio.h>
#include <cstdio>
#include<vector>
#include <string>
#include <cmath>
#include "opencv2/highgui/highgui.hpp"
using namespace std;

int matrixInversion(vector<vector<double>> &a, int n)
{
    int *is = new int[n];
    int *js = new int[n];
    int i,j,k;
    double d,p;
    for ( k = 0; k < n; k++)
    {
        d = 0.0;
        for (i=k; i<=n-1; i++)
            for (j=k; j<=n-1; j++)
            {
                p=fabs(a[i][j]);
                if (p>d) { d=p; is[k]=i; js[k]=j;}
            }
            if ( 0.0 == d )
            {
                free(is); free(js); printf("err**not inv\n");
                return(0);
            }
            if (is[k]!=k)
                for (j=0; j<=n-1; j++)
                {
                    p=a[k][j];
                    a[k][j]=a[is[k]][j];
                    a[is[k]][j]=p;
                }
            if (js[k]!=k)
                for (i=0; i<=n-1; i++)
                {
                    p=a[i][k];
                    a[i][k]=a[i][js[k]];
                    a[i][js[k]]=p;
                }
            a[k][k] = 1.0/a[k][k];
            for (j=0; j<=n-1; j++)
                if (j!=k)
                {
                    a[k][j] *= a[k][k];
                }
            for (i=0; i<=n-1; i++)
                if (i!=k)
                    for (j=0; j<=n-1; j++)
                        if (j!=k)
                        {
                            a[i][j] -= a[i][k]*a[k][j];
                        }
            for (i=0; i<=n-1; i++)
                if (i!=k)
                {
                    a[i][k] = -a[i][k]*a[k][k];
                }
    }
    for ( k = n-1; k >= 0; k--)
    {
        if (js[k]!=k)
            for (j=0; j<=n-1; j++)
            {
                p = a[k][j];
                a[k][j] = a[js[k]][j];
                a[js[k]][j]=p;
            }
            if (is[k]!=k)
                for (i=0; i<=n-1; i++)
                {
                    p = a[i][k];
                    a[i][k]=a[i][is[k]];
                    a[i][is[k]] = p;
                }
    }
    free(is); free(js);
    return(1);
}


int main()
{
    //prior_testMatrix();	// 这个是测试原来的方法的函数，内容和原来的main函数一致
    int i,j;
    static vector<vector<double>> a(4,vector<double>(4));
    a={{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}};
    static vector<vector<double>> b(4,vector<double>(4));	// 拷贝a

    for (i=0; i< 4; i++){
        for (j=0; j< 4; j++){
            b[i][j]=a[i][j];
            }// 拷贝a
    }

    cout << (i=matrixInversion(b,4)) << endl;	// 计算逆矩阵，结果在b中
    if (i!=0){
        cout<<"矩阵A"<<endl;
        for (i=0; i<=3; i++){
            for (j=0; j<4; j++)
            cout<<a[i][j]<<endl;
        }
        cout<<"矩阵A的逆矩阵"<<endl;
        for (i=0; i<4; i++){
            for (j=0; j<=3; j++)
            cout<<b[i][j]<<endl;
        }
    }
}

方法2：伴随矩阵求逆（耗时）

#include <iostream>
#include <stdio.h>
#include <cstdio>
#include<vector>
#include <string>
#include <cmath>
#include "opencv2/highgui/highgui.hpp"
#define N 3    //测试矩阵维数定义
using namespace std;

//按第一行展开计算|A|
double getA(vector<vector<double> > &arcs,int n)
{
    if(n==1){
        return arcs[0][0];
    }
    double ans = 0;
    vector<vector<double> > temp(N,vector<double>(N));
    temp={{0,0,0},{0,0,0},{0,0,0}};
    int i,j,k;
    for(i=0;i<n;i++){
        for(j=0;j<n-1;j++){
            for(k=0;k<n-1;k++){
                temp[j][k] = arcs[j+1][(k>=i)?k+1:k];
            }
        }
        double t = getA(temp,n-1);
        if(i%2==0){
            ans += arcs[0][i]*t;
        }
        else{
            ans -=  arcs[0][i]*t;
        }
    }
    return ans;
}

//计算每一行每一列的每个元素所对应的余子式，组成A*
void  getAStart(vector<vector<double> > &arcs,int n,vector<vector<double> > &ans)
{
    if(n==1){
        ans[0][0] = 1;
        return;
    }
    int i,j,k,t;
    vector<vector<double> > temp(N,vector<double>(N));
    temp={{0,0,0},{0,0,0},{0,0,0}};

    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            for(k=0;k<n-1;k++){
                for(t=0;t<n-1;t++){
                    temp[k][t] = arcs[k>=i?k+1:k][t>=j?t+1:t];
                }
            }
            ans[j][i] = getA(temp,n-1);
//此处顺便进行了转置
           if((i+j)%2 == 1){
                ans[j][i] = - ans[j][i];
            }
        }
    }
}
//得到给定矩阵src的逆矩阵保存到des中。
bool GetMatrixInverse(vector<vector<double> > &src,int n,vector<vector<double> > &des)
{
    double flag=getA(src,n);
    cout<<"flag "<<flag<<endl;
    vector<vector<double> > t(N,vector<double>(N));
    t={{0,0,0},{0,0,0},{0,0,0}};
    if(0==flag){
        cout<< "原矩阵行列式为0，无法求逆。请重新运行" <<endl;
        return false;//如果算出矩阵的行列式为0，则不往下进行
    }
    else{
        getAStart(src,n,t);
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                des[i][j]=t[i][j]/flag;
            }
        }
    }
    return true;
}

int main()
{
    bool flag;//标志位，如果行列式为0，则结束程序
    vector<vector<double> > matrix_before(N,vector<double>(N));
    //随机数据，可替换
    srand((unsigned)time(0));
    for(int i=0; i<N ;i++)
    {
        for(int j=0; j<N;j++)
        {
            matrix_before[i][j]=rand()%100 *0.01;
        }
    }
    cout<<"原矩阵："<<endl;
    for(int i=0; i<N ;i++){
        for(int j=0; j<N;j++){
            cout <<matrix_before[i][j]<<" ";
        }
        cout<<endl;
    }
    vector<vector<double> > matrix_after(N,vector<double>(N));
    flag=GetMatrixInverse(matrix_before,N,matrix_after);
    if(false==flag)
        return 0;
    cout<<"逆矩阵："<<endl;
    for(int i=0; i<N ;i++){
        for(int j=0; j<N;j++){
            cout <<matrix_after[i][j] <<" ";
            //cout << *(*(matrix_after+i)+j)<<" ";
        }
        cout<<endl;
    }
    GetMatrixInverse(matrix_after,N,matrix_before);
    cout<<"反算的原矩阵："<<endl;//为了验证程序的精度
    for(int i=0; i<N ;i++){
        for(int j=0; j<N;j++){
            cout << matrix_before[i][j]<<" ";
        }
        cout<<endl;
    }
    return 0;
}

将其整理为类实现：

#include <iostream>
#include <stdio.h>
#include <cstdio>
#include <vector>
#include <string>
#include "opencv2/highgui/highgui.hpp"

#define N 3    //测试矩阵维数定义
using namespace std;

class Inv{
public:
    Inv(){};//构造函数
    ~Inv(){};//析构函数
    double getA(vector<vector<double>> &arcs,int n);
    void  getAStart(vector<vector<double>> &arcs,int n,vector<vector<double>> &ans);
    bool GetMatrixInverse(vector<vector<double>> &src,int n,vector<vector<double>> &des);
};

int main(void)
{
    bool flag;//标志位，如果行列式为0，则结束程序
    vector<vector<double>> matrix_before(3,vector<double>(3));
    matrix_before={{1,0,0},{0,1,0},{0,0,1}};
    vector<vector<double> > matrix_after(N,vector<double>(N));
    Inv I;
    flag=I.GetMatrixInverse(matrix_before,N,matrix_after);
    if(false==flag)
        return 0;
    cout<<"逆矩阵："<<endl;
    for(int i=0; i<N ;i++){
        for(int j=0; j<N;j++){
            cout <<matrix_after[i][j] <<" ";
        }
        cout<<endl;
    }
  /*cout<<"反算的原矩阵："<<endl;//为了验证程序的精度
    for(int i=0; i<N ;i++){
        for(int j=0; j<N;j++){
            cout << matrix_before[i][j]<<" ";
        }
        cout<<endl;
    }*/
    return 0;
}

//按第一行展开计算|A|
double Inv::getA(vector<vector<double>> &arcs,int n)
{
    if(n==1){
        return arcs[0][0];
    }
    double ans = 0;
    vector<vector<double>> temp(N,vector<double>(N));
    temp={{0,0,0},{0,0,0},{0,0,0}};
    int i,j,k;
    for(i=0;i<n;i++){
        for(j=0;j<n-1;j++){
            for(k=0;k<n-1;k++){
                temp[j][k] = arcs[j+1][(k>=i)?k+1:k];
            }
        }
        double t = getA(temp,n-1);
        if(i%2==0){
            ans += arcs[0][i]*t;
        }
        else{
            ans -=  arcs[0][i]*t;
        }
    }
    return ans;
}
//计算每一行每一列的每个元素所对应的余子式，组成A*
void  Inv::getAStart(vector<vector<double>> &arcs,int n,vector<vector<double>> &ans)
{
    if(n==1){
        ans[0][0] = 1;
        return;
    }
    int i,j,k,t;
    vector<vector<double>> temp(N,vector<double>(N));
    temp={{0,0,0},{0,0,0},{0,0,0}};

    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            for(k=0;k<n-1;k++){
                for(t=0;t<n-1;t++){
                    temp[k][t] = arcs[k>=i?k+1:k][t>=j?t+1:t];
                }
            }
            ans[j][i] = getA(temp,n-1);
//此处顺便进行了转置
           if((i+j)%2 == 1){
                ans[j][i] = - ans[j][i];
            }
        }
    }
}
//得到给定矩阵src的逆矩阵保存到des中。
bool Inv::GetMatrixInverse(vector<vector<double>> &src,int n,vector<vector<double>> &des)
{
    double flag=getA(src,n);
    cout<<"flag "<<flag<<endl;
    vector<vector<double>> t(N,vector<double>(N));
    t={{0,0,0},{0,0,0},{0,0,0}};
    if(0==flag){
        cout<< "原矩阵行列式为0，无法求逆。请重新运行" <<endl;
        return false;//如果算出矩阵的行列式为0，则不往下进行
    }
    else{
        getAStart(src,n,t);
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                des[i][j]=t[i][j]/flag;
            }
        }
    }
    return true;
}

方法3：高斯-约旦消元法求逆

#include<iostream>
#include<iomanip>
using namespace std;

//左下角变换0 若对角不为1元素下方元素为1交换两行
//对角变1
//右上角变换0
void LD_triangle(float temp[6][12],int n);
void RU_triangle(float temp[6][12],int n);
void Aii_Cto_1(float temp[6][12],int n);
void Display(float A[6][12],int m,int n);
int Find(float temp[6][12],int i,int n);
void Swap(float temp[6][12],int i,int j,int n);

class Matrix{
public:
    Matrix(){
        n=9;
        for(int i=0;i<6;i++)
            for(int j=0;j<12;j++)
                A[i][j]=0;//初始化矩阵A 6*12       存放的是目标矩阵|单位矩阵
        Det_A=0;
    }
    void Conversion();
    void Input();
    bool Cal_DetA();
    friend void Save_To_Result(float temp[6][12],Matrix &D);
    void Cheak_IsSolvable();
    void Out_Result();
private:
    int n;                  //阶
    float A[6][12];         //矩阵
    float result[6][6];     //逆矩阵结果
    float Det_A;            //|A|值：是否存在逆矩阵
};


int main(){
    cout<<"\n\t\t-------求逆矩阵（初等变换法）-------"<<endl<<endl;
    while(1){
        Matrix D;
        D.Input();
        D.Cheak_IsSolvable();   //可解
        D.Conversion();         //变换
        D.Out_Result();         //输出
    }
    return 1;
}
//输入要处理的矩阵以及单位矩阵，可更改输入方式
void Matrix::Input(){
    cout<<"\n---行列数：";
    cin>>n;
    cout<<endl<<"输入行列式："<<endl;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            cin>>A[i][j];
    int j=0;
    for(int i=n;i<2*n;i++){
        A[j][i]=1;
        j++;
    }
}
//检验是否可解
void Matrix::Cheak_IsSolvable(){
    while(1){
        if(Cal_DetA()){
        cout<<"\n|A|="<<Det_A<<endl<<endl;
            cout<<"(A|E)="<<endl;
            Display(A,n,2*n);
            break;
        }
        else{
            cout<<"此矩阵的行列式为0，无解"<<endl;
            Input();
        }
    }
}
void Matrix::Conversion(){
    float temp[6][12];
    for(int i=0;i<n;i++)//替身
        for(int j=0;j<2*n;j++)
            temp[i][j]=A[i][j];
    LD_triangle(temp,n);//左下角变换0
    Aii_Cto_1(temp,n);   //对角变1
    RU_triangle(temp,n);//右上角变换
    Save_To_Result(temp,*this);
}
void LD_triangle(float temp[6][12],int n){
    for(int i=0;i<n-1;i++){
        if(temp[i][i]){//对角线元素不为0
            if(temp[i][i]!=1.0)
                if(int j=Find(temp,i,n)) Swap(temp,i,j,n);
            for(int m=i+1;m<n;m++){
                if(temp[m][i]){
                    float t=-(temp[m][i])/temp[i][i];
                    for(int p=0;p<2*n;p++){//R(m)+tRi
                        temp[m][p]=temp[m][p]+t*temp[i][p];
                    }
                    if(t>0)
                        cout<<"R"<<m+1<<"+"<<t<<"R"<<i+1<<":"<<endl;
                    else
                        cout<<"R"<<m+1<<t<<"R"<<i+1<<":"<<endl;
                    Display(temp,n,2*n);
                }
                else continue;
            }
        }

        else {//若对角线元素为0
            int m;          for(m=i+1;m<n;m++){//
                if(temp[m][i]){//使对角线元素非0
                    for(int p=0;p<n;p++)//Ri+Rm
                        temp[i][p]=temp[m][p]+temp[i][p];
                    break;
                }
                else continue;
            }
            if(m==n){   //可逆矩阵 不会出现此情况（对角元素为零，且此元素以下全为零）
                //此情况 矩阵 行列式值为0
            }
            cout<<"R"<<i+1<<"+"<<"R"<<m+1<<endl;
            Display(temp,n,2*n);
            i--;
        }
    }
}
int Find(float temp[6][12],int i,int n){
    for(int m=i+1;m<n;m++)
        if(temp[m][i]==1.0)
            return m;
    return 0;
}
void Swap(float temp[6][12],int i,int j,int n){
    cout<<"R"<<i+1<<"<-->R"<<j+1<<endl;
    for(int m=0;m<2*n;m++){
        temp[i][m]+=temp[j][m];
        temp[j][m]=temp[i][m]-temp[j][m];
        temp[i][m]-=temp[j][m];
    }
    Display(temp,n,2*n);
}
void Aii_Cto_1(float temp[6][12],int n){
    for(int i=0;i<n;i++){
        if(temp[i][i]!=1.0){
            float t=1.0/temp[i][i];
            cout<<"R"<<i+1<<"/"<<temp[i][i]<<endl;
            for(int p=0;p<2*n;p++){//R(m)+tRi
                temp[i][p]*=t;
            }
        Display(temp,n,2*n);
        }
    }
}
void RU_triangle(float temp[6][12],int n){
    for(int i=n-1;i>0;i--){
        for(int m=i-1;m>=0;m--){
            if(temp[m][i]){
                float t=-(temp[m][i])/(temp[i][i]);
                for(int p=i;p<2*n;p++){ //R(m)+tempRi
                    temp[m][p]=temp[m][p]+t*temp[i][p];
                }
                if(t>0)
                    cout<<"R"<<m+1<<"+"<<t<<"R"<<i+1<<":"<<endl;
                else
                    cout<<"R"<<m+1<<t<<"R"<<i+1<<":"<<endl;
                Display(temp,n,2*n);
            }
            else continue;
        }
    }
}
void Display(float A[6][12],int m,int n){
    for(int i=0;i<m;i++){
        cout<<"|";
        for(int j=0;j<n;j++){
            cout<<setw(8)<<A[i][j];
            if((j+1)%m==0)
                cout<<"  |";
        }
        cout<<endl;
    }
    cout<<endl;
}
void Save_To_Result(float temp[6][12],Matrix &D){
    for(int i=0;i<D.n;i++)
        for(int j=D.n;j<2*D.n;j++)
            D.result[i][j-D.n]=temp[i][j];
}
bool Matrix::Cal_DetA(){
    float temp[6][6];
    for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                temp[i][j]=A[i][j];
    for(int i=0;i<n-1;i++){
        if(temp[i][i]){//对角线元素不为0
            for(int m=i+1;m<n;m++){
                if(temp[m][i]){
                    float tem=-(temp[m][i])/temp[i][i];
                    for(int p=0;p<n;p++)//R(m)+tempRi
                        temp[m][p]=temp[m][p]+tem*temp[i][p];
                }
                else continue;
            }
        }
        else {//若对角线元素为0
            int m;
            for(m=i+1;m<n;m++){//
                if(temp[m][i]){//使对角线元素非0
                    for(int p=0;p<n;p++)//Ri+Rm
                        temp[i][p]=temp[m][p]+temp[i][p];
                    break;
                }
                else continue;
            }
            if(m==n){
                Det_A=0;
                return false;
            }
            i--;
        }
    }
    Det_A=temp[0][0];//求对角积
    for(int i=1;i<n;i++)
        Det_A=Det_A*temp[i][i];

    if(Det_A) return true;
    else return false;
}
void Matrix::Out_Result(){
    cout<<"Result:"<<endl;
    for(int i=0;i<n;i++){
        cout<<"|";
        for(int j=0;j<n;j++){
            cout<<setw(8)<<result[i][j];
        }
        cout<<"  |"<<endl;
    }
    cout<<endl;
}

第一次修改结果：
为便于自己的使用，将矩阵修改为动态矩阵的形式，并在类中实现。
Invmatrix.h

#ifndef INVMATRIX_H
#define INVMATRIX_H
#include<iostream>
#include<iomanip>
#include<vector>
using namespace std;
#define N 3
//左下角变换0 若对角不为1元素下方元素为1交换两行
//对角变1
//右上角变换0

void LD_triangle(vector<vector<double>> &temp,int n);
void RU_triangle(vector<vector<double>> &temp,int n);
void Aii_Cto_1(vector<vector<double>> &temp,int n);
void Display(vector<vector<double>> &A,int m,int n);
int Find(vector<vector<double>> &temp,int i,int n);
void Swap(vector<vector<double>> &temp,int i,int j,int n);

class Matrix{
public:
    Matrix(int _n);
    void Conversion();
    void Input(vector<vector<double>> &B);
    bool Cal_DetA();
    void Cheak_IsSolvable();
    vector<vector<double>> Out_Result();
private:
    int n;                  //阶
    float Det_A;            //|A|值 ：是否存在逆矩阵
};
#endif // INVMATRIX_H

Invmatrix.cpp

#include <Invmatrix.h>
vector<vector<double>> A(3,vector<double>(6));            //全局矩阵;
vector<vector<double>> D(3,vector<double>(3));
vector<vector<double>> temp(3,vector<double>(6));
Matrix::Matrix(int _n)
{
        n=_n;
        Det_A=0;
}

void Matrix::Input(vector<vector<double>> &B){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            A[i][j]=B[i][j];
        }
    }
    int j=0;
    for(int i=n;i<2*n;i++){
        A[j][i]=1;
        j++;
    }
}

void Matrix::Cheak_IsSolvable(){
    while(1){
        if(Cal_DetA()){
        cout<<"\n|A|="<<Det_A<<endl<<endl;
            cout<<"(A|E)="<<endl;
            Display(A,n,2*n);
            break;
        }
        else{
            cout<<"此矩阵的行列式为0，无解"<<endl;
        }
    }
}
void Matrix::Conversion(){
    //vector<vector<double>> temp(n,vector<double>(2*n));
    for(int i=0;i<n;i++)//替身
        for(int j=0;j<2*n;j++)
            temp[i][j]=A[i][j];
    LD_triangle(temp,n);//左下角变换0
    Aii_Cto_1(temp,n);   //对角变1
    RU_triangle(temp,n);//右上角变换
}
void LD_triangle(vector<vector<double>> &temp,int n){
    for(int i=0;i<n-1;i++){
        if(temp[i][i]){//对角线元素不为0
            if(temp[i][i]!=1.0)
                if(int j=Find(temp,i,n))
                    Swap(temp,i,j,n);
            for(int m=i+1;m<n;m++){
                if(temp[m][i]){
                    float t=-(temp[m][i])/temp[i][i];
                    for(int p=0;p<2*n;p++){//R(m)+tRi
                        temp[m][p]=temp[m][p]+t*temp[i][p];
                    }
                    if(t>0)
                        cout<<"R"<<m+1<<"+"<<t<<"R"<<i+1<<":"<<endl;
                    else
                        cout<<"R"<<m+1<<t<<"R"<<i+1<<":"<<endl;

                    Display(temp,n,2*n);
                }
                else continue;
            }
        }

        else {//若对角线元素为0
            int m;
            for(m=i+1;m<n;m++){//
                if(temp[m][i]){//使对角线元素非0
                    for(int p=0;p<n;p++)//Ri+Rm
                        temp[i][p]=temp[m][p]+temp[i][p];
                    break;
                }
                else continue;
            }
            if(m==n){   //可逆矩阵 不会出现此情况（对角元素为零，且此元素以下全为零）
                //此情况 矩阵 行列式值为0
            }
            cout<<"R"<<i+1<<"+"<<"R"<<m+1<<endl;
            Display(temp,n,2*n);
            i--;
        }
    }
}

int Find(vector<vector<double>> &temp,int i,int n){
    for(int m=i+1;m<n;m++)
        if(temp[m][i]==1.0)
            return m;
    return 0;
}

void Swap(vector<vector<double>> &temp,int i,int j,int n){
    cout<<"R"<<i+1<<"<-->R"<<j+1<<endl;
    for(int m=0;m<2*n;m++){
        temp[i][m]+=temp[j][m];
        temp[j][m]=temp[i][m]-temp[j][m];
        temp[i][m]-=temp[j][m];
    }
    Display(temp,n,2*n);
}

void Aii_Cto_1(vector<vector<double>> &temp,int n){
    for(int i=0;i<n;i++){
        if(temp[i][i]!=1.0){
            float t=1.0/temp[i][i];
            cout<<"R"<<i+1<<"/"<<temp[i][i]<<endl;
            for(int p=0;p<2*n;p++){//R(m)+tRi
                temp[i][p]*=t;
            }
        Display(temp,n,2*n);
        }
    }
}
void RU_triangle(vector<vector<double>> &temp,int n){
    for(int i=n-1;i>0;i--){
        for(int m=i-1;m>=0;m--){
            if(temp[m][i]){
                float t=-(temp[m][i])/(temp[i][i]);
                for(int p=i;p<2*n;p++){ //R(m)+tempRi
                    temp[m][p]=temp[m][p]+t*temp[i][p];
                }
                if(t>0)
                    cout<<"R"<<m+1<<"+"<<t<<"R"<<i+1<<":"<<endl;
                else
                    cout<<"R"<<m+1<<t<<"R"<<i+1<<":"<<endl;
                Display(temp,n,2*n);
            }
            else continue;
        }
    }
}
void Display(vector<vector<double>> &A,int m,int n){
    for(int i=0;i<m;i++){
        cout<<"|";
        for(int j=0;j<n;j++){
            cout<<setw(8)<<A[i][j];
            if((j+1)%m==0)
                cout<<"  |";
        }
        cout<<endl;
    }
    cout<<endl;
}
bool Matrix::Cal_DetA()
{
  //  vector<vector<double>> temp(n,vector<double>(n));
    for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                temp[i][j]=A[i][j];//复制矩阵A左侧（原始矩阵）
            }
    }
    for(int i=0;i<n-1;i++){
        if(temp[i][i]){//对角线元素不为0
            for(int m=i+1;m<n;m++){
                if(temp[m][i]){
                    float tem=-(temp[m][i])/temp[i][i];
                    for(int p=0;p<n;p++)//R(m)+tempRi
                        temp[m][p]=temp[m][p]+tem*temp[i][p];
                }
                else continue;
            }
        }
        else {//若对角线元素为0
            int m;
            for(m=i+1;m<n;m++){//如果对角元素i为0;，那么判断下一行的此列元素是否为0，若不为零，则将该行元素分别与下一行的元素相加。
                if(temp[m][i]){//使对角线元素非0
                    for(int p=0;p<n;p++)//Ri+Rm
                        temp[i][p]=temp[m][p]+temp[i][p];
                    break;
                }
                else continue;
            }
            if(m==n){
                Det_A=0;//因为如果该对角元素为0，表示其上行的该列元素也为0，那么若果其下所有行在此列元素都为0，表名该列所有元素均为0，则该矩阵无解；
                return false;
            }
            i--;//执行完该语句后，表示else语句执行完成，进入 for(int i=0;i<n-1;i++)循环
        }
    }
    Det_A=temp[0][0];//求对角积
    for(int i=1;i<n;i++)
        Det_A=Det_A*temp[i][i];

    if(Det_A) return true;
    else return false;
}

vector<vector<double>> Matrix::Out_Result(){

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            D[i][j]=0;
        }
    }

    for(int i=0;i<n;i++){
        int p=0;
        for(int j=n;j<2*n;j++){
            D[i][p]=temp[i][j];
            p++;
        }
    }
    return D;
}

main.cpp

#include <Invmatrix.h>
#include<iostream>
#include<iomanip>
#include<vector>
using namespace std;


int main(){
    cout<<"\n\t\t-------求逆矩阵（初等变换法）-------"<<endl<<endl;
        vector<vector<double>> a(3,vector<double>(3));
        a={{1,2,0},{1,1,0},{0,0,1}};
        Matrix D(3);
        D.Input(a);
        D.Cheak_IsSolvable();   //可解
        D.Conversion();         //变换
        vector<vector<double>> M(3,vector<double>(3));
        M=D.Out_Result();         //输出
        for(int i=0;i<3;i++){
            for(int j=0;j<3;j++){
                cout<<M[i][j]<<endl;
            }
        }
        return 0;
}

第二次修改：
删掉了一些不必要的成员函数，如friend void Save_To_Result()。因为要用到最终的可逆矩阵，因此将结果输出函数设置成vector<vector> Out_Result()，此外，使用该程序不需要进行每一个步骤的显示，只需得到最终的结果，因此去掉一系列的显示，精简程序：
Invmatrix.h

#ifndef INVMATRIX_H
#define INVMATRIX_H
#include<iostream>
#include<iomanip>
#include<vector>
using namespace std;
#define N 3
//左下角变换0 若对角不为1元素下方元素为1交换两行
//对角变1
//右上角变换0
void LD_triangle(vector<vector<double>> &temp,int n);
void RU_triangle(vector<vector<double>> &temp,int n);
void Aii_Cto_1(vector<vector<double>> &temp,int n);
int Find(vector<vector<double>> &temp,int i,int n);
void Swap(vector<vector<double>> &temp,int i,int j,int n);

class Matrix{
public:
    Matrix(int _n);
    void Conversion();
    void Input(vector<vector<double>> &B);
    bool Cal_DetA();
    void Cheak_IsSolvable();
    vector<vector<double>> Out_Result();
private:
    int n;                  //阶
    float Det_A;            //|A|值 ：是否存在逆矩阵
};
#endif // INVMATRIX_H

Invmatrix.cpp

扫描二维码关注公众号，回复： 11157613 查看本文章 
#include <Invmatrix.h>
vector<vector<double>> A(3,vector<double>(6));            //全局矩阵;
vector<vector<double>> D(3,vector<double>(3));
vector<vector<double>> temp(3,vector<double>(6));
Matrix::Matrix(int _n)
{
        n=_n;
        Det_A=0;
}

void Matrix::Input(vector<vector<double>> &B){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            A[i][j]=B[i][j];
        }
    }
    int j=0;
    for(int i=n;i<2*n;i++){
        A[j][i]=1;
        j++;
    }
}

void Matrix::Cheak_IsSolvable(){
    while(1){
        if(Cal_DetA()){
        cout<<"\n|A|="<<Det_A<<endl<<endl;
            break;
        }
        else{
            cout<<"此矩阵的行列式为0，无解"<<endl;
        }
    }
}
void Matrix::Conversion(){
    //vector<vector<double>> temp(n,vector<double>(2*n));
    for(int i=0;i<n;i++)//替身
        for(int j=0;j<2*n;j++)
            temp[i][j]=A[i][j];
    LD_triangle(temp,n);//左下角变换0
    Aii_Cto_1(temp,n);   //对角变1
    RU_triangle(temp,n);//右上角变换
}
void LD_triangle(vector<vector<double>> &temp,int n){
    for(int i=0;i<n-1;i++){
        if(temp[i][i]){//对角线元素不为0
            if(temp[i][i]!=1.0)
                if(int j=Find(temp,i,n))
                    Swap(temp,i,j,n);
            for(int m=i+1;m<n;m++){
                if(temp[m][i]){
                    float t=-(temp[m][i])/temp[i][i];
                    for(int p=0;p<2*n;p++){//R(m)+tRi
                        temp[m][p]=temp[m][p]+t*temp[i][p];
                    }
                }
                else continue;
            }
        }
        else {//若对角线元素为0
            int m;
            for(m=i+1;m<n;m++){//
                if(temp[m][i]){//使对角线元素非0
                    for(int p=0;p<n;p++)//Ri+Rm
                        temp[i][p]=temp[m][p]+temp[i][p];
                    break;
                }
                else continue;
            }
            if(m==n){   //可逆矩阵 不会出现此情况（对角元素为零，且此元素以下全为零）
                //此情况 矩阵 行列式值为0
            }
            i--;
        }
    }
}

int Find(vector<vector<double>> &temp,int i,int n){
    for(int m=i+1;m<n;m++)
        if(temp[m][i]==1.0)
            return m;
    return 0;
}

void Swap(vector<vector<double>> &temp,int i,int j,int n){
    for(int m=0;m<2*n;m++){
        temp[i][m]+=temp[j][m];
        temp[j][m]=temp[i][m]-temp[j][m];
        temp[i][m]-=temp[j][m];
    }
}

void Aii_Cto_1(vector<vector<double>> &temp,int n){
    for(int i=0;i<n;i++){
        if(temp[i][i]!=1.0){
            float t=1.0/temp[i][i];
            for(int p=0;p<2*n;p++){//R(m)+tRi
                temp[i][p]*=t;
            }
        }
    }
}

void RU_triangle(vector<vector<double>> &temp,int n){
    for(int i=n-1;i>0;i--){
        for(int m=i-1;m>=0;m--){
            if(temp[m][i]){
                float t=-(temp[m][i])/(temp[i][i]);
                for(int p=i;p<2*n;p++){ //R(m)+tempRi
                    temp[m][p]=temp[m][p]+t*temp[i][p];
                }
            }
            else continue;
        }
    }
}

bool Matrix::Cal_DetA()
{
  //  vector<vector<double>> temp(n,vector<double>(n));
    for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                temp[i][j]=A[i][j];//复制矩阵A左侧（原始矩阵）
            }
    }
    for(int i=0;i<n-1;i++){
        if(temp[i][i]){//对角线元素不为0
            for(int m=i+1;m<n;m++){
                if(temp[m][i]){
                    float tem=-(temp[m][i])/temp[i][i];
                    for(int p=0;p<n;p++)//R(m)+tempRi
                        temp[m][p]=temp[m][p]+tem*temp[i][p];
                }
                else continue;
            }
        }
        else {//若对角线元素为0
            int m;
            for(m=i+1;m<n;m++){//如果对角元素i为0;，那么判断下一行的此列元素是否为0，若不为零，则将该行元素分别与下一行的元素相加。
                if(temp[m][i]){//使对角线元素非0
                    for(int p=0;p<n;p++)//Ri+Rm
                        temp[i][p]=temp[m][p]+temp[i][p];
                    break;
                }
                else continue;
            }
            if(m==n){
                Det_A=0;//因为如果该对角元素为0，表示其上行的该列元素也为0，那么若果其下所有行在此列元素都为0，表名该列所有元素均为0，则该矩阵无解；
                return false;
            }
            i--;//执行完该语句后，表示else语句执行完成，进入 for(int i=0;i<n-1;i++)循环
        }
    }
    Det_A=temp[0][0];//求对角积
    for(int i=1;i<n;i++)
        Det_A=Det_A*temp[i][i];

    if(Det_A) return true;
    else return false;
}
vector<vector<double>> Matrix::Out_Result(){

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            D[i][j]=0;
        }
    }
    for(int i=0;i<n;i++){
        int p=0;
        for(int j=n;j<2*n;j++){
            D[i][p]=temp[i][j];
            p++;
        }
    }
    return D;
}

main.cpp

#include <Invmatrix.h>
#include<iostream>
#include<iomanip>
#include<vector>
using namespace std;
int main(){
    cout<<"\n\t\t-------求逆矩阵（初等变换法）-------"<<endl<<endl;
        vector<vector<double>> a(3,vector<double>(3));
        a={{1,2,0},{1,1,0},{0,0,1}};
        Matrix D(3);
        D.Input(a);
        D.Cheak_IsSolvable();   //可解
        D.Conversion();         //变换
        vector<vector<double>> M(3,vector<double>(3));
        M=D.Out_Result();         //输出
        for(int i=0;i<3;i++){
            for(int j=0;j<3;j++){
                cout<<M[i][j]<<endl;
            }
        }
        return 0;
}
罗宾酱
原创文章 15 获赞 8 访问量 916
关注 私信

猜你喜欢

转载自blog.csdn.net/weixin_39652282/article/details/104181625
今日推荐
周排行
Leetcode简单题61~80
解决zookeeper磁盘IO高的问题
多线程相关方法详解
Maven-setting.xml文件详解
Maven 项目的 classpath 理解
渊亭科技大数据笔试题
配置JVM内存分配
计算机网络个人学习笔记 （三）网络层 ：第三部分 连载
js中两个等号(==)和三个等号(===)的区别
用C程序自动打开电脑上的程序
每日归档
更多
2024-09-18(0)
2024-09-17(0)
2024-09-16(0)
2024-09-15(0)
2024-09-14(0)
2024-09-13(0)
2024-09-12(0)
2024-09-11(0)
2024-09-10(0)
2024-09-09(0)

代码天地 © 2018 codetd.com

境外服务器   最新电视剧2022