数值分析实验报告 (2).doc

北京建筑工程学院理学院数值分析课程 实验报告课程名称数值分析 实验名称 实验地点 基础C-423 日期_姓名 班级 学号 指导教师 张健、王恒友 成 绩 【实验目的】1、掌握用MATLAB计算拉格朗日插值的方法，改变节点的数目，对插值结果进行初步分析。 2、 通过实例学习用插值方法解决实际问题。【实验要求】 1、按题目要求完成实验内容2、写出相应的Matlab程序3、给出实验结果4、对实验结果进行分析讨论5、写出相应的实验报告【实验内容】1、考虑函数在-5,5上取n+1个等距节点，构造Lagrange插值多项式（1） 分别画出n=2,4,6,8,10,12,14,16,18,20的Lagrange插值函数的图形和函数y=f(x)的图形（要求画在一张图上）（2） 分析观察的现象，大体说出Lagrange插值收敛的区间2、机翼加工问题已知机翼轮廓上的数据如表所示表：机翼轮廓数据表x/m035791112131415y/m01.21.72.02.12.01.81.21.01.6加工时需要x每改变0.1m时的y值，并画出相应的轮廓曲线，选用Newton插值方法。【实验步骤及结果】1. (1)n=2,4,6,8,10,12,14,16,18,20的Lagrange插值函数的图形和函数y=f(x)的图形(2) Lagrange插值收敛的区间为-1,1。一般认为差值多项式的次数越高，逼近f(x)的精度越高。实际上并非如此，由上图可以看出，在x=0附近当i值越大时Li(x)（i=2,4,6，20）与f(x)的近似程度比较好，但当x远离零点时，近似效果越差，以致完全失真。我们可以缩短差值区间，分段插值来减少误差。2. 机翼加工问题差商值的计算公式：fx1,x2,xn=i=0nf(xi)j=0,jin(xi-xj) 差商形式的差值公式fx=fx1+fx1,x2x-x1+fx1,x2,x3x-x1x-x2+fx1,x2,xnx-x1x-x2x-xn-1+fx1,x2,xn,xx-x1x-x2x-xn易知：Rx=fx1,x2,xn,xx-x1x-x2x-xn=0。用递归的方法编写newton差值公式：Newton(i,n)Ci+Newtoni+1,n ,in0 ,i=n+1通过matlab画出如下图所示结果。利用三次样条插值：【附录】1. （1）chazhi.m文件clear;clc;W=-5:1:5;%Xi(k),Yi(k)矩阵for i=1:10 n=2*i; for j=0:n X(i,j+1)=-5+10*(j/n); z=X(i,j+1); Y(i,j+1)=1./(1+z.2); end H(i,:)=L(W,X(i,:),Y(i,:); endplot(W,H(1,:),y-+,W,H(2,:),g-o,W,H(3,:),c-.,W,H(4,:),b-,W,H(6,:),y:+ ,W,H(7,:),g:o,W,H(8,:),c:.,W,H(9,:),b:);hold onplot(W,1./(1+W.2),r*),gtext(f1),gtext(f2),gtext(f3),gtext(f4),gtext(f5),gtext(f6),gtext(f7),gtext(f8),gtext(f9),gtext(f10),gtext(f11);clear i j n z;（2）Lagrange插值函数function y=L(x,xi,yi)n=0;m1=length(x);m2=length(xi);for i=1:m2 n=n+1; if xi(i)=5 break; endend for i=1:m1 s=0; for k=1:n t=1; for j=1:n if j=k t=t*(x(i)-xi(j)/(xi(k)-xi(j); end end s=s+t*yi(k); end y(i)=s;endend2. （1）Tmain.m文件clear;clc;W=0:0.1:15;m=length(W);x0=0 3 5 7 9 11 12 13 14 15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;n=length(x0);for i=1:m P(i)=Newton(x0,y0,W(i),1,n);endplot(W,P);clear i n m;（2）newton差值函数：function N=Newton(x0,y0,x,i,n)if i=n c=ChaShang(x0,y0,i); p(i)=c+Newton(x0,y0,x,i+1,n )*(x-x0(i);else p(i)=0;endN=p(i);end（3）差值函数：function sum=ChaShang(x,y0,n)sum=0;for i=1:n t=1; for j=1:n if j=i t=t*(x(i)-x(j); end end sum=sum+y0(i)/t;%f(x(i)/t;endsum;end三次样条插值代码：x0=0 3 5 7 9 11 12 13 14 15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;plot(x0