如何对无限范围内的一维积分进行数值积分（使用什么数值方法，以及使用什么技巧），其中被积函数中的一个或多个函数是一维量子谐振子波函数。除其他外，我想在谐振子基础上计算某些函数的矩阵元素：
phi n (x) = N n H n (x) exp(-x 2 /2)
其中 H n (x) 是Hermite 多项式V m,n = \int_{-infinity}^{infinity} phi m (x) V(x) phi n (x) dx
在存在宽度不同的量子谐波波函数的情况下也是如此。
问题是波函数 phi n (x) 具有振荡行为，这对于较大的n来说是个问题，并且来自 GSL（GNU 科学图书馆）的自适应 Gauss-Kronrod 正交等算法需要很长时间来计算，并且误差很大。
一个不完整的答案，因为我现在的时间有点短；如果其他人无法完成图片，我可以稍后提供更多详细信息。
随时随地应用波函数的正交性。这应该会显着减少计算量。
尽你所能进行分析。提升常数，按部分拆分积分，等等。隔离感兴趣的区域；大多数波函数是带限的，减少感兴趣的区域将大大节省工作。
对于正交本身，您可能希望将波函数分成三部分并分别积分：中心的振荡位加上两侧的指数衰减尾部。如果波函数是奇数，你会很幸运，尾巴会相互抵消，这意味着你只需要担心中心。对于偶数波函数，您只需要积分一个并将其加倍（万岁对称！）。否则，使用高阶 Gauss-Laguerre 求积法则对尾部进行积分。您可能必须自己计算规则；我不知道表格是否列出了好的 Gauss-Laguerre 规则，因为它们不经常使用。随着规则中节点数量的增加，您可能还想检查错误行为；很久没有使用 Gauss-Laguerre 规则了 不记得他们是否表现出龙格现象。使用您喜欢的任何方法整合中心部分；当然，Gauss-Kronrod 是一个可靠的选择，但也有 Fejer 正交（有时可以更好地缩放到大量节点，这可能在振荡被积函数上工作得更好）甚至梯形规则（它在某些振荡函数上表现出惊人的精度）。挑一个试试看；如果结果很差，请尝试另一种方法。挑一个试试看；如果结果很差，请尝试另一种方法。挑一个试试看；如果结果很差，请尝试另一种方法。
有史以来最难的问题？几乎不 ：）
我推荐一些其他的东西：
WKB近似值？
我现在不打算解释或限定其中的任何内容。此代码按原样编写，可能不正确。我什至不确定它是否是我正在寻找的代码，我只记得几年前我遇到了这个问题，在搜索我的档案时我发现了这个。您需要自己绘制输出，提供了一些说明。我会说无限范围内的积分是我解决的一个问题，并且在执行代码时，它会在“无穷大”处指出舍入误差（这在数字上只是意味着大）。
// 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 ,pos_infinity]。在第一个区间，你可以使用任何你想要的gsl积分函数，在第二个区间，你可以使用qawo或qawf。
2.或者你可以把函数积分到一个上限b，也就是积分在[0, b ]中。因此可以使用高斯传奇方法计算积分，这在 gsl 中提供。虽然实际值和计算值之间可能存在一些差异，但如果b设置得当，差异可以忽略不计。只要差异小于您想要的精度。而这个使用gsl函数的方法只调用一次，可以多次使用，因为返回值是point及其对应的权重，积分只是f(xi)*wi之和，更多细节可以搜索gauss legendre维基百科上的正交。乘法和加法运算比积分快得多。
3.还有一个计算无穷大面积积分的函数——qagi，你可以在gsl-user's guide中搜索。但是每次你需要计算积分时都会调用它，这可能会导致一些耗时，但我不确定它会在你的程序中使用多长时间。
我建议我提供的 NO.2 选择。
如果您要使用小于 n = 100 的谐波振荡器函数，您可能想尝试：
http://www.mymathlib.com/quadrature/gauss_hermite.html
该程序通过具有 100 个零和权重（H_100 的零）的 gauss-hermite 求积计算积分。一旦你超过 Hermite_100，积分就不那么准确了。
使用这种集成方法，我编写了一个程序来精确计算您想要计算的内容，并且效果很好。此外，可能有一种方法可以通过使用 Hermite 多项式零点的渐近形式来超越 n=100，但我还没有研究过。