Wie numerische Integration mit Quanten harmonischen Oszillators Wellenfunktion zu tun?
https://stackoverflow.com/questions/599619
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
Wie zu tun numerische Integration (was für numerische Methode und welche Tricks zu verwenden) für eindimensionale Integration über unendlichen Bereich, wo eine oder mehr Funktionen im Integra sind 1D Quanten harmonischen Oszillators Wellenfunktionen. Unter anderem möchte ich Matrixelemente einer Funktion in der harmonischen Oszillator Basis berechnen:
phi n (x) = N n H n (x) exp (-x 2 / 2)
wobei H n (x) ist Hermite-PolynomV m, n = \ int _ {- Unendlich} ^ {Unendlich} phi M (x) V (x) phi n (x) dx
Auch in dem Fall, wo es Quanten harmonische Wellenfunktionen mit unterschiedlichen Breiten.
Das Problem besteht darin, dass Wellenfunktionen phi n (x) hat Schwingungsverhalten, was ein Problem für großes ist n und Algorithmus wie adaptive Gauss-Kronrod Quadratur von GSL ( GNU Scientific Library) lange dauern, zu berechnen und hat große Fehler.
Lösung
Eine unvollständige Antwort, da ich ein wenig zu kurz bin pünktlich zur Zeit; wenn andere nicht das Bild abzurunden, kann ich später mehr Details liefern.
-
Übernehmen Orthogonalität der Wellenfunktionen, wann immer und wo immer möglich. Dies sollte deutlich die Menge der Berechnung abgeholzt.
-
Sie analytisch, was Sie können. Heben Sie Konstanten, geteilt Integrale von Teilen, was auch immer. Isolieren Sie die Region von Interesse; die meisten sind Wellenfunktionen bandbegrenzt und den Bereich von Interesse zu reduzieren wird eine Menge tun, um Arbeit zu speichern.
-
Für die Quadratur selbst, haben Sie wahrscheinlich die Wellenfunktionen in drei Teile geteilt werden sollen und integrieren, jeweils getrennt: die oszillierende Bit in der Mitte und den exponentiell abklingende Schwanz auf beiden Seiten. Wenn die Wellenfunktion ungerade ist, hat man Glück und der Schwanz gegenseitig aufheben, so dass Sie nur um das Zentrum zu kümmern. Für eine noch Wellenfunktionen, haben Sie nur eine integrieren und verdoppeln (hurra für Symmetrie!). Andernfalls integriert den Schwanz eine hohe Ordnung Gauß-Laguerre Quadratur Regel. Sie könnten die Regeln selbst berechnen müssen; Ich weiß nicht, ob Tabellen gute Gauss-Laguerre Regeln aufzulisten, da sie nicht zu oft verwendet. Sie wollen wahrscheinlich auch das Fehlerverhalten als die Anzahl der Knoten überprüfen in der Regel nach oben geht; es ist lange her, seit ich Gauß-Laguerre-Regeln verwendet und ich erinnere mich nicht, ob sie Runges Phänomen aufweisen. Integrieren Sie das Mittelteil mit welcher Methode Sie mögen; Gauss-Kronrod ist eine gute Wahl, natürlich, aber es gibt auch Fejer Quadratur (die manchmal zu hohen Anzahl von Knoten besser skaliert, die auf einer oszillierenden Integra funktionieren könnten schöner) und sogar die Trapezregel (die mit bestimmten Schwingungsfunktionen erstaunliche Präzision zeigt ). Wählen Sie ein und probieren Sie es aus; wenn die Ergebnisse schlecht sind, gibt eine andere Methode einen Schuss.
Hardest Frage immer auf SO? Kaum:)
Andere Tipps
Ich würde ein paar andere Dinge empfehlen:
- Versuchen Sie, die Funktion auf einen endlichen Bereich Umwandlung der Integration besser handhabbar zu machen.
- Verwenden Sie Symmetrie, wo möglich - in die Summe zweier Integrale von minus unendlich auf Null und Null bis unendlich brechen und sehen, ob die Funktion Symmetrie oder antisymmetrisch ist. Es könnte Ihre Berechnung erleichtern.
- Schauen Sie in Gauß-Laguerre Quadratur und sehen, ob es Ihnen helfen können.
Die WKB Annäherung?
Ich werde nicht erklären, oder jetzt jeden dieses Recht zu qualifizieren. Dieser Code wird geschrieben wie er ist und wahrscheinlich nicht korrekt. Ich bin nicht einmal sicher, ob es der Code ich suchte, ich erinnere mich nur, dass vor Jahren ich dieses Problem habe und meine Archive auf der Suche fand ich diese. Sie wird die Ausgabe selbst, einige Anweisung plotten müssen bereitgestellt. Ich werde sagen, dass die Integration über unendliche Reihe ist ein Problem, das ich angesprochen und bei der Ausführung des Codes besagt es die Abrundungsfehler bei ‚Unendlichkeit‘ (die numerisch nur bedeuten, groß).
// 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);
}
Wenn dieser Code scheint richtig, falsch, interessant oder Sie tun, konkrete Fragen stellen und ich werde sie beantworten.
Ich bin ein Student in Physik mit Schwerpunkt, und ich stieß auch das Problem. In diesen Tagen Ich denke immer über diese Frage und meine eigene Antwort bekommen. Ich denke, es kann Ihnen helfen, diese Frage zu lösen.
1.In GSL gibt es Funktionen können Ihnen helfen, die oszillierende Funktion integrieren - qawo & qawf. Vielleicht können Sie einen Wert gesetzt, a . Und die Integration kann in Schlepptau Teile getrennt werden, [0, a ] und [ a , pos_infinity]. Im ersten Intervall können Sie jede GSL Integrationsfunktion verwenden Sie wollen, und in dem zweiten Intervall, können Sie qawo oder qawf verwenden.
2.Or Sie die Funktion einer Obergrenze integrieren können, b , wird diese integriert in [0, b ]. So kann die Integration mit einer gauss legendry Methode berechnet werden, und dies wird in GSL zur Verfügung gestellt. Obwohl es vielleicht ein Unterschied zwischen dem tatsächlichen Wert und dem berechneten Wert, aber wenn Sie b richtig, kann der Unterschied vernachlässigt werden. Solange die Differenz kleiner ist als die Genauigkeit, die Sie wollen. Und diese Methode die GSL-Funktion wird nur einmal aufgerufen und viele Male verwenden können, weil der Rückgabewert Punkt und dem entsprechenden Gewicht, und Integration ist nur die Summe von f (xi) * wi, für weitere Informationen können Sie gauss Legendre suchen Quadratur auf Wikipedia. Multiple und Additionsoperation ist viel schneller als Integration.
3.There ist auch eine Funktion, die die Unendlichkeit Bereich Integration berechnen kann - qagi, können Sie es in der GSL-Benutzerhandbuch suchen. Aber das heißt jedes Mal, müssen Sie die Integration berechnen, und dies kann einige Zeit verursachen raubend, aber ich bin mir nicht sicher, wie lange es in Ihnen Programm verwenden.
Ich schlage vor, NO.2 Wahl, die ich angeboten.
Wenn Sie vorhaben, mit Harmonischer Oszillator Funktionen weniger arbeiten als n = 100 Sie könnten versuchen wollen:
http://www.mymathlib.com/quadrature/gauss_hermite.html
Das Programm berechnet ein integrale über Hermite-Gaußschen Quadratur mit 100 Nullen und Gewichten (die Nullstellen von H_100). Sobald Sie über Hermite_100 gehen, sind die Integrale nicht so genau.
Mit dieser Integrationsmethode ich ein Programm geschrieben Berechnung genau das, was Sie berechnen wollen und es funktioniert ziemlich gut. Auch könnte es eine Möglichkeit geben, über n = 100 unter Verwendung der asymptotische Form der Hermite-Polynom Nullen zu gehen, aber ich habe nicht in sie aussehen.
