用C++求矩阵的特征值

怎么用求矩阵行列式的函数求矩阵特征值啊
我知道求特征特是要用原矩阵减特征值乘单位矩阵的行列式等于0
如果这个函数是输入一个矩阵输出它的行列式,那么求特征值能用到这样一个函数么
用程序求特征值是什么思路啊

当使用C++来求解矩阵的特征值时,可以通过以下步骤进行:

1、准备矩阵:首先,根据输入的矩阵数据构建一个 n x n 的矩阵 A。你可以使用数组或矩阵类库(如Eigen)来表示和操作矩阵。
2、计算特征多项式:特征多项式是一个关于 λ 的多项式,其中 λ 是特征值。你可以使用矩阵减去 λ 乘以单位矩阵,然后计算该矩阵的行列式来得到特征多项式。
3、求解特征值:一般来说,特征多项式是一个高阶多项式,需要通过求解多项式的根来得到特征值。你可以使用数值计算方法(如牛顿法、二分法等)来近似地找到多项式的根,这些根即为特征值。
具体地,你可以使用线性代数库(如Eigen)来简化矩阵操作和特征值计算的过程。下面是一个使用Eigen库计算矩阵特征值的示例代码:

在上述代码中,我们使用Eigen库的EigenSolver类来计算矩阵的特征值。通过solver.eigenvalues()可以获取到计算得到的特征值。我们将结果以易读的方式打印出来。

请注意,这只是一个基本的示例代码,你可以根据具体的需求和矩阵的大小进行相应的调整。同时,对于大型矩阵或特殊类型的矩阵,可能需要使用更高效的特征值计算方法来求解特征值。但是,使用库函数能够提供更好的性能和精度,并且避免了手动实现复杂的算法。

  • img

img

【1】C++矩阵处理工具——Eigen实现:
https://blog.csdn.net/abcjennifer/article/details/7781936

#include <iostream>
#include <Eigen/Dense>
 
//using Eigen::MatrixXd;
using namespace Eigen;
using namespace Eigen::internal;
using namespace Eigen::Architecture;
 
using namespace std;
 
 
int main()
{
 
#pragma region one_d_object
 
    cout<<"*******************1D-object****************"<<endl;
 
    Vector4d v1;
    v1<< 1,2,3,4;
    cout<<"v1=\n"<<v1<<endl;
 
    VectorXd v2(3);
    v2<<1,2,3;
    cout<<"v2=\n"<<v2<<endl;
 
    Array4i v3;
    v3<<1,2,3,4;
    cout<<"v3=\n"<<v3<<endl;
 
    ArrayXf v4(3);
    v4<<1,2,3;
    cout<<"v4=\n"<<v4<<endl;
 
#pragma endregion
 
#pragma region two_d_object
    
    cout<<"*******************2D-object****************"<<endl;
    //2D objects:
    MatrixXd m(2,2);
 
    //method 1
    m(0,0) = 3;
    m(1,0) = 2.5;
    m(0,1) = -1;
    m(1,1) = m(1,0) + m(0,1);
 
    //method 2
    m<<3,-1,
        2.5,-1.5;
    cout <<"m=\n"<< m << endl;
 
#pragma endregion
 
#pragma region Comma_initializer
 
    cout<<"*******************Initialization****************"<<endl;
 
    int rows=5;
    int cols=5;
    MatrixXf m1(rows,cols);
    m1<<( Matrix3f()<<1,2,3,4,5,6,7,8,9 ).finished(),
        MatrixXf::Zero(3,cols-3),
        MatrixXf::Zero(rows-3,3),
        MatrixXf::Identity(rows-3,cols-3);
    cout<<"m1=\n"<<m1<<endl;
 
#pragma endregion
 
#pragma region Runtime_info
    
    cout<<"*******************Runtime Info****************"<<endl;
 
    MatrixXf m2(5,4);
    m2<<MatrixXf::Identity(5,4);
    cout<<"m2=\n"<<m2<<endl;
 
    MatrixXf m3;
    m3=m1*m2;
    cout<<"m3.rows()="<<m3.rows()<<"  ;  "
             <<"m3.cols()="<< m3.cols()<<endl;
    
    cout<<"m3=\n"<<m3<<endl;
 
#pragma endregion
    
#pragma region Resizing
    
    cout<<"*******************Resizing****************"<<endl;
 
    //1D-resize    
    v1.resize(4);
    cout<<"Recover v1 to 4*1 array : v1=\n"<<v1<<endl;
 
    //2D-resize
    m.resize(2,3);
    m.resize(Eigen::NoChange, 3);
    m.resizeLike(m2);
    m.resize(2,2);
    
#pragma endregion
 
#pragma region Coeff_access
    
    cout<<"*******************Coefficient access****************"<<endl;
 
    float tx=v1(1);
    tx=m1(1,1);
    cout<<endl;
 
#pragma endregion
 
#pragma  region Predefined_matrix
 
    cout<<"*******************Predefined Matrix****************"<<endl;
 
    //1D-object
    typedef  Matrix3f   FixedXD;
    FixedXD x;
    
    x=FixedXD::Zero();
    x=FixedXD::Ones();
    x=FixedXD::Constant(tx);//tx is the value
    x=FixedXD::Random();
    cout<<"x=\n"<<x<<endl;
 
    typedef ArrayXf Dynamic1D;
    //或者 typedef VectorXf Dynamic1D
    int size=3;
    Dynamic1D xx;
    xx=Dynamic1D::Zero(size);
    xx=Dynamic1D::Ones(size);
    xx=Dynamic1D::Constant(size,tx);
    xx=Dynamic1D::Random(size);
    cout<<"xx=\n"<<x<<endl;
 
    //2D-object
    typedef MatrixXf Dynamic2D;
    Dynamic2D y;
    y=Dynamic2D::Zero(rows,cols);
    y=Dynamic2D::Ones(rows,cols);
    y=Dynamic2D::Constant(rows,cols,tx);//tx is the value
    y=Dynamic2D::Random(rows,cols);
 
#pragma endregion
 
#pragma region Arithmetic_Operators
 
    cout<<"******************* Arithmetic_Operators****************"<<endl;
 
    //add & sub
    MatrixXf m4(5,4);
    MatrixXf m5;
    m4=m2+m3;
    m3-=m2;
 
    //product
    m3=m1*m2;
 
    //transposition
    m5=m4.transpose();
    //m5=m.adjoint();//伴随矩阵    
    
    //dot product
    double xtt;
    cout<<"v1=\n"<<v1<<endl;
    v2.resize(4);
    v2<<VectorXd::Ones(4);
    cout<<"v2=\n"<<v2<<endl;
 
    cout<<"*************dot product*************"<<endl;
    xtt=v1.dot(v2);
    cout<<"v1.*v2="<<xtt<<endl;
 
    //vector norm
 
    cout<<"*************matrix norm*************"<<endl;
    xtt=v1.norm();
    cout<<"norm of v1="<<xtt<<endl;
    xtt=v1.squaredNorm();
    cout<<"SquareNorm of v1="<<xtt<<endl;
 
#pragma endregion
 
cout<<endl;
}

【2】参考这个,C++求矩阵特征值:
http://blog.163.com/shikang999@126/blog/static/17262489620114208458265/@126/blog/static/17262489620114208458265/


bool Matrix_EigenValue(double *K1,int n,int LoopNumber,double Error1,double *Ret)

 {

         int i,j,k,t,m,Loop1;

         double b,c,d,g,xy,p,q,r,x,s,e,f,z,y,temp,*A;

         A=new double[n*n];

         Matrix_Hessenberg(K1,n,A);

         m=n;

         Loop1=LoopNumber;

         while(m!=0)

         {

                   t=m-1;

                   while(t>0)

                   {

                            temp=abs(A[(t-1)*n+t-1]);

                            temp+=abs(A[t*n+t]);

                            temp=temp*Error1;

                            if (abs(A[t*n+t-1])>temp)

                            {

                                     t--;

                            }

                            else

                            {

                                     break;

                            }

                   }

                   if (t==m-1)

                  {

                            Ret[(m-1)*2]=A[(m-1)*n+m-1];

                            Ret[(m-1)*2+1]=0;

                            m-=1;

                            Loop1=LoopNumber;

                   }

                   else if(t==m-2)

                   {

                            b=-A[(m-1)*n+m-1]-A[(m-2)*n+m-2];

                            c=A[(m-1)*n+m-1]*A[(m-2)*n+m-2]-A[(m-1)*n+m-2]*A[(m-2)*n+m-1];

                            d=b*b-4*c;

                            y=sqrt(abs(d));

                            if (d>0)

                            {

                                     xy=1;

                                     if (b<0)

                                     {

                                               xy=-1;

                                     }

                                     Ret[(m-1)*2]=-(b+xy*y)/2;

                                     Ret[(m-1)*2+1]=0;

                                     Ret[(m-2)*2]=c/Ret[(m-1)*2];

                                     Ret[(m-2)*2+1]=0;

                            }

                            else

                            {

                                     Ret[(m-1)*2]=-b/2;

                                     Ret[(m-2)*2]=-b/2;

                                     Ret[(m-1)*2+1]=y/2;

                                      Ret[(m-2)*2+1]=-y/2;

                            }

                            m-=2;

                            Loop1=LoopNumber;

                   }

                   else

                   {

                            if (Loop1<1)

                            {

                                     return false;

                            }

                            Loop1--;

                            j=t+2;

                            while (j<m)

                            {

                                     A[j*n+j-2]=0;

                                     j++;

                            }

                            j=t+3;

                            while (j<m)

                            {

                                     A[j*n+j-3]=0;

                                     j++;

                            }

                            k=t;

                            while (k<m-1)

                            {

                                     if (k!=t)

                                     {

                                               p=A[k*n+k-1];

                                               q=A[(k+1)*n+k-1];

                                               if (k!=m-2)

                                               {

                                                        r=A[(k+2)*n+k-1];

                                               }

                                               else

                                               {

                                                        r=0;

                                               }

                                     }

                                     else

                                     {

                                               b=A[(m-1)*n+m-1];

                                               c=A[(m-2)*n+m-2];

                                               x=b+c;

                                               y=b*c-A[(m-2)*n+m-1]*A[(m-1)*n+m-2];

                                               p=A[t*n+t]*(A[t*n+t]-x)+A[t*n+t+1]*A[(t+1)*n+t]+y;

                                               q=A[(t+1)*n+t]*(A[t*n+t]+A[(t+1)*n+t+1]-x);

                                               r=A[(t+1)*n+t]*A[(t+2)*n+t+1];

                                     }

                                     if (p!=0 || q!=0 || r!=0)

                                     {

                                               if (p<0)

                                               {

                                                        xy=-1;

                                               }

                                               else

                                               {

                                                        xy=1;

                                               }

                                               s=xy*sqrt(p*p+q*q+r*r);

                                               if (k!=t)

                                               {

                                                        A[k*n+k-1]=-s;

                                               }

                                               e=-q/s;

                                               f=-r/s;

                                               x=-p/s;

                                               y=-x-f*r/(p+s);

                                               g=e*r/(p+s);

                                               z=-x-e*q/(p+s);

                                               for (j=k;j<m;j++)

                                               {

                                                        b=A[k*n+j];

                                                        c=A[(k+1)*n+j];

                                                        p=x*b+e*c;

                                                        q=e*b+y*c;

                                                        r=f*b+g*c;

                                                        if (k!=m-2)

                                                        {

                                                                 b=A[(k+2)*n+j];

                                                                 p+=f*b;

                                                                 q+=g*b;

                                                                 r+=z*b;

                                                                 A[(k+2)*n+j]=r;

                                                        }

                                                        A[(k+1)*n+j]=q;

                                                        A[k*n+j]=p;

                                               }

                                               j=k+3;

                                               if (j>m-2)

                                               {

                                                        j=m-1;

                                               }

                                               for (i=t;i<j+1;i++)

                                               {

                                                        b=A[i*n+k];

                                                        c=A[i*n+k+1];

                                                        p=x*b+e*c;

                                                        q=e*b+y*c;

                                                        r=f*b+g*c;

                                                        if (k!=m-2)

                                                        {

                                                                 b=A[i*n+k+2];

                                                                 p+=f*b;

                                                                 q+=g*b;

                                                                 r+=z*b;

                                                                 A[i*n+k+2]=r;

                                                        }

                                                        A[i*n+k+1]=q;

                                                        A[i*n+k]=p;

                                               }

                                     }

                                     k++;

                            }

                   }

         }

         delete []A;

         return true;

}

来自GPT的回答:

要用C++求矩阵的特征值,可以借助线性代数库或数值计算库来进行计算。以下是一种常见的做法:

首先,你需要定义一个矩阵并将其填充为你的实际数据。可以使用C++的数组或向量容器来表示矩阵。
cpp
Copy code
#include <iostream>
#include <vector>

using namespace std;

// 定义矩阵类型
typedef vector<vector<double>> Matrix;

int main() {
    // 创建一个3x3的矩阵
    Matrix matrix = {{1.0, 2.0, 3.0},
                     {4.0, 5.0, 6.0},
                     {7.0, 8.0, 9.0}};
    
    // 其他代码...
    
    return 0;
}
接下来,你需要使用线性代数库或数值计算库中的函数来计算矩阵的特征值。在C++中,你可以使用一些库,如Eigen、Armadillo或GSL等,这些库提供了特征值计算的功能。
以下是使用Eigen库来计算矩阵特征值的示例:

cpp
Copy code
#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main() {
    Matrix3d matrix;
    matrix << 1.0, 2.0, 3.0,
              4.0, 5.0, 6.0,
              7.0, 8.0, 9.0;

    EigenSolver<Matrix3d> solver(matrix);
    Vector3cd eigenvalues = solver.eigenvalues();

    cout << "Eigenvalues:" << endl;
    cout << eigenvalues << endl;

    return 0;
}
在这个示例中,我们使用Eigen库创建了一个3x3的矩阵,并使用EigenSolver计算其特征值。计算结果存储在eigenvalues向量中,并通过cout打印出来。

请注意,这只是一个简单的示例,你可以根据实际情况进行调整和扩展。确保你在项目中包含了相应的库头文件,并按照库的文档使用正确的命名空间和函数来执行特征值计算。

希望这可以帮助你开始使用C++来求解矩阵的特征值。

  • 这个问题的回答你可以参考下: https://ask.csdn.net/questions/266484
  • 你也可以参考下这篇文章:C++遍历文件夹获取各文件名称并筛选指定格式类型的文件或具有特定名称的文件
  • 除此之外, 这篇博客: 设计模式之工厂模式,以C++语言举例。中的 2.1 使用虚函数与抽象类方法 部分也许能够解决你的问题, 你可以仔细阅读以下内容或跳转源博客中阅读:
  •         这样,你就可以使用接口或抽象类的指针或引用来操作产品,而无需关心具体的产品类型。  使用虚方法,你就可以在子类中重写工厂方法,以便创建不同类型的产品。

    // 定义产品接口
    class Product {
     public:
      virtual ~Product() {}
      virtual void DoSomething() = 0;
    };
    
    // 定义具体产品类
    class ConcreteProductA : public Product {
     public:
      void DoSomething() override {
        // 在这里实现具体的操作
      }
    };
    
    class ConcreteProductB : public Product {
     public:
      void DoSomething() override {
        // 在这里实现具体的操作
      }
    };
    
    // 定义工厂类
    class Factory {
     public:
      static Product* CreateProduct(int type) {
        if (type == 1) {
          return new ConcreteProductA;
        } else if (type == 2) {
          return new ConcreteProductB;
        } else {
          return nullptr;
        }
      }
    };
    
    int main() {
      // 使用接口或抽象类的指针或引用来操作产品
      Product* product = Factory::CreateProduct(1);
      product->DoSomething();
      delete product;
      return 0;
    }
    
  • 您还可以看一下 林男老师的小学生c++趣味编程入门视频教程 少儿C十十信息学奥赛竞赛网课课程中的 数组越界——捉迷藏小节, 巩固相关知识点