1、第 二 章 插 值 与 拟 合2.1插值与拟合的基本概念2.1.插 值与 插值 函数已 知由 ()gx(可 能未 知或 非常 复杂 )产 生的 一批 离散 数据 (,),0,1,iixyi n=， 且 n+1个 互异 插值 节点 01 1n naxx xxb=， 在 插值 区间 内寻 找一 个相 对简 单的 函数 ()fx， 使其 满足 下列 插值 条件 ：再 利 用 已 求 得 的 ()fx计 算 任 一 非 插 值 节 点*x的 近 似 值 * *()yfx=， 这 就 是 插 值 。 其 中()fx称 为插 值函 数， ()gx称 为被 插函 数。下 面介 绍几 种常 用的 而且 有现
2、成 MATLB命 令的 插值 方法 的数 学原 理。1.分 段线 性插 值 将 两 个 相 邻 节 点 用 直 线 连 起 来 ,如 此 形 成 的 一 条 折 线 就 是 分 段 线 性 插 值 函 数 ， 记 作()nIx， 它 满足 ()nj jIxy=， 且 ()nIx在 每个 小区 间 1, , 0,1,j jxxj n+=， 上 是线 性函 数 。()nIx可 以表 示为0() ()nn jjjIxylx=1 111 11, , 0() , ,j j jj jjj j jj jxx xxxjxxxxlx xxxjnxx + + = = 舍 去 ,舍 去 ,0 其 他()nIx有 良
3、好 的收 敛性 ，即 对于 ,xab时 ，有 lim ()()nnIxgx=。用 ()nIx计 算 x点 的插 值时 ,只 用到 x左 右的 两个 节点 ,计 算量 与节 点个 数 n无 关 。 但 是 n越 大， 分段 越多 ，插 值误 差越 小。 MATLB中 有现 成的 分段 线性 插值 命令 ：y=interp1(x0,y,x)其 中 x0,y为 节点 数组 （同 长度 ）， x为 插值 点数 组， y为 插值 数组 。2.三 次样 条插 值 三 次样 条函 数记 作 ()( )Sxaxb， 要求 它满 足以 下条 件：（ 1） 在每 个小 区间 1,(1,)i ixxi n =上 是
4、三 次多 项式 ；（ 2） 在 axb上 二阶 导数 连续 ；（ 3） (),(0,1,)i iSxyi n=。由 条件 （ 1） ，不 妨将 ()Sx记 为1()(),(1,)i i iSxSxxxi n=， 3 2()i i i i iSxaxbxcxd=+，其 中 ,iiiiabcd为 待定 系数 ，共 4n个 。由 条件 （ 2） ，111() (),() (), 1,2,1() ().ii i iii i iii i iSxSxSxSxi nSxSx+= = = =上 式 与 条 件 （ 3） 共 有 4n-2个 方 程 ， 为 确 定 ()Sx的 4n个 待 定 参 数 ， 还 需
5、 两 个 条 件 。 在 实际 应用 中通 常有 以下 三种 类型 的端 点条 件作 为附 加条 件。第 一类 给 定两 端点 的一 阶导 数 0(),()nSxSx ；第 二类 给 定两 端点 的二 阶导 数 0(),()nSxSx ， 最常 用的 是所 谓的 自然 边界 条件 ：0()()0nSxSx =；第 三类 对 于周 期函 数， 即两 端点 已经 满足0()()nSxSx=时 ，令 它们 的一 阶导 数及 二阶 导数 分别 相等 ，即 0()()nSxSx =， 0()()nSxSx =， 称为 周期 条件 。这 样， 就构 成了 4n元 线性 方程 组， 可以 证明 有唯 一解
6、，即 ()Sx被 唯一 确定 。对 于三 次样 条插 值， MATLB中 有现 成的 命令y=interp1(x0,y,x,spline)或 者y=spline(x0,y,x)其 中 x0,y为 节点 数组 （ 同 长度 ） ， x为 插值 点数 组 ， y为 插值 数组 ， 端 点为 自然 边界 条件 。Spline命 令还 可以 处理 上述 第一 类端 点条 件 ， 只 需将 输入 数 组 y0改 为 y0=ay0b，其 中 a， b分 别为 0(),()nSxSx 。3.高 维插 值 MATLB还 给 出 了 高 维 插 值 函 数 interpN()， 其 中 N可 以 为 2， 3，
7、 ， 例 如 N=2时 ，为 二维 插值 ，调 用格 式为 Zi=interp2(x,y,z,xi,yi,method)其 中 x,yz为 插 值 节 点 ， zi为 被 插 值 点 （ xi,yi） 处 的 插 值 结 果 。 method表 示 采 用 的 插 值 方法 ： nearst表 示 最 邻 近 插 值 ； linear表 示 线 性 插 值 ； cubic表 示 三 次 插 值 。 所 有 插 值 方 法都 要求 x,y是 单调 的网 格。2.1.2最 小二 乘拟 合已 知 一 批 离 散 的 数 据 (,),0,1,iixyi n=， ix互 不 相 同 ， 寻 求 一 个
8、拟 合 函 数()yfx=， 使 ()ifx与 iy的 误 差 平 方 和 在 最 小 二 乘 意 义 下 最 小 。 在 最 小 二 乘 意 义 下 确 定 的()fx称 为最 小二 乘拟 合函 数。1.一 元最 小二 乘法给 定 平 面 上 的 点 (,),0,1,iixyi n=， 进 行 曲 线 拟 合 有 多 种 方 法 ， 最 小 二 乘 法 是 解 决曲 线拟 合最 常用 的一 种方 法。 最小 二乘 法的 提法 是： 求 ()fx， 使 2 21 1()n ni i ii i fxy= = 达 到最 小。拟 合时 选用 一定 形式 的拟 合函 数 ， 拟 合函 数可 由一 些简
9、 单的 “ 基 函数 ” （ 例 如幂 函数 ，三 角函 数等 等） 0 1(),(),()mxx x来 线性 表示 ：00 1() () () ()mfxcxcxcx=+现 在 要 确 定 系 数 01,mc c， 使 达 到 极 小 。 为 此 ， 将 ()fx的 表 达 式 代 入 中 ， 就 成为01,mc c的 函数 ， 令 对 (1,2,)ici m=的 偏导 数等 于零 ， 于 是得 到 m +1个 方程 组 ，由 此求 出 (1,2,)ici m=。 通 常取 基函 数为 231,mxxx， 这 时拟 合函 数 ()fx为 多项式 函数 。特 别地 ，当 m =1时 ， ()f
10、xabx=+， 称为 一元 线性 拟合 函数 。MATLBA中 提供 了多 项式 拟合 的语 句：a=polyfit(x,y,m)x,y为 要拟 合的 数据 ，是 长度 自定 义的 数组 ， m 为 拟合 多项 式的 次数 。输 出参 数为 拟合多 项式 的系 数。 多项 式为 降幂 形式 ： 11 2 1mm mmyaxax axa +=+输 出系 数次 序 12 1, maaa+=。多 项式 在 x处 的拟 合值 y可 用下 面程 序计 算：y=polyval(,x)2.如 何选 择拟 合函 数 已 知 一 组 数 据 (,),0,1,iixyi n=， 选 择 什 么 样 的 函 数 (
11、)fx呢 ？ 一 是 根 据 机 理 分 析来 确定 函数 形式 ， 二 是根 据散 点图 直观 判断 函数 ()fx的 形式 。 要 比 较两 个模 型哪 个拟 合效果 更佳 ，则 比较 两个 模型 的残 差平 方和 。残 差平 方和 较小 者更 佳。设 iy为 拟合 函数 的值 ， iy为 测量 值， 则残 差 2( )i iieyy=3.曲 面拟 合简 介 实 际 问 题 中 可 能 遇 到 曲 面 拟 合 问 题 ， 可 以 将 一 元 最 小 二 乘 方 法 的 有 关 概 念 和 结 论 推 广到 多元 最小 二乘 方法 。已 知 m 个 自变 量1(,)mxx和 一个 因变 量
12、y的 一组 观测 值12(,),1,2,i i miixxxyi n=， 要确 定函 数 1(,)myfxx=， 使得2121m in(,)n i i mi iiJfxxxy= 一 般地 ， 首 先通 过机 理分 析或 数据 的直 观判 断 ， 去 确定 函数1(,)myfxx=的 结构 ， 假 定函 数中 含有 未知 参数 12,kaa， 然后 ，通 过最 小二 乘原 理具 体确 定参 数 12,kaa。2.1.3温 度预 测问 题在 12小 时内 ， 每 隔 1小 时测 量一 次温 度 。 温 度依 次为 ： 5， 8， 9， 15， 25， 29， 31， 30，2， 25， 27， 2
13、4。 (单 位 : )（ 1） 试 分 别 用 分 段 线 性 插 值 、 三 次 样 条 插 值 方 法 估 计 在 3.2h， 6.5h， 7.1h， 1.7h的 温 度值 ，每 隔 1/0h估 计一 次温 度值 并画 出其 图形 。（ 2） 用 多项 式拟 合， 估计 在 3.2h， 6.5h， 7.1h， 1.7h的 温度 值。hours=1:2;tem ps5,8,9,15,2,9,31,0,2,5,27,4;t=interp1(hours,tem ps,3.2,6.5,7.1,.7) %线 性插 值Tinterp1(hours,tem ps,3.2,6.5,7.1,.7,splin
14、e) 三 次样 条插 值计 算结 果为 t=10.2030.030.9024.90T9.673430.42731.7525.3820每 隔 1/0h估 计一 次温 度值 并画 出其 图形 ：hours=1:2;tem ps5,8,9,15,2,9,31,0,2,5,27,4;h=1:0.1:2;tinterp1(hours,tem ps,h,spline);plot(hours,tem ps,+,h,tours,tem ps,r:)xlabel(时 间 ),ylabel(温 度 )图 2.1温度插值图形用 三次 多项 式拟 合： hours=1:2;tem ps5,8,9,15,2,9,31,
15、0,2,5,27,4;a=polyfit(hours,tem ps,3)tem ps1=polyval(,hours);plot(hours,tem ps,ro,hours,tem ps1,b.)得 到 3 20.00650.32837.12814.4343y x x x= +， 图形 如下 ：计 算在 3.2h， 6.5h， 7.1h， 1.7h的 温度 值为 ： 14.8， 26.4， 27.30， 23.61。2.1.4山区地貌图 在 某山 区 （ 平 面区 域 （ 0， 280） （ 0， 240） 内 ， 单 位 ： m ） 测 得一 些地 点的 高度 （单 位： m ） 如表 1.
16、所 示， 试作 出该 山区 的地 貌图 。 表 1.某 山区 一些 地点 的高 度X 04080120160202402802401430145014701320128012010894020 1450148015015015014301301201601460150150160150160160160120137015012010150160150138080 1270150120101350145012015040 1230139015015014090101060 18013201450142014013070901.用 二 维 三 次 插 值 方 法 ， 可 以 得 到 不 同 方 向
17、 上 的 山 体 的 高 度 （ 插 值 结 果 ） ， 然 后 在 m atlb中 建立 M-file程 序如 下：y=0:40:240;x0:40:280;z=1801320145014201401307090;123013901501501409010106;12701501201013501450120150;1370150120101501601501380;1460150150160150160160160;145014801501501501430130120;14301450147013201280120108940;xi,yi=m eshgrid(0:2:80,:20:40)
18、;zi=interp2(x,y,zxi,yi,cubic);surf(xi,yi,zi);xlabel(x轴 ),ylabel(y轴 ),zlabel(z轴 );axis(028024020);tile(山 区地 貌图 );2.2行驶汽车车距问题美 国的 某些 司机 培训 课程 中有 这样 的规 则： 正常 驾驶 条件 下车 速每 增加 16公 里 /小 时 ，后 面与 前面 一辆 车的 距离 应增 加一 个车 身的 长度 。 又 云 ， 实 现这 个规 则的 一种 简便 办法 是所谓 2秒 准则 。 即 后车 司机 从前 车经 过某 一标 志开 始默 数 2秒 钟后 到达 同一 标志 ， 而
19、 不管 车速如 何。试 判簖 2秒 准则 与上 述规 则是 一样 的么 ？这 个规 则的 合理 性如 何 ， 是 否有 更好 的规 则 。 下 面表 2-1是 测得 的 车 速和 刹车 距离 的一 组数 据表 2-1车 速和 刹车 距离 的一 组数 据车 速 （ km /h） 2040608010120140刹 车距 离 (m )6.517.83.657.183.418.0153.1.问 题分 析 制 订这 样的 规则 是为 了在 后车 急刹 车情 况下 不致 撞上 前车 ，即 要确 定汽 车的 刹车 距离 。刹 车距 离显 然与 车速 有关 ， 先 计算 出汽 车以 16km /h的 车速
20、在 2秒 钟行 驶的 距离 ， 这 个距 离为 ： 16公 里 /小 时 （ 1小 时 /360秒 ） 2秒 =8.9米 ， 而 汽 车 平 均 长 度 4.6m， 所 以 “ 2秒准 则 ” 与 上述 规则 并不 一致 。要 判断 规则 是否 合理 ，则 要对 刹 车距 离进 行分 析。刹 车 距 离 由 反 应 距 离 和 制 动 距 离 两 部 分 组 成 ， 反 应 距 离 指 司 机 看 到 需 要 刹 车 的 情 况 到汽 车制 动器 开始 起作 用汽 车行 驶的 距离 ， 制 动距 离指 制动 器开 始起 作用 到汽 车完 全停 止行 驶的 距离 。反 应距 离由 反应 时间 和
21、车 速决 定， 反应 时间 取决 于司 机个 人状 况（ 灵敏 、机 警等 ）和 制动 系统 的灵 敏性 ，由 于很 难对 反应 时间 进行 区别 ，因 此， 通常 认为 反应 时间 为常 数 ， 而且 在这 段时 间内 车速 不变 。制 动距 离与 制动 作用 力、 车重 、车 速以 及路 面状 况等 因素 有关 。由 能量 守恒 制动 力所 作的 功被 汽车 动能 的改 变所 抵消 。 设 计制 动器 的一 个合 理原 则是 ， 最 大制 动力 大体 上与 汽车 质量 成正 比 ， 使 汽车 的减 速度 基本 上是 常数 。 路 面状 况可 认为 是固 定的 。2.模 型假 设基 于以 上
22、分 析， 作如 下假 设： （ 1） 刹 车距 离 d等 于反 应距 离 1d和 制动 距离 2d之 和。（ 2） 反 应距 离1d与 车速 v成 正比 ，比 例系 数为 反应 时间 1t。（ 3） 刹 车 时 使 用 最 大 制 动 力 F， F作 的 功 等 于 汽 车 动 能 的 改 变 ， 且 F与 车 的 质 量 m 成 正比 。3.模 型的 建立 由 假设 211dtv=由 假 设 3， 在 F作 用 下 行 驶 距 离 2d作 的 功 F2d使 车 速 从 v变 成 0， 动 能 的 变 化 为 2/mv，有 F2d=2/mv， 又 Fm， 按 照牛 顿第 二定 理 Fma=可
23、知 ， 刹 车时 的减 速度 a为 常数 ，于 是 22dkv=其 中 k为 比例 系数 ，实 际上 1/2ka=。 由假 设 1， 刹车 距离 为21dtvkv=+下 面对 参数 1t和 k进 行估 计 。 通 常有 经验 估计 和数 据拟 合两 种方 法 ， 这 里我 们采 用数 据拟 合的 方法 。如 果没 有首 先推 导出 模型 形式 的话 ，也 可直 接用 多项 式拟 合。 由 上述 推导 得出 ，实 际刹 车距 离的 拟合 多项 式为 d=k1v2+k2v。其 M-file程 序如 下：v=20:20:140/3.6;v2=v.2;x=v;v2;d=6.5,17.8,33.6,57
24、.1,83.4,118.0,153.5;a=xddd=*addd=6.5,17.8,33.6,57.1,83.4,118.0,153.5;b=polyfit(v,ddd,2)y=polal(b,v)plot(v,ddd,ro,v,dd,b)运 行后 得到 1 20.0853, 0.6522k k= =， 即 20.6520.853d v v=+如 果不 知道 函数 的形 式可 以作 出实 际距 离的 2次 多项 式拟 合拟 合 ，d=k1v2+k2v+k3其 M-file程 序如 下：v20:140/3.6;d=6.5,17.8,3.6,57.1,83.4,18.0,153.;k=polyfi
25、t(v,d,2)d1polyval(k,v)plot(v,d1,b,v,d,r*);图 2.函数 d=k1v2+k2v拟合图 二次多项式 d=k1v2+k2v+k3拟合图图 2.给 出了 实际 刹车 距离 与拟 合曲 线的 计算 刹车 距离 的比 较， 对应 数据 列在 表 2-中 。4.模 型应 用由 20.6520.853d v v=+， 根 据 假 设 及 模 型 的 建 立 知 道 ， 反 应 时 间 1t=0.652， 如 果设 制动 时间 为 2t， 如 果已 知速 度则 由物 理学 匀变 速运 动的 运动 方程 可得 到相 应速 度的 制动 距离 ，21222dvt at=， 其
26、中 12ka=5.86。再 考察 2秒 准则 的合 理性 ， 根 据车 距和 汽车 车速 ， 即 可计 算出 汽车 行驶 一段 车距 路程 的时 间 ， 而 车距 就取 为拟 合的 刹车 距离 （ 当 然如 果有 最大 刹车 距离 数据 ， 则 利用 最大 刹车 距离计 算会 更好 ） 。 这个 时间 我们 称为 “ 车 距时 间 ” ， 部分 结果 列在 表 2-中 。根 据计 算结 果， 可以 将 “ 车 距时 间 ” 2秒 规则 改为 ：车 速 019（ km /h） ,“ 车 距时 间 ”为 1秒 ； 车 速 2059（ km /h） ,“ 车 距 时 间 ” 为 2秒 ； 车 速 6
27、09（ km /h） ,“ 车 距 时 间 ”为 3秒 ；车 速 10140（ km /h） ,“ 车 距时 间 ” 为 4秒 。表 2-车 速与 刹车 距离 和计 算刹 车距 离及 车距 时间车 速（ km /h） 2040608010120140实 际刹 车距 离 (m ) 6.517.83.657.183.418.0153.计 算刹 车距 离 （ m ） 6.217.834.5656.283.916.52154.37车 距时 间 (s) 1. 1.62.12.53.03.53.72.3国土面积的数值计算1.问 题的 提出 已 知欧 洲一 个国 家的 的地 图 ， 为 了算 出它 的国 土
28、面 积 ， 首 先对 地图 作如 下测 量 ： 以 由西向 东方 向为 x轴 ， 由 南向 北方 向为 y轴 ， 选 择方 便的 原点 ， 并 将从 最西 边界 点到 最东 边界 点在 x轴 上的 区间 适当 的分 为若 干段 ， 在 每个 分点 的 y方 向测 出南 边界 点和 北边 界点 的 y坐 标1y和 2y， 这样 就得 到了 表 2-3的 数据 ( )m单位： 。根 据 地 图 的 比 例 我 们 知 道 18m相 当 于 40k， 试 由 测 量 数 据 计 算 该 国 国 土 的 近 似 面积 ，与 它的 精确 值 241288km比 较。表 2-3某 国国 土地 图边 界测
29、量值 （ 单位 ： m）x7.010.513.017.534.040.54.548.056.0y4 4547505038303034y4 5970729310101010x61.068.576.580.591.096.010.0104.0106.5y363441454643373 28y1718161818121241212x1.518.0123.5136.5142.0146.0150.157.0158.0y32655 5452506 6 68y1212168381828685682.模 型的 假设我 们作 如下 假设 （ 1） 假 设测 量的 地图 和数 据准 确 ， 由 最 西边 界点 与
30、最 东边 界点 分为 上下 两条 连续 的边 界曲线 ，边 界内 的所 有土 地均 为该 国国 土。 （ 2） 假 设我 们 从 最西 边界 点到 最东 边界 点 ， 变 量 ,xab， 划 分 ,ab为 n小 段1,i ixx，并 由此 将国 土分 为 n小 块 ， 设 每一 小块 均为 X型 区域 。 即 作垂 直于 x轴 的直 线穿 过该 区域 ，直 线与 边界 曲线 最多 只有 两个 交点 。 3.数 值积 分方 法计 算国 土面 积 根 据测 量数 据我 们利 用 Matlb软 件对 上下 边界 进行 三次 多项 式插 值， 得到 图 2.3。数 值积 分法 的基 本思 想是 将上
31、边界 点与 下边 界点 分别 利用 插值 函数 求出 两条 曲线 ， 则 曲线 所围 面积 即为 国土 面积 （地 图上 的国 土面 积 ） ， 然后 根据 比例 缩放 关系 求出 国土 面积 的近似 解 。 在 求国 土面 积时 ， 利 用求 平面 图形 面积 的数 值积 分方 法 将 该面 积分 成若 干个 小长方 形， 分别 求出 长方 形的 面积 后相 加即 为该 面积 的近 似解 。 设 上边 界函 数为2()fx， 下 边界 函数 为 1()fx， 由 定积 分定 义可 知曲 线所 围区 域面 积为 ：2 11()lim ()()nb i i ia nifxdxf f x= 其 中
32、 1,i i ixx。图 2.3某国家国土区域插值边界%三 次多 项式 插值 及面 积计 算源 程序 ：clearlx=7.01.513.017.534.04.54.548.056.061.068.576.580.591.096.01.014.016.51.518.0123.5136.5142.0146.0150.157.0158.0;y1=454750380343641456437328326554250668;y2459702931010101017181618181212412121212168318286568;newx=7:0.1:58;Llength(newx);newy1=int
33、erp1(x,y1,newx,linear);ney2interp1(x,y2,nex,linear);Area=sum (newy2-newy1)*0.1/82*160fil(nexnex(L-1:-:2),newy1newy2(L-1:-:2),gren)holdonplot(x,y1,2,r*)xlabel(东 西距 离（ 单位 ： m m ） ),ylabel(南 北距 离（ 单位 ： m m ） )tile(国 土面 积计 算 -三 次插 值（ 比例 尺为 9:2000） )采 用三 次多 项式 插值 计算 所的 面积 为 42142km， 与其 准确 值 24128km， 只相 差 2.73%。
