4

实际上是否可以在 c / c++ 中计算复数矩阵的矩阵指数？

我已经设法使用 GNU 科学库中的 blas 函数得到两个复杂矩阵的乘积。对于 matC = matA * matB：

gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, matA, matB, GSL_COMPLEX_ZERO, matC);

我已经设法通过使用未记录的

gsl_linalg_exponential_ss(&m.matrix, &em.matrix, .01);

但这似乎不接受复杂的论点。

有没有办法做到这一点？我曾经认为 c++ 可以做任何事情。现在我认为它过时和神秘......

4

4 回答 4

3

几个选项：

  1. 修改gsl_linalg_exponential_ss代码以接受复杂的矩阵

  2. 将复杂的 NxN 矩阵写为实数 2N x 2N 矩阵

  3. 对矩阵进行对角化，取特征值的指数，然后将矩阵旋转回原来的基

  4. 使用可用的复矩阵乘积，根据其定义实现矩阵指数：exp(A) = sum_(n=0)^(n=infinity) A^n/(n!)

您必须检查哪些方法适合您的问题。

C++ 是一种通用语言。如上所述，如果您需要特定功能，则必须找到可以执行此操作的库或自己实现它。或者，您可以使用 MatLab 和 Mathematica 等软件。如果这太贵了，还有开源替代品，例如 Sage 和 Octave。

于 2012-04-12T22:19:54.773 回答
1

“我曾经认为 c++ 可以做任何事情”——如果通用语言的核心内置了复杂的数学，那么该语言就有问题。

对于这些非常具体的任务，有一个广为接受的解决方案：库。要么自己写，要么更好，使用已经存在的。

我自己很少需要 C++ 中的复杂矩阵，我总是使用 Matlab 和类似的工具。但是，如果您了解 Matlab，您可能会对这个http://www.mathtools.net/C_C__/Mathematics/index.html感兴趣。

还有几个其他库可能会有所帮助：

于 2012-04-06T21:04:23.690 回答
0

我也在考虑这样做，将复杂的 NxN 矩阵写为真正的 2N x 2N 矩阵是解决问题的最佳方法，然后使用gsl_linalg_exponential_ss().

假设A=Ar+i*Ai，其中A是复矩阵，Ar和Ai是实矩阵。然后写出新的矩阵B=[Ar Ai ;-Ai Ar]（这里的矩阵是用matlab符号写的）。现在计算 的指数B，即 eB=[eB1 eB2 ;eB3 eB4]，然后 的指数由A下式给出（对矩阵和 求和）。eA=eB1+1i.*eB2
eB11i.*eB2

于 2013-02-04T17:48:37.183 回答
0

我编写了一个代码来使用 gsl 函数 gsl_linalg_exponential_ss(&m.matrix, &em.matrix, .01); 计算复矩阵的矩阵指数。这里有完整的代码和编译结果。我用 Matlab 检查了结果，结果一致。

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
void my_gsl_complex_matrix_exponential(gsl_matrix_complex *eA, gsl_matrix_complex *A, int dimx)
{
    int j,k=0;
    gsl_complex temp;
    gsl_matrix *matreal =gsl_matrix_alloc(2*dimx,2*dimx);
    gsl_matrix *expmatreal =gsl_matrix_alloc(2*dimx,2*dimx);
    //Converting the complex matrix into real one using A=[Areal, Aimag;-Aimag,Areal]
    for (j = 0; j < dimx;j++)
        for (k = 0; k < dimx;k++)
        {
            temp=gsl_matrix_complex_get(A,j,k);
            gsl_matrix_set(matreal,j,k,GSL_REAL(temp));
            gsl_matrix_set(matreal,dimx+j,dimx+k,GSL_REAL(temp));
            gsl_matrix_set(matreal,j,dimx+k,GSL_IMAG(temp));
            gsl_matrix_set(matreal,dimx+j,k,-GSL_IMAG(temp));
        }

    gsl_linalg_exponential_ss(matreal,expmatreal,.01);

    double realp;
    double imagp;
    for (j = 0; j < dimx;j++)
        for (k = 0; k < dimx;k++)
        {
            realp=gsl_matrix_get(expmatreal,j,k);
            imagp=gsl_matrix_get(expmatreal,j,dimx+k);
            gsl_matrix_complex_set(eA,j,k,gsl_complex_rect(realp,imagp));
        }
    gsl_matrix_free(matreal);
    gsl_matrix_free(expmatreal);
}

int main()
{
    int dimx=4;
    int i, j ;
    gsl_matrix_complex *A = gsl_matrix_complex_alloc (dimx, dimx);
    gsl_matrix_complex *eA = gsl_matrix_complex_alloc (dimx, dimx);

    for (i = 0; i < dimx;i++)
    {
        for (j = 0; j < dimx;j++)
        {
            gsl_matrix_complex_set(A,i,j,gsl_complex_rect(i+j,i-j));
            if ((i-j)>=0)
            printf("%d+%di ",i+j,i-j);
            else
            printf("%d%di  ",i+j,i-j);

        }
        printf(";\n");
    }
    my_gsl_complex_matrix_exponential(eA,A,dimx);
    printf("\n Printing the complex matrix exponential\n");
    gsl_complex compnum;
    for (i = 0; i < dimx;i++)
    {
        for (j = 0; j < dimx;j++)
        {
            compnum=gsl_matrix_complex_get(eA,i,j);
            if (GSL_IMAG(compnum)>=0)
                printf("%f+%fi\t ",GSL_REAL(compnum),GSL_IMAG(compnum));
            else
                printf("%f%fi\t ",GSL_REAL(compnum),GSL_IMAG(compnum));
        }
        printf("\n");
    }
    return(0);
}
于 2013-02-05T22:49:34.923 回答