计算物理第三次作业一维薛定谔方程定态解张楚珩 ( ) May 28, 问题描述对于定态薛定谔方程 [ d2 dx 2 + 2m (V (x) E)]ψ(x) = 0 h2 ψ(x min ) = ψ(x max ) = 0 (1) 在给定势阱的情况下, 求使这个问题

Size: px

Start display at page:

Download "计算物理第三次作业一维薛定谔方程定态解张楚珩 ( ) May 28, 问题描述对于定态薛定谔方程 [ d2 dx 2 + 2m (V (x) E)]ψ(x) = 0 h2 ψ(x min ) = ψ(x max ) = 0 (1) 在给定势阱的情况下, 求使这个问题"

深滑杜
5 years ago
Views:

1 计算物理第三次作业一维薛定谔方程定态解张楚珩 ( ) May 28, 问题描述对于定态薛定谔方程 [ d2 dx 2 + 2m (V (x) E)]ψ(x) = 0 h2 ψ(x min ) = ψ(x max ) = 0 (1) 在给定势阱的情况下, 求使这个问题有非零解的本征值 E 0 及其相应的波函数在本次作业中, 我们将要求解如下几种势阱 : V (x) = (x 1) 2 (x + 2) 2 (2) { < x < 1 有限高方势阱 V (x) = (3) 0 x < 1 或 x > 1 抛物势阱 V (x) = x 2 (4) 2 问题分析及算法设计我们利用打靶法来求解该本征值问题可以看到, 本征值问题 ( d2 dx 2 + k2 )ψ = 0 ψ(x min ) = ψ(x max ) = 0 (5) 相比于定解问题, 本征值问题多了一个待定的参数 k, 因此先猜测一个试验本征值 k, 同时任设一个非零参数 δ, 把微分方程变为初值问题 ( d2 dx 2 + k2 )ψ = 0 (6) ψ(x min ) = ψ(x max ) = 0, ψ (x min ) = ψ (x max ) = δ 对于本题要求解的定态薛定谔方程, 我们可以从 x min 出发先去积分, 可以产生一个数值解 ψ < 它在经典禁戒的区指数增长, 转过左折点之后进入经典容许的区域, 并在此区域内振荡 1

2 Figure 1: 波函数在势阱下的状态如 Figure 1 所示如果越过右节点再继续积分下去, 那么这个数值将变得不稳定因为即使在一个精确的能量本征值上也可能混入一个不想要的指数增长成分, 这将导致进入经典禁戒区的积分很可能是不准确的因此, 我们采取如下做法在每一个试探的本征值上, 由 x min 出发向前积分产生一个数值解 ψ <, 同时, 由 x max 出发向前积分产生一个数值解 ψ > 为了判断这个试探能量值是不是能量本征值, 我们选择一个接合点 x m 上比较 ψ < (x m ) 和 ψ > (x m ) 是否满足如下数学表达式 ψ < (x m ) = ψ > (x m ) ψ < dx x m = ψ > dx (7) x m ψ < (x m ) 在实际的算法中, 我们先将右边的波函数 ψ > 乘上一个系数得到一个新的波函数 ψ > (x m ) ψ < (x) x < x m ψ new (x) = ψ < (x m ) ψ > (x m ) ψ (8) >(x) x > x m 同时, 我们使用表达式 f 1 ψ(x m ) [ψ <(x m h) ψ > (x m h)] < ERROR_T OLERANCE (9) 来判定, 所试探的能量值是否为符合精度要求的本征值 3 代码下面给出了 C 语言编写的代码 : 1 // 2 // main.c 3 // Schrodinger 2

3 4 // 5 // Created by ZHANG CH on // Copyright (c) 年 2014 NJU. All rights reserved. 7 // 8 9 // 问题陈述 : 10 // 利用打靶法解一维薛定谔方程本征值与本征波函数 : 11 // d2u / dx2 + k2(x)u = 0 12 // 其中 k2(x) = 2m/h^2 [E V(x)] 13 // 14 // 二阶差分公式 : 15 // d2u / dx2 = (u_i 1 2 u_i + u_i +1)/( h ^2) = k2(x)u 16 // 推得 : u[i +1] = (2 h^2 k2(x)) u[i] u[i 1] # include <stdio.h> 19 # include <math.h> # define X_MIN ( 4.0) // 的下界 x 22 # define X_MAX (4.0) // 的上界 x 23 # define DELTA (1.0) // 边界处的导数 24 # define h (0.001) // 步长 25 # define MAXN (8000) // 数组大小应等于 (X_MAX X_MIN )/h 26 # define MIDN (4000) // 处的数组下标 x_m 即 x_m =0 27 # define COEFFIENT (1.0) // 归一化系数 2m/h^2 28 # define E_START ( 10.0) // 起始处的能量试探值需比真实本征值小 29 # define DELTA_E_START (0.1) // 最开始的时候能量试探值的增量 30 # define ERROR_TOLERANCE ( ) // 最好的评估函数的误差容忍 f double u[ MAXN ] = {0}; // 储存波函数值的数组, 即波函数 psi 33 double E = E_START ; // 试探能量值 double function_ k2 ( double E, double x); 36 double function_ f ( double E); 37 double function_ V ( double x); 38 void run (); 39 void print (); int main ( int argc, const char * argv []) 42 { 43 run (); 44 print (); 45 return 0; 46 } void run () 49 { 3

4 50 double de = DELTA_ E_ START ; 51 double f_ old = function_ f ( E); 52 double f_ new = function_ f ( E += de ); 53 while (1) { 54 if ( f_old * f_new <0) 55 { 56 de = (de /2.0); 57 } 58 f_old = f_new ; 59 f_new = function_ f ( E += de ); 60 if ( fabs ( f_new ) < ERROR_ TOLERANCE ) { 61 break ; 62 } 63 } 64 } void print () 67 { 68 FILE * fp; 69 int i; 70 fp = fopen ( output.txt, w+ ); 71 fprintf (fp, E = %lf\n, E); 72 fprintf (fp, [ ); 73 for (i=0; i< MAXN ; i++) { 74 if (i!= 0) fprintf (fp,, ); 75 fprintf (fp, %lf, u[i ]); 76 } 77 fprintf (fp, ] ); 78 fclose (fp ); 79 } double function_ f ( double E) 82 { 83 int i; 84 double Umr, Uml; 85 u[0] = 0; u[1] = DELTA ; u[ MAXN 1] = 0; u[ MAXN 2] = DELTA ; 86 // 使用二阶差分计算关于的定解问题 u 87 // 正向 88 for (i=maxn 3; i >= MIDN ; i ) { 89 u[i] =(2 h*h* function_k2 (E, X_MIN +(i*h )))* u[i +1] u[i +2]; 90 } 91 Umr = u[ MIDN ]; 92 // 逆向 93 for (i=2; i <= MIDN ; i++) { 94 u[i] = (2 h*h* function_k2 (E, X_MIN +(i*h )))* u[i 1] u[i 2]; 95 } 4

5 96 Uml = u[ MIDN ]; 97 for (i= MIDN +1; i< MAXN ; i ++) { 98 u[i] = u[i]* Uml/Umr; 99 } 100 return (1/ Uml) * (2* u[ MIDN ] u[ MIDN +1] u[ MIDN 1]); 101 } double function_ k2 ( double E, double x) 104 { 105 return COEFFIENT *( E function_ V ( x )); 106 } // 改变一下函数的值可以求出不同势垒下的本征值和本征波函数 109 double function_ V ( double x) 110 { 111 // 作业题目中给定的势垒 112 return (x 1)*(x 1)*(x +2)*( x +2); 113 // 有限高方势垒 114 // return ((x <1) && (x> 1))? 100: 0; 115 // 抛物势垒 116 // return x* x; 117 } 另外给出了 MATLAB 下的代码, 代码调用了 MATLAB 内置的四阶和五阶混合龙格库塔法 ode45 ode45 用 4 阶方法提供候选解,5 阶方法控制误差, 是一种自适应步长 ( 变步长 ) 的常微分方程数值解法, 其整体截断误差为 x 5 1 % 在文件 main.m 中 2 function main 3 global E 4 clc ; 5 n=9; 6 eold =0; 7 tol =1e 2; 8 xturn = 0.5; 9 kk =0.1; 10 xmin = 4; 11 % xmax =3; for k=1:n 14 E= eold +abs ( eold )/70; de =1; [x1,u1 ]= ode45 (@fun,[ xmin xturn ],[0 kk ]); 19 fold = u1( length (x1 ),2); 20 5

6 21 while (abs ( fold ) > tol) 22 E = E + de; 23 [x1,u1 ]= ode45 (@fun,[ xmin xturn ],[0 kk ]); 24 fnew = u1( length (x1 ),2); 25 if ( fold * fnew < 0) 26 de = (de./2); 27 end 28 fold = fnew ; 29 end 30 eold = E; 31 subplot (3,3,k); 32 plot (x1,u1 (:,1), b,( 1).*ones ( size (x1)) x1,u1 (:,1), b ); 33 axis tight ; 34 title ([ 第,num2str (k), 个本征态, E=,num2str ( eold )]); 35 end 1 % 在文件 fun.m 中 2 function yy=fun(x,psi) 3 global E 4 yy =[ psi(2);( E+((x 1)^2)*((x +2)^2))* psi (1)]; 4 计算结果 4.1 三个不同势阱下的波函数对于双峰抛物线势阱 V (x) = (x 1) 2 (x + 2) 2 可以得到相应的本征值 E 0 = , 相应的波函数如图 Figure2 对于有限高方势阱可以得到相应的本征值 E 0 = , 相应的波函数如图 Figure3 对于抛物势阱可以得到相应的本征值 E 0 = , 相应的波函数如图 Figure4 4.2 双峰抛物线下的不同本征值下的波函数对于双峰抛物线势阱 V (x) = (x 1) 2 (x + 2) 2 可以得到不同的本征值 E, 由此, 图 5 给出了程序计算得到的前九个能级的波函数 6

7 Figure 2: 双峰抛物线势阱下的波函数 Figure 3: 有限高方势阱下的波函数 7

8 Figure 4: 抛物势阱下的波函数 8

9 Figure 5: 双峰抛物线势阱下前九个能级的波函数 9

新・解きながら学ぶＣ言語

新・解きながら学ぶＣ言語 330!... 67!=... 42 "... 215 " "... 6, 77, 222 #define... 114, 194 #include... 145 %... 21 %... 21 %%... 21 %f... 26 %ld... 162 %lf... 26 %lu... 162 %o... 180 %p... 248 %s... 223, 224 %u... 162 %x... 180

计算物理第三次作业一维薛定谔方程定态解 张楚珩 ( ) May 28, 问题描述 对于定态薛定谔方程 [ d2 dx 2 + 2m (V (x) E)]ψ(x) = 0 h2 ψ(x min ) = ψ(x max ) = 0 (1) 在给定势阱的情况下, 求使这个问题

计算物理第三次作业一维薛定谔方程定态解张楚珩 ( ) May 28, 问题描述对于定态薛定谔方程 [ d2 dx 2 + 2m (V (x) E)]ψ(x) = 0 h2 ψ(x min ) = ψ(x max ) = 0 (1) 在给定势阱的情况下, 求使这个问题