FORTRAN数值方法及其在物理学中应用7节.ppt

1,第七章 实验物理学中的插值和数据拟合,2,学习意义,例如在大学物理实验中根据观测或实验得到一系列的 数据，确定了与自变量的某些点相应的函数值，而要 计算未测量到的点的函数值，为此可根据测量到的数 据构造一个适当的函数近似地去代替要寻求的函数， 这就是本章要介绍的插值和数据拟合在有关科学研究和工程计算中，常常会遇到计算函数 值等一类问题，然而函数关系往往是很复杂的，甚 至没有明显的解析表达式3,§7.1 实验数据的拉格朗日插值法,如已知,利用线性插值有：,若有一组实验观测结果,4,拉格朗日插值公式,如进行四点插值（三次拉格朗日插值多项式），有,5,通常采用三点插值（二次拉格朗日抛物线）插值中的位置6,,7,,8,例: 在用外接电源给电容器充电时，电容器两端的电压V将会随着充电时间t发生变化，已知在某一次实验时，通过测量得到下列观测值，试用拉格朗日插值法画出V随着时间t变化的曲线图，并计算当时间t=7s时，电容器两端电压的近似值利用上题计算程序，t=7.0代入计算可得：,V=9.708000,9,电容器电压随时间变化曲线,红色表示原始测量数据点，白色为插值结果10,解：取,代入拉格朗日（三点）二次插值多项式,11,对于某些函数，适当地提高插值多项式的次数，会提高 计算精度。

当函数 是连续函数时，加密插值节点虽然使插值函数 与被插值函数在更多节点上的值相等，但由于插值多项式 函数在某些非节点处的振荡可能加大，因而可能使在非节 节点处的误差变得很大讨论,节点加密会增加插值多项式计算次数，不利于控制舍入误差12,在接近区间两端点附近，两函数差别较大并非插值多项式的次数越高，其精度就越高高次多项式插值的Runge反例,13,作 业,14,§7.2 差商与牛顿插值公式,15,、差商概念,差商的写法与点 的排列顺序无关，具有对称性，即,16,差商的计算可列差商表：,17,,18,二、牛顿插值多项式,数多项式，设为：,由有关代数理论知，通过这,19,由方程组中第一式可得：,依此最后求得：,20,牛顿一次插值多项式为：,牛顿二次插值多项式为：,增加节点，以前计算结果仍可用，可以推出：,编程求牛顿插值多项式时，先造差商表，再按以下公式计算：,21,例5：给出通过例4中数据点,的牛顿三次插值多项式解：利用例4中的差商表及牛顿插值多项式公式得：,,22,作 业,23,§7.3 Hermite插值,插值函数曲线的光滑度可能很差，有时就不能满足实际,需要1),前面的拉格朗日插值具体实现时为分段低次插值法，它,无法保证分段插值多项式在节点处导数的连续性，因而,24,满足条件(1)的Hermite插值是最简单的情形。

值，并要求插值多项式满足,,分析,(2),性和稳定性都不能得以保证，所以采取分段插值法25,一、两点三次Hermite插值,(3),26,两点三次Hermite插值多项式,Step 1 :,27,Step 2 :,,(4)－(7)式得到,28,(9),(10),Step 3 :,29,类似地，设,(12),,（两点三次Hermite多项式）,30,解：设,因此有：,31,的Hermite插值公式为,(15),(17),Hermite插值有时只要求插值函数在个别节点上的导数取给定,的值，这时需要根据具体问题建立Hermite插值公式，而不是,说明,简单地套用公式(14)或(15)32,导数值如下表：,解：首先计算各次拉格朗日插值基函数及其导数：,,,,,,33,,,,,34,,,35,作 业,,,,,,,,,,,,,36,二、分段两点三次Hermite插值,37,,38,(19),(20),39,(21),(22),(23),(24),这种性质为局部非零性40,基函数的线性组合：,(25),41,用分段三次Hermite插值多项式计算,的近似值，并与准确值相比较42,可以证明分段三次Hermite插值多项式函数在插值区间上,计算结果表明用分段三次Hermite插值多项式可较好地,近似,43,作 业,,,,,,,,,,,,44,§7.4 三次样条插值,一、三次样条函数,实际问题中提出的插值问题，有一些要求插值函数曲线,例：根据风洞实验，可以构造出一种具有所需要特征的机翼，,具有较高的光滑性。

45,说明,为构造具有二阶光滑度的插值函数，引进三次样条,只有一阶光滑度定义4,满足：,函数概念46,二、三次样条插值多项式,及节点上的函数值：,(1),(2),(3),47,(4),表示导数阶数即所谓满足边界条件通常使用的边界条件有以下三类：,,实际问题对三次样条插值函数在端点处的状态也有要求，,第一类边界条件:,(5),48,第二类边界条件：,(6),时，样条函数在两端点不受力，呈自然状态，,称为自然边界条件第三类边界条件：周期性条件,(7),条件,49,,分析,个数50,(9),51,(11),(12),分析,然后求得：,得到：,52,(13),(15),的方程组(13)，其中每一个方程都联系着,在相邻三个节点上一阶导数值，称三对角方程,讨论,各种边界条件的情形－基本方程组,53,,(16),系数矩阵都是严格对角占优的三对角线矩阵，追赶法可求唯一解,54,,得,55,用矩阵表示为,(20),(15),系数矩阵都是严格对角占优的三对角线矩阵，追赶法可求唯一解,56,由条件(7)式，,由(10)和(11)式得,,,(23),57,将方程(23)与(13)联立，,,(26),,,系数矩阵也是严格对角占优的，所以非奇异，方程组有唯一解。

58,例9：给定函数表,满足第二类边界条件，相应的方程组为(20)未,该问题需,59,(2).计算,计算结果如下表：,(3).求解方程组,将上面求得数据代入方程组，求得,60,,61,表示为分段三次多项式，即,,,62,曲线图形,63,作 业,64,§7.5 数值微分,一、插值型求导公式,导数值,的近似值：,65,插值型求导公式,,时应特别注意误差分析故上式对应的导数误差为：,误差为：,66,,是无法估计的等号右边第二项变为零，此时余项公式：,67,下面仅考察节点处的导数值为简化讨论，假定所给的,1. 两点公式,即,截断误差,68,2. 三点公式,69,,还可以建立高阶数值微分方程公式,,70,3.五点公式,同三点公式类似，可以导出五点公式：,相邻的节点，若一侧的节点不是两个（即一个或没有），,依次为,则用另一侧节点补足71,解：,节点及其对应的导数值的计算结果列表如下：,72,计算程序,real*8 x(6),y(6),d(6) data x/0.1,0.2,0.3,0.4,0.5,0.6/ data y/1.2051709,1.4214028, 1.6498588, * 1.8918247, 2.1487213, 2.4221188/ h=0.1 N=6 d(1)=(-3.*y(1)+4*y(2)-y(3))/(2*h) write(*,*) x(1),d(1) do 10 i=2,N-1 d(i)=(-y(i-1)+y(i+1))/(2*h) 10 write(*,*) x(i),d(i) d(N)=(y(N-2)-4*y(N-1)+3*y(N))/(2*h) write(*,*) x(N),d(N) end,73,解：解析解为：,数值解为：,。

计算结果比较列于下表：,74,二、样条求导公式,利用拉格朗日插值多项式导出的数值微分公式只能求节点上,问题：,非节点处的导数为解决这一问题，可采用三次样条插值多项式代替,75,有时将三点或五点公式与三次样条求导公式配合使用，可以,先用比较简单的三点或五点公式求出等间距节点处一阶导数，,然后再利用三次样条求导公式求出非节点处的一阶或二阶数,值导数，这样可以回避解三对角方程说明,76,,作 业,EX7-8: 写出三点求导公式的误差,77,78,§7.6 最小二乘曲线拟合,一、最小二乘法的一般原理,１、数据量较大目的,数据特点,２、数据是通过测量得到，数据本身有一定的误差79,分析,这样可抵消原数据组中测量误差若用插值法，通过这几个已知点所求得的插值多项式必定,是高次插值多项式，而高次插值是数值不稳定的由于数据本身存在的误差，利用插值所得到的插值多项式,必保留了所有测量误差，所得结果与实际问题误差较大对数据拟合问题，一般地讲，并不要求所得到的近似解析,表达式通过所有已知点，而只要求尽可能通过它们附近，,80,问题,最小二乘法（最小二乘曲线拟合法）81,拟合函数,合而成，即,82,二、用最小二乘法求解矛盾方程组,求解线性方程组时，通常要求未知数的个数与方程式个,(1),数相等。

若方程式的个数多于未知数的个数，方程无解，,称为矛盾（超定）方程组该线性方程组在实验中对应于,即,83,矛盾方程组求解过程,寻求各未知数的一组值，使方程组中各式能近似相等最小二乘法的基本思想：矛盾方程组无精确解，只能寻求某,种意义下的近似解，这种近似解并非对精确解之近似，而是,按最小二乘法原则，求各个方程式误差平方和,这组值是矛盾方程组的最优近似解令,84,而,,,（极值条件）,85,方程组(1)的正规方程组显然，(2)的解是(1)的最优近似解记,，,(3),(2),,86,令,(2)式可写为：,为对称阵 2)式又可写为：,(4),(2),87,２、求解正规方程组，得矛盾方程组的最优近似解最小二乘法解矛盾方程组步骤,88,例12：用最小二乘法求矛盾方程组,的最优近似解因此有：,，,解正规方程组：,解：由题知：,89,作 业,90,三、用多项式作最小二乘曲线拟合,取基函数为：,则拟合多项式为：,分析,91,(5),即,它对应的正规方程组为：,92,求得矛盾方程组的唯一一组最优近似解，应使,取极小，从而求得所给数据的最小二乘拟合多项式对称正定）,93,只需计算,和,1. 计算正规方程组的系数矩阵和常数项各元素,2. 利用迭代法求正规方程组的解,最小二乘数据拟合步骤：,,94,例13：通过实验获得如下数据,用最小二乘法求多项式曲线，使与此数据相拟合。

解：取,建立正规方程组，首先计算,,,95,实验数据和拟合多项式曲线,96,dimension x(7),y(7),a(3) data x/1,2,3,4,6,7,8/ data y/2,3,6,7,5,3,2/ N=7 M=3 call LSM(x,y,a,N,M) write(*,*) write(*,20) (I,a(I),I=1,M) 20 format(1X,‘a(',I2,' )=',D15.6) write(*,*) end,最小二乘法求解正规方程组系数c,d ! subroutine LSM(x,y,a,N,M) dimension x(N),y(N),a(M),c(M,M),d(M) do 15 i=1,m do 10 j=1,m c(i,j)=0.0 d(i)=0.0 do 40 k=1,n c(i,j)=c(i,j)+x(k)**(i+j-2) 40 d(i)=d(i)+y(k)*x(k)**(i-1) 10 continue 15 continue write(*,*) c write(*,*) d call agaus(c,d,m,a) return end,计算程序,97,试用最小二乘法求最佳数据拟合函数。

解：将数据经过绘制，它近似为一条指数曲线,两边取常用对数得,即,因此取指数曲线：,测量数据,98,得正规方程组：,,测量数据和拟合曲线,99,,,作 业,100,作 业,,101,本章小结,一、拉格朗日插值法,点插值公式,三点公式及计算程序,102,二、牛顿插值法,构造并计算差商表,103,三、Hermite插值,两点三次Hermite插值,分段两点三次Hermite插值,104。