1、实 验 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:mz=x(i);s=0.0;for k=1:np=1.0;f
2、or j=1:nif j=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(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
3、=y_lagr1(x0,y0,x);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,spline); 5 1 5 2 3 4 5 6 8 9 1 2). 0 1 2 3 4 5 7 80 0 0 0 0 0 0 1 2 3 4 5 6 7 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),运行后,得到各节点和插值点的值如下:X y y1 y2 y3 00.06670.13330.200
4、0 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
5、.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,产生一组数据(x i,yi,i=1,2,n) ,再在 yi 上添加随机干扰(可用 rand 产生(0,1)均匀分布随机
6、数,或用 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.6679 1 2 4 5 6 7 8 0 0 0 0 m=3 ,y_2_3f
7、= 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) 与 成指数变化关系,所以在0()exp/)vtVt通过曲线拟合的时候,使用指数曲线 ya 1ea2x。(注:这是一种非线性拟合). 首先进行变量代换(
8、为了程序运行的方便,在程序中用 v1代替 v(t),v 2是拟合后的曲线方程):对 变 形 后 取 对 数 , 有0(1()exp(/)vtvt 0ln(0)ln()(/tt。令 y=ln(10- ) ,a1=ln(10- ) ,a2= -1/,则(tv0,= -1/ a20exp(a(拟合图形与节点已显示在下图中) 。程序如下: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/a 1v0=10-exp(a2)v2=1
9、0-(10-v0)*exp(-t/);plot(t,v1,k+,t,v2,r)结果: = 3.5269v0 = 5.6221第 8 题分析:1)由已知点的排列可见,选前五个点作直线拟合比较合理。为了达到拟合为正比例函数的目的,在 x 和 F 的数组中增加一组(0,0) ,尽管如此,程序运行后仍然有 a2 。02)后五个点的曲线可用二次曲线拟合。程序如下: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=po
10、lyfit(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) K T( C) P(103KN/m2) K68 97981 00848 106 97918 0069668 13324 00897 106 14277 00
11、75387 90078 00762 140 96563 0061187 13355 00807 140 12463 00651试求 T=99(C)和 P=10.3(10 3KN/m2)下的 K。解: 首先,找出温度 T相等时,导热系数 K与压力 P的关系。由于在每一温度时,仅给出两个 K、P 值,因此采用线性近似,把 K、P 看作是线性关系。要求的是温度为 99时的导热系数,99介于 87和 106之间,所以先分别求出 P=10.3(10 3KN/m2),T=87和T=106时的导热系数 K。运行如下程序:p1=9.7981,13.324; k1=0.0848,0.0897; %T=68p2=
12、9.0078,13.355; k2=0.0762,0.0807; %T=87p3=9.7918,14.277; k3=0.0696,0.0753; %T=106p4=9.6563,12.463; k4=0.0611,0.0651; %T=140a2=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)xl
13、abel (丙烷压力 P(*103KN/m2)ylabel (丙烷导热系数 K)title(在不同温度下丙烷导热系数与压力的关系图)gtext(T=68),gtext(T=87),gtext(T=106),gtext(T=140)运行后得到下图:图中所标点即 P=10.3(10 3KN/m2)时,T=87和 T=106对应的导热系数K值。在 T=87和 T=106之间,仍采用线性近似来求T=99时的导热系数 K。运行如下程序:x=87,106;y=x1,x2;a=polyfit(x,y,1);z=polyval(a,99)z=0.0729plot(99,z,k+,x,y)gridxlabel(丙烷温度 T()ylabel(丙烷导热系数 K)title(压力 10.3*103KN/m2时丙烷导热系数与温度的关系图)结论:T=99()和 P=10.3(10 3KN/m2)下的 K=0.0729。?疑问:能否用二维分段插值函数 interp2 来做?怎么看不懂系统的出错提示? 3 4 5 6 8 ) 0