第4讲 插值与拟合

x0 , x1 , x2 , x3 为节点作一个三次 Newton 插值多项式 N a ( x ) ，以 xn , xn 1, xn 2 , xn 3 作一个三 次 Newton 插值多项式 N b ( x ) ，要求 S ' (a ) N ' a (a ), S ' (b) N ' b (b). 由这种边界条件建立的三次样条称为 f ( x ) 的 Lagrange 三次样条插值函数。 （ii） S ' ' (a ) y ' ' 0 , S ' ' (b) y ' ' n 。特别地 y ' ' 0 y ' ' n 0 时，称为自然边界条件。 （iii） S ' (a 0) S ' (b 0), S ' ' (a 0) S ' ' (b 0) ，此条件称为周期条件。
线性插值 立方样条插值
'cubic' 立方插值 所有的插值方法要求 x0 是单调的。 当 x0 为等距时可以用快速插值法， 使用快速插值法的格式为 '*nearest'、 '*linear'、 '*spline'、 '*cubic'。 1.3.2 三次样条插值 对于三次样条插值，提倡使用函数 csape， csape 的返回值是 pp 形式，要求插值点的函 数值，必须调用函数 fnval。 pp=csape(x0,y0)使用默认的边界条件，即 Lagrange 边界条件。 pp=csape(x0,y0,conds,valconds) 中的 conds 指定插值的边界条件，其值可为 'complete' 边界为一阶导数，一阶导数的值在 valconds 参数中给出，若忽略 valconds 参数，则按缺省情况处理。 'not-a-knot' 非扭结条件。 'periodic' 周期条件。 'second' 边界为二阶导数， 二阶导数的值在 valconds 参数中给出， 若忽略 valconds 参数，二阶导数的缺省值为 [0, 0]。 'variational' 设置边界的二阶导数值为[0,0]。 对于一些特殊的边界条件，可以通过 conds 的一个 1 2 矩阵来表示，conds 元素的取值 为 0， 1， 2。 conds(i)=j 的含义是给定端点 i 的 j 阶导数，即 conds 的第一个元素表示左边界的条件， 第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶导数，对 应的值由 valconds 给出。 详细情况请使用帮助 doc csape。 1.4 举例 例 1 已知速度曲线 v(t ) 上的三个数据点如表 1 所示。用三次样条函数求速度函数，并 求位移 S
: a x0 x1 xn 1 xn b
1
如果函数 S ( x ) 满足： （ i）在每个小区间 [ xi , xi 1 ](i 0,1,, n 1) 上 S ( x ) 是 m 次多项式； （ii） S ( x ) 在 [a, b] 上具有 m 1 阶连续导数。 则称 S ( x ) 为关于分划 的 m 次样条函数，其图形为 m 次样条曲线。 显然，折线是一次样条曲线。 1.2.2 三次样条插值 利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。例如分段线性插 值是一次样条插值。 我们只介绍三次样条插值， 即已知函数 y f ( x ) 在区间 [a, b] 上的 n 1 个节点
1.3 Matlab 插值工具箱 1.3.1 一维插值函数 Matlab 中有现成的一维插值函数 interp1，语法为 y=interp1(x0,y0,x,'method') 其中 method 指定插值的方法，默认为线性插值。其值可为 'nearest' 最近项插值
2
'linear' 'spline'
x y1 y2
根据地图的比例我们知道18mm相当于40km， 试由测量数据计算该国国土的近似面积和边 2 界的近似长度，并与国土面积的精确值41288km 比较。 解 该地区的示意图见图 1。
140 120
100
80
60
40
20
0
20
40
60
80
100
120
140
160
图 1 区域边界示意图 若区域的下边界和上边界曲线的方程分别为 y1 y1 ( x) ， y2 y2 ( x) ， a x b ，则 该地区的边界线长为
17.5
34.0
4
40.5
44.5
48.0
56.0
44 45 47 50 50 38 30 30 34 44 59 70 72 93 100 110 110 110 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 36 34 41 45 46 43 37 33 28 117 118 116 118 118 121 124 121 121 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0 32 65 55 54 52 50 66 66 68 121 122 116 83 81 82 86 85 68]; x0=a(1:3:end,:); x0=x0'; x0=x0(:); % 提取点的横坐标 y1=a(2:3:end,:); y1=y1'; y1=y1(:); % 提出下边界的纵坐标 y2=a(3:3:end,:); y2=y2'; y2=y2(:); % 提出上边界的纵坐标 pp1=csape(x0,y1); pp2=csape(x0,y2); % 计算三次样条插值函数 dp1=fnder(pp1); dp2=fnder(pp2); % 求三次样条插值函数的导数 L1=quad(@(x)sqrt(1+fnval(dp1,x).^2)+sqrt(1+fnval(dp2,x).^2),x0(1),x0(end)) L2=L1/18*40 % 换算成边界的实际长度 S1=quad(@(x)fnval(pp2,x)-fnval(pp1,x),x0(1),x0(end)) % 计算地图上面积 S2=S1/18^2*1600 delta=(S2-41288)/41288*100 注：为了熟悉一些Matlab数值积分的命令，下面我们给出根据表2的数据直接进行数值 积分计算边界长度和国土面积的程序 clc, clear a=[7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 44 45 47 50 50 38 30 30 34 44 59 70 72 93 100 110 110 110 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 36 34 41 45 46 43 37 33 28 117 118 116 118 118 121 124 121 121 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0 32 65 55 54 52 50 66 66 68 121 122 116 83 81 82 86 85 68]; x0=a(1:3:end,:); x0=x0'; x0=x0(:); % 提取点的横坐标 y1=a(2:3:end,:); y1=y1'; y1=y1(:); % 提出下边界的纵坐标 y2=a(3:3:end,:); y2=y2'; y2=y2(:); plot(x0,y1,'*-') %画下边界曲线 hold on plot(x0,y2,'.-') %画上边界曲线 L1=trapz(x0,sqrt(1+gradient(y1,x0).^2)) % 计算下边界的长度 L2=trapz(x0,sqrt(1+gradient(y2,x0).^2)) % 计算上边界的长度 L=L1+L2; %计算地图上边界的长度 LL=L/18*40 %计算实际的边界长度 S=trapz(x0,y2-y1); %计算地图上的近似面积 SS=S/18^2*1600 例3 已知平面区域 0 x 1400 ， 0 y 1200 的高程数据见表 3（单位： m） 。求点
n
用 I n ( x ) 计算 x 点的插值时， 只用到 x 左右的两个节点， 计算量与节点个数 n 无关。 但n 越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够 了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。 1.2 样条插值 许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形， 内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续，而且要有 连续的曲率，这就导致了样条插值的产生。 1.2.1 样条函数的概念 所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或 细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线） ，并使连接点 处有连续的曲率。三次样条插值就是由此抽象出来的。 数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间 [a, b] 的一 个分划
v(t )dt 。
0
2
t v(t )
表 1 速度的三个观测值 0 1 1 3
2 4
解 求解的 Matlab 程序如下 clc, clear t0=0:2; v0=[1 3 4]; pp=csape(t0,v0) xishu=pp.coefs %显示三次多项式的系数 s1=quadl(@(t)fnval(pp,t),0,2) %求数值积分 s2=integral(@(t)fnval(pp,t),0,2) % 用第 2 种方法求数值积分 求出的速度三次样条函数为
a x0 x1 xn 1 xn b
上的值 yi f ( xi )(i 0,1,, n) ，求插值函数 S ( x ) ，使得 （i） S ( xi ) yi
(i 0,1,, n) ；
（ 1）
（ii）在每个小区间 [ x j , x j 1 ]( j 0,1,, n 1) 上 S ( x ) 是三次多项式，记为 S j ( x ) ； （iii） S ( x ) 在 [a, b] 上二阶连续可微。 函数 S ( x ) 称为 f ( x ) 的三次样条插值函数。 由条件（ii） ，不妨将 S ( x ) 记为

b
a
1 y1 ' ( x)2 dx 1 y2 ' ( x)2 dx ，
a
b
计算时用数值积分即可。 计算该区域的面积，可以把该区域看成是上、下两个边界为曲边的曲边四边形，则区域 的面积
S ( y2 ( x) y1 ( x))dx .
a
b
计算相应的数值积分就可求出面积。 为了提高计算的精度，可以把上、下边界曲线分别进行三次样条插值，利用三次样条函 数计算相应的弧长和曲边四边形的面积。 利用三次样条插值计算时，得到边界长度的近似值为 1161km ，区域面积的近似值为 2 2 42476km ，与其准确值 41288km 只相差 2.8768%。 计算的Matlab程序如下 clc, clear a=[7.0 10.5 13.0
第 4 讲 插值与拟合
司守奎 烟台市，海军航空工程学院数学教研室 Email：sishoukui@163.com 插值：求过已知有限个数据点的近似函数。 拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它 在这些点上的总偏差最小。 插值和拟合都是要根据一组数据构造一个函数作为近似， 由于近似的要求不同， 二者的 数学方法上是完全不同的。 而面对一个实际问题， 究竟应该用插值还是拟合， 有时容易确定， 有时则并不明显。 1 插值方法 插值中使用较多的是分段线性插值和三次样条插值。 1.1 分段线性插值 简单地说， 将每两个相邻的节点用直线连起来， 如此形成的一条折线就是分段线性插值 函数，记作 I n ( x ) ，它满足 I n ( xi ) yi ，且 I n ( x ) 在每个小区间 [ xi , xi 1 ] 上是线性函数
2 0.5t 2.5t 1, t [0,1], v(t ) 2 0.5(t 1) 1.5(t 1) 3, t [1,2].
位移 S
v(t )dt 5.6667 。
0
2
3
例 2 已知欧洲一个国家的地图，为了算出它的国土面积和边界长度，首先对地图作如 下测量：以由西向东方向为 x 轴，由南向北方向为 y 轴，选择方便的原点，并将从最西边界 点到最东边界点在 x 轴上的区间适当地分为若干段， 在每个分点的 y 方向测出南边界点和北 边界点的 y 坐标 y1 和 y 2 ，这样就得到了表 2 的数据（单位：mm）。 表 2 某国国土地图边界测量值（ 单位：mm） 10.5 13.0 17.5 34.0 40.5 45 47 50 50 38 59 68.5 34 118 118.0 65 122 70 76.5 41 116 123.5 55 116 72 80.5 45 118 136.5 54 83 93 91.0 46 118 142.0 52 81 100 96.0 43 121 146.0 50 82
x y1 y2 x y1 y2
7.0 44 44 61.0 36 117 111.5 32 121
44.5 30 110 101.0 37 124 150.0 66 86
48.0 30 110 104.0 33 121 157.0 66 85
56.0 34 110 106.5 28 121 158.0 68 68
(i 0,1,, n) 。 I n ( x ) 可以表示为
I n ( x ) y i li ( x )Hale Waihona Puke Baidu
i 0
n
x xi 1 x x , x [ xi 1 , xi ] (i 0时舍去） i 1 i x xi 1 li ( x ) ， x [ xi , xi 1 ] （i n时舍去） x x i i 1 0， 其它 I n ( x ) 有良好的收敛性，即对于 x [a, b] 有， lim I n ( x ) f ( x ) 。
S ( x) {S j ( x), x [ x j , x j 1 ], j 0,1,n 1}
S j ( x) a j x 3 b j x 2 c j x d j
其中 a j , b j , c j , d j 为待定系数，共 4n 个。由条件（ iii）
S j ( x j 1 ) S j 1 ( x j 1 ) （ 2） j 0,1,, n 1 S ' j ( x j 1 ) S ' j 1 ( x j 1 ) S ' ' ( x ) S ' ' ( x ) j 1 j 1 j j 1 容易看出， （1 ） 、 （ 2）式共含有 4n 2 个方程，为确定 S ( x ) 的 4n 个待定参数，尚需再给出
2 个条件。 常用的三次样条函数的边界条件有 3 种类型： （i） S ' (a ) y ' 0 , S ' (b) y ' n 。由这种边界条件建立的样条插值函数称为 f ( x ) 的完备 三次样条插值函数。 特别地， y ' 0 y ' n 0 时，样条曲线在端点处呈水平状态。 如果 f ' ( x ) 不知道，我们可以要求 S ' ( x ) 与 f ' ( x ) 在端点处近似相等。这时以