数学建模讲座(五)插值和拟合.ppt

MATLAB使用之二插值问题与拟合问题 应用数学学院高崇山 插值与拟合的关系 插值和拟合都是函数逼近或者数值逼近的重要组成部分 他们的共同点都是通过已知一些离散点集M上的约束 求取一个定义在连续集合S M包含于S 的未知连续函数 从而达到获取整体规律的目的 即通过 窥几斑 来达到 知全豹 简单的讲 所谓拟合是指已知某函数的若干离散函数值 f1 f2 fn 通过调整该函数中若干待定系数f 1 2 n 使得该函数与已知点集的差别 最小二乘意义 最小 如果待定函数是线性 就叫线性拟合或者线性回归 主要在统计中 否则叫作非线性拟合或者非线性回归 表达式也可以是分段函数 这种情况下叫作样条拟合 而插值是指已知某函数的在若干离散点上的函数值或者导数信息 通过求解该函数中待定形式的插值函数以及待定系数 使得该函数在给定离散点上满足约束 插值函数又叫作基函数 如果该基函数定义在整个定义域上 叫作全域基 否则叫作分域基 如果约束条件中只有函数值的约束 叫作Lagrange插值 否则叫作Hermite插值 从几何意义上将 拟合是给定了空间中的一些点 找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点 而插值是找到一个 或几个分片光滑的 连续曲面来穿过这些点 插值的主要内容 1 各种插值方法1 1Lagrange插值法1 2分段插值法1 3三次样条插值法1 4二维插值2 插值的Matlab实现2 1一维插值2 2二维插值3 建模实例 水塔流量的估计 1 1Lagrange插值 已知y f x 该函数未知 在互异的n 1个点x0 x1 x2 xn处的函数值y0 y1 y2 yn 则构造一个过n 1个点 xk yk k 0 1 2 n的次数不超过n的多项式y Ln x 称为插值多项式 使其满足Ln xk yk 称为插值条件 然后用y Ln x 作为准确函数y f x 的近似值 此方法称为插值法 Theorem 满足插值条件的次数不超过n的多项式是唯一存在的 Lagrange插值多项式的构造 显然 Ln x 就是满足插值条件的n次多项式 上式称为基函数 Lagrange插值程序 functiony lagr x0 y0 x lagrange插值法的程序 x0 y0 表示已知的n个节点 x表示m个插值点y表示对应于x的m个插值n length x0 m length x fori 1 mz x i s 0 0 fork 1 np 1 0 forj 1 nifj kp p z x0 j x0 k x0 j endends p y0 k s endy i s end Lagrange插值法的缺点 多数情况下 Lagrange插值法效果是不错的 但随着节点数n的增大 Lagrange多项式的次数也会升高 可能造成插值函数的收敛性和稳定性变差 如龙格 Runge 现象 在 1 1 上用n 1个等距节点作插值多项式Ln x 使得它在节点处的值与函数y 1 1 25x2 在对应节点的值相等 当n增大时 插值多项式在区间的中间部分趋于y x 但对于满足条件0 728 x 1的x Ln x 并不趋于y x 在对应点的值 产生了Runge现象 Runge现象的程序 1 clc clf clearall m 21 x 1 1 m 1 1 y 1 1 25 x 2 z 0 x n 3 x0 1 1 n 1 1 y0 1 1 25 x0 2 y1 lagr x0 y0 x subplot 2 2 1 plot x z r x y m holdon 原曲线plot x y1 b gtext L4 x FontSize 12 pause Lagrange曲线n 5 x0 1 1 n 1 1 y0 1 1 25 x0 2 y1 lagr x0 y0 x subplot 2 2 2 plot x z r x y m holdon 原曲线plot x y1 b gtext L8 x FontSize 12 pause Lagrange曲线 Runge现象的程序 2 n 7 x0 1 1 n 1 1 y0 1 1 25 x0 2 y1 lagr x0 y0 x subplot 2 2 3 plot x z r x y m holdon 原曲线plot x y1 b gtext L12 x FontSize 12 pause Lagrange曲线n 9 x0 1 1 n 1 1 y0 1 1 25 x0 2 y1 lagr x0 y0 x subplot 2 2 4 plot x z r x y m holdon 原曲线plot x y1 b gtext L16 x FontSize 12 Lagrange曲线 1 2分段线性插值 可以证明 In x f x 1 3三次样条 设在区间 a b 上 已给n 1个互不相同的节点a x0 x1 xn b而函数y f x 在这些节点的值f xi yi i 0 1 n 如果分段函数S x 满足下列条件 就称S x 为f x 在点x0 x1 xn的三次样条插值函数 1 S x 在子区间 xi xi 1 的表达式Si x 都是次数为3的多项式 2 S xi yi 3 S x 在区间 a b 上有连续的二阶导数 三次样条 即Si x aix3 bix2 cix dii 0 1 nxi 1 x xi 4n个变量 需要4n个方程S xi yii 0 1 n n 1个方程 Si xi Si 1 xi i 1 n 1在xi连续 n 1个方程 Si xi Si 1 xi i 1 n 1在xi连续 n 1个方程 Si xi Si 1 xi i 1 n 1在xi连续 n 1个方程 再加两个条件S x0 S xn 0自然边界条件 2个方程 可以证明 满足上述4n个线性方程组有唯一解 1 4二维插值 1 4 1网格节点插值法已知m n个节点 xi yj zij i 1 2 m j 1 2 n 一般设a x1 x2 xm b c y1 y2 yn d 求任意一点 x y xi yj 处的插值z 1 最临近点插值 2 分片线性插值 3 双线性插值1 4 2散乱点插值法在T a b c d 上散乱分布n个点 一般采用反距离加权平均法 2 插值的Matlab实现 2 1一维插值2 2二维插值 2 1一维插值的实现 基本格式 yc interp1 x y cx method x y分别表示已知数据点的横 纵坐标向量 x必须单调 cx为需要插值的横坐标数据 method为插值方法 有 nearest 最临近点插值 linear 线性插值 默认 spline 三次样条插值 cubic 三次插值注 Lagrange插值法需自编程序 见前面 2 2二维插值的实现 2 2 1插值节点为网格节点 即x y向量是单调的 调用格式 zi interp2 x y z xi yi method Method的4种情况 nearest 最临近点插值 linear 线性插值 默认 spline 三次样条插值 cubic 三次插值说明 这里x和y是两个独立的向量 它们必须是单调的 z是矩阵 是由x和y确定的点上的值 z和x y之间的关系是z i f x y i z j f x j y 即 当x变化时 z的第i行与y的第i个元素相关 当y变化时z的第j列与x的第j个元素相关 如果没有对x y赋值 则默认x 1 n y 1 m n和m分别是矩阵z的行数和列数 二维插值的实现 2 2 2插值点为散乱节点格式 cz griddata x y z cx cy method Method有 linear线性插值 默认 bilinear双线性插值cubic三次插值bicubic双三次插值nearest最近邻域插值注意 cy必须为列向量 3 应用实例 3 1数控机床加工零件3 2山区地形地貌图3 3海底曲面图 3 1数控机床加工零件 图1零件的轮廓线 x间隔0 2 表1x间隔0 2的加工坐标x y 图1右半部的数据 加工时需要x每改变0 05时的y值 模型将图1逆时针方向转90度 轮廓线上下对称 只需对上半部计算一个函数在插值点的值 图2逆时针方向转90度的结果 数控机床加工零件程序 按照表1输入原始数据x 0 0 2 5 4 8 0 2 0 y 54 714 313 683 052 52 051 691 41 1810 860 740 640 570 5 0 440 40 360 320 290 260 240 20 150 1 4 1 96 2 37 2 71 3 3 25 3 47 3 67 3 84 4 4 14 4 27 4 39 4 49 4 58 4 66 4 74 4 8 4 85 4 9 4 94 4 96 4 98 4 99 5 逆时针方向转90度 节点 x y 变为 u v v0 x u0 y 按0 05的间隔在u方向产生插值点u 5 0 05 5 在v方向计算分段线性插值v1 interp1 u0 v0 u 在v方向计算三次样条插值v2 spline u0 v0 u 在 x y 坐标系输出结果 v1 v2 u subplot 1 3 1 plot x y axis 05 55 gtext 原轮廓线 FontSize 12 subplot 1 3 2 plot v1 u axis 05 55 gtext 分段线性插值 FontSize 12 subplot 1 3 3 plot v2 u axis 05 55 gtext 三次样条插值 FontSize 12 数控机床加工零件运行结果 3 2山区地形地貌图 已知某处山区地形选点测量坐标数据为 x 00 511 522 533 544 55y 00 511 522 533 544 555 56海拔高度数据为 z 8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287 山区地形地貌图程序 原始地貌图程序 x 0 5 5 y 0 5 6 xx yy meshgrid x y z 8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287 mesh xx yy z 加密后的地貌图x 0 5 5 y 0 5 6 z 8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287 xi linspace 0 5 50 加密横坐标数据到50个yi linspace 0 6 80 加密纵坐标数据到60个 xii yii meshgrid xi yi 生成网格数据zii interp2 x y z xii yii cubic 插值mesh xii yii zii 加密后的地貌图 山区地形地貌图结果 3 3海底曲面图 例 在某海域测得一些点 x y 处的水深z由下表给出 在矩形区域 75 200 50 150 内画出海底曲面的图形 海底曲面图程序 clc clf clearall x 129140103 588185 5195105157 5107 57781162162117 5 y 7 5141 52314722 5137 585 5 6 5 81356 5 66 584 33 5 z 48686889988949 plot3 x y z o holdon 原始数据点 插值cx 75 0 5 200 cy 70 0 5 150 cz griddata x y z cx cy cubic 三次插值meshz cx cy cz 海底曲面图结果 曲线拟合问题 各种拟合方法 1 线性拟合函数2 多项式曲线拟合函数3 稳健回归函数4 向自定义函数拟合5 非线性曲线拟合 1 线性拟合 调用格式 b regress y X b bint r rint stats regress y X b bint r rint stats regress y X alpha 说明 b regress y X 返回X处y的最小二乘拟合值 该函数求解线性模型 y X 是p 1的参数向量 是服从标准正态分布的随机干扰的n 1的向量 y为n 1的向量 X为n p矩阵 bint返回 的95 的置信区间 r中为形状残差 rint中返回每一个残差的95 置信区间 Stats向量包含R2统计量 回归的F值和p值 线性拟合 例1 设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到 即y 10 x 求线性拟合方程系数 程序 x ones 10 1 1 10 y x 10 1 normrnd 0 0 1 10 1 b bint regress y x 0 05 回归方程为 y 9 9213 1 0143x 2 多项式曲线拟合函数 调用格式 p polyfit x y n p s polyfit x y n 说明 x y为数据点 n为多项式阶数 返回p为幂次从高到低的多项式系数向量p 矩阵s用于生成预测值的误差估计 多项式曲线拟合函数 例2 由离散数据x0 1 2 3 4 5 6 7 8 91y 3 511 41 61 9 6 4 81 52拟合出多项式 程序 x 0 1 1 y 3 511 41 61 9 6 4 81 52 n 3 p polyfit x y n xi linspace 0 1 100 z polyval p xi 多项式求值plot x y o xi z k x y b legend 原始数据 3阶曲线 多项式曲线拟合函数 也可由函数给出数据 例3 x 1 20 y x 3 sin x 程序 x 1 20 y x 3 sin x p polyfit x y 6 xi 1inspace 1 20 100 z poyval p xi 多项式求值函数plot x y o xi z k x y b legend 原始数据 6阶曲线 多项式曲线拟合函数 再用10阶多项式拟合程序 x 1 20 y x 3 sin x p polyfit x y 10 xi linspace 1 20 100 z polyval p xi plot x y o xi z k x y b legend 原始数据 10阶多项式 多项式曲线求值函数 调用格式 y polyval p x y DELTA polyval p x s 说明 y polyval p x 为返回对应自变量x在给定系数p的多项式的值 y DELTA polyval p x s 使用polyfit函数的选项输出s得出误差估计YDELTA 它假设polyfit函数数据输入的误差是独立正态的 并且方差为常数 则YDELTA将至少包含50 的预测值 多项式曲线拟合的评价和置信区间函数调用格式 Y DELTA polyconf p x s Y DELTA polyconf p x s alpha 说明 Y DELTA polyconf p x s 使用polyfit函数的选项输出s给出Y的95 置信区间YDELTA 它假设polyfit函数数据输入的误差是独立正态的 并且方差为常数 1 alpha为置信度 3 稳健回归函数 稳健回归是指此回归方法相对于其他回归方法而言 受异常值的影响较小 调用格式 b robustfit x y b stats robustfit x y 说明 b返回系数估计向量 stats返回各种参数估计 稳健回归函数 例5 演示一个异常数据点如何影响最小二乘拟合值与稳健拟合 首先利用函数y 10 2x加上一些随机干扰的项生成数据集 然后改变一个y的值形成异常值 调用不同的拟合函数 通过图形观查影响程度 程序 x 1 10 y 10 2 x randn 10 1 y 10 0 bls regress y ones 10 1 x 线性拟合brob robustfit x y 稳健拟合scatter x y holdonplot x bls 1 bls 2 x plot x brob 1 brob 2 x r 稳健回归函数 分析 稳健拟合 实线 对数据的拟合程度好些 忽略了异常值 最小二乘拟合 点线 则受到异常值的影响 向异常值偏移 4 自定义函数拟合 对于给定的数据 根据经验拟合为带有待定常数的自定义函数 所用函数 nlinfit 调用格式 beta r J nlinfit X y fun betao 说明 beta返回函数 fun 中的待定常数 r表示残差 J表示雅可比矩阵 X y为数据 fun 自定义函数 beta0待定常数初值 向自定义函数拟合 在化工生产中获得的氯气的级分y随生产时间x下降 假定在x 8时 y与x之间有如下形式的非线性模型 现收集了44组数据 利用该数据通过拟合确定非线性模型中的待定常数 xyxyxy80 49160 43280 4180 49180 46280 40100 48180 45300 40100 47200 42300 40100 48200 42300 38100 47200 43320 41120 46200 41320 40120 46220 41340 40120 45220 40360 41120 43240 42360 36140 45240 40380 40140 43240 40380 40140 43260 41400 36160 44260 40420 39160 43260 41 自定义函数拟合 首先定义非线性函数的m文件 model mfunctionyy model beta0 x a beta0 1 b beta0 2 yy a 0 49 a exp b x 8 程序 x 8 008 0010 0010 0010 0010 0012 0012 0012 0014 0014 0014 0016 0016 0016 0018 0018 0020 0020 0020 0020 0022 0022 0