Wie kann man eine numerische Integration mit der Wellenfunktion eines harmonischen Quantenoszillators durchführen?

Question

Wie zu tun ist numerische Integration (welche numerische Methode und welche Tricks) für eine eindimensionale Integration über einen unendlichen Bereich, bei der eine oder mehrere Funktionen im Integranden sind 1d quantenharmonischer Oszillator Wellenfunktionen. Unter anderem möchte ich die Matrixelemente einer Funktion auf der Basis des harmonischen Oszillators berechnen:

phi _n (x) = N _n H _n (x) exp(-x 2 /2)
_wobei H _n (x) ist Hermite-Polynom_

V _m,n \= \int_ {Unendlichkeit}^{Unendlichkeit} phi _m (x) V(x) phi _n (x) dx

Das gilt auch für den Fall, dass es quantenharmonische Wellenfunktionen mit unterschiedlicher Breite gibt.

Das Problem ist, dass die Wellenfunktionen phi _n (x) haben ein oszillierendes Verhalten, was ein Problem für große n und Algorithmen wie die adaptive Gauß-Kronrod-Quadratur von GSL (GNU Scientific Library) benötigen viel Zeit für die Berechnung und weisen große Fehler auf.

Answer 1

Eine unvollständige Antwort, da ich im Moment etwas wenig Zeit habe; wenn andere das Bild nicht vervollständigen können, kann ich später mehr Details liefern.

1. Anwendung der Orthogonalität der Wellenfunktionen, wann immer und wo immer möglich. Dadurch dürfte sich der Rechenaufwand erheblich verringern.

2. Tun Sie analytisch, was Sie können. Heben Sie Konstanten an, teilen Sie Integrale in Teile auf, was auch immer. Isolieren Sie den interessierenden Bereich; die meisten Wellenfunktionen sind bandbegrenzt, und eine Reduzierung des interessierenden Bereichs spart viel Arbeit.

3. Für die Quadratur selbst möchten Sie wahrscheinlich die Wellenfunktionen in drei Teile aufspalten und jeden separat integrieren: den oszillierenden Teil in der Mitte und die exponentiell abklingenden Schwänze auf beiden Seiten. Wenn die Wellenfunktion ungerade ist, haben Sie Glück und die Schwänze heben sich gegenseitig auf, so dass Sie sich nur um das Zentrum kümmern müssen. Bei geraden Wellenfunktionen müssen Sie nur eine integrieren und diese verdoppeln (ein Hoch auf die Symmetrie!). Andernfalls integrieren Sie die Schwänze mit einer Gauß-Laguerre-Quadraturregel hoher Ordnung. Möglicherweise müssen Sie die Regeln selbst berechnen; ich weiß nicht, ob in Tabellen gute Gauß-Laguerre-Regeln aufgeführt sind, da sie nicht allzu oft verwendet werden. Wahrscheinlich möchten Sie auch das Fehlerverhalten mit zunehmender Anzahl von Knoten in der Regel überprüfen; es ist lange her, dass ich Gauß-Laguerre-Regeln verwendet habe, und ich kann mich nicht erinnern, ob sie das Runge-Phänomen aufweisen. Integrieren Sie den Mittelteil mit einer beliebigen Methode; Gauß-Kronrod ist natürlich eine solide Wahl, aber es gibt auch die Fejer-Quadratur (die manchmal bei einer hohen Anzahl von Knoten besser skaliert, was bei einem oszillierenden Integranden besser funktionieren könnte) und sogar die Trapezregel (die bei bestimmten oszillierenden Funktionen eine erstaunliche Genauigkeit aufweist). Wählen Sie eine Methode aus und probieren Sie sie aus; wenn die Ergebnisse schlecht sind, versuchen Sie es mit einer anderen Methode.

Die schwierigste Frage überhaupt bei SO? Wohl kaum :)

Answer 2

Ich würde ein paar andere Dinge empfehlen:

1. Versuchen Sie, die Funktion in einen endlichen Bereich zu transformieren, um die Integration handhabbarer zu machen.
2. Verwenden Sie nach Möglichkeit Symmetrie - zerlegen Sie die Funktion in die Summe zweier Integrale von negativem Unendlichen bis Null und von Null bis Unendlich und prüfen Sie, ob die Funktion symmetrisch oder antisymmetrisch ist. Das könnte Ihre Berechnungen erleichtern.
3. Blick in Gauß-Laguerre-Quadratur und sehen Sie, ob es Ihnen helfen kann.

Answer 3

El WKB Annäherung?

Answer 4

Ich werde das jetzt nicht erklären oder qualifizieren. Dieser Code ist so geschrieben wie er ist und wahrscheinlich falsch. Ich bin mir nicht einmal sicher, ob es der Code ist, nach dem ich gesucht habe. Ich erinnere mich nur daran, dass ich vor Jahren dieses Problem gelöst habe, und bei der Suche in meinen Archiven fand ich diesen Code. Sie müssen die Ausgabe selbst zeichnen, eine Anleitung ist vorhanden. Ich möchte sagen, dass die Integration über einen unendlichen Bereich ein Problem ist, mit dem ich mich befasst habe, und bei der Ausführung des Codes wird der Rundungsfehler bei "unendlich" angegeben (was numerisch einfach groß bedeutet).

// 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 Ihnen dieser Code richtig, falsch oder interessant erscheint oder Sie spezielle Fragen haben, fragen Sie mich und ich werde sie beantworten.

Answer 5

Ich bin Student im Hauptfach Physik und bin auch auf dieses Problem gestoßen. In diesen Tagen denke ich immer wieder über diese Frage nach und komme zu meiner eigenen Antwort. Ich denke, sie kann Ihnen helfen, diese Frage zu lösen.

1. in gsl gibt es Funktionen, die Ihnen helfen können, die oszillierende Funktion zu integrieren - qawo und qawf. Vielleicht können Sie einen Wert einstellen, a . Und die Integration kann in zwei Teile aufgeteilt werden, [0, a ] und [ a ,pos_infinity]. Im ersten Intervall können Sie jede beliebige gsl-Integrationsfunktion verwenden, im zweiten Intervall qawo oder qawf.

2. oder Sie können die Funktion bis zu einer Obergrenze integrieren, b , die in [0, b ]. Die Integration kann also mit einer Gauß-Legenden-Methode berechnet werden, die in gsl enthalten ist. Es kann zwar ein gewisser Unterschied zwischen dem realen Wert und dem berechneten Wert bestehen, aber wenn Sie die b kann der Unterschied vernachlässigt werden. Solange die Differenz geringer ist als die gewünschte Genauigkeit. Und diese Methode mit der gsl-Funktion wird nur einmal aufgerufen und kann viele Male verwendet werden, weil der Rückgabewert Punkt und sein entsprechendes Gewicht ist, und die Integration ist nur die Summe von f(xi)*wi, für weitere Details können Sie Gauß-Legende-Quadratur auf wikipedia suchen. Multiple und Addition Betrieb ist viel schneller als Integration.

Es gibt auch eine Funktion, die die Integration der Unendlichkeitsfläche berechnen kann - Qagi, Sie können sie im gsl-Benutzerhandbuch suchen. Aber diese Funktion wird jedes Mal aufgerufen, wenn Sie die Integration berechnen müssen, und das kann einige Zeit in Anspruch nehmen, aber ich bin nicht sicher, wie lange es in Ihrem Programm verwendet wird.

Ich schlage die von mir angebotene Option Nr. 2 vor.