数值计算方法教案-插值方法.doc

复习：1.数值计算方法的含义2误差及误差限3.误差与有效数字4.数值计算中应注意的问题第二章 插值方法一插值的含义问题提出：已知函数在n+1个点上的函数值，求任意一点的函数值。说明：函数可能是未知的；也可能是已知的，但它比较复杂，很难计算其函数值。解决方法：构造一个简单函数来替代未知（或复杂）函数，则用作为函数值的近似值。二、泰勒（Taylor）插值1.问题提出：已知复杂函数在点的函数值，求附近另一点的函数值。2.解决方法：构造一个代数多项式函数，使得与在点充分逼近。泰勒多项式为：显然，与在点，具有相同的i阶导数值（i=0,1,n）。3.几何意义为：与都过点；与在点处的切线重合；与在点处具有相同的凹凸性；其几何意义可以由下图描述，显然函数能相对较好地在点逼近。4.误差分析（泰勒余项定理）：，其中在与之间。5.举例：已知函数，求。分析：本题理解为，已知“复杂”函数在=100点的函数值为，求的附近一点+15的函数值。解：（1）构造1次泰勒多项式函数：。其中，则有：故有误差分析：函数在100,115区间绝对值的极大值为，则有：于是近似值10.75有三位有效数字。几何意义：显然，也过点（100,10），且就是函数在点（100,10）处的切线，如下图所示。（2）构造2次泰勒多项式函数：。把，及代入，有。分析误差函数在100,115区间绝对值的极大值为，则有于是近似值10.721875有四位有效数字。运行文件taylor.m：%已知函数f(x)=x(1/2)，求f(115)%一次泰勒插值subplot(1,2,1);f=inline(x(1/2);p1=inline(5+0.05*x);fplot(f,-50,300);hold onfplot(p1,-50,300);plot(115,10.75,*)line(115,115,0,10.75)%二次泰勒插值subplot(1,2,2);p2=inline(10+1/20*(x-100)-1/4000/2*(x-100)2);fplot(f,-30,300);hold onfplot(p2,-30,300);plot(115,10.72,*)line(115,115,0,10.72)可以得到以下图形：6.泰勒插值存在的问题：1.函数必须存在n+1阶导函数，即使存在n+1阶导数，计算的工作量也比较大；2.要求h为个小量，若h较大，则计算的误差就很大。三拉格朗日（Lagrange）插值1.问题提出：已知函数在n+1个点上的函数值，求任意一点的函数值。说明：函数可能是未知的；也可能是已知的，但它比较复杂，很难计算其函数值。2.解决方法：构造一个n次代数多项式函数来替代未知（或复杂）函数，则用作为函数值的近似值。设，构造即是确定n+1个多项式的系数。3.构造的依据：当多项式函数也同时过已知的n+1个点时，我们可以认为多项式函数逼近于原来的函数。根据这个条件，可以写出非齐次线性方程组：其系数矩阵的行列式D为范德萌行列式：故当n+1个点的横坐标各不相同时，方程组系数矩阵的行列式D不等于零，故方程组有唯一解。即有以下结论。结论：当已知的n+1个点的横坐标各不相同时，则总能够构造唯一的n次多项式函数，使也过这n+1个点。4.几何意义5举例：已知函数，求。分析：本题理解为，已知“复杂”函数，当x=81,100,121,144时，其对应的函数值为：y=9,10,11,12，当x=115时，求函数值。解：（1）线性插值：过已知的（100,10）和（121,11）两个点，构造1次多项式函数，于是有则。（2）抛物插值：构造2次多项式函数，使得它过已知的（100,10）、（121,11）和（144,12）三个点。于是有2次拉格朗日插值多项式： 则有10.722755505364206.拉格朗日n次插值多项式公式：其中称为基函数（k=0,1,.,n），每一个基函数都是关于x的n次多项式，其表达式为：拉格朗日公式特点：1把每一点的纵坐标单独组成一项；2.每一项中的分子是关于x的n次多项式，分母是一个常数；3.每一项的分子和分母的形式非常相似，不同的是：分子是，而分母是7.误差分析（拉格朗日余项定理），其中在所界定的范围内。针对以上例题的线性插值，有函数在100,115区间绝对值的极大值为，则有：于是近似值有三位有效数字。针对以上例题的抛物线插值，有函数在100,115区间绝对值的极大值为，则有于是近似值10.72275550536420有四位有效数字。8.拉格朗日插值公式的优点公式有较强的规律性，容易编写程序利用计算机进行数值计算。9. 拉格朗日插值通用程序程序流程图如下：文件lagrange.m如下：%拉格朗日插值close alln=input(已知的坐标点数n=?);x=input(x1,x2,.,xn=?);y=input(y1,y2,.,yn=?);xx=input(插值点=？);syms t %定义t为符号量p=0;for k=1:n l=1; for j=1:k-1 l=l*(t-x(j)/(x(k)-x(j); end for j=k+1:n l=l*(t-x(j)/(x(k)-x(j); end p=p+l*y(k);endp=inline(p); %把符号算式p变为函数形式fplot(p,min(min(x),xx)-1,max(max(x),xx)+1); %画多项式函数hold onp(xx) %显示插值点plot(x,y,o,xx,p(xx),*); %画已知点和插值点在MATLAB命令窗口输入：lagrange然后有以下对话过程和结果，已知的坐标点数n=?6x1,x2,.,xn=?1,3,5,7,9,11y1,y2,.,yn=?-1,20,0,-1,12,3插值点=？8ans = 5.67187500000000有以下图形：10.作业1，已知函数sin(x)过以下数据点：x0.791.01.6sin(x)0.7103530.8414710.999574请用线性插值和抛物插值，计算sin(0.63)的值，并分析误差。四牛顿（Newton）插值复习：（1）问题提出：已知函数在n+1个点的值(x0,y0),(x1,y1),.(xn,yn),求当x=x时，y的值。（2）解决方法：构造n次多项式函数，使它也过已知的n+1个点。（3）拉格朗日公式：，（4）拉格朗日公式的优点：结构规律性强，便于编写程序。（5）拉格朗日插值的缺点：无承袭性（继承性）若算出3点的抛物插值精度不够，再进行4点的3次多项式插值时，必须从头算起，前面算出的3点抛物插值的计算结果不能利用。而泰勒插值却是具有承袭性的，如线性插值的结果不精确，那么再加上一项，就变成了泰勒抛物插值，如：泰勒1次插值：泰勒2次插值：。而牛顿插值就是具有承袭性的插值公式1.差商的概念设n+1个点互不相等，则定义：和两点的一阶差商为：,三点的二阶差商为：,四点的三阶差商为：n+1个点的n阶差商为：差商具有对称性：；2.牛顿插值解决的问题与拉格朗日插值解决的问题相同只是表述 n次多项式的公式不同。3.牛顿插公式的推导根据差商的概念，有：是两点的一阶差商；是三点的二阶差商；把以上各式从后向前逐次代入，可以得到：其中以上的表达式称为牛顿插值公式，可以证明，n次牛顿插值多项式与n次拉格朗日插值多项式完全相同，只是表达形式不同。故，拉格朗日余项定理与牛顿余项定理相同：，其中在所界定的范围内。则有公式：4.牛顿插值差商表xiyi一阶差商二阶差商n阶差商*x0y01x1y1fx0,x1(x-x0)x2y2fx1,x2fx0,x1,x2(x-x0)(x-x1)x3y3fx2,x3fx1,x2,x3(x-x0)(x-x2)xn-1yn-1xnynfxn-1,xnfxn-2,xn-1,xnfx0,xn(x-x0)(x-xn-1)5.举例例1：已知函数f(x)当x=-2,-1,0,1时，其对应函数值为f(x)=13,-8,-1,4。求f(0.5)的值。解：根据已知点，填写以下差商表：xiyi一阶差商二阶差商三阶差商*-2131-1-8-21(x+2)0-1714(x+2)(x+1)145-1-5(x+2)(x+1)x则四点三次牛顿插值多项式为：故，=3.625可以在MATLAB下运行程序newton01.m：p3=inline(13-21*(x+2)+14*(x+2)*(x+1)-5*(x+2)*(x+1)*x);fplot(p3,-2.5,2.5);hold onxi=-2,-1,0,1;yi=13,-8,-1,4;plot(xi,yi, *);plot(0.5,p3(0.5),o);可以得到以下图形：例2：已知函数f(x)当x=-2,-1,0,1,2时，其对应函数值为f(x)=13,-8,-1,4,1。求f(0.5)的值。解：该题目与例1相比，就是多了一个点，所以和例1的差商表相比，只需多一列，多一行： xiyi一阶差商二阶差商三阶差商四阶差商*-2131-1-8-21(x+2)0-1714(x+2)(x+1)145-1-5(x+2)(x+1)x21-3-4-11(x+2)(x+1)x(x-1)而5个点的4次牛顿插值多项式是在的基础上多增加1项：则可以在MATLAB下运行程序newton02.m：p4=inline(13-21*(x+2)+14*(x+2)*(x+1)-5*(x+2)*(x+1)*x+(x+2)*(x+1)*x*(x-1);fplot(p4,-2.5,2.5,r);hold onxi=-2,-1,0,1,2;yi=13,-8,-1,4,1;plot(xi,yi, *);plot(0.5,p4(0.5),o);可以得到以下图形：6.牛顿插值的优点（1）具有承袭性质（2）利用差商表，计算多点插值，比拉格朗日公式计算方便。7.牛顿插值算法的通用程序以下是程序流程图：MATLAB的通用程序newton.m为：%牛顿插值close alln=input(已知的坐标点数n=?);x=input(x1,x2,.,xn=?);y=input(y1,y2,.,yn=?);xx=input(插值点=？);% 计算差商：fx1,x2,fx1,x2,x3,.,fx1,x2,.,xnf=y;for i=1:n-1 % 计算第i阶差商 for k=n:-1:i+1 f(k)=(f(k)-f(k-1)/(x(k)-x(k-i); endendsyms t %定义t为符号量p=f(1);for k=2:n l=1; for j=1:k-1 l=l*(t-x(j); end p=p+l*f(k);endp=inline(p); %把符号算式p变为函数形式fplot(p,min(min(x),xx)-1,max(max(x),xx)+1); %画多项式函数hold onp(xx) %显示插值点plot(x,y,o,xx,p(xx),*); %画已知点和插值点在MATLAB命令窗口输入：newton然后有以下对话过程和结果，已知的坐标点数n=?6x1,x2,.,xn=?1,3,5,7,9,11y1,y2,.,yn=?-1,20,0,-1,12,3插值点=？8ans = 5.67187500000000有以下图形：8.作业（1）过(0,6)，(1,7)，(2,20)，(3,81)，(4,250)五个点做多项式函数p(x)，并求p(-2)的值。（2）给出下列函数表，已知函数f(x)是一个多项式函数，试求其次数及x的最高幂的系数。x012345f(x)-7-452665128（3）请写出下面数列中？的值 2,5,9,15,23,？ 2,8,15,29,50,？,125五 埃尔米特（Hermite）插值1问题提出已知函数在n+1个点上的函数值及一阶导函数值，求任意一点的函数值。2.解决方法：构造一个2n+1次代数多项式函数，使得 即，多项式函数也过这n+1个点，且函数f(x)和在这n+1个点上具有相同的切线。3. 埃尔米特插值公式：当节点横坐标各不相同时，存在唯一的n+1次代数多项式函数：其中，4. 举例例1求满足下列条件的埃尔米特插值多项式。12231-1解：根据埃尔米特插值公式有：把表中值代入，得：例2已知函数满足下列数据表：121/30.20-0.14构造3次埃尔米特插值多项式。解：根据埃尔米特插值公式可以构造为：在Matlab命令窗口输入：f=inline(x/(x3+2);p3=inline(19/150*x3-16/25*x2+9/10*x-4/75);fplot(f,0,3);hold onfplot(p3,0,3,r);plot(1,2,1/3,0.2,*);绘出如下图形例3.求二次多项式满足，。其中为已知常数。解：设，根据已知条件有，于是基函数一定含有因子，基函数一定含有因子，基函数一定含有因子。设，则有解得：，则有：，六 分段插值1龙格(Runge)现象（高次多项式插值的缺陷）针对函数选取6个节点:xi：-5,-3,-1,1,3,5；yi：1/26,0.1,0.5,0.5,0.1,1/26可以构造5次多项式函数若选项11个节点xi：-5,-4,-3,-2,-1,0,1,2,3,4,5；yi：1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26可以构造10次多项式函数利用用拉格朗日插值的通用程序（或牛顿插值的通用程序）可以画出f(x),P5(x)和P10(x)的图形。程序runge.m如下：%用拉格朗日插值公式分析龙格现象close alln=6;x=-5,-3,-1,1,3,5;y=1/26,0.1,0.5,0.5,0.1,1/26;x6=x;y6=y;syms t %定义t为符号量p=0;for k=1:n l=1; for j=1:k-1 l=l*(t-x(j)/(x(k)-x(j); end for j=k+1:n l=l*(t-x(j)/(x(k)-x(j); end p=p+l*y(k);endp5=inline(p); %把符号算式p变为函数形式f=inline(1/(x2+1);fplot(f,-5,5); %画原来的函数hold onfplot(p5,-5,5,g); %画5次多项式函数n=11;x=-5,-4,-3,-2,-1,0,1,2,3,4,5;y=1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26;syms t %定义t为符号量p=0;for k=1:n l=1; for j=1:k-1 l=l*(t-x(j)/(x(k)-x(j); end for j=k+1:n l=l*(t-x(j)/(x(k)-x(j); end p=p+l*y(k);endp10=inline(p); %把符号算式p变为函数形式fplot(p10,-5,5,r); %画10次多项式函数legend(f (x),P_5(x),P_1_0(x)plot(x6,y6,*); %画6个已知节点plot(x,y,o); %画10个已知节点plot(-5,5,0,0,k); %画坐标轴plot(0,0,-0.5,2,k);运行该程序，可以绘制出如下图形：从图中可以看出，随着节点的增加，采用高次多项式插值，可以在某些区域较好的逼近原来的函数（如在-2,2区间）；但在高次多项式的两端出现了激烈震荡的现象，这就是所谓的龙格现象。从该图可以看出，在附近时，与f(x)偏离很远。例如，而。这就说明用高次插值多项式来近似f(x)的效果并不好，因而通常不用高次插值。2分段线性插值当节点较多时，可以采用分段线性插值，公式如下：，以上公式即为两点的线性拉格朗日插值公式。例如针对龙格现在的函数选取11个节点:xi：-5,-4,-3,-2,-1,0,1,2,3,4,5；yi：1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26可以构造10个1次多项式函数，即分段线性函数。在MATLAB命令窗口输入：f=inline(1/(x2+1);fplot(f,-5,5); %画原来的函数hold onx=-5,-4,-3,-2,-1,0,1,2,3,4,5;y=1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26;plot(x,y,o)plot(x,y,r)可以得到以下图形：显然和龙格现象相比，分段线性插值函数比和都能更好的逼近原函数f(x)。3三次样条插值（1）样条函数的概念分段线性插值在节点处没有连续的一阶导函数，其光滑性较差。对于飞机的机翼的型线及船舶型往往要求有二阶光滑度（即在节点处要求二阶导函数连续）。样条函数的概念来源于工程设计的实践。所谓“样条”（spline）是早期工程设计中的一种绘图工具，它是富有弹性的细长条。绘图时，用压铁迫使样条通过指定的型值点，并保证样条的光滑外形。在绕度不大的情况下，样条的曲线即为三次样条函数。（2）几何意义（3）构造三次样条函数的理论分析如上图所示，通过已知的六个点，构造5个三次多项式函数分别是：红色、蓝色、黑色、紫色和绿色5根曲线。为确定一根曲线，就需要确定4个待定系数，所以总共需要4*5=20个待定系数。另外，分析需要的约束条件。每一根函数都要过已知的左右两个点，则有5*2=10个约束条件。此外，每两个相邻曲线在相邻点处要求充分光滑，即在连接点处左右两个函数在该点具有1次和2次的导函数连续，图中有4个“中间点”，故又有4*2=8个约束条件。若在整个图形的两端在加2个约束条件，整个3次样条函数就确定了。如：左右两端点上的1阶导函数已知；左右两端点上的2阶导函数已知，如（称为自然边界条件）；若原来的函数f(x)是以xn-x0为周期的周期函数，则y0=yn，且。（4）用MATLAB函数interp1进行三次样条函数的插值例1.对龙格现象中的函数进行11个点的三次样条插值：x=-5,-4,-3,-2,-1,0,1,2,3,4,5;y=1/26,1/17,0.1,0.2,0.5,1,0.5,0.2,0.1,1/17,1/26;xi=-5:0.01:5;yi=interp1(x,y,xi,spline);plot(xi,yi,x,y,o)hold onf=inline(1/(x2+1);fplot(f,-5,5,r)可以绘出下图：红色为原来的函数f(x)，蓝色为通过曲线f(x)上的11个点而构成的三次样条函数。从图中可以发现三次样条函数很好的地描述了函数f(x)。例2.随机构造15个点，用牛顿法（或拉格朗日法）构造14次代数多项式函数，然后再根据这15个点构造3次样条函数。现在在MATLAB命令窗口输入：x=1:15;y=round(10*randn(1,15);然后运行程序newton，有以下对话过程：已知的坐标点数n=?15x1,x2,.,xn=?xy1,y2,.,yn=?y插值点=？2ans = 0xi=1:0.01:15;yi=interp1(x,y,xi,spline);plot(xi,yi,r)可以绘出如下图形显然图中，可以看到龙格现象，如果，在另一个图中重新画3次样条函数：closeplot(xi,yi,x,y,o)可以得到：七、小结本章学习了1泰勒插值2拉格朗日插值3.牛顿插值4.埃尔米特插值5.龙格现象6.分段线性插值7.分段三次插值(3次样条函数)作业：（1）已知单调连续函数y=f(x)的下列数据xi-1.10.01.22.1yi-2.2-1.l1.02.1用插值法计算，x为多少时，f(x)=0。提示：把xi和yi“颠倒”理解。（2）用插值法计算矩阵的特征多项式。提示：为3次多项式函数，故让分别取0,1,2,3时，求出的函数值，再构造3次多项式函数。