Komplexe Matrix Exponential in C++

Question

Ist es tatsächlich möglich, die Matrix Exponential einer komplexen Matrix in C/C++ zu berechnen?

Ich habe es geschafft, das Produkt von zwei komplexen Matrizen mit Hilfe von Blasfunktionen aus der GNU Science Library zu nehmen. für matC = matA * matB:

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

Und ich habe es geschafft, die Exponentialmatrix einer Matrix zu erhalten, indem die undokumentierte

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

Verwendung Aber das scheint nicht komplexe Argumente zu akzeptieren.

Gibt es trotzdem etwas zu tun? Ich dachte früher, dass C++ zu allem fähig ist. Jetzt denke ich, es ist veraltet und kryptisch ...

Answer 1

"Ich dachte früher, dass C++ zu allem fähig ist" - wenn eine allgemeine Sprache komplexe Mathematik in ihrem Kern eingebaut hat, dann stimmt irgendetwas mit dieser Sprache nicht.

Für solche sehr spezifische Aufgaben gibt es eine gut akzeptierte Lösung: Bibliotheken. Entweder schreiben Sie Ihre eigenen, oder viel besser, verwenden Sie eine bereits bestehende.

Ich selbst brauche selten komplexe Matrizen in C++, dafür habe ich immer Matlab und ähnliche Tools verwendet. Jedoch könnte diese http://www.mathtools.net/C_C__/Mathematics/index.html Sie interessieren, wenn Sie Matlab kennen.

Es gibt ein paar andere Bibliotheken, die hilfreich sein könnten:

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

Answer 2

Mehrere Optionen:

1. ändern Sie den gsl_linalg_exponential_ss Code-Komplex zu akzeptieren Matrizen

2. schreiben Sie Ihre komplexen NxN Matrix als echte 2N x 2N Matrix

3. Diagonalisierung der Matrix nehmen die exponentielle der Eigenwerte, und drehen Sie die Matrix zurück in die ursprüngliche Basis

4. die komplexe Matrix Produkt, das verfügbar ist, implementieren die Exponentialmatrix entsprechend seiner Definition: exp(A) = sum_(n=0)^(n=infinity) A^n/(n!)

Sie müssen prüfen, welche Methoden geeignet sind für Ihr Problem.

C++ ist eine allgemeine Sprache. Wie oben erwähnt, müssen Sie, wenn Sie eine bestimmte Funktionalität benötigen, eine Bibliothek finden, die das kann oder selbst implementiert. Alternativ können Sie Software wie MatLab und Mathematica verwenden. Wenn das zu teuer ist, gibt es Open-Source-Alternativen, z.B. Salbei und Oktave.

Answer 3

Ich dachte auch, das Gleiche zu tun, schreiben Ihre komplexe NxN-Matrix als echte 2N x 2N Matrix ist der beste Weg, um das Problem zu lösen, dann verwenden Sie gsl_linalg_exponential_ss().

Angenommen A=Ar+i*Ai, wo A die komplexe Matrix ist und Ar und Ai sind die reellen Matrizen. Dann schreibe die neue Matrix B=[Ar Ai ;-Ai Ar] (Hier ist die Matrix in Matlab-Schreibweise geschrieben). Nun berechnen die exponentielle von B, nämlich eB=[eB1 eB2 ;eB3 eB4], dann exponential von A ist gegeben durch, eA=eB1+1i.*eB2
(Summieren der Matrizen eB1 und 1i.*eB2).

Answer 4

Ich habe einen Code geschrieben, um die Matrix-Exponentialfunktion der komplexen Matrizen mit der GSL-Funktion, gsl_linalg_exponential_ss zu berechnen (& m.matrix, & em.matrix, .01); Hier haben Sie den vollständigen Code und die Kompilierungsergebnisse. Ich habe das Ergebnis mit dem Matlab überprüft und das Ergebnis stimmt zu.

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