Implementation of the Numerov Method for the 1D square well

[evamaster]

I want to solve the Schrodinger via the Numerov Method but I had some troubles. I'm programing in C++, so here is my code:

#include<cstdlib>
#include<iostream>
#include<cmath>

using namespace std;

double x_min=-4.0 , x_max=4.0;
int N=2000;             
double r=(x_max-x_min)/(1.0*N); 
double d=2.0;
double p=0.4829;    // 2m/(hbar^2)
double Vo=20.0; // Altura del pozo

double x_m=0.1;                                     //Matching point
int i_x_m=(x_m-x_min)/r;

double Control=-123456789;

double SlopeLeft,SlopeRight;

double PAR;

double K2(double x, double E);
double NumerovL(int i, double k21, double k22, double k23, double Y[]);
double NumerovR(int i, double k21, double k22, double k23, double Y[]);
double FuncLeft(double E, double Y[]);
double FuncRight(double E, double Y[]);
void PrintFunc(double Y[]);
void Normalizar(double Y[]);
double f(double E, double Y[]);
double Biseccion(double a, double b, double Y[]);

//=========================MAIN===============================

int main(int argc, char **argv)
{  

double Y[N+1];       // Función de Onda

double paso=0.02;   // Escala   en la que se varia la energía
double Eo=0;

for(double E=0 ; E<=Vo ; E+=paso) //Cálculo de las funciones IMPARES
{
    PAR=-1;
    Eo=Biseccion(E,E+paso,Y);

    if(Eo != Control && SlopeRight*SlopeLeft<0.) 
    {
        Y[i_x_m]=FuncRight(Eo,Y);
        Y[i_x_m]=FuncLeft(Eo,Y);


        Normalizar(Y);

        PrintFunc(Y);   
    }

}


for(double E=0 ; E<=Vo ; E+=paso)   //Cálculo de las funciones PARES
{
    PAR=1;
    Eo=Biseccion(E,E+paso,Y);

    if(Eo != Control && SlopeRight*SlopeLeft>0.)
    {
        Y[i_x_m]=FuncRight(Eo,Y);
        Y[i_x_m]=FuncLeft(Eo,Y);

        Normalizar(Y);

        PrintFunc(Y);   
    }
}


  return 0;
}

 //=========================FUNCIONES===============================

 double K2(double x, double E) 
 {
  double k2;

if(fabs(x)<=d)
{
    k2=p*E;
    return k2;
}
else
{
    k2=p*(E-Vo);
    return k2;
}
}

double NumerovL(int i, double k21, double k22, double k23, double Y[]) 
{ // Para la función de Onda Izquierda
double A1,B1,C1,N;
A1=2.0*(1.0-(5.0/12.0)*r*r*k21)*Y[i-1];
B1=(1.0+(1.0/12.0)*r*r*k22)*Y[i-2];
C1=1.0+(1.0/12.0)*r*r*k23;
N=(A1-B1)/(C1);
return N;
}

double NumerovR(int i, double k21, double k22, double k23, double Y[]) 
{ // Para la función de Onda Derecha
double A1,B1,C1,N;
A1=2.0*(1.0-(5.0/12.0)*r*r*k21)*Y[i+1];
B1=(1.0+(1.0/12.0)*r*r*k22)*Y[i+2];
C1=1.0+(1.0/12.0)*r*r*k23;
N=PAR*(A1-B1)/(C1);
return N;
}

double FuncLeft(double E, double Y[])
{
double k21,k22,k23,Yleft,b;

b=sqrt(p*(Vo-E));

Y[0]=exp(b*x_min);
Y[1]=exp(b*(x_min+r));


for(int i=2 ; i<i_x_m ; i++) // Se calcula la función de Onda Izquierda
{
    k21=K2(x_min+(i-1)*r,E);
    k22=K2(x_min+(i-2)*r,E);
    k23=K2(x_min+i*r,E);

    Y[i]=NumerovL(i,k21,k22,k23,Y);

    if(i==i_x_m-1) //Función de Onda Izquierda en el Matching point
    {
        k21=K2(x_min+(i)*r,E);
        k22=K2(x_min+(i-1)*r,E);
        k23=K2(x_min+(i+1)*r,E);

        Yleft=NumerovL(i+1,k21,k22,k23,Y);
    }
}

SlopeLeft=(Yleft-Y[i_x_m-1])/r;

return Yleft;
}

double FuncRight(double E, double Y[])
{
double k21,k22,k23,Yright,b;

b=sqrt(p*(Vo-E));

Y[N]=PAR*exp(-b*(x_min+N*r));   
Y[N-1]=PAR*exp(-b*(x_min+(N-1)*r));

for(int i=N-2 ; i>i_x_m; i--) // Se calcula la función de Onda Derecha
{
    k21=K2(x_min+(i+1)*r,E);
    k22=K2(x_min+(i+2)*r,E);
    k23=K2(x_min+i*r,E);

    Y[i]=PAR*NumerovR(i,k21,k22,k23,Y);

    if(i==i_x_m+1) //Función de Onda Derecha en el Matching point
    {
        k21=K2(x_min+(i)*r,E);
        k22=K2(x_min+(i+1)*r,E);
        k23=K2(x_min+(i-1)*r,E);

        Yright=NumerovR(i-1,k21,k22,k23,Y);
    }
}

SlopeRight=PAR*(Y[i_x_m+1]-Yright)/r;

return Yright;
}

void PrintFunc(double Y[])
{
  for(int i=0 ; i<=N+1 ; i++)
 {
    cout << x_min+i*r << "\t" << Y[i] << endl;
 }
}

void Normalizar(double Y[])
{
  double S=0;

 for(int i=0 ; i<=N+1 ; i++)
 {
     S += Y[i]*Y[i]*r;
 }  

S=sqrt(S);

  for (int i=0 ; i<=N+1 ; i++)
 {
     Y[i]=Y[i]/S;
 }

}

double f(double E, double Y[])
{
double F;

F=FuncLeft(E,Y)-PAR*FuncRight(E,Y);

return F;
}

 double Biseccion(double a, double b, double Y[])
 {

  double Tol=0.00001; //Tolerancia para encontrar la raiz

 double RET=-123456789;

 if(f(a,Y)*f(b,Y)<0)
 {
     while(fabs(a-b)>Tol)
    {
         double x_m,fa,fm;

        fa=f(a,Y);
        x_m=(a+b)/2.0;
        fm=f(x_m,Y);
        //fb=f(b);

        if(fa*fm<0)
        {
            b=x_m;
            //RET=b;
        }
        else
        {
            a=x_m;
            //RET=a;
        }
    }
    RET=a;
}   

return RET;
}

The code compiles perfectly but the problem arises when I want to plot the odd solutions. I obtain two satisfactory solutions but another two that it's function is continuous but not it's derivative.

I would be very thankful if somebody can help me with this problem.

 

[maajdl]

You should give more explanation.
Without it would take too much time to guess what you did based only on your code.
I was surprised by this:

double Y[N+1]; // Función de Onda

I would have expected the wave function to be a complex number.

 

[DrClaude]

I would have expected the wave function to be a complex number.

Eigenstates of a real potential can be chosen to be purely real.

 

[evamaster]

You should give more explanation.
Without it would take too much time to guess what you did based only on your code.
I was surprised by this:

double Y[N+1]; // Función de Onda

I would have expected the wave function to be a complex number.

Basically the code takes all the energies, i.e. $0<E<Vo$ and the fuction "Biseccion" applies the Bisection algorithm between an energy $E$ and $E+step$. So the fuction finds the eigen-energy for which the left and right (from the Numerov Method wave functions matches.


Here is an example of the plot that I obtain:

As you can see, there are two graphs that are not a satisfactory solution to the problem.

Eigenstates of a real potential can be chosen to be purely real.

Yes, http://i.imgur.com/vB9B5.gif

 

[DrClaude]

It seems like the absolute value of the derivative is the same at thepoint where it is discontinuous. Check your PAR factors to see if you don't have too many/too few.

And if I may comment: try to avoid global variables. And why C++? You don't seem to be using any of its functionality.

 

[evamaster]

It seems like the absolute value of the derivative is the same at thepoint where it is discontinuous. Check your PAR factors to see if you don't have too many/too few.

And if I may comment: try to avoid global variables. And why C++? You don't seem to be using any of its functionality.

Why you recommend to avoid global variables? And I use C++ because i have to.