基本运算计算物理

基本运算计算物理第1页，课件共73页，创作于2023年2月本章内容插值14数值积分3数值微分2拟合45方程求根分子振动的半经典量子化46第2页，课件共73页，创作于2023年2月1.1插值xx0

x1

x2……

xnyy0

y1

y2

……

yn给定一组离散的数据寻找一个解析形式的函数φ(x)，满足φ(xi)=yi,i=1,2,…,n问题的提出最常用的函数形式是多项式，称为多项式插值。Y(x0,y0)(x1,y1)(xn,yn)x0x1xnX第3页，课件共73页，创作于2023年2月xx0

x1

x2yy0

y1

y2最直观的方法——以二阶为例设插值函数为带入数据，得离散数据为解这个方程，即可得到相应系数a0,a1,a2。第4页，课件共73页，创作于2023年2月N阶呢设插值函数为n阶多项式xx0

x1

……

xnyy0

y1

……

yn离散数据为带入数据，得当n很大时，求解这个方程计算量太大，需另寻它法第5页，课件共73页，创作于2023年2月拉格朗日插值——从一阶说起xx0

x1yy0

y1插值方程为一直线方程，可表示为对离散数据第6页，课件共73页，创作于2023年2月为了推广到高阶，将其写成更对称的形式x0x1A0

(x)10A1

(x)01函数值函数节点满足其中第7页，课件共73页，创作于2023年2月x0x1x2A0

(x)100A1

(x)010A2

(x)001函数值函数节点更进一步——二阶插值y1y2x0x1XYOy=f(x)x2y0第8页，课件共73页，创作于2023年2月这称之为拉格朗日多项式插值。一般的——N阶多项式插值xx0

x1

……

xnyy0

y1

……

yn已知离散数据插值多项式为其中满足第9页，课件共73页，创作于2023年2月过高阶的插值可能导致严重的振荡行为，即Runge现象。是否阶数越高，效果越好？怎样改进？-55例：连续函数可以看出，L10(x)的误差在区间两端非常大在区间[-5,5]上取等距插值节点讨论第10页，课件共73页，创作于2023年2月用分段低次插值，最简单的就是分段线性插值不光滑！解决的办法——分段插值第11页，课件共73页，创作于2023年2月区间[a,b]有离散点：a=x0<x1<…<xn=b，及其对应的函数值yi(i=0,1,…,n)，如果函数S(x)满足条件：(2)在每个子区间[xi

,xi+1](i=0,1,…,n-1)上是三次多项式则称S(x)是y=f(x)的三次样条插值函数确定S(x)需要4n个条件,而我们只给出了4n-2个条件端点函数值：2个内节点函数值及连续条件：2(n-1)

个内节点一、二阶导数连续条件：2(n-1)

个解决之道——在每个区间上，不是用线性函数而是用三阶多项式进行插值。(1)在[a,b]上具有二阶连续导数，即(3)样条插值第12页，课件共73页，创作于2023年2月需补充2个边界条件，根据实际情况有不同的取法例如，可取周期边界条件：第13页，课件共73页，创作于2023年2月样条插值的例子第14页，课件共73页，创作于2023年2月Matlab指令interp1(X,Y,Xi，method)X,Y 数据点及其函数值Xi 待求的数据点method 插值方法，可取’linear’,’spline’等返回值

数据点Xi的函数值第15页，课件共73页，创作于2023年2月1.2拟合给定一组离散数据：要求一个函数p(x)，按照最小二乘原理，使得达到最小值，这一方法被称为拟合，其中p(x)称为拟合曲线。xx1

……

xnyy1

……

yn第16页，课件共73页，创作于2023年2月则问题变为求ak(k=0,1,2,…,m),使得一种常见的拟合曲线为多项式拟合，即选取拟合曲线p(x)为m次多项式极小。第17页，课件共73页，创作于2023年2月由多元函数取极值条件，得这就是正则方程。解这个方程即可得到拟合函数的形式。第18页，课件共73页，创作于2023年2月当最高幂次为1时，即为最小二乘法的拟合直线，正则方程简化为第19页，课件共73页，创作于2023年2月Matlab指令polyfit(X,Y,N)X,Y 数据点及其函数值N 拟合多项式的最高幂次返回值

为按降幂排列的多项次系数第20页，课件共73页，创作于2023年2月例子i123456789xi-1-0.75-0.5-0.2500.250.50.751yi-0.22090.32950.88261.43922.00032.56453.13343.76014.2836设所求的最小二乘二次拟合多项式为相应的正则方程为第21页，课件共73页，创作于2023年2月其解为a0=2:0034;a1=2:2625;a2=0:0378，所求多项式为第22页，课件共73页，创作于2023年2月实例——密立根油滴实验第23页，课件共73页，创作于2023年2月指数拟合如果数据点的分布近似于指数曲线，则可考虑用指数函数去拟合数据。同时对数据yi

也取对数得Inyi,利用数据组(xi,Inyi)求出最小二乘拟合直线。再取指数即得到原数据组的最小拟合指数对上式取对数第24页，课件共73页，创作于2023年2月例子ti0.20.30.40.50.60.70.8Ii3.162.381.751.341.000.740.56已知一发射源的发射强度具有指数形式I=I0e-αt,现有一组观测数据如下第25页，课件共73页，创作于2023年2月WithfourparametersIcanfitanelephant,andwithfiveIcanmakehimwigglehistrunk.关于拟合，要提醒的第26页，课件共73页，创作于2023年2月1.3数值微分为什么要学习数值微分？我们碰到的函数往往没有解析形式，例如可能是如下数表形式x0.10.20.30.40.5f(x)44.5688.5这就需要借助于数值微分更重要的，数值微分是其它很多数值方法的基础第27页，课件共73页，创作于2023年2月求X=0处的导数x…-2h-h0h2h…f(x)…f-2f-1f0f1f2…已知第28页，课件共73页，创作于2023年2月方法：泰勒展开x=h处的函数值x=-h处的函数值第29页，课件共73页，创作于2023年2月两点公式向前差分向后差分第30页，课件共73页，创作于2023年2月三点公式（中心差分公式）第31页，课件共73页，创作于2023年2月进一步的改进——5点公式虽然精度更大，但是计算量也更大！通常的计算三点公式已足够。第32页，课件共73页，创作于2023年2月例子：不同微分方法的精度比较不同方法计算dsinx/dx|x=1=0.540302的误差第33页，课件共73页，创作于2023年2月不同公式的精度：五点>三点>两点注意误差随步长h的减小先减小再增大注意舍入误差微分公式涉及两个很接近的f值相减。当步长过小时，计算机的舍入误差会使导数计算不准确。第34页，课件共73页，创作于2023年2月高阶导数二阶导数更高阶的阶导数以此类推这也容易从二阶导数的定义直接得出第35页，课件共73页，创作于2023年2月1.4数值积分牛顿-莱布尼兹公式被积函数f(x)有解析表达式f(x)的原函数F(x)为初等函数但是这要求为什么要数值积分？第36页，课件共73页，创作于2023年2月1)f(x)没有解析表达式2)f(x)有解析表达式，但原函数不是初等函数，例如我们面临的的问题x0.10.20.30.40.5f(x)44.5688.5它们的原函数都不是初等函数第37页，课件共73页，创作于2023年2月数值积分的出发点积分转化为求和将区间[a,b]分割为n等份，每个小区间的宽度为h=(b-a)/n第38页，课件共73页，创作于2023年2月在每个区间[xi,xi+1]进行线性插值线性近似——梯形法则第39页，课件共73页，创作于2023年2月二阶多项式近似——辛普生法则在区间[xi-1,xi+1]上对f(x)进行二阶插值。第40页，课件共73页，创作于2023年2月不同积分方法的结果比较注意，此时误差随h减小而减小，舍入误差并不重要，这是因为积分公式中，所有f的值的符号都一样。第41页，课件共73页，创作于2023年2月选取更多的点，进行更高阶的插值可以得到更高阶的算法，如Bode规则（四阶插值）高阶的算法simpson3/8算法（三阶插值）第42页，课件共73页，创作于2023年2月过高阶的插值可能导致严重的振荡行为，从而给出被积函数不准确的插值。所以为了得到更高精度，往往用低阶插值，同时减小h。第43页，课件共73页，创作于2023年2月反常积分分为两类：积分区间有限,在积分区间内被积函数有奇点积分区间为无限反常积分的处理通过积分变量的变换，将反常积分变换为普通积分策略第44页，课件共73页，创作于2023年2月（1）

积分区间内含有奇点的积分在x=1

处有一个奇点，假设函数g

在这点的值有限，则积分为有限值。做变换t=(1-x)1/2，积分变为第45页，课件共73页，创作于2023年2月（2）无限区间的积分g(x)在x很大时趋于常数，积分为有限值。做变换t=x-1，积分变为第46页，课件共73页，创作于2023年2月quad

用自适应辛普森算法。根据积分精度的需要，自动调节积分取点的数目。调用格式为quad(函数句柄,积分下限，积分上限)quad(@(x)sin(x),0.5,0.6)Matlab自带的积分指令第47页，课件共73页，创作于2023年2月1.5数值求根高于四次方程的一般代数方程没有一般形式的代数解更不用说更为复杂的方程阿贝尔(1802—1829)第48页，课件共73页，创作于2023年2月求f(x)=0的根原理：若f

C[a,b]，且f(a)·f(b)<0，则f(x)=0在(a,b)上必有一根。yxbaf(x)x*1.5.1二分法（搜索法）第49页，课件共73页，创作于2023年2月给定有根区间[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*第50页，课件共73页，创作于2023年2月

简单易用

稳妥，对f(x)要求不高，只要连续即可收敛

收敛速度慢二分法的优缺点第51页，课件共73页，创作于2023年2月

例1:用二分法求方程在区间(1,2)内的实根,要求误差限为。二分法——例题第52页，课件共73页，创作于2023年2月

解：令

f(1)<0,f(2)>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.6758二分法——例题

例1:用二分法求方程在区间(1,2)内的实根,要求误差限为。第53页，课件共73页，创作于2023年2月xyx*x0几何意义1.5.2牛顿法迭代形式为第54页，课件共73页，创作于2023年2月[牛顿迭代法]1:初始化

x0,δ,置i:=02:如果|f(xi)|≤δ

，则停止.3:计算xi+1:=xi-f(xi)/f'(xi)4:如果|f(xi+1)|≤δ

，则停止.5:i:=i+1,转至3.牛顿法的算法构造第55页，课件共73页，创作于2023年2月例1：利用牛顿迭代法求解f(x)=ex-1.5-tan-1x的零点。初始点x0=-7.0

第56页，课件共73页，创作于2023年2月解：f(x0)=-0.702×10-1，f'(x)=ex-(1+x2)-1

计算迭代格式:

计算结果如下表：(取|f(x)|<=10-10)kxf(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-012例1：利用牛顿迭代法求解f(x)=ex-1.5-tan-1x的零点。初始点x0=-7.0

第57页，课件共73页，创作于2023年2月注：Newton’sMethod收敛性依赖于x0

的选取。x*x0x0x0算法说明第58页，课件共73页，创作于2023年2月x0x1切线

割线

切线斜率割线斜率任意2个初值x0

和x1可以启动这个递推关系。1.5.3弦割法第59页，课件共73页，创作于2023年2月三种求根算法的比较第60页，课件共73页，创作于2023年2月二分法最稳妥，但是效率最低。牛顿法效率最高，但是要求函数解析，容易计算导数。弦割法是前面两种方法的折衷，既有较高的效率，又不必像牛顿法那样必须计算函数f的导数。

如果初始猜测比较接近待求解，则其收敛速度几乎与牛顿算法一样快。三种求根算法的比较如果待求解附近，函数的行为不好，则自动的牛顿法和弦割法都可能无法收敛或收敛到错误的结果。保险的做法是先用二分法初步的定出解的位置，再用两个自动方法中的一个定出解的精确位置。第61页，课件共73页，创作于2023年2月fzero(@(x)sin(x),[3,4]）Matlab自带的求根指令X=fzero(函数句柄，猜测的初始值或搜寻的区间)fzero：求单变量函数的零点使用zeroin算法（结合了二分法、弦割法以及其它方法的一种综合方法）。fzero(@(x)sin(x),3）第62页，课件共73页，创作于2023年2月1.6分子振动的半经典量子化原子的相互作用势为Lennard-Jones或6-12形式势能最低处为rmin=21/6a，深度为V0第63页，课件共73页，创作于2023年2月能量为En的相对运动的振动态可以用一维薛定谔方程的束缚态解Ψn(r)来描述目标：对给定的势，求得能量En标准办法：数值求解常微分方程本征问题我们这里采用的方法：半经典量子化第64页，课件共73页，创作于2023年2月通过考虑原