4

c / c ++で複素行列の行列指数を計算することは実際に可能ですか?

GNU Science Library の blas 関数を使用して、2 つの複雑な行列の積を得ることができました。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 に答える