MATLAB数据拟合实用教程.ppt

1 用Matlab进行数据拟合 1 多项式曲线拟合 polyfit y0 polyval p x0 p polyfit x y m 其中 x y为已知数据点向量 分别表示横 纵坐标 m为拟合多项式的次数 结果返回m次拟合多项式系数 从高次到低次存放在向量p中 可求得多项式在x0处的值y0 2 例1已知观测数据点如表所示 分别用3次和6次多项式曲线拟合这些数据点 x 0 0 1 1y 0 447 1 978 3 28 6 16 7 08 7 34 7 66 9 56 9 48 9 3 11 2 plot x y k markersize 25 axis 01 3 216 p3 polyfit x y 3 p6 polyfit x y 6 编写Matlab程序如下 3 t 0 0 1 1 2s polyval p3 t s1 polyval p6 t holdonplot t s r linewidth 2 plot t s b linewidth 2 grid x 0 0 1 1y 0 447 1 978 3 28 6 16 7 08 7 34 7 66 9 56 9 48 9 3 11 2 plot x y k markersize 25 axis 01 3 216 p3 polyfit x y 3 p6 polyfit x y 6 4 例2用切削机床进行金属品加工时 为了适当地调整机床 需要测定刀具的磨损速度 在一定的时间测量刀具的厚度 得数据如表所示 5 解 描出散点图 在命令窗口输入 t 0 1 16 y 30 029 128 428 128 027 727 527 227 026 826 526 326 125 725 324 824 0 plot t y 6 解 描出散点图 在命令窗口输入 t 0 1 16 y 30 029 128 428 128 027 727 527 227 026 826 526 326 125 725 324 824 0 plot t y a 0 301229 3804 holdon plot t y1 holdoff a polyfit t y 1 y1 0 3012 t 29 3804 7 例2用切削机床进行金属品加工时 为了适当地调整机床 需要测定刀具的磨损速度 在一定的时间测量刀具的厚度 得数据如表所示 切削时间t h 0 30 0 1 29 1 2 28 4 3 28 1 4 28 0 5 27 7 6 27 5 7 27 2 8 27 0 刀具厚度y cm 切削时间t h 9 26 8 10 26 5 11 26 3 12 26 1 13 25 7 14 25 3 15 24 8 16 24 0 刀具厚度y cm 拟合曲线为 y 0 3012t 29 3804 8 例3一个15 4cm 30 48cm的混凝土柱在加压实验中的应力 应变关系测试点的数据如表所示 1 55 2 47 2 93 3 03 已知应力 应变关系可以用一条指数曲线来描述 即假设 式中 表示应力 单位是N m2 表示应变 2 89 9 已知应力 应变关系可以用一条指数曲线来描述 即假设 式中 表示应力 单位是N m2 表示应变 解选取指数函数作拟合时 在拟合前需作变量代换 化为k1 k2的线性函数 于是 令 即 10 在命令窗口输入 x 500 1 0e 61000 1 0e 61500 1 0e 62000 1 0e 62375 1 0e 6 y 3 103 1 0e 32 465 1 0e 31 953 1 0e 31 517 1 0e 31 219 1 0e 3 z log y a polyfit x z 1 k1 exp 8 3009 w 1 552 472 933 032 89 plot x w y1 exp 8 3009 x exp 494 5209 x plot x w x y1 r 11 已知应力 应变关系可以用一条指数曲线来描述 即假设 式中 表示应力 单位是N m2 表示应变 拟合曲线为 令 则 求得 于是 12 在实际应用中常见的拟合曲线有 直线 多项式 一般n 2 3 不宜过高 双曲线 一支 指数曲线 13 2 非线性曲线拟合 lsqcurvefit 功能 x lsqcurvefit fun x0 xdata ydata x resnorm lsqcurvefit fun x0 xdata ydata 根据给定的数据xdata ydata 对应点的横 纵坐标 按函数文件fun给定的函数 以x0为初值作最小二乘拟合 返回函数fun中的系数向量x和残差的平方和resnorm 14 例4已知观测数据点如表所示 求三个参数a b c的值 使得曲线f x aex bx2 cx3与已知数据点在最小二乘意义上充分接近 首先编写存储拟合函数的函数文件 functionf nihehanshu x xdata f x 1 exp xdata x 2 xdata 2 x 3 xdata 3 保存为文件nihehanshu m 15 例4已知观测数据点如表所示 x y 0 3 1 0 1 3 27 0 2 3 81 0 3 4 5 0 4 5 18 0 5 6 0 6 7 05 0 7 8 56 0 8 9 69 0 9 11 25 1 13 17 求三个参数a b c的值 使得曲线f x aex bx2 cx3与已知数据点在最小二乘意义上充分接近 编写下面的程序调用拟合函数 xdata 0 0 1 1 ydata 3 1 3 27 3 81 4 5 5 18 6 7 05 8 56 9 69 11 25 13 17 x0 0 0 0 x resnorm lsqcurvefit nihehanshu x0 xdata ydata 16 编写下面的程序调用拟合函数 xdata 0 0 1 1 ydata 3 1 3 27 3 81 4 5 5 18 6 7 05 8 56 9 69 11 25 13 17 x0 0 0 0 x resnorm lsqcurvefit nihehanshu x0 xdata ydata 程序运行后显示 x 3 00224 03040 9404 resnorm 0 0912 17 例4已知观测数据点如表所示 x y 0 3 1 0 1 3 27 0 2 3 81 0 3 4 5 0 4 5 18 0 5 6 0 6 7 05 0 7 8 56 0 8 9 69 0 9 11 25 1 13 17 求三个参数a b c的值 使得曲线f x aex bx2 cx3与已知数据点在最小二乘意义上充分接近 说明 最小二乘意义上的最佳拟合函数为 f x 3ex 4 03x2 0 94x3 此时的残差是 0 0912 18 f x 3ex 4 03x2 0 94x3 拟合函数为 19 练习 1 已知观测数据点如表所示 求用三次多项式进行拟合的曲线方程 2 已知观测数据点如表所示 求a b c的值 使得曲线f x aex bsinx clnx与已知数据点在最小二乘意义上充分接近 20 插值问题 g表达式复杂 甚至无表达式 21 1 分段线性插值 实用插值方法 2 三次样条插值 细木条 样条 22 输入 节点x0 y0 插值点x 均为数组 长度自定义 输出 插值y 与x同长度数组 1 分段线性插值 已有程序y interp1 x0 y0 x y interp1 x0 y0 x linear 2 三次样条插值 已有程序y interp1 x0 y0 x spline 或y spline x0 y0 x 用Matlab作插值计算 23 例5对在 1 1 上 用n 20的等距分点进行分段线性插值 绘制f x 及插值函数的图形 解在命令窗口输入 x 1 0 1 1y 1 1 9 x 2 xi 1 0 1 1yi interp1 x y xi plot x y r xi yi 24 例6对在 5 5 上 用n 11个等距分点作分段线性插值和三次样条插值 用m 21个插值点作图 比较结果 解在命令窗口输入 n 11 m 21x 5 10 m 1 5y 1 1 x 2 z 0 xx0 5 10 n 1 5y0 1 1 x0 2 y1 interp1 x0 y0 x y2 interp1 x0 y0 x spline x y y1 y2 plot x z r x y k x y1 b x y2 g gtext Piece linear gtext Spline gtext y 1 1 x 2 25 01 00001 00001 00000 50000 80000 75000 82051 00000 50000 50000 50001 50000 30770 35000 29732 00000 20000 20000 20002 50000 13790 15000 14013 00000 10000 10000 10003 50000 07550 07940 07454 00000 05880 05880 05884 50000 04710 04860 04845 00000 03850 03850 0385 例6对在 5 5 上 用n 11个等距分点作分段线性插值和三次样条插值 用m 21个插值点作图 比较结果 x y y1 y2 26 解在命令窗口输入 例7在一天24h内 从零点开始每间隔2h测得的环境温度为 12 9 9 10 18 24 28 27 25 20 18 15 13 单位 推测在每1s时的温度 并描绘温度曲线 t 0 2 24T 129910182428272520181513 plot t T ti 0 1 3600 24T1i interp1 t T ti plot t T ti T1i r T2i interp1 t T ti spline plot t T ti T1i r ti T2i g 27 例8在飞机的机翼加工时 由于机翼尺寸很大 通常在图纸上只能标出部分关键点的数据 某型号飞机的机翼上缘轮廓线的部分数据如下 x04 749 051938577695114133 y05 238 111 9716 1517 116 3414 6312 166 69 x152171190 y7 033 990 28 例8在飞机的机翼加工时 由于机翼尺寸很大 通常在图纸上只能标出部分关键点的数据 某型号飞机的机翼上缘轮廓线的部分数据如下 x 04 749 051938577695114133152171190 y 05 238 111 9716 1517 116 3414 6312 169 697 033 990 xi 0 0 001 190 yi interp1 x y xi spline plot xi yi 29 例9天文学家在1914年8月份的7次观测中 测得地球与金星之间距离 单位 m 并取其常用对数值与日期的一组历史数据如下所示 试推断何时金星与地球的距离 单位 m 的对数值为9 9352 日期 18202224262830 距离对数 9 96189 95449 94689 93919 93129 92329 9150 解由于对数值9 9352位于24和26两天所对应的对数值之间 所以对上述数据用三次样条插值加细为步长为1的数据 30 解由于对数值9 9352位于24和26两天所对应的对数值之间 所以对上述数据用三次样条插值加细为步长为1的数据 x 18 2 30 y 9 96189 95449 94689 93919 93129 92329 9150 xi 18 1 30 yi interp1 x y xi spline A xi yi A 18 000019 000020 000021 000022 000023 000024 000025 000026 000027 000028 000029 000030 00009 96189 95819 95449 95069 94689 94309 93919 93529 93129 92729 92329 91919 9150 31 练习 1 设在区间 2 2 上用10等分点作为节点 分别用三种插值方法 1 计算并输出在该区间的20等分点的函数值 2 输出这个函数及两个插值函数的图形 3 对输出的数据和图形进行分析 32 1 设在区间 2 2 上用10等分点作为节点 分别用三种插值方法 1 计算并输出在该区间的20等分点的函数值 zi 0 01830 03870 07730 14110 23690 36850 52730 69800 85210 95991 00000 95990 85210 69800 52730 36850 23690 14110 07730 03870 0183 33 1 设在区间 2 2 上用10等分点作为节点 分别用两种插值方法 2 输出这个函数及两个插值函数的图形 34 练习 2 已知某型号飞机的机翼断面下缘轮廓线上的部分数据如表所示 假设需要得到x坐标每改变0 1时的y坐标 分别用两种插值方法对机翼断面下缘轮廓线上的部分数据加细 并作出插值函数的图形 35 例5给药方案 一种新药用于临床之前 必须设计给药方案 在快速静脉注射的给药方式下 所谓给药方案是指 每次注射剂量多大 间隔时间多长 药物进入机体后随血液输送到全身 在这个过程中不断地被吸收 分布 代谢 最终排除体外 药物在血液中的浓度 即单位体积血液中的药物含量 称血药浓度 在最简单的一室模型中 将整个机体看作一个房室 称中心室 室内的血药浓度是均匀的 快速静脉注射后 浓度立即上升 然后逐渐下降 当浓度太低时 达不到预期的治疗效果 血药浓度太高 又可能导致药物中毒或副作用太强 临床上 每种药物有一个最小有效浓度c1和一个最大治疗浓度c2 设计给药方案时 要使血药浓度保持在c1 c2之间 设本题所研究药物的最小有效浓度c1 10 最大治疗浓度c2 25 36 例5给药方案 显然 要设计给药方案 必须知道给药后血药浓度随时间变化的规律 为此 从实验和理论两方面着手 在实验方面 对某人用快速静脉注射方式一次注入该药物300mg后 在一定时刻t 小时 采集血样 测得血药浓度c 如表 血药浓度c t 的测试数据 37 例5给药方案