插值及其误差.doc

插值及其误差
x
sin x
cos x
tan x
1.567
0.999 992 8
0.003 796 3
263.411 25
1.568
0.999 996 1
0.002 796 3
357.611 06
1.569
0.999 998 4
0.001 796 3
556.690 98
1.570
0.999 999 7
0.000 796 3
1255.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;
end
sin 1.5695=0.9999991749999999
cos 1.5695=0.001296300000000056
tan 1.5695=819.0342874999274
1.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的一个1×2矩阵来表示，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 源程序
clc
clear all
close all
x0=[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插值结果
sin 1.5695 0.9999991749999999 0.15231474667511273738e-7
cos 1.5695 0.001296300000000056 -0.20389791089232194786e-4
sin 1.5695/ cos 1.5695 771.42573092644965982 0.20405438626494428004e-4
tan 1.5695 819.0342874999274 0.61736687561832218063e-1
3.2分段线性插值结果
sin 1.5695 0.9999990500000000 -0.10976863026795959047e-6
cos 1.5695 0.001296300000000074 -0.20389791075348535345e-4
sin 1.5695/ cos 1.5695 771.42563449814315391 0.20280435958963498322e-4
tan 1.5695 906.2282849999482 0.17476866619063291530
3.3 边界为一阶导数的三次样条插值
sin 1.5695 0.9999991749999999 0.15231474667511273738e-7
cos 1.5695 0.001296300000000055 -0.20389791089901286439e-4
sin 1.5695/ cos 1.5695 771.42573092645022825 0.20405438627231302622e-4
tan 1.5695 819.0342874999272 0.61736687561832072346e-1
3.4 边界为一阶导数的三次样条插值（使用csape）
sin 1.5695 0.9999991749999999 0.15231474667511273738e-7
cos 1.5695 0.001296300000000056 -0.20389791089232194786e-4
sin 1.5695/ cos 1.5695 771.42573092644965982 0.20405438626494428004e-4
tan 1.5695 819.0342874999274 0.61736687561832218063e-1
3.5 边界为二阶导数的三次样条插值
sin 1.5695 0.9999991250000000 -0.34768567262268111256e-7
cos 1.5695 0.001296300000000083 -0.20389791068323067898e-4
sin 1.5695/ cos 1.5695 771.42569235511530223 0.20355437544243442131e-4
tan 1.5695 858.8508187499194 0.11335195281356016950
3.6 MATLAB计算的近似值（作为准确值）
sin 1.5695 0.99999915976853804
cos 1.5695 0.0012963264318251843
tan 1.5695 771.40998996724351855