计算科学 - 一维方形井的 Numerov 方法的实现 - 吾爱随笔录

我想通过 Numerov 方法解决薛定谔，但我遇到了一些麻烦。我正在用 C++ 编程，所以这是我的代码：

#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;
}

基本上代码需要所有的精力，即 $0<E<Vo$ 并且函数“Biseccion”在能量之间应用 Bisection 算法 $E$ 和 $E+step$ . 因此，该函数找到左右（来自Numerov 方法）波函数匹配的特征能量。

代码编译完美，但是当我想绘制奇怪的解决方案时问题就出现了。我得到了两个令人满意的解决方案，但另外两个它的功能是连续的但不是它的导数。这是我获得的情节示例：

程序计算的特征函数

如您所见，有两个图表不能令人满意地解决问题。

如果有人可以帮助我解决这个问题，我将非常感激。