6

Cách thực hiện tích hợp số (phương pháp số và thủ thuật nào) để tích hợp một chiều qua phạm vi vô hạn, trong đó một hoặc nhiều hàm trong tích phân là 1d quantum harmonic oscillator hàm sóng. Trong số những người khác tôi muốn để tính toán các yếu tố ma trận của một số chức năng trong dao động điều hòa cơ sở:Làm thế nào để thực hiện tích hợp số với hàm sóng dao động điều hòa lượng tử?

phi n (x) = N n H n (x) exp (-x/2)
nơi H n (x) là Hermite polynomial

V m, n = \ int _ {- vô cực}^{} vô phi 012.m (x) V (x) phi n (x) dx

Cũng trong trường hợp có các hàm số sóng điều hòa lượng tử với độ rộng khác nhau.

Vấn đề là hàm số sóng phi n (x) có hành vi dao động, mà là một vấn đề đối với lớn n, và thuật toán như adaptive Gauss-Kronrod vuông góc từ (Thư viện Khoa học GNU) GSL mất nhiều thời gian để tính toán, và có lỗi lớn.

+1

Đây có phải là câu hỏi khó nhất trên SO? – MrTelly

+4

Không chỉ đơn thuần là mối quan tâm của một miền mà hầu hết không quen thuộc, bí truyền! = Khó. – Saem

+0

Gauss-Laguerre đã được giới thiệu vào [GSL2.3] (https://www.gnu.org/software/gsl/doc/html/integration.html) – zmwang

Trả lời

8

Câu trả lời chưa hoàn chỉnh, vì tôi hơi ít thời gian vào lúc này; nếu những người khác không thể hoàn thành bức tranh, tôi có thể cung cấp thêm chi tiết sau.

  1. Áp dụng trực giao các hàm sóng bất cứ khi nào và bất cứ nơi nào có thể. Điều này sẽ làm giảm đáng kể lượng tính toán.

  2. Thực hiện phân tích mọi thứ bạn có thể. Nâng các hằng số, phân tách tích phân theo từng phần, bất cứ thứ gì. Cô lập khu vực quan tâm; hầu hết các hàm sóng đều bị hạn chế về băng tần và việc giảm diện tích quan tâm sẽ làm rất nhiều việc để tiết kiệm.

  3. Đối với chính cầu phương, bạn có thể muốn chia các hàm sóng thành ba phần và tích hợp riêng từng phần: bit dao động ở giữa cộng với đuôi phân rã theo cấp số nhân ở hai bên. Nếu hàm sóng là lẻ, bạn sẽ may mắn và đuôi sẽ hủy nhau, nghĩa là bạn chỉ phải lo lắng về trung tâm. Đối với các hàm sóng, ngay cả bạn chỉ phải tích hợp một và tăng gấp đôi nó (hooray cho đối xứng!). Nếu không, hãy tích hợp đuôi bằng cách sử dụng quy tắc cầu phương Gauss-Laguerre thứ tự cao. Bạn có thể phải tự tính toán các quy tắc; Tôi không biết nếu các bảng liệt kê các quy tắc Gauss-Laguerre tốt, vì chúng không được sử dụng quá thường xuyên. Bạn cũng có thể muốn kiểm tra hành vi lỗi vì số lượng nút trong quy tắc tăng lên; đã lâu rồi kể từ khi tôi sử dụng các quy tắc Gauss-Laguerre và tôi không nhớ liệu chúng có hiện tượng của Runge hay không. Tích hợp phần trung tâm bằng cách sử dụng bất kỳ phương pháp nào bạn thích; Gauss-Kronrod là một sự lựa chọn vững chắc, nhưng cũng có Fejer (đôi khi có vảy tốt hơn với số lượng nút cao, có thể hoạt động tốt hơn trên một tích phân dao động) và thậm chí quy tắc hình thang (thể hiện độ chính xác tuyệt đối với một số hàm dao động nhất định)). Chọn một và dùng thử; nếu kết quả kém, đưa ra một phương pháp khác một shot.

Câu hỏi khó nhất trên SO? Hardly :)

1

Giá trị xấp xỉ WKB?

0

Tôi sẽ không giải thích hoặc hội đủ điều kiện bất kỳ điều này ngay bây giờ. Mã này được viết như là và có lẽ không chính xác. Tôi thậm chí không chắc chắn nếu nó là mã tôi đang tìm kiếm, tôi chỉ nhớ rằng năm trước tôi đã làm vấn đề này và khi tìm kiếm kho lưu trữ của tôi, tôi tìm thấy điều này. Bạn sẽ cần phải tự vẽ đầu ra, một số hướng dẫn được cung cấp. Tôi sẽ nói rằng sự tích hợp trên phạm vi vô hạn là một vấn đề mà tôi đã giải quyết và khi thực hiện mã nó nói ra lỗi vòng tròn tại 'vô cùng' (mà số lượng chỉ có nghĩa là lớn).

// 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); 
     } 

Nếu mã này có vẻ đúng, sai, thú vị hoặc bạn có câu hỏi cụ thể, hãy hỏi và tôi sẽ trả lời.

4

tôi khuyên bạn nên một vài điều khác:

  1. Hãy thử chuyển đổi chức năng lên một miền hữu hạn để thực hiện việc tích hợp dễ quản lý hơn.
  2. Sử dụng đối xứng khi có thể - chia nhỏ thành tổng của hai tích phân từ vô cực âm đến 0 và 0 thành vô cùng và xem hàm có đối xứng hay đối xứng. Nó có thể làm cho máy tính của bạn dễ dàng hơn.
  3. Nhìn vào Gauss-Laguerre quadrature và xem nó có thể giúp bạn không.
0

Tôi là sinh viên chuyên ngành vật lý và tôi cũng gặp phải sự cố. Những ngày này tôi tiếp tục suy nghĩ về câu hỏi này và nhận được câu trả lời của riêng tôi. Tôi nghĩ rằng nó có thể giúp bạn giải quyết câu hỏi này.

1.Trong gsl, có các chức năng có thể giúp bạn tích hợp chức năng dao động - qawo & qawf. Có thể bạn có thể đặt giá trị, a. Và tích hợp có thể được tách thành các phần kéo, [0, a] và [a, pos_infinity]. Trong khoảng thời gian đầu tiên, bạn có thể sử dụng bất kỳ hàm tích hợp gsl nào bạn muốn và trong khoảng thứ hai, bạn có thể sử dụng qawo hoặc qawf.

2. Hoặc bạn có thể tích hợp hàm vào giới hạn trên, b, được tích hợp trong [0, b]. Vì vậy, tích hợp có thể được tính toán bằng cách sử dụng một phương pháp lồng tiếng gauss, và điều này được cung cấp trong gsl. Mặc dù có thể có một số khác biệt giữa giá trị thực và giá trị được tính toán, nhưng nếu bạn đặt b đúng cách, sự khác biệt có thể bị bỏ qua. Miễn là sự khác biệt nhỏ hơn độ chính xác bạn muốn. Và phương pháp này sử dụng hàm gsl chỉ được gọi một lần và có thể sử dụng nhiều lần, vì giá trị trả về là điểm và trọng số tương ứng của nó, và tích hợp chỉ là tổng f (xi) * wi, để biết thêm chi tiết bạn có thể tìm kiếm quadrature trên wikipedia. Thao tác nhiều và bổ sung nhanh hơn nhiều so với tích hợp.

3.Có một chức năng có thể tính tích hợp vùng vô cực - qagi, bạn có thể tìm kiếm trong hướng dẫn của người dùng gsl. Nhưng điều này được gọi là mọi lúc bạn cần tính toán tích hợp, và điều này có thể gây tốn thời gian, nhưng tôi không chắc nó sẽ sử dụng bao lâu trong chương trình của bạn.

Tôi đề xuất lựa chọn số 2 mà tôi đã cung cấp.

0

Nếu bạn đang đi để làm việc với chức năng dao động điều hòa ít hơn n = 100 bạn có thể muốn thử:

http://www.mymathlib.com/quadrature/gauss_hermite.html

Chương trình tính một tích phân qua cầu phương gauss-Hermite 100 zero và trọng lượng (các số 0 của H_100). Một khi bạn đi qua Hermite_100 tích phân không chính xác.

Sử dụng phương pháp tích hợp này Tôi đã viết một chương trình tính toán chính xác những gì bạn muốn tính và nó hoạt động khá tốt. Ngoài ra, có thể có một cách để vượt xa n = 100 bằng cách sử dụng hình thức tiệm cận của các số không đa thức Hermite nhưng tôi chưa xem xét nó.

Các vấn đề liên quan