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

1、章绍辉版数学建模第五章作业第五章作业第一题：（1） 第三种边界条件第三种边界条件要求给定三次样条s(x)在区间x0,xn的左右端点的一阶导数和。我们参考了例5.1.4在210页计算出了左右端点的一阶导数，其中= -3.3667、= 2.3333MATLAB程序如下：x=0,1,3,6,8,9;y=-3.3667,3,1,2,0,2,4,2.33333;pp=csape(x,y,complete)pp.coefss=(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、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 onendtitle(第三种边界条件及其一、二、三阶导函数的图像)legend(三次样条（第三种边界条件）

3、,样条的一阶导函数,样条的二阶导函数,样条的三阶导函数)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: 1ans = -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 00.0678 0.1977 1.7345 2.0000所绘

4、制的图形如下：结果说明：计算结果说明该三次样条的分段多项式为：S(x)= 执行以下命令可以验算该三次样条在区间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的左右端点的二阶导数和。我们分别取左右端点的二阶导数为=-2, =3MATLAB程序如下：x=0,1,3,6,8,9;y=-2,3,1,2,0,2,4,3;p

5、p=csape(x,y,second)pp.coefss=(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

6、(k),c);v3=d3s(u,x(k),c); plot(u,v,k,u,v1,k-.,u,v2,k-,u,v3,k:),hold onendtitle(第四种边界条件及其一、二、三阶导函数的图像)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: 1ans = 0.9088 -1.0000 -1.9088

7、 3.0000 -0.4427 1.7265 -1.1823 1.0000 0.1901 -0.9296 0.4116 2.0000 -0.1319 0.7809 -0.0344 00.5034 -0.0103 1.5069 2.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.coefss=(t,tj,c)c(1).*(t-tj).3+c(2).*(t-tj).2+c(3).*(t-tj)+c(4);d1s=

8、(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 onendtitle(第五种边界条件及其一、

9、二、三阶导函数的图像)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: 1ans = -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

10、 0 -0.0700 0.2633 1.8066 2.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.coefss=(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,

11、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 onendtitle(第六种边界条件及其一、二、三阶导函数的图像)legend(三次样条（第六种边界条件）,样条的一阶导函数,样条的二阶导函数,样条的三阶导函数)plot(0,1,3,6,8,9,3,1,2,0,2,4,ko

12、),hold off运行结果为：pp = form: pp breaks: 0 1 3 6 8 9 coefs: 5x4 double pieces: 5 order: 4 dim: 1ans = 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所绘制的图形如下：结果说明：计算结果说明该三次样条的分段多项式为：S(x)= 第二题：本题要求通过给出的机翼断面轮廓线的

13、一些点的坐标绘制机翼断面轮廓线的图像并计算机翼断面的面积，我们分别运用了拉格朗日插值法、分段插值法和三次样条插值法来得到机翼断面上下轮廓线的函数并绘制出图像，然后再运用梯形法则对机翼断面上下轮廓线的函数进行积分，上轮廓线的积分减去下轮廓线的积分就得到机翼断面的面积。（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=zer

14、os(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 scali

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

16、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所绘制的图形如下：（分段插值

17、机翼横断面轮廓线）（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

18、(1).*ones(size(t);pp=csape(x,y1),pp.coefs;s1=0;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 onends2=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 onendplot(x,y1

19、,ko),plot(x,y2,ko),hold offxlabel(x);ylabel(y); title(三次样条插值法机翼断面轮廓线)S=s1-s2运行结果为：S = 11.2917所绘制的图形为：（三次样条插值法机翼断面轮廓线）结果分析：综上，我们得到用拉格朗日插值法得到的机翼断面面积为40.3426，用分段插值法得到的机翼断面面积为10.75，用三次样条插值法得到的机翼断面面积为11.2917。拉格朗日插值是比较常见的一种函数插值，但这里拉格朗日插值会发生Runge现象，得到的机翼断面面积明显偏大，不可采用。分段插值可以避免高次插值可能出现的大幅度波动现象，克服了拉格朗日插值可能发生不

20、收敛的缺点，但分段插值存在基点处不光滑、插值精度低的缺点。而三次样条插值不但克服了拉格朗日插值的不收敛性，也提高了分段插值函数在节点处的光滑性，所以说三次样条插值是较低次数的多项式而达到较高阶光滑性的方法。第四题：我们以时间为横坐标，车辆数为纵坐标，得到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) 利用了三次样条插值的方法（选用第二种边界条件）

21、，再用梯形法则计算每一分段函数的积分进行循环叠加，最后得出一天过桥的车辆数。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)

22、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 onendplot(x,y,ko),hold offs=Sxlabel(时间)ylabel(车辆数)title(一天内各时刻过桥车辆数据图)运行结果为：pp = form: pp breaks: 0 2 4 5 6 7 8 9 10.5000 12.5000 14 16 17

23、 18 19 20 21 22 23 24 coefs: 19x4 double pieces: 19 order: 4 dim: 1s = 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的函数)图像，进一步求出该车的最大速度和超速时间

24、段。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.co

25、efss=(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:4u=t(k):.01:t(k+1);p=s(u,t(k),c(k,:);plot(u,p,k),hold onendplot

26、(t,y,ko),hold offxlabel(时间/s)ylabel(位置/m)title(该汽车在0到13秒内位置数据图)figure(2)for k=1:4u=t(k):.01:t(k+1);p=d1s(u,t(k),c(k,:);plot(u,p,k),hold onendplot(t,v,ko)plot(0,13,80/3.6,80/3.6,k-),hold offxlabel(时间/s)ylabel(速度m/s)title(该汽车在0到13秒内速度数据图)运行结果为：c = 0.2963 -0.3333 20.0000 0 -0.7500 2.5000 26.0000 65.000

27、0 0.2593 -1.6667 27.0000 121.0000 -0.1440 0.6800 24.0000 194.0000y_10 = 243.5680v_10 = 24.9920所绘制的图形如下：(1) 当t=10s时，这辆汽车的位置大约为243.56m，速度大约为24.99m/s(2) 由上图可以直观的看到，与速度上限v=22.22m/s的交点在第一段和第四段。有前面计算结果可以得到这两段的三次样条函数为和 则由，代入数据化简得到两个方程，通过解方程 和 ，可以求出交点。执行代码：s=solve(0.8889*x2-0.6666*x-2.2222=0)s1=solve(-0.432

28、0*x2+8.272*x-36.7502=0)运行得到：s = -1.2500151443640650147734304226318 1.9999307704187393313444732845961s1 = 7.0063928305296322024009166077546 12.141755317618515945747231540394结果分析：这辆汽车在2s12.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.7778tmax = 11112tmax1 = 4.1112结果分析：在观测时间段内，这辆车的最高速度为vmax =28.7778m/s,此时t=4.11s。