Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Implementation of the Numerov Method for the 1D square well

  1. Jan 8, 2014 #1
    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:

    Code (Text):

    #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 continous but not it's derivative.

    I would be very thankful if somebody can help me with this problem.
     
  2. jcsd
  3. Jan 9, 2014 #2

    maajdl

    User Avatar
    Gold Member

    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.
     
  4. Jan 9, 2014 #3

    DrClaude

    User Avatar

    Staff: Mentor

    Eigenstates of a real potential can be chosen to be purely real.
     
  5. Jan 9, 2014 #4


    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 wich 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.

    Yes, http://i.imgur.com/vB9B5.gif
     
  6. Jan 10, 2014 #5

    DrClaude

    User Avatar

    Staff: Mentor

    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.
     
  7. Jan 11, 2014 #6
    Why you recomend to avoid global variables? And I use C++ because i have to.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Implementation of the Numerov Method for the 1D square well
  1. Program Implementation (Replies: 21)

Loading...