宁波大学-数值方法 (3).ppt

《宁波大学-数值方法 (3).ppt》由会员分享，可在线阅读，更多相关《宁波大学-数值方法 (3).ppt（70页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、第六章第六章插值插值 /*Interpolation*/当精确函数当精确函数y=f(x)非常复杂或未知时，在一系列节非常复杂或未知时，在一系列节点点x0 xn处测得函数值处测得函数值y0=f(x0),yn=f(xn)，由由此构造一个简单易算的近似函数来逼近此构造一个简单易算的近似函数来逼近f(x)即即g(x)f(x)，自然地，希望g(x)通过所有的离散点，即满满足条件足条件g(xi)=f(xi)(i=0,n)。x0 x1x2x3x4xg(x)f(x)问题是：问题是：v是否存在唯一v如何构造v误差估计应用：应用：查三角函数表、对数表、平方根表等机械加工中计算机图形学、图形图像编解码中的应用也很广

2、泛多项式多项式!这里的这里的g(x)称为称为f(x)的的插值函数插值函数。最常用的插值函数是最常用的插值函数是?简单易算简单易算?1拉格朗日多项式拉格朗日多项式/*Lagrange Polynomial*/niyxPiin,.,0,)(=求求n次多项式次多项式使得使得条件：条件：无重合节点，即无重合节点，即n=1已知已知x0,x1;y0,y1，求求使得使得111001)(,)(yxPyxP=可见可见P1(x)是过是过(x0,y0)和和(x1,y1)两点的直线。两点的直线。101xxxx 010 xxxx =y0+y1l0(x)l1(x)=10)(iiiyxl称为称为拉氏基函数拉氏基函数/*La

3、grange Basis*/，满足条件满足条件li(xj)=ij/*Kronecker Delta*/101xxxx 010 xxxx l0(x)=l1(x)=例例1:已知已知,利用插值一次多项式求利用插值一次多项式求的近似值。的近似值。解：设,由 lg10 和 lg20 两个值的线性插值得到lg12,且具有两位有效数字。于是,拉格朗日型一次插值多项式为:则插值基本多项式为：2、二次插值x0 x1x2y0y1y2f(x)p2(x)要求构造一个不超过二次的代数多项式：使满足：?不妨令由条件得于是得二次插值或是抛物线插值函数：若记：你来给我想一个好记忆的你来给我想一个好记忆的Li(x)的统一形式：

4、的统一形式：呵呵，这下我记住了。例例2:已知：已知：xi101520yi=lgxi11.17611.3010利用此三值的二次插值多项式求利用此三值的二次插值多项式求lg12的近似值。的近似值。解：设则故所以利用三个点进行抛物插值得到的利用三个点进行抛物插值得到的lg12的值的值,与精确值相比与精确值相比,具有具有3位有效数字，精度提高了位有效数字，精度提高了 高次插值通常优于高次插值通常优于低次插值低次插值但但绝对不是次数越绝对不是次数越高就越好，嘿嘿高就越好，嘿嘿1LagrangePolynomialn 1希望找到希望找到li(x)，i=0,n使得使得 li(xj)=ij；然后然后令令=ni

5、iinyxlxP0)()(，则显然有，则显然有Pn(xi)=yi。li(x)每个每个li 有有n个根个根x0 xi xn=njj i jiniiixxCxxxxxxCxl00)().().()(=j i jiiiixxCxl)(11)(LagrangePolynomial与与有关，而与有关，而与无关无关节点节点f但绝对不是次数越高就越好，嘿嘿但绝对不是次数越高就越好，嘿嘿例：例：在在 5,5上考察上考察的的Ln(x)。取取-5-4-3-2-1012345-0.500.511.522.5n越大，越大，端点附近抖动端点附近抖动越大，称为越大，称为Runge现象现象Ln(x)f(x)分段分段低次低次

6、插值插值 function inter_y=lag_func(x,y,xi)%其中x,y为已知的数据点及其函数值%xi为待插值数据点n=length(xi);inter_y=0;end for j=1:n if j=i a=a*(xi-x(j)/(x(i)-x(j);a=1;end inter_y=inter_y+a*y(i);for i=1:n end提问：提问：v如果当等插值数据点不只是单个数据点，而是一个数组呢？function inter_y=lag_func(x,y,xi)%其中x,y为已知的数据点及其函数值 xi为待插值数据点n=length(x);inter_y=0;for j=

7、1:n if j=i a=1;end讲解编程，加深对算法的理解，掌握算法编程讲解编程，加深对算法的理解，掌握算法编程 inter_y=inter_y+a*y(i);a=a*(xi-x(j)/(x(i)-x(j);v分层设计v课后思考：如果当等插值数据点不只是课后思考：如果当等插值数据点不只是单个数据点，而是一个数组呢？单个数据点，而是一个数组呢？end for i=1:n end for j=1:nv课后思考：如果当等插值数据点不只是课后思考：如果当等插值数据点不只是单个数据点，而是一个数组呢？单个数据点，而是一个数组呢？2牛顿插值牛顿插值/*Newtons Interpolation*/La

8、grange插值虽然易算，但若要增加一个节点时，插值虽然易算，但若要增加一个节点时，全部基函数全部基函数li(x)都需重新算过。都需重新算过。将将Ln(x)改写成改写成的形式，希望每加一个节点时，的形式，希望每加一个节点时，只附加一项只附加一项上去即可。上去即可。?差商差商(亦称均差亦称均差)/*divided difference*/1阶差商阶差商/*the 1st divided difference of f w.r.t.xiand xj*/2阶差商阶差商2NewtonsInterpolation10111010,.,.,.,+=kkkkkxxxxxfxxxfxxf(k+1)阶差商：阶差

9、商：Warning:myheadisexplodingWhatisthepointofthisformula?差商的值与差商的值与xi的顺序无关！的顺序无关！2NewtonsInterpolation 牛顿插值牛顿插值/*Newtons Interpolation*/12n+11+(x x0)2+(x x0)(x xn 1)n+1Nn(x)Rn(x)ai=f x0,xi2NewtonsInterpolation实际计算过程为实际计算过程为f(x0)f(x1)f(x2)f(xn 1)f(xn)f x0,x1f x1,x2 f xn 1,xnf x0,x1,x2 f xn 2,xn 1,xnf x

10、0,xn f(xn+1)f xn,xn+1 f xn 1,xn,xn+1 f x1,xn+1f x0,xn+1给定数值：x13/202F(x)313/435/3构造出函数f的均差表，并写出3次牛顿插值多项式：3stdivideddifference2stdivideddifference1stdivideddifference f(f(xkxk)xkxk3次牛顿插值多项式：x210131y5616224376例2求证：存在三次多项式 P3(x)满足下面函数表证:由于表中给出了六个点的函数值，根据Newton插值公式，一般情况下可以构造出一个5次插值多项式，这个5次多项式必然满足上面函数表.构造

11、差商表如下:由下表知，由于四阶差商以上均为0，所以这个5次多项式实际上是3次多项式故存在三次多项式 满足是所给的函数表。xy一阶差商二阶差商三阶差商四阶差商25656 1164040 02141313 12072 2 3431207376931520 一般情况下，给定n+1个节点，可构造一个n次插值多项式，若得到低于n次的插值多项式，称为“退化”情况，利用Newton插值，很容易检查出是否为“退化”情况，因为利用差商（差分）表，当某一阶差商（差分）为常数时，则下一阶差商（差商）必定为0，此时必会出现“退化”情况。2NewtonsInterpolation 等距节点公式等距节点公式 /*Form

12、ulae with Equal Spacing*/向前差分向前差分/*forward difference*/iiifff=+1ikikiikffff111)(+=k1k1=向后差分向后差分/*backward difference*/111 =ikikikfffi 1iifff=当当节点节点等距等距分布时分布时:Moregivenonp.113-114.xk fk fk 2 fk 3 fk 4 fk x0 f0 x1 f1 x2 f2 x3 f3 x4 f4 f0 f1 f2 f32 f02 f12 f23 f03 f14 f0前向差分表：2NewtonsInterpolation 差分的重

13、要性质差分的重要性质线性：线性：若若f(x)是是m次多项式，则次多项式，则是是次多项式，次多项式，而而差分值可由函数值算出差分值可由函数值算出：=+=njjknjknfjnf0)1(=+=njnjkjnknfjnf0)1(其中其中/*二项式系数二项式系数*/函数值可由差分值算出函数值可由差分值算出：kjnjknfjnf =+=0kkkhkfxxf!,.,00=knkknnnhkfxxxf!,.,1=kkkhff0)()(=由由Rn表达式表达式前前后后2NewtonsInterpolation牛顿公式牛顿公式牛顿前差公式牛顿前差公式/*Newtons forward-difference for

14、mula*/已知f(xi)(i=0,1,2,n)xi-xi-1=h用差分代替差商当x(x0,x1),作变换x=x0+th,0t0andmrealnumbers.Thenumbersareseparatedbyspacesandnewlines.1LagrangePolynomialOutput(represents a space)represents a space)Foreachai,youaresupposedtoprintthefunctionvalueatthispointinthefollowingformat:fprintf(outfile,f(%6.3f)=%12.8en,a,

15、f);Theoutputsoftwotestcasesmustbeseperatedbyablankline.Sample Input Sample Input 70.50.70.91.11.31.51.71.90.480.640.780.890.961.000.990.9550.741.60.551.21.8510.01.00.01.010.51Sample OutputSample Outputf(0.740)=6.68735821e001f(1.600)=1.00341309e+000f(0.550)=5.26938820e001f(1.200)=9.29135742e001f(1.85

16、0)=9.52078056e001f(0.500)=5.00000000e001 When you start writing the program,you will find how easy it is to calculate the Lagrange polynomial.Oh yeah?What if I find the current interpolation not accurate enough?Then you might want to take more interpolating points into account.Right.Then all the Lag

17、range basis,li(x),will have to be re-calculated.Excellent point!We will come to discuss this problemnext time.4分段低次插值分段低次插值/*piecewise polynomial approximation*/RememberwhatIhavesaid?IncreasingthedegreeofinterpolatingpolynomialwillNOTguaranteeagoodresult,sincehigh-degreepolynomialsareoscillating.例：例

18、：在在 5,5上考察上考察的的Ln(x)。取取-5-4-3-2-1012345-0.500.511.522.5n越大，越大，端点附近抖动端点附近抖动越大，称为越大，称为Runge现象现象Ln(x)f(x)分段分段低次低次插值插值3厄米插值厄米插值/*Hermite Interpolation*/不仅要求函数值重合，而且要求若干阶不仅要求函数值重合，而且要求若干阶导数导数也重合。也重合。即：要求插值函数即：要求插值函数(x)满足满足(xi)=f(xi),(xi)=f(xi),(mi)(xi)=f(mi)(xi).注：注：N个条件可以确定个条件可以确定阶多项式阶多项式(N个系数个系数)。N 1要求

19、在要求在1个节点个节点x0处直到处直到m0阶导数都重合的插阶导数都重合的插值多项式即为值多项式即为Taylor多项式多项式其余其余项为项为一般只考虑一般只考虑f与与f 的值。的值。3HermiteInterpolation例：例：设设x0 x1 x2,已知已知f(x0)、f(x1)、f(x2)和和f(x1),求多项式求多项式P(x)满足满足P(xi)=f(xi)，i=0,1,2，且且P(x1)=f(x1),并估计误并估计误差。差。模仿模仿Lagrange多项式的思想，设多项式的思想，设解：解：首先，首先，P 的阶的阶数数=3+=213)()()()()(=0iiixgx1f xhxfxPh0(

20、x)有根有根x1,x2，且且h0(x1)=0 x1是重根。是重根。)()()(22100 xxxxCxh =又又:h0(x0)=1C0h2(x)h1(x)有根有根x0,x2)()()(201xxxxBAxxh +=由余下条件由余下条件h1(x1)=1和和h1(x1)=0可可解。解。与与h0(x)完全类似完全类似。(x)g1有根有根x0,x1,x2g1)()()(2101xxxxxxCx =g1又又:(x1)=1C1可解。可解。其中其中hi(xj)=ij,hi(x1)=0,(xi)=0,(x1)=1g1g1与与Lagrange分析分析完全类似完全类似3HermiteInterpolation一般

21、地，已知一般地，已知x0,xn 处有处有y0,yn 和和y0,yn，求求H2n+1(x)满足满足H2n+1(xi)=yi，H2n+1(xi)=yi。解：解：设设+=ni)()()(=0iixhxhyixH2n+1 n=0iyi其中其中hi(xj)=ij,hi(xj)=0,(xj)=0,(xj)=ijgigihi(x)有根有根x0,xi,xn且都是且都是2重根重根)()()(2xlBxAxhiiii+=由余下条件由余下条件hi(xi)=1和和hi(xi)=0可解可解Ai和和Bi (x)gi有根有根x0,xn,除了除了xi外都是外都是2重根重根gi)()(iili2(x)xxCx=gi又又:(xi

22、)=1Ci=1 hi)(x)(ili2(x)xx=设设则则这样的这样的Hermite插值唯插值唯一一 仿照Lagrange插值的做法，首先确定多项式插值空间的维数，注意到，我们的条件共有2(n+1)个条件，所以，最高次数为2n+13HermiteInterpolationQuiz:给定给定xi=i+1,i=0,1,2,3,4,5.下面哪个是下面哪个是h2(x)的图像的图像？x0-10.5123456yxy0-10.5123456斜率斜率=1 求求Hermite多项式的基本步骤：多项式的基本步骤：写出相应于条件的写出相应于条件的hi(x)、gi(x)的组合形式；的组合形式；对每一个对每一个hi(

23、x)、gi(x)找出尽可能多的条件给出的根；找出尽可能多的条件给出的根；根据多项式的总阶数和根的个数写出表达式；根据多项式的总阶数和根的个数写出表达式；根据尚未利用的条件解出表达式中的待定系数；根据尚未利用的条件解出表达式中的待定系数；最后完整写出最后完整写出H(x)。4PiecewisePolynomialApproximation 分段线性插值分段线性插值 /*piecewise linear interpolation*/在每个区间在每个区间上，用上，用1阶多项式阶多项式(直线直线)逼近逼近f(x):记记，易证：当，易证：当时，时，一致一致失去了原函数的光滑性。失去了原函数的光滑性。分段

24、分段Hermite插值插值 /*Hermite piecewise polynomials*/给定给定在在上利用两点的上利用两点的y 及及y 构造构造3次次Hermite函数函数导数一般不易得到。导数一般不易得到。Howcanwemakeasmoothinterpolationwithoutaskingtoomuchfromf?Headache4分段低次插值分段低次插值/*piecewise polynomial approximation*/RememberwhatIhavesaid?IncreasingthedegreeofinterpolatingpolynomialwillNOTgua

25、ranteeagoodresult,sincehigh-degreepolynomialsareoscillating.例：例：在在 5,5上考察上考察的的Ln(x)。取取-5-4-3-2-1012345-0.500.511.522.5n越大，越大，端点附近抖动端点附近抖动越大，称为越大，称为Runge现象现象Ln(x)f(x)分段分段低次低次插值插值4PiecewisePolynomialApproximation 分段线性插值分段线性插值 /*piecewise linear interpolation*/在每个区间在每个区间上，用上，用1阶多项式阶多项式(直线直线)逼近逼近f(x):记记

26、，易证：当，易证：当时，时，一致一致失去了原函数的光滑性。失去了原函数的光滑性。分段分段Hermite插值插值 /*Hermite piecewise polynomials*/给定给定在在上利用两点的上利用两点的y 及及y 构造构造3次次Hermite函数函数导数一般不易得到。导数一般不易得到。Howcanwemakeasmoothinterpolationwithoutaskingtoomuchfromf?Headache5样条多项式样条多项式半截多项式：半截多项式：样条多项式样条多项式性质：性质：1.小区间上是k次多项式；2.k-1次连续可导；3.不存在k阶导数；三次样条三次样条/*Cu

27、bic Spline*/设设。三次样条函数三次样条函数,且且在在每每个个上上为为三三次次多多项项式式/*cubic polynomial*/。若若它它同同时时还还满满足足，则则称称为为f 的的三三次次样样条条插插值值函函数数/*cubic spline interpolant*/.注：注：三次样条与分段三次样条与分段Hermite插值的根本区别在于插值的根本区别在于S(x)自自身光滑身光滑，不需要知道，不需要知道f 的导数值（除了在的导数值（除了在2个端点可能需要）个端点可能需要）；而；而Hermite插值依赖于插值依赖于f 在所有插值点的导数值。在所有插值点的导数值。f(x)H(x)S(x)

28、5CubicSpline 6.5.3 6.5.3 构造三次样条插值函数的构造三次样条插值函数的三弯矩法三弯矩法 /*method of bending moment*/在在上，记上，记,for)()(1jjjxxxxSxS =对每个对每个j,此为此为3次多项式次多项式则则Sj”(x)为为次多项式，需次多项式，需个点的值确定之。个点的值确定之。12设设Sj”(xj 1)=Mj 1，Sj”(xj)=Mj对于对于x xj 1,xj可可得到得到积分积分1次，可得次，可得Sj(x)：jjjjjjjAhxxMhxxM+2)(2)(2121Sj(x)=jjjjjjjjBxAhxxMhxxM+6)(6)(31

29、31Sj(x)=利用已知利用已知Sj(xj 1)=yj 1Sj(xj)=yj可解可解积分积分2次，可得次，可得Sj(x)：Sj”(x)=jjjjjjhxxMhxxM11 +5CubicSplinejjjjjjjhMMhyyA611 =jjjjjjjjjjjjhxxhMyhxxhMyBxA12211)6()6(+=+下面解决下面解决Mj：利用利用S 在在xj 的的连续性连续性xj 1,xj:Sj(x)=jjjjjjjjjjjhMMxxfhxxMhxxM6,2)(2)(112121 +1111211216,2)(2)(+jjjjjjjjjjjhMMxxfhxxMhxxMxj,xj+1:Sj+1(x

30、)=jjjjjjjAhxxMhxxM+2)(2)(2121Sj(x)=利用利用Sj(xj)=Sj+1(xj)xj 1,xj:Sj(x)=jjjjjjjjjjjhMMxxfhxxMhxxM6,2)(2)(112121 +1111211216,2)(2)(+jjjjjjjjjjjhMMxxfhxxMhxxMxj,xj+1:Sj+1(x)=5CubicSpline利用利用Sj(xj)=Sj+1(xj)，合并关于合并关于Mj 1、Mj、Mj+1的同类的同类项，项，记记,整理：整理：11jjjjhhh+=l l1jj-=l lm m),(6111jjjjjjjxxfxxfhhg-+-+=211gMMMj

31、jjjjj=+-l lm m j 1n 1有有个未知数，个未知数，有有个方程。个方程。n 1n+1解方程还需解方程还需2个个边界条件边界条件/*boundary conditions*/5CubicSpline第第1类边条件类边条件/*clamped boundary*/：S(a)=y0，S(b)=yna,x1:S1(x)=1011012112106,2)(2)(hMMxxfhaxMhxxM +010110),(62gy0 xxfhMM=+nnnnnngxxfynhMM=+),(6211类似地利用类似地利用xn 1,b 上的上的Sn(x)有有n+1个方程个方程(加了两行加了两行)：211gMM

32、Mjjjjjj=+-l lm m1 j n 1010110),(62gy0 xxfhMM=+nnnnnngxxfynhMM=+),(62115CubicSpline第第2类边条件：类边条件：S”(a)=y0”=M0，S”(b)=yn”=Mn特别地，特别地，M0=Mn=0称为称为自由边界自由边界/*free boundary*/，对应的对应的样条函数称为样条函数称为自然样条自然样条/*Natural Spline*/。有有n-1个方程（减了两行两列）：个方程（减了两行两列）：211gMMMjjjjjj=+-l lm m1 j n 15CubicSpline第第3类边条件类边条件/*periodi

33、c boundary*/：当当f为为周期函数周期函数时，时，yn=y0，S(a+)=S(b)，S”(a+)=S”(b)有有n个方程：个方程：M0=Mn（0行和行和n行相同）行相同）(加了一行加了一行)：yn=y0，S(a+)=S(b)，S”(a+)=S”(b)a=x0,b=xnspline插值的解题步骤：插值的解题步骤：1.2.3.4.5.第第2类边条件：类边条件：S”(a)=y0”=M0，S”(b)=yn”=Mn这时：这时：例题：给定以下数据表，求三次样条插值函数，使其满足边界条件：例题：给定以下数据表，求三次样条插值函数，使其满足边界条件：y”(27.7)=y”(30)=0X27.7282

34、930F(x)4.14.34.13.0第二类边界条件下：这里n=3,节点为0,1，2，3jjjjjjjjjjjjhxxhMyhxxhMyBxA12211)6()6(+=+jjjjjjjjBxAhxxMhxxM+6)(6)(3131Sj(x)=如果如果y(27.7)=3.0y(30)=-4.0,那又如何？那又如何？课堂作业：课堂作业：例题：给定以下数据表，求三次样条插值函数，使其满例题：给定以下数据表，求三次样条插值函数，使其满足边界条件：足边界条件：X27282930F(x)4.14.34.04.2如果如果y”(27)=y”(30)=1,如何？如何？5CubicSpline插值法小结插值法小结 Lagrange:给出给出y0 yn，选基函数选基函数li(x)，其次数为其次数为节点数节点数1。Newton Ln(x)，只是形式不同；节点等距或渐增节点只是形式不同；节点等距或渐增节点时方便处理。时方便处理。Hermite：给出给出yi 及及yi，选选hi(x)及及hi(x)。Spline：分段低次分段低次,自身光滑自身光滑,f 的导数只在边界给出。的导数只在边界给出。HW:p.131#5-2