插值法及程序.doc

数值计算方法第2章大作业 题目： 插值逼近 姓名： 冯 婷 婷 学号： 104753140674 专业：运筹学与控制论 插值逼近一、插值问题设是区间上的一个实函数,是上个互异实数,已知在的值, 求一个次数不超过的多项式使其满足这就是多项式插值问题。其中称为的次插值多项式, 称为被插函数, 称为插值节点, 称为插值点,称为插值区间, 式称为插值条件。 从几何意义来看,上述问题就是要求一条多项式曲线, 使它通过已知的个点,并用近似表示.即其中为实数，就称为插值多项式，相应的插值法称为多项式插值。二、拉格朗日插值法在求满足插值条件次插值多项式之前，先考虑一个简单的插值问题：对节点中任一点，作一n次多项式，使它在该点上取值为1，而在其余点上取值为零，即 上式表明个点都是次多项式的零点，故可设 其中，为待定系数。由条件立即可得 故 由上式可以写出个次插值多项式。我们称它们为在个节点上的次基本插值多项式或次插值基函数。三、牛顿插值法任何一个不高于次多项式，都可以表示成函数的线性组合。既可以把满足插值条件的次插值多项式写成如下形式其中，为待定系数。这种形式的插值多项式称为牛顿插值多项式，记为，即 因此，牛顿插值多项式是插值多项式的另一种表示形式。设函数在等距节点处的函数值为已知，其中是正常数，称步长。我们称两个相邻点和处函数之差为函数在点处以为步长的一阶向前差分，记作，即于是，函数在各节点处的一阶差分依次为又称一阶差分的差分为二阶差分。 一般的，定义函数在点处的阶差分为。在等距节点情况下，可以利用差分表示牛顿插值多项式的系数。 四、三次样条插值定义：给定区间的一个划分,如果函数满足： (1); (2)在每个小区间上是次数不超过的多项式; (3)在每个内节点上具有二阶连续导数,则称为关于上述划分的一个三次多项式样条函数，简称三次样条。在每个小区间上是一个次数不超过的多项式, 因此需确定四个待定常数, 一共有个小区间,故应确定个系数, 在个内节点上具有二阶连续导数，应满足条件即有个连续条件，再加上 满足的插值条件个，共计个，因此还需要个条件才能确定，通常可在区间的端点上各加一个条件（称为边界条件）。可根据实际问题的要求给定，常见有以下三种： (1)已知两端的一阶导数值，即。 (2)两端的二阶导数已知，即。(3)当是以为周期函数时，则要求也是周期函数。这是边界条件应满足五、例题展示分别用拉格朗日插值法，牛顿插值法和三次样条插值法算法计算例 已知函数表如下：0.10.20.30.40.50.60.099830.198670.295520.389420.479430.56464计算的值。1、利用拉格朗日插值法计算过程如下： （1）利用线性插值法求近似值因为0.12位于0.1与0.2之间，故取节点（2）利用抛物插值所求的近似值为因为0.12位于0.1，0.2与0.3之间，故取节点%建立M文件function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s;end%在命令行窗口输入x=0.1:0.1:0.6;y=sin(x);x1=0.1 0.2;y1=sin(x1);z1=lagrange(x1,y1,x);x2=0.1 0.2 0.3;y2=sin(x2);z2=lagrange(x2,y2,x);a=0.12;b1=lagrange(x1,y1,a)b2=lagrange(x2,y2,a)plot(x,y,+k,x,z1,-g,x,z2,-r);xlabel(X)ylabel(Y)title(lagrange插值)%输出结果b1 = 0.1196b2 = 0.11982、利用牛顿插值法计算过程如下：构造差分表如下：0.10.20.30.40.099830.198670.295520.389420.098840.096850.09390-0.00199-0.00295-0.00096（1）利用线性插值所求的近似值为（1）利用线性插值法求近似值因为0.12位于0.1与0.2之间，故取节点（2）利用抛物插值所求的近似值为因为0.12位于0.1，0.2与0.3之间，故取节点%建立M文件function y=newton(x0,y0,x)n=length(x0);m=length(x);y1=zeros(m);for j=1:mfor i=n:(-1):2y1(j)=(x(j)-x0(i-1)*(cs(x0(1:i),y0(1:i)+y1(j);endy(j)=y1(j)+y0(1);endEnd%在命令行窗口输入x=0.1:0.1:0.6;y=sin(x);x1=0.1 0.2;y1=sin(x1);z1=newton(x1,y1,x);x2=0.1 0.2 0.3;y2=sin(x2);z2=newton(x2,y2,x);a=0.12;b1=newton(x1,y1,a)b2=newton(x2,y2,a)plot(x,y,+k,x,z1,-g,x,z2,-r);xlabel(X)ylabel(Y)title(Newton插值)%输出结果b1 = 0.1196b2 = 0.11983、利用三次样条插值法计算过程如下： 0.10.20.30.40.099830.198670.295520.389421.000000.99998%建立M文件function y,M=scyt(X,Y,D0,Dn,x) n=length(X);m=length(x);S=-diag(ones(n,1)+diag(ones(n-1,1),-1);h=X*S;H=h(1:n-1);L=zeros(1,n-2);for i=1:n-2L(i)=H(i+1)/(H(i)+H(i+1);endmu=ones(1,n-2)-L;L=1,L;mu1=mu,1;A=2*diag(ones(n,1)+diag(mu1,-1)+diag(L,1);d0=6/H(1)*(Y(2)-Y(1)/H(1)-D0);dn=6/H(n-1)*(Dn-(Y(n)-Y(n-1)/H(n-1);d=zeros(1,n-2);for j=2:n-1d(j-1)=6*cs(X(j-1:j+1),Y(j-1:j+1);end d1=d0,d,dn;M=inv(A)*d1;y=zeros(1,m);for i=1:mif(x(i)=X(n)for j=1:n-1if(x(i)=X(j)&x(i)X(j+1) y(i)=(X(j+1)-x(i)3*M(j)+(-X(j)+x(i)3*M(j+1)/6/H(j)+(X(j+1)-x(i)*Y(j)+(x(i)-X(j)*Y(j+1)/H(j)-H(j)/6*(X(j+1)-x(i)*M(j)+(-X(j)+x(i)*M(j+1); endendelsey(i)=Y(n);endendEnd%在命令行窗口输入x0=0.1:0.1:0.4;y0=0.099830.198670.295520.389