¿Cómo hacer la integración numérica con la función de onda del oscilador armónico cuántico?

Question

Cómo hacer integración numérica (qué método numérico, y qué trucos usar) para la integración unidimensional sobre el rango infinito, donde una o más funciones en el integrando son 1d quantum harmonic oscillator funciones de onda. Entre otros Quiero calcular elementos de matriz de alguna función en la base oscilador armónico:

phi _n (x) = N _n H _n (x) exp (-x/2)
donde H _n (x) es Hermite polynomial

V _{m, n} = \ int _ {- infinito}^{infinito} phi _m (x) V (x) phi _n (x) dx

También en el caso donde hay funciones de onda de armónicos cuánticos con anchuras diferentes.

El problema es que las funciones de onda phi _n (x) tienen un comportamiento oscilatorio, que es un problema para gran n, y el algoritmo como adaptativa en cuadratura de Gauss-Kronrod de GSL (Biblioteca GNU Scientific) toma mucho tiempo para calcular, y tiene grandes errores

Answer 1

Una respuesta incompleta, ya que estoy un poco corto de tiempo en este momento; si otros no pueden completar la imagen, puedo proporcionar más detalles más adelante.

1. Aplica la ortogonalidad de las funciones de onda siempre que sea posible. Esto debería reducir significativamente la cantidad de cálculos.

2. Analice lo que pueda. Levanta las constantes, divide las integrales por partes, lo que sea. Aislar la región de interés; la mayoría de las funciones de onda son de banda limitada, y la reducción del área de interés hará mucho para ahorrar trabajo.

3. Para la cuadratura en sí misma, es probable que desee dividir las funciones de onda en tres partes e integrarlas por separado: la parte oscilatoria en el centro más las colas exponencialmente decrecientes en cada lado. Si la función de onda es extraña, tienes suerte y las colas se cancelarán entre sí, lo que significa que solo tienes que preocuparte por el centro. Para las funciones de onda uniformes, solo tiene que integrar una y duplicarla (¡hurra por la simetría!). De lo contrario, integre las colas usando una regla de cuadratura de Gauss-Laguerre de orden superior. Puede que tenga que calcular las reglas usted mismo; No sé si las tablas enumeran las buenas reglas de Gauss-Laguerre, ya que no se usan con demasiada frecuencia. Probablemente también desee comprobar el comportamiento del error a medida que aumenta la cantidad de nodos en la regla; Ha pasado mucho tiempo desde que utilicé las reglas de Gauss-Laguerre y no recuerdo si exhiben el fenómeno de Runge. Integre la parte central utilizando el método que desee; Gauss-Kronrod es una opción sólida, por supuesto, pero también hay cuadratura Fejer (que a veces se escala mejor a un gran número de nodos, que podría funcionar mejor en un integrando oscilatorio) e incluso la regla trapezoidal (que exhibe una precisión asombrosa con ciertas funciones oscilatorias) Elige uno y pruébalo; si los resultados son pobres, dale una oportunidad a otro método.

Pregunta más difícil en SO? Apenas :)

Answer 2

¿Aproximación al WKB?

Answer 3

No voy a explicar o calificar nada de esto en este momento. Este código está escrito tal cual y probablemente incorrecto. Ni siquiera estoy seguro de si es el código que estaba buscando, solo recuerdo que hace años hice este problema y al buscar en mis archivos encontré esto. Deberá trazar la salida usted mismo, se proporciona alguna instrucción. Diré que la integración en un rango infinito es un problema que abordé y, tras la ejecución del código, indica el error de redondeo en 'infinito' (que numéricamente simplemente significa grande).

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

Si este código parece correcto, incorrecto, interesante o tiene preguntas específicas, las voy a responder y las responderé.

Answer 4

lo recomiendo algunas otras cosas:

1. Trate transformación de la función en un dominio finito para que la integración sea más manejable.
2. Utilice la simetría siempre que sea posible: divídala en la suma de dos integrales desde infinito negativo a cero y de cero a infinito y vea si la función es simetría o antisimétrica. Podría hacer su computación más fácil.
3. Mire en Gauss-Laguerre quadrature y vea si puede ayudarlo.

Answer 5

Soy un estudiante con especialización en física, y también encontré el problema. Estos días sigo pensando en esta pregunta y obtengo mi propia respuesta. Creo que puede ayudarte a resolver esta pregunta.

1.En gsl, hay funciones que pueden ayudarlo a integrar la función oscilatoria - qawo & qawf. Tal vez pueda establecer un valor, a. Y la integración puede separarse en partes de remolque, [0, a] y [a, pos_infinity]. En el primer intervalo, puede usar cualquier función de integración gsl que desee, y en el segundo intervalo, puede usar qawo o qawf.

2.O puede integrar la función a un límite superior, b , que está integrado en [0, b ]. Entonces, la integración puede calcularse usando un método de gauss legendry, y esto se proporciona en gsl. Aunque puede haber alguna diferencia entre el valor real y el valor calculado, pero si establece b correctamente, la diferencia puede omitirse. Siempre que la diferencia sea menor que la precisión que desee. Y este método que utiliza la función gsl solo se llama una vez y puede usarse muchas veces, porque el valor de retorno es el punto y su peso correspondiente, y la integración es solo la suma de f (xi) * wi, para más detalles puede buscar gauss legendre cuadratura en wikipedia. La operación múltiple y adicional es mucho más rápida que la integración.

3. También hay una función que puede calcular la integración del área de infinito: qagi, puede buscarla en la guía del usuario de gsl. Pero esto se llama cada vez que necesita calcular la integración, y esto puede consumir mucho tiempo, pero no estoy seguro de cuánto tiempo usará en su programa.

Sugiero la opción NO.2 que ofrecí.

Answer 6

Si va a trabajar con funciones oscilador armónico menor que n = 100 es posible que desee probar:

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

El programa calcula una integral a través de cuadratura de Gauss-Hermite con 100 ceros y pesos (los ceros de H_100). Una vez que revisas Hermite_100, las integrales no son tan precisas.

Usando este método de integración, escribí un programa calculando exactamente lo que quiere calcular y funciona bastante bien. Además, podría haber una manera de ir más allá de n = 100 mediante el uso de la forma asintótica de los ceros de polinomios de Hermite, pero no los he investigado.