章绍辉版数学建模第五章作业

来源：个人技术集锦

第五章作业

第一题：

（1）第三种边界条件

第三种边界条件要求给定三次样条s(x)在区间[x0,xn]的左右端点的一阶导数S'(x0)和S'(xa)。我们参考了例5.1.4在210页计算出了左右端点的一阶导数，其中S'(x0)= -3.3667、S'(xa)= 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)); for k=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:'),hold on end

title('第三种边界条件及其一、二、三阶导函数的图像')

legend('三次样条（第三种边界条件）','样条的一阶导函数','样条的二阶导函数','样条的三阶导函数')

plot([0,1,3,6,8,9],[3,1,2,0,2,4],'ko'),hold off

运行结果为： pp =

form: 'pp'

breaks: [0 1 3 6 8 9] coefs: [5x4 double] pieces: 5 order: 4 dim: 1 ans =

-0.0389 1.4056 -3.3667 3.0000 -0.3514 1.2890 -0.6722 1.0000 0.1696 -0.8197 0.2664 2.0000 -0.0848 0.7064 -0.0736 0

0.0678 0.1977 1.7345 2.0000

所绘制的图形如下：

第三种边界条件及其一、二、三阶导函数的图像43210-1-2-3-4三次样条（第三种边界条件）样条的一阶导函数样条的二阶导函数样条的三阶导函数0123456789

结果说明：

计算结果说明该三次样条的分段多项式为：

0.0389x31.4056x23.3667x3,0x10.3514(x1)31.289(x1)20.6722(x1)1,1x3S(x)= 0.1696(x3)30.8197(x3)20.2664(x3)2,3x6

0.0848(x6)30.7064(x6)20.0736(x6),6x80.0678(x8)30.1977(x8)21.7345(x8)2,8x9执行以下命令可以验算该三次样条在区间[0，9]的左端点x=0和右端点x=9的一阶导数分别为-3.3667和2.3333：

在MATLAB的command window运行：[1.*pp.coefs(1,3),d1s(9,8,pp.coefs(5,:))] 运行结果为： ans =

-3.3667 2.3333

（2）第四种边界条件

第四种边界条件要求给定三次样条s(x)在区间[x0,xn]的左右端点的二阶导数s''(x0)和s''(xn)。我们分别取左右端点的二阶导数为s''(x0)=-2,s''(xn)=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)); for k=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:'),hold on end

title('第四种边界条件及其一、二、三阶导函数的图像')

legend('三次样条（第四种边界条件）','样条的一阶导函数','样条的二阶导函数','样条的三阶导函数')

plot([0,1,3,6,8,9],[3,1,2,0,2,4],'ko'),hold off 运行结果为： pp =

form: 'pp'

breaks: [0 1 3 6 8 9] coefs: [5x4 double] pieces: 5 order: 4 dim: 1 ans =

0.9088 -1.0000 -1.9088 3.0000 -0.4427 1.7265 -1.1823 1.0000 0.1901 -0.9296 0.4116 2.0000 -0.1319 0.7809 -0.0344 0

0.5034 -0.0103 1.5069 2.0000 所绘制的图形为：

第四种边界条件及其一、二、三阶导函数的图像6三次样条（第四种边界条件）543210-1-2-3样条的一阶导函数样条的二阶导函数样条的三阶导函数0123456789

结果说明：

计算结果说明该三次样条的分段多项式为：

0.9088x3x21.9088x3,0x10.4427(x1)31.7265(x1)21.1823(x1)1,1x3S(x)= 0.1901(x3)30.9296(x3)20.4116(x3)2,3x6

0.1319(x6)30.7809(x6)20.0344(x6),6x80.5034(x8)30.0103(x8)21.5069(x8)2,8x9

（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

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:'),hold on end

title('第五种边界条件及其一、二、三阶导函数的图像')

legend('三次样条（第五种边界条件）','样条的一阶导函数','样条的二阶导函数','样条的三阶导函数')

plot([0,1,3,6,8,9],[3,1,2,0,2,4],'ko'),hold off

运行结果为： pp =

form: 'pp'

breaks: [0 1 3 6 8 9] coefs: [5x4 double] pieces: 5 order: 4 dim: 1 ans =

-0.3240 2.1291 -3.8052 3.0000 -0.3240 1.1573 -0.5188 1.0000 0.1633 -0.7864 0.2229 2.0000 -0.0700 0.6833 -0.0866 0 -0.0700 0.2633 1.8066 2.0000

所绘制的图形为：

第五种边界条件及其一、二、三阶导函数的图像5三次样条（第五种边界条件）43210-1-2-3-4样条的一阶导函数样条的二阶导函数样条的三阶导函数0123456789

结果说明：

计算结果说明该三次样条的分段多项式为：

0.324x32.1291x23.8052x3,0x10.324(x1)31.1573(x1)20.5188(x1)1,1x3S(x)= 0.1633(x3)30.7864(x3)20.2229(x3)2,3x6

0.07(x6)30.6833(x6)20.0866(x6),6x80.07(x8)30.2633(x8)21.8066(x8)2,8x9（4）第六种边界条件

MATLAB程序如下：

x=[0,1,3,6,8,9];y=[3,1,2,0,2,4]; pp=csape(x,y,'periodic') pp.coefs

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:'),hold on end

title('第六种边界条件及其一、二、三阶导函数的图像')

legend('三次样条（第六种边界条件）','样条的一阶导函数','样条的二阶导函数','样条

的三阶导函数')

plot([0,1,3,6,8,9],[3,1,2,0,2,4],'ko'),hold off 运行结果为： pp =

form: 'pp'

breaks: [0 1 3 6 8 9] coefs: [5x4 double] pieces: 5 order: 4 dim: 1 ans =

1.9961 -3.7833 -0.2127 3.0000 -0.5296 2.2048 -1.7912 1.0000 0.1754 -0.9728 0.6728 2.0000 0.0537 0.6061 -0.4272 0 -1.5706 0.9285 2.6421 2.0000 所绘制的图形如下：

第六种边界条件及其一、二、三阶导函数的图像15三次样条（第六种边界条件）样条的一阶导函数样条的二阶导函数样条的三阶导函数1050-5-100123456789

结果说明：

计算结果说明该三次样条的分段多项式为：

1.9961x33.7833x20.2127x3,0x10.5296(x1)32.2048(x1)21.7912(x1)1,1x3S(x)= 0.1754(x3)0.9728(x3)0.6728(x3)2,3x6

320.0537(x6)30.6061(x6)20.4272(x6),6x81.5706(x8)30.9285(x8)22.6421(x8)2,8x9第二题：

本题要求通过给出的机翼断面轮廓线的一些点的坐标绘制机翼断面轮廓线的图像并计算机翼断面的面积，我们分别运用了拉格朗日插值法、分段插值法和三次样条插值法来得到机翼断面上下轮廓线的函数并绘制出图像，然后再运用梯形法则对机翼断面上下轮廓线的函数进行积分，上轮廓线的积分减去下轮廓线的积分就得到机翼断面的面积。

（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: Polynomial is badly conditioned. Remove repeated data points or try centering and scaling as described in HELP POLYFIT. > In polyfit at 81

Warning: Polynomial is badly conditioned. Remove repeated data points or try centering and scaling as described in HELP POLYFIT. > In polyfit at 81 s =

40.3426

所绘制的图形如下：（多项式插值机翼横断面轮廓线）

多项式插值机翼横断面轮廓线420-2-4-6-8-10-12-14-1605x1015y

（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

所绘制的图形如下：（分段插值机翼横断面轮廓线）

分段插值机翼横断面轮廓线5432y10-105x1015

（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];

for k=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'),hold on end s2=0;

pp=csape(x,y2);pp.coefs; for k=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'),hold on end

plot(x,y1,'ko'),plot(x,y2,'ko'),hold off xlabel('x');ylabel('y');

title('三次样条插值法机翼断面轮廓线') S=s1-s2

运行结果为： S =

11.2917

所绘制的图形为：（三次样条插值法机翼断面轮廓线）

三次样条插值法机翼断面轮廓线3.532.52y1.510.5005x1015

结果分析：综上，我们得到用拉格朗日插值法得到的机翼断面面积为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;

for k=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'),hold on end

plot(x,y,'ko'),hold off s=S

xlabel('时间') ylabel('车辆数')

title('一天内各时刻过桥车辆数据图')

运行结果为： pp =

form: 'pp'

breaks: [0 2 4 5 6 7 8 9 10.5000 12.5000 14 16 17 18 19 20 21 22 23 24] coefs: [19x4 double] pieces: 19 order: 4 dim: 1 s =

219.6689

所绘制的图形为：

一天内各时刻过桥车辆数据图302520车辆数151050-50510时间152025

结果分析：所以一天大约有220辆汽车过桥。

第五题：

我们以时间为横坐标，位置为纵坐标，得到5个结点：x1(0,0)、x2(3,65)、x3(5,121)、x4(8,194)、x5(13,313).那么位置x关于时间t的函数在这5个节点处的导数为：

x'(0)20,x'(3)26,x'(5)27,x'(8)24,x'(13)20.因此，我们分四段

x1x2,x2x3,x3x4,x4x5进行三次样条差值，得到每一段的三次样条函数，并画出x(t)和

进一步求出该车的最大速度和超速时间v(t)x'(t)(该车的速度v关于时间t的函数)图像，

段。

MATLAB程序如下：

t=[0 3 5 8 13];y=[0 65 121 194 313];v=[20 26 27 24 20]; 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) for k=1:4

u=t(k):.01:t(k+1);p=s(u,t(k),c(k,:)); plot(u,p,'k'),hold on end

plot(t,y,'ko'),hold off xlabel('时间/s') ylabel('位置/m')

title('该汽车在0到13秒内位置数据图') figure(2) for k=1:4

u=t(k):.01:t(k+1);p=d1s(u,t(k),c(k,:)); plot(u,p,'k'),hold on end

plot(t,v,'ko')

plot([0,13],[80/3.6,80/3.6],'k--'),hold off xlabel('时间/s') ylabel('速度m/s')

title('该汽车在0到13秒内速度数据图') 运行结果为： c =

0.2963 -0.3333 20.0000 0 -0.7500 2.5000 26.0000 65.0000 0.2593 -1.6667 27.0000 121.0000 -0.1440 0.6800 24.0000 194.0000 y_10 = 243.5680 v_10 = 24.9920

所绘制的图形如下：

该汽车在0到13秒内位置数据图350300250位置/m2001501005000246时间/s8101214

该汽车在0到13秒内速度数据图2928272625速度m/s2423222120190246时间/s8101214X: 13Y: 22.22

(1) 当t=10s时，这辆汽车的位置大约为243.56m，速度大约为24.99m/s

(2) 由上图可以直观的看到，与速度上限v=22.22m/s的交点在第一段和第四段。有前面计

算结果可以得到这两段的三次样条函数为

v(t)3c1t22c2tc3,其中[c1,c2,c3][0.2963,0.3333,20.0000],

和 v(t)3c1(t8)2c2(t8)c3,其中[c1,c2,c3][0.144,0.68,24],

2则由v(t)800，代入数据化简得到两个方程，通过解方程 3.60.8889t20.6666t2.22220 和-0.432t28.272t36.75020 ，可以求出交点。

执行代码：

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。

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文

全部栏目

章绍辉版数学建模第五章作业