Como fazer integração numérica com quantum wavefunction oscilador harmônico?
https://stackoverflow.com/questions/599619
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianPergunta
Como fazer integração numérica (o método numérico, e que truques para uso) para a integração de uma dimensão ao longo gama infinita, onde uma ou mais funções no integrando são quantum 1d oscilador harmônico onda funções. Entre outros Quero elementos de matriz calcular de alguma função na base de oscilador harmônico:
Phi n (x) = N n H n (x) exp (-X 2 / 2)
em que H n (x) é Hermite polinomialV m, n = \ int _ {- infinito} ^ {infinito} Phi m (x) V (x) Phi n (x) dx
Além disso, no caso onde há wavefunctions harmônicos quânticos com diferentes larguras.
O problema é que as funções de onda Phi n (x) têm um comportamento oscilatório, que é um problema de grande n , e algoritmo como adaptativa em quadratura de Gauss-Kronrod de GSL ( GNU Scientific Library) leva muito tempo para calcular, e têm grandes erros.
Solução
Uma resposta incompleta, desde que eu sou um pouco curto no tempo no momento; Se os outros não podem completar o quadro, eu posso fornecer mais detalhes mais tarde.
-
Aplicar ortogonalidade das funções de onda quando e sempre que possível. Isso deve reduzir significativamente a quantidade de computação.
-
Do analiticamente o que puder. constantes elevador, integrais divididas por partes, qualquer que seja. Isolar a região de interesse; a maioria das funções de onda são limitados-band, e reduzindo a área de interesse vai fazer muito para salvar o trabalho.
-
Para o próprio quadratura, você provavelmente vai querer dividir as funções de onda em três pedaços e integrar cada um separadamente: o bit oscilatório no centro mais as caudas exponencialmente em decomposição em ambos os lados. Se a função de onda é estranho, você tem sorte e as caudas irá cancelar o outro, ou seja, você só tem que se preocupar com o centro. Pois até wavefunctions, você só tem de integrar um e dobrá-lo (hooray para a simetria!). Caso contrário, integrar as caudas usando uma alta ordem da regra de quadratura de Gauss-Laguerre. Você pode ter que calcular as regras si mesmo; Eu não sei se a lista tabelas boa Gauss-Laguerre regras, como eles não são usados ??com muita freqüência. Você provavelmente também quer verificar o comportamento de erro como o número de nós na regra sobe; tem sido um longo tempo desde que eu usei regras Gauss-Laguerre e eu não me lembro se eles apresentam o fenômeno de Runge. Integrar a parte central usando qualquer método que você gosta; Gauss-Kronrod é uma escolha sólida, é claro, mas há também Fejer quadratura (que às vezes escala melhor a altos números de nós, o que mais agradável de trabalho poder em um integrando oscilatório) e até mesmo a regra trapezoidal (que exibe impressionante precisão com certas funções oscilatórias ). Escolha uma e experimentá-lo; Se os resultados são pobres, dar outro método um tiro.
pergunta mais difícil já no SO? Dificilmente:)
Outras dicas
Eu recomendo algumas outras coisas:
- Tente transformar a função em um domínio finito para tornar a integração mais administrável.
- Use simetria onde for possível - dividi-lo na soma de duas integrais de infinito negativo a zero e zero a infinito e ver se a função é a simetria ou anti-simétrica. Ele poderia fazer sua computação mais fácil.
- Gauss-Laguerre quadratura e ver se ele pode ajudá-lo.
O WKB aproximação?
Eu não vou explicar ou qualificar nada disso agora. Este código é escrito como é, e provavelmente incorreta. Não estou mesmo certo se é o código que eu estava procurando, eu só lembro que anos atrás eu fiz este problema e sobre buscando meus arquivos eu encontrei este. Você vai precisar para traçar a saída mesmo, alguma instrução é fornecido. Posso dizer que a integração ao longo gama infinita é um problema que eu abordados e após a execução do código afirma o erro de arredondamento em 'infinito' (que numericamente apenas 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);
}
Se este código parece correto, errado, interessante ou você tem perguntas específicas perguntar e eu vou respondê-las.
Eu sou um estudante com especialização em física, e eu também encontrou o problema. Esses dias eu fico pensando sobre esta questão e obter a minha própria resposta. Eu acho que pode ajudar a resolver esta questão.
1.In GSL, existem funções podem ajudá-lo a integrar a função de oscilação - QAWO & QAWF. Talvez você pode definir um valor, a . E a integração pode ser separado em partes de reboque, [0, a ] e [ a , pos_infinity]. No primeiro intervalo, você pode usar qualquer função de integração gsl você quer, e no segundo intervalo, você pode usar QAWO ou QAWF.
2.Ou você pode integrar a função de um limite superior, b , que é integrado em [0, b ]. Assim, a integração pode ser calculado usando um método legendry Gauss, e isso é fornecido na GSL. Embora talvez haja alguma diferença entre o valor real eo valor calculado, mas se você definir b corretamente, a diferença pode ser negligenciada. Enquanto a diferença é menor do que a precisão que você deseja. E este método usando a função GSL só é chamado uma vez e pode usar muitas vezes, porque o valor de retorno é ponto e seu peso correspondente, e integração é apenas a soma de f (xi) * wi, para mais informações você pode pesquisar Gauss Legendre quadratura na wikipedia. operação múltipla e disso é muito mais rápido do que a integração.
3.There também é uma função que pode calcular a integração zona infinity - QAGI, você pode pesquisá-lo no guia do utilizador GSL. Mas isso é chamado toda vez que você precisa para calcular a integração, e isso pode causar algum demorado, mas não tenho certeza quanto tempo vai usar em você programa.
Eu sugiro escolha NO.2 eu ofereci.
Se você estiver indo para o trabalho com funções de oscilador harmônico menos de n = 100 que você pode querer tentar:
http://www.mymathlib.com/quadrature/gauss_hermite.html
O programa calcula um integrante através de quadratura de Gauss-hermite com 100 zeros e pesos (os zeros de H_100). Uma vez que você passar por cima Hermite_100 as integrais não são tão precisos.
Usando este método de integração Eu escrevi um programa calcular exatamente o que você deseja calcular e ele funciona muito bem. Além disso, pode haver uma maneira de ir além n = 100 usando a forma assintótica dos zeros Hermite-polinomiais mas eu não olhei para ele.
