3、解一个简单的一阶常微分方程。
![]()
先来分析一下这个方程，ψ是要求的函数，
是一阶导数，这个表达式的含义就是：所求的函数上每个点的斜率等于这个点在ψ轴上的值乘以(-k)。
我们给定一个初始点，再给定一个x轴的步长h，下一个点的ψ值等于初始点ψ值加上下一个点的一阶导数乘以步长h，这样一步步延伸，延伸的线段就是所求函数的近似值。步长越短，解越精确。
计算第(i)个点的一阶导数。
![]()
计算第(i)个点的函数值ψ。
函数初始值设为0.5，自变量x初始值设为 -0.8，常数K初始值设为1，线段步长设为0.01，程序得到如下图像：
(3)_解一个简单的一阶常微分方程.zip
(7.72 KB, 下载次数 Times of downloads: 3)
用鼠标滚轮可以对视图进行缩放操作；按键盘
A键曲线加粗；按
Z键曲线变细；按
↑键起始点上移；按
↓键起始点下移；按
←键起始点左移；按
→键起始点右移。越精确。
所有得到的这些曲线都是方程的解。实验结果是动态的，非常直观，这是无法通过任何理论公式体验到的。
4、解一个简单的二阶常微分方程。
这是一维薛定谔方程：
![]()
我们把一维定态薛定谔方程的二阶导数放在左边，把右边
这一坨看作常数k，就得到如下简单的二阶常微分方程。
先来看一下这个方程的解法，把这个方程的解法搞清楚后，解完整的一维方程就简单多了。
是二阶导数，那么一阶导数是多少呢？我们知道，二阶导数是一阶导数的导数，一阶导数也是所求的函数，也就是说我们在求解二个相互关联的函数：原函数和一阶导数。好了，我们依旧用上一节的方法来解一下这个方程。
计算第(i)个点的二阶导数。
计算第(i)个点的一阶导数。
计算第(i)个点的函数值ψ。
函数初始值设为0，一阶导数初始值设为1，常数K初始值设为1，线段步长设为0.1，程序得到如下图像：
(4)_解一个简单的二阶常微分方程.zip
(7.75 KB, 下载次数 Times of downloads: 2)
可以看出，借助计算机做重复复杂的计算工作，可以让我们把精力集中到现象的观察中，这样能大大提高新方法的发明及理解。
原子单位制下的所有物理量全部用同一单位a.u.，使用原子单位后，薛定谔方程大大简化，下面是一些常用原子单位。
长度1a.u. = 玻尔半径 = 0.5291772083D-10m
时间1a.u. = 2.418884326505(16)D-17s
速度1a.u. = 2187691.2633(72)m/s
质量1a.u. = 电子质量 = 9.10938188D-31Kg
质子质量 = 1836.1527a.u.
原子量单位(amu) = 碳十二质量的十二分之一 = 1822.88848a.u.
动量1a.u. = 1.99285166(34)D-24Kg*m/s
能量1a.u. = 4.35974417(75)D-18J
角动量1a.u. = 1.05457160D-34J*s
电荷1a.u. = 1.602176462D-19C
力1a.u. = 8.2387225(14)D-8N
真空电容率 = 0.07957747154594767(e^2a0me/h^2)
约化真空电容率1a.u. =
电势能 =
=
约化普朗克常数1a.u. =
=
=
=
这是一维自由粒子薛定谔方程：
假设我们所求的粒子是电子，质量m就是1a.u.，约化普朗克常数也是1a.u.，上面的方程就化为下面的简单方程。
可以看出来这个方程跟上一节的方程差别不大，但是为了提高计算精度，本节采用3步计算法。
第1步-预测
预测第(i)个点的二阶导数。
预测第(i)个点的一阶导数。
预测第(i)个点的函数值ψ。
第2步-校正
校正第(i)个点的二阶导数。
校正第(i)个点的一阶导数。
校正第(i)个点的函数值ψ。
第3步-计算
计算第(i)个点的二阶导数。
计算第(i)个点的一阶导数。
计算第(i)个点的函数值ψ。
计算第(i)个点的函数值ρ。
函数初始值设为0，一阶导数初始值设为1，能量E初始值设为1，线段步长设为0.1，程序得到如下图像：
(5)_原子单位制，一维自由粒子.zip
(7.94 KB, 下载次数 Times of downloads: 1)
用鼠标滚轮可以对视图进行缩放操作；按键盘A键曲线加粗；按Z键曲线变细；按↑键步长线段加长；按↓键步长线段减短；按←键整条曲线左右加长；按→键整条曲线左右减短；按˂键能量E值增大；按˃键能量E值减小；按PageUp键函数曲线纵向拉长；按PageDown键函数曲线纵向缩短；按
Home键波函数ψ与概率密度函数ρ相切换。
绿色的是一阶导数，黑色的是函数。能量E的值越大，函数周期越短，反之越长。
我们可以看出，数值解法不需要复杂的数学计算，不需要很深的数学功底，通过简单的程序设计就能解决问题。
一维谐振子是量子力学的简谐运动，粒子受中心力场的引力，粒子离中心力场越远引力越大，引力大小跟距离的二次方成正比。
这是一维谐振子薛定谔方程：
假设我们所求的粒子是电子，应用原子单位制，质量m就是1a.u.，约化普朗克常数file:///C:\Temp\ksohtml\wpsED4C.tmp.png也是1a.u.，中心力场设为原点，常数k设为1，上面的方程就化为下面的简单方程。
本节依然采用3步计算法。
第1步-预测
预测第(i)个点的二阶导数。
预测第(i)个点的一阶导数。
预测第(i)个点的函数值ψ。
第2步-校正
校正第(i)个点的二阶导数。
校正第(i)个点的一阶导数。
校正第(i)个点的函数值ψ。
第3步-计算
计算第(i)个点的二阶导数。
![]()
计算第(i)个点的一阶导数。
计算第(i)个点的函数值ψ。
计算第(i)个点的函数值ρ。
用鼠标滚轮可以对视图进行缩放操作；按键盘A键曲线加粗；按Z键曲线变细；按↑键步长线段加长；按↓键步长线段减短；按←键整条曲线左右加长；按→键整条曲线左右减短；按˂键能量E值增大；按˃键能量E值减小；按PageUp键函数曲线纵向拉长；按PageDown键函数曲线纵向缩短；按Home键波函数ψ与概率密度函数ρ相切换。
函数初始值设为0，一阶导数初始值设为0.1，能量E初始值设为0.5，线段步长设为0.05，曲线纵向缩放系数设为 -37，从理论上讲，基态能量E为0.5a.u.，但数值解是有误差的，无论多精确，函数右边都会拉起来，按˂键或˃键，微调能量大小，使曲线右边拉到最长，这时的能量值就是数值计算所得到的能量值，跟理论值0.5是有误差的。基态概率密度函数ρ的图像如下图所示：
按Home键，切换到波函数ψ图像，这时函数图像成为一条直线，连续按PageUp键，曲线纵向拉长，得到如下基态波函数ψ的图像：
按Home键，切换回概率密度函数ρ图像，连续按PageDown键，曲线纵向缩回原始状态，连续按˂键，能量不断调大，函数曲线右侧逐渐拉起来，当能量值接近1.5时函数曲线右侧又逐渐落下来，曲线由单峰变为双峰，或者在左上角输入框内填1.5，按能量按钮，直接得到第一激发态概率密度函数ρ的图像，按˂键或˃键，微调能量大小，使曲线右边拉到最长，这时的能量值就是数值计算所得到的能量值，跟理论值1.5也是有误差的。第一激发态概率密度函数ρ的图像如下图所示：
重复上述步骤，得到下面的图像。
(6)_一维谐振子.zip
(8.66 KB, 下载次数 Times of downloads: 3)
一维谐振子的计算结果是比较完美的，跟理论计算结果吻合得比较好。
把氢原子用圆管束缚起来，当圆管截面半径趋于零时，就得到了一维氢原子。
一维氢原子和一维谐振子运动状态相似，粒子都受中心力场的引力，所不同的是氢原子电子离原子核越远引力越小，引力大小跟距离的倒数成正比。
如果把原子核看做质点，就得到下面的一维氢原子薛定谔方程：
应用原子单位制，电子质量m就是1a.u.，电子电荷e是 - 1a.u.，质子电荷e是1a.u.，约化普朗克常数file:///C:\Temp\ksohtml\wpsF115.tmp.png也是1a.u.，中心力场设为原点，约化真空电容率file:///C:\Temp\ksohtml\wpsF116.tmp.png是1a.u.，上面的方程就化为下面的简单方程。
本节依然采用3步计算法。
第1步-预测
预测第(i)个点的二阶导数。
预测第(i)个点的一阶导数。
预测第(i)个点的函数值ψ。
第2步-校正
校正第(i)个点的二阶导数。
校正第(i)个点的一阶导数。
校正第(i)个点的函数值ψ。
第3步-计算
计算第(i)个点的二阶导数。
![]()
计算第(i)个点的一阶导数。
计算第(i)个点的函数值ψ。
计算第(i)个点的函数值ρ。
用鼠标滚轮可以对视图进行缩放操作；按键盘A键曲线加粗；按Z键曲线变细；按↑键步长线段加长；按↓键步长线段减短；按←键整条曲线左右加长；按→键整条曲线左右减短；按˂键能量E值增大；按˃键能量E值减小；按PageUp键函数曲线纵向拉长；按PageDown键函数曲线纵向缩短；按Home键波函数ψ与概率密度函数ρ相切换。
这是理论计算的能级公式（原子单位制）。
一维氢原子的求解遇到了很大的困难，我得到了很多理论计算结果中没有的解，量子数n是非整数，但是这些解又是符合方程条件的，我很困惑，先来展示没有问题的解吧。
函数初始值设为0，一阶导数初始值设为0.1，能量E初始值设为 -0.5，线段步长设为0.01，曲线纵向缩放系数设为 -120，从理论上讲，基态能量E为 -0.5a.u.，但数值解是有误差的，无论多精确，函数右边都会拉起来，按˂键或˃键，微调能量大小，使曲线右边拉到最长，这时的能量值就是数值计算所得到的能量值，跟理论值 -0.5是有误差的。基态概率密度函数ρ的图像如下图所示：
按Home键，切换到波函数ψ图像，这时函数图像成为一条直线，连续按PageUp键，曲线纵向拉长，得到如下基态波函数ψ的图像：
按Home键，切换回概率密度函数ρ图像，连续按PageDown键，曲线纵向缩回原始状态，在左上角输入框内填-0.125，按能量按钮，曲线变为第一激发态概率密度函数ρ的图像，按˂键，微调能量大小，使曲线右边拉到最长，如下图所示：
重复上述步骤，得到下面的图像。
(7)_一维氢原子.zip
(8.69 KB, 下载次数 Times of downloads: 4)
这些解是跟理论计算结果一致的，再来看看那些不正常的解。
(7)_一维氢原子_第1个不正常解.zip
(8.69 KB, 下载次数 Times of downloads: 2)
(7)_一维氢原子_第2个不正常解.zip
(8.7 KB, 下载次数 Times of downloads: 1)
这些解的量子数n都不是整数，却符合方程给定的条件，那么这些解的物理意义是什么呢？我认为数学解法属于科学理论，需要严格论证，并且要跟实验结果吻合；数值解法则属于工程技术，需要积累大量的经验，并且要跟理论结果相一致。既然这些解是不正常的，如何在数值计算中避免，还望各位大神赐教。
发表于 Post on 2018-2-26 00:09:36
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author