插值及其误差.doc

插值及其误差xsin xcos xtan x1.5670.999 992 80.003 796 3263.411 251.5680.999 996 10.002 796 3357.611 061.5690.999 998 40.001 796 3556.690 981.5700.999 999 70.000 796 31255.765 59用表中的数据和任一插值公式求：（1）用tan x表格直接计算tan 1.569 5。（2）用sin 1.569 5和cos 1.569 5来计算tan 1.569 5。并讨论这两个结果中误差变化的原因。插值：求过已知有限个数据点的近似函数。1 插值方法 下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite 插值和三次样条插值。1.1 拉格朗日多项式插值1.1.1 插值多项式 用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数在区间上个不同点处的函数值，求一个至多次多项式（1）使其在给定点处与同值，即满足插值条件（2）称为插值多项式，称为插值节点，简称节点，称为插值区间。从几何上看，次多项式插值就是过个点，作一条多项式曲线近似曲线。 次多项式（1）有个待定系数，由插值条件（2）恰好给出个方程（3）记此方程组的系数矩阵为，则是范德蒙特(Vandermonde)行列式。当互不相同时，此行列式值不为零。因此方程组（3）有唯一解。这表明，只要个节点互不相同，满足插值要求（2）的插值多项式（1）是唯一的。 插值多项式与被插函数之间的差称为截断误差，又称为插值余项。当充分光滑时，其中。1.1.2 拉格朗日插值多项式实际上比较方便的作法不是解方程（3）求待定系数，而是先构造一组基函数是次多项式，满足令（4）上式称为次 Lagrange 插值多项式，由方程（3）解的唯一性，个节点的次Lagrange 插值多项式存在唯一。1.1.3 用Matlab作Lagrange插值 Matlab中没有现成的Lagrange插值函数，必须编写一个M文件实现Lagrange插值。 设个节点数据以数组输入，个插值点以数组输入，输出数组为个插值。编写一个名为lagrange.m的M文件：function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;endsin 1.5695=0.9999991749999999cos 1.5695=0.001296300000000056tan 1.5695=819.03428749992741.2 分段线性插值1.2.1 插值多项式的振荡 用Lagrange插值多项式近似，虽然随着节点个数的增加，的次数变大，多数情况下误差会变小。但是增大时，的光滑性变坏，有时会出现很大的振荡。理论上，当，在内并不能保证处处收敛于。Runge给出了一个有名的例子：对于较大的，随着的增大，振荡越来越大，事实上可以证明，仅当时，才有，而在此区间外，是发散的。 高次插值多项式的这些缺陷，促使人们转而寻求简单的低次多项式插值。1.2.2 分段线性插值 简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作，它满足，且在每个小区间上是线性函数。 可以表示为 有良好的收敛性，即对于， 用计算点的插值时，只用到左右的两个节点，计算量与节点个数无关。但越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。 1.2.3 用Matlab实现分段线性插值 用 Matlab 实现分段线性插值不需要编制函数程序，Matlab 中有现成的一维插值函数interp1。y=interp1(x0,y0,x,method)method指定插值的方法，默认为线性插值。其值可为：nearest最近项插值linear线性插值spline逐段3次样条插值cubic保凹凸性3次插值。 所有的插值方法要求x0是单调的。 当 x0 为等距时可以用快速插值法，使用快速插值法的格式为*nearest、*linear、*spline、*cubic。1.3 样条插值 许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续，而且要有连续的曲率，这就导致了样条插值的产生。1.3.1 样条函数的概念 所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。 数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间的一个分划如果函数满足： （1）在每个小区间上是次多项式； （2）在上具有阶连续导数。则称为关于分划的次样条函数，其图形称为次样条曲线。称为样条节点，称为内节点，称为边界点，这类样条函数的全体记做，称为次样条函数空间。 显然，折线是一次样条曲线。 若，则是关于分划的次多项式样条函数。次多项式样条函数的一般形式为其中和均为任意常数，而在实际中最常用的是k =2和3的情况，即为二次样条函数和三次样条函数。二次样条函数：对于上的分划，则（5）其中。三次样条函数：对于上的分划，则（6）其中。利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。例如分段线性插值是一次样条插值。下面我们介绍二次、三次样条插值。1.3.2 二次样条函数插值 首先，我们注意到中含有个特定常数，故应需要个插值条件，因此，二次样条插值问题可分为两类： 问题（1）： 已知插值节点和相应的函数值以及端点（或）处的导数值（或），求使得（7） 事实上，可以证明这两类插值问题都是唯一可解的。 对于问题（1），由条件（7） 引入记号为未知向量，为已知向量。 于是，问题转化为求方程组的解的问题，即可得到二次样条函数的表达式。1.3.3 三次样条函数插值 由于中含有个特定常数，故应需要个插值条件，已知插值节点和相应的函数值，这里提供了个条件，还需要2个边界条件。 常用的三次样条函数的边界条件有3种类型： （1）。由这种边界条件建立的样条插值函数称为的完备三次样条插值函数。 特别地，时，样条曲线在端点处呈水平状态。 如果不知道，我们可以要求与在端点处近似相等。这时以为节点作一个三次Newton插值多项式，以作一个三次 Newton插值多项式，要求由这种边界条件建立的三次样条称为的Lagrange三次样条插值函数。（2）。特别地时，称为自然边界条件。（3），（这里要求）此条件称为周期条件。1.3.4 三次样条插值在Matlab中的实现 在Matlab中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。这个条件强迫第1个和第2个三次多项式的三阶导数相等。对最后一个和倒数第2个三次多项式也做同样地处理。 Matlab中三次样条插值也有现成的函数： y=interp1(x0,y0,x,spline)； y=spline(x0,y0,x)； pp=csape(x0,y0,conds)，y=ppval(pp,x)。其中x0，y0是已知数据点，x是插值点，y是插值点的函数值。 对于三次样条插值，提倡使用函数csape，csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppval。pp=csape(x0,y0)：使用默认的边界条件，即Lagrange边界条件。pp=csape(x0,y0,conds)中的conds指定插值的边界条件，其值可为：complete 边界为一阶导数，即默认的边界条件not-a-knot 非扭结条件periodic 周期条件second 边界为二阶导数，二阶导数的值0, 0。variational 设置边界的二阶导数值为0,0。 对于一些特殊的边界条件，可以通过conds的一个12矩阵来表示，conds元素的取值为1，2。此时，使用命令pp=csape(x0,y0_ext,conds)其中y0_ext=left, y0, right，这里left表示左边界的取值，right表示右边界的取值。 conds(i)=j的含义是给定端点i 的j 阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=2,1表示左边界是二阶导数，右边界是一阶导数，对应的值由left和right给出。2 源程序clcclear allclose allx0=1.567 1.568 1.569 1.570;% y0=0.9999928 0.9999961 0.9999984 0.9999997;y00=0.0037963 0.0027963 0.0017963 0.0007963;y0=263.41125 357.61106 556.69098 1255.76559;x=1.5695;digits(16);y1=vpa(lagrange(x0,y0,x) %调用前面编写的Lagrange插值函数y2=vpa(interp1(x0,y0,x) %分段线性插值y3=vpa(interp1(x0,y0,x,spline) %边界为一阶导数的三次样条插值pp1=csape(x0,y0); y4=vpa(ppval(pp1,x) %边界为一阶导数的三次样条插值pp2=csape(x0,y0,second); y5=vpa(ppval(pp2,x) %边界为二阶导数的三次样条插值y1=vpa(lagrange(x0,y00,x) y2=vpa(interp1(x0,y00,x)y3=vpa(interp1(x0,y00,x,spline)pp1=csape(x0,y00); y4=vpa(ppval(pp1,x)pp2=csape(x0,y00,second); y5=vpa(ppval(pp2,x)3 结果3.1 Lagrange插值结果sin1.56950.99999917499999990.15231474667511273738e-7cos1.56950.001296300000000056-0.20389791089232194786e-4sin1.5695/ cos 1.5695771.425730926449659820.20405438626494428004e-4tan1.5695819.03428749992740.61736687561832218063e-13.2分段线性插值结果sin1.56950.9999990500000000-0.10976863026795959047e-6cos1.56950.001296300000000074-0.20389791075348535345e-4sin1.5695/ cos 1.5695771.425634498143153910.20280435958963498322e-4tan1.5695906.22828499994820.174768666190632915303.3 边界为一阶导数的三次样条插值sin1.56950.99999917499999990.15231474667511273738e-7cos1.56950.001296300000000055-0.20389791089901286439e-4sin1.5695/ cos 1.5695771.425730926450228250.20405438627231302622e-4tan1.5695819.03428749992720.61736687561832072346e-13.4 边界为一阶导数的三次样条插值（使用csape）sin1.56950.99999917499999990.15231474667511273738e-7cos1.56950.001296300000000056-0.20389791089232194786e-4sin1.5695/ cos 1.5695771.425730926449659820.20405438626494428004e-4tan1.5695819.03428749992740.61736687561832218063e-13.5 边界为二阶导数的三次样条插值sin1.56950.9999991250000000-0.34768567262268111256e-7cos1.56950.001296300000000083-0.20389791068323067898e-4sin1.5695/ cos 1.5695771.425692355115302230.20355437544243442131e-4tan1.5695858.85081874991940.113351952813560169503.6 MATLAB计算的近似值（作为准确值）sin1.56950.99999915976853804cos1.56950.0012963264318251843tan1.5695771.40998996724351855