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

Question

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

фи _п (х) = N _п Н _п (х) ехр (-х/2)
где Н _п (х) Hermite polynomial

в _{м, п} = \ Int _ {- бесконечность}^{} бесконечности фи _м (х) V (х) фи _п (х) дй

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

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

Answer 1

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

1. Применить ортогональность волновых функций всякий раз, когда это возможно. Это должно значительно сократить объем вычислений.

2. Аналитически, что бы вы ни делали. Подвижные константы, разделенные интегралы по частям, что угодно. Изолировать интересующую область; большинство волновых функций ограничены по диапазону, и сокращение области интереса будет делать многое, чтобы сохранить работу.

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

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

Answer 2

WKB аппроксимация?

Answer 3

Я не собираюсь объяснять никому прямо сейчас. Этот код написан как есть и, вероятно, неверен. Я даже не уверен, что это код, который я искал, я просто помню, что много лет назад я сделал эту проблему, и, найдя в своих архивах, я нашел это. Вам нужно будет распечатать вывод самостоятельно, предоставляется некоторая инструкция. Я скажу, что интеграция по бесконечному диапазону - это проблема, к которой я обратился, и после выполнения кода она объявляет ошибку округления на «бесконечности» (которая численно просто означает большой).

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

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

Answer 4

Я рекомендую несколько других вещей:

1. Попробуйте трансформации функции на конечной области, чтобы сделать интеграцию более управляемым.
2. Используйте симметрию, где это возможно, разложите ее на сумму двух интегралов от отрицательной бесконечности до нуля и от нуля до бесконечности и посмотрите, симметрична или антисимметрична. Это может сделать ваши вычисления проще.
3. Посмотрите на Gauss-Laguerre quadrature и посмотрите, не поможет ли он вам.

Answer 5

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

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

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

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

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

Answer 6

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

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

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

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