W06 插值 拟合

函数逼近 插值与拟合 Matlab第四讲 第4章插值与拟合 插值方法Lagrange法Hermite法分段三次样条法Matlab插值函数interp1pchipspline 拟合方法多项式拟合多元线性回归非线性拟合Matlab拟合函数polyfitregressnlinfit 函数逼近 一般采用函数表示自然规律中各参数的相互关系 解析式 离散点 x 具有简单的形式 插值和拟合是两种重要的函数逼近方法 4 1函数插值 表4 1列表函数 已知某熔盐在423 473K的密度如表4 1所示 估计450K时的密度 如果可以寻找出与已测得的实验数据相适应的近似解析函数式 则可根据近似解析表达式求出未列点的函数值 函数插值的定义 图4 1n次插值多项式的几何表示 插值就是定义一个在特定点取给定值的函数的过程 当选用不同的插值函数时 相应的得到不同的插值方法 多项式插值法 代数多项式插值法 代数多项式具有形式简单 计算方便 存在各阶导数等良好性质 是最常用的插值函数 代数多项式 插值原则 问题 如何确定多项式的系数 以使多项式成为插值多项式 构造多次插值多项式 代数多项式插值基本原理 如果以上的多项式满足插值原则 就称P x 为F x 的插值多项式 已知 对于函数y F x 给出在n 1个互异插值节点上的函数值 可构造一个n次的多项式 代数多项式插值基本原理 确定 次多项式系数 由于节点互异 行列式不等于0 方程组存在唯一解 由于满足插值原则 则有 以上方程组的系数行列式为范得蒙行列式 多项式插值法 线性插值 n 1 对a1和a0进行不同的变换 可得Lagrange或Newton形式的插值函数 求解方程组得到插值多项式系数的方法很多 以拉格朗日 Lagrange 法和牛顿 Newton 法最为常见 插值多项式 直线的两点式表达式 拉格朗日线性插值 分别称为节点x0和x1的一次插值基函数插值函数为基函数的线性组合 组合系数就是对应节点上的函数值 二次插值多项式 拉格朗日二次插值 对于i 0 可假定L0 x A x x1 x x2 基函数满足 基函数如何确定 A 1 x0 x1 x0 x2 拉格朗日插值法求多项式模型为 拉格朗日n次插值 解 按拉格朗日方法 有 例题 拉格朗日插值法的特点 Lagrange插值多项式的优点在于不要求数据点是等间隔的 其缺点是数据点数不宜过大 通常不超过7个 否则计算工作量大且误差大 计算不稳定 拉格朗日插值多项式次数越高 利用被插函数节点信息越多 理应误差越小 可事实并非如此 龙格 Runge 就给出了一个例子 取等矩节点 作拉格朗日插值多项式 当时 函数及插值多项式的图形如4 2所示 龙格现象 1 1 x 0 5 1 0 1 5 y 0 0 2 0 2 由图可见 在区间 0 2 0 2 上比较接近 但在区间 1 1 两端则误差很大 当n继续增大时 部分区间上插值多项式误差偏大的现象更严重 这种现象称 龙格现象 龙格现象 分段插值法 研究表明 分段低次插值多项式的误差小于高次多项式 分段插值法 为避免拉格朗日插值法的不稳定现象 实际进行插值计算时 通常将插值范围分为若干段 然后在每个分段上使用低次插值 这就是分段插值法 我们称由 4 2 式构成的插值多项式为分段线性Lagrange插值多项式 在每个分段插值区间 xi xi 1 内 4 2 分段线性插值法 分段线性插值法的实现 Matlab中的函数interp1可实现线性插值法运算 调用格式为 yi interp1 x y xi method 含义 输入离散数据x y xi 输出xi对应的插值yi 其中 1 x y为样本点 yi为插值点自变量值xi对应的函数值 2 method共有6种参数可供选择 当省略method时 即默认为linear线性插值 注意 向量x为单调 若y为矩阵 则对y的每一列进行插值xi可以为向量 例题 例3 已知某转子流量计在100 1000mL min流量范围内 刻度值与校正值有如下关系 见表4 3 试用线性插值法计算流量计的刻度值为785时 实际流量为多少 表4 3例3数据 例题 解 Matlab计算程序如下 X 100 200 300 400 500 600 700 800 900 1000 Y 105 3 207 2 308 1 406 9 507 5 605 8 707 4 806 7 908 0 1007 9 Xk 780 Yk interp1 X Y Xk 执行结果 Yk 786 8400 这里 X和Y分别表示样本点的刻度值和校正值 Xk和Yk分别表示插值点的刻度值和校正值 线性插值可以保证函数的连续性 但对导数没有要求 分段三次Hermite插值 如果已知插值点的函数值和一阶导数值 利用它们可以构造一个三次函数 使其在插值点的函数值和一阶导数值与已知值相等 这样的分段三次函数称为分段三次Hermite插值 插值原则 许多最有效的插值技术都是基于分段三次多项式 分段三次Hermite插值 以第i个插值区间为例 实际上就是要求取如下方程 依次解得各个区间的系数 最后得分段函数 保形分段三次Hermite插值 函数pchip piecewisecubicHermiteinterpolatingpolynomial 函数pchip可以保证插值函数的一阶导数连续 是一种保形 shape preserving 插值函数 保证函数值不会过度偏离给定数据 问题是如何确定节点处的斜率 节点两侧异号 则取0两侧同号则取平均值端点处另计 保形分段三次Hermite插值 pchip的调用格式为 yi pchip x y xi 或者用yi interp1 x y xi pchip 含义 输入离散数据x y xi 输出xi对应的插值yi 其中 x y为样本点 yi为插值点自变量值xi对应的函数值 三次样条插值函数 三次样条插值原理 在n个小区间构造S x 共有n个三次多项式 需确定4n个参数 在所有节点上 n 1个方程 在除端点外的节点上 3 n 1 个方程 边界条件 已知S x0 S xn 已知S x0 S xn S x0 S xn 0 自然边界条件已知S x0 S xn S x0 S xn 2个方程 三次样条插值法的实现 Matlab中设有专用函数实现三次样条插值法 可调用interp1和spline两个函数实现 interp1函数调用格式为 yi interp1 x y xi spline spline函数调用格式为 yi spline x y xi 两种函数的插值效果等价 其中参数x y xi yi的意义及要求与线性插值interp1中的完全一致 例5 为了绘制乙醇 水二元溶液的气液平衡曲线 查到24点实验数据如表4 4所示 表中x和y分别表示液相和气相中乙醇摩尔分数 但表中缺少低浓度点的数据 且数据分布不均匀 试用三次样条插值法对以下插值点插值并画出乙醇 水二元溶液的气液平衡曲线 例题 表4 4乙醇 水二元溶液的气液平衡曲线实验值 解 Matlab计算程序如下 x 0 0 018 0 054 0 078 0 111 0 147 0 172 0 203 0 210 0 241 0 297 0 320 0 381 0 391 0 423 0 425 0 497 0 507 0 567 0 685 0 767 0 838 0 882 0 894 y 0 0 163 0 328 0 386 0 452 0 491 0 512 0 524 0 520 0 546 0 573 0 577 0 600 0 608 0 621 0 623 0 653 0 656 0 684 0 743 0 796 0 846 0 884 0 894 xk 0 010 0 018 0 020 0 050 0 078 0 100 0 118 0 200 0 300 0 400 0 425 0 500 0 600 0 700 0 800 0 838 0 866 0 899 yk spline x y xk plot x y xk yk r grid legend 样本点 三次样条插值 这里 x和y分别表示样本点的液相和气相中乙醇摩尔分数 xk和yk分别表示插值点的液相和气相中乙醇摩尔分数 注 用yk interp1 x y xk spline 求插值可得到同样的结果 例题 图4 3三次样条插值乙醇 水二元溶液的气液平衡曲线 例题 4 2最小二乘法曲线拟合 化工实验和工程实践中 可测得许多离散的实验数据和工业数据 通常需要寻找一条连续光滑曲线来近似反映已知数据组间存在的某种关系的一般趋势 所得近似函数可以很好地逼近离散数据 但不要求函数完全通过数据点 以排除实验误差的影响这个函数逼近的过程称为曲线拟合或经验建模 0导入 常用最小二乘法曲线拟合 区别 求解误差上的区别 拟合 考虑观察值的误差 误差不可避免时 以偏差的某种最小为拟合标准 如最小二乘法插值 在观测点处无误差运算目的的区别 拟合 是将数据点用最恰当的曲线描述出来 以反映问题的规律 插值 根据某些已知点求近似函数 再根据近似函数求未知点的函数值 联系 二者都是函数逼近的主要方法 都需要求取函数表达式 插值与拟合的联系 最小二乘法曲线拟合的原理 最小二乘法曲线拟合的原理 如果观测数据存在较大误差 通常采用 近似函数在各实验点的计算结果与实验结果的偏差平方和最小 的原则建立近似函数 最小 若 称此曲线拟合法为最小二乘法曲线拟合 定义 经验建模经验建模又分为两种情况 一是无任何理论依据 但有经验公式可供选择 例如很多物性数据 热容 密度 饱和蒸气压 与温度的关系常表示为 二是没有任何经验可循的情况 只能将实验数据画出图形与已知函数图形进行比较 选择图形接近的函数形式作拟合模型 半经验建模如果建模过程中先由一定的理论依据写出模型结构 再由实验数据估计模型参数 这时建立的模型为半经验模型 例如 描述反应速率常数与温度的关系可用阿仑纽斯方程 即这种情况下 工作要点在于如何确定函数中的各未知系数 最小二乘法建模 在选定关联函数的形式之后 问题变换为如何确定所选关联函数中的待定系数 最小二乘法按计算方法特点又分为线性和非线性最小二乘法 对于一元线性函数 测定了m个自变量值 和m个应变量值 计算出m个应变量值 定义误差 1一元线性最小二乘法的方法概述 线性最小二乘法 线性最小二乘法是常用的曲线拟合方法 线性最小二乘法又分为一元和多元等不同情况 欲使Q最小 按极值的必要条件 要满足 由最小二乘法 设 线性最小二乘法 可推导出 解此方程组 就可求出参数a b 上式称为一元线性最小二乘法的法方程组 线性最小二乘法 多元线性最小二乘法的方法概述 设系统共有n个影响因子 得到m次实验数据 若可用多元线性函数拟合时 形式如下 若k代表第k次实验的数据 则相应的预测值表示为 由最小二乘法设 欲使Q最小 按极值的必要条件 要满足 在应用最小二乘法曲线拟合时 通常遇到更多的是非线性函数 对比线性模型拟合 非线性模型拟合要困难的多 最好设法使模型转化为线性形式 有些非线性模型是不能变换成线性模型的 这时应该用直接非线性最小二乘法进行处理 非线性模型拟合的二个途径 直接采用非线性拟合 通过代换转化为线性关系 非线性最小二乘法 非线性函数的线性化处理 化工中常见的函数双曲线幂函数指数函数负指数函数对数函数S型曲线n次多项式 均可以线性化处理 7 n次多项式 令 代换方程为 非线性函数的线性化处理 最小二乘法曲线拟合的Matlab实现 最小二乘法曲线拟合的Matlab实现 1一元多项式拟合函数polyfit是多项式函数拟合命令 它只能拟合一元多项式 可化为线性的非线性函数拟合 必须变换变量后 然后才可以用此函数 调用格式 p polyfit x y n 其中 x y为输入的拟合数据 通常用数组方式输入 n表示多项式的最高阶数 当n 1时不可省略 输出参数p为拟合多项式的系数 在Matlab中有专用的拟合指令 可以方便地给出求解运算 函数polyval 常常与polyfit 联合使用 计算拟合多项式在指定点的值 其调用格式为 y polyval p x y deltay polyval p xi s 返回xi处的拟合函数值 deltay是50 置信区间的y的误差 某化学反应其反应产物的浓度随时间变化的数据如下 试用最小二乘法关联y与t的关系 例题 由曲线形状特点 可用拟合实验数据 例题 解 1 建立数学模型 首先对实验点绘图 得 2 线性化变换 令 则 3 计算参数 程序如下 t 5 10 15 20 25 30 35 40 45 50 55 y 1 27 2 16 2 86 3 44 3 87 4 15 4 37 4 51 4 58 4 62 4 64 Y log y X 1 t p polyfit X Y 1 例题 即 b 7 4962 lna 1 6516 则 a exp 1 6516 5 2153所以y与t的关系可表示为 执行结果 p 7 49621 6516 多项式拟合效果的判断与检验 1 图形直观判断数据点是否有异常 对异常数据点应检查是否为输入有误 或实验不正常所致拟合曲线是否符合原始数据点趋势数据点是否均匀分布与拟合曲线两侧 多项式拟合效果的判断与检验 2 相关系数相关系数越接近1越好 r 0 9正线性相关 r 0 9负线性相关 r 0线性无关 r 0 09非线性相关 多项式拟合效果的判断与检验 3 置信区间置信区间越小越好 至少不能大于拟合所得参数的数量级 残差的标准偏差 如果试验的随机误差服从正态分布 则实验值落在区间 yi 2 之内的概率为95 多项式拟合效果的判断与检验 4 F检验F检验实际就是方差分析 即分析自变量和因变量相互关系 在给定的显著性水平 下 从F分布表查得F 1 n 2 一般取0 01和0 05 当FF0 01时 x y有十分显著的关系 polyfit拟合效果的判断与检验 1 相关系数 以上一例题为例 计算polyfit拟合的效果 p s polyfit X Y 1 alpha 0 05 y delta polyconf p x s alpha R sqrt sum polyval p X mean Y 2 sum Y mean Y 2 2 置信区间 3 F检验 ybar mean Y yt1 polyval p X n length X SSR1 sum yt1 ybar 2 MSR1 SSR1 1SSE1 sum Y yt1 2 MSE1 SSE1 n 2 ftest1 MSR1 MSE1fvalue1 finv 0 95 1 n 2 多元线性拟合函数 Matlab统计工具箱中使用命令regress 实现多元线性拟合 常用调用格式为 b regress y x 其中 x y代表自变量数据矩阵x和因变量数据向量y b为多元线性拟合的参数结果向量 注 因变量数据向量y和自变量数据矩阵x按以下排列方式输入 这里 n为变量个数 m为实验次数 设某气体反应可表示为 其反应动力学方程可用下列非线性方程表示 式中 V为反应速度 K为反应速度常数 依次为A B C的分压 下表为实验测定的不同分压下的V值 试确定此气体反应动力学方程 例题 解 1 首先将动力学方程的两边取对数 令 则上式变换为 2 编程计算参数 程序如下 PA 8 998 8 199 7 901 7 001 7 001 3 310 6 501 7 997 PB 8 298 7 001 6 203 4 302 3 900 3 401 2 601 2 199 PC 2 699 4 402 5 900 8 199 9 895 9 796 10 903 17 797 V 8 58 6 05 4 73 3 35 2 18 2 11 1 88 1 04 x1 log PA x2 log PB x3 log PC y log V x ones 8 1 x1 x2 x3 将x表示成要求的x形式b regress y x 例题 执行结果 b 1 49140 00970 6454 0 6611则 1 4914 K exp 1 4914 4 4433 0 0097 0 6454 0 6611 所以此气体反应动力学方程可表示为 例题 regress回归效果的判断与检验 以上一例题为例 计算regress拟合的效果 alpha 0 05 b bint r rint stats regress y x alpha b bint为回归系数估计值及他们的置信区间 r rint为残差及其置信区间 stats为检验回归模型的统计量 包括四个量为复相关系数R2 为F统计量的值 为相对于所得F统计量值的概率p 当p 是认为回归有意义 误差方差的估计 直接非线性拟合函数 Matlab求解非线性拟合的函数较多 多采用最优化方法解决 其中nlinfit函数的简单调用格式为 beta nlinfit x y fun beta0 这里 输入数据x y分别为m n矩阵和m维列向量 n为变量个数 m为实验次数 因变量数据向量y和自变量数据矩阵x排列方式同regress函数 fun是事先用m 文件定义的非线性函数 beta0是回归系数的初值 beta是估计出的回归系数 nlinfit函数拟合的检验