第五章薛定谔方程数值解法1

1、第五章第五章 薛定谔方程数值解法薛定谔方程数值解法量子力学的基本方程是量子力学的基本方程是薛定谔方程薛定谔方程 (5.0.1) 其中其中 为普朗克常数；为普朗克常数；H为粒子的为粒子的哈密顿量；哈密顿量； 为波函数，用来描述粒子的微观运为波函数，用来描述粒子的微观运动状态，一般是空间位置和时间的函数。动状态，一般是空间位置和时间的函数。iHthh,2/nH可表示为可表示为UmH222(5.0.2) 其中其中 是拉普斯算符，是拉普斯算符，U是势能。薛定谔方是势能。薛定谔方程是一个偏微分方程，我们要同时求程是一个偏微分方程，我们要同时求本征值本征值和和本征波函数本征波函数。现在我们考虑几种特殊情况

2、。现在我们考虑几种特殊情况下，对薛定谔方程进行计算机求解。下，对薛定谔方程进行计算机求解。25.1 一维方势阱的计算机求解一维方势阱的计算机求解n考虑一维空间的粒子运动，它的势能考虑一维空间的粒子运动，它的势能U具有具有如下性质如下性质WXXWXVU, 0000(5.1.1) 如图如图5.1.1所示所示:V0A=0X=W图图5.1.1 方势阱方势阱 n由于势能和时间无关，属于定态问题。考虑一由于势能和时间无关，属于定态问题。考虑一维问题，故薛定谔方程维问题，故薛定谔方程(5.0.1)式可简化。为此式可简化。为此设设( ) ( )x f t(5.1.2) 代入代入(5.0.1式，用分离变数法，可

3、得式，用分离变数法，可得Efdtdfi和和222( )2dU xEm dx(5.1.3) (5.1.4) 其中其中E为波函数的本征值。为波函数的本征值。n由由(5.1.3)式，直接可得式，直接可得/( )iEtf tce其中其中c为任意常数。粒子的波函数为任意常数。粒子的波函数 可表示成可表示成/( , )( )iEtx tx e(5.1.5) 这就是这就是定态波函数定态波函数，其中常数，其中常数c已经包括在已经包括在 中。中。 ( )xn几率密度为几率密度为22|( , )|( )|x tx与时间无关。与时间无关。222( ( )2dU xEm dx(5.1.6) 上式与时间无关，称上式与时

4、间无关，称定态薛定谔方程定态薛定谔方程。所以，求解薛定谔方程所以，求解薛定谔方程(5.0.1)式变为求解式变为求解(5.1.4)式，式，即即 n为简单起见，我们令为简单起见，我们令 ，于是，于是(5.1.6)式变为式变为21, 1m22( )dU xEdx(5.1.7) 所以，现在解薛定谔方程是不难的。所以，现在解薛定谔方程是不难的。先来考察波函数的一般性质。先来考察波函数的一般性质。n在方势阱内在方势阱内：2212drdx (5.1.8) 其中其中 。因为。因为 。所。所以以 是实数。是实数。01VEr00 EV1r1111( )sincosxArxBrx(5.1.9) 其中其中A1和和B1

5、是待定常数。是待定常数。方程方程(5.1.8)式的解是式的解是n在方势阱外在方势阱外：222drdx(5.1.10) 其中其中 。Er222( )r xr xjjxC eD e(5.1.11) 它的一般解为它的一般解为 n其中其中j = 0，表示势阱左边的波函数，表示势阱左边的波函数，j = 1，表示势阱右边的波函数。表示势阱右边的波函数。 2200( )r xr xxC eD e2211( )r xr xxC eD e可以利用边界条件来确定系数可以利用边界条件来确定系数A1，B1，C 0，C 1，D0 和和D 1。左边有：左边有： 右边有：右边有： n因为在因为在 时，波函数要有物理意时，波

6、函数要有物理意义必须有义必须有 。所以，可以得到：。所以，可以得到： 和和 。x0)(x10C 00D 10C 00D n注意注意：对于任意的：对于任意的E值，值， 和和 的要求的要求 不可能同时得到满足。不可能同时得到满足。n下面下面求势阱内的波函数求势阱内的波函数。 在在 处（即左边界），波函数是连续的，处（即左边界），波函数是连续的，根据根据(5.1.9)式和式和(5.1.11)式，有式，有1)0(1)0(10BC(5.1.12) 同时，在同时，在 处，波函数的一次导数是连处，波函数的一次导数是连续的，对续的，对(5.1.9)式和式和(5.1.11)式求导数后，式求导数后，有有11202

7、)0()0(ArrCr0X0 xn可解得可解得211rAr(5.1.13) 同理，利用在势阱的另一边同理，利用在势阱的另一边 处，波函数处，波函数和它的一次导数都连续，得到和它的一次导数都连续，得到2111()sincosrWWrWrr2211()r Wr WWC eDe(5.1.14-1) (5.1.14-2) WX n令令(5.1.14-1)=(5.1.14-2)，建立，建立 和和 满满足的一个方程式。足的一个方程式。2111()cossinWrWrrWr222121()r Wr WWr C er De(5.1.15-1) (5.1.15-2) 令令(5.1.15-1)=(5.1.15-2

8、)，建立，建立 和和 满足的另一个方程式。满足的另一个方程式。1C1D1C1Dn由此解得由此解得 和和 为为1C1D2121)()(21)()(2122rWWeDrWWeCWrWr(5.1.16) n由上面过程可见，由上面过程可见，C 0和和D 0的值完全确定了的值完全确定了A1，B1，C1和和D1的值，结果完全确定了波的值，结果完全确定了波函数。一般情况下函数。一般情况下C1的值不一定为零。为了的值不一定为零。为了满足满足 时，时， 为零，必须要为零，必须要求求C1=0。这就对。这就对E的取值进行了限制，不能的取值进行了限制，不能取任意值，只能取某些确定值，才能保证取任意值，只能取某些确定值

9、，才能保证C1=0的要求。的要求。x)(xn我们用图解来说明上面的情况。我们用图解来说明上面的情况。 设设E1是能量本征值。对能量是能量本征值。对能量E可能有三种情况：可能有三种情况：E E1。我们画出波函数，。我们画出波函数，如图如图5.1.2所表示。由图可见，这所表示。由图可见，这3个波函数都个波函数都满足当满足当 时，时， 。但是，当。但是，当 时，只有对于时，只有对于 的波函的波函 数，数， ，即即C1 = 0。这个波函数称为方程。这个波函数称为方程(5.1.7)式的式的本征波函数。本征波函数。x ( )0 xx1EE0)(x能量本征值的确定：能量本征值的确定：图图5.2图图5.1.2

10、 不同能量不同能量E对应的波函数对应的波函数n对对E E1的情况。由的情况。由C1的表达式的表达式(5.1.14)和和(5.1.15)式可知，式可知，C10，因此在，因此在 时，时，波函数向下发散。对波函数向下发散。对E 0，当当 时，波函数向上发散。时，波函数向上发散。xxn 对第二个本征值对第二个本征值E2，同样可存在，同样可存在E值大于、值大于、 等于和小于等于和小于E2的三种情况，如图的三种情况，如图5.1.3所示。所示。图图5.1.3 本征值为本征值为E2时不同能量时不同能量E对应的波函数对应的波函数n由上述讨论可知，粒子的能量只能取得某些由上述讨论可知，粒子的能量只能取得某些 分裂

11、的值，如图分裂的值，如图5.1.4所示。所示。E的值与势阱的参数的值与势阱的参数 V0和和W有关，我们的中心问题是求解薛定谔方程有关，我们的中心问题是求解薛定谔方程(5.1.7)式的本征波函数和能量本征值。式的本征波函数和能量本征值。 V0A=0X=WE1E2E3En图图5-4 一维方势阱内粒子的能级一维方势阱内粒子的能级n我们进一步分析图我们进一步分析图5.1.2和图和图5.1.3中的波函中的波函数。波函数通过数。波函数通过 x 轴的交点称为轴的交点称为结点结点。在。在图图5.1.2中能量中能量E小于小于E1的的波函数没有结的的波函数没有结点，能量点，能量E大于大于E1的波函数有一个结点。在

12、的波函数有一个结点。在图图5.1.3中能量中能量E小于小于E2的波函数有一个结的波函数有一个结点，能量点，能量E大于大于E2的波函数有两个结点。的波函数有两个结点。n由此推广到一般规律：若由此推广到一般规律：若 和和 是薛定谔方程是薛定谔方程的两个解（不一定是本征波函数），而相应的能的两个解（不一定是本征波函数），而相应的能量是量是 和和 ， 对应着对应着n-1个结点，个结点， 对应对应着着n个结点，则第个结点，则第n个能量本征值必然处在个能量本征值必然处在 和和 之间，即之间，即 。更一般说，若。更一般说，若 和和 分别有分别有n-1个和个和n+k个结点，则能量本征个结点，则能量本征值值 必

13、然处在必然处在 和和 之间。之间。 abaEbEabaEbE2/ )(banEEEabknnnEEE,1aEbEn根据上面分析，我们得到结论：根据上面分析，我们得到结论：E1的的本征波函数没有结点本征波函数没有结点E2的本征波函数的本征波函数有一个结点，有一个结点，E3的本征波函数有二个的本征波函数有二个结点等等。对一维方势阱情况，我们用结点等等。对一维方势阱情况，我们用计算机通过确定结点来求解全部本征值计算机通过确定结点来求解全部本征值和本征波函数。和本征波函数。n第一步，给定参数第一步，给定参数。设。设 V=势阱的深度势阱的深度 W=势阱的宽度势阱的宽度 Emax为猜测的能量本征值的上限。

14、在本问题中，为猜测的能量本征值的上限。在本问题中，其值为零。由图其值为零。由图5.1.4可见，若可见，若Emax大于零，势大于零，势阱不起作用，变为自由粒子的运动。阱不起作用，变为自由粒子的运动。 Emin为猜测的能量本征值的下限。在现在的问题为猜测的能量本征值的下限。在现在的问题中，其值为中，其值为 。因为由图。因为由图5.1.4可见，粒子在势可见，粒子在势阱内运动，最小的可能能量为阱内运动，最小的可能能量为 。0V0V计算步骤：计算步骤：n第二步，将能量下限和能量上限之间分成第二步，将能量下限和能量上限之间分成M个能个能量。量。即从能量即从能量Emin开始，其增量为开始，其增量为maxmi

15、n1EEDEEM 对每一个能量，用对每一个能量，用113.141593rWGWN求出方势阱内波函数和求出方势阱内波函数和x 轴相交的结点数。轴相交的结点数。 n因为因为 是半波长的数目。如果它小于是半波长的数目。如果它小于1，取值，取值为零，没有结点。若大于为零，没有结点。若大于1而小于而小于2，取值为，取值为1，有，有一个结点。若大于一个结点。若大于2而小于而小于3，取值为，取值为2，有两个结，有两个结点 。 如 此 下 去 ， 即 可 求 得 全 部 结 点 情 况 。 在点 。 如 此 下 去 ， 即 可 求 得 全 部 结 点 情 况 。 在FORTRAN语言中，实数可以自动转化为整数

16、。在语言中，实数可以自动转化为整数。在方势阱右边，波函数是否与方势阱右边，波函数是否与x轴相交，引起对结点轴相交，引起对结点数贡献，这需要计算波函数，由波函数的系数数贡献，这需要计算波函数，由波函数的系数C1和和D1决定。若决定。若 ，则结点数增加，则结点数增加1，若，若 ，则对结点数没有影响。则对结点数没有影响。/1Wr0/11DC0/11DCn第三步，由结点数计算结果，定出能级第三步，由结点数计算结果，定出能级E1，E2，。 根据前面分析，相邻的两个能量根据前面分析，相邻的两个能量 和和 分别对分别对应应0个和个和1个结点，则能量本征值个结点，则能量本征值E1必处在必处在 和和 之间，取之

17、间，取 ，若相邻的两个能，若相邻的两个能量分别对应一个和两个结点，则能量本征值量分别对应一个和两个结点，则能量本征值E2必必处在这两个能量之间。如此等等，就可以求得一处在这两个能量之间。如此等等，就可以求得一系列的能量本征值。由能量本征值，根据前面的系列的能量本征值。由能量本征值，根据前面的公式，求得波函数的系数公式，求得波函数的系数A，B，C，D，由此确，由此确定了相应的本征波函数。用计算机可画出本征波定了相应的本征波函数。用计算机可画出本征波函数的图形。函数的图形。aEbEaEbE2/ )(1baEEE图图5.5 一维方势阱的一维方势阱的薛定谔方程的本征值薛定谔方程的本征值和本征波函数计和

18、本征波函数计算流程图算流程图5.2 粒子在辏力场中的运动粒子在辏力场中的运动n用计算机求解方势阱的问题是求解一维的薛定用计算机求解方势阱的问题是求解一维的薛定谔方程定态问题。谔方程定态问题。 ( )( )U rU rn 讨论粒子在辏力场中的运动是求解三维薛定谔方程讨论粒子在辏力场中的运动是求解三维薛定谔方程 的定态解。辏力场的势能为的定态解。辏力场的势能为表示势能和方向无关，只是粒子到力心距离表示势能和方向无关，只是粒子到力心距离 r 的函数。的函数。这种情况在实际问题中也有常遇到的情况。这种情况在实际问题中也有常遇到的情况。 n例如谐振子的势能是例如谐振子的势能是2221)(rmrU(5.2

19、.1) 又如带电粒子在一个固定点电荷所产生的电又如带电粒子在一个固定点电荷所产生的电场中运动，这是库仑场的情况，设运动粒子场中运动，这是库仑场的情况，设运动粒子的电荷为的电荷为 ，固定点的电荷为，固定点的电荷为 ，则，则粒子的势能为粒子的势能为eZerZerU2)(5.2.2) 其中其中 m 是振子的质量，是振子的质量， 是振子的频率。是振子的频率。n如果考虑原子的外层电子运动，这时内层电子的如果考虑原子的外层电子运动，这时内层电子的作用可近似考虑成电子云，它的密度为作用可近似考虑成电子云，它的密度为 ，这，这时的势能由两部分组成，即时的势能由两部分组成，即)(rrrrZerU)()(2(5.

20、2.3)第一项是核子的贡献，第一项是核子的贡献， 为带正电的核子数，第为带正电的核子数，第二项为电子云的贡献。以上三种情况都属于二项为电子云的贡献。以上三种情况都属于辏力辏力场场的情况。的情况。Zn对辏力场的情况，定态的薛定谔方程为对辏力场的情况，定态的薛定谔方程为22( )2U rEm (5.2.4) 由于由于 是球对称的是球对称的, 与方向无关与方向无关.)(rU),()(Yrru(5.2.5) 其中其中 是是 的函数，的函数， 仅仅是角度仅仅是角度 、 的函数。的函数。 rru/ )(r),(Y可用分离变数法进行求解，波函数可用分离变数法进行求解，波函数 可分解为可分解为n此时拉普拉斯算

21、符此时拉普拉斯算符 为为222222222sin1sinsin11rrrrrr将拉普法斯算符和（将拉普法斯算符和（5.2.5）式代入（）式代入（5.2.4）式，）式，可得可得22211sinsinsinYYY 0)()(2)(22 rurErUmru和和（5.2.6） （5.2.7） n方程（方程（5.2.6）式可按）式可按 和和 再分离为两个方程，再分离为两个方程，求得解析解，求得解析解， 是球谐函数。同时得到是球谐函数。同时得到 的的 值为值为),(Y, 2, 1, 0) 1(lll方程（方程（5.2.7）式是）式是径向波动径向波动方程，在一般的方程，在一般的 情况时，通常很难求得（情况时

22、，通常很难求得（5.2.7）式的解析解。这）式的解析解。这是一个二阶常微分方程，可以用其它数值解法来是一个二阶常微分方程，可以用其它数值解法来求解求解.（留作练习题，可以用氢原子为例）。（留作练习题，可以用氢原子为例）。)(rUn现在用计算机对（现在用计算机对（5.2.7）式进行数值解。为）式进行数值解。为书写简单，设书写简单，设22) 1()(2)(rllErUmrf于是（于是（5.2.7）式变为）式变为)()()(rurfru (5.2.8) n为了近似进行数值计算，考虑为了近似进行数值计算，考虑 函数在函数在 处作泰勒展开处作泰勒展开 )(! 4)(! 3)(! 2)()()(432ru

23、hruhruhruhruhru )(! 4)(! 3)(! 2)()()(432ruhruhruhruhruhru(5.2.9)(5.2.10) )(hrur将（将（5.2.9）式和（）式和（5.2.10）式相加除）式相加除2，得以，得以)(! 6)(! 4)(! 2)()()(21) 6(6) 4(42ruhruhruhruhruhru (5.2.11) n对（对（5.2.11）式两边进行二次求导数，得到）式两边进行二次求导数，得到 )(! 4)(! 2)()()(21)6(4)4(2ruhruhruhruhru(5.2.12) 将方程（将方程（5.2.12）式乘）式乘 与（与（5.2.11

24、）式相减，）式相减，消去消去 的项，结果为的项，结果为 12/2h4h12)()(21)()(212hhruhruhruhru 121301)(! 4)(2)(12)()6(622ruhruhruhru(5.2.13) n忽略大于等于忽略大于等于 的项。同时利用薛定谔方程的项。同时利用薛定谔方程(5.2.8)6h)()()()()()()()()(hruhrfhruhruhrfhrururfru 代入（代入（5.2.13）式，得到）式，得到)(121)(21)(121)(2122hrfhhruhrfhhru)()(126)(121)(22rurfhrfhru(5.2.14) n令令 ，于是（，

25、于是（5.2.14）式变为）式变为)(12)(2rfhrT)(102)()(1)()(1)(rTruhrThruhrThru(5.2.15) （5.2.15）式是我们要求的结果。通过上述运）式是我们要求的结果。通过上述运算，我们将二阶的微分方程变成了（算，我们将二阶的微分方程变成了（5.2.15）式的主方程。这个方法称为式的主方程。这个方法称为Numorovs方法。方法。在主方程中舍去了在主方程中舍去了 的项，因此这个方法在理的项，因此这个方法在理论上达到五级精确度，误差数量级为论上达到五级精确度，误差数量级为 。 6h)(6hOn由舍去的部分由舍去的部分)(240)6(6ruh误差(5.2.

26、16) 主方程（主方程（5.2.15）说明径向波函数）说明径向波函数 、 和和 的三点数值之间的关系，由前面二点的三点数值之间的关系，由前面二点的值可求得第三点的值。但是其中的的值可求得第三点的值。但是其中的 函数没有函数没有确定，确定， 为为 )(hru)(ru)(hruff22) 1()(2)(rllErUmrf(5.2.17) n其中其中 可由具体物理问题的势能模型确定下来，可由具体物理问题的势能模型确定下来，而而 是未知的。解主方程（是未知的。解主方程（5.2.15）式必须同）式必须同时求本征值和本征波函数。时求本征值和本征波函数。n在具体求解以前，我们首先分析一下波函数可在具体求解以

27、前，我们首先分析一下波函数可能具有的形式。能具有的形式。n为此，考察为此，考察 函数的性质，由（函数的性质，由（5.2.17）式）式 的表示及的表示及 的形式，对库仑场情况（的形式，对库仑场情况（5.2.2式式和和5.2.3式），式），将将 分成三个区域分成三个区域：)(rf)(rUEf)(rUrn第一个区域，第一个区域， 很小，这时很小，这时 为为rf2)(rArf第二个区域，第二个区域， 取一般数值时取一般数值时 rErrrBrf)()(第三个区域，第三个区域，1rErf)(n定性图形如图定性图形如图5.2.1所示。所示。r1f(r)图图 5-6 原子的原子的f(r) 定性图定性图 n对应三个区域对应三个区域 第一区域第一区域， ，设波函数，设波函数0fu于是于是 areu araeu fuuaeauar 22因此，因此， ，可得，可得2/1fa rfeu2/1(5.2.18) n第二区域第二区域， ，设，设0fikreu 即取即取 ，于是有，于是有krusinikrikeu fuukekuikr 22因此，因此， ，可得，可得2/1fk )sin(2/1rfu (5.2.19) n第三区域第三区域， ，设，设0fbreu于是于是 brbeufuebubr 2因此，因此， ，可得，可得2/1fb rfeu2/1(5.2.20) n波函数的定性