第一章基本运算

1、第一章 基本数学运算本章要学会计算的物理问题：双原子分子的振动能级本章内容插值14数值积分3数值微分2拟合45方程求根分子振动的半经典量子化461.1 插值x x0 x1 x2 xny y0 y1 y2 yn给定一组离散的数据寻找一个解析形式的函数 (x)， 满足 (xi) = yi, i=1,2, n 问题的提出最常用的函数形式是多项式，称为多项式插值。Y(x0,y0)(x1,y1)(xn,yn)x0 x1xnXx x0 x1 x2y y0 y1 y222102)(xaxaax222221012121100202010yxaxaayxaxaayxaxaa最直观的方法以二阶为例设插值函数为带入

2、数据，得离散数据为解这个方程，即可得到相应系数 a0, a1, a2。N阶呢nnnnnnnnnnyxaxaxaayxaxaxaayxaxaxaa22101121211000202010 nnnxaxaxaax2210)(设插值函数为n阶多项式x x0 x1 xny y0 y1 yn离散数据为带入数据，得当n很大时，求解这个方程计算量太大，需另寻它法拉格朗日插值从一阶说起x x0 x1y y0 y1插值方程为一直线方程，可表示为对离散数据为了推广到高阶，将其写成更对称的形式x0 x1A0 (x)10A1 (x)01函数值函数值函数函数节点节点1010 xxxxA0101xxxxA满足其中0122

3、021()()()()xxxxAxxxx221100)()()()(yxAyxAyxAxyx0 x1x2A0 (x)100A1 (x)010A2 (x)001函数值函数值函数函数节点节点更进一步二阶插值y1y2x0 x1XYOy=f(x)x2y01200102()()()()xxxxAxxxx2011210()()()()xxxxAxxxxjnjjyxAxy)()(0njiiijijxxxxxA0)(这称之为拉格朗日多项式插值。一般的N阶多项式插值x x0 x1 xny y0 y1 yn已知离散数据插值多项式为其中满足过高阶的插值可能导致严重的振荡行为，即Runge现象。是否阶数越高，效果越好

4、？怎样改进？)(xf)(10 xL-5 5例：连续函数215)(xxf可以看出，L10(x)的误差在区间两端非常大在区间-5,5上取等距插值节点讨论用分段低次插值，最简单的就是分段线性插值不光滑！解决的办法分段插值区间a,b 有离散点：a= x0 x1 三点两点 注意误差随步长 h 的减小先减小再增大注意舍入误差微分公式涉及两个很接近的 f 值相减。当步长过小时，计算机的舍入误差会使导数计算不准确。高阶导数二阶导数更高阶的阶导数以此类推这也容易从二阶导数的定义直接得出1.4 数值积分牛顿-莱布尼兹公式被积函数 f(x) 有解析表达式f(x) 的原函数 F(x) 为初等函数)()()(aFbFd

5、xxfIab但是这要求为什么要数值积分？1) f(x) 没有解析表达式2) f(x)有解析表达式，但原函数不是初等函数 ，例如我们面临的的问题x0.10.20.30.40.5f(x)44.5688.5它们的原函数都不是初等函数数值积分的出发点积分转化为求和将区间 a,b 分割为 n 等份，每个小区间的宽度为 h=(b-a)/n在每个区间 xi, xi+1 进行线性插值线性近似梯形法则二阶多项式近似辛普生法则在区间 xi-1, xi+1 上对 f(x) 进行二阶插值。不同积分方法的结果比较注意，此时误差随h减小而减小，舍入误差并不重要，这是因为积分公式中，所有 f 的值的符号都一样。选取更多的点

6、，进行更高阶的插值可以得到更高阶的算法，如 Bode规则（四阶插值） 高阶的算法simpson 3/8 算法（三阶插值）过高阶的插值可能导致严重的振荡行为，从而给出被积函数不准确的插值。所以为了得到更高精度，往往用低阶插值，同时减小 h。反常积分分为两类： 积分区间有限,在积分区间内被积函数有奇点 积分区间为无限反常积分的处理通过积分变量的变换，将反常积分变换为普通积分策略（1） 积分区间内含有奇点的积分积分区间内含有奇点的积分在 x=1 处有一个奇点，假设函数 g 在这点的值有限，则积分为有限值。做变换 t=(1-x)1/2，积分变为（2）无限区间的积分）无限区间的积分g(x) 在 x 很大

7、时趋于常数，积分为有限值。做变换 t=x-1，积分变为 quad 用自适应辛普森算法。根据积分精度的需要，自动调节积分取点的数目。 调用格式为 quad(函数句柄, 积分下限，积分上限)quad(x)sin(x), 0.5, 0.6)Matlab自带的积分指令1.5 数值求根高于四次方程的一般代数方程没有一般形式的代数解更不用说更为复杂的方程1110( )(0)nnnnnf xa xaxa xaa( )( ) 10cos x cosh x 阿贝尔阿贝尔(18021829)求 f (x) = 0 的根原理：若 f Ca, b，且 f (a) f (b) 0，则 f (x)=0 在 (a, b)

8、上必有一根。yxbaf (x)x*1.5.1 二分法（搜索法）给定有根区间 a, b ( f(a) f(b) 0) 和 精度 1. 令 x = (a+b)/22. 如果 b a ， 停机，输出 x3. 如果 f (a) f (x) 0 , 则令 b = x，否则令 a = x, 返回第1步二分法的算法实现abx1x2abx* 简单易用 稳妥，对 f (x) 要求不高，只要连续即可收敛 收敛速度慢 二分法的优缺点 例1: 用二分法求方程 在区间 (1,2)内 的实根, 要求误差限为 。01523 xx21011.11.21.31.41.51.61.71.81.92-4-3-2-10123453(

9、 )251f xxx二分法例题 解：令 f (1)0 记 I0=1,2 , x0 =(1+2)/2=1.5 因为 f (x0) f (1)0 得 I1=1.5, 2 , x1 =(1.5+2)/2=1.75 f (x1) f (1.5)0 得 I2=1.5, 1.75 , x2 =(1.5+1.75)/2=1.625 . I6=1.681875, 1.6875, I7=1.671875, 1.679688 b7 - a7=0.781310-2 10-2 x*x7 =1.6758152)(3xxxf二分法例题 例1: 用二分法求方程 在区间 (1,2)内 的实根, 要求误差限为 。01523 x

10、x210 xyx*x0)()(1kkkkxfxfxx 0100()()f xxxfx 几何意义1.5.2 牛顿法迭代形式为牛顿迭代法1: 初始化 x0 , ,置 i:=02: 如果| f(xi ) | ，则停止. 3: 计算 xi+1:=xi - f (xi) / f (xi)4: 如果 | f (xi+1) | ，则停止.5: i:=i+1, 转至3.1()()iiiif xxxfx 牛顿法的算法构造例1： 利用牛顿迭代法求解 f(x)=ex-1.5-tan-1x 的零点。初始点 x0=-7.0 解： f (x0)=-0.70210-1，f (x)=ex-(1+x2)-1 计算迭代格式: 计

11、算结果如下表：(取|f(x) |=10-10)k x f (x) 0 -7.0000 -0.07018881 -10.6771 -0.02256662 -13.2792 -0.004366023 -14.0537 -0.000239024 -14.1011 -7.99585e-0075 -14.1013 -9.00833e-0121()()kkkkf xxxfx 例1： 利用牛顿迭代法求解 f(x)=ex-1.5-tan-1x 的零点。初始点 x0=-7.0 Newtons Method 收敛性依赖于x0 的选取。x*x0 x0 x0算法说明x0 x1切线 割线 切线斜率割线斜率11)()()

12、( kkkkkxxxfxfxf)()()(111 kkkkkkkxfxfxxxfxx任意2个初值 x0 和 x1可以启动这个递推关系。1.5.3 弦割法三种求根算法的比较 二分法最稳妥，但是效率最低。 牛顿法效率最高，但是要求函数解析，容易计算导数。 弦割法是前面两种方法的折衷，既有较高的效率，又不必像牛顿法那样必须计算函数 f 的导数。 如果初始猜测比较接近待求解，则其收敛速度几乎与牛顿算法一样快。 三种求根算法的比较如果待求解附近，函数的行为不好，则自动的牛顿法和弦割法都可能无法收敛或收敛到错误的结果。保险的做法是先用二分法初步的定出解的位置，再用两个自动方法中的一个定出解的精确位置。fz

13、ero(x)sin(x),3,4）Matlab自带的求根指令X=fzero(函数句柄，猜测的初始值或搜寻的区间)fzero：求单变量函数的零点使用zeroin算法（结合了二分法、弦割法以及其它方法的一种综合方法）。fzero(x)sin(x),3）1.6 分子振动的半经典量子化原子的相互作用势为 Lennard-Jones 或6-12形式势能最低处为rmin=21/6a，深度为V0能量为En的相对运动的振动态可以用一维薛定谔方程的束缚态解n(r)来描述目标：对给定的势，求得能量En标准办法：数值求解常微分方程本征问题我们这里采用的方法：半经典量子化通过考虑原子核在势场中做经典运动，然后应用量子

14、化规则，可以定出其振动能 En。这就是所谓的Bohr-Sommerfeld-Wilson量子化法则。半经典量子化Bohr 根据对应原理的思想得出了一个角动量子化的条件，即电子运动的角动量 J 只能是 的整数倍Bohr的量子化法则Somemerfeld 等为处理多自由度体系的周期运动的能量量子化，给出了推广的量子化条件其中其中 qk, pk 代表一对共轭的正则坐标与动量， 代表对周期运动积分一个周期。推广的量子化法则在相空间轨迹为势场 V 中，原子核间距的经典的运动，可以在能量 -V0E0上发生。原子核间距在 rin 和 rout 之间周期性的振动，对应相空间中的一条封闭轨道，轨道方程为其中n为非负整数量子化条件为即其中为了对 Lennard-Jones势确定量子化条件，定义几个无量纲量量子化条件变为是标度化后的位势。无量纲化分子分子H2HDO221.724.8150若知道分子的转动惯量（可以从分子转动的能量得知）和离解能（把分子分解为组成它的两个原子所需的能量），就能从实验观测测定参量 a 和 V0，从而定出 。问题1假设原子间相互作用势为谐振子势，利用半经典量子化条件，解析求系统能级。其中问题2证明：L-J势下，小幅振动的频率之预期值为问题3对氢分子更适用的势函数为调