Cómo hacer la integración numérica con la función de onda cuántica oscilador armónico?
https://stackoverflow.com/questions/599619
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPregunta
¿Cómo hacer integración numérica (qué método numérico, y qué trucos para usar) para unidimensional integración sobre rango infinito, donde una o más funciones en el integrando son 1d oscilador armónico cuántico funciones de onda. Entre otros que quiero para calcular los elementos de matriz de alguna función en la base oscilador armónico:
phi n (x) = N n H n (x) exp (-x 2 / 2)
donde H n (x) es Hermite polinomioV 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) tener un comportamiento oscilatorio, que es un problema para la gran n , y el algoritmo como adaptativa en cuadratura de Gauss-Kronrod de GSL ( GNU Scientific Library) tomará mucho tiempo para calcular, y tienen grandes errores.
Solución
Una respuesta incompleta, ya que soy un poco corto de tiempo en el momento; Si otros no pueden completar el cuadro, puedo suministrar más detalles más adelante.
-
Aplicar ortogonalidad de las funciones de onda siempre y cuando sea posible. Esto debería reducir significativamente la cantidad de cálculos.
-
Haga analíticamente lo que pueda. Levante constantes, integrales divididas 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 van a hacer mucho para guardar el trabajo.
-
Para el propio cuadratura, es probable que desee dividir las funciones de onda en tres piezas e integrar cada uno por separado: el bit oscilatoria en el centro además de las colas de manera exponencial-descomposición de cualquier lado. Si la función de onda es impar, tienes suerte y las colas se anulan entre sí, lo que significa que sólo tiene que preocuparse por el centro. Para funciones de onda incluso, es suficiente para integrar uno y duplicarlo (hurra por la simetría!). De lo contrario, integrar las colas utilizando una regla de cuadratura de Gauss-Laguerre alto orden. Puede que tenga que calcular las reglas a sí mismo; No sé si tablas muestran buenas reglas de Gauss-Laguerre, ya que no se utilizan con mucha frecuencia. Es probable que también desee comprobar el comportamiento de error como el número de nodos en la regla sube; que ha sido un largo tiempo desde que utiliza reglas de Gauss-Laguerre y no me acuerdo si presentan el fenómeno de Runge. Integrar la parte central utilizando cualquier método que te gusta; Gauss-Kronrod es una opción sólida, por supuesto, pero también hay Fejer cuadratura (que a veces escala mejor a un gran número de nodos, que podrían funcionar mejor en un integrando oscilatoria) e incluso la regla trapezoidal (que exhibe una precisión impresionante con ciertas funciones oscilatorias ). Elegir uno y probarlo; si los resultados son pobres, dar otro método de un solo tiro.
pregunta más difícil nunca en SO? Apenas:)
Otros consejos
Me gustaría recomendar algunas otras cosas:
- Intenta transformar la función en un dominio finito para que la integración sea más manejable.
- Uso simetría siempre que sea posible - dividirla en la suma de dos integrales de menos infinito a cero y de cero a infinito y ver si la función es la simetría o anti-simétrica. Podría hacer que su cómputo más fácil.
- de Gauss-Laguerre cuadratura y ver si se le puede ayudar.
El WKB aproximación?
No voy a explicar o calificar algo de esto en este momento. Este código se escribe como es y probablemente incorrecta. Ni siquiera estoy seguro de si es el código que estaba buscando, solo recuerdo que hace años que hice este problema y buscar en mis archivos encontré esto. Usted tendrá que trazar la salida de sí mismo, se proporciona alguna instrucción. Me gustaría decir que la integración sobre el rango infinito es un problema que me ha dirigido y después de la ejecución del código Afirma el error de redondeo al 'infinito' (que numéricamente sólo 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 usted tiene preguntas específicas preguntar y voy a responder a ellos.
Soy un estudiante con especialización en física, y también me encontré con el problema. En estos días no dejo de pensar acerca de esta cuestión y obtener mi propia respuesta. Creo que puede ayudar a resolver esta cuestión.
1.In GSL, hay funciones pueden ayudar a integrar la función oscilatoria - QAWO y qawf. Tal vez se puede establecer un valor, a . Y la integración se puede separar en partes de remolque, [0, a ] y [ a , pos_infinity]. En el primer intervalo, se puede utilizar cualquier función de integración GSL desea, y en el segundo intervalo, puede utilizar QAWO o qawf.
2.O puede integrar la función a un límite superior, b , que está integrado en [0, b ]. Por lo que la integración se puede calcular utilizando un método legendry gauss, y esto se proporciona en GSL. Aunque puede haber alguna diferencia entre el valor real y el valor calculado, pero si define b correctamente, la diferencia puede ser descuidado. Siempre y cuando la diferencia es menor que la precisión que desee. Y este método utiliza la función de GSL sólo se le llama una vez y se puede utilizar en muchas ocasiones, debido a que el valor de retorno es punto y su peso correspondiente, y la integración es más que la suma de f (x) * wi, para más información, se puede buscar gauss Legendre cuadratura en la wikipedia. operación múltiple y además es mucho más rápido que la integración.
3. Hay también una función que permite calcular la integración del área del infinito - qagi, usted puede buscar en la guía de GSL-usuario. Pero esto es llamado todo lo que necesita para calcular la integración, y esto puede causar el consumo de algún tiempo, pero no estoy seguro de cuánto tiempo va a utilizar en que programa.
Sugiero elección NO.2 ofrecí.
Si se 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 se pasa de Hermite_100 las integrales no son tan precisos.
El uso de este método de integración escribí un programa para calcular exactamente lo que desea calcular y funciona bastante bien. Además, puede ser una manera de ir más allá de n = 100 mediante el uso de la forma asintótica de los ceros Hermite-polinómicas pero no han mirado en él.
