插值与拟合综合作业.doc

1.画出由参数方程确定的空间曲线图形
（1）新建一个M-file
（2）输入命令：

（3） 结果为：

2.绘出旋转抛物面的图形
（1）新建一个M-file
（2）输入命令：

（3） 结果为：

3.从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24．试估计每隔1/10小时的温度值．
（1）新建一个M-file
（2）输入命令：

（3） 结果为：

4.已知飞机轮廓线的数据如下表所列，求每改变0.1时的的值。
表1 轮廓数据线
x
0
3
5

7
9
11
12
13
14
15
y
0
1.2
1.7
2.0
2.1
2.0
1.8
1.2
1.0
1.6
1） 新建一个M-file
2） 输入命令：

3） 结果为：

最后，在Matlab中输入要求的值，比如：第五个值（也就是点x=0.5）就输入，yi(5)
结果为ans =

0.1915

5. 例：测得平板表面3×5网格点处的温度分别为：
82
81
80
82
84
79
63
61
65
81
84
84
82
85
86
试作出平板表面的温度分布曲面z=f(x,y)的图形．
1.先在三维坐标画出原始数据，画出粗糙的温度分布曲线图.
输入以下命令：
x=1:5;
y=1:3;
temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];
mesh(x,y,temps)

再输入以下命令:
xi=1:0.01:5;
yi=1:0.01:3;
zi=interp2(x,y,temps,xi',yi,'cubic');
surf(xi,yi,zi)
shading flat
画出插值后的温度分布曲面图.

6.要在一山区修建公路，首先测得一些地点的高程，数据如表所列：
4800
1350
1370
1390
1400
1410
960
940
880
800
690
570
430
290
210
150
4400
1370
1390
1410
1430
1440
1140
1110
1050
950
820
690
540
380
300
210
4000
1380
1410
1430
1450
1470
1320
1280
1200
1080
940
780
620
460
370
350
3600
1420
1430
1450
1480
1500
1550
1510
1430
1300
1200
980
850
750
550
500
3200
1430
1450
1460
1500
1550
1600
1550
1600
1600
1600
1550
1500
1500
1550
1500
2800
950
1190
1370
1500
1200
1100
1550
1600
1550
1380
1070
900
1050
1150
1200
2400
910
1090
1270
1500
1200
1100
1350
1450
1200
1150
1010
880
1000
1050
1100
2000
880
1060
1230
1390
1500
1500
1400
900
1100
1060
950
870
900
930
950
1600
830
980
1180
1320
1450
1420
1400
1300
700
900
850
840
380
780
750
1200
740
880
1080
1130
1250
1280
1230
1040
900
500
700
780
750
650
550
800
650
760
880
970
1020
1050
1020
830
800
700
300
500
550
480
350
400
510
620
730
800
850
870
850
780
720

650
500
200
300
350
320
0
370
470
550
600
670
690
670
620
580
450
400
300
100
150
250
y/x
0
400
800
1200
1600
2000
2400
2800
3200
3600
4000
4400
4800
5200
5600
试利用表中数据，绘制这一山区的地貌网格图、等高线图，并用几种插值方法进行比较

输入：
x=0:400:5600; %x轴坐标的值
y=4800:-400:0; %y轴坐标的值

%%——————————录入每行（x，y）的高程——————————————————
x1=[1350,1370,1390,1400,1410,960,940,880,800,690,570,430,290,210,150];
x2=[1370,1390,1410,1430,1440,1140,1110,1050,950,820,690,540,380,300,210];
x3=[1380,1410,1430,1450,1470,1320,1280,1200,1080,940,780,620,460,370,350];
x4=[1420,1430,1450,1480,1500,1550,1510,1430,1300,1200,980,850,750,550,500];
x5=[1430,1450,1460,1500,1550,1600,1550,1600,1600,1600,1550,1500,1500,1550,1500];
x6=[950,1190,1370,1500,1200,1100,1550,1600,1550,1380,1070,900,1050,1150,1200];
x7=[910,1090,1270,1500,1200,1100,1350,1450,1200,1150,1010,880,1000,1050,1100];
x8[880,1060,1230,1390,1500,1500,1400,900,1100,1060,950,870,900,930,950];
x9=[830,980,1180,1320,1450,1420,1400,1300,700,900,850,840,380,780,750];
x10=[740,880,1080,1130,1250,1280,1230,1040,900,500,700,780,750,650,550];
x11=[650,760,880,970,1020,1050,1020,830,800,700,300,500,550,480,350];
x12=[510,620,730,800,850,870,850,780,720,650,500,200,300,350,320];
x13=[370,470,550,600,670,690,670,620,580,450,400,300,100,150,250];
Heigh=[x1;x2;x3;x4;x5;x6;x7;x8;x9;x10;x11;x12;x13];
surf(x,y,Heigh),title('地貌网格图') %地貌网格图

figure(2); %新建一个图形窗口
[X,Y]=meshgrid(x,y);
subplot(1,2,1)
[c,h]=contour(X,Y,Heigh,20);clabel(c,h),title('平面等高线') %画出平面等高线
subplot(1,2,2)
contour3(x,y,Heigh,800),title('地貌图') %绘制地貌图

xi=linspace(0,5600,50);yi=linspace(0,4800,50); %给出新的插值坐标
[XI,YI]=meshgrid(xi,yi);

%%-------------比较四种插值后的表面图——-----------------------------
figure(3); %新建一个图形窗口
%——————对数据（xi,yi,hi）使用三次插值————————————----
HI1=interp2(X,Y,Heigh,XI,YI,'cubic'); %对数据（xi,yi,hi）使用三次多项式在网格{X,Y}上插值
subplot(2,2,1),surf(XI,YI,HI1); %画出三次多项式插值后的地貌图
title('三次多项式插值后的地貌图')
shading flat
%——————对数据（xi,yi,hi）使用三次样条插值————————————
HI2=interp2(X,Y,Heigh,XI,YI,'spline'); %对数据（xi,yi,hi）使用三次样条在网格{X,Y}上插值
subplot(2,2,2),surf(XI,YI,HI2); %画出三次样条插值后的地貌图
title('三次样条插值后的地貌图')
shading flat
%——————对数据（xi,yi,hi）使用最近邻点插值————————————
HI3=interp2(X,Y,Heigh,XI,YI,'nearest'); %对数据（xi,yi,hi）使用最近邻点在网格{X,Y}上插值
subplot(2,2,3),surf(XI,YI,HI3); %画出最近邻点插值后的地貌图
title('最近邻点插值后的地貌图')
shading flat
%——————对数据（xi,yi,hi）使用线性插值————————————----
HI4=interp2(X,Y,Heigh,XI,YI,'linear'); %对数据（xi,yi,hi）使用线性在网格{X,Y}上插值
subplot(2,2,4),surf(XI,YI,HI4); %画出线性插值后的地貌图
title('线性插值后的地貌图')
shading flat

%%-------------比较四种插值后的等高线图——-------------------------
figure(4); %新建一个图形窗口
subplot(2,2,1),contour(XI,YI,HI1,10,'r');title('三次多项式插值后的平面等高线')
subplot(2,2,2),contour(XI,YI,HI2,10,'r');title('三次样条插值后的平面等高线')
subplot(2,2,3),contour(XI,YI,HI3,10,'r');title('最近邻点插值后的平面等高线')
subplot(2,2,4),contour(XI,YI,HI4,10,'r');title('线性插值后的平面等高线')
输出：

7.2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功．整个试验期为20多天，小浪底从6月19日开始预泄放水，直到7月13日恢复正常供水结束．小浪底水利工程按设计拦沙量为75.5亿立方米，在这之前，小浪底共积泥沙达14.15亿吨．这次调水调试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙．在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2700立方米/每秒，使小浪底水库的排沙量也不断地增加．下面是由小浪底观测站从6月29日到7月10日检测到的试验数据：
表1: 试验观测数据 单位：水流为立方米 / 秒，含沙量为公斤 / 立方米
日期
6.29
6.30
7.1
7.2
7.3
7.4
时间
8:00
20:00
8:00
20:00
8:00
20:00
8:00
20:00
8:00
20:00
8:00
20:00
水流量
1800
1900
2100
2200

2300
2400
2500
2600
2650
2700
2720
2650
含沙量
32
60
75
85
90
98
100
102
108
112
115
116
日期
7.5
7.6
7.7
7.8
7.9
7.10
时间
8:00
20:00
8:00
20:00
8:00
20:00
8:00
20:00
8:00
20:00
8:00
20:00
水流量
2600
2500
2300
2200
2000
1850
1820
1800
1750
1500
1000
900
含沙量
118
120
118
105
80
60
50
30
26
20
8
5
现在，根据试验数据建立数学模型研究下面的问题：
(1) 给出估算任意时刻的排沙量及总排沙量的方法；
(2) 确定排沙量与水流量的变化关系。

(1) A.某一时刻的排沙量，其中为时刻的水流量，而为时刻的含沙量。因为排沙量与时间的散点图基本符合正态曲线，如图一所示。所以，排沙量的对数与时间的函数关系就应该符合二次函数关系（曲线见图二），因而排沙量取对数后，在与时间进行二次回归。
B.假设,两边取对数后先由表二做出排沙量的自然对数与时间的散点图见图一，并用Matlab软件进行拟合，得到排沙量的自然对数与时间的回归方程为：

所以，任意时刻排沙量与时间之间的关系为：

所以，任意时刻总排沙量与时间之间的关系为：

在M文件中输入：
t=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24]; %时间
v(t)=[1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720,2650,2600,2500,2300,2200,2000,1850,1820,1800,1750,1500,1000,900,]; %水流量
S(t)=[32,60,75,85,90,98,100,102,108,112,115,116,118,120,118,105,80,60,50,30,26,20,8,5]; %含沙量
y=v(t).*S(t); %排沙量
s=log(y); %排沙量的对数
plot(t,y,'*'),title('图一');
figure(2)
plot(t,s,'*')，title（'图二');
p=polyfit(t,s,2) %对(x,y)做二阶线性拟合
z=polyval(p,t);
figure（3）
plot(t,s,'r*',t,z,'b')
syms t
V=int(exp(-0.0209*t^2+0.4289*t+10.6312),t)

输出结果为：

p =

-0.0209 0.4289 10.6312

V =

50/209*pi^(1/2)*exp(107272353/8360000)*209^(1/2)erf(1/100*209^(1/2)*t-4289/41800*209^(1/2))

(2) 确定排沙量与水流量的变化关系，从实验数据可以看出，开始排沙量随水流量的增加而增加，而后是随水流量的减少而减少的。可以得知：变化规律不是呈线性关系，所以，把这个问题分两部分（第一部分是从开始水流量增加到最大值，第二部分是从水流量的最大值到结束）分别进行研究，从而得出水流量与排沙量的关系。具体数据如下表：

第一部分实验观测数据
单位：水流量为 ，含沙量为
序号
1
2
3
4
5
6
7
8
9
10
11
水流量
1800
1900
2100

2200
2300
2400
2500
2600
2650
2700
2720
含沙量
32
60
75
85
90
98
100
102
108
112
115

第二部分实验观测数据
单位：水流量为 ，含沙量为
序号
1
2
3
4
5
6
7
8
9
10
11
12
13
水流量
2650
2600
2500
2300
2200
2000
1850
1820
1800
1750
1500
1000
900
含沙量
116
118
120
118
105
80
60
50
30
26
20
8
5

对于第一部分，我们用Matlab作图可以看出其变化趋势，同时，我们用多项式作最小拟合。
程序如下：
x=[1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720];
y=[32,60,75,85,90,98,100,102,108,112,115];
plot(x,y,'r')
如图所示：

利用已知数据对其三次多项式拟合，编写Matlab命令如下：
x=[1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720];
y=[32,60,75,85,90,98,100,102,108,112,115];
A=polyfit(x,y,3)
z=polyval(A,x);
plot(x,y,'*',x,y,'r',x,z,'b')

结果为：A =

1.0e+003 *

0.0000 -0.0000 0.0032 -2.4929
所以，可得拟合多项式为：

利用已知数据对其四次多项式拟合，编写Matlab命令如下：
x=[1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720];
y=[32,60,75,85,90,98,100,102,108,112,115];
A=polyfit(x,y,4)
z=polyval(A,x);
plot(x,y,'*',x,y,'r',x,z,'b')

结果为：A =

1.0e+003 *

-0.0000 0.0000 -0.0000 0.0121 -7.4347

所以，可得拟合多项式为：

从上面的三次多项式拟合和四次多项式拟合效果来看，差别不是很大，基本可以看出排沙量与水流量的关系。

对于第二部分，可以用类似的处理。用最小二乘法作三次和四次多项式拟合，拟合效果与Matlab程序如下：
第二部分三次多项式拟合效果图：

x=[2650,2600,2500,2300,2200,2000,1850,1820,1800,1750,1500,1000,900];
y=[116,118,120,118,105,80,60,50,30,26,20,8,5];
p=polyfit(x,y,3)
z=polyval(p,x);
plot(x,y,'*',x,y,'r',x,z,'b')

结果为：p =

-0.0000 0.0006 -1.0075 497.8489

所以，可得拟合多项式为：

三次多项式拟合效果图：

x=[2650,2600,2500,2300,2200,2000,1850,1820,1800,1750,1500,1000,900];
y=[116,118,120,118,105,80,60,50,30,26,20,8,5];
p=polyfit(x,y,4)
z=polyval(p,x);
plot(x,y,'*',x,y,'r',x,z,'b')

结果为：p =

-0.0000 0.0000 -0.0016 1.4334 -468.8030
所以，可得拟合多项式为：

通过比较可知：四次多项式拟合更合适
关系：开始排沙量随水流量的增加而增加，而后是随水流量的减少而减少的