matlab拟合数据应用.ppt

1 实验目的 实验内容 2 掌握用数学软件求解拟合问题 1 直观了解拟合基本内容 1 拟合问题引例及基本理论 4 实验作业 2 用数学软件求解拟合问题 3 应用实例 2 拟合 2 拟合的基本原理 1 拟合问题引例 3 拟合问题引例1 求600C时的电阻R 设R at ba b为待定系数 4 拟合问题引例2 求血药浓度随时间的变化规律c t 作半对数坐标系 semilogy 下的图形 MATLAB aa1 5 曲线拟合问题的提法 已知一组 二维 数据 即平面上n个点 xi yi i 1 n 寻求一个函数 曲线 y f x 使f x 在某种准则下与所有数据点最为接近 即曲线拟合得最好 y f x i为点 xi yi 与曲线y f x 的距离 6 拟合与插值的关系 函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似 由于近似的要求不同 二者的数学方法上是完全不同的 实例 下面数据是某次实验所得 希望得到X和f之间的关系 MATLAB cn 问题 给定一批数据点 需确定满足特定要求的曲线或曲面 解决方案 若不要求曲线 面 通过所有数据点 而是要求它反映对象整体的变化趋势 这就是数据拟合 又称曲线拟合或曲面拟合 若要求所求曲线 面 通过所给所有数据点 就是插值问题 7 最临近插值 线性插值 样条插值与曲线拟合结果 8 曲线拟合问题最常用的解法 线性最小二乘法的基本思路 第一步 先选定一组函数r1 x r2 x rm x m n 令f x a1r1 x a2r2 x amrm x 1 其中a1 a2 am为待定系数 第二步 确定a1 a2 am的准则 最小二乘准则 使n个点 xi yi 与曲线y f x 的距离 i的平方和最小 记 问题归结为 求a1 a2 am使J a1 a2 am 最小 9 线性最小二乘法的求解 预备知识 超定方程组 方程个数大于未知量个数的方程组 超定方程一般是不存在解的矛盾方程组 如果有向量a使得达到最小 则称a为上述超定方程的最小二乘解 10 线性最小二乘法的求解 定理 当RTR可逆时 超定方程组 3 存在最小二乘解 且即为方程组RTRa RTy的解 a RTR 1RTy 所以 曲线拟合的最小二乘法要解决的问题 实际上就是求以下超定方程组的最小二乘解的问题 11 线性最小二乘拟合f x a1r1 x amrm x 中函数 r1 x rm x 的选取 1 通过机理分析建立数学模型来确定f x 2 将数据 xi yi i 1 n作图 通过直观判断确定f x 12 用MATLAB解拟合问题 1 线性最小二乘拟合 2 非线性最小二乘拟合 13 用MATLAB作线性最小二乘拟合 1 作多项式f x a1xm amx am 1拟合 可利用已有程序 a polyfit x y m 2 对超定方程组 3 多项式在x处的值y可用以下命令计算 y polyval a x 14 例对下面一组数据作二次多项式拟合 15 1 输入以下命令 x 0 0 1 1 y 0 4471 9783 286 167 087 347 669 569 489 3011 2 R x 2 x ones 11 1 A R y MATLAB zxec1 解法1 用解超定方程的方法 2 计算结果 9 810820 1293 0 0317 16 1 输入以下命令 x 0 0 1 1 y 0 4471 9783 286 167 087 347 669 569 489 3011 2 A polyfit x y 2 z polyval A x plot x y k x z r 作出数据点和拟合曲线的图形 2 计算结果 9 810820 1293 0 0317 解法2 用多项式拟合的命令 MATLAB zxec2 17 1 lsqcurvefit已知数据点 xdata xdata1 xdata2 xdatan ydata ydata1 ydata2 ydatan 用MATLAB作非线性最小二乘拟合 Matlab的提供了两个求非线性最小二乘拟合的函数 lsqcurvefit和lsqnonlin 两个命令都要先建立M 文件fun m 在其中定义函数f x 但两者定义f x 的方式是不同的 可参考例题 lsqcurvefit用以求含参量x 向量 的向量值函数F x xdata F x xdata1 F x xdatan T中的参变量x 向量 使得 18 输入格式为 1 x lsqcurvefit fun x0 xdata ydata 2 x lsqcurvefit fun x0 xdata ydata options 3 x lsqcurvefit fun x0 xdata ydata options grad 4 x options lsqcurvefit fun x0 xdata ydata 5 x options funval lsqcurvefit fun x0 xdata ydata 6 x options funval Jacob lsqcurvefit fun x0 xdata ydata 说明 x lsqcurvefit fun x0 xdata ydata options 19 lsqnonlin用以求含参量x 向量 的向量值函数f x f1 x f2 x fn x T中的参量x 使得最小 其中fi x f x xdatai ydatai F x xdatai ydatai 2 lsqnonlin 已知数据点 xdata xdata1 xdata2 xdatan ydata ydata1 ydata2 ydatan 20 输入格式为 1 x lsqnonlin fun x0 2 x lsqnonlin fun x0 options 3 x lsqnonlin fun x0 options grad 4 x options lsqnonlin fun x0 5 x options funval lsqnonlin fun x0 说明 x lsqnonlin fun x0 options 21 例2用下面一组数据拟合中的参数a b k 该问题即解最优化问题 22 MATLAB fzxec1 1 编写M 文件curvefun1 mfunctionf curvefun1 x tdata f x 1 x 2 exp 0 02 x 3 tdata 其中x 1 a x 2 b x 3 k 2 输入命令tdata 100 100 1000cdata 1e 03 4 54 4 99 5 35 5 65 5 90 6 10 6 26 6 39 6 50 6 59 x0 0 2 0 05 0 05 x lsqcurvefit curvefun1 x0 tdata cdata f curvefun1 x tdata F x tdata x a b k 解法1 用命令lsqcurvefit 23 3 运算结果为 f 0 00430 00510 00560 00590 00610 00620 00620 00630 00630 0063x 0 0063 0 00340 2542 4 结论 a 0 0063 b 0 0034 k 0 2542 24 MATLAB fzxec2 解法2用命令lsqnonlinf x F x tdata ctada x a b k 1 编写M 文件curvefun2 mfunctionf curvefun2 x tdata 100 100 1000 cdata 1e 03 4 54 4 99 5 35 5 65 5 90 6 10 6 26 6 39 6 50 6 59 f x 1 x 2 exp 0 02 x 3 tdata cdata 2 输入命令 x0 0 2 0 05 0 05 x lsqnonlin curvefun2 x0 f curvefun2 x 函数curvefun2的自变量是x cdata和tdata是已知参数 故应将cdatatdata的值写在curvefun2 m中 25 3 运算结果为f 1 0e 003 0 2322 0 1243 0 2495 0 2413 0 1668 0 07240 02410 11590 20300 2792x 0 0063 0 00340 2542 可以看出 两个命令的计算结果是相同的 4 结论 即拟合得a 0 0063b 0 0034k 0 2542 26 MATLAB解应用问题实例 1 电阻问题 2 给药方案问题 3 水塔流量估计问题 27 MATLAB dianzu1 电阻问题 得到a1 3 3940 a2 702 4918 方法2 直接用 结果相同 MATLAB dianzu2 28 一室模型 将整个机体看作一个房室 称中心室 室内血药浓度是均匀的 快速静脉注射后 浓度立即上升 然后迅速下降 当浓度太低时 达不到预期的治疗效果 当浓度太高 又可能导致药物中毒或副作用太强 临床上 每种药物有一个最小有效浓度c1和一个最大有效浓度c2 设计给药方案时 要使血药浓度保持在c1 c2之间 本题设c1 10 c2 25 ug ml 一种新药用于临床之前 必须设计给药方案 药物进入机体后血液输送到全身 在这个过程中不断地被吸收 分布 代谢 最终排出体外 药物在血液中的浓度 即单位体积血液中的药物含量 称为血药浓度 29 要设计给药方案 必须知道给药后血药浓度随时间变化的规律 从实验和理论两方面着手 给药方案 1 在快速静脉注射的给药方式下 研究血药浓度 单位体积血液中的药物含量 的变化规律 t 问题 2 给定药物的最小有效浓度和最大治疗浓度 设计给药方案 每次注射剂量多大 间隔时间多长 分析 理论 用一室模型研究血药浓度变化规律 实验 对血药浓度数据作拟合 符合负指数变化规律 3 血液容积v t 0注射剂量d 血药浓度立即为d v 2 药物排除速率与血药浓度成正比 比例系数k 0 模型假设 1 机体看作一个房室 室内血药浓度均匀 一室模型 模型建立 在此 d 300mg t及c t 在某些点处的值见前表 需经拟合求出参数k v 用线性最小二乘拟合c t MATLAB lihe1 计算结果 用非线性最小二乘拟合c t 给药方案设计 设每次注射剂量D 间隔时间 血药浓度c t 应c1 c t c2 初次剂量D0应加大 给药方案记为 2 1 计算结果 给药方案 c1 10 c2 25k 0 2347v 15 02 34 故可制定给药方案 即 首次注射375mg 其余每次注射225mg 注射的间隔时间为4小时 35 估计水塔的流量 2 解题思路 3 算法设计与编程 1 问题 36 某居民区有一供居民用水的园柱形水塔 一般可以通过测量其水位来估计水的流量 但面临的困难是 当水塔水位下降到设定的最低水位时 水泵自动启动向水塔供水 到设定的最高水位时停止供水 这段时间无法测量水塔的水位和水泵的供水量 通常水泵每天供水一两次 每次约两小时 水塔是一个高12 2米 直径17 4米的正园柱 按照设计 水塔水位降至约8 2米时 水泵自动启动 水位升到约10 8米时水泵停止工作 表1是某一天的水位测量记录 试估计任何时刻 包括水泵正供水时 从水塔流出的水流量 及一天的总用水量 37 38 流量估计的解题思路 拟合水位 时间函数 确定流量 时间函数 估计一天总用水量 39 拟合水位 时间函数测量记录看 一天有两个供水时段 以下称第1供水时段和第2供水时段 和3个水泵不工作时段 以下称第1时段t 0到t 8 97 第2次时段t 10 95到t 20 84和第3时段t 23以后 对第1 2时段的测量数据直接分别作多项式拟合 得到水位函数 为使拟合曲线比较光滑 多项式次数不要太高 一般在3 6 由于第3时段只有3个测量记录 无法对这一时段的水位作出较好的拟合 40 2 确定流量 时间函数对于第1 2时段只需将水位函数求导数即可 对于两个供水时段的流量 则用供水时段前后 水泵不工作时段 的流量拟合得到 并且将拟合得到的第2供水时段流量外推 将第3时段流量包含在第2供水时段内 41 3 一天总用水量的估计总用水量等于两个水泵不工作时段和两个供水时段用水量之和 它们都可以由流量对时间的积分得到 42 算法设计与编程 1 拟合第1 2时段的水位 并导出流量 2 拟合供水时段的流量 3 估计一天总用水量 4 流量及总用水量的检验 43 1 拟合第1时段的水位 并导出流量设t h为已输入的时刻和水位测量记录 水泵启动的4个时刻不输入 第1时段各时刻的流量可如下得 1 c1 polyfit t 1 10 h 1 10 3 用3次多项式拟合第1时段水位 c1输出3次多项式的系数2 a1 polyder c1 a1输出多项式 系数为c1 导数的系数3 tp1 0 0 1 9 x1 polyval a1 tp1 x1输出多项式 系数为a1 在tp1点的函数值 取负后边为正值 即tp1时刻的流量 MATLAB llgj1 4 流量函数为 44 2 拟合第2时段的水位 并导出流量设t h为已输入的时刻和水位测量记录 水泵启动的4个时刻不输入 第2时段各时刻的流量可如下得 1 c2 polyfit t 10 9 21 h 10 9 21 3 用3次多项式拟合第2时段水位 c2输出3次多项式的系数2 a2 polyder c2 a2输出多项式 系数为c2 导数的系数3 tp2 10 9 0 1 21 x2 polyval a2 tp2 x2输出多项式 系数为a2 在tp2点的函数值 取负后边为正值 即tp2时刻的流量 MATLAB llgj2 4 流量函数为 45 3 拟合供水时段的流量在第1供水时段 t 9 11 之前 即第1时段 和之后 即第2时段 各取几点 其流量已经得到 用它们拟合第1供水时段的流量 为使流量函数在t 9和t 11连续 我们简单地只取4个点 拟合3次多项式 即曲线必过这4个点 实现如下 xx1 polyval a1 89 取第1时段在t 8 9的流量xx2 polyval a2 1112 取第2时段在t 11 12的流量xx12 xx1xx2 c12 polyfit 891112 xx12 3 拟合3次多项式tp12 9 0 1 11 x12 polyval c12 tp12 x12输出第1供水时段各时刻的流量 MATLAB llgj3 拟合的流量函数为 46 在第2供水时段之前取t 20 20 8两点的流水量 在该时刻之后 第3时段 仅有3个水位记录 我们用差分得到流量 然后用这4个数值拟合第2供水时段的流量如下 dt3 diff t 22 24 最后3个时刻的两两之差dh3 diff h 22 24 最后3个水位的两两之差dht3 dh3 dt3 t 22 和t 23 的流量t3 2020 8t 22 t 23 xx3 polyval a2 t3 1 2 dht3 取t3各时刻的流量c3 polyfit t3 xx3 3 拟合3次多项式t3 20 8 0 1 24 x3 polyval c3 tp3 x3输出第2供水时段 外推至t 24 各时刻的流量 MATLAB llgj4 拟合的流量函数为 47 3 一天总用水量的估计第1 2时段和第1 2供水时段流量的积分之和 就是一天总用水量 虽然诸时段的流量已表为多项式函数 积分可以解析地算出 这里仍用数值积分计算如下 y1 0 1 trapz x1 第1时段用水量 仍按高度计 0 1为积分步长y2 0 1 trapz x2 第2时段用水量y12 0 1 trapz x12 第1供水时段用水量y3 0 1 trapz x3 第2供水时段用水量y y1 y2 y12 y3 237 8 0 01 一天总用水量 计算结果 y1 146 2 y2 266 8 y12 47 4 y3 77 3 y 1250 4 MATLAB llgjz 48 4 流量及总用水量的检验计算出的各时刻的流量可用水位记录的数值微分来检验 用水量y1可用第1时段水位测量记录中下降高度968 822 146来检验 类似地 y2用1082 822 260检验 供水时段流量的一种检验方法如下 供水时段的用水量加上水位上升值260是该时段泵入的水量 除以时段长度得到水泵的功率 单位时间泵入的水量 而两个供水时段水泵的功率应大致相等 第1 2时段水泵的功率可计算如下 p1 y12 260 2 第1供水时段水泵的功率 水量仍以高度计 tp4 20 8 0 1 23 xp2 polyval c3 tp4 xp2输出第2供水时段各时刻的流量p2 0 1 trapz xp2 260 2 2 第2供水时段水泵的功率 水量仍以高度计 计算结果 p1 154 5 p2 140 1 MATLAB ll 49 计算结果 流量函数为 50 流量曲线见图 n 3 4 n 5 6 51 练习1用给定的多项式 如y x3 6x2 5x 3 产生一组数据 xi yi i 1 2 n 再在yi上添加随机干扰 可用rand产生 0 1 均匀分布随机数 或用rands产生N 0 1 分布随机数 然后用xi和添加了随机干扰的yi作的3次多项式拟合 与原系数比较 如果作2或4次多项式拟合 结果如何 52 练习2 用电压V 10伏的电池给电容器充电 电容器上t时刻的电压为 其中V0是电容器的初始电压 是充电常数 试由下面一组t V数据确定V0 53 用非线性最小二乘拟合c t 用lsqcurvefit 2 主程序lihe2 m如下cleartdata 0 250 511 523468 cdata 19 2118 1515 3614 1012 899 327 455 243 01 x0 10 0 5 x lsqcurvefit curvefun3 x0 tdata cdata f curvefun3 x tdata x MATLAB lihe2 1 用M 文件curvefun3 m定义函数functionf curvefun3 x tdata d 300f x 1 d exp x 2 tdata x 1 v x 2 k