/* ===================================================================== 調和振動子ポテンシャルの波動関数を規格化 Last modified: Dec 19 2003 ===================================================================== */ #include #include #include "cpgplot.h" /* 物理設定 */ static double E = 2000.0; /* 固有値 */ static double dE = -1.0; static int parity = 1; /* ルンゲ・クッタ法のきざみ幅 */ static double h = 0.001; /* --------------------------------------------------------------------- ポテンシャル関数 --------------------------------------------------------------------- */ double V(double x) { /* V(x)=x^2 */ return pow(x,2); } /* --------------------------------------------------------------------- 導関数の定義 --------------------------------------------------------------------- */ double dpsi(double x, double psi, double z) { /* \frac{d\psi}{dx} = z */ return(z); } double dz(double x, double psi, double z) { /* \frac{dz}{dx} = (V(x)-E)\psi */ return( (V(x) - E) * psi ); } /* --------------------------------------------------------------------- 固有値を求める・波動関数を規格化する --------------------------------------------------------------------- */ void ravelAndNormalize(double x0, double psi0, double z0, double h){ double x, psi, z, k[4], q[4], psi_old, psi_tmp; double is_converge_start = 4.0; /* 発散評価の開始点 */ double is_converge_value = 0.00001; /* 発散評価のしきい値 */ double converge_point; /* 発散点 */ double integral; /* 半区間での (波動関数)^2 の積分値 */ double A; /* 規格化因子 */ int decimal = 14; /* 計算する小数点以下の桁数 */ int i, j, n; int loop = decimal/h; int color = 9; /* 固有値Eの決定 */ for (n=1 ; n1) && ((psi_old > 0) && (psi < 0)) ) break; if ( (i>1) && ((psi_old < 0) && (psi > 0)) ) break; printf("E = %2.12lf\n", E); /* モニタ */ /* 固有値変更 */ E += dE; /* 波動関数初期化 */ x = x0; psi = psi0, z = z0; } } /* 規格化する前の (波動関数)^2 の描画・発散点の決定 */ x = x0; psi = psi0, z = z0; for (i=1 ; i<=loop ; i++) { k[0] = h * dpsi(x, psi, z); q[0] = h * dz (x, psi, z); k[1] = h * dpsi(x+h/2, psi+k[0]/2, z+q[0]/2); q[1] = h * dz (x+h/2, psi+k[0]/2, z+q[0]/2); k[2] = h * dpsi(x+h/2, psi+k[1]/2, z+q[1]/2); q[2] = h * dz (x+h/2, psi+k[1]/2, z+q[1]/2); k[3] = h * dpsi(x+h, psi+k[2], z+q[2]); q[3] = h * dz (x+h, psi+k[2], z+q[2]); /* PGPLOT */ cpgsci(color); cpgmove(x, pow(psi,2)); cpgdraw(x+h/6, pow(psi+k[0]/6,2)); cpgsci(color); cpgmove(x+h/6, pow(psi+k[0]/6,2)); cpgdraw(x+h/2, pow(psi+k[1]/2,2)); cpgsci(color); cpgmove(x+h/2, pow(psi+k[1]/2,2)); cpgdraw(x+h*5/6, pow(psi+k[2]*5/6,2)); cpgsci(color); cpgmove(x+h*5/6, pow(psi+k[2]*5/6,2)); cpgdraw(x+h, pow(psi+k[3],2)); psi += k[0]/6 + k[1]/3 + k[2]/3 + k[3]/6; z += q[0]/6 + q[1]/3 + q[2]/3 + q[3]/6; /* 次のきざみ幅へ */ x = x0 + i * h; /* 発散点の決定 */ if ( (x>is_converge_start) && (fabs(psi)converge_point) break; } /* 規格化因子 */ A = 1 / sqrt(2*integral); /* 規格化した (波動関数)^2 の描画 */ x = x0; psi = psi0, z = z0, color=2; for (i=1 ; i<=loop ; i++) { k[0] = h * dpsi(x, psi, z); q[0] = h * dz (x, psi, z); k[1] = h * dpsi(x+h/2, psi+k[0]/2, z+q[0]/2); q[1] = h * dz (x+h/2, psi+k[0]/2, z+q[0]/2); k[2] = h * dpsi(x+h/2, psi+k[1]/2, z+q[1]/2); q[2] = h * dz (x+h/2, psi+k[1]/2, z+q[1]/2); k[3] = h * dpsi(x+h, psi+k[2], z+q[2]); q[3] = h * dz (x+h, psi+k[2], z+q[2]); /* PGPLOT */ cpgsci(color); cpgmove(x, pow(A*psi,2)); cpgdraw(x+h/6, pow(A*(psi+k[0]/6),2)); cpgsci(color); cpgmove(x+h/6, pow(A*(psi+k[0]/6),2)); cpgdraw(x+h/2, pow(A*(psi+k[1]/2),2)); cpgsci(color); cpgmove(x+h/2, pow(A*(psi+k[1]/2),2)); cpgdraw(x+h*5/6, pow(A*(psi+k[2]*5/6),2)); cpgsci(color); cpgmove(x+h*5/6, pow(A*(psi+k[2]*5/6),2)); cpgdraw(x+h, pow(A*(psi+k[3]),2)); psi += k[0]/6 + k[1]/3 + k[2]/3 + k[3]/6; z += q[0]/6 + q[1]/3 + q[2]/3 + q[3]/6; /* 次のきざみ幅へ */ x = x0 + i * h; /* 発散点まできたら終了 */ if (x>converge_point) break; } /* 計算結果表示 */ printf(" E = %2.12lf\n 収束点 = %2.12lf\n 積分値 = %2.12lf\n |psi(0)|^2 = %2.12lf\n", E, converge_point, 2*integral, pow(A*psi0,2)); } /* --------------------------------------------------------------------- 調和振動子ポテンシャルの描画(青線) V(x)=x^2 --------------------------------------------------------------------- */ void DrawPotential(double xmax) { double i; cpgsci(4); cpgslw(1); cpgmove(-xmax, pow(xmax,2)); for (i=-xmax; i