当前位置： 文江博客话题详情

如何用量子谐振子波函数进行数值积分？

发布于 2024-07-14 12:21:05 字数 748 浏览 10 评论 0原文

如何对无限范围内的一维积分进行数值积分（使用什么数值方法以及使用什么技巧），其中被积函数中的一个或多个函数是 1d 量子谐振子 波函数。 除此之外，我想计算谐振子基础上某个函数的矩阵元素：

phin(x) = Nn Hn(x) exp(-x2/ 2)
其中 Hn(x) 是 Hermite 多项式

Vm,n = \int_{-无穷大}^{无穷大} phim(x) V(x) phin (x) dx

同样在存在不同宽度的量子谐波波函数的情况下。

问题是波函数 phin(x) 具有振荡行为，这对于较大的 n 和来自 GSL（GNU 科学图书馆）的自适应高斯-克朗罗德求积等算法来说是一个问题）计算时间长，误差大。

原文

How to do numerical integration (what numerical method, and what tricks to use) for one-dimensional integration over infinite range, where one or more functions in the integrand are 1d quantum harmonic oscillator wave functions. Among others I want to calculate matrix elements of some function in the harmonic oscillator basis:

phin(x) = Nn Hn(x) exp(-x2/2)
where Hn(x) is Hermite polynomial

Vm,n = \int_{-infinity}^{infinity} phim(x) V(x) phin(x) dx

Also in the case where there are quantum harmonic wavefunctions with different widths.

The problem is that wavefunctions phin(x) have oscillatory behaviour, which is a problem for large n, and algorithm like adaptive Gauss-Kronrod quadrature from GSL (GNU Scientific Library) take long to calculate, and have large errors.

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（6）

奈何桥上唱咆哮 2024-07-21 12:21:05

一个不完整的答案，因为我现在时间有点短； 如果其他人不能完成图片，我可以稍后提供更多细节。

  1. 随时随地应用波函数的正交性。 这应该会显着减少计算量。

  2. 尽你所能进行分析。 升力常数、按部分分割积分等等。 隔离感兴趣的区域； 大多数波函数都是带限的，减少感兴趣的区域将大大节省工作量。

  3. 对于求积本身，您可能希望将波函数分成三部分并分别积分：中心的振荡位加上两侧的指数衰减尾部。 如果波函数是奇数，那么你很幸运，尾部会相互抵消，这意味着你只需要担心中心。 对于偶数波函数，您只需对一积分并加倍（对称性万岁！）。 否则，使用高阶高斯-拉盖尔求积规则对尾部进行积分。 你可能需要自己计算规则； 我不知道表格是否列出了好的高斯-拉盖尔规则，因为它们不经常使用。 您可能还想检查规则中节点数量增加时的错误行为； 我已经很久没有使用高斯-拉盖尔规则了，我不记得它们是否表现出龙格现象。 使用您喜欢的任何方法整合中心部分； 当然，Gauss-Kronrod 是一个可靠的选择，但还有 Fejer 求积（有时可以更好地扩展到大量节点，这可能在振荡被积函数上效果更好），甚至梯形规则（对于某些振荡函数表现出惊人的准确性） ）。 选择一个并尝试一下； 如果结果很差，请尝试另一种方法。

SO 有史以来最难的问题？ 几乎不 ：）

An incomplete answer, since I'm a little short on time at the moment; if others can't complete the picture, I can supply more details later.

  1. Apply orthogonality of the wavefunctions whenever and wherever possible. This should significantly cut down the amount of computation.

  2. Do analytically whatever you can. Lift constants, split integrals by parts, whatever. Isolate the region of interest; most wavefunctions are band-limited, and reducing the area of interest will do a lot to save work.

  3. For the quadrature itself, you probably want to split the wavefunctions into three pieces and integrate each separately: the oscillatory bit in the center plus the exponentially-decaying tails on either side. If the wavefunction is odd, you get lucky and the tails will cancel each other, meaning you only have to worry about the center. For even wavefunctions, you only have to integrate one and double it (hooray for symmetry!). Otherwise, integrate the tails using a high order Gauss-Laguerre quadrature rule. You might have to calculate the rules yourself; I don't know if tables list good Gauss-Laguerre rules, as they're not used too often. You probably also want to check the error behavior as the number of nodes in the rule goes up; it's been a long time since I used Gauss-Laguerre rules and I don't remember if they exhibit Runge's phenomenon. Integrate the center part using whatever method you like; Gauss-Kronrod is a solid choice, of course, but there's also Fejer quadrature (which sometimes scales better to high numbers of nodes, which might work nicer on an oscillatory integrand) and even the trapezoidal rule (which exhibits stunning accuracy with certain oscillatory functions). Pick one and try it out; if results are poor, give another method a shot.

Hardest question ever on SO? Hardly :)

回复 原文
像极了他 2024-07-21 12:21:05

我还建议其他一些事情：

  1. 尝试将函数转换到有限域以使集成更易于管理。
  2. 尽可能使用对称性 - 将其分解为从负无穷大到零和从零到无穷大的两个积分之和，并查看该函数是对称的还是反对称的。 它可以让你的计算变得更容易。
  3. 查看高斯-拉盖尔求积，看看它是否可以帮助您。

I'd recommend a few other things:

  1. Try transforming the function onto a finite domain to make the integration more manageable.
  2. Use symmetry where possible - break it up into the sum of two integrals from negative infinity to zero and zero to infinity and see if the function is symmetry or anti-symmetric. It could make your computing easier.
  3. Look into Gauss-Laguerre quadrature and see if it can help you.
回复 原文
荒路情人 2024-07-21 12:21:05

WKB 近似值？

The WKB approximation?

回复 原文
一江春梦 2024-07-21 12:21:05

我现在不打算解释或限定其中的任何内容。 这段代码是按原样编写的，可能不正确。 我什至不确定这是否是我正在寻找的代码，我只记得几年前我做过这个问题，在搜索我的档案时我发现了这个。 您需要自己绘制输出，提供了一些说明。 我想说的是，无限范围内的积分是我解决的一个问题，并且在执行代码时，它指出“无穷大”的舍入误差（在数字上仅意味着大）。

// 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);
        }

如果这段代码看起来正确、错误、有趣，或者您确实有具体问题，请提出，我会回答。

I am not going to explain or qualify any of this right now. This code is written as is and probably incorrect. I am not even sure if it is the code I was looking for, I just remember that years ago I did this problem and upon searching my archives I found this. You will need to plot the output yourself, some instruction is provided. I will say that the integration over infinite range is a problem that I addressed and upon execution of the code it states the round off error at 'infinity' (which numerically just means large).

// 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);
        }

If this code seems correct, wrong, interesting or you do have specific questions ask and I will answer them.

回复 原文
离旧人 2024-07-21 12:21:05

我是一名物理专业的学生，​​我也遇到了这个问题。 这些天我一直在思考这个问题并得到自己的答案。 我想它可以帮助你解决这个问题。

1.gsl中有函数可以帮助你集成振荡函数--qawo & qawf。 也许您可以设置一个值，a。 积分可以分为两部分：[0,a] 和 [a,pos_infinity]。 在第一个间隔中，您可以使用任何您想要的 gsl 集成函数，在第二个间隔中，您可以使用 qawo 或 qawf。

2.或者您可以将函数积分到上限 b，即积分在 [0,b] 中。 因此可以使用高斯传说方法来计算积分，这在 gsl 中提供。 虽然实际值和计算值之间可能存在一些差异，但如果设置正确，则差异可以忽略不计。 只要差异小于您想要的精度即可。 而这个使用gsl函数的方法只调用一次，可以使用多次，因为返回值是点及其对应的权重，而积分只是f(xi)*wi的和，更多详情可以搜索gauss legendre维基百科上的正交。 乘法和加法运算比积分快得多。

3.还有一个可以计算无穷远面积积分的函数--qagi，你可以在gsl-user's Guide中搜索到。 但是每次你需要计算积分时都会调用这个，这可能会导致一些耗时，但我不确定它在你的程序中会使用多长时间。

我建议我提供的第二个选择。

I am a student majoring in physics, and I also encountered the problem. These days I keep thinking about this question and get my own answer. I think it may help you solve this question.

1.In gsl, there are functions can help you integrate the oscillatory function--qawo & qawf. Maybe you can set a value, a. And the integration can be separated into tow parts, [0,a] and [a,pos_infinity]. In the first interval, you can use any gsl integration function you want, and in the second interval, you can use qawo or qawf.

2.Or you can integrate the function to a upper limit, b, that is integrated in [0,b]. So the integration can be calculated using a gauss legendry method, and this is provided in gsl. Although there maybe some difference between the real value and the computed value, but if you set b properly, the difference can be neglected. As long as the difference is less than the accuracy you want. And this method using the gsl function is only called once and can use many times, because the return value is point and its corresponding weight, and integration is only the sum of f(xi)*wi, for more details you can search gauss legendre quadrature on wikipedia. Multiple and addition operation is much faster than integration.

3.There is also a function which can calculate the infinity area integration--qagi, you can search it in the gsl-user's guide. But this is called everytime you need to calculate the integration, and this may cause some time consuming, but I'm not sure how long will it use in you program.

I suggest NO.2 choice I offered.

回复 原文
梦幻的味道 2024-07-21 12:21:05

如果您要使用小于 n = 100 的谐波振荡器函数，您可能需要尝试：

http ://www.mymathlib.com/quadrature/gauss_hermite.html

该程序通过具有 100 个零和权重（H_100 的零）的高斯-hermite 求积计算积分。 一旦超过 Hermite_100，积分就不那么准确了。

使用这种集成方法，我编写了一个程序来准确计算您想要计算的内容，并且它运行得相当好。 另外，可能有一种方法可以通过使用 Hermite 多项式零点的渐近形式来超越 n=100，但我还没有研究过。

If you are going to work with Harmonic oscillator functions less than n = 100 you might want to try:

http://www.mymathlib.com/quadrature/gauss_hermite.html

The program computes an integral via gauss-hermite quadrature with 100 zeroes and weights (the zeroes of H_100). Once you go over Hermite_100 the integrals are not as accurate.

Using this integration method I wrote a program calculating exactly what you want to calculate and it works fairly well. Also, there might be a way to go beyond n=100 by using the asymptotic form of the Hermite-polynomial zeroes but I haven't looked into it.

回复 原文
~没有更多了~

关于作者

陌若浮生

暂无简介

0 文章
0 评论
638 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

1CH1MKgiKxn9p

文章 0 评论 0

ゞ记忆︶ㄣ

文章 0 评论 0

JackDx

文章 0 评论 0

信远

文章 0 评论 0

yaoduoduo1995

文章 0 评论 0

霞映澄塘

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文