phin(x) = Nn Hn(x) exp(-x2/2)
其中 Hn(x) 是 Hermite polynomialVm,n = \int_{-infinity}^{infinity} phim(x) V(x) phin(x) dx
当存在具有不同宽度的量子谐波函数时也适用。
问题在于波函数 phin(x) 具有振荡行为，这对于较大的 n 来说是一个问题，而像 GSL（GNU Scientific Library）中的自适应 Gauss-Kronrod 积分算法需要很长时间才能计算，并且误差很大。尽可能地应用波函数的正交性。这应该会显著减少计算量。
尽可能地进行分析处理。提取常数，按部就班地拆分积分。隔离感兴趣的区域;大多数波函数是带限的，减少感兴趣的区域将大大节省工作量。
对于积分本身，您可能需要将波函数分成三个部分分别进行积分：中心的振荡部分加上两侧指数衰减的尾巴。如果波函数是奇数，则运气好，尾巴将互相抵消，意味着您只需要担心中心部分。对于偶数波函数，您只需要积分一个并将其加倍（对称性万岁！）。否则，使用高阶Gauss-Laguerre积分规则积分尾巴。您可能需要自己计算规则；我不知道表格是否列出了良好的Gauss-Laguerre规则，因为它们并不经常使用。您还可能想检查规则节点数增加时的误差行为；我很久没有使用Gauss-Laguerre规则了，不记得它们是否表现出Runge现象。使用任何您喜欢的方法积分中心部分；当然，Gauss-Kronrod是一个可靠的选择，但还有Fejer积分（有时对高节点数缩放更好，这可能在振荡积分因子上效果更好），甚至梯形规则（对于某些振荡函数表现出惊人的精度）。选择一种并尝试一下；如果结果不好，再尝试另一种方法。
SO上最难的问题吗？根本不是 ：）
我现在不打算解释或证明任何东西。这段代码是按原样编写的，可能是错误的。我甚至不确定它是否是我正在寻找的代码，我只记得多年前我做过这个问题，并在搜索我的档案时找到了这个。你需要自己绘制输出，提供了一些指令。我会说，无限范围内的积分是一个我解决的问题，在执行代码时它会在“无穷大”处报告舍入误差（数值上意味着大）。
// compile g++ base.cc -lm
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <math.h>
using namespace std;
int main ()
{
double xmax,dfx,dx,x,hbar,k,dE,E,E_0,m,psi_0,psi_1,psi_2;
double w,num;
int n,temp,parity,order;
double last;
double propogator(double E,int parity);
double eigen(double E,int parity);
double f(double x, double psi, double dpsi);
double g(double x, double psi, double dpsi);
double rk4(double x, double psi, double dpsi, double E);
ofstream datas ("test.dat");
E_0= 1.602189*pow(10.0,-19.0);// ev joules conversion
dE=E_0*.001;
//w^2=k/m v=1/2 k x^2 V=??? = E_0/xmax x^2 k-->
//w=sqrt( (2*E_0)/(m*xmax) );
//E=(0+.5)*hbar*w;
cout << "Enter what energy level your looking for, as an (0,1,2...) INTEGER: ";
cin >> order;
E=0;
for (n=0; n<=order; n++)
{
parity=0;
//if its even parity is 1 (true)
temp=n;
if ( (n%2)==0 ) {parity=1; }
cout << "Energy " << n << " has these parameters: ";
E=eigen(E,parity);
if (n==order)
{
propogator(E,parity);
cout <<" The postive values of the wave function were written to sho.dat \n";
cout <<" In order to plot the data should be reflected about the y-axis \n";
cout <<" evenly for even energy levels and oddly for odd energy levels\n";
}
E=E+dE;
}
}
double propogator(double E,int parity)
{
ofstream datas ("sho.dat") ;
double hbar =1.054*pow(10.0,-34.0);
double m =9.109534*pow(10.0,-31.0);
double E_0= 1.602189*pow(10.0,-19.0);
double dx =pow(10.0,-10);
double xmax= 100*pow(10.0,-10.0)+dx;
double dE=E_0*.001;
double last=1;
double x=dx;
double psi_2=0.0;
double psi_0=0.0;
double psi_1=1.0;
// cout <<parity << " parity passsed \n";
psi_0=0.0;
psi_1=1.0;
if (parity==1)
{
psi_0=1.0;
psi_1=m*(1.0/(hbar*hbar))* dx*dx*(0-E)+1 ;
}
do
{
datas << x << "\t" << psi_0 << "\n";
psi_2=(2.0*m*(dx/hbar)*(dx/hbar)*(E_0*(x/xmax)*(x/xmax)-E)+2.0)*psi_1-psi_0;
//cout << psi_1 << "=psi_1\n";
psi_0=psi_1;
psi_1=psi_2;
x=x+dx;
} while ( x<= xmax);
//I return 666 as a dummy value sometimes to check the function has run
return 666;
}
double eigen(double E,int parity)
{
double hbar =1.054*pow(10.0,-34.0);
double m =9.109534*pow(10.0,-31.0);
double E_0= 1.602189*pow(10.0,-19.0);
double dx =pow(10.0,-10);
double xmax= 100*pow(10.0,-10.0)+dx;
double dE=E_0*.001;
double last=1;
double x=dx;
double psi_2=0.0;
double psi_0=0.0;
double psi_1=1.0;
do
{
psi_0=0.0;
psi_1=1.0;
if (parity==1)
{double psi_0=1.0; double psi_1=m*(1.0/(hbar*hbar))* dx*dx*(0-E)+1 ;}
x=dx;
do
{
psi_2=(2.0*m*(dx/hbar)*(dx/hbar)*(E_0*(x/xmax)*(x/xmax)-E)+2.0)*psi_1-psi_0;
psi_0=psi_1;
psi_1=psi_2;
x=x+dx;
} while ( x<= xmax);
if ( sqrt(psi_2*psi_2)<=1.0*pow(10.0,-3.0))
{
cout << E << " is an eigen energy and " << psi_2 << " is psi of 'infinity' \n";
return E;
}
else
{
if ( (last >0.0 && psi_2<0.0) ||( psi_2>0.0 && last<0.0) )
{
E=E-dE;
dE=dE/10.0;
}
}
last=psi_2;
E=E+dE;
} while (E<=E_0);
}
我是一名物理专业的学生，也遇到了这个问题。这些天我一直在思考这个问题，并得出了自己的答案。我认为它可能会帮助你解决这个问题。
1. 在 gsl 中，有一些函数可以帮助你积分振荡函数——qawo 和 qawf。也许你可以设置一个值 a。然后将积分分成两部分，[0，a] 和 [a，正无穷]。在第一个区间，你可以使用任何 gsl 积分函数，而在第二个区间，你可以使用 qawo 或 qawf。
2. 或者你可以将函数集成到上限 b 中，即在 [0，b] 中进行积分。因此，可以使用高斯勒让德方法计算积分，这在 gsl 中提供。虽然实际值和计算值之间可能存在一些差异，但如果您适当设置了 b，则可以忽略差异。只要差异小于所需的精度即可。而且，使用 gsl 函数的这种方法仅被调用一次，并且可以多次使用，因为返回值是点及其相应的权重，积分仅是 f(xi)*wi 的总和。有关更多详细信息，您可以在维基百科上搜索高斯勒让德积分。乘法和加法运算比积分快得多。
3. 还有一个可以计算无穷区域积分的函数——qagi，您可以在 gsl 用户指南中搜索它。但是，每次需要计算积分时都会调用此函数，这可能会耗费一些时间，但我不确定在您的程序中需要多长时间。
我建议选择第二个选项。
http://www.mymathlib.com/quadrature/gauss_hermite.html
该程序通过高斯-埃尔米特积分法计算积分，使用了100个零点和权重（H_100的零点）。一旦超过Hermite_100，积分精度就不如之前那么准确。