数值分析实验报告1.doc

《数值分析实验报告1.doc》由会员分享，可在线阅读，更多相关《数值分析实验报告1.doc（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、实验报告实验项目名称 插值法 实 验 室 数 学 实 验 室 所属课程名称 数值逼近 实 验 类 型 算法设计 实 验 日 期 班 级 学 号 姓 名 成 绩 实验概述：【实验目的及要求】本次实验的目的是熟练数值分析第二章“插值法”的相关内容，掌握三种插值方法：牛顿多项式插值，三次样条插值，拉格朗日插值，并比较三种插 值方法的优劣。本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并在 MATLAB软件中去实现。【实验原理】数值分析第二章“插值法”的相关内容，包括：牛顿多项式插值，三次样条插值，拉格朗日插值的相应算法和相关性质。【实验环境】（使用的软硬件）软件：MATLAB

2、 2012a硬件：电脑型号：联想 Lenovo 昭阳 E46A笔记本电脑操作系统：Windows 8 专业版 处理器：Intel（R）Core（TM）i3 CPU M 350 2.27GHz 2.27GHz实验内容：【实验方案设计】第一步，将书上关于三种插值方法的内容转化成程序语言，用 MATLAB 实现；第二步，分别用牛顿多项式插值，三次样条插值，拉格朗日插值求解不同的问题。【实验过程】（实验步骤、记录、数据、分析）实验的主要步骤是：首先分析问题，根据分析设计 MATLAB 程序，利用程序算出问题答案，分析所得答案结果，再得出最后结论。实验一：已知函数在下列各点的值为xi 0.2 0.4 0

3、.6 .0.8 1.0f（x i）0.98 0.92 0.81 0.64 0.38试用 4 次牛顿插值多项式 P4（x）及三次样条函数 S（x）（自然边界条件）对数据进行插值。用图给出（x i，y i），x i=0.2+0.08i，i=0，1, 11, 10，P 4（x）及 S（x）。(1)首先我们先求牛顿插值多项式，此处要用 4 次牛顿插值多项式处理数据。已知 n 次牛顿插值多项式如下：Pn=f(x0)+fx0,x1(x-x0)+ fx0,x1,x2(x-x0) (x-x1)+ fx0,x1，x n(x-x0) (x-xn-1)我们要知道牛顿插值多项式的系数，即均差表中得部分均差。在 MAT

4、LAB 的 Editor 中输入程序代码，计算牛顿插值中多项式系数的程序如下：function varargout=newtonliu(varargin)clear,clcx=0.2 0.4 0.6 0.8 1.0;fx=0.98 0.92 0.81 0.64 0.38;newtonchzh(x,fx);function newtonchzh(x,fx)%由此函数可得差分表n=length(x);fprintf(*差分表*n);FF=ones(n,n);FF(:,1)=fx;for i=2:nfor j=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i

5、+1);endendfor i=1:nfprintf(%4.2f,x(i);for j=1:ifprintf(%10.5f,FF(i,j);endfprintf(n);end由 MATLAB 计算得：xi f(xi) 一阶差商 二阶差商 三阶差商 四阶差商0.20 0.9800.40 0.920 -0.300000.60 0.810 -0.55000 -0.625000.80 0.640 -0.85000 -0.75000 -0.208331.00 0.380 -1.30000 -1.12500 -0.62500 -0.52083所以有四次插值牛顿多项式为：P4（x）=0.98-0.3(x-0

6、.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.75-4302.500.50.2.50 210三次样条插值函数计算的程序如下：function tgsanci(n,s,t) %n 代表元素数，s,t 代表端点的一阶导。x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;n=5for j=1:1:n-1h(j)=x(j

7、+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);endd(1)=0d(n)=0a=zeros(n,n);for j=1:1:na(j,j)=2;endr(1)=0;u(n)=0;for j=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*dp=zeros(n-1,

8、4); %p 矩阵为 S(x)函数的系数矩阵for j=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.1071 0)T 即 M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0则根据三次样条函数定义，可得：S(x)= 0.1,8x).0(91. 306.80-x.125893- .,6368574-5

9、9.20 4.75. . .6 ,2)(4333 ，）（）（）（ ） ，（）（）（）（ ） ，（）（）（）（ ，）（）（）（ x接着，在 Command Window 里输入画图的程序代码，下面是画牛顿插值以及三次样条插值图形的程序：x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;plot(x,y)hold on for i=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*(

10、x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=0 1 10 11x0=0.2+0.08*kfor i=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 )hold ony1=spline(x,y,x0)plot(x0,y1,o)hold ons=c

11、sape(x,y,variational)fnplt(s,r)hold ongtext(三次样条自然边界)gtext(原图像)gtext(4 次牛顿插值)运行上述程序可知：给出的（x i，y i），x i=0.2+0.08i，i=0，1, 11, 10点，S（x）及 P4（x）图形如下所示：0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.20.20.30.40.50.60.70.80.91三三三三三三三三三三三4三三三三三实验二：在区间-1,1上分别取 用两组等距节点对龙格函数 作多项式插值及三次样条插10,2n 21()5fx值，对每个 值，分别画出插值函数即

12、 的图形。()fx我们先求多项式插值：在 MATLAB 的 Editor 中建立一个多项式的 M-file,输入如下的命令（如牛顿插值公式）：function C,D=newpoly(X,Y)n=length(X);D=zeros(n,n) D(:,1)=Y for j=2:nfor k=j:nD(k,j)=(D(k,j-1)- D(k-1,j-1)/(X(k)-X(k-j+1);endendC=D(n,n);for k=(n-1):-1:1C=conv(C,poly(X(k) m=length(C);C(m)= C(m)+D(k,k);end当 n=10 时，我们在 Command Wind

13、ow 中输入以下的命令：clear,clf,hold on;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,.);grid on;xp=-1:0.2:1;z=1./(1+25*xp.2);plot(xp,z,r)得到插值函数和 f（x）图形：当 n=20 时，我们在 Command Window 中输入以下的命令：clear,clf,hold on;X=-1:0.1:1;Y=1./(1+25*X.2); C,D=newpoly(X,Y);x=-1:0.01:1;y=poly

14、val(C,x);plot(x,y,X,Y,.); grid on;xp=-1:0.1:1;z=1./(1+25*xp.2);plot(xp,z,r)得到插值函数和 f（x）图形：下面再求三次样条插值函数，在 MATLAB 的 Editor 中建立一个多项式的 M-file,输入下列程序代码：function S=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*(

15、D(1);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N); for k=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); end M(N)=U(N-1)/B(N-1);for k=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;for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1

16、);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 时，我们在 Command Window 中输入以下的命令：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(3,:),x3-X(3);x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,.)结果如图：