GSL комплексная матрица - собственные значения / собственные векторы - PullRequest
Новая
       7

GSL комплексная матрица - собственные значения / собственные векторы

0 голосов
/ 26 ноября 2018

Я написал программу, которая вычисляет собственные значения и собственные векторы эрмитовой матрицы.

Кто-нибудь знает, как это правильно делается в GSL?Вот то, что у меня уже есть. 

//hermitian matrix
0 1 0 -i
1 0 -i 0
0 i 0 1
i 0 1 0


#include <iostream>
#include <stdio.h>
#include <cmath>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
using namespace std;
const int N = 4;

int main(){

    gsl_eigen_hermv_workspace *workN = gsl_eigen_hermv_alloc(N);
    gsl_matrix_complex *A = gsl_matrix_complex_alloc(N, N);
    gsl_complex i = gsl_complex_rect(0.0,1.0);
    gsl_complex ii = gsl_complex_rect(0.0,-1.0);
    gsl_vector *eval = gsl_vector_alloc(N);
    gsl_matrix_complex *evec = gsl_matrix_complex_alloc(N, N);


    double mTab[] = { 
    0, 1, 0, 5, 
    1, 0, 5, 0, 
    0, 5, 0, 1,
    5, 0, 1, 0
    };

    gsl_matrix_complex_view tmpM = gsl_matrix_complex_view_array(mTab, N, N);

   gsl_matrix_complex_memcpy(A, &tmpM.matrix);
   gsl_matrix_complex_set(A, 0, 3, ii);
   gsl_matrix_complex_set(A, 1, 2, ii);
   gsl_matrix_complex_set(A, 2, 1, i);
   gsl_matrix_complex_set(A, 3, 0, i);
   gsl_eigen_hermv(A, eval, evec, workN);

   for(int i=0; i < N; i++){
       for(int j=0; j < N; j++){
           gsl_complex z = gsl_matrix_complex_get(A, i, j);
           cout << GSL_REAL(z) << "+ i" << GSL_IMAG(z) << " ";
       }
       cout << "\n";
   }

   cout << "\n";
   for(int i=0; i < N; i++){
       cout << gsl_vector_get(eval, i) << " ";
   }

   return 0;
}

Вот так я вывожу собственные векторы 

  for(int i=0; i < N; i++){
       for(int j=0; j < N; j++){
           gsl_complex z = gsl_matrix_complex_get(A, i, j);
           cout << GSL_REAL(z) << "+ i" << GSL_IMAG(z) << " ";
       }
       cout << "\n";
   }

Наконец, вот как я объявил эту матрицу. 

  double mTab[] = { 
    0, 1, 0, 5, 
    1, 0, 5, 0, 
    0, 5, 0, 1,
    5, 0, 1, 0
    };

Позже я добавил комплексные числа.

Мне удалось распечатать собственные векторы, но я не знаю, как это сделать для собственных значений.Любая помощь с этим ценится?.

1 Ответ

0 голосов
/ 26 ноября 2018

Ошибка при использовании double mTab для gsl_matrix_complex_view_array.Это предполагает массив комплексных чисел, представленных в виде действительных частей, за которыми следуют мнимые части в одном большом массиве double значений.Вы можете изменить свое определение на: 

double mTab[] = {
0, 0, 1, 0, 0, 0, 5, 0,
1, 0, 0, 0, 5, 0, 0, 0,
0, 0, 5, 0, 0, 0, 1, 0,
5, 0, 0, 0, 1, 0, 0, 0,
};

(что также означает, что вам не нужно использовать «фиктивную» переменную 5, чтобы просто переписать ее на ± i позже.) Затем код для печати собственных значенийработает хорошо.

Также у вас есть опечатка в цикле печати eigenvector: она должна быть 

gsl_complex z = gsl_matrix_complex_get(evec, i, j);

not 

gsl_complex z = gsl_matrix_complex_get(A, i, j);
...