Implementation of the Numerov Method for the 1D square well

Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
Sign up Log in
5 replies · 3K views
evamaster
Messages
5
Reaction score
0
I want to solve the Schrödinger via the Numerov Method but I had some troubles. I'm programing in C++, so here is my code:

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.
 
Physics news on Phys.org
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.
 
maajdl said:
I would have expected the wave function to be a complex number.
Eigenstates of a real potential can be chosen to be purely real.
 
maajdl said:
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:

pjwj1.png


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

DrClaude said:
Eigenstates of a real potential can be chosen to be purely real.

Yes, http://i.imgur.com/vB9B5.gif
 
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.
 
DrClaude said:
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.
 

Similar threads

Numerov method on a exponential grid mesh
  • Telemachus
  • · Replies 6 ·
  • Programming and Computer Science
Replies
6
Probability Density in an infinite 1D square well
  • gv3
  • · Replies 9 ·
  • Advanced Physics Homework Help
Replies
9
10 Electrons in an Infinitely Deep 1D Square Well
  • t00dles23
  • · Replies 1 ·
  • Advanced Physics Homework Help
Replies
1