清华大学_计算方法(数学实验)实验2插值与拟合.doc

实 验 2 插 值 与 拟 合 系 班 姓名 学号 【实验目的】1、 掌握用MATLAB计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点的数目，对三种插值结果进行初步分析。2、 掌握用MATLAB作线性最小二乘的方法。3、 通过实例学习如何用插值方法与拟合方法解决实际问题，注意二者的联系和区别。【实验内容】预备：编制计算拉格朗日插值的M文件：以下是拉格朗日插值的名为y_lagrl的M文件：function y=y_lagr1(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s; end第1题（d）选择函数y=exp(-x2) (-2x2),在n个节点上（n不要太大，如511）用拉格朗日、分段线性、三次样条三种插值方法，计算m个插值点的函数值（m要适中，如50100）。通过数值和图形输出，将三种插值结果与精确值进行比较。适当增加n，在作比较，由此作初步分析。 运行如下程序：n=7;m=61;x=-2:4/(m-1):2;y=exp(-x.2);z=0*x;x0=-2:4/(n-1):2;y0=exp(-x0.2);y1=y_lagr1(x0,y0,x);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,spline);xyy1y2y3 plot(x,z,w,x,y,r-,x,y1,b:,x,y2,m,x,y3,b) gtext(y=exp(-x2),gtext(Lagr.),gtext(Piece.-linear.),gtext(Spline),运行后，得到各节点和插值点的值如下：Xy y1 y2 y3 00.06670.13330.2000 0.86670.93331.0000 1.60001.6667 1.93332.00001.0000 1.0000 1.0000 1.00000.9956 0.9958 0.9641 0.99470.9824 0.9831 0.9282 0.97970.9608 0.9624 0.8924 0.9561 0.4718 0.4626 0.4995 0.48360.4185 0.4062 0.4523 0.43290.3679 0.3531 0.4051 0.3837 0.0773 0.1292 0.1087 0.05340.0622 0.1271 0.0937 0.0347 0.0238 0.0685 0.0334 0.01040.0183 0.0183 0.0183 0.0183将三种插值结果y1,y2,y3与精确值y项比较，显然y1在节点处不光滑,拉格朗日插值出现较大的振荡,样条插值得结果是最好的.增加n值（使n=11），再运行以上程序，得到的图形如右图所示，比较这两个图可发现，节点增加后，三种插值方法结果的准确度均有所提高，因此可近似地认为：增加节点个数可以提高插值结果的准确程度。第3题用给定的多项式，如y=x3-6x2+5x-3，产生一组数据（xi,yi,i=1,2,n），再在yi上添加随机干扰（可用rand产生（0，1）均匀分布随机数，或用randn产生N（0，1）分布随机数），然后用xi和添加了随机干扰的yi作3次多项式拟合，与原系数比较。如果作2或4次多项式拟合，结果如何？解：2编制y_2_3.m文件n=15;x=0:8/(n-1):8;y=x.3-6*x.2+5*x-3;z=0*x;y0=y+rand(1,15);f=polyfit(x,y0,m);r=polyval(f,x)pl2ot(x,z,k,x,y,r:,x,r,b) 程序及运行结果如下： m=2 ,y_2_3f = 5.9888 -31.9916 17.6679m=3 ,y_2_3f = 0.9850 -5.8326 4.5549 -2.4273m=4 ,y_2_3f = 0.0006 0.9950 -6.0302 5.2801 -2.8196运行后，比较拟合后多项式和原式的系数，发现三次多项式系数与原系数比较接近，四次多项式的四次项系数很小。作图后，发现二次多项式与原函数的差别比较大，而三次多项式和四次多项式符合得比较好。第6题解: 分析:电容器充电的数学模型已经建立。由题目给的方程(已知V=10)可见，(t)与成指数变化关系，所以在通过曲线拟合的时候，使用指数曲线ya1ea2x。(注：这是一种非线性拟合). 首先进行变量代换(为了程序运行的方便，在程序中用v1代替v(t)，v2是拟合后的曲线方程)：对变形后取对数，有。令y=ln(10-) ,a1=ln(10-) ,a2= -1/,则 ，= -1/ a2（拟合图形与节点已显示在下图中）。程序如下：t=0.5 1 2 3 4 5 7 9;v1=6.36 6.48 7.26 8.22 8.66 8.99 9.43 9.63;y=log(10- v1);a=polyfit(t,y,1);a1=a(1);a2=a(2);=-1/a1v0=10-exp(a2)v2=10-(10-v0)*exp(-t/);plot(t,v1,k+,t,v2,r)结果： = 3.5269v0 = 5.6221第8题分析：1）由已知点的排列可见，选前五个点作直线拟合比较合理。为了达到拟合为正比例函数的目的，在x和F的数组中增加一组（0，0），尽管如此，程序运行后仍然有a2。2）后五个点的曲线可用二次曲线拟合。程序如下：x1=0 1 2 4 7 9;F1=0 1.5 3.9 6.6 11.7 15.6;k=polyfit(x1,F1,1)y1=k(1).* x1;x2=9 12 13 15 17;F2=15.6 18.8 19.6 20.6 21.1;s=polyfit(x2,F2,2)y2=s(1).* x2.2+s(2).* x2+s(3);plot(x1,y1,g, x1, F1,k+, x2,y2,r, x2,F2,k+),axis(0 20 0 22)结果：k = 1.7085 0.0008s = -0.0764 2.6728 -2.2613第9题在化工生产中常常需要知道丙烷在各种温度T和压力P下的导热系数K。下面是实验得到的一组数据：T（C）P(103KN/m2)KT（C）P(103KN/m2)K6897981008481069791800696681332400897106142770075387900780076214096563006118713355008071401246300651试求T=99（C）和P=10.3（103KN/m2）下的K。解： 首先，找出温度T相等时，导热系数K与压力P的关系。由于在每一温度时，仅给出两个K、P值，因此采用线性近似，把K、P看作是线性关系。要求的是温度为99时的导热系数，99介于87和106之间，所以先分别求出P=10.3（103KN/m2）,T=87和T=106时的导热系数K。运行如下程序： p1=9.7981,13.324; k1=0.0848,0.0897; %T=68 p2=9.0078,13.355; k2=0.0762,0.0807; %T=87 p3=9.7918,14.277; k3=0.0696,0.0753; %T=106 p4=9.6563,12.463; k4=0.0611,0.0651; %T=140 a2=polyfit(p2,k2,1); a3=polyfit(p3,k3,1); x1=polyval(a2,10.3); x2=polyval(a3,10.3); %x1,x2分别是P=10.3（103KN/m2）下 87和106时的k值 plot(10.3,x1,k+,10.3,x2,k+,p1,k1,p2,k2,p3,k3,p4,k4) xlabel (丙烷压力P(*103KN/m2) ylabel (丙烷导热系数K) title(在不同温度下丙烷导热系数与压力的关系图) gtext(T=68),gtext(T=87),gtext(T=106),gtext(T=140)运行后得到下图：图中所标点即P=10.3（103KN/m2）时,T=87和T=106对应的导热系数K值。 在T=87和T=106之间，仍采用线性近似来求T=99时的导热系数K。运行如下程序：x=87,106;y=x1,x2;a=polyfit(x,y,