清华大学 - 计算方法（数学实验）实验2插值与拟合

实 验 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）

2

选择函数y=exp(-x) (-2≤x≤2),在n个节点上（n不要太大，如5~11）用拉格朗日、分段线性、三次样条三种插值方法，计算m个插值点的函数值（m要适中，如50~100）。通过数值和图形输出，将三种插值结果与精确值进行比较。适当增加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');

[x'y'y1'y2'y3']

plot(x,z,'w',x,y,'r--',x,y1,'b:',x,y2,'m',x,y3,'b') gtext('y=exp(-x^2)'),gtext('Lagr.'),gtext('Piece.-linear.'),gtext('Spline'),

运行后，得到各节点和插值点的值如下： X y y1 y2 y3 1y=exp(-x2)0 1.0000 1.0000 1.0000 1.0000 0.9Piece.-linear.0.0667 0.9956 0.9958 0.9641 0.9947 0.80.1333 0.9824 0.9831 0.9282 0.9797 0.70.2000 0.9608 0.9624 0.8924 0.9561 … … … … 0.60.8667 0.4718 0.4626 0.4995 0.4836 0.50.9333 0.4185 0.4062 0.4523 0.4329 0.41.0000 0.3679 0.3531 0.4051 0.3837 Lagr.… … … … 0.31.6000 0.0773 0.1292 0.1087 0.0534 0.21.6667 0.0622 0.1271 0.0937 0.0347 Spline0.1… … … … 1.9333 0.0238 0.0685 0.0334 0.0104 0-2-1.5-1-0.500.511.522.0000 0.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')

200150100500-50012345678180160140120100806040200 程序及运行结果如下： m=2 ,y_2_3

f = 5.9888 -31.9916 17.6679 m=3 ,y_2_3

-20012345678

f = 0.9850 -5.8326 4.5549 -2.4273 m=4 ,y_2_3

f = 0.0006 0.9950 -6.0302 5.2801 -2.8196

18016014012010080604020运行后，比较拟合后多项式和原式 的系数，发现三次多项式系数与原 系数比较接近，四次多项式的四次项系数很小。作图后，发现二次多项式与原函数的差别比较大，而三次多项式和四次多项式符合得比较好。

第6题

0-20012345678解: 分析:电容器充电的数学模型已经建立。由题目给的方程

v(t)?V?(V?V0)exp(?t/?)(已知V=10)可见，(t)与τ成指数变化关系，所以在通过曲线拟合的时候，使用指数曲线y＝a1ea2x。(注：这是一种非线性拟合). 首先进行变量代换(为了程序运行的方便，在程序中用v1代替v(t)，v2是拟合后的曲线方程)：

对v(t)?10?(10?v0)exp(?t/?)变形后取对数，有

ln(10?v(t))?ln(10?v0)?(?t/?)。

令y=ln(10-v(t)) ,a1=ln(10-?0) ,a2= -1/τ,则

v0?10?exp(a1) ，τ= -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/a1

v0=10-exp(a2)

v2=10-(10-v0)*exp(-t/τ); plot(t,v1,'k+',t,v2,'r')

结果：

τ = 3.5269 v0 = 5.6221

第8题

分析：1）由已知点的排列可见，选前五个点作直线拟合比较合理。为了达到拟合为正比例函数的目的，在x和F的数组中增加一组（0，0），尽管如此，程序运行后仍然有a2?0。

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.0008

s = -0.0764 2.6728 -2.2613

第9题

在化工生产中常常需要知道丙烷在各种温度T和压力P下的导热系数K。下面是实验得到的一组数据： 32P(103KN/m2) K K T（°C） T（°C） P(10KN/m) 68 106 9．7981 0．0848 9．7918 0．0696 68 106 13．324 0．0897 14．277 0．0753 87 87 9．0078 13．355 0．0762 0．0807 140 140 9．6563 12．463 0．0611 0．0651 试求T=99（°C）和P=10.3（103KN/m2）下的K。

解： 首先，找出温度T相等时，导热系数K与压力P的关系。由于在每一温度时，仅给出两个K、P值，因此采用线性近似，把K、P看作是线性关系。要求的是温度为99℃时的导

32

热系数，99℃介于87℃和106℃之间，所以先分别求出P=10.3（10KN/m）,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(*10^3KN/m^2)') ylabel ('丙烷导热系数K')

title('在不同温度下丙烷导热系数与压力的关系图')

gtext('T=68℃'),gtext('T=87℃'),gtext('T=106℃'),gtext('T=140℃')

32

运行后得到下图：图中所标点即P=10.3（10KN/m）时,T=87℃和T=106℃对应的导热系数K值。 ?ú2?í????è??±?íéμ?èè?μêyó??1á|μ?1??μí?0.09 T=68? ? C 在T=87℃和T=106℃之间， 0.085仍采用线性近似来求 T=87??C0.08T=99℃时的导热系数K。 0.075T=106??C运行如下程序：

0.07

T=140??Cx=[87,106];y=[x1,x2]; 0.065a=polyfit(x,y,1);

0.069101112131415z=polyval(a,99) ±?íé?1á|P(*103KN/m2)z=0.0729 ?1á|10.3*10KN/mê±±?íéμ?èè?μêyó????èμ?1??μí?0.079plot(99,z,'k+',x,y) 0.078grid 0.0770.076xlabel('丙烷温度T(℃)')

0.075ylabel('丙烷导热系数K')

0.074title('压力10.3*10^3KN/m^2

0.073时丙烷导热系数与温度的关系图') 0.072 0.071μêyKéμ?èè?±?í32结论：T=99（℃）和P=10.3 （10KN/m）下的 K=0.0729。 3

2

μêyKéμ?èè?±?í0.0786889092949698100102104106±?íé???èT(??)

？疑问：能否用二维分段插值函数interp2来做？怎么看不懂系统的出错提示？

本文来源：https://www.bwwdw.com/article/089f.html