《计算物理物理学中常用的数值方法.pptx》由会员分享，可在线阅读，更多相关《计算物理物理学中常用的数值方法.pptx（41页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 物理学中经常遇到求单变量函数方程：f f（x x）0 0的根的问题（1 1）可严格求解（解析求解）（2 2）不可严格求解第1页/共41页可严格求解可严格求解可严格求解可严格求解 例例1 1：在热力学中，求最可几速率（最概然速率）已知气体分子的速率分布，即单位体积内速率 内的分子数：其中：n 单位体积内的分子数 m 分子质量 k 玻尔兹曼常数 T 温度 麦克斯韦速率分布式第2页/共41页可严格求解可严格求解可严格求解可严
2、格求解最可几速率（最概然速率）：速率 附近的单位速率间隔内的分子数最多。其可如下求出：化简可得此方程可以严格求解：第3页/共41页 例例2 2：单缝衍射问题 已知夫琅和费单缝衍射的光强分布为：除了中央极大外，另外的光强极大的位置由下式得出：化简可得：只能用数值法求解不可严格求解不可严格求解不可严格求解不可严格求解第4页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 实际生活中的二分法：在一个风雨交加的夜里，从某水库闸房到防洪指挥部的电话线路发生了故障，这是
3、一条10km长的线路，如何迅速查出故障所在？要把故障可能发生的范围缩小到50100m左右，即一两根电线杆附近，要检查多少次？方法分析：方法分析：定义定义：每次取中点，将区间一分为二，再经比较，按需要留下其中一个小区间的方法叫二分法，也叫对分法。常用于：查找线路电线、水管、气管等管道线路故障；实验设计、资料查询；是方程求根的常用方法！七次七次第5页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步基本概念：根的存在定理 假设函数y=f(x),xa,b,且f(a)
4、f(b)0,函数图象如图则至少存在一点x a,b 使得f(x)=0,这就是根的存在定理。二分法就是将方程的有根取间对分,然后在选择比原间缩小一半的有根区间,如此继续下去,直到得到满足精度要求的根为止的一种简单的区间方法。yxba第6页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步原理分析：给定方程f(x)=0,设f(x)在区间a,b连续,且f(a)f(b)0,则方程f(x)在(a,b)内至少有一根,为便于讨论,不妨设方程f(x)=0在(a,b)内只有一实根
5、。采取使有根区间逐步缩小,从而得到满足精度要求的实根的近似值 。取a,b区间二等分的中点x0=(a+b)/2,若f(x0)=0,则x0是f(x)=0的实根 若f(a)f(x0)0 成立,则实根必在区间(a,x0)内,取a1=a,b1=x0;否则 必在区间(x0,b)内,取a1=x0,b1=b,这样,得到新区间(a1,b1),其长度为a,b的一半,如此继续下去,进行k次等分后,得到一组不断缩小的区间,a,b,a1,b1,.ak,bk.第7页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortr
6、an程序设计初步程序设计初步其中每一个区间都是前一个长度的一半,从而 ,的长度为 如此继续下去,则有这些区间将收敛于一点 ,该点即为所求的根.当做到第k步时，有 为给定精度 此时 即为所求方程的近似解.以上方法就是用于求非线性方程实根近似值的二分法第8页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 基本步骤：给定精确度 ，则（1 1）先猜a a，b b两个值以及精确度，使得f(a)*f(b)f(a)*f(b)小于0 0，也就是f(af(a）和 f(b)f
7、(b)必须异号。若f(x)是连续函数，则可保证在a,b 区间存在一个c值，使得f(c)=0。（2）令c=(a+b)/2，如果f(c)=0就找到了一个解，工作完成。（3）f(c)不为0时，如果f(a)，f(c)异号，则以a，c为新的两个猜测值来重复步骤（2）；如果f(b)，f(c)异号，b，c为新的值来重复步骤（2）。根据步骤编写程序 第9页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 收敛性分析：现在来研究用二分法求 函数的根时的精确性。假定f(x)是连
8、续函数,并且它在区间a0.b0的两端点 所取的值有相反的符号,于是在a0.b0中有一个根r,如果用中点c0=(a0+b0)/2作为对r的估计则有 r-c0(b0-a0)/2.如图所示:r-c0 a0 r c0 b0现在应用二分算法,并将算出的量用a0,b0,c0,a1,b1,c1 等来表示,则由同样的推理,r-cn(bn-an)/2 n=0,1,2,.由于这些区间的宽度在每一步中都除以2,所以可断定 r-cn 可见经过n步后将算出一个近似根其误差至多为 （b-a）/2n+1 第10页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截
9、法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 优缺点：u二分法计算过程简单,程序容易实现.可在大范围内求根。但该方法收敛较慢,且不能求偶数重根和复根,一般用于求根的初始近似值,而后在使用其它的求根方法。二分法收敛速度不快,其收敛速度仅与一个以 1/2为比值的等比级数相同。第11页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 牛顿法基本思想：牛顿法是求解单变量函数方程的一种重要方法。这种方法的基本思想是将非线性方程f(x)=0
10、f(x)=0近似化为线性方程求解 具体推导过程第12页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 牛顿法几何意义：由f(xf(xn+1n+1)=f(x)=f(xn n)+f)+f(x(xn n)(x-x)(x-xn n)=0)=0知，x xn+1n+1是（x xn n,f(x,f(xn n)）的切线 与X X轴的交点的横坐标（如图）。也就是说，新的近似值x xn+1n+1是用代替曲线y=f(x)y=f(x)的切线与x x 轴相交得到的。继续取点（x x
11、n n1 1,f(x,f(xn n1 1)）,再做切线与x x轴相交，又可得x xn n2 2 。由图可见，只要初值取的充分靠近 ,这个序列就会很快收敛于与x x的交点。uNewtonNewton迭代法又称切线法第13页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 牛顿法基本步骤：1.准备。选定初始近似值x0，计算出f0=f(x0);f0=f(x0)2.迭代。按公式 x1=x0-f0/f0 迭代一次，得到新的近似值x1，计算f1=f(x1);f1=f(x
12、1)。3.控制。如果x1满足 ，则终止迭代，以 x1 作为所求的根；否则，转下一步。4.修改。如果迭代次数达到预先指定的次数N，或者f0=0,则方法失败，否则以（x1,f1,f1）代替（x0,f0,f0）转步骤2继续进行迭代。程序示例(exp(x)lnx-x2=0)第14页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 牛顿法优缺点：优点：牛顿迭代法具有平方收敛的速度，所以在迭代过程中只要迭代几次就会得到很精确的解。这是牛顿迭代法比简单迭代法优越的地方 缺
13、点：选定的初值要接近方程的解，否则有可能的不到收敛的结果。再者，牛顿迭代法计算量比较大。因每次迭代除计算函数值外还要计算微商值如果f f0 0=0 =0，应该如何处理？第15页/共41页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 弦截法：弦截法的迭代公式：其几何意义为：设曲线y=f(x)上横坐标为x0和xk的点分别为P0 和Pk，则差商 表示弦P0Pk的斜率，弦的方程为：因而此法可形象地称为弦截法Ox*xk+1xkPkP0 x0 xyy=f(x)第16页/共41
14、页洛阳师范学院物理系2.1 2.1 函数方程求根函数方程求根u引言引言u二分法二分法u牛顿法牛顿法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 弦截法：弦截法的收敛速度比牛顿法慢得多。为了加快收敛速，我们改用差商 来代替牛顿迭代格式中的一阶导数 于是我们就可以得到快速弦截法的迭代公式：快速弦截法也是迭代法。不过，前面讨论的迭代法在计算xk+1时只用到上一步的结果xk（一步迭代）。而快速弦截法则是一种多步迭代（两步迭代：在计算xk+1时要用到前两步的结果xk和xk-1，因此在使用快速弦截法时，必须给出两个初始近似根x0和x1）。程序示例第17页/共41页洛阳师
15、范学院物理系2.2 2.2 定积分计算定积分计算u问题提出问题提出u梯形法梯形法u辛普森法辛普森法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初步 问题的提出：计算定积分的方法：计算定积分的方法：(1)求原函数；问题：问题：(1)被积函数的原函数不能用初等函数表示；(2)被积函数难于用公式表示，而是用图形或表格给出的；(3)被积函数虽然能用公式表示，但计算其原函数很困难(2)利用牛顿莱布尼茨公式得结果第18页/共41页洛阳师范学院物理系u问题提出问题提出u梯形法梯形法u辛普森法辛普森法u弦截法弦截法第一章第一章 FortranFortran程序设计初步程序设计初
16、步解决办法解决办法：建立定积分的近似计算方法常用方法常用方法：矩形法、梯形法、辛普森法思路：思路：2.2 2.2 定积分计算定积分计算第19页/共41页洛阳师范学院物理系u问题提出问题提出u矩形法矩形法u梯形法梯形法u辛普森法辛普森法第一章第一章 FortranFortran程序设计初步程序设计初步则有矩形法（1）2.2 2.2 定积分计算定积分计算第20页/共41页洛阳师范学院物理系2.2 2.2 定积分计算定积分计算u问题提出问题提出u矩形法矩形法u梯形法梯形法u辛普森法辛普森法第一章第一章 FortranFortran程序设计初步程序设计初步矩形法（2）则有第21页/共41页洛阳师范学院
17、物理系2.2 2.2 定积分计算定积分计算u问题提出问题提出u矩形法矩形法u梯形法梯形法u辛普森法辛普森法第一章第一章 FortranFortran程序设计初步程序设计初步梯形法梯形法就是在每个小区间上，以窄梯形的面积近似代替窄曲边梯形的面积，如图在实际的计算中，公式可以写为第22页/共41页洛阳师范学院物理系2.2 2.2 定积分计算定积分计算u问题提出问题提出u矩形法矩形法u梯形法梯形法u辛普森法辛普森法第一章第一章 FortranFortran程序设计初步程序设计初步辛普森法：所谓的辛普森方法，就是用一个二次曲线来近似一个子区间上的被积函数。把积分空间a,b分为两个子区间，分点为(a,a
18、+h,a+2h=b),h=(b-a)/2，在此区间上，把被积函数用二次函数来近似，即取：具体推导公式第23页/共41页洛阳师范学院物理系2.3 2.3 特殊函数计算特殊函数计算u勒让德勒让德u拉盖尔拉盖尔第一章第一章 FortranFortran程序设计初步程序设计初步Legendre Polynomials：勒让德多项式 Pl(x)满足递推公式：而且已知：P0(x)=1 P1(x)=x程序编写第24页/共41页洛阳师范学院物理系2.3 2.3 特殊函数计算特殊函数计算u勒让德勒让德u拉盖尔拉盖尔第一章第一章 FortranFortran程序设计初步程序设计初步Laguerre Polynom
19、ials：勒让德多项式 Pl(x)满足递推公式：而且已知：L0(a)=1L1(a)=a+1-x程序编写第25页/共41页subroutine glp(alpha,n,x,gl)!generalized Laguerre polynomial!implicit real*8(a-h,o-z)dimension gl(0:n)gl(0)=1.d0gl(1)=alpha+1.d0-xdo 10 i=2,ngl(i)=(2*i+alpha-1.d0-x)*gl(i-1)-(i+alpha-1.d0)*gl(i-2)/dble(i)10continuereturnend 第26页/共41页Program
20、 LegendreDIMENSION P(0:20)lmax=5x=0.5CALL land(lmax,x,P)WRITE(*,*)lmax,x,P(lmax)ENDSUBROUTINE lgnd(lmax,x,P)DIMENSION P(0:lmax)P(0)=1.0P(1)=xDO L=1,lmax P(L+1)=(2.0*L+1)*x*P(L)-L*P(L-1)/(L+1)END DOEND 第27页/共41页把三点的坐标带入二次函数式中：由此可得：辛普森求积公式。由于其使用二次曲线近似，所以它具有二次代数精度，即二阶精度。第28页/共41页 在实际的计算中，我们会把区间a,b分成N个子
21、区间，然后在每个子区间中用辛普森求积公式计算，最后再把结果加起来 这是一个在实际计算中常常被采用的方法，它对于一般有限区间的积分都能给出较好的正确结果。第29页/共41页梯形法实际上是把每个子区间上的被积函数近似为线性函数然后计算其积分值，再把结果加起来。在这个子区间线性函数与被积函数有两个交点，所以这种方法具有一阶精度如果我们需要更高的精度，怎么办？第30页/共41页 PROGRAM cha_3 dl=1.0E-06;a=1.0;b=2.0;dx=(b-a)/10.0 x0=(a+b)/2.0 CALL secant(dl,x0,dx,istep)WRITE(*,999)istep,x0,d
22、x999 FORMAT(I4,2F16.8)END SUBROUTINE secant(dl,x0,dx,istep)Istep=0 x1=x0+dx DO WHILE(ABS(dx).GT.dl)d=f(x1)-f(x0)x2=x1-f(x1)*(x1-x0)/d x0=x1 x1=x2 dx=x1-x0 istep=istep+1 END DO END FUNCTION f(x)f=exp(x)*alog(x)-x*x END 把二分法和牛顿法（弦截法）结合起来，会扩大算法的适用范围第31页/共41页注：注：注：注：Newton法的收敛性依赖于法的收敛性依赖于x0 的选取。的选取。x*x0
23、 x0 x0第32页/共41页 PROGRAM Newton dl=1.0E-06;a=1.0;b=2.0;dx=b-a x0=(a+b)/2.0 Istep=0 DO WHILE(ABS(dx).GT.dl)x1=x0-f(x0)/df(x0)dx=x1-x0 x0=x1 istep=istep+1 END DO WRITE(*,999)istep,x0,dx999 FORMAT(I4,2F16.8)END FUNCTION f(x)f=exp(x)*alog(x)-x*x END FUNCTION df（x）df=exp(x)*(alog(x)+1.0/x)-2.0*xexp(x)*(al
24、og(x)+1.0/x)-2.0*x END END第33页/共41页是允许的误差C是取绝对值或者相对误差的控制常数，一般可取C1第34页/共41页第35页/共41页设 是f(x)=0的一个近似根，把f(x)在 处作泰勒展开若取前两项来近似代替f(x)(称为f(x)的线性化)，则得近似的线性方程设 ,则其解为 当然这只是一个近似解，因为我们只取了级数展开的前两项。可以把该x作为新的近似根，记为x1，根据同样的过程可以找到更精确的近似解x2，一般有迭代公式：只要f(xn)与零足够接近，便认为xn为函数方程的解牛顿迭代公式第36页/共41页 二分法程序示例二分法程序示例二分法程序示例二分法程序示例
25、PROGRAM ch3_1USE numericalIMPLICIT noneREAL a,b !两个猜值REAL ans !算出的值DO WHILE(.TRUE.)WRITE(*,*)”输入两个猜测值”READ(*,*)a,b IF (f(a)*f(b)0)EXIT !f(a)*f(b)zero)!f(c)小于zero时，就视f(c)=0，结束循环 fa=func(a)fb=func(b)IF(fa*fc0)THEN !F(a)*f(c)0，以a,c值为新的区间 b=c c=(a+b)/2.0ELSE a=c c=(a+b)/2.0END IF fc=func(c)END DObisect=cRETURNEND function第39页/共41页REAL function f(x)IMPLICIT noneREAL xf=(x+3)*(x-3)RETURNEND functionEND MODULE第40页/共41页洛阳师范学院物理系感谢您的观看。第41页/共41页