浙江大学【物理学】毕业论文基于FlexPDE 软件实现薛定谔方程求解的可视化.doc

浙江大学【物理学】毕业论文：基于FlexPDE 软件实现薛定谔方程求解的可视化 本科生毕业论文 ( 2012 )届论文题目:基于 FlexPDE软件实现薛定谔方程求解的可视化 学 院: 物理科学与技术学院专 业:物理学 学 号: 200810800094姓 名:韦 贝 佩指导老师姓名及职称: 王 宁教授2012年 5月 目 录内容摘要.- 2 - 1 引言- 3 - 2 FLEXPDE软件.- 6 - 2.1 FLEXPDE软件简介.- 6 - 2.2 FLEXPDE软件初始界面.- 6 - 2.3 新建FLEXPDE脚本.- 8 - 2.3.1新建FlexPDE脚本操作.- 8 - 2.3.2新建FlexPDE脚本界面.- 8 - 2.3.3新建FlexPDE脚本程序介绍.- 9 - 2.3.4程序编写.- 11 - 2.4 打开和运行FLEXPDE脚本- 12 - 3 运用FLEXPDE解薛定谔方程- 14 - 3.1 一维线性谐振子.- 14 - 3.2 三维各项同性谐振子- 18 - 3.4 三维WOODS-SAXON势.- 21 - 4 总结与展望.- 22 - 参考文献.- 23 - ABSTRACT.- 24 - 致 谢- 25 - 1 -基于FlexPDE软件实现薛定谔方程求解的可视化专业:物理学 学号:200810800094 学生姓名:韦贝佩 指导老师姓名:王宁 内容摘要 量子力学是现代物理学各分支学科的基础。而要学习好量子力学课程,首先要解决的关键问题就是求解薛定谔方程。对于复杂系统,运用数学方法求解薛定谔方程是比较复杂的;甚至,薛定谔方程的精确求解依然是量子力学中的一大难题。FlexPDE是一个基于有限元方法求解偏微分方程的商业软件,对求解复杂系统的薛定谔方程有着重要的应用价值。 本文主要论述了FlexPDE的一般操作以及运用FlexPDE求解薛定谔方程的具体实例介绍。首先我们介绍了FlexPDE的界面和基本程序结构。然后基于该软件,解决了在复杂系统中求解薛定谔方程的困难,实现薛定谔方程求解的可视化。通过对该软件的介绍和运用该软件进行求解薛定谔方程的具体实例操作,为处理量子力学中更为复杂的情况打下基础,从而有利于我们更好地学习量子力学和进行相关教学。关键词: FlexPDE; 偏微分方程;薛定谔方程;- 2 -1 引言 对于微观世界的研究,离不开量子力学。量子力学是物理学专业学生的一门必修课,更是现代物理学的各个分支学科及相关的边缘学科的基础。量子力学中最核心的问题就是解决波函数如何随时间演化以及在各种具体情况下找出描述体系状态的各种可能的波函数【1】。薛定谔在1926年提出的波动方程(薛定谔方程)成功的解决了这个问题,因而我们学习量子力学要解决的一个关键问题就是求解薛定谔方程。 薛定谔方程是量子力学的一个基本假定,并不能从什么比它更根本的假定来证明它,它的正确性只能靠实验来检验。薛定谔方程是将物质波的概念和波动方程相结合建立的二阶偏微分方程,揭示了原子世界中物质运动的基本规律,每个微观系统都有一个相应的薛定谔方程式,通过解方程可得到波函数的具体形式以及对应的能量,从而了解微观系统的性质。 定态薛定谔方程的表达式为: 2h2- j +V xj Ej2m这是一个二阶线性偏微分方程,x,y,z是待求函数,它是 x,y,z三个变量的复数函数(就是说函数值不一定是实数,也可能是复数)。式子最左边的倒三角是一个二阶偏微分算符,意思是分别对x,y,z的x,y,z坐标求偏导的平方和。式中的 V(r)是代表势场,视实际问题而定的。例如在一维无限深势井中,势场为: 0,0xa Vx ?,其他 要解薛定谔方程,仅仅给出势场的值是不够的,还需给出一定的边界条件。而实际问题中的边界条件是由波函数的性质(连续性、单值性、有限性)【2】得出。例如在一维无限深势井中,在x0和xa处,我们可以认为波函数为零即粒子出现在势井外的概率几乎为零。若要用数学的方法解薛定谔方程,则不同类型的问题要用到不同的方法去解决。对于在一维无限深势井中求解薛定谔方程还是比较简便的,但是对于在其他势井中求解薛定谔方程,甚至增加到多维时,则求解过- 3 -程将会比较复杂。例如在一维线性谐振子中求解此薛定谔方程时,需要用到Hermite方程,甚至 Bessel方程等特殊方程;不仅如此,相应的图像我们根本不会画。另外,对于复杂系统,薛定谔方程的精确求解依然是量子力学中的一大难题。 从而,对于初学者来说,学习量子力学时,将会感到非常困难和抽象。因为:用数学方法求解薛定谔方程过程的复杂性;不能结合图像来了解其中原理或进行实际分析;不能预见性的知道求解结果如何。对于老师来说,在教授量子力学课程时,也只能按照教材的原理和加上个人的理解进行口头教学,而不能结合具体图像(特别是三维图像)来进行教学。 同样,在其他学科领域中,我们也经常会遇到偏微分方程的求解和画图的问题。这往往是一项复杂的工程,也是我们学习和教学中的一个难点。所以我们需要寻找一个一般的高效的解决方案。 FlexPDE软件是使用有限元方法给出线性和非线性偏微分方程数值解的一个工具软件。它可以有效地解决上述问题(能很好的帮助我们解决这一问题)。运用 FlexPDE软件,求解偏微分方程可以自动运行,可以得到想要的各种图形与结果,有利于我们进行学习和教学。 另外,运用 FlexPDE解偏微分方程还具备以下优点: 1.求解速度快。当输入了完整的解决方案后,计算机可自动运行,整个过程可在短时间内完成。 2.程序编写简单。自动生成程序的主要部分,按照生成程序提示修改完整即可。程序语言也简单易懂。 3.程序的可移植性强。对于同类型问题,只要把程序稍作修改,问题即可迎刃而解,另外软件本身还提供了一些常用的运算函数。 4.求解结果明了易懂。计算机自动运行结果会以文字和图像的方式显示,让人容易获得其信息,且能帮助进行实际分析,实现可视化。 5.求解精确,精度可控性强。FlexPDE可提供误差分析,通过 SELECT控件可以控制解决方案的误差范围。 可是,目前,我们可以容易找到英文版的 FlexPDE说明书,却无中文版,而且关于 FlexPDE的文章也比较少。另外,在我找到的为数不多的几篇相关领- 4 -域研究人员发表的基于 FlexPDE 的相关研究的论文或期刊中,虽涉及到FlexPDE,但文章中都只稍微介绍了一下 FlexPDE,更多的是阐述具体应用方面的知识,而且大多数并没有说明怎样运用 FlexPDE软件来实现相关研究的详细过程。基于以上现状,我们可以知道,目前会应用 FlexPDE软件的人员不多,仅限于高知识水平的研究人员在使用。而对于一般知识水平,特别是英文不太好的人员来说,FlexPDE软件是比较难学的,所以不能普遍的应用到我们的学习和教学中。 基于以上困难,本文首先通过参阅了英文版的教程和相关文章,旨在使用中文的形式,说明 FlexPDE软件的使用方法。然后重点阐明运用 FlexPDE软件去求解量子力学中的薛定谔方程的具体过程及结合图像进行结果分析,从而使得省去了解微分方程的繁杂过程。目的是帮助大家学习这一软件,并很好的应用于学习和教学中,以及对薛定谔方程有进一步的深刻了解。- 5 -2 FlexPDE软件 2.1 FlexPDE软件简介 FlexPDE软件的全称为Finite Element Solution Environment for Partial Differential Equations即偏微分方程有限元软件。FlexPDE是 PDE Solution公司专门开发的使用有限元方法给出线性和非线性偏微分方程数值解的一个工具软件。该软件可以在一维到三维区域内求解微分方程和求解含时方程,并将其结果以图片形式输出。FlexPDE软件的工作原理就是把那些可以用偏微分方程描述的系统转化为有限元模型并求解 。其中,网格的数量和密度可根据误差自动调整 。 软件推出后 ,就被世界上许多大学和实验室采用,其中最常用于物理、电机、电子、通讯、土木、机械、化工、化学、生物学、地质学、数学等需运用到偏微分方程的科学领域。因为多数的物理和工程问题在某种程度上可以被描述成为部分微分方程。许多领域的学者可以使用FlexPDE来将他们的试验或设备建模,进行重要的预测或检验;工程师也可以使用FlexPDE来设计最优化研究,适合度研究以及概念分析等。 2.2 FlexPDE软件初始界面 当我们运行FlexPDE软件时我们将看到如图 1界面图 1- 6 -上图的中文注释如图2图 2当用鼠标点击各项菜单时,会出现下拉菜单。界面和注释如下: 一、文件菜单如图3图 3 二、控制菜单如图4图 4三、停止菜单如图5(当一个问题正在运行时,有时候需要用到这个菜单栏) 图4图 5四、编辑菜单如图6图 6 - 7 -五、帮助菜单如图7图 72.3 新建 FlexPDE脚本 2.3.1新建 FlexPDE脚本操作 当我们要新建一个脚本时,可以点击 File(文件)?New Script(新建一个脚本),则会出现如图8界面图 8然后在文件名栏输入文件名,保存类型选择.PDE。然后点击保存。则新建的脚本将保存在FlexPDE目录下,生成的文件是以“.pde ”为后缀扩展名的文件。2.3.2新建 FlexPDE脚本界面 新建FlexPDE脚本的主界面如图 9所示。- 8 - 图 9 2.3.3新建 FlexPDE脚本程序介绍 FlexPDE程序称为问题描述符。它是通过一些分段的关联语句部分使问题系统化。每个部分都是由一个恰当的名字开头,后面跟一个或多个用来定义问题的语句。问题描述符如上图右边窗口所示。红色字体部分为描述符的主要结构,蓝色字体和黑色字体部分为我们可以需要根据具体问题改写部分。而绿色字体部分为注释,编写时可留可不留。 - 9 -以上是单个问题的描述符,而当遇到同时要用多个问题时。问题描述符的格式为: section 1statementsection 2statement 1statement 2*section 3statement 1statement 2 下面我们再详细的介绍一下问题描述符的几个重要部分 1. SELECT 在SELECT部分,它是可选的,用于必要时,覆盖默认的一些内部的选择方案。SELECT部分是用来控制流程来解决问题的。下面将阐述两个从属SELECT节的解决方案控制。ERRLIM这是主要的精确控制。无论是空间误差控制XERRLIM的时间误差控制TERRLIM。 注:ERRLIM是一个相对误差在依赖估计变量。该解决方案不能保证尽在此错误。这可能有必要调整或手动强制更大ERRLIM网密度达到预期的解决方案的准确性。 MODES 选择特征值求解,并指定所需的数模式。默认情况下不运行一个特征值问题。 2. Variables 变量部分是用来定义变量名的,描述符的各部分变量都要按照定义的变量名书写。在指派变量名时要遵守以下规则: 1、 变量名称必须以字母开头。不能以一个数字或符号开头。 2、 变量名称可以用字符t以外的单个字符,因为这是留给时间变量。 3、 变量名称可为任何长度和任何组合的字符、数字及/或符号,但除了描述符中本来自带的词。 4、 变量名称不得有分隔符。复合名称可用_表示 例如temperature_celsius 5、 变量名称不得使用 - ,这被认为是减法符号。 3. BOUNDARIES 边界部分是用来描述方程的求解区域。每个问题描述符都必修包含边界部- 10 - 分。问题边界是由在二维笛卡尔空间内行走的周边路径构成。以这种方式,物理区域分解为REGION ,FEATURE 和 EXCLUDE 小分部。每个问题描述符至少必须有一个REGION分部,FEATURE 和 EXCLUDE分部可以选择。下面我们主要介绍REGION分部。 REGIONS 一个区域是一个二维问题域或三维投影的问题域部分,由路径构成封闭区域但一维区域除外。例如: REGION 1 an outer box START0,0LINE TO 10,0 TO 10,10 TO 0,10 TO CLOSE这里表示的区域是如图10所示 (0,10) (10,10) (10,0) (0,0) 图 102.3.4程序编写 我们已经了解了FlexPDE程序,下面我们开始编写程序(求解正弦函数的偏微分问题)。 我们先看看这个函数fcosx的导数,按照数学知识,我们可以清楚的知道,cosx的导数是 sinx。这里我们可以用 FlexPDE尝试运算一下。具体程序编写如下: TITLE定义问题的标题 cosx COORDINATES cartesian1 这是一维应用 DEFINITIONS 定义各辅助参量 Lx10 fcosx 定义函数 fxdxf f的导数 - 11 - BOUNDARIES region domain 描述的区域 start-Lx line to Lx 描述的区域为-Lx到 Lx PLOTS elevationf, fx from -Lx to Lx 显示 f和 f的导数的图像,从-Lx到 Lx END 程序结束的标志 把程序以cos(x).pde保存到 FlexPDE5中。 2.4 打开和运行 FlexPDE脚本 而当我们要打开一个脚本时,可以点击File(文件)?Open Script(打开一个脚本),这里我们以打开刚才保存的 cos(x).pde为例。则会出现如图 11界面图 11当我们需要知道求解结果时,可以点击 (运行),计算机将开始进行运算,运算结束后,结果以图片的形式输出,如图12- 12 -图 12图像 a表示函数 f的图像,图像 b表示 f的导数,图像 b刚好是我们熟知的函数 sinx的图像。 另外,我们还可以进行三种操作:一、在状态窗口中右击鼠标选择字体。二、在图像窗口中双击图像至最大化,同样双击最大化图像还原。三、可以通过 和工具切换编辑窗口和图像窗口。 - 13 - 3 运用 FlexPDE解薛定谔方程 在学习量子力学时,我们经常需要在一定的势场中求解定态薛定谔方程。但在有些势场中运用数学方法求解定态薛定谔方程,是比较复杂和困难的。下面具体介绍几种运用 FlexPDE软件求解这一类问题。 3.1 一维线性谐振子 2 2一维线性谐振子的势场为 Vxm x /2。以10为例,参考程序如下: TITLE Schrodinger equation with 1D harmonic oscillator potential SELECT ERRLIM1E-4精度控制 modes15 讨论前十五级 COORDINATES cartesian1 一维问题 VARIABLES Phi定义波函数 Phi DEFINITIONS K-197.3*197.3/2*938定义辅助量 R0Pi*2 V0.5*938*100/197.3/197.3*x*x Vxm2*x2/2 shift 0Elambda+shift能级 EQUATIONS Divk*gradPhi +V*Phi LAMBDA*Phi + shift*Phi薛定谔方程 BOUNDARIES REGION 1问题边界的区域 1 START -15POINT valuePhi0 波函数的有限性 line to 15 POINT valuePhi0 波函数的有限性- 14 - PLOTS ELEVATION20*Phi+E,V,E FROM -12 to 12输出 20+E,V,E的正视图 ELEVATION20*Phi*Phi+E,V, E FROM -12 to 12 输出 20*+E,V,E的正视图SUMMARY REPORTE数值输出能级 E END 程序结束标志运行上述程序,可得到能级:为了方便观察,我们把波函数增大 20倍。所以输出在 E的基础上波函数增大 20倍20+E,势场V和能级 E的组合图如图 13至图 16。分别表示 n0至 n3的图形。 图 13 图 14 - 15 -图 15 图 16以及在 E的基础上波函数的平方增大 20倍20*+E,势场V和能级 E的组合图如图 17至图 20。分别表示 n0至 n3的图形。图 21 图 22图 17 图 18 图 19 图 20- 16 - n0至 n9的波函数和对应的 E以及 V放到同一个图上是如图 21;同样把 n0至 n9的波函数的平方*和对应的 E以及 V放到同一个图上是如图22。 图 21图 21 图 22 我们可以看到以上这两个图与曾谨言量子力学卷第四版第 84页的图形是基本一致的,说明我们用 FlexPDE软件解出的波函数与实际波函数相符。我们再来看看解出的能级与实际能级之间的差别。如表 1 表 1 FlexPDE解出的能级E2与实际能级E1比较表 绝对误差n E1n+0.5*hbar* E2 相对误差(|E2-E1|/E1)*100% |E2-E1| 0 5 5.000000 0.00000 0.00000% 1 15 15.00000 0.00000 0.00000% 2 25 25.00001 0.00001 0.00004% 3 35 35.00002 0.00002 0.00006% 4 45 45.00005 0.00005 0.00011% 5 55 55.00009 0.00009 0.00016% 6 65 65.00015 0.00015 0.00023% 7 75 75.00024 0.00024 0.00032% 8 85 85.00036 0.00036 0.00042% 9 95 95.00050 0.00050 0.00053% 10 105 105.0007 0.00070 0.00067% 11 115 115.0010 0.00100 0.00087% 12 125 125.0015 0.00150 0.00120% - 17 - 13 135 135.0020 0.00200 0.00148% 14 145 145.0022 0.00220 0.00152% 从上表可以看出解出的能级与实际能级相差很小,说明 FlexPDE软件可以较精确的解出薛定谔方程的波函数以及相应能级,且以图形或数据的形式输出结果,实现了薛定谔方程求解的可视化。 另外 FlexPDE软件还能帮助我们把复杂的抽象的三维薛定谔方程的求解可视化,下面我们以三维各项同性谐振子为例。 3.2 三维各项同性谐振子 2 2三维各项同性谐振子的势能表达式为 Vrm r /2。参考程序如下: TITLE 3D_harmonicCOORDINATEScartesian3 三维问题 VARIABLESPhi 定义变量 Phi SELECTMODES20 讨论的级数 ERRLIM1E-4 精度控制 DEFINITIONSR115 定义辅助量 K197.3*197.3/2*938 hbar2/2*mu V0.5*938*100/197.3/197.3* x2+y2+z2 Vxm2*r2/2 ZspheresqrtR12-x2-y2,0定义辅助量shift0 Elambda+shiftEQUATIONS Div-K*gradPhi +V*Phi LAMBDA*Phi + shift*Phi 薛定谔方程Extrusionsurface Sphere BottomZ-Zsphere 下半区域 layer can - 18 - surface Sphere TopZZsphere 上半区域BOUNDARIES REGION 1surface Sphere Bottom ValuePhi0波函数连续性surface Sphere Top ValuePhi0 波函数连续性STARTR1,0 从(R1,0)开始ValuePhi0波函数连续性ARC center0,0 Angle360 TO CLOSE 以(0,0)为中心旋转 360度! TIME 0 TO 1 if time dependent MONITORSPLOTSCONTOURPHI on z0.0 等高线输出 PhiSURFACEPHI on z0.0 输出 Phi的曲面图 ELEVATIONPhi From -12,0,0 to 12,0,0输出 Phi的正视图SUMMARY export mergeREPORTE数值输出能级 EEND 程序结束标志 运行程序,得到的能级如下: 图 23- 19 - 从上图可以很容易看出,三维各向同性谐振子的能谱与一维谐振子的能谱相似,都是均匀分布;但也有不同之处,三维各向同性谐振子的能级一般是简并的。三维各向同性谐振子的能级能级表达式为 EN+3/2*hbar*,其中 N 2n + l。r因此,对于能级 EN, l和n有多种组合形式,且分析可知三维各向同性谐振子 ENr1能级的简并度为 f N +1N + 2。所以当 N0时 f 1, N1时 f 3等等。n n n2以上三维谐振子的能级数据结构正反映出了量子力学中的能级简并的特点。 而波函数切面图,三维图和曲线图如图 24至 26,为了简洁,只给出 Mode 2的图形。图 24 图 25 图 26运用 FlexPDE,我们可以轻松的解出三维谐振子的能级以及得出一到三维的图形。以至于帮助我们在学习和教学中进行形象的分析。 前面我们已经论证了运用 FlexPDE求解薛定谔方程的准确性。所以我们还可以运用 FlexPDE预见性地求解在较复杂势场中的薛定谔方程。Woods-Saxon势是更接近于实际原子核中的单粒子势,它的性质介于三维各向同性谐振子势和球方势阱之间【1】。下面我们以求解在三维的 Woods-Saxon势中的薛定谔方程为例。看看在此场中的波函数是怎样的图形,是否与在三维各向同性谐振子势中的波函数相似。- 20 - 3.3 三维 Woods-Saxon势 Woods-Saxon势表达式为:V r -Vo/ 1 + exprR/a(Vo,R,a0)参考程序只要在三维线性谐振子的基础上稍微修改即可,增加定义一个R22*Pi和修改为 V-50/1+expsqrtx2+y2+z2-R2/0.5。 运行程序得到波函数的切面图,三维图,曲线图如图 27至图 29(同样只给出 Mode 2)。 图 27 图 28 图 29以上我们可以看出 Woods-Saxon势的波函数图像与三维各向同性谐振子的图像相似。这样我们就可以把在 Woods-Saxon势下的薛定谔方程的求解可视化了。同样我们也可以运用 FlexPDE求解其他势场下的薛定谔方程,只要稍微改动程序即可。- 21 - 4 总结与展望 总结 本文在引言中,首先提出了量子力学中求解薛定谔方程的困难,而这一问题可用 FlexPDE软件解决;其实提出了 FlexPDE软件的相关中文资料比较难找的问题。 基于 FlexPDE软件的相关中文资料比较难找的问题,本文在第二章中,介绍了 FlexPDE软件的基本功能和简单操作,并结合简单的微分当成进行讲解,帮助初学者轻松的学习这一软件。第三章是本文的重点,详细介绍了应用FlexPDE求解薛定谔方程以及结合求解结果进行教学分析。 展望 本论文研究进展中,发现一些问题可作进一步研究,比如以下几方面: 1. FlexPDE的英文版说明书比较好找,而无中文版。展望以后能有 FlexPDE的中文版说明书,或者汉化的 FlexPDE软件。 2. FlexPDE的相关资料比较少,展望能有更多的研究者发表相关资料,以帮助后人学习。 3. FlexPDE是一个非常有用的偏微分方程求解软件,而目前只有少数的研究者在用。展望 FlexPDE能普及到我们的学习与教学中。利于我们的学习和教学水平。 4.FlexPDE不仅可以解决量子力学问题,在其它学科有关偏微分方程问题中应该也可用到。 5.目前 FlexPDE软件可以图像形式显示结果,但缺乏数值解的功能,展望以后可以有这方面的改善。- 22 -参考文献 【1】 曾谨言 量子力学卷 第四版 科学出版社 2007 【2】 周世勋 量子力学教程 第二版高等教育出版社 2009 【3】