CODEWITHSUNDEEP
lapackeigenvector
mwa1
mwa1

Reputation: 151

lapack tridiag - eigenvectors wrong sign

I'm using the lapack function tridiag to find the eigenvectors of the hamiltonian in a simulation of the infinite square well, and I seem to randomly get a wrong sign (for p=4 and 6 in this case) :

red is the theory and green the simulation

enter image description here

(they're also not normalized, but I think it's normal ?)

Here is the code :

#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;


extern "C"
void dsteqr_( char * , int *, double [], double [], double [], int *, double [],
int *  ) ; 

void  tridiag( int n, double d[], double e[], double z[] ) 
/* Appel du sous-programme dsteqr de Lapack : calcul des valeurs et vecteurs
propres d'une matrice tridiagonale symmetrique   */
{
    int info ; // diagnostic
    double work[2*n-2] ; // espace de travail
    char compz='v' ; // mettre 'n' pour n'avoir que les valeurs propres

    dsteqr_( &compz, &n, d, e, z, &n, work, &info ) ; // scalaires pointeurs
    if ( info == 0 ) cout << "La diagonalisation s'est bien passée" << endl ;
    else cout << "Attention ! SSTEQR : info =" << info << endl ;
}

double V(double x){
    return 0;
}

int main() {
    int n=100;
    double d[n], e[n-1], z[n*n], x[n], v[n], L = 5, delta_x = L/n;
    ofstream of1("eigenvalues.dat"), of2("eigenvectors.dat");
    for(int i=0; i<n; i++){
        for (int j=0; j<n; j++){
            if (i == j) z[i+n*j] = 1;
            else z[i+n*j] = 0;
        }
        x[i] = (i+0.5)*delta_x - 0.5*L;
        v[i] = V(x[i]);
        d[i] = 2/pow(delta_x, 2) + v[i];
        if (i<n-2){ 
            e[i] = -1/pow(delta_x, 2);
        }
    }
    tridiag(n,d,e,z);
    for (int i=0; i<n; i++){
        of1 << d[i] << endl;
        for (int j=0; j<n; j++){
            of2 << z[i+n*j] << " ";
        }
        of2 << x[i] << endl;
    }
    of1.close();
    of2.close();
    return 0;
}

So does anyone know what is happening ?

Upvotes: 1

Views: 136

Answers (1)

francis
francis

Reputation: 9817

The magnitude of your numerical result being larger than the magnitude of the reference result is due to the fact that L=5. The normalization of LAPACK would correspond to L=1. To get the correct vectors, divide all values by L !

Regarding the sign... If you change the signs of vectors in an orthonormal basis, the result is still an orthonormal basis. Moreover if A is an Hermitian matrix and if x is a vector so that A.x=l.x where l is a scalar, then A.(-x)=l.(-x). So the sign of the vector you get is not very significant... Moreover, if the matrix features many times the same eigenvalue, the result can contain any of the orthonormal basis of the subspace associated to this eigenvalue. For instance, if the matrix A is the identity, any orthonormal basis is an acceptable result for dsteqr().

Upvotes: 1

Related Questions