插值和拟合

插值和拟合都是函数逼近或者数值逼近的重要组成部分他们的共同点都是通过已知一些离散点集M上的约束，求取一个定义在连续集合S(M包含于S)的未知连续函数，从而达到获取整体规律的目的，即通过窥几斑来达到知全豹。简单的讲，所谓拟合是指已知某函数的若干离散函数值{f1,f2,…,fn}，通过调整该函数中若干待定系数f(λ1,λ2,…,λ3),使得该函数与已知点集的差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫作非线性拟合或者非线性回归。表达式也可以是分段函数，这种情况下叫作样条拟合。而插值是指已知某函数的在若干离散点上的函数值或者导数信息，通过求解该函数中待定形式的插值函数以及待定系数，使得该函数在给定离散点上满足约束。插值函数又叫作基函数，如果该基函数定义在整个定义域上，叫作全域基，否则叫作分域基。如果约束条件中只有函数值的约束，叫作Lagrange插值，否则叫作Hermite插值。从几何意义上将，拟合是给定了空间中的一些点，找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点；而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。一、概念的引入1.插值与拟合在现实生活中的应用l机械制造：汽车外观设计l采样数据的重新建构：电脑游戏中场景的显示，地质勘探，医学领域（CT）2.概念的定义l插值：基于[a,b]区间上的n个互异点，给定函数f(x)，寻找某个函数去逼近f(x)。若要求φ(x)在xi处与f(xi)相等，这类的函数逼近问题称为插值问题，xi即是插值点l逼近：当取值点过多时，构造通过所有点的难度非常大。此时选择一个次数较低的函数最佳逼近这些点，一般采用最小二乘法l光顾：曲线的拐点不能太多，条件：①二阶几何连续②不存在多余拐点③曲率变化较小l拟合：曲线设计过程中用插值或通过逼近方法是生成的曲线光滑(切变量连续)光顾二、插值理论设函数y=f(x)在区间[a,b]上连续，在[a,b]上有互异点x0,x1,…,xn处取值y0,y1,…,yn。如果函数φ(x)在点xi上满足φ(xi)=yi(i=0,1,2,…,n)，则称φ(x)是函数y=f(x)的插值函数，x0,x1,…,xn是插值节点。若此时φ(x)是代数多项式P(x)，则称P(x)为插值多项式。显然f(x)≈φ(x)，x∈[a,b]1.拉格朗日插值构造n次多项式Pn(x)=yklk(x)=y0l0(x)+y1l1(x)+…+ynln(x)，这是不超过n次的多项式，其中基函数lk(x)=显然lk(x)满足lk(xi)=此时Pn(x)≈f(x),误差Rn(x)=f(x)-Pn(x)=其中∈(a,b)且依赖于x，=（x-x0）(x-x1)…(x-xn)很显然，当n=1、插值节点只有两个xk,xk+1时P1(x)=yklk(x)+yk+1lk+1(x)其中基函数lk(x)=lk+1(x)=2.牛顿插值构造n次多项式Nn(x)=f(x0)+f(x0,x1)(x-x0)+f(x0,x1,x2)(x-x0)(x-x1)+…+f(x0,x1,x2,…,xn)(x-x0)(x-x1)…(x-xn)称为牛顿插值多项式，其中(二个节点，一阶差商)(三个节点，二阶差商)(n+1个节点，n阶差商)注意：由于插值多项式的唯一性，有时为了避免拉格朗日余项Rn(x)中n+1阶导数的运算，用牛顿插值公式Rn(x)=f(x)-Nn(x)=f(x,x0,…,xn)ωn+1(x),其中ωn+1(x)=(x-x0)(x-x1)…(x-xn)3.分段插值------子区间内，避免函数在某些区间失真1）线性插值已知n+1个不同节点x0,x1,…,xn,构造分段一次线性多项式P(x)，使之满足lP(x)在[a,b]上连续lP(xk)=yklP(x)在[xi,xi+1]上是线性函数，P(x)=2）两点带导数插值---避免尖点、一阶连续区间[a,b]上两个互异节点xi,xi+1，已知实数yi,yi+1,mi,mi+1，为了构造次数不大于3的多项式满足条件引入,使之满足可以求出此时=+，其中4.三次样条插值------二阶可导对于给定n+1个不同节点x0,x1,…,xn及函数值y0,y1,…,yn，其中a=x01…n=b。构造三次样条插值函数S(x)。S(x)称为三次样条函数时需满足：lS(x)在[a,b]上二阶导数连续lS(xk)=yk(k=0,1,…,n)l每个子区间[xk,xk+1]上S(x)是三次多项式(k=0,1,…,n)5.例题已知函数y=f(x)的观测值X1234Y0-5-63求三次插值多项式φ(x)及φ(2.5).解：(1)拉格朗日插值P3(x)=y0l0(x)+y1l1(x)+y2l2(x)+y3l3(x)=(-5)+(-6)+3=x3-4x2+3φ(x)≈P3(x)=x3-4x2+3φ(2.5)=2.53-4*2.52+3=-6.375(2)牛顿插值xif(xi)一阶差商二阶差商三阶差商102–5-53–6-1243951N3(x)=f(x0)+f(x0,x1)(x-x0)+f(x0,x1,x2)(x-x0)(x-x1)+f(x0,x1,x2,x3)(x-x0)(x-x1)(x-x2)=0+(-5)(x-1)+2*(x-1)(x-2)+1*(x-1)(x-2)(x-3)=x3-4x2+3三、Matlab在插值中的应用1.Lagrange插值１)方法回顾对给定的n个节点x1,x2,…,xn及对应的函数值y1,y2,…,yn，利用n次Lagrange插值多项式公式，插值区间内任意x的函数值y可以通过下式求出：2)Matlab实现函数Lagrange.mfunctiony=lagrange(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end3)例题的Matlab实现x=[1234];y=[0–5–63];lagrange(x,y,2.5)2.Runge现象及分段线性插值１）Runge现象Runge在本世纪初发现：在[-1,1]上用n+1个等距结点作插值多项式Pn(x)，使其在个结点的值与函数y(x)=1/(1+25x2)在结点的值相等。但在n→∞时，插值多项式Pn(x)在区间中部趋于y(x)。但对于0.726≤∣x∣≤1的x，Pn(x)严重发散。通过下面的例子，以图形的方式体会Runge现象（f(x)=1/(1+x2)）x=[-5:1:5];y=1./(1+x.^2);x0=[-5:0.1:5];y0=lagrange(x,y,x0);y1=1./(1+x0.^2);%绘制图形plot(x0,y0,'--r')holdonplot(x0,y1,'-b')2）Matlab实现分段插值------一维插值interp1lyi=interp1(x,y,xi)对(x,y)进行插值，计算插值点xi的函数值lyi=interp1(y,xi)默认x=1:n，n是向量y的元素个数lyi=interp1(x,y,xi,’method’)指定特定算法插值，method可以是如下字符串Ølinear线性插值Øspline三次样条插值Øcubic三次插值要求：x是单调，但不要求连续等距。如果x连续等距，可以选用快速插值法。调用函数时只需在method前加”*”，如”*spline”3）例题①用一维线性插值解决Runge现象y2=interp1(x,y,x0);plot(x0,y2,'*m')②正弦曲线的插值示例x=0:0.1:10;y=sin(x);xi=0:0.25:10;yi=interp1(x,y,xi);plot(x,y,’o’,xi,yi)3.Hermite插值1)方法介绍已知n个插值节点x1,x2,…,xn及起对应的函数值y1,y2,…,yn和一阶导数值y1’,y2’,…,yn’，则计算插值区域内任意x的函数值y的Hermite插值公式其中：,2)Matlab实现Hermite.mfunctiony=hermite(x0,y0,y1,x)n=length(x0);m=length(x);fork=1:myy=0.0;fori=1:nh=1.0;a=0.0;forj=1:nifj~=ih=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;a=1/(x0(i)-x0(j))+a;endendyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));endy(k)=yy;end3）例题如下数据表，构造Hermite多项式，并求出Sin0.34的近似值x0.300.320.35Sin(x)0.295520.314570.34290Cos(x)0.955340.949240.93937x0=[0.30.320.35];y0=[0.295520.314570.34290];y1=[0.955340.949240.93937];x=[0.3:0.005:0.35];y=hermite(x0,y0,y1,x);plot(x,y)y=hermite(x0,y0,y1,0.34);ysin(0.34)y2=sin(x);holdonplot(x,y2,'--r')4.三次样条插值y3=interp1(x,y,x0,'*spline');y3=spline(x,y,x0);plot(x0,y3,'-g')四、数据拟合1.方法介绍在实际生活中，往往需要从一组实验数据（xi,yi）中寻找出变量x,y之间的函数关系。由于观测数据不可避免出现误差，因此并不需要y=f(x)一定要经过所有的点，而只要求在给定点xi上误差Δi=f(xi)-yi按某种标准达到最小。通常用欧式范数║Δ║2作为误差量度的标准。这就是所谓的最小二乘法。注意：数据拟合与插值的最大区别在于拟合需要给出一个曲线方程的具体解析形式，而插值只需求出该点的内插数值。1)线性拟合线性拟合以最简单的一次线性方程f(x)=a1x+a0拟合数据。按最小二乘法，a1,a0需满足最小，因此可以通过求出此时的a1,a02)超定方程的解法求解方程Ax=b，其中A为m*n阶矩阵，x(1:n)，b(1:m)均为列向量，且mn。由于该超定方程个数多于未知数个数，当增广矩阵的秩大于系数矩阵的秩时无解。现在求其最小二乘解，它就是使余向量rx=b-Ax的谱范数║rx║2=(rxTrx)1/2最小的n维向量。具体解法可以通过求解该方程组的法方程组ATAx=ATb获得。2.Matlab的实现1)线性拟合及多项式拟合ployfit(x,y,i)以最高次为i的多项式拟合数据点(x,y)例1x=[012345];y=[021627077110];coef=polyfit(x,y,1);a1=coef(1),a0=coef(2);ybest=a1*x+a0;s=sum((y-ybest).^2);axis([-1,6,-20,120]);plot(x,y,'*')holdonplot(x,ybest)例2如下给出从二阶到十阶多项式拟合曲线的比较程序，并给出拟合曲线x=[012345];y=[021627077110];xi=0:0.2:5;forn=2:10bb=polyfit(x,y,n);yi=polyval(bb,xi);plot(xi,yi,x,y,'*')title([int2str(n),'次多项式拟合曲线'])gridonpauseend例3在某个实验中得到如下一组数据：x1234567y0.31010.49000.64000.80000.92001.05001.2000已知x,y满足y=kxn，求参数k与n。提示:y=kxnàlny=lnk+nlnxLOG(x)EXP(x)**可线性化的非线性模型模型形式变换后形式变量和参数的变化YXa1a22)超定方程的解法例：用最小二乘法求