Comment faire l'intégration numérique avec fonction d'onde quantique oscillateur harmonique?
https://stackoverflow.com/questions/599619
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
Comment faire intégration numérique (quelle méthode numérique, et quels trucs à utiliser) pour l'intégration d'une dimension sur la gamme infinie, où dans la integrand une ou plusieurs fonctions sont oscillateur harmonique quantique 1d onde fonctions. Entre autres, je veux calculer des éléments de la matrice d'une fonction dans la base de l'oscillateur harmonique:
phi n (x) = N n H n (x) exp (-x 2 / 2)
où H n (x) est Hermite polynômeV m, n = \ int _ {- infini} ^ {infini} phi m (x) V (x) phi n (x) dx
En outre, dans le cas où il y a des fonctions d'onde harmoniques quantiques avec des largeurs différentes.
Le problème est que les fonctions d'onde phi n (x) ont un comportement oscillatoire, ce qui est un problème pour les grandes n , et l'algorithme comme quadrature adaptatif Gauss-Kronrod de GSL ( GNU Scientific Library) prendre de temps pour calculer, et ont de grandes erreurs.
La solution
Une réponse incomplète, puisque je suis un peu à court de temps au moment; si les autres ne peuvent pas compléter le tableau, je peux fournir plus de détails plus tard.
-
Appliquer orthogonalité des fonctions d'onde chaque fois que possible. Cela devrait réduire significativement la quantité de calcul.
-
Faites analytiquement ce que vous pouvez. Soulevez constantes, par des pièces séparées Intégrales, peu importe. Isoler la région d'intérêt; la plupart sont limités dans la bande fonction d'onde, et en réduisant la zone d'intérêt va faire beaucoup pour sauver le travail.
-
Pour la quadrature elle-même, vous voulez probablement diviser les fonctions d'onde en trois morceaux et d'intégrer chacun séparément: le bit oscillatoire dans le centre ainsi que les queues de façon exponentielle pourrissant de chaque côté. Si la fonction d'onde est impair, vous avez de la chance et les queues s'annulent, ce qui signifie que vous avez seulement à vous soucier du centre. Car même fonction d'onde, il suffit d'intégrer un et doubler (hourra pour la symétrie!). Dans le cas contraire, intégrer les queues à l'aide d'un ordre élevé règle en quadrature de Gauss-Laguerre. Vous pourriez avoir à calculer les règles vous-même; Je ne sais pas si la liste des tables de bonnes règles Gauss-Laguerre, car ils ne sont pas utilisés trop souvent. Vous voulez sans doute aussi de vérifier le comportement d'erreur que le nombre de noeuds dans la règle monte; il a été longtemps que je les règles Gauss-Laguerre et je ne me souviens pas si elles présentent un phénomène de Runge. Intégrer la partie centrale à l'aide quelle que soit la méthode que vous voulez; Gauss-Kronrod est un bon choix, bien sûr, mais il y a aussi quadrature Fejer (qui adapte parfois mieux à un nombre élevé de noeuds, ce qui pourrait travailler plus agréable sur un integrand oscillatoire) et même la règle trapézoïdale (qui présente une précision étonnante avec certaines fonctions oscillatoires ). Choisissez un et l'essayer; si les résultats sont mauvais, donner une autre méthode un coup de feu.
Hardest question jamais sur le SO? A peine:)
Autres conseils
Je recommande quelques autres:
- Essayez de transformer la fonction sur un domaine fini de rendre l'intégration plus facile à gérer.
- Utilisation symétrie lorsque cela est possible - décomposer en la somme de deux intégrales de l'infini négatif à zéro et zéro à l'infini et de voir si la fonction est symétrie ou anti-symétrique. Il pourrait rendre votre calcul plus facile.
- Regardez dans Gauss-Laguerre quadrature et voir si elle peut vous aider.
Le WKB approximation?
Je ne vais pas expliquer ou de qualifier tout cela en ce moment. Ce code est écrit est probablement incorrect. Je ne suis même pas sûr que ce soit le code que je cherchais, je me souviens juste qu'il ya des années, je l'ai fait ce problème et à la recherche de mes archives, j'ai trouvé cela. Vous aurez besoin de tracer vous-même sortie, une instruction est fournie. Je dirai que l'intégration sur la gamme infinie est un problème que j'adressais et lors de l'exécution du code qu'il établit le tour de l'erreur à « l'infini » (qui signifie simplement numérique grand).
// 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 ce code semble correct, mauvais, intéressant ou vous avez des questions précises poser et je vais y répondre.
Je suis un étudiant qui se spécialise en physique, et j'ai aussi rencontré le problème. Aujourd'hui, je continue à penser à cette question et d'obtenir ma propre réponse. Je pense que cela peut vous aider à résoudre cette question.
1.In GSL, il y a des fonctions peuvent vous aider à intégrer la fonction oscillatoire - qawo & qawf. Peut-être vous pouvez définir une valeur, a . Et l'intégration peut être séparée en plusieurs parties de remorquage, [0, a ] et [ a , pos_infinity]. Dans le premier intervalle, vous pouvez utiliser une fonction d'intégration de GSL que vous voulez, et dans le second intervalle, vous pouvez utiliser qawo ou qawf.
2.Or vous pouvez intégrer la fonction à une limite supérieure, b , qui est intégré dans [0, b ]. Ainsi, l'intégration peut être calculée en utilisant une méthode legendry gauss, et cela est prévu dans GSL. Bien qu'il y ait peut-être une différence entre la valeur réelle et la valeur calculée, mais si vous définissez b correctement, la différence peut être négligée. Tant que la différence est inférieure à la précision que vous voulez. Et cette méthode en utilisant la fonction GSL est appelée une seule fois et peut utiliser plusieurs fois, parce que la valeur de retour est le point et son poids correspondant, et l'intégration est que la somme de f (xi) * wi, pour plus de détails vous pouvez rechercher gauss legendre quadrature sur wikipedia. opération multiples et plus est beaucoup plus rapide que l'intégration.
3. Il est également une fonction qui permet de calculer l'intégration de la zone de l'infini - qagi, vous pouvez rechercher dans le guide de l'utilisateur GSL. Mais cela est appelé à chaque fois que vous devez calculer l'intégration, et cela peut causer un certain temps, mais je ne sais pas combien de temps il vous utiliser dans le programme.
Je suggère le choix NO.2 je.
Si vous allez travailler avec des fonctions d'oscillateur harmonique inférieur à n = 100, vous pouvez essayer:
http://www.mymathlib.com/quadrature/gauss_hermite.html
Le programme calcule une intégrale par l'intermédiaire de quadrature de Gauss-Hermite avec 100 zéros et de poids (les zéros de H_100). Une fois que vous allez sur Hermite_100 les intégrales ne sont pas aussi précis.
En utilisant cette méthode d'intégration j'ai écrit un programme de calcul exactement ce que vous voulez calculer et il fonctionne assez bien. En outre, il pourrait y avoir un moyen d'aller au-delà n = 100 en utilisant la forme asymptotique des zéros Hermite-polynôme, mais je n'ai pas étudié la question.
