数值分析论文(作业)

数值分析--代数插值法的论述姓名：蔺孝宝学号：12023316班级：1203学院：商洛学院数计学院数学与计算科学系日期2014.12.29商洛学院-1-代数插值法1.摘要插值法是函数逼近的重要方法之一，有着广泛的应用。在生产和实验中，函数f(x)或者其表达式不便于计算复杂或者无表达式而只有函数在给定点的函数值(或其导数值)，此时我们希望建立一个简单的而便于计算的函数(x)，使其近似的代替f(x)，有很多种插值法,其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点,常用的插值还有Hermit插值,分段插值和样条插值.这里主要介绍拉格朗日(Lagrange)插值和牛顿(Newton)插值并在MATLAB中的应用操作。【关键字】插值法拉格朗日插值牛顿插值MATLAB正文：一、调用MATLAB内带函数插值1、MATLAB内带插值函数列举如下：interp1interpftinterp2interp3interpnsplinemeshgridndgridgriddata一维数据内插(查表法)使用FFT方法的一维数据内插二维数据内插(查表法)三维数据内插(查表法)多维数据内插(查表法)三次样条内插为三维绘图产生X和Y阵为多维函数和内插产生阵列数据网格2、取其中的一维数据内插函数（interp1）为例，程序如下：其调用格式为：yi=interp1(x,y,xi)yi=interp1(x,y,xi,method)举例如下：x=0:10:100y=[40444652657680828892110];xi=0:1:100yi=interp1(x,y,xi,'spline')3、其他内带函数调用格式为：Interpft函数：y=interpft(x,n)y=interpft(x,n,dim)interp2函数：商洛学院-2-ZI=interp2(X,Y,Z,XI,YI)，ZI=imerp2(Z,ntimes)ZI=interp2(Z,XI,YI)，ZI=interp2(X,Y,Z,XI,YI,method)interp3函数：VI=interp3(X,Y,Z,V,XI,YI,ZI)VI=interp3(V,ntimes)VI=interp3(V,XI,YI,ZI)VI=interp3(…,method)Interpn函数：VI=interpn(X1,X2,X3,…,V,Y1,Y2,Y3,…)VI=interpn(V,ntimes)VI=interpn(V,Yl,Y2,Y3,…)VI=interpn(…,method)Spline函数：yi=spline(x,y,xi)pp=spline(x,y)meshgrid函数：[X,Y]=meshgrid(x,y)[X,Y]=meshgrid(x)[X,Y,Z]=meshgrid(x,y,z)Ndgrid函数：[X1,X2,X3,…]=ndgrid(x1,x2,x3,…)[X1,X2,X3,…]=ndgrid(x)Griddata函数：ZI=griddata(x,y,z,XI,YI)[XI,YI,ZI]=griddata(x,y,z,xi,yi)[…]=griddata(…method)二、两种插值法分析与MATLAB应用1.1拉格朗日插值1.1.1基本原理构造n次多项式Pn(x)=yklk(x)=y0l0(x)+y1l1(x)+…+ynln(x)，这是不超过n次的多项式，其中基函数lk(x)=)...()()...()(()...()()...()(()1110)1110nkkkkkkknkkxxxxxxxxxxxxxxxxxxxx显然lk(x)满足lk(xi)=)(0)(1kiki此时Pn(x)≈f(x),误差Rn(x)=f(x)-Pn(x)=(x))!1()(1)1(nnnf其中∈(a,b)且依赖于x，(x)1n=（x-x0）(x-x1)…(x-xn)很显然，当n=1、插值节点只有两个xk,xk+1时P1(x)=yklk(x)+yk+1lk+1(x)商洛学院-3-其中基函数lk(x)=11kkkxxxxlk+1(x)=kkkxxxx11.1.2优缺点可对插值函数选择多种不同的函数类型，由于代数多项式具有简单和一些良好的特性，故常选用代数多项式作为插值函数。利用插值基函数很容易得到拉格朗日插值多项式，公式结构紧凑，在理论分析中甚为方便，但当插值节点增减时全部插值基函数Lk(x)(k=0,1,…,n)均要随之变化，整个公式也将发生变化，这在实际计算中是很不方便的，为了克服这一缺点，提出了牛顿插值可以克服这一缺点。1.1.3数值实验建立M文件：functionf=Language(x,y,x0)symstl;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等！');return;%检错endh=sym(0);for(i=1:n)l=sym(y(i));for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));end;h=h+l;endsimplify(h);商洛学院-4-if(nargin==3)f=subs(h,'t',x0);%计算插值点的函数值elsef=collect(h);f=vpa(f,6);%将插值多项式的系数化成6位精度的小数end在MATLAB中输入：x=[18316668707270;]y=[23335251434046];f=Language(x,y)plot(x,y)结果为：f=Inf+(-t)*Inf-54329.8*t^2+1503.75*t^3-22.2065*t^4+0.16789*t^5-0.000512106*t^6图形如下：MATLAB实现拉格朗日插值建立如下拉格朗日插值函数：functiony=lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:m商洛学院-5-z=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;end画图程序如下：x=[-5:1:5];y=1./(1+x.^2);x0=[-5:0.001:5];y0=lagrange(x,y,x0);y1=1./(1+x0.^2);plot(x0,y0,'r')holdonplot(x0,y1,'g')注：画出的图形为n=10的图形得到图形如下：n=10的图像商洛学院-6-1.2牛顿插值1.2.1基本原理构造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)称为牛顿插值多项式，其中101010)()(),(xxxfxfxxf(二个节点，一阶差商)202110210),(),(),,(xxxxfxxfxxxf(三个节点，二阶差商)nnnnxxxxxfxxxfxxxf02111010),...,,(),...,,(),...,,((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)1.2.2优缺点牛顿插值法具有承袭性和易变性的特点，当增加一个节点时，只要再增加一项就可以了即],,,[)())(()()(10101kkkkkxxxfxxxxxxxNxN而拉格朗日插值若要增加一个节点时全部基函数都需要重新算过。牛顿插值法既适合于用来计算函数值，也适合于做理论推导，比如说可用来推导微分方程的数值求解公式。1.2.3数值实验1、差商相关公式00100120101012011011011,,,...,,...,,,...,,,,...,,,...,...nnnkkkkkknnNxfxfxxxxfxxxxxxxfxxxxfxxxxfxxxfxxxxxxxxxxxx商洛学院-7-2、程序思想根据实际计算理论，利用Newton插值多项式计算，事先应知道其节点个数。所以，本程序设置在100个节点即最高次可达99次的多项式计算，通过输入节点个数来控制每次要计算的多项式的次数。首先，计算各阶的函数值。根据节点的个数与每阶阶商具有的阶商个数循环计算各阶商的值，并根据相应的节点下标控制输出各阶商值。其次，输出最高阶表达式。由于最高阶表达式是从0位开始计算的，所以，在知道节点下标的基础上，每次以二维数组两下标是否相等来输出相应的阶商，及x的值与节点的起止下标输出相应的表达式。最后，根据输入的节点数计算其对应的函数值。3阶计算表达式，任意输入节点值，根据输入的节点值的大小判断其所在的范围，以输出表达式同样的思想计算节点对应的函数值。最高阶表达式，亦是通过循环计算获得相应的函数值。3、牛顿插值法基本思路与计算步骤：给定插值点序列（))(,iixfx,,,1,0,ni。构造牛顿插值多项式)(uNn。输入要计算的函数点,x并计算)(xNn的值，利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面)(xNn的各项系数恰好又是各阶均差，而各阶均差可用均差公式来计算。为的一阶均差。为的k阶均差。均差表：零阶均差一阶均差二阶均差三阶均差X0f(X0)X1f(X1)f[X0,X1]X2f(X2)f[X1,X2]f[X0,X1,X2]X3f(X3)f[X2,X3]f[X1,X2,X3]f[X0,X1,X2X3]MMMMM牛顿插值法计算步骤：（1）输入n值及（))(,iixfx,,,1,0,ni；要计算的函数点x。（2）对给定的,x由kx商洛学院-8-00010101201101()()(),()(),,()()(),,nnnNxfxxxfxxxxxxfxxxxxxxxxfxxx计算()nNx的值。（3）输出()nNx。建立M文件：function[c,d]=newpoly(x,y)%牛顿插值的MATLAB实现%这里x为n个节点的横坐标所组成的向量，y为纵坐标所组成的向量。%c为所求的牛顿插值多项式的系数构成的向量。n=length(x);%取x的个数。d=zeros(n,n);%构造nXn的空数组。d(:,1)=y';forj=2:nfork=j:nd(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));endendc=d(n,n);fork=(n-1):-1:1c=conv(c,poly(x(k)));%conv求积，poly(x)将该多项式的系数赋给向量。m=length(c);c(m)=c(m)+d(k,k);end（4）测试数据与结果：测试数据：（第三章习题第三题第2题）f(x)=lnx的数值如表所示，构造牛顿插值多项式并求ln0.53的值。X0.40.50.60.70.8lnx-0.916291-0.693147-0.510826-0.357765-0.223144解：由表可知x0=0.4,x1=0.5,x2=0.6,x3=0.7,x4=0.7，函数值：Y0=-0.916291,y1=-0.693147,y2=-0.510826,y3=-0.357765,y4=-0.223144建立一个主程序np.mclcclearnewpoly([0.4,0.5,0.6,0.7,0.8],[-0.916