調和振動子ポテンシャルの数値解

• 調和振動子ポテンシャル
• コーディング
  • アルゴリズム
  • ソースコード
• 実行結果
• 理論値との比較

解析的に固有値が求められる調和振動子ポテンシャルで数値解を求め、 解析解と比較してみます。一次元調和振動子のシュレディンガー方程式は

です。この調和振動子の解析解はよく知られていて、次のようになっています。

なるべく式を簡単にしたいので、 =1, =1 として計算します。するとするとポテンシャルは

となります。また、 は

より となります。だから固有値は

と求まります。数値計算した結果がこの解析解の固有値とどれくらい差があるかを調べて、 数値計算の精度を求めてみます。

アルゴリズム

微分方程式を解くアルゴリズムには ルンゲ・クッタ法 を使っています。 まず適当な E を与たシュレディンガー方程式(微分方程式)を解き、 波動関数 φ(x) が x→∞ で +∞ に発散する場合は E を少し増加(E + dE)させてもう一度微分方程式を解きます。

E を増加させながら微分方程式を解くことを繰り返し、 発散が +∞ から -∞ に変わったら E を減少させます。 このとき、減少の度合はいままでの 10分の1 にします(E - dE/10)。

E を減少させながらまた微分方程式を解いていき、 発散が -∞ から +∞ に変わったら E を増加させます。 増加の度合はさらに 10分の1 にします(E + dE/100)。 この操作を繰り返し、 +∞になるEと、 -∞になるEとのはさみ撃ちで 固有関数 E を決めていきます。つぎのようなイメージです．

固有関数 E が収束していく様子です．発散の符号が変化したところで色を変えています． だんだん x→∞ で 0 に近づいていることに注目してください． ただ，誤差がたまってくるのでどうしても途中で発散してしまいます．

ソースコード

[harmonic-simu.c]

     1	/* =====================================================================
     2	   ルンゲ・クッタ法で一次元調和振動子ポテンシャルの固有値を求める
     3	   Last modified: Oct 20 2003
     4	   ===================================================================== */
     5	
     6	
     7	
     8	#include <stdio.h>
     9	#include <math.h>
    10	#include "cpgplot.h"
    11	
    12	
    13	
    14	// 物理設定
    15	double E    =  0.0;
    16	double dE   =  1.0;
    17	int parity  =  1;
    18	
    19	// きざみ幅
    20	double h = 0.001;
    21	
    22	
    23	
    24	/* --------------------------------------------------------------------- 
    25	   potential function 
    26	   --------------------------------------------------------------------- */
    27	
    28	double V(double x)  // V(x)=1/2*x^2
    29	{
    30	    return (pow(x,2)/2);
    31	}
    32	
    33	
    34	
    35	/* ---------------------------------------------------------------------
    36	   導関数の定義
    37	   --------------------------------------------------------------------- */
    38	
    39	double dphi(double x, double phi, double z)  // dphi/dx = z
    40	{
    41	  return(z);
    42	}
    43	double dz(double x, double phi, double z)    // dz/dx = (V(x)-E)phi
    44	{
    45	  return( (V(x) - E) * phi );
    46	}
    47	
    48	
    49	
    50	/* ---------------------------------------------------------------------
    51	   ルンゲ・クッタ法
    52	   --------------------------------------------------------------------- */
    53	
    54	void rungeKuttaMethod(double x0, double phi0, double z0, double h){
    55	  int i, j, n, loop, color;
    56	  double x, phi, z, k[4], q[4], phi_old[2];
    57	  loop = 15/h; color = 9;
    58	  
    59	  for(n=1 ; n<=14 ; n++){
    60	    dE = -1 * dE/10; 
    61	    x = x0; phi = phi0, z = z0;  // 波動関数初期化
    62	
    63	    while(1){
    64	      for(i=1 ; i<=loop ; i++){ 
    65		// 4次のRunge-kutta法 の連立
    66		k[0] = h * dphi(x, phi, z);
    67		q[0] = h * dz(x, phi, z);
    68		k[1] = h * dphi(x + h/2, phi + k[0]/2, z + q[0]/2);
    69		q[1] = h * dz(x + h/2, phi + k[0]/2, z + q[0]/2);
    70		k[2] = h * dphi(x + h/2, phi + k[1]/2, z + q[1]/2);
    71		q[2] = h * dz(x + h/2, phi + k[1]/2, z + q[1]/2);
    72		k[3] = h * dphi(x + h, phi + k[2], z + q[2]);
    73		q[3] = h * dz(x + h, phi + k[2], z + q[2]);
    74		phi += k[0]/6 + k[1]/3 + k[2]/3 + k[3]/6;
    75		z   += q[0]/6 + q[1]/3 + q[2]/3 + q[3]/6;
    76		x = x0 + i * h;
    77	      }
    78	
    79	      // 今回と前回のphiを保存
    80	      phi_old[0] = phi_old[1]; phi_old[1] = phi;
    81	
    82	      // 今回と前回のphiが異符号ならwhile(1)を抜ける
    83	      // i=1ではphi_oldの値が無意味なので比較しない
    84	      if ( (i>=2) && ((phi_old[0] > 0) && (phi < 0)) ) break;
    85	      if ( (i>=2) && ((phi_old[0] < 0) && (phi > 0)) ) break;
    86	
    87	      printf("%2.15lf\n", E);  // モニタ
    88	      E += dE;  // 固有値変更
    89	      x = x0; phi = phi0, z = z0;
    90	    }
    91	  }
    92	
    93	  // 求まった固有値を表示
    94	  printf("E = %2.15lf\n", E);
    95	
    96	  // 波動関数の描画
    97	  x = x0; phi = phi0, z = z0;
    98	  for(i=1 ; i<=loop ; i++){ 
    99	    k[0] = h * dphi(x, phi, z);
   100	    q[0] = h * dz(x, phi, z);
   101	    k[1] = h * dphi(x + h/2, phi + k[0]/2, z + q[0]/2);
   102	    q[1] = h * dz(x + h/2, phi + k[0]/2, z + q[0]/2);
   103	    k[2] = h * dphi(x + h/2, phi + k[1]/2, z + q[1]/2);
   104	    q[2] = h * dz(x + h/2, phi + k[1]/2, z + q[1]/2);
   105	    k[3] = h * dphi(x + h, phi + k[2], z + q[2]);
   106	    q[3] = h * dz(x + h, phi + k[2], z + q[2]);
   107	    cpgsci(color); cpgmove(x, phi); cpgdraw(x + h/6, phi + k[0]/6);
   108	    cpgsci(color); cpgmove(x + h/6, phi + k[0]/6); cpgdraw(x + h/2, phi + k[1]/2);
   109	    cpgsci(color); cpgmove(x + h/2, phi + k[1]/2); cpgdraw(x + h*5/6, phi + k[2]*5/6);
   110	    cpgsci(color); cpgmove(x + h*5/6, phi + k[2]*5/6);  cpgdraw(x + h, phi + k[3]);
   111	    phi += k[0]/6 + k[1]/3 + k[2]/3 + k[3]/6;
   112	    z   += q[0]/6 + q[1]/3 + q[2]/3 + q[3]/6;
   113	    x = x0 + i * h;
   114	  }
   115	}
   116	
   117	
   118	
   119	/* ---------------------------------------------------------------------
   120	   調和振動子ポテンシャルの描画(青線)  V(x)=1/2(x)
   121	   --------------------------------------------------------------------- */
   122	
   123	void DrawPotential(double xmax)
   124	{
   125	  double i;
   126	  cpgsci(4); cpgslw(1);
   127	  cpgmove(-xmax, pow(xmax,2)/2);
   128	  for(i = -xmax ; i < xmax ; i += 0.1){
   129	    cpgdraw( i, pow(i,2)/2 );
   130	  }
   131	}
   132	
   133	
   134	
   135	/* ---------------------------------------------------------------------
   136	   main
   137	   --------------------------------------------------------------------- */
   138	
   139	int main(void)
   140	{
   141	  double phi0, z0;
   142	  double xmax, bottom, top;
   143	
   144	  // グラフ描画設定
   145	  xmax   = 12.0;
   146	  bottom = 2.0;
   147	  top    = 2.0;
   148	
   149	  cpgopen( "/xserv" );
   150	  cpgpap( 8.0, .6 );
   151	  cpgenv( 0, xmax, -bottom, top, 0, 0 );
   152	
   153	  if(parity == -1){  // 奇関数
   154	    phi0 = 0.0; z0 = 1.0;
   155	  }
   156	  else{              // 偶関数
   157	    phi0 = 1.0; z0 = 0.0; 
   158	  }
   159	
   160	  DrawPotential(xmax);
   161	  rungeKuttaMethod(0.0, phi0, z0, h);
   162	
   163	  cpgclos();
   164	
   165	  return 0;
   166	}