数学建模第五章.docx
《数学建模第五章.docx》由会员分享，可在线阅读，更多相关《数学建模第五章.docx（20页珍藏版）》请在冰豆网上搜索。
数学建模第五章
第五章作业
1,
（3）边界条件:
给定三次样条在左右端点的一阶导数
仿照例题5.1.4计算和绘图
程序：
x=[0,1,3,6,8,9];y=[-1,3,1,2,0,2,4,1];
pp=csape（x,y,'complete'）,pp.coefs
s=@（t,tj,c）c
（1）.*（t-tj）.^3+c
（2）.*（t-tj）.^2+c（3）.*（t-tj）+c（4）;
d1s=@（t,tj,c）3.*c
（1）.*（t-tj）.^2+2.*c
（2）.*（t-tj）+c（3）;
d2s=@（t,tj,c）6.*c
（1）.*（t-tj）+2.*c
（2）;
d3s=@（t,tj,c）6.*c
（1）.*ones（size（t））;
fork=1:
pp.pieces
c=pp.coefs（k,:
）;u=x（k）:
.01:
x（k+1）;v=s（u,x（k）,c）;
v1=d1s（u,x（k）,c）;v2=d2s（u,x（k）,c）;v3=d3s（u,x（k）,c）;
plot（u,v,'k',u,v1,'k:
',u,v2,'k-.',u,v3,'k--'）,holdon
end
legend（'三次样条（给定样条在区间左右端点的一阶导数）','样条的一阶导函数',...
'样条的二阶导函数','样条的三阶导函数'）
y1=[1,4,3,2,3,1];
plot（x,y1,'ko'）,holdoff
运行结果：
pp=
form:
'pp'
breaks:
[013689]
coefs:
[5x4double]
pieces:
5
order:
4
dim:
1
ans=
1.4903-2.4903-1.00003.0000
-0.48791.9807-1.50971.0000
0.1795-0.94690.55802.0000
-0.01570.6691-0.27540
-0.78740.57492.21262.0000
计算结果说明该三次样条的分段多项式形式为
2,
（1）多项式插值
计算多项式插值的MATLAB的函数M文件为：
functionyi=polyinterp（x,y,xi）
n=length（x）;m=length（xi）;
fork=1:
m
z=xi（k）;s=0;
fori=1:
n
p=1;
forj=1:
n
ifj~=i
p=p*（z-x（j））/（x（i）-x（j））;
end
end
s=p*y（i）+s;
end
yi（k）=s;
end
以下为用来绘图的MATLAB脚本：
x=[0,3,5,7,9,11,12,13,14,15];y=[0,1.8,2.2,2.7,3.0,3.1,2.9,2.5,2.0,1.6];
s=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
xi=-1:
.01:
16;yi=zeros（size（xi））;
figure
（1）
fork=1:
10
w=zeros（1,10）;w（k）=w（k）+1;
wi=polyinterp（x,w,xi）;
yi=yi+y（k）.*wi;
subplot（5,2,k）,plot（x,w,'k*',xi,wi,'k'）,axis（[-1,16,-1,3]）
end
si=zeros（size（xi））;
figure
（2）
fork=1:
10
w=zeros（1,10）;w（k）=w（k）+1;
wi=polyinterp（x,w,xi）;
si=si+s（k）.*wi;
subplot（5,2,k）,plot（x,w,'k*',xi,wi,'k'）,axis（[-1,16,-1,3]）
end
figure（3）,
plot（x,y,'k*',xi,yi,'k'）,holdon
plot（x,s,'k*',xi,si,'k'）,holdoff,
axis（[-1,16,-1,5]）
title（'Lagrange插值多项式'）
gtext（'L_16（x）=1.8*I1+2.2*I2+2.7*I3+3*I4+3.1*I5+2.9*I6+2.5*I7+2*I8+1.6*I9'）;
gtext（'R_16（x）=1.2*I1+1.7*I2+2.0*I3+2.1*I4+2*I5+1.8*I6+1.2*I7+1*I8+1.6*I9'）
图像为：
figure1：
figure2：
figure3：
由图可知，高次的多项式插值会发生震荡，不可求面积。
（2）分段线性插值法
程序：
x=[0,3,5,7,9,11,12,13,14,15];y=[0,1.8,2.2,2.7,3.0,3.1,2.9,2.5,2.0,1.6];
s=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
xi=-1:
.1:
16;yi=zeros（size（xi））;
figure
（1）
fork=1:
10
w=zeros（1,10）;w（k）=w（k）+1;
wi=interp1（x,w,xi）;
yi=yi+y（k）.*wi;
subplot（5,2,k）,plot（x,w,'k*',xi,wi,'k'）,axis（[-1,16,-1,3]）;
end
b1=trapz（x,y）
si=zeros（size（xi））;
figure
（2）
fork=1:
10
w=zeros（1,10）;w（k）=w（k）+1;
wi=interp1（x,w,xi）;
si=si+s（k）.*wi;
subplot（5,2,k）,plot（x,w,'k+',xi,wi,'k'）,axis（[-1,16,-1,3]）
end
b2=trapz（x,s）
figure（3）,
plot（x,y,'ko',xi,yi,'-k'）,holdon
plot（x,s,'ko',xi,si,'--k'）,holdoff,
title（'Lagrange插值多项式'）,
axisequal
b1=
33.1500
b2=
22.4000
图像：
Figure1：
Figure2：
Figure3：
计算：
在这种方法下，截面面积为10.75
（3）三次样条插值
分析:
通过计算上,下轮廓的数值积分,然后作差得到截面面积
程序:
计算上轮廓的程序
x=[0,3,5,7,9,11,12,13,14,15];
y=[0,1.8,2.2,2.7,3.0,3.1,2.9,2.5,2.0,1.6];
pp=csape（x,y）,pp.coefs;
s=@（t,tj,c）c
（1）.*（t-tj）.^3+c
（2）.*（t-tj）.^2+c（3）.*（t-tj）+c（4）;
fork=1:
pp.pieces
c=pp.coefs（k,:
）;u=x（k）:
.1:
x（k+1）;v=s（u,x（k）,c）;
v1=d1s（u,x（k）,c）;v2=d2s（u,x（k）,c）;v3=d3s（u,x（k）,c）;
plot（u,v,'m'）,holdon
end
legend（'三次样条（边界条件2）'）
plot（x,y,'ko'）,holdoff
a=trapz（0:
.01:
15,ppval（pp,0:
.01:
15））,
axisequal
运行结果：
pp=
form:
'pp'
breaks:
[035791112131415]
coefs:
[9x4double]
pieces:
9
order:
4
dim:
1
a=
33.8639
图像：
计算下轮廓的程序
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];
pp=csape（x,y）,pp.coefs;
s=@（t,tj,c）c
（1）.*（t-tj）.^3+c
（2）.*（t-tj）.^2+c（3）.*（t-tj）+c（4）;
d1s=@（t,tj,c）3.*c
（1）.*（t-tj）.^2+2.*c
（2）.*（t-tj）+c（3）;
d2s=@（t,tj,c）6.*c
（1）.*（t-tj）+2.*c
（2）;
d3s=@（t,tj,c）6.*c
（1）.*ones（size（t））;
fork=1:
pp.pieces
c=pp.coefs（k,:
）;u=x（k）:
.1:
x（k+1）;v=s（u,x（k）,c）;
v1=d1s（u,x（k）,c）;v2=d2s（u,x（k）,c）;v3=d3s（u,x（k）,c）;
plot（u,v,'m',u,v1,'k--',u,v2,'k-',u,v3,'k:
'）,holdon
end
legend（'三次样条（边界条件2）','样条一阶导函数','样条二阶导函数','样条三阶导函数',4）
plot（x,y,'ko'）,holdoff
b=trapz（0:
.01:
15,ppval（pp,0:
.01:
15））,
axisequal
运行结果：
pp=
form:
'pp'
breaks:
[035791112131415]
coefs:
[9x4double]
pieces:
9
order:
4
dim:
1
b=
22.5722
图像：
结论：
截面面积（=a-b=33.8639-22.5722=）11.2917
4.解：
问题分析：
根据表5.5中测量的过桥车辆数据记录，然后估计一天过桥的车辆总数。
由题设可知一天过桥的车辆总数是每一时刻过桥的车辆数对时间积分得到。
可使用多项式、分段线性、三次样条三种一维插值的方法来进行求解，再用复化梯形求积公式来求面积，从而得出一天过桥辆车的总数。
显然，多项式插值一般适合于结点数目较少即多项式次数较低的情况，故不采用此种方法来进行插值。
所以我们选择用分段线性和三次样条插值方法进行插值，然后再进一步比较优劣，以更好地估计一天过桥的车辆总数。
模型的建立和求解：
先用Matlab编程进行线性分段插值和三次样条插值，再用复化梯形求积公式求得一天过桥的车的数量，同时得到对应的图像。
程序如下：
x=[0245678910.512.514161718192021222324];
y=[2202582512101279282210911893];
x0=x+（1/60）*（1/2）;
y0=y/（1/60）;
u=（1/60）*（1/2）:
（1/60）:
24+（1/60）*（1/2）;
y1=interp1（x0,y0,u）;
sum1=trapz（y1）*（1/60）
plot（u,y1,'k'）,grid,xlabel（'time/h'）,
ylabel（'一天过桥的车的数量'）,title（'线性分段插值'）
运行结果：
sum1=13425
y2=spline（x0,y0,u）;
sum2=trapz（y2）*（1/60）
plot（u,y2,'k'）,grid,xlabel（'time/h'）,
ylabel（'一天过桥的车的数量'）,title（'三次样条插值'）
运行结果：
sum2=1.3192e+004
结果分析：
用分段线性插值方法和三次样条插值方法进行插值进行插值，用复化梯形求积公式求面积，得到一天通过桥梁的车总量分别为13425辆和13192辆。
比较图形，可以知道分段线性插值方法得出的结果更符合实际情况，拟合得更好更合理，于是得出一天过桥的车的数量是13425辆。
5.解：
问题分析：
设这辆汽车在t时刻的位置为x（t），则x（t）关于t的一阶导数
是这辆汽车在t时刻的速度，由表5.6可知x（t）在区间[0,13]的左、右端点的一阶导数就是在t=0s和t=13s时汽车的速度，均为20m/s，即
。
由上面的分析，首先求出x（t）的函数表达式，再对其求导得到一阶导数
，结合
和
的图像得到所要求解的问题的答案即可。
模型的建立和求解：
用Matlab样条工具箱函数scape计算给定左右端点一阶导数的三次样条，得到位置x关于时间t的函数。
Matlab程序如下：
x=[0,3,5,8,13];y=[20,0,65,121,194,313,20];
pp=csape（x,y,'complete'）,pp.coefs
s=@（t,tj,c）c
（1）.*（t-tj）.^3+c
（2）.*（t-tj）.^2+c（3）.*（t-tj）+c（4）;
d1s=@（t,tj,c）3.*c
（1）.*（t-tj）.^2+2.*c
（2）.*（t-tj）+c（3）;
d2s=@（t,tj,c）6.*c
（1）.*（t-tj）+2.*c
（2）;
fork=1:
pp.pieces
c=pp.coefs（k,:
）;u=x（k）:
0.01:
x（k+1）;v=s（u,x（k）,c）;
v1=d1s（u,x（k）,c）;v2=d2s（u,x（k）,c）;plot（u,v,'k',u,v1,'k-',u,v2,'k--'）,holdon
end
legend（'三次样条','样条的一阶导函数',...
'样条的二阶导函数'）
y1=[0,65,121,194,313];
plot（x,y1,'ko'）,holdoff
输出的结果：
pp=
form:
'pp'
breaks:
[035813]
coefs:
[4x4double]
pieces:
4
order:
4
dim:
1
ans=
0.3008-0.346920.00000
-0.69042.360526.040765.0000
0.2757-1.782027.1976121.0000
-0.14600.699723.9507194.0000
由此可以得到这辆汽车在t时刻的位置的分段函数x（t）如下：
对x（t）求导数得到导函数
也即速度关于t的函数如下：
（1）求t=10s时，这辆车的位置和速度
此题只需将t=10s代入
得到
也即在
时的位置是
，速度是
。
（2）由题意汽车在这段公路上行驶时，最大速度不能超过
，下面用matlab编程画出速度函数v（t）关于时间t的图像，matlab程序如下：
f1=@（x）0.9024*x.^2-0.6938*x+20;
f2=@（x）-2.0712*（x-3）.^2+4.721*（x-3）+26.0407;
f3=@（x）0.8271*（x-5）.^2-3.564*（x-5）+27.1976;
f4=@（x）-0.438*（x-8）.^2+1.3994*（x-8）+23.9507;
plot（0:
0.01:
3,f1（0:
0.01:
3）,3:
0.01:
5,f2（3:
0.01:
5）,5:
0.01:
8,f3（5:
0.01:
8）,8:
0.01:
13,f4（8:
0.01:
13）,0:
0.01:
14,22.2）,holdon
legend（'0≤t≤3','3≤t≤5','5≤t≤8','8≤t≤13'）
xlabel（'时间（s）'）,ylabel（'速度（m/s）'）
title（'速度v（t）关于时间t的函数图5.2'）
运行结果为：
由图像可以看出在
和
这两个时间段里面，汽车的速度已经超出了
，要求出汽车开始超速和结束超速的时间，只需要求出在在
和
这两个时间段里速度函数与
的交点的横坐标即可。
Matlab程序如下：
h1=@（x）f1（x）-22.2;
h2=@（x）f4（x）-22.2;
x1=fzero（h1,[0,3]）
x2=fzero（h2,[8,13]）
命令窗口输出结果是：
x1=1.9924，x2=12.1566。
因此，当
时，汽车开始超速，当
时，汽车结束超速，即当
时，汽车超速。
（3）由第
（2）题易知这辆汽车的最高速度在
这个时间段里面出现。
因此，只需求出在
内速度的最大值即可。
下面用matlab程序算出-f2的最小值点坐标（就可以求出f2的最大值点）
h3=@（x）-f2（x）;
[u,v]=fminbnd（h3,3,5）
命令窗口输出如下：
u=4.1397，v=-28.7309
即在
时，速度达到最大值，为
。
升级会员 