章绍辉版数学建模第五章作业.docx
- 文档编号:11912387
- 上传时间:2023-04-16
- 格式:DOCX
- 页数:22
- 大小:166.50KB
章绍辉版数学建模第五章作业.docx
《章绍辉版数学建模第五章作业.docx》由会员分享,可在线阅读,更多相关《章绍辉版数学建模第五章作业.docx(22页珍藏版)》请在冰豆网上搜索。
章绍辉版数学建模第五章作业
第五章作业
第一题:
(1)第三种边界条件
第三种边界条件要求给定三次样条s(x)在区间[x0,xn]的左右端点的一阶导数
和
。
我们参考了例5.1.4在210页计算出了左右端点的一阶导数,其中
=-3.3667、
=2.3333
MATLAB程序如下:
x=[0,1,3,6,8,9];y=[-3.3667,3,1,2,0,2,4,2.33333];
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
title('第三种边界条件及其一、二、三阶导函数的图像')
legend('三次样条(第三种边界条件)','样条的一阶导函数','样条的二阶导函数','样条的三阶导函数')
plot([0,1,3,6,8,9],[3,1,2,0,2,4],'ko'),holdoff
运行结果为:
pp=
form:
'pp'
breaks:
[013689]
coefs:
[5x4double]
pieces:
5
order:
4
dim:
1
ans=
-0.03891.4056-3.36673.0000
-0.35141.2890-0.67221.0000
0.1696-0.81970.26642.0000
-0.08480.7064-0.07360
0.06780.19771.73452.0000
所绘制的图形如下:
结果说明:
计算结果说明该三次样条的分段多项式为:
S(x)=
执行以下命令可以验算该三次样条在区间[0,9]的左端点x=0和右端点x=9的一阶导数分别为-3.3667和2.3333:
在MATLAB的commandwindow运行:
[1.*pp.coefs(1,3),d1s(9,8,pp.coefs(5,:
))]
运行结果为:
ans=
-3.36672.3333
(2)第四种边界条件
第四种边界条件要求给定三次样条s(x)在区间[x0,xn]的左右端点的二阶导数
和
。
我们分别取左右端点的二阶导数为
=-2,
=3
MATLAB程序如下:
x=[0,1,3,6,8,9];y=[-2,3,1,2,0,2,4,3];
pp=csape(x,y,'second')
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
title('第四种边界条件及其一、二、三阶导函数的图像')
legend('三次样条(第四种边界条件)','样条的一阶导函数','样条的二阶导函数','样条的三阶导函数')
plot([0,1,3,6,8,9],[3,1,2,0,2,4],'ko'),holdoff
运行结果为:
pp=
form:
'pp'
breaks:
[013689]
coefs:
[5x4double]
pieces:
5
order:
4
dim:
1
ans=
0.9088-1.0000-1.90883.0000
-0.44271.7265-1.18231.0000
0.1901-0.92960.41162.0000
-0.13190.7809-0.03440
0.5034-0.01031.50692.0000
所绘制的图形为:
结果说明:
计算结果说明该三次样条的分段多项式为:
S(x)=
(3)第五种边界条件
MATLAB程序如下:
x=[0,1,3,6,8,9];y=[3,1,2,0,2,4];
pp=csape(x,y,'not-a-kont')
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
title('第五种边界条件及其一、二、三阶导函数的图像')
legend('三次样条(第五种边界条件)','样条的一阶导函数','样条的二阶导函数','样条的三阶导函数')
plot([0,1,3,6,8,9],[3,1,2,0,2,4],'ko'),holdoff
运行结果为:
pp=
form:
'pp'
breaks:
[013689]
coefs:
[5x4double]
pieces:
5
order:
4
dim:
1
ans=
-0.32402.1291-3.80523.0000
-0.32401.1573-0.51881.0000
0.1633-0.78640.22292.0000
-0.07000.6833-0.08660
-0.07000.26331.80662.0000
所绘制的图形为:
结果说明:
计算结果说明该三次样条的分段多项式为:
S(x)=
(4)第六种边界条件
MATLAB程序如下:
x=[0,1,3,6,8,9];y=[3,1,2,0,2,4];
pp=csape(x,y,'periodic')
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
title('第六种边界条件及其一、二、三阶导函数的图像')
legend('三次样条(第六种边界条件)','样条的一阶导函数','样条的二阶导函数','样条的三阶导函数')
plot([0,1,3,6,8,9],[3,1,2,0,2,4],'ko'),holdoff
运行结果为:
pp=
form:
'pp'
breaks:
[013689]
coefs:
[5x4double]
pieces:
5
order:
4
dim:
1
ans=
1.9961-3.7833-0.21273.0000
-0.52962.2048-1.79121.0000
0.1754-0.97280.67282.0000
0.05370.6061-0.42720
-1.57060.92852.64212.0000
所绘制的图形如下:
结果说明:
计算结果说明该三次样条的分段多项式为:
S(x)=
第二题:
本题要求通过给出的机翼断面轮廓线的一些点的坐标绘制机翼断面轮廓线的图像并计算机翼断面的面积,我们分别运用了拉格朗日插值法、分段插值法和三次样条插值法来得到机翼断面上下轮廓线的函数并绘制出图像,然后再运用梯形法则对机翼断面上下轮廓线的函数进行积分,上轮廓线的积分减去下轮廓线的积分就得到机翼断面的面积。
(1)拉格朗日插值法
MATLAB程序如下:
x=[0,3,5,7,9,11,12,13,14,15];
y1=[0,1.8,2.2,2.7,3.0,3.1,2.9,2.5,2.0,1.6];
y2=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
xi=0:
.05:
15;yi=zeros(size(xi));ti=zeros(size(xi));
yi=polyval(polyfit(x,y1,9),xi);
ti=polyval(polyfit(x,y2,9),xi);
plot(x,y1,'k.',x,y2,'k.',xi,yi,'k-',xi,ti,'k-');
axis([0,15,-17,5]);
title('多项式插值机翼横断面轮廓线');
xlabel('x');ylabel('y');
s=trapz(xi,yi)-trapz(xi,ti)
运行结果:
Warning:
Polynomialisbadlyconditioned.Removerepeateddatapoints
ortrycenteringandscalingasdescribedinHELPPOLYFIT.
>Inpolyfitat81
Warning:
Polynomialisbadlyconditioned.Removerepeateddatapoints
ortrycenteringandscalingasdescribedinHELPPOLYFIT.
>Inpolyfitat81
s=
40.3426
所绘制的图形如下:
(多项式插值机翼横断面轮廓线)
(2)分段线性差值
MATLAB程序如下:
x=[0,3,5,7,9,11,12,13,14,15];
y1=[0,1.8,2.2,2.7,3.0,3.1,2.9,2.5,2.0,1.6];
y2=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
xi=0:
.05:
15;yi=zeros(size(xi));ti=zeros(size(xi));
yi=interp1(x,y1,xi);
ti=interp1(x,y2,xi);
plot(x,y1,'k*',x,y2,'k*',xi,yi,'k-',xi,ti,'k-'),
axis([0,15,-1,5]),xlabel('x');ylabel('y');
title('分段插值机翼横断面轮廓线')
s=trapz(xi,yi)-trapz(xi,ti)
运行结果为:
s=
10.7500
所绘制的图形如下:
(分段插值机翼横断面轮廓线)
(3)三次样条差值
MATLAB程序如下:
x=[0,3,5,7,9,11,12,13,14,15];
y1=[0,1.8,2.2,2.7,3.0,3.1,2.9,2.5,2.0,1.6];
y2=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
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));
pp=csape(x,y1),pp.coefs;
s1=0;
fork=1:
pp.pieces,
c=pp.coefs(k,:
);u=x(k):
.01:
x(k+1);v=s(u,x(k),c);
s1=s1+trapz(u,v);
plot(u,v,'k'),holdon
end
s2=0;
pp=csape(x,y2);pp.coefs;
fork=1:
pp.pieces
c=pp.coefs(k,:
);u=x(k):
.01:
x(k+1);v=s(u,x(k),c);
s2=s2+trapz(u,v);
plot(u,v,'k'),holdon
end
plot(x,y1,'ko'),plot(x,y2,'ko'),holdoff
xlabel('x');ylabel('y');
title('三次样条插值法机翼断面轮廓线')
S=s1-s2
运行结果为:
S=
11.2917
所绘制的图形为:
(三次样条插值法机翼断面轮廓线)
结果分析:
综上,我们得到用拉格朗日插值法得到的机翼断面面积为40.3426,用分段插值法得到的机翼断面面积为10.75,用三次样条插值法得到的机翼断面面积为11.2917。
拉格朗日插值是比较常见的一种函数插值,但这里拉格朗日插值会发生Runge现象,得到的机翼断面面积明显偏大,不可采用。
分段插值可以避免高次插值可能出现的大幅度波动现象,克服了拉格朗日插值可能发生不收敛的缺点,但分段插值存在基点处不光滑、插值精度低的缺点。
而三次样条插值不但克服了拉格朗日插值的不收敛性,也提高了分段插值函数在节点处的光滑性,所以说三次样条插值是较低次数的多项式而达到较高阶光滑性的方法。
第四题:
我们以时间为横坐标,车辆数为纵坐标,得到20个结点:
(0,2)、(2,2)、(4,0)、(5,2)、(6,5)、(7,8)、(8,25)、(9,12)、(10.5,10)、(12.5,12)、(14,7)、(16,9)、(17,28)、(18,22)、(19,10)、(20,9)、(21,11)、(22,8)、(23,9)、(24,3)
利用了三次样条插值的方法(选用第二种边界条件),再用梯形法则计算每一分段函数的积分进行循环叠加,最后得出一天过桥的车辆数。
MATLAB程序如下:
x=[0,2,4,5,6,7,8,9,10.5,12.5,14,16,17,18,19,20,21,22,23,24];
y=[2,2,0,2,5,8,25,12,10,12,7,9,28,22,10,9,11,8,9,3];
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));
S=0;
fork=1:
pp.pieces
c=pp.coefs(k,:
);u=x(k):
.01:
x(k+1);v=s(u,x(k),c);
S=S+trapz(u,v);
plot(u,v,'k'),holdon
end
plot(x,y,'ko'),holdoff
s=S
xlabel('时间')
ylabel('车辆数')
title('一天内各时刻过桥车辆数据图')
运行结果为:
pp=
form:
'pp'
breaks:
[0245678910.500012.500014161718192021222324]
coefs:
[19x4double]
pieces:
19
order:
4
dim:
1
s=
219.6689
所绘制的图形为:
结果分析:
所以一天大约有220辆汽车过桥。
第五题:
我们以时间为横坐标,位置为纵坐标,得到5个结点:
x1(0,0)、x2(3,65)、x3(5,121)、x4(8,194)、x5(13,313).那么位置x关于时间t的函数在这5个节点处的导数为:
.因此,我们分四段x1x2,x2x3,x3x4,x4x5进行三次样条差值,得到每一段的三次样条函数,并画出
和
(该车的速度v关于时间t的函数)图像,进一步求出该车的最大速度和超速时间段。
MATLAB程序如下:
t=[035813];y=[065121194313];v=[2026272420];
pp1=csape([t
(1),t
(2)],[v
(1),y
(1),y
(2),v
(2)],'complete');
pp2=csape([t
(2),t(3)],[v
(2),y
(2),y(3),v(3)],'complete');
pp3=csape([t(3),t(4)],[v(3),y(3),y(4),v(4)],'complete');
pp4=csape([t(4),t(5)],[v(4),y(4),y(5),v(5)],'complete');
c=[pp1.coefs;pp2.coefs;pp3.coefs;pp4.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));
y_10=s(10,8,c(4,:
))
v_10=d1s(10,8,c(4,:
))
figure
(1)
fork=1:
4
u=t(k):
.01:
t(k+1);p=s(u,t(k),c(k,:
));
plot(u,p,'k'),holdon
end
plot(t,y,'ko'),holdoff
xlabel('时间/s')
ylabel('位置/m')
title('该汽车在0到13秒内位置数据图')
figure
(2)
fork=1:
4
u=t(k):
.01:
t(k+1);p=d1s(u,t(k),c(k,:
));
plot(u,p,'k'),holdon
end
plot(t,v,'ko')
plot([0,13],[80/3.6,80/3.6],'k--'),holdoff
xlabel('时间/s')
ylabel('速度m/s')
title('该汽车在0到13秒内速度数据图')
运行结果为:
c=
0.2963-0.333320.00000
-0.75002.500026.000065.0000
0.2593-1.666727.0000121.0000
-0.14400.680024.0000194.0000
y_10=
243.5680
v_10=
24.9920
所绘制的图形如下:
(1)当t=10s时,这辆汽车的位置大约为243.56m,速度大约为24.99m/s
(2)由上图可以直观的看到,与速度上限v=22.22m/s的交点在第一段和第四段。
有前面计算结果可以得到这两段的三次样条函数为
和
则由
,代入数据化简得到两个方程,通过解方程
和
,可以求出交点。
执行代码:
s=solve('0.8889*x^2-0.6666*x-2.2222=0')
s1=solve('-0.4320*x^2+8.272*x-36.7502=0')
运行得到:
s=
-1.2500151443640650147734304226318
1.9999307704187393313444732845961
s1=
7.0063928305296322024009166077546
12.141755317618515945747231540394
结果分析:
这辆汽车在2s—12.14s间超速了。
(3)由该车的速度图像可知,速度最大值在第二段取到。
执行代码:
d1s=@(t,c)3.*c
(1).*(t-3).^2+2.*c
(2).*(t-3)+c(3);
t=3:
0.0001:
5;
c=[-0.7500,2.5000,26.0000,65.0000];
[vmax,tmax]=max(d1s(t,c))
tmax1=3+tmax*0.0001
运行得到:
vmax=
28.7778
tmax=
11112
tmax1=
4.1112
结果分析:
在观测时间段内,这辆车的最高速度为vmax=28.7778m/s,此时t=4.11s。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 章绍辉版 数学 建模 第五 作业