拉格朗日插值法

拉格朗日插值法的一些讨论学院：班级：姓名：学号:引言——在数值分析中，拉格朗日插值法是以法国十八世纪数学家约瑟夫·路易斯·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律，而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。拉格朗日插值法最早被英国数学家爱德华·华林于1779年发现，不久后（1783年）由莱昂哈德·欧拉再次发现。1795年，拉格朗日在其著作《师范学校数学基础教程》中发表了这个插值方法，从此他的名字就和这个方法联系在一起。数据建模有两大方法：一类是插值方法，另一类是拟合函数一般的说，插值法比较适合数据准确或数据量小的情形。然而Lagrange插值有很多种，1阶，2阶,…n阶。我们可以利用拉格朗日插值求方程，根据它的程序求原方程的图像。下面我具体介绍分析一下拉格朗日插值的算法设计及应用。具体算法1、基本概念已知函数y=f(x)在若干点ix的函数值iy=ixf（i=0,1,,n）一个差值问题就是求一“简单”的函数p(x)：p(ix)=iy,i=0,1,,n,(1)则p(x)为f(x)的插值函数，而f(x)为被插值函数会插值原函数，0x，1x，2x，...,nx为插值节点，式（1）为插值条件，如果对固定点x求f(x)数值解，我们称x为一个插值节点，f(x)p(x)称为x点的插值，当x[min(0x，1x，2x，...,nx)，max(0x，1x，2x，...,nx)]时，称为内插，否则称为外插式外推，特别地，当p(x)为不超过n次多项式时称为n阶Lagrange插值。2、Lagrange插值公式（1）线性插值)1(1L设已知0x，1x及0y=f(0x),1y=f(1x),)(1xL为不超过一次多项式且满足)(01xL=0y,)(11xL=1y,几何上，)(1xL为过（0x，0y），（1x，1y）的直线，从而得到)(1xL=0y+0101xxyy（x-0x）.（2）为了推广到高阶问题，我们将式（2）变成对称式)(1xL=0l（x）0y+1l(x)1y.其中，0l（x）=101xxxx,1l(x)=010xxxx。均为1次多项式且满足0l（x）=1且1l(x)=0。或0l（x）=0且1l(x)=1。两关系式可统一写成)(iixl=jiji01。（3）（2）n阶Lagrange插值)(xLn设已知0x，1x，2x，...,nx及iy=f(ix)(i=0,1,.....,n),)(xLn为不超过n次多项式且满足iinyxL)(（i=0,1,...n）.易知)(xLn=0l（x）0y+....+)(xlnny.其中，)(xli均为n次多项式且满足式（3）（i,j=0,1,...,n）,再由jx（ji）为n次多项式)(xli的n个根知)(xli=cniijjxx0.最后，由1)()(0nijjjijixxcxlc=nijjjixx0)(1,i=0,1,...,n.总之，)(xLn=iniiyxl0)(，)(xli=.0nijjjijxxxx式为n阶Lagrange插值公式，其中，)(xli（i=0,1,...n）称为n阶Lagrange插值的基函数。3，Lagrange插值余项设0x，1x，2x，...,nx[a,b],f(x)在[a,b]上有连续的n+1阶导数，)(xLn为f(x)关于节点0x，1x，2x，...,nx的n阶Lagrange插值多项式，则对任意x[a,b],).()!1()()()()()1(xnfxLxfxRnnn其中，位于0x，1x，2x，...,nx及x之间（依赖于x），(x)=njjxx0).(Eg1:已知函数表sin6=0.5000,sin4=0.7071,sin3=0.8660,分别由线性插值与抛物插值求sin92的数值解，并由余项公式估计计算结果的精度。解：（1）这里有三个节点，线性插值需要两个节点，根据余项公式，我们选取前两个节点，易知：sin921L(92)=0.5000+645000.07071.0(92-6)=0.5000+0.207132=0.6381截断误差，)92(1R=)492)(692(2)(sinx310615.7361821，得.105.010615.713知结果至少有1位有效数字。（2）易知sin925000.0)3-6)(4-6()33-92)(4-92()92(2L））（（））（（34643926920.7071+8660.04363492692））（（））（（=7071.0985000.092910.8660=0.6434截断误差为：xxR)492)(492)(692(6)(sin)92(2210861..09361861得.105.010861.824知结果至少有两位数字。比较本题精确解sin92=0.642787609...,实际误差限分别为0.0047和0.00062。4，拉格朗日插值算法和程序functionyy=nalagr(x,y,xx)m=length(x);n=length(y);ifm~=n,error('向量x与y的长度必须一致');ends=0;fori=1:nt=ones(1,length(xx));forj=1:nifj~=it=t.*(xx-x(i))/(x(i)-x(j));endends=s+t*y(i);endyy=s;用以上程序的Eg1的结果为x=pi*[1/61/4];y=[0.50.7071];xx=2*pi/9;yy1=nalagr(x,y,xx)yy1=-0.5690x=pi*[1/61/41/3];y=[0.50.70710.866];yy2=nalagr(x,y,xx)yy2=0.8023fplot('sin',[pi/6,pi/3]);holdon;plot(x,y,'o',xx,0.6381,'g^',xx,0.6434,'rv');holdoff;图形为具体例题问题：1、已知916291.0)4.0(f，693147.0)5.0(f，510826.0)6.0(f，分别用线性插值法和二次拉格朗日插值法分别计算)54.0(f，并且说明哪个结果误差更小。2、给出节点03.37)15.3(f，24.7)00.1(f，05.1)01.0(f，03.2)02.1(f，06.17)03.2(f，05.23)25.3(f，作出五次拉格朗日插值多项式，并且求出)09.1(f及其误差。三、实验结果1．（1）线性插值：程序：formatshortx=[0.40.5];y=[-0.916291-0.693147];L0=poly(x(2))/(x(1)-x(2));L1=poly(x(1))/(x(2)-x(1));z=y(1)*L0+y(2)*L1;a=poly2sym(z)b=polyval(z,0.54)结果：a=27893/12500*x-4073206373581395/2251799813685248b=-0.6039（2）二次拉格朗日程序：clcformatshortx=[0.40.50.6];y=[-0.916291-0.693147-0.510826];L0=poly([x(2)x(3)])/((x(1)-x(2))*(x(1)-x(3)));L1=poly([x(1)x(3)])/((x(2)-x(1))*(x(2)-x(3)));L2=poly([x(1)x(2)])/((x(3)-x(1))*(x(3)-x(2)));z=y(1)*L0+y(2)*L1+y(3)*L2;a=poly2sym(z)b=polyval(z,0.54)结果：a=-287266324356479/140737488355328*x^2+1145173905872887/281474976710656*x-1248114652880531/562949953421312b=-0.6153误差：function[R]=wucha(M,x,X)n=length(X);r=1;q=1;fori=1:nr=i*r;%q=(x-X(i))*q;%endR=M*abs(q)/r;结果：（1）x=0.54x=0.5400X=[0.40.5]X=0.40000.5000symsM[R]=wucha(M,x,X)R=7/2500*M（2）x=0.54x=0.5400X=[0.40.50.6]X=0.40000.50000.6000symsM[R]=wucha(M,x,X)R=7/125000*M答：线性插值的误差为：R=7/2500*M；二次拉格朗日插值的误差为：R=7/125000*M显然，二次拉格朗日的误差较小。2．五次拉格朗日clcx=[-3.15-10.011.022.033.25];y=[37.037.241.052.0317.0623.05];L0=poly([x(2)x(3)x(4)x(5)x(6)])/((x(1)-x(2))*(x(1)-x(3))*(x(1)-x(4))*(x(1)-x(5))*(x(1)-x(6)));L1=poly([x(1)x(3)x(4)x(5)x(6)])/((x(2)-x(1))*(x(2)-x(3))*(x(2)-x(4))*(x(2)-x(5))*(x(2)-x(6)));L2=poly([x(1)x(2)x(4)x(5)x(6)])/((x(3)-x(1))*(x(3)-x(2))*(x(3)-x(4))*(x(3)-x(5))*(x(3)-x(6)));L3=poly([x(1)x(2)x(3)x(5)x(6)])/((x(4)-x(1))*(x(4)-x(2))*(x(4)-x(3))*(x(4)-x(5))*(x(4)-x(6)));L4=poly([x(1)x(2)x(3)x(4)x(6)])/((x(5)-x(1))*(x(5)-x(2))*(x(5)-x(3))*(x(5)-x(4))*(x(5)-x(6)));L5=poly([x(1)x(2)x(3)x(4)x(5)])/((x(6)-x(1))*(x(6)-x(2))*(x(6)-x(3))*(x(6)-x(4))*(x(6)-x(5)));z=y(1)*L0+y(2)*L1+y(3)*L2+y(4)*L3+y(5)*L4+y(6)*L5;a=poly2sym(z)b=polyval(z,1.09)结果：a=-1762131072886105/9007199254740992*x^5-1725132155832521/36028797018963968*x^4+1255027041024201/562949953421312*x^3+7895561387872097/2251799813685248*x^2-5312509994841133/1125899906842624*x+4939690858283287/4503599627370496b=2.6381for法：formatshortfunction[C,L]=lagranl(X,Y)n=length(X);fork=1:nV=1;fori=1:nifk~=iV=conv(V,poly(X(i)))/(X(k)-X(i));endendL1(k,:)=V;end