matlab数学实验作业.doc

深圳大学考试答题纸(以论文、报告等形式考核专用)二一 二一 一 学年度第 一 学期课程编号23160008课程名称Matlab与数学实验主讲教师仇高新评分学 号 2008160120 姓名黄日才专业年级 电子学院微电子二班教师评语：题目：Matlab与数学实验期末作业1. 用数值积分公式计算 (结果保留小数点后8位)：(1) 取积分步长, 用梯形公式计算S= 6.247641317417333 。 x=0:pi/2:2*pi; trapz(0:pi/2:2*pi, sqrt(1-(0.15.*sin(x).2)(2) 要求相对误差为10-6, 用Simpson公式S= 6.247691887569109 ，Matlab命令是_ quad(sqrt(1-(0.15.*sin(x).2),0,2*pi,1e-6)_. 2. 设用数值解法算出 y(1)= 1.163536246222507 ,你用的方法是 Runga-Kutta 方法 ，调用的 Matlab命令t,x=ode45(verderpol2,0:0.1:1.5, 1;0) ，算法精度为 4阶 。解：先编写函数文件verderpol2.m：function xprime=verderpol2(t,x)xprime=x(2); x(1)*sin(t);调用命令：t,x=ode45(verderpol2,0:0.1:1.5, 1;0)；3. 设 用数值解法算出y(1)= 0.2714 （精确到4位小数）, 你用的方法是 Runga-Kutta 方法 ，调用的 Matlab命令是 t,x=ode45(verderpol3, 0:0.1:1.5, 1; 0) ，算法精度为 4-5阶 。解：先编写函数文件verderpol3.m：function xprime=verderpol3(t,x)xprime=x(2); x(2)*sin(t)-x(1)*exp(t);调用命令：t,x=ode45(verderpol3, 0:0.1:1.5, 1; 0) 4. 用电压V=14伏的电池给电容器充电，电容器上t时刻的电压满足：，其中是电容器的初始电压，是充电常数。试用下列数据确定和。t（秒）0.30.51.02.04.07.0v(t)5.68736.14347.16338.862611.032812.6962 你用的方法是 最小二乘法 ，结果是= 4.971112 ，= 3.586875 。程序answer4：t=0.3,0.5,1.0,2.0,4.0,7.0;v=5.6873,6.1434,7.1633,8.8626,11.0328,12.6962;fun=inline(14-(14-k(1)*exp(-t/k(2),k,t);k0=1,1;k = lsqcurvefit(fun,k0, t, v);disp(V0=,mat2str(k(1) )disp(=,mat2str(k(2) )5. 小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力。当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）。重力加速度取9.8米/秒2.A. 建立火箭升空过程的数学模型（微分方程）；B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度。解：火箭上升可以分为两个过程：1、 有燃料产生推力的过程；2、燃料耗尽之后上升的过程。第一个过程：持续的时间为： 向上加速过程火箭的质量为： (1) t40s空气阻力与速度平方成正比，即，v为火箭的速度，方向是竖直向下根据牛顿第二定律得 a为加速度；在加速过程有： (2) 有 (3 ) 和 (4) y为高度；联立（1）、（2）、（3）和（4）带入数字得微分方程： :初始条件为：y(0)=0 ; y(0)=0第二个过程：t40火箭只受到重力和阻力的作用由于燃料已经用完,则剩余质量为：m=900-600=300kg对物体进行受力分析得： (5)联立（3）（4）（5） 式得： (6)代入数值得微分方程 ： t40 （7） 初始条件由第一个过程的终值给出。 运行answer5.m 结果：引擎关闭瞬间火箭的高度:8322.96171214975m引擎关闭瞬间火箭的速度:258.982232154038m/s引擎关闭前瞬间火箭的速度:0.770937904682857m/s2火箭到达的最高点高度:9191.95731955837m火箭到达最高点的时间:51s引擎关闭后瞬间火箭的高度:-99.2413333333333m/s2程序先编写两个函数如下：verderpol5_1.mfunction xprime=verderpol5(t,x)xprime=x(2); (30000-0.4*x(2)2-(900-15*t)*9.8)/(900-15*t);verderpol5_2.mfunction xprime=verderpol5_2(t,x)xprime=x(2); (-0.4*x(2)2)/300-9.8;answer5.m如下：t,x=ode45(verderpol5_1, 0:10:40, 0; 0);a_before=(30000-0.4*x(10)2-(900-15*t(5)*9.8)/(900-15*t(5);y=x(5);v=x(10);disp(引擎关闭瞬间火箭的高度:,mat2str(y) m)disp(引擎关闭瞬间火箭的速度:,mat2str(v) m/s)disp(引擎关闭前瞬间火箭的速度:,mat2str(a_before) m/s2)t,x=ode45(verderpol5_2, 40:1:60, 8323; 259);h=max(x);disp(火箭到达的最高点高度:,mat2str(h(1) m)for k=1:20 if x(k)=h(1); t_max=t(k); break; endenddisp(火箭到达最高点的时间:,mat2str(t_max) s);a_after=(-0.4*h(2)2)/300-9.8; disp(引擎关闭后瞬间火箭的高度:,mat2str(a_after) m/s2);6. 冰淇淋的下部为椎体，上部为半球。设它由锥面和球面围成，用蒙特卡罗方法计算它的体积。解：方程是一个以（0，0，1）为球心，半径为1，椎体是与半球切面相接的高度为1的椎体，椎体地面和半球地面面积相等。所求的椎体和球面围成的体积包含在球体里面，而且椎体的顶点刚好落在球面上。设计一个与球相切的正方体，边长为2，正方体的体积为8. 在正方体内随机的投点，分别判断点是否落在半圆内或者落在椎体内。则点落在圆锥和球面围城的体积里面的概率应该等于圆锥和球面围城的体积和正方体的体积之比。假说总投点数为n,落在圆锥和球面围城的体积内的点数为M，则应有： V/ V正=m/n 在matlab里面编写程序answe6.m运行得：冰淇淋的体积:3.144992 ，程序answe6.m如下：clear;n=100000;m=0;for i=1:n x=rand(1)*2-1;y=rand(1)*2-1;z=rand(1)*2; if(z=1) if(x2+y2+(z-1)2=1) m=m+1; end else if(x2+y2)0.5=z) m=m+1; end endendV=8*m/n;disp(冰激淋的体积:,mat2str(V);7. 容器盛满水后，低端直径为的小孔开启。根据水力学知识，当水面高度为时，水从小孔中流出的速度（为重力加速度，0.6为孔口收缩系数）。若容器为倒圆锥形，现测得容器高和上底面直径均为1.2m，小孔直径为3cm，问水从小孔中流完需要多少时间？2分钟时水面高度是多少？解：射容器中总共的水的体积为：V,经过t时间后，从容器中流出的水的体积为Vout,剩余水的体 （1）, 而： （2） ;设剩余水体高度为h 则： （3）将（2）、（3）代入（1）得 （4）对式（4）两边求导得： （5） 即： （6） 代入数值进一步化简得： (7)在matlab里面编写程序answe7.m运行得：水从小孔中流完需要时间:264s2分钟时水面高度是:0.944182010688867m可知水从小孔中流完需要时间:264s；2分钟时水面高度是:0.9442m程序answe7.m如下：fun=inline(-(5.4*10-4)*sqrt(2*9.8*h)/(h2),t,h);t,h=ode23(fun,0:1:300,1.2);a=min(h);for k=1:300 if h(k)=a t_over=t(k); endenddisp(水从小孔中流完需要时间:,mat2str(t_over) s);h_120=h(120);disp(2分钟时什么的高度是:,mat2str(h_120) m);8. 某公司投资2000万元建设一条生产线。投产后，在t时刻（单位：年）所需的原材料、设备维护、人员经费等所带来的追加成本为（百万元），而产品销售所带来的收入为（百万元），试确定该生产线在何时停产可获得最大利润？最大利润是多少？（注：在t时刻的利润 = 销售收入-追加成本）解：因为在t时刻的利润 = 销售收入-追加成本所以构造在t时刻的利润函数f(t)=H(t)-G(t)=13-t-3t2/3; 令t1/3=x,则f(x)=-x3-3x2+13 令f(x)=-x3-3x2+13=0,求解得到t的三个极值。 可得矩阵P=-1,-3,0,13求最佳生产时间的源程序如下：p=-1,-3,0,13;x=roots(p);t=x.3 再分别以t三个值为上限对函数f(t)求积分，比较大小后，得到最大利润与最佳生产时间。求最大利润的程序代码如下： x=0:0.01:t（1）;y=13-x-3*x.(2/3);trapz(x,y) x=0:0.01:t（2）;y=13-x-3*x.(2/3);trapz(x,y) x=0:0.01:t（3）;y=13-x-3*x.(2/3);trapz(x,y)运行程序answer8.m结果：最佳生产时间:4.64647020645183年最大利润:6.32077370357982百万元可知最佳生产时间t=4.6465年，可获得的最大利润 6.3208百万元。程序answer8.m：clear;p=-1,-3,0,13;x=roots(p);t=x.3;x=0:0.01:t(1);y0=13-x-3*x.(2/3);a=trapz(x,y0); x=0:0.01:t(2);y1=13-x-3*x.(2/3);b=trapz(x,y1); x=0:0.01:t(3);y2=13-x-3*x.(2/3);c=trapz(x,y2);d=a,b,c; sum=max(d);maxsum=sum-20;if sum=a t=t(1);elseif sum=b t=t(2);elseif sum=c t=t(3);end disp(最佳生产时间:,mat2str(t) 年) disp(最大利润:,mat2str(maxsum) 百万元)9. 编写能实现如下功能的主程序与子函数文件。屏幕提示从键盘输入一个四位整数，输入的整数按如下规则加密后输出。加密规则：每位数字都加上7，然后除以10所得的余数取代该数字。解：设四位数的y千,百，十,个位上数字分别为a ,b,c,d,即y=a*1000+b*100+c*10+d每位数都加7之后为a+7 ,b+7,c+7,d+7，然后四个数字分别除以10之后去余数，再取代原来的相应位置上的数字在matl