Как выполнить численное интегрирование с волновой функцией квантового гармонического осциллятора?

https://stackoverflow.com/questions/599619

11-09-2019
|

Вопрос

Как это сделать численное интегрирование (какой численный метод и какие приемы использовать) для одномерного интегрирования в бесконечном диапазоне, где одна или несколько функций в подынтегральном выражении являются 1d квантовый гармонический генератор волновые функции.Среди прочего, я хочу вычислить матричные элементы некоторой функции в базисе гармонического осциллятора:

фи_n(x) = N_n H_n(x) опыт(-x²/2)
где H_n(x) является Многочлен Эрмита

V_{m, n} = \int_{-бесконечность}^{бесконечность} phi_m(x) V(x) фи_n(x) dx

Также в случае, когда существуют квантовые гармонические волновые функции разной ширины.

Проблема в том, что волновые функции phi_n(x) имеют колебательное поведение, что является проблемой для больших n, а алгоритм, подобный адаптивной квадратуре Гаусса-Кронрода из GSL (GNU Scientific Library), требует много времени для вычисления и имеет большие ошибки.

Решение

Неполный ответ, поскольку в данный момент у меня немного не хватает времени;если другие не смогут дополнить картину, я могу предоставить более подробную информацию позже.

Применяйте ортогональность волновых функций всякий раз, когда и где это возможно.Это должно значительно сократить объем вычислений.
Делайте аналитически все, что в ваших силах.Поднимите константы, разделите интегралы на части, что угодно.Изолируйте область, представляющую интерес;большинство волновых функций ограничены полосой пропускания, и уменьшение области интереса значительно сэкономит работу.
Для самой квадратуры вы, вероятно, захотите разделить волновые функции на три части и интегрировать каждую отдельно:колеблющийся наконечник в центре плюс экспоненциально убывающие хвосты с обеих сторон.Если волновая функция нечетная, вам повезет, и хвосты будут компенсировать друг друга, а это значит, что вам нужно беспокоиться только о центре.Для четных волновых функций вам нужно только интегрировать одну и удвоить ее (ура симметрии!).В противном случае интегрируйте хвосты, используя квадратурное правило Гаусса-Лагерра высокого порядка.Возможно, вам придется рассчитать правила самостоятельно;Я не знаю, перечислены ли в таблицах хорошие правила Гаусса-Лагерра, поскольку они используются не слишком часто.Вероятно, вы также захотите проверить поведение ошибки по мере увеличения количества узлов в правиле;прошло много времени с тех пор, как я использовал правила Гаусса-Лагерра, и я не помню, проявляют ли они феномен Рунге.Интегрируйте центральную часть любым удобным вам способом;Гаусса-Кронрода, конечно, хороший выбор, но есть также квадратура Фейера (которая иногда лучше масштабируется до большого числа узлов, что может лучше работать с колебательным подынтегральным выражением) и даже правило трапецеидальности (которое демонстрирует потрясающую точность с определенными колебательными функциями).Выбери один и попробуй его;если результаты плохие, попробуйте другой метод.

Самый сложный вопрос, когда-либо задававшийся на SO?Вряд ли :)

Другие советы

Я бы порекомендовал несколько других вещей:

Попробуйте преобразовать функцию в конечную область, чтобы сделать интеграцию более управляемой.
Используйте симметрию там, где это возможно - разбейте ее на сумму двух интегралов от отрицательной бесконечности до нуля и от нуля до бесконечности и посмотрите, является ли функция симметричной или антисимметричной.Это могло бы упростить ваши вычисления.
Загляни в Квадратура Гаусса-Лагерра и посмотрим, сможет ли это вам помочь.

В ВКБ приближение?

Я не собираюсь объяснять или квалифицировать что-либо из этого прямо сейчас.Этот код написан как есть и, вероятно, неверен.Я даже не уверен, тот ли это код, который я искал, я просто помню, что много лет назад я решал эту проблему и, поискав в своих архивах, нашел это.Вам нужно будет вывести выходные данные самостоятельно, некоторые инструкции прилагаются.Я скажу, что интегрирование в бесконечном диапазоне - это проблема, которую я рассматривал, и при выполнении кода он выдает ошибку округления в 'infinity' (что численно просто означает 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);
        }

Если этот код кажется правильным, неправильным, интересным или у вас действительно есть конкретные вопросы, задавайте, и я отвечу на них.

Я студент, специализирующийся на физике, и я тоже столкнулся с этой проблемой.В эти дни я продолжаю думать над этим вопросом и получаю свой собственный ответ.Я думаю, это может помочь вам решить этот вопрос.

1.In gsl, есть функции, которые могут помочь вам интегрировать колебательную функцию - qawo и qawf.Может быть, вы сможете установить значение, a.И интеграция может быть разделена на две части, [0,a] и [a,pos_infinity].В первом интервале вы можете использовать любую функцию интеграции gsl, которую вы хотите, а во втором интервале вы можете использовать qawo или qawf.

2. Или вы можете интегрировать функцию до верхнего предела, b, который интегрирован в [0,b].Таким образом, интегрирование может быть вычислено с использованием метода Гаусса лежандри, и это предусмотрено в gsl.Хотя, возможно, существует некоторая разница между реальным значением и вычисленным значением, но если вы установите b собственно, этой разницей можно пренебречь.До тех пор, пока разница меньше точности, которую вы хотите.И этот метод, использующий функцию gsl, вызывается только один раз и может использоваться много раз, потому что возвращаемым значением является точка и соответствующий ей вес, а интегрирование - это всего лишь сумма f(xi) * wi, более подробную информацию вы можете найти в квадратуре Гаусса Лежандра в Википедии.Операция умножения и сложения выполняется намного быстрее, чем интеграция.

3. Существует также функция, которая может рассчитать интегрирование области бесконечности - qagi, вы можете найти ее в руководстве пользователя gsl.Но это вызывается каждый раз, когда вам нужно рассчитать интеграцию, и это может занять некоторое время, но я не уверен, как долго это будет использоваться в вашей программе.

Я предлагаю вариант № 2, который я предложил.

Если вы собираетесь работать с функциями гармонического осциллятора, меньшими, чем n = 100, вы можете попробовать:

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

Программа вычисляет интеграл через квадратуру Гаусса-Эрмита со 100 нулями и весами (нулями H_100).Как только вы переходите к Hermite_100, интегралы становятся не такими точными.

Используя этот метод интеграции, я написал программу, вычисляющую именно то, что вы хотите вычислить, и она работает довольно хорошо.Кроме того, возможно, есть способ выйти за пределы n = 100, используя асимптотическую форму нулей полинома Эрмита, но я не рассматривал это.

Лицензировано под: CC-BY-SA с атрибуция

Не связан с StackOverflow