C++の指数関数複素数行列

Question

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

私は、GNUサイエンスライブラリのblas関数を使用して2つの複雑な行列の積を取ることができました。 MATC =マタ* 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++は何でもできると思っていました。今、私は古くて暗いと思う...

Answer 1

"私は、C++は何でもできると思っていました" - 一般的な目的の言語には複雑な数学が組み込まれていれば、その言語に何か間違いがあります。

このような非常に特殊なタスクには、十分に受け入れられた解決策があります：ライブラリ。あなた自身で書くか、もっと良いものを作るか、既存のものを使うかのどちらかです。

私はC++で複雑な行列を必要とすることはめったになく、Matlabとそれに類するツールを常に使用していました。しかし、Matlabを知っているなら、このhttp://www.mathtools.net/C_C__/Mathematics/index.htmlが興味深いかもしれません。

• http://eigen.tuxfamily.org/index.php?title=Main_Page
• http://math.nist.gov/lapack++/

Answer 2

いくつかのオプション：

1. が複合体を受け入れるようにgsl_linalg_exponential_ssコードを変更

   助けになるかもしれないカップルの他のライブラリがあります。マトリックス

2. が本当の2Nのx 2N行列

3. 対角化行列としてあなたの複雑なN×N個の行列を書き、固有値の指数を取り、複雑なマトリックスを使用して、元の基礎

4. に行列を回転させますご利用可能な製品は、その者の定義に従って指数行列を実装：exp(A) = sum_(n=0)^(n=infinity) A^n/(n!)

あなたは、あなたの問題のために適切である方法で確認する必要があります。

C++は汎用言語です。上記のように、特定の機能が必要な場合は、それを実行できるライブラリを見つけるか、自分で実装する必要があります。 MatLabやMathematicaのようなソフトウェアを使うこともできます。それが高すぎると、オープンソースの代替案があります。セージとオクターブ。

Answer 3

私は同じことをやろうとしていましたが、あなたの複雑なNxNマトリックスを実際の2N×2Nマトリックスを書くことが問題を解決する最善の方法であると考えて、gsl_linalg_exponential_ss()を使用してください。

A=Ar+i*Aiと仮定、Aは複素行列とArとAiあるは実数行列です。次に、新しい行列を書くB=[Ar Ai ;-Ai Ar]（行列はmatlab表記で書かれています）。今Bの指数関数を計算し、それが（行列eB1と1i.*eB2を合計する）によって与えられる。Aの指数次いでeB=[eB1 eB2 ;eB3 eB4]、eA=eB1+1i.*eB2
あります。

Answer 4

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); 
}