被積分関数の 1 つまたは複数の関数が1d 量子調和振動子波動関数である場合の、無限範囲にわたる 1 次元積分の数値積分の実行方法 (数値積分法と使用するトリック) 。とりわけ、調和振動子基底でいくつかの関数の行列要素を計算したい:
phi n (x) = N n H n (x) exp(-x 2 /2)
ここで、H n (x) はエルミート多項式です。V m,n = \int_{-infinity}^{infinity} phi m (x) V(x) phi n (x) dx
また、幅の異なる量子調和波動関数が存在する場合も同様です。
問題は、波動関数 phi n (x) が振動動作をすることです。これはnが大きい場合の問題であり、GSL (GNU Scientific Library) の適応ガウスクロンロッド求積法のようなアルゴリズムは計算に時間がかかり、大きな誤差があります。
現時点では少し時間が足りないため、不完全な回答です。他の人が写真を完成できない場合は、後で詳細を提供できます。
可能な限り波動関数の直交性を適用します。これにより、計算量が大幅に削減されるはずです。
できることは何でも分析的に行います。定数を持ち上げたり、部分ごとに積分を分割したり、何でも。関心領域を分離します。ほとんどの波動関数は帯域制限されており、関心のある領域を減らすと、作業を大幅に節約できます。
直交自体については、おそらく波動関数を 3 つの部分に分割し、それぞれを個別に統合することをお勧めします。中央の振動ビットと、両側の指数関数的に減衰するテールです。波動関数が奇数の場合は、運が良く、裾が互いに打ち消し合うため、中心だけを気にする必要があります。偶数の波動関数の場合、1 つを積分して 2 倍にするだけで済みます (対称性に万歳!)。それ以外の場合は、高次のガウス・ラゲール求積規則を使用してテールを統合します。ルールを自分で計算する必要がある場合があります。あまり頻繁に使用されないため、表に適切なガウス・ラゲール規則がリストされているかどうかはわかりません。また、ルール内のノード数が増加したときのエラー動作も確認する必要があります。Gauss-Laguerre ルールを使用してから長い時間が経ちましたが、使用していません。それらがルンゲの現象を示すかどうかを思い出してください。好きな方法で中央部分を統合します。もちろん、Gauss-Kronrod は堅実な選択ですが、Fejer 求積法 (振動被積分関数でより適切に機能する可能性のある多数のノードにスケーリングすることがある) や台形則 (特定の振動関数で驚異的な精度を示す) もあります。 ）。1 つ選んで試してみてください。結果が悪い場合は、別の方法を試してください。1 つ選んで試してみてください。結果が悪い場合は、別の方法を試してください。1 つ選んで試してみてください。結果が悪い場合は、別の方法を試してください。
SOでこれまでで最も難しい質問は? しそうにない ：）
他にもいくつかお勧めします：
WKB近似？
私は物理学を専攻している学生で、私もこの問題に遭遇しました。最近、私はこの質問について考え続け、自分なりの答えを見つけています。この質問を解決するのに役立つと思います。
1.gslには、振動関数を統合するのに役立つ関数があります-qawoとqawf。たぶん、値を設定できます。また、積分は [0, a ] と [ a ,pos_infinity]の 2 つの部分に分けることができます。最初の間隔では任意の gsl 統合関数を使用でき、2 番目の間隔では qawo または qawf を使用できます。
2.または、関数を[0, b ]に統合された上限bに統合することもできます。したがって、積分はガウスの伝説的な方法を使用して計算でき、これは gsl で提供されます。実際の値と計算値には多少の違いがあるかもしれませんが、bを適切に設定すれば、その違いは無視できます。差が希望する精度よりも小さい限り。また、gsl 関数を使用するこのメソッドは 1 回だけ呼び出され、何度でも使用できます。戻り値はポイントとそれに対応する重みであり、積分は f(xi)*wi の合計のみであるため、詳細については gauss legendre を検索できます。ウィキペディアの四角形。乗算と追加の演算は、積分よりもはるかに高速です。
3.無限領域積分を計算できる関数qagiもあり、gsl-user's guideで検索できます。ただし、これは積分を計算する必要があるたびに呼び出されるため、時間がかかる場合がありますが、プログラムでどれくらいの時間使用されるかはわかりません。
私が提供したNO.2の選択肢を提案します。
現時点では、これについて説明したり、限定したりするつもりはありません。このコードはそのまま書かれており、おそらく間違っています。それが私が探していたコードであるかどうかさえわかりません。何年も前にこの問題を解決し、アーカイブを検索したときにこれを見つけたことを覚えています。自分で出力をプロットする必要があります。いくつかの指示が提供されています。無限範囲での積分は私が対処した問題であり、コードを実行すると、「無限大」での丸め誤差が示されます (数値的には単に大きいことを意味します)。
// 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);
}
このコードが正しい、間違っている、興味深いと思われる場合、または特定の質問がある場合は、それらにお答えします.