地球物理计算方法课件：第一章_插值方法 4

地球物理计算方法地球物理与信息技术学院地球物理与信息技术学院随着节点的增加，采用高次多项式插值，可以在某些区域较好的逼近原来的函数（如在-2,2区间）；但在高次多项式的两端出现了激烈震荡的现象，这就是所谓的龙格现象。104.81.8044P4.80.0416f高次高次多项式插值多项式插值的的病态病态性性 复习复习3分段线性插值公式：11111iiiiiiiixxxxSxyyxxxx，1iixxx再把所有区间上的插值函数联合起来。在几何上就是用折线替代曲线,如右图所示复习复习4分段三次插值函数S(x)整体上具有一阶连续导数，满足下述条件：()iiS xy()iiS xy),2,1,0(ni(1)各区间插值节点上(2)在每个小区间xi,xi+1上，都是一个三次多项式：230123()iiiiiS xaa xa xa x复习复习5由埃尔米特插值公式，可直接写出分段三次Hermite插值函数的分段表达式：2211111112211111()1212()()iiiiiiiiiiiiiiiiiiiiiiiixxxxxxxxS xyyxxxxxxxxxxxxxxyxxyxxxx复习复习1130()0S(),S(),hR xxC a bxC a b 当当时时，。但但分分段段低低次次插插值值整整体体光光滑滑性性低低，整整体体连连续续，但但一一阶阶导导数数不不连连续续，整整体体一一阶阶导导数数连连续续，但但二二阶阶导导数数不不连连续续.212201,1|()|()(),8max,ma|()1)|(xii nxa bR xf xSxM hxa bhhMfx 分分段段线线性性插插值值的的误误差差4344(4)401,(2)1|()|()()|4!2max,max|()|ii nxa bR xf xSxM hhhMHerftexmi 分分段段插插值值的的误误差差复习复习7对于区间a,b上n+1个节点样本点数据),(,),(),(1100nnyxyxyx构造一个插值函数S(x)满足:iiyxS)(),2,1,0(ni,(1)(2)在每个小区间xi,xi+1上，都是一个三次多项式：332210)(xaxaxaaxSiiiii(3)(),(),(xSxSxS 在a,b上连续；S(x)是一个光滑的分段函数，这样的函数称为三次样条（Spline）插值函数。333333S(0)S(0)S(0)S(0),1,2,1S(0)S(0)iiiiiixxxxinxx 复习复习8为使方程组有解，在区间a,b的两个端点上各加一个条件，即称之为边界条件。常用的边界条件有以下三种：边界条件：边界条件：（1）给定两端点处的导数值导数值,（3）如果f(x)是以b-a为周期的周期函数，则S(x)也应是具有同样的周期的周期函数，在端点处需要满足:（2）给定两端处的二阶导数值二阶导数值,)0()0(),0()0(bSaSbSaSnnyxSyxS)(,)(00nnyxSyxS )(,)(00固定边界条件自然边界条件周期边界条件复习复习构造方法：构造方法：n 用一阶导数表示的三次样条 n 用二阶导数表示的三次样条 3(),(0,1,)iiS xmin3()(0,1,)iiSxMin根据样条函数特性（光滑性）采用导数或二阶导数待定的方法三转角插值三弯矩插值复习复习101iiihxx1,.,1in11iiiihhh0112112122321112111()2)2)2(1)2niiiiiinnnnnyymmmmmmmmmm三三转角方程（固定边界条件）：转角方程（固定边界条件）：1113(1)iiiiiiiiiyyyyhh复习复习11或者写成下式，通过解方程求出mi，最后带入样条函数即可。12233222121212nn 1.1102233221111(1)2nnnnnnnmymmmmy 三对角型系数矩阵复习复习对于自然边界条件与周期边界条件均可得到类似的三对角型系数矩阵。对于用二阶导数表示的三次样条函数（三弯矩插值）同样可得类似的三对角型系数矩阵（三弯矩方程）。线性方程组的数值解法)0()0(),0()0(bSaSbSaSnnyxSyxS )(,)(003()iiSxM复习复习高次插值出现龙格现象高次插值出现龙格现象Taylor/Lagrange/Newton插值插值Hermite插值插值 分段线性分段线性插值插值 分段分段线性插值在节点线性插值在节点处光滑性差处光滑性差分段分段Hermite插值插值 导数导数值不值不容易确定容易确定三次样条插值（先由三次样条插值（先由函数值函数值确定确定导数值导数值，再由，再由分段分段Hermite插值解决问题插值解决问题）实际观测的数据通常含有误差，此时往往只需要求实际观测的数据通常含有误差，此时往往只需要求出一个能够符合大部分数据分布趋势的近似函数即出一个能够符合大部分数据分布趋势的近似函数即可，并不要求近似函数在所有的观测点上与实际观可，并不要求近似函数在所有的观测点上与实际观测值严格相等，通过插值方法来解决此类问题并不测值严格相等，通过插值方法来解决此类问题并不合适。合适。函数逼近问题函数逼近问题小结小结1.9 曲线拟合的最小二乘法曲线拟合的最小二乘法第第1 1章章 插值方法插值方法15问题引入 在工程实践和科学实验中，我们经常需要对实际问题建立函数关系，即y=f(x)。但是不知道f(x)的真实表达式，只得到一系列离散的观测数据（有误差）：xi x1 x2 xnf(xi)f(x1)f(x2)f(xn)u插值带来较大误差；u想能用一个经验函数（规律性）经验函数（规律性）y=s(x)对真实函数y=f(x)作近似，从而得到在其它离散点处的函数值。16 y=s(x)oxy(x0,y0)x1 x2 xn(x1,y1)x0(xn,yn)(x2,y2)拟合的基本思想求一个经验函数 ，()s x曲线拟合问题也就是函数逼近问题。17火山岩储层声波速度与孔隙度的关系火山岩储层声波速度与孔隙度的关系PPOR(6130V)/10790；一般根据岩石测量参数建立岩石物理模型，钻孔提供孔隙度数据可以与岩石储层声波速度建立关系，如下图所示，其关系为实际拟合问题中，已知点个数通常远大于拟合多项式的次数。岩石物理应用18oy ynx y0 x1 x2 xn y1 y2 x0pn(x)s(x)拟合与插值的区别19拟合/逼近问题与插值问题的比较相同点：两类问题都是用简单函数近似未知函数；不同点：插值问题要求近似函数与未知函数在已知节点函数值相等；拟合/逼近问题只要求近似函数能够符合大部分已知数据的分布趋势即可。根据近似函数的形式不同，可分为：代数插值代数逼近（代数多项式）三角插值三角逼近（三角函数）有理插值有理逼近（有理函数）20这些点的分布接近一条真线，yabx这是一个直线方程，无论怎么选择a,b，直线都不可能同时过全部数据点。怎样选取a,b,才能使直线“最好”地反映数据点的基本趋势？首先要建立好坏标准；1、直线拟合21残差（误差）衡量残差（误差）衡量 ：如果某个直线能很好逼近数据规律，即a,b确定，(1,2,3,4)iiyabx i所以根据方程可以计算近似值，他与实测值yi之差(1,2,3,4)iiiiieyyyabxi称为残差残差。残差的大小可作为衡量近似函数好坏的标准。22准则（判断标准）：准则（判断标准）：（2）使残差的绝对值之和最小，即（1）使残差的最大绝对值最小，即（3）使残差的平方和最小，即miniie maxminiie 2miniie 函数的最佳一致逼近函数最佳平方逼近/最小二乘拟合23数据拟合的最小二乘法问题数据拟合的最小二乘法问题根据给定的数据组（xi,yi)(i=1,2,N),选取近似函数形式，即给定函数类H,求函数 ，使Hx)(2211()minnniiiiiQeyx 该函数称为这组数据的最小二乘函数。通常H取为一些比较简单函数的集合，一般为低次多项式，如1次多项式（a,b待定）：()xa bx 24据函数极值方法，下式成立0;0QQab将Q代入可以得到线性方程组2iiiiiiaNbxyaxbxx y25例例:给定一组实验数据如下i1234xi2468yi1.12.84.97.2求y=a+bx的一次拟合曲线4201620120100.4abab解：N=4,带入公式1.1,1.02ab 1.102.1xy得到：2iiiiiiaNbxyaxbxx y262、多项式拟合根据给定的数据组（xi,yi)(i=1,2,N),求一个m 次多项式（mN)010()mmjmmjjypxaa xa xa x使所有采样点的误差满足：220111()(,)minnniimimiieypxQ a aa Q是一个多元函数，所以求这个多元函数的极值问题。27由多元函数取极值的必要条件，得方程组201011(,)()njmmiijimiiQ a aayaa xa xa x1020nmkjikiiikjQya xxa 求导移项得011mNNkjjkiiikiiaxy x),2,1,0(mj28201211112310121111112201211111NNNNmiimiiiiiiNNNNNmiiimiiiiiiiiNNNNNmmmmmiiimiiiiiiiia Naxaxaxyaxaxaxaxy xaxaxaxaxy x写成方程组：这是最小二乘拟合多项式的系数aj所满足的方程组，称为正则方程组或法方程组。29正则方程组解的存在性(P39)30正则方程组解的存在性1、反证法（P39).2、证明系数行列式不为零（数学归纳法）.3、由函数组1,x,x2xm的的线性无关性也可以证明方程组存在唯一解。31是否满足极小值条件(P39)32 有时根据给定数据图形，其拟合函数y=f(x)表面上不是线性或多项式的形式，但通过变换仍可化为线性模型。例如，bxaexS)(,ln)(lnbxaxS若两边取对数得这样就变成了形如的线性模型.此时，若令 则,ln),(ln)(bBaAxSxS,)(BxAxS非线性拟合模型非线性问题线性化33例例 求数据表的最小二乘法拟合的二次多项式函数xi-1-0.75-0.5-0.25 00.250.50.751yi504025201821355666x=-1:0.25:1;y=50,40,25,20,18,21,35,56,66;p=polyfit(x,y,2);xi=-1:0.01:1;yi=polyval(p,xi);plot(xi,yi,x,y,o);在MATLAB命令窗口实现作业作业：P57P57，习题一，习题一35,3735,371010月月2020号（号（下周四下周四）上课时交）上课时交