Implementation of the Numerov Method for the 1D square well

Click For Summary

Discussion Overview

The discussion revolves around the implementation of the Numerov Method to solve the Schrödinger equation for a one-dimensional square well potential. Participants are sharing code snippets, troubleshooting issues related to the continuity of wave functions, and discussing the properties of eigenstates in quantum mechanics.

Discussion Character

  • Technical explanation
  • Debate/contested
  • Mathematical reasoning

Main Points Raised

  • One participant describes their C++ code for solving the Schrödinger equation using the Numerov Method and notes issues with the continuity of the derivative of the wave function.
  • Another participant requests more explanation about the code, particularly expressing surprise that the wave function is represented as a real number instead of a complex number.
  • Some participants mention that eigenstates of a real potential can be purely real, suggesting that the observed behavior in the wave functions may be acceptable under certain conditions.
  • There is a suggestion to check the factors used in the code, particularly the global variable PAR, to ensure they are correctly implemented.
  • Concerns are raised about the use of global variables in the code, with a recommendation to avoid them for better programming practices.
  • One participant defends their choice of C++ for the implementation, stating it is a requirement.

Areas of Agreement / Disagreement

Participants express differing views on the representation of the wave function and the implications of its continuity properties. There is no consensus on the best approach to resolve the issues presented in the code.

Contextual Notes

Participants note that the code takes a range of energies and applies the Bisection algorithm to find eigen-energies, but the specific conditions leading to the observed discontinuities in the wave function derivatives remain unresolved.

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

Debugging a Backprop Neural Network: Issues & Solutions
  • · Replies 3 ·
Replies
3
Views
2K
C/C++ How to Implement Polynomial Operations in C++?
  • · Replies 5 ·
Replies
5
Views
3K
Java Runge-Kutta method for Simple pendulum in java
  • · Replies 1 ·
Replies
1
Views
3K
Trying to compile a source code in C
  • · Replies 5 ·
Replies
5
Views
2K
C/C++ Solve ODEs Backwards in Time with C++ Euler Method
  • · Replies 1 ·
Replies
1
Views
8K
Generating Random Numbers with the Acceptance-Rejection Method
  • · Replies 22 ·
Replies
22
Views
5K
The code work in windows system but does not work in linux
  • · Replies 2 ·
Replies
2
Views
3K
C/C++ Two Body problem with RK4 in C++
  • · Replies 13 ·
Replies
13
Views
6K
C/C++ C++ simpson's method -loop doesn't stop
  • · Replies 8 ·
Replies
8
Views
4K
Is this the correct implementation of Fermat's Primality Test?
  • · Replies 8 ·
Replies
8
Views
2K