数值分析实验报告1

实验报告实验项目名称插值法实验室数学实验室所属课程名称数值逼近实验类型算法设计实验日期班级学号姓名成绩实验概述：【实验目的及要求】本次实验的目的是熟练《数值分析》第二章“插值法”的相关内容，掌握三种插值方法：牛顿多项式插值，三次样条插值，拉格朗日插值，并比较三种插值方法的优劣。本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并在MATLAB软件中去实现。【实验原理】《数值分析》第二章“插值法”的相关内容，包括：牛顿多项式插值，三次样条插值，拉格朗日插值的相应算法和相关性质。【实验环境】（使用的软硬件）软件：MATLAB2012a硬件：电脑型号：联想Lenovo昭阳E46A笔记本电脑操作系统：Windows8专业版处理器：Intel（R）Core（TM）i3CPUM350@2.27GHz2.27GHz实验内容：【实验方案设计】第一步，将书上关于三种插值方法的内容转化成程序语言，用MATLAB实现；第二步，分别用牛顿多项式插值，三次样条插值，拉格朗日插值求解不同的问题。【实验过程】（实验步骤、记录、数据、分析）实验的主要步骤是：首先分析问题，根据分析设计MATLAB程序，利用程序算出问题答案，分析所得答案结果，再得出最后结论。实验一：已知函数在下列各点的值为xi0.20.40.6.0.81.0f（xi）0.980.920.810.640.38试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值。用图给出{（xi，yi），xi=0.2+0.08i，i=0，1,11,10}，P4（x）及S（x）。(1)首先我们先求牛顿插值多项式，此处要用4次牛顿插值多项式处理数据。已知n次牛顿插值多项式如下：Pn=f(x0)+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+···+f[x0,x1，···xn](x-x0)···(x-xn-1)我们要知道牛顿插值多项式的系数，即均差表中得部分均差。在MATLAB的Editor中输入程序代码，计算牛顿插值中多项式系数的程序如下：functionvarargout=newtonliu(varargin)clear,clcx=[0.20.40.60.81.0];fx=[0.980.920.810.640.38];newtonchzh(x,fx);functionnewtonchzh(x,fx)%由此函数可得差分表n=length(x);fprintf('*****************差分表*****************************\n');FF=ones(n,n);FF(:,1)=fx';fori=2:nforj=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1));endendfori=1:nfprintf('%4.2f',x(i));forj=1:ifprintf('%10.5f',FF(i,j));endfprintf('\n');end由MATLAB计算得：xif(xi)一阶差商二阶差商三阶差商四阶差商0.200.9800.400.920-0.300000.600.810-0.55000-0.625000.800.640-0.85000-0.75000-0.208331.000.380-1.30000-1.12500-0.62500-0.52083所以有四次插值牛顿多项式为：P4（x）=0.98-0.3(x-0.2)-0.62500(x-0.2)(x-0.4)-0.20833(x-0.2)(x-0.4)(x-0.6)-0.52083(x-0.2)(x-0.4)(x-0.6)(x-0.8）(2)接下来我们求三次样条插值函数。用三次样条插值函数由上题分析知,要求各点的M值：06.7500-4.5000-3.7500-0MMMMM2.5000000020.50000000.50002.500000.50002.5000000243210三次样条插值函数计算的程序如下：functiontgsanci(n,s,t)%n代表元素数，s,t代表端点的一阶导。x=[0.20.40.60.81.0];y=[0.980.920.810.640.38];n=5forj=1:1:n-1h(j)=x(j+1)-x(j);endforj=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endforj=1:1:n-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=0d(n)=0a=zeros(n,n);forj=1:1:na(j,j)=2;endr(1)=0;u(n)=0;forj=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*d'p=zeros(n-1,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);endp得到m=(0-1.6071-1.0714-3.10710)T即M0=0;M1=-1.6071;M2=-1.0714;M3=-3.1071;M4=0则根据三次样条函数定义，可得：S(x)=]0.1,8.0[x)8.0(9.10.13036.38.00-x0.12.5893-]8.0,6.0[x6.0x3036.3x8.00857.4.60-x5893.2-x.800.8929-]6.0,4.0[x4.0x7508.4x6.04.6536.40-x8929.0x.601.3393-]4.0,2.0[x)2.0(6536.44.0900.42.01.33934.003333333，）（）（）（），（）（）（）（），（）（）（）（，）（）（）（xxxxxxx接着，在CommandWindow里输入画图的程序代码，下面是画牛顿插值以及三次样条插值图形的程序：x=[0.20.40.60.81.0];y=[0.980.920.810.640.38];plot(x,y)holdonfori=1:1:5y(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=[011011]x0=0.2+0.08*kfori=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endplot(x0,y0,'o',x0,y0)holdony1=spline(x,y,x0)plot(x0,y1,'o')holdons=csape(x,y,'variational')fnplt(s,'r')holdongtext('三次样条自然边界')gtext('原图像')gtext('4次牛顿插值')运行上述程序可知：给出的{（xi，yi），xi=0.2+0.08i，i=0，1,11,10}点，S（x）及P4（x）图形如下所示：实验二：在区间[-1,1]上分别取10,20n用两组等距节点对龙格函数21()125fxx作多项式插值及三次样条插0.20.30.40.50.60.70.80.911.11.20.20.30.40.50.60.70.80.91三次样条自然边界原图像4次牛顿插值值，对每个n值，分别画出插值函数即()fx的图形。我们先求多项式插值：在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令（如牛顿插值公式）：function[C,D]=newpoly(X,Y)n=length(X);D=zeros(n,n)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)))m=length(C);C(m)=C(m)+D(k,k);end当n=10时，我们在CommandWindow中输入以下的命令：clear,clf,holdon;X=-1:0.2:1;Y=1./(1+25*X.^2);[C,D]=newpoly(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,'.');gridon;xp=-1:0.2:1;z=1./(1+25*xp.^2);plot(xp,z,'r')得到插值函数和f（x）图形：当n=20时，我们在CommandWindow中输入以下的命令：clear,clf,holdon;X=-1:0.1:1;Y=1./(1+25*X.^2);[C,D]=newpoly(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,'.');gridon;xp=-1:0.1:1;z=1./(1+25*xp.^2);plot(xp,z,'r')得到插值函数和f（x）图形：下面再求三次样条插值函数，在MATLAB的Editor中建立一个多项式的M-file,输入下列程序代码：functionS=csfit(X,Y,dx0,dxn)N=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N));C=H(2:N);U=6*diff(D);B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1));B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N));fork=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=U(N-1)/B(N-1);fork=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;fork=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;S(k+1,4)=Y(k+1);end当n=10时，我们在CommandWindow中输入以下的命令：clear,clcX=-1:0.2:1;Y=1./(25*X.^2+1);dx0=0.0739644970414201;dxn=-0.0739644970414201;S=csfit(X,Y,dx0,dxn)x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1));x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2));x3=0:0.01:0.5;y3=polyval(S