曲线曲面的插值与拟合方法次课.pptx

《曲线曲面的插值与拟合方法次课.pptx》由会员分享，可在线阅读，更多相关《曲线曲面的插值与拟合方法次课.pptx（43页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、引 例 1 函 数 查 表 问 题：已 知 标 准 正 态 分 布 函 数 表，求 表 中 没 有 的 值 (2.3 4)=0.9 9 0 3 6 (2.3 5)=0.9 9 0 6 1求 (2.3 4 5 7)(2.3 5-2.3 4 5 7)/(2.3 5-2.3 4)*(2.3 4)+(2.3 4 5 7-2.3 4)/(2.3 5-2.3 4)*(2.3 5)引 例 2 地 图 绘 制 问 题：假 如 我 们 在 地 图 边 界 获 取 了 一 些 边 界 点 的 坐 标，连 接 这 些 边 界 点 形 成 闭 合 曲 线，可 以 用 来 近 似 表 示 真 实 边 界 线，如 何 更

2、 准 确 地 逼 近 真 实 边 界 线？函数查表与地图边界线绘制(2.3457)=？第1页/共43页如何更准确地逼近真实边界线？第2页/共43页插值在数码图像放大中的应用引例3 图像插值放大：数码相机运用插值的方法可以创造出比传感器实际像素更多的图像，这种处理称为“数码变焦”。106*40原始图像：原始图像：左边：左边：最近邻插值最近邻插值放大放大450%右边：右边：双三次插值双三次插值放大放大450%第3页/共43页插值在图像三维重建中的应用Surface recostruction from scattered points cloud 第4页/共43页分段线性插值和拉格朗日插值分段线性

3、插值：分段线性插值：用直线用直线(线性线性)连接数据点列上相邻的两点。连接数据点列上相邻的两点。比如比如在两点在两点xi-1,xi上线性插值函数为上线性插值函数为拉格朗日插值：拉格朗日插值：用用n次拉格朗日插值多项式次拉格朗日插值多项式连接数据点列上相邻的连接数据点列上相邻的n+1个点。个点。Pszjs71第5页/共43页拉格朗日插值基函数的构造比如 在三个点x0,x1,x2上lagrange插值函数为(线性插值是拉格朗日插值最简单的情形)第6页/共43页分段三次埃尔米特插值条件数分 段 三 次 埃 尔 米 特 插 值：线 性 插 值 在 每 一 小 段 上(两 点 之 间)，用 到 2 个

4、条 件 q(xi)=yi，所 以 确 定 了 一 个 线 性 插 值 函 数；三 次 埃 尔 米 特 插 值 在 每 一 小 段 上，用 到 4 个 条 件 q(xi)=yi，q (xi)=yi,所 以 确 定 一 个 3 次 多 项 式 插 值 函 数。分 段 插 值 主 要 是 为 了 避 免 高 次 插 值 可 能 出 现 的 大 幅 度 振 荡 现 象，在 实 际 应 用 中 通 常 采 用 分 段 低 次 插 值 来 提 高 近 似 程 度，比 如 可 用 分 段 线 性 插 值 或 分 段 三 次 埃 尔 米 特 插 值 来 逼 近 已 知 函 数，但 它 们 的 总 体 光 滑

5、性较 差，为 了 克 服 这 一 缺 点，三 次 样 条 插 值 成 为 比 较 理 想 的 工 具。第7页/共43页三次样条(spline)插值的概念 样条的概念出自工程设计和机械加工(飞机、船舶外形曲线设计)中的绘图工具(曲线尺)，简单说就是具有连续二阶导数的三次插值多项式函数。第8页/共43页三次样条(spline)插值的条件数 首先从段数n=2分析：我们知道在每一小段的三次多项式有4个系数，所以如下图，总共需要有4*2=8个方程来确定；由q(xi)=yi可以确定2*2=4个方程，又由内部节点q1(xi)=q2(xi)和q1(xi)=q2(xi)可以确定2*(2-1)=2个方程，看来剩下

6、的8-(4+2)=2个方程只有靠外部给定(边界条件)了q1q2x0 x1x2第9页/共43页一维曲线等距插值函数interpinterps syntaxOne-dimensional r times longer data interpolation y=interp(y,r)题例题例 在原始数据点中增倍插值在原始数据点中增倍插值x=0:0.001:1;y=sin(2*pi*30*x)+sin(2*pi*60*x);yi=interp(y,4);subplot(1,2,1);stem(y(1:30);title(Original Points);subplot(1,2,2);stem(yi(1

7、:120);title(Interpolated Points);第10页/共43页一维曲线等距插值函数interp1interp1s syntaxOne-dimensional data interpolation yi=interp1(x,y,xi,method)nearest Nearest neighbor interpolationlinear Linear interpolation(default)spline Cubic spline interpolationcubic Piecewise cubic Hermite interpolation题例题例 在一天在一天24小时内

8、小时内,从零点开始每间隔从零点开始每间隔2小时测小时测得的环境温度，推测在得的环境温度，推测在15点点6分的的温度分的的温度x=0:2:24;y=12,9,9,10,18,24,28,27,25,20,18,15,13;plot(x,y,-ro);hold on;xi=15.1;yi=interp1(x,y,xi,spline),xi=0:1/3600:24;yi=interp1(x,y,xi,spline);plot(xi,yi,b-);第11页/共43页二维曲面等距插值函数interp2interp2s syntaxTwo-dimensional data interpolation ZI

9、=interp2(X,Y,Z,XI,YI,method)nearest Nearest neighbor interpolationlinear Bilinear interpolation(default)spline Cubic spline interpolationcubic Bicubuc interpolation第12页/共43页二维曲面等距插值函数interp2动画展示：三维空间中的曲面等距格点动画展示：三维空间中的曲面等距格点第13页/共43页二维曲面等距插值函数interp2题例题例 粗糙山顶曲面的平滑处理粗糙山顶曲面的平滑处理(等距情形等距情形)load mountain

10、.mat%载入山顶地形数据载入山顶地形数据mesh(x,y,z)%绘制原始山顶地形图绘制原始山顶地形图xi=linspace(0,5,50);yi=linspace(0,6,80);xii,yii=meshgrid(xi,yi);zii=interp2(x,y,z,xii,yii,spline);%三次样条插值三次样条插值figure;surf(xii,yii,zii)%绘制平滑处理后的山顶曲面绘制平滑处理后的山顶曲面hold on;xx,yy=meshgrid(x,y);plot3(xx,yy,z+0.1,ob);第14页/共43页二维曲面等距插值函数interp2题例题例 粗糙山顶曲面的平

11、滑处理粗糙山顶曲面的平滑处理(等距情形等距情形)第15页/共43页二维曲面散乱插值函数griddatagriddatas syntaxData interpolation for scattered points ZI=griddata(x,y,z,XI,YI)XI,YI,ZI=griddata(x,y,z,xi,yi).=griddata(.,method)linear Triangle-based linear interpolationcubic Triangle-based cubic(default)nearest Nearest neighbor v4 MATLAB 4 gridd

12、ata methodMATLAB二维插值函数二维插值函数griddata，可以将平面或曲面上的可以将平面或曲面上的散乱点散乱点插值为插值为规则网格规则网格第16页/共43页二维曲面散乱插值函数griddata题例题例 粗糙山顶曲面的平滑处理粗糙山顶曲面的平滑处理(散乱情形散乱情形)rand(seed,0)x=rand(100,1)*4-2;y=rand(100,1)*4-2;z=x.*exp(-x.2-y.2);plot3(x,y,z,o);hold onti=-2:.25:2;XI,YI=meshgrid(ti,ti);ZI=griddata(x,y,z,XI,YI);mesh(XI,YI,

13、ZI);第17页/共43页二维曲面散乱插值函数griddata题例题例 墨西哥草帽的平滑处理墨西哥草帽的平滑处理(散乱情形散乱情形)x=rand(100,1)*16-8;y=rand(100,1)*16-8;r=sqrt(x.2+y.2)+eps;z=sin(r)./r;plot3(x,y,z,.,MarkerSize,15)hold onxlin=linspace(min(x),max(x),33);ylin=linspace(min(y),max(y),33);X,Y=meshgrid(xlin,ylin);Z=griddata(x,y,z,X,Y,cubic);mesh(X,Y,Z);a

14、xis tight;第18页/共43页南半球气旋变化的可视图形第19页/共43页山区地貌的可视化图形第20页/共43页水塔用水量估计通用程序通用程序通用程序tbp69.m可近似计算时间段内的用水量可近似计算时间段内的用水量格式为：格式为：tbp69(ts,tf)其中其中ts为起点时间，为起点时间，tf为终点时间为终点时间第21页/共43页实验一：水塔用水量估计实验一：水塔用水量估计Thats all3Q!第22页/共43页第四讲 插值与拟合之拟合(下)内容：拟合是离散函数逼近的重要方法，利用它 可通过函数在有限个点处的取值状况，拟 合出近似替代函数，进而估算出函数在其 他点处的近似值。目的：学

15、习拟合的基本思想和方法，掌握Matlab 的多项式/一般拟合函数/曲线拟合工具箱要求：掌握Matlab拟合函数，处理拟合应用问题了解基于最小二乘法则拟合的基本思想掌握拟合函数 polyfit lsqcur vefit cur vefit掌握cf tool曲线拟合工具箱(多目标函数多法则)第23页/共43页关于数据拟合的两个要素.在 工 程 实 践 和 科 学 计 算 中，用 某 种 经 验 函 数 解 析 式 y=f(x)来 近 似 刻 画 采 集 数 据(x,y)之 间 的 关 系 的 方 法 就 叫 拟 合，所 谓“拟 合”有 “最 贴 近”之 意 。与 插 值 不 同，拟 合 的 主 要

16、 目标 是 要 离 散 点 尽 量 靠 近 拟 合 函数。一 般 过 程 是，我 们 首 先 根据 采 样 点 的 散 点 分 布 图，大 致推 测 x 与 y 之 间 的 经 验 函 数 形 式(比 如 多 项 式、指 数 函 数 等)，然 后 依 据 某 种 法 则(比 如 最 常 用 的 最 小 二 乘 法 则)，确 定 出 的 经 验 函 数 解 析 式 中 的 待 定 参 数。其 中 经 验 函 数 和 拟 合 法 则 是 拟 合 的 两 个 关 键 要 素！第24页/共43页引 例 1 化 合 物 浓 度 随 时 间 变 化 的 规 律：与 插 值 面 临 的 问 题 相 似，我

17、们 被 要 求 去 求 解 或 预 测 表 格 中 没 有 的 因 变 量 对 应 值，与 插 值 的 解 决 思 路 不 同，我 们 试 图 获 得 比 较 完 备 的 解 决 方 案：设 计 并 求 出 离 散 数 据 点 的 近 似 替 代 函 数，有 了 近 似 函 数 解 析 式，就 可 以 进 一 步 代 值 计 算 或 作 图 分 析。为 了 揭 示 浓 度 y 与 时 间 t 之 间 呈 现 的 函 数 规 律，我 们 首 先 作 出 散 点 图，帮 助 分 析 和 设 计 经 验 函 数化合物浓度随时间变化的规律第25页/共43页化合物浓度随时间变化的规律如图,化合物浓度y随

18、时间t大致呈抛物线状（二次函数）变化，这种分析和判断来自已有经验.t=1:16;c=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6;plot(t,c,-ro)第26页/共43页化合物浓度随时间变化的规律经验函数形式：经验函数形式：已经拟定为多项式函数：已经拟定为多项式函数：y=at2+bt+c剩下的工作是确定拟合原则：剩下的工作是确定拟合原则：可选的法则很多，其中最常用的是最小二乘法则可选的法则很多，其中最常用的是最小二乘法则(method of Least Squares)，即，即各点偏差平方和最

19、小各点偏差平方和最小高斯和勒让德关于最小二乘法的发明权第27页/共43页化合物浓度随时间变化的规律对 经 验 函 数 形 式 确 定 的 补 充 说 明：拟 合 函 数 解 析 式 选 用 什 么 形 式(用 多 项 式，还 是 用 幂 函 数?)，主 要 取 决 于 采 样 点 的 分 布 无 疑，那 么 如 何 求 出 这 些 含 有 待 定 参 数 的 解 析 式 呢？把 各 点 偏 差 的 平 方 和 最 小 作 为 一 个 目 标 函 数，实 际 上 考 虑 为 极 值 问 题，极 值 点 导 数 为 零。具 体 计 算 时，我 们 在 把经 验 函 数 用 一 系 列 拟 合 基

20、函 数 线 性 表 出 同 时，在 j 个 采 样 点 对 待 定 参 数 C j 求 偏 导(=0),获 取 j 个 方 程，进 而 解 出 C j，具 体 参 见 数 值 计 算 S Z J S p 9 0 9 1 在 本 例 中，已 经 拟 定 拟 合 的 目 标 函 数 为 多 项 式 函 数：y=a t2+b t+c ，所 以 只 要 解 出 三 个 待 定 参 数 a,b,c，问 题 即 获 解 决 第28页/共43页基于最小二乘的多项式拟合函数基 于 最 小 二 乘 的 多 项 式 拟 合 函 数 p o l y f i t：P o l y n o m i a l c u r v

21、 e f i t t i n g .S y n t a x：p =p o l y f i t(x,y,n)其 中 n 是 拟 合 多 项 式 的 阶 数，不 能 超 过（散 点 数 据 对 数-1）下 面 回 到 化 合 物 浓 度 随 时 间 变 化 的 引 例：t=1:1 6;c=4 6.4 8 8.4 9.2 8 9.5 9.7 9.8 6 1 0 1 0.2 1 0.3 2 1 0.4 2 1 0.5 1 0.5 5 1 0.5 8 1 0.6 ;p l o t(t,c,k o);h o l d o n;%作 散 点 图p 2=p o l y f i t(t,c,2);y 2=p o

22、l y 2 s t r(p 2,t ),%作 多 次 拟 合 比 较p 5=p o l y f i t(t,c,5);y 5=p o l y 2 s y m(p 5,t ),f=i n l i n e(y 5)t i=0:.0 0 1:2 0;p l o t(t i,p o l y v a l(p 2,t i),b-,t i,f(t i),r-);d i s p(化 合 物 在 刻 度 1 1.2 的 浓 度 近 似 值 为 ,n u m 2 s t r(f(1 1.2)d i s p(化 合 物 在 刻 度 1 7.8 的 浓 度 预 测 值 为 ,n u m 2 s t r(f(1 7.8

23、)s t e m(1 1.2 1 7.8 ,f(1 1.2)f(1 7.8),r );x l a b e l(时 间 t );y l a b e l(化 合 物 浓 度 c );t i t l e(化 合 物 浓 度 随 时 间 变 化 的 规 律 )第29页/共43页引 例 2 确 定 医 用 薄 膜 渗 透 率 的 数 学 模 型：某 种 医 用 薄 膜 允 许 一 种 物 质 分 子从 高 浓 度 溶 液 V B 穿 过 薄 膜 向 低 浓 度溶 液 V A 中 扩 散。通 过 单 位 面 积 膜 S分 子 扩 散 的 速 度 与 膜 两 侧 溶 液 的 浓度 差 成 正 比，比 例 系

24、 数 K 表 示 薄 膜被 该 物 质 分 子 穿 透 的 能 力，称 为 渗 透 率，定 时 测 量 薄 膜 V B 侧 的 溶 液 浓 度 值 C B，以 此 确 定 K 的 值 V A=V B=1 0 0 0 c m3,S=1 0 c m2,容 器 的 B 部 分 溶 液 浓 度 C B 的 测 试 结 果 如 下 表：(C B 单 位 为 m g/c m3 )确定医用薄膜渗透率的数学模型第30页/共43页确定医用薄膜渗透率的数学模型由质量守恒考察由质量守恒考察 t,t+t 时间段时间段B向向A中渗透物质：中渗透物质：VA*CA(t+t)-VA*CA(t)=SKCB(t)-CA(t)t

25、推出推出dCA(t)/dt=SK/VA*CB(t)-CA(t)两边除以两边除以 t,令令 t0又由质量守恒考察整个容器中物质总量始终不变：又由质量守恒考察整个容器中物质总量始终不变：VA*CA(t)+VB*CB(t)=VA*aA+VB*aB 推出推出CA(t)=aA+VB/VA*aB-VB/VA*CB(t)代入上式代入上式2 推出推出dCB(t)/dt=-SK(1/VA+1/VB)CB(t)+SK(aA/VB+aB/VA)CB(0)=aB 初值条件初值条件 此带初值微分方程可由此带初值微分方程可由dsolve求解求解在上式中，已知的包括在上式中，已知的包括VA，VB，S以及一组以及一组t和和C

26、B(t)值值未知的包括未知的包括aA，aB，K，下面通过数据拟合确定渗透率，下面通过数据拟合确定渗透率K第31页/共43页确定医用薄膜渗透率的数学模型在上式中，代入已知值在上式中，代入已知值VA=VB=1000cm3，S=10cm2令令a=(aA*VA+aB*VB)/(VA+VB)，b=VA(aB-aA)/(VA+VB)简化之后的表达式为：简化之后的表达式为：CB(t)=a+b*exp(-0.02*k*t)编写被调编写被调M文件文件 tbp79.m function CB=tbp79(x,t)CB=x(1)+x(2)*exp(-0.02*x(3)*t);编写主调编写主调M文件文件 fittbp

27、79.m（片段）（片段）x=curvefit(tbp79,x0,t,CB)%curvefit拟合及图像拟合及图像x=lsqcurvefit(tbp79,x0,t,CB)%lsqcurvefit拟合及图像拟合及图像求解结果：求解结果：a=x(1)=0.0070;b=x(2)=-0.0030;k=x(3)=0.1012进一步求解：进一步求解：aA=0.01；aB=0.004最终数学模型：最终数学模型：CB(t)=0.007-0.003*exp(-0.002*t)第32页/共43页基于最小二乘的一般拟合函数基 于 最 小 二 乘 的 一 般 拟 合 函 数 l s q c u r v e f i t

28、：S o l v e n o n l i n e a r c u r v e-f i t t i n g (d a t a-f i t t i n g)p r o b l e m s i n t h e l e a s t-s q u a r e s s e n s e.S y n t a x：x =l s q c u r v e f i t(f u n,x 0,x d a t a,y d a t a)x,r e s n o r m =l s q c u r v e f i t(.)范 例:f u n c t i o n F=m y f u n(x,x d a t a)F=x(1)*x d a

29、 t a.2+x(2)*s i n(x d a t a)+x(3)*x d a t a.3;%下 面 是 主 调 函 数 f i t m y f u n.m 的 部 分x d a t a =3.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 1 0.0 5.4 ;y d a t a =1 6.5 1 5 0.6 2 6 3.1 2 4.7 2 0 8.5 9.9 2.7 1 6 3.9 3 2 5.0 5 4.3 ;x 0 =1 0,1 0,1 0 ;%S t a r t i n g g u e s s x,r e s n o r m =l s q c u r v e f i t(

30、m y f u n,x 0,x d a t a,y d a t a)第33页/共43页MATLAB工具箱的版本更新.同名升级同名升级换名升级换名升级第34页/共43页CFTool曲线拟合工具箱简介 基 于 M AT L A B 的 曲 线 拟 合 问 题，已 经 提 供 独 立 的 t o o l b o x 供 调 用，该 t o o l b o x 采 用 G U I 界 面，功 能 强 大，下 面 简 单 介 绍 如 何 使 用 该 To o l b o x 解 决 一 般 曲 线 拟 合 问 题。在 c o m m a n d w i n d o w 中 键 入 指 令 c f t o

31、 o l 即 可 启 动 曲 线 拟 合 工 具 箱。在 该 集 成 环 境 里 面，可 以 实 现 多 种 经 验 函 数，多 种 法 则 的 曲 线 拟 合，实 时 绘 制 图 像 并 进 行 误 差 分 析。需 要 注 意 的 是：在 进 入 C u r v e F i t t i n g To o l b o x 环 境 进 行 曲 线 拟 合 之 前，需 要 预 先 在 w o r k s p a c e 输 入 或 载 入 供 拟 合 的 数 据 源第35页/共43页CFTool-选择Data导入数据下面还是以引例的采样数据为例，进行演示：下面还是以引例的采样数据为例，进行演示：t

32、=1:16;y=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6;cftool导入数据导入数据绘制散点图绘制散点图第36页/共43页CFTool-选择Fitting拟合数据进行拟合进行拟合拟合方法拟合方法结果和误差分析结果和误差分析 这里可供选择的这里可供选择的拟合类型拟合类型和和可选参可选参数数比较多比较多,包括多项式函数包括多项式函数,指数函数指数函数,幂函数等幂函数等,如何确定如何确定最优的方案最优的方案？第37页/共43页CFTool-拟合效果评价指标SSE-The sum of square

33、s due to error.This statistic measures the deviation of the responses from the fitted values of the responses.A value closer to 0 indicates a better fit.R-square-The coefficient of multiple determination.This statistic measures how successful the fit is in explaining the variation of the data.A valu

34、e closer to 1 indicates a better fit.Adjusted R-square-The degree of freedom adjusted R-square.A value closer to 1 indicates a better fit.It is generally the best indicator of the fit quality when you add additional coefficients to your model.RMSE-The root mean squared error.A value closer to 0 indi

35、cates a better fit.第38页/共43页CFTool-选择Analysis分析数据分析数据分析数据 这里可以计算拟合函数在采样点的这里可以计算拟合函数在采样点的一阶二阶导函数值，并绘制相应的函一阶二阶导函数值，并绘制相应的函数图像数图像第39页/共43页CFTool-拟合薄膜渗透数学模型导入数据：导入数据：load bmstl.matcftool拟合数据：拟合数据：syms a b c d tf=a*exp(b*t)+c*exp(d*t)cbvst=subs(f,a,b,c,d,0.008839,-0.0001336,-0.004818,-0.001422);f=inline(cbvst)figuret=100:1:1000;plot(t,f(t),b-)第40页/共43页CFTool-拟合薄膜渗透数学模型分析数据：分析数据：第41页/共43页实验一：实验一：P87 MalthusP87 Malthus人口模型人口模型实验二：实验二：P87 P87 旧车价格预测旧车价格预测Thats all3Q!第42页/共43页感谢您的观看。第43页/共43页