第七讲MATLAB在插值与逼近中的应用

1第七讲MATLAB在插值与逼近中的应用1插值与逼近1.1为什么要逼近数学上来讲，逼近就是在精度要求的范围内对要研究函数给出近似的函数值，甚至函数表达式。为什么我们不直接计算要研究的函数或函值本身？理由如下：用给定函数表达式计算函值很困难甚至根本不可能。如，sinx、tgx、Inx等。由实验与测量得到的变量间对应关系常常是一函数值表（今后我们也称为表列函数）。但表所表示函在表某个中间位置的函数值却是无法知道的。函数可能被隐含地定义，而事实上又不能用一个直接规律给出。例如，由方程ey+y+sinx=0确定的隐含数。计算逼近函数的值往往比计算函值本身更快。特别地，当原来函数以无穷级数的形式给出，只能如此。计算机存储量有限，而其计算量相对来说却很大，从某种意义上来讲，逼近实际上也是为了取长补短。如，我们不可能将所有的sinx的值都存在计算机内，但我们将会看到，利用琏近我们的却可以很方便地算出任一点的函数值。实际应用中，只要函数值符合某一个精度要求也就够了。1.2逼近的分类逼近函数是为了更方便地计算函数，更简单地表达函数。因此，常用一些简单函数或这些简单函数的线性组合来逼近。通常的逼近形式有：我们称ф（x）,i＝0，1，2，…，m为逼近函数，f(x)称为逼近函数。1.3逼近的原则已知函数f(x)在n+1个点xi（i＝0，1，2，…，n）的函数值为f(xi)（i＝0，1，2，…，n）。要求出f(x)的逼近函数g(x)，则要选定逼近基函数，确定上式中的常数ai（i＝0，1，2，…，m）。基函数选定往往跟实际问题有关；而确定常数ai（I＝0，1，2，…，m）以保证逼近函数g(x)能更近似地表示函数f(x)，则是我们这里要解决的问题。为此，就要首先给出一个准则，来描述“更近似”。;)(,)()(/)(4]cossin[32;p(x)1000m0kmjjjmniiinmnmkkkkkxbxpxaxpxpxpkxbkxaxa，其中式之比．有理分式：两个多项；．三角多项式：项式来逼近；的区间上。用不同的多．分段多项式：在不同．多项式：miiixmkkkxaxfeInxxxxba00)()(6321sin6;)exp(5为中逼近形式可统一表示．．．上述等的线性组合。，，，．其它形式，如，．指数函数：2定义距离：其中，p0为一实数。则“更近似“即指“e更小“。因此，确定ai（i＝0，1，2，…，m）使得e取得最小即可。e称为逼近误差。若p=1，称为一致逼近，p=2，称为平方逼近。从上式不难看出，就此式而言，e最好的最小值为零，此时，g(xi)=f(xi)，逼近函数g(x)恰好经过所有n+1个已知点（xi，f(xi)），（i＝0，1，2，…，n）。1.4什么叫插值给定n+1个数据点（x0,y0），(x1,y1)，…，(xn,yn)，若逼近函数经过n+1个数据点，即在已知数据点上的逼近误差为零，则称逼近函数为插值函数，简称为插值。若存在P(xi)=yi(i=0,1,…,n)称P(x)为y=f(x)的插值函数，求插值函数P(x)的方法称为插值法。主要算法有Lagrange插值、Newton插值、分段线性插值、Hermite插值及三次样条插值等。1.5Lagrange插值1.5.1线性插值过函数y=f(x)上的两点(x0,y0)(x1,y1)作一直线p1(x)近似地替代f(x)即：p1(x0)=y0p1(x1)=y1由点斜式1.5.2抛物插值过函数y=f(x)上的三点(x0,y0),(x1,y1),(x2,y2)作一抛物线p2(x)近似地替代f(x)即：p2(x0)=y0p2(x1)=y1p2(x2)=y2作二次式l0(x),使其满足l0(x0)=1,l0(x1)=0,l0(x2)=0,易推出：10100101101001000)0(01010)(1yxxxxyxxxxyxxxxxxxxyyxxxxyyyxp1)1(1,0)0(10)1(0,1)0(0010)(1,101)(0xlxlxlxlxxxxxlxxxxxl则令1)(10)(0)(1yxlyxlxppmipixfxge10|)()(|＝3同理则1.5.3Lagrange插值设函数y=f(x)在给定的两两互异的节点x0,x1,…,xn上的函数值为y0,y1,…,yn,求作一个次数≤n的多项式使它满足这就是Lagrange插值多项式1.5.4Lagrange插值的流程图))(())(()(1)(l))(()(l201021000210xxxxxxxxxlxxxxxcx得由假设))(())(()(2101201xxxxxxxxxl))(())(()(1202102xxxxxxxxxlnnnxaxaxaaxp...)(2210n0,1,2,...,i)(iinyxpninijjjijinijjjijnijjjniiijniiixxxxyxxxxxxxxcxxxxxxxxcxlxyxlx00n0i01100in)(p)(l)())...()()...(()(ji1ji0)(l)()(p所以其中假设2)(21)(10)(0)(2yxlyxlyxlxp41.5.5Lagrange编程functiony=lagrange(x0,y0,x)%lagrangeinsertn=length(x0);p=0;fori=1:nl=1.0;forj=1:nifj~=il=l*(x-x0(j))/(x0(i)-x0(j));end5endp=p+l*y0(i);endy=p;例例给出f(x)=e-x的数值表，用lagrange插值计算e-0.2的近似值x0.100.150.250.30ln(x)0.9048370.8607080.7788010.740818x0=[0.100.150.250.3];y0=[0.9048370.8607080.7788010.740818];y=lagrange(x0,y0,0.2)y=0.818730166666667例给出f(x)=lnx的数值表，用lagrange插值计算ln0.54的近似值x0.40.50.60.70.8ln(x)-0.916291-0.693147-0.510826-0.356675-0.223144x0=[0.4:0.1:0.8];y0=[-0.916291-0.693147-0.510826-0.356675-0.223144];y=lagrange(x0,y0,0.54)y=-0.6161427152y=log(0.54)y=-0.6161861394238171.6分段插值（将插值区间划分小区间[xi,xi+1]，在每个小区间构造插值多项式pi(x)，将每个小区间插值多项式pi(x)拼接）1.6.1高次插值的Runge现象对于函数x=[-5:1:5];y=1./(1+x.^2);x0=[-5:0.1:5];y0=1./(1+x0.^2);y1=lagrange(x,y,x0);plot(x0,y1,'--or')holdonplot(x0,y0,'-b')211)(xxf6-505-0.500.511.52clfx1=[-5:1:0];y1=1./(1+x1.^2);x2=[0:1:5];y2=1./(1+x2.^2);x3=[-5:0.1:0];x4=[0:0.1:5];x0=[-5:0.1:5];y0=1./(1+x0.^2);y3=lagrange(x1,y1,x3);y4=lagrange(x2,y2,x4);plot(x3,y3,'or',x4,y4,'+r')holdonplot(x0,y0,'-b')-50500.10.20.30.40.50.60.70.80.911.7MATLAB插值函数插值函数及其功能函数功能描述interp1一维插值interp2二维插值interp3三维插值7interpft一维快速傅立叶插值interpnN-D维插值spline三次样条数据插值pchip保形分段三次插值（分段三次Hermit插值多项式）yi=interp1(x,y,xi,method)methodsare:'nearest'-nearestneighborinterpolation'linear'-linearinterpolation'spline'-piecewisecubicsplineinterpolation(SPLINE)'pchip'-shape-preservingpiecewisecubicinterpolation'cubic'-sameas'pchip''v5cubic'-thecubicinterpolationfromMATLAB5,whichdoesnotextrapolateanduses'spline'ifXisnotequallyspaced.x=0:8;y=sin(x);xi=0:0.1:8;ni=interp1(x,y,xi,'nearest');li=interp1(x,y,xi,'linear');si=interp1(x,y,xi,'spline');ci=interp1(x,y,xi,'cublic');plot(x,y,'o',xi,ni,'ms',xi,li,'b*',xi,ci,'k+',xi,si,'r-')legend('原数据','nearest','linear','cublic','spline')8012345678-1-0.500.511.5原原原nearestlinearcublicspline1.8样条曲线的MATLAB实现1.8.1样条曲线在工程实践与科学中的应用样条曲线拟合问题：由试验或观测得到了一批数据点，要求用一个函数近似地表明数据点间的函数关系，并画出函数的样条曲线。样条曲线插值问题：由试验、观测或计算得到了若干离散点组成的点列，要求用光滑的样条把这些离散点联结起来。样条曲线逼近问题：在样条曲线形状设计中，给定了折线轮廓，要求用样条曲线逼近这个折线轮廓。MATLAB专门提供了样条工具箱，以便能够方便地处理各种数据，生成样条曲线。1.8.2Splinetool的使用功能：用一系列方法生成各种样条曲线语法：splinetool(x,y)(回车)根据数组x,y在图形用户界面下生成样条曲线。splinetool(回车)启动图形用户界面9选择数据输入的方式（如选第一项Provideyouowndata）x=linspace(0,2*pi,31)y=cos(x)10选择各种生成样条曲线的方法：三次样条插值法，可选择各种约束条件平滑样条法，可选择函数的阶数(4或6),可改变精度的允许值最小二乘法，可选择函数的阶数(1到14)样条插值法，可选择函数的阶数(1到14)选择方法的辅助图形（通过View菜单）显示样条函数的一阶导数曲线图显示样条函数的二阶导数曲线图显示误差曲线图11第一个点的三次导数和第二点的三次导数一样；最后一个点的三次导数和倒数第二个点一样'complete'or'clamped'第一类边界条件（缺省边界条件）第二类边界条件周期（第三类）边界条件使最后一个点的二阶导数等于零自然边界条件用所给定四个点通过Lagrange插值法生成三次样条曲线第一类边界条件：给定函数在x0,xn端点处的一阶导数，第二类边界条件：给定函数在x0,xn端点处的二阶导数，当函数在x0,xn端点处的二阶导数等于0时，称为自然边界条件，此时的样条函数称为自然样条函数第三类边界条件：设f(x)是周期函数，并设x0,xn是一个周期，给定函数在x0,xn端点处的一阶导数和二阶导数相等通过菜单中各项选择，得到所需样条曲线。1.9最小二乘拟合对于有限区间[a,b]上的连续函数f(x),找到一个多项式p(x),使得f(x)-p(x)的欧几里德