《2022年插值与拟合》由会员分享,可在线阅读,更多相关《2022年插值与拟合(17页珍藏版)》请在金锄头文库上搜索。
1、山东科技大学毕业论文-17 - 第四章插值与拟合插值与拟合是来源于实际, 又广泛应用于实际的两种重要的方法。随着计算机的不断发展及计算水平的不断提高,他们已在国民生产和科学研究等方面扮演着越来越重要的角色。在测量上观测值不等间隔和后续值的预测及观测值的曲线拟合,都需要插值与拟合来完成。本节将对插值法中重要的LAGRANGE 插值,分段线性插值, HERMITE 插值及三次洋条插值做介绍, 说明利用 MATLAB 处理此类问题的方法。 另外还对拟合中最为重要的最小二乘法拟合和快速FURIER 变换加以介绍。第一节插值运算(一) LAGRANGE 插值1方法介绍对给定的N 个插值节点x1,x2,x
2、3, xn 及对应的函数值y1,y2,y3,yn,利用(n-1)次 LAGRANGE 插值多项式公式 ,则对插值区间内的任意 X 的函数值 Y 可通过下式求得 . Y(x)=nkYk1( nkjj,1xxxx) 2MATLAB 的实现Lagrange.m Function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - -
3、 - 第 1 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-18 - for j=1;n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s; end3题目举例例: 给出 F(x)=ln(x) 的数值表 ,用 LAGRANGE 插值计算 ln0.54 的近似值解: 在 MATLAB 命令窗口中输入x=0.4:0.1:0.8; y=-0.916291 -0.693147 -0.510826 -0.356675 -0.223144; lagrange(x,y,0.54) ans= -0.6
4、16143 说明:同精确解 ln(0.54)=-0.616186比较起来 ,误差还是可以接受的 ,特别是在工程应用之中 . (二) RUNGE 现象的产生和分段线性插值1方法介绍x 0.4 0.5 0.6 0.7 0.8 lnx -0.916291 -0.693147 -0.510826 -0.35667 -0.223144 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-19 - 上面根据区间 a,b上
5、给出的节点做插值多项式ln(x)近似值,一般总认为 ln(x)的次数越高逼近f(x)的精度就越好。但事实并非如此,本世纪Runge 就给出了一个等距节点插值多项式不收敛的例子。他给出的函数为:f(x)=211x它在区间 -5,5上的各阶导数存在,但在此区间上取n 个节点所构造 Lagrange的插值多项式在全区间内并非都收敛的, 而且分散的很厉害。看如下的例子。例:取 n=10,用 Lagrange插值法进行插值计算。在 MATLAB 命令窗口中输入x=-5:1:5; y=1./(1+x.2); x0=-5:0.1:5; y0=lagrange(x,y,x0); y1=1./(1+x0.2);
6、 plot(x0,y0, r )hold on plot(x0,y1-b)得到的插值曲线已经严重的偏离了原曲线为了解决 Runge 问题,引入分段线性插值。2算法分析所谓分段线性插值就是通过插值点用折线连接起来逼近原曲线,这也是计算机绘制图形的基本原理。3MATLAB 的实现实现分段线性插值不需要编制函数程序,MATLAB 自身提供了内部名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-20 - 的功能函数
7、 interp1. Interp1 一维插值yi=interp1(x,y,xi) 对节点 (x,y)插值,求插值点的函数值x 为节点的向量值, y 为对应的节点的函数值。如果为矩阵,则插值对的每一列进行,若y和 y1的长度超出 x或 xi 的长度,则返回 NaN yi=interp1(y,xi) 此格式默认为 x=1:n, n 为向量 y 的长度值或 y 矩阵的 size(y,1) yi=interp1(x,y,xi,method)Method 指定插值的算法,默认为线性算法。其值可为nearest 线性最近项插值linear线性插值spline 立方样条插值cubic 立方插值所有的插值方法
8、要求x 是单调的。 x 也可以是非连续x 等距且连续时可以用快速插值法。例:x=0:0.1:10 y=sin(x); xi=0:0.25:10; yi=interp1(x,y,xi); plot(x,y,o,xi,yi)名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-21 - (三)三次样条插值1方法介绍样条函数可以给出光滑的插值曲线,因此在数值逼近,常微分方程和偏微分方程的数值及科学和工程的计算中均起到
9、重要的作用。下面着重讨论三次样条函数。设区间 a,b上给定有关划分 a=x0x1=x0(j)&9x(k)x0(j+1) l(k)=j; end end end for k=1:m sum=(3*(x0(l(k)+1)-x(k)2/h(l(k)2- 2*(x0(l(k)+1)-x(k)3/h(l(k)3*y0(l(k); sum=sum+(3*(x(k)-x0(l(k)2/h(l(k)2- 2*(x(k)-x0(l(k)3/h(l(k)3)*m(l(k); sum=sum+h(l(k)*(x0(l(k)+1)-x(k)2/h(l(k)2- (x0(l(k)+1-x(k)3/h(l(k)3)*m(
10、l(k)0; s(k)=sum-h(l(k)*(x(k)-x0(l(k)2/h(l(k)2- (x(k)-x0(l(k)3/h(l(k)3)*m(l(k)+1); end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-24 - 例:给定数据表x 27.7 28 29 30 f(x) 4.1 4.3 4.1 3.0 试求三次样条插值函数满足边界条件s(27.7)=s (30)=0解:在 MATLAB 命令
11、窗口中输入:x=27.7 28 29 30; y=4.1 4.3 4.1 3.0; x0=27.7:0.15:30; y=spline2(x,y,0,0,x0); plot(x0,y) 得到的图形为(四)二维插值二维插值基于与一维插值同样的思想,但是针对两个自变量的函数进行插值。二维插值简单地可以理解为连续三维空间地取值运算,如求解随平面位置变化地温度,湿度,气压等。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页,共 17 页 - - - - - - - - - 山东科技大学
12、毕业论文-25 - 二维插值命令形式为zi=interp2(x,y,z,x0,y0,method)式中, x,y 为自变量数组, z 为测量数组, x0,y0 为指定地自变量插值计算点数组, method 是二维插值使用的方法,包括nearest 最近点插值linear 双线性插值(缺省)cubic 三次函数插值二维插值要求所有自变量取值都是单调地。例如,考虑一个平面范围内地温度变化,在多个离散点上测量出当地的温度,数据如下:x= 1 2 3 4 5 6 y= 1 2 3 4 t= 12 10 11 11 13 15 16 22 28 35 27 20 18 21 26 32 28 25 20
13、 25 30 33 32 30 t 以矩阵形式给出 , t 的列数与x 的维数相同,行数与y 的维数相同。用三维绘图命令mesh(x,y,t)名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-26 - 可以画出温度变化的三维空间网格图形,如下图用直线连接相邻两点。如果只考虑一个方向上的温度变化,如取:x01:0.1:6;细化 x 的范围y0= 指定 y 为常值2 t1=intep2(x,y,t,x0,y0)
14、; 线性插值t2=interp2(x,y,t,x0,y0,cubic );三次插值plot(x0,t1,k,x0,t2,k)比较两个插值结果名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-27 - 如果考虑两个方向上的插值,可以运行下列程序:x0=1:0.1:6; y0=1:0.06:4; 编制 M 函数 f1: functiont3=f1(x,y,t,x0,y0) n=length(y0); for
15、i=1:n tt=interp2(x,y,t,x0,y0(i); for j=1:length(tt) t3(i,j)=tt(j); end end 运行该 M 函数,结果存在t3(x0,y0)矩阵,会图:t3=f1(x,y,t,x0,y0); mesh(x0,y0,t3) 绘出的图形如下图。与上图比较可以看出,在x,y 范围细化以后,采用线性二维插值,结果比一维插值更为精确。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页,共 17 页 - - - - - - - - -
16、山东科技大学毕业论文-28 - 第二节曲线拟合在许多应用领域中,人们经常需要从一系列已知离散点上的数据集(x1,y1),(x2,y2), (xn,yn)得到一个解析函数 y=f(x) 。 得到的解析函数 f(x) 应当在原离散点 xi上尽可能接近给定的yi的值。这一过程称为曲线拟合。最常用的曲线拟合是最小二乘法曲线拟合。拟合结果可使误差的平方和最小,即找出使ni1yixif)(2 最小的f(x).如果自己去编程序, 求解最小二乘解, 过程是比较复杂的。 MATLAB提供的 polyfit 函数,根据给定的自变量数组x 和函数数组 y,按照拟合的阶数要求自动的求解满足最小二乘意义的一阶或高阶解析
17、函数f(x),使用很方便。为了说明这个问题,我们取以下函数为例:y=0.5-2*x.2 令自变量 x=0:0.1:1 键入下列命令,得到y 的准确解for i=1:length(x) y(i)=0.5-2*x(i)2; end 键入y 得到y= Columns 1 though 7 0.5000 0.4800 0.4200 0.3200 0.1800 0 -0.2200 Columns 8 through 11 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 12 页,共 17 页
18、- - - - - - - - - 山东科技大学毕业论文-29 - -0.4000 -0.8100 -1.1000 -1.5000 令 y 的每一项数据都有一些偏差,构成待拟合的数据集:y0.52 0.45 0.40 0.35 0.18 0.02 0.25 0.40 0.81 1.1 1.5 分别进行一维和二维曲线拟合:键入命令m=1; fxy1=polyfit(x,y,m) fxy1= -1.9873 0.7991 拟合出的多项式为fxy1=-1.9873x+0.7991,是一个线性方程键入命令:m=2; fxy2=polyfit(x,y,m) fxy2= -2.0350 0.0477 0.
19、4938 拟合出的多项式为fxy2=-2.035x2+0.0477x+0.4938,是一个二阶方程。很明显,二维曲线拟合得到的结果与拟合得到的结果与给定的多项式y=0.5-2*x2 更为接近,一维拟合效果较差。利用多项式估计命令分别计算一.二维拟合结果的数值,并与原函数值绘在一张图上。y1=polyval(fxy1,x); y2=polyval(fxy2,x); plot(x,y,o,x,y1,k: ,x,y2.k)名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 13 页,共 17
20、页 - - - - - - - - - 山东科技大学毕业论文-30 - 由图可以看出,二维拟合结果与原函数值基本重合。为了使用命令,我们必须给出自变量数据组和期望的最佳拟合数据的多项式的阶数。选择不同的阶数,会得到不同的拟合结果,这一点从上面的举例中可以明显看出,如何选择阶数,需要在系统辨识方面有更多的知识。最小二乘法拟合介绍1. 方法介绍在科学实验的统计方法研究中,往往要从一组实验数据(xi,yi )中寻找出自变量 x 和因变量 Y 之间的函数关系y=f(x) 。由于观测数据往往不够准确,因此并不要求y=f(x)经过所有的点( xi,yi),而只要求在给定的点 xi 上误差 ei=f(xi)
21、-yi按照某种标准最小,通常就采用欧氏范数作为误差的量度标准。这就是最小二乘法。 2. MATLAB 实现名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 14 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-31 - 在 MATLAB 中实现最小二乘法拟合通常可以采用两种途径:1利用 polyfit功能函数进行多项式拟合。2利用常用的矩阵除法解决复杂型函数的拟合最小二乘法拟合如表所示的数据:此例用功能函数进行拟合,在MATLAB 命令窗口中输入:x=0.5
22、 1.0 1.52.0 2.5 3.0;y=1.75 2.45 3.81 4.80 7.00 8.60 ;a=polyfitfit (X,Y,2)a= 0.5614 0.8287 1.1560 x1=0.5:0.05:3.0;y1=a(3)+a(2)*x1+a(1)*x2; plot(x,y,*)hold on plot(x1,y1,-r ) 得到下图x 0.5 1.0 1.5 2.0 2.5 3.0 y 1.75 2.45 3.81 4.80 7.00 8.60 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 -
23、 - - - - - - 第 15 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-32 - 用最小二乘法求一个形如Y=A+BX2的经验公式,它与下表所示的数据拟合。下面用另一种方法来求解此拟合问题。用求解矩阵的方法来解, 把 A,B看成未知量,已知Xi,Yi ,求解一超定方程。在 MATLAB 中实现为:x=19 25 31 38 44; y=19.0 32.3 49.0 73.3 97.8; x1=x.2 x1= 361625 961 1444 1963 x1=ones(5,1),x1 x1= 1361 1625 xi 19 25 31 38 44 yi 19.
24、0 32.3 49.0 73.3 97.8 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 16 页,共 17 页 - - - - - - - - - 山东科技大学毕业论文-33 - 1961 11444 11963 ab=x1yab= 0.9726 0.0500 x0=19:0.2:44; y0=ab(1)+ab(2)*0.2; clf plot(x,y, o ) hold on plot(x0,y0, -r )名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 17 页,共 17 页 - - - - - - - - -