常微分方程的数值解教程文件.doc
《常微分方程的数值解教程文件.doc》由会员分享,可在线阅读,更多相关《常微分方程的数值解教程文件.doc(26页珍藏版)》请在咨信网上搜索。
1、常微分方程的数值解精品文档实验4 常微分方程的数值解【实验目的】1掌握用MATLAB软件求微分方程初值问题数值解的方法;2通过实例用微分方程模型解决简化的实际问题;3了解欧拉方法和龙格-库塔方法的基本思想和计算公式,及稳定性等概念。【实验内容】题3小型火箭初始重量为1400kg,其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s,由此产生32000N的推力,火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方,比例系数为0.4kg/m,求引擎关闭瞬间火箭的高度、速度、加速度,及火箭到达最高点的时的高度和加速度,并画出高度、速度、加速度随时间变化的图形。模型及其求解火
2、箭在上升的过程可分为两个阶段,在全过程中假设重力加速度始终保持不变,g=9.8m/s2。在第一个过程中,火箭通过燃烧燃料产生向上的推力,同时它还受到自身重力(包括自重和该时刻剩余燃料的重量)以及与速度平方成正比的空气阻力的作用,根据牛顿第二定律,三个力的合力产生加速度,方向竖直向上。因此有如下二式:a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8dh/dt=v又知初始时刻t=0,v=0,h=0。记x(1)=h,x(2)=v,根据MATLAB可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下:function dx = r
3、ocket( t,x )a=(32000-0.4*x(2)2)/(1400-18*t)-9.8;dx=x(2);a;endts=0:1:60;x0=0,0;t,x=ode45(rocket,ts,x0);h=x(:,1);v=x(:,2);a=(32000-0.4*(v.2)./(1400-18*t)-9.8;t,h,v,a;数据如下:thva00013.061.006.5713.1913.302.0026.4426.5813.453.0059.7640.0613.504.00106.5753.5413.435.00166.7966.8913.266.00240.2780.0212.997.0
4、0326.7292.8312.618.00425.79105.2212.159.00536.99117.1111.6210.00659.80128.4311.0211.00793.63139.1410.3812.00937.85149.189.7113.001091.79158.559.0214.001254.71167.238.3315.001425.93175.227.6516.001604.83182.556.9917.001790.78189.226.3618.001983.13195.275.7619.002181.24200.755.2120.002384.47205.704.69
5、21.002592.36210.184.2222.002804.52214.193.7923.003020.56217.793.4124.003240.08221.013.0725.003462.65223.922.7726.003687.88226.562.5027.003915.58228.972.2728.004145.60231.142.0629.004377.76233.111.8930.004611.86234.911.7431.004847.68236.571.6232.005085.02238.141.5133.005323.85239.611.4134.005564.1124
6、0.991.3335.005805.77242.281.2736.006048.72243.501.2137.006292.87244.681.1738.006538.11245.831.1339.006784.48246.961.0940.007031.96248.051.0741.007280.54249.101.0542.007530.19250.121.0343.007780.85251.141.0244.008032.49252.151.0045.008285.12253.160.9946.008538.75254.150.9847.008793.39255.120.9748.009
7、049.01256.070.9749.009305.58257.030.9650.009563.08257.990.9551.009821.52258.950.9452.0010080.93259.900.9353.0010341.30260.830.9354.0010602.62261.750.9455.0010864.86262.670.9456.0011127.98263.610.9357.0011392.04264.540.9158.0011657.03265.460.9159.0011922.96266.350.9260.0012189.78267.260.92因此,在引擎关闭的瞬间
8、,火箭的速度为267.26m/s,高度为12189.78m,加速度为0.92m/s2。(2)在第二个阶段,火箭的重量保持不变,没有向上的推力,只收到重力和空气阻力。因此有如下关系式:a=dv/dt=(-mg-0.4v2)/m=-0.4v2/320-9.8dh/dt=v假设在80秒之前达到最高点,以60秒时的速度、高度、加速度为初始值进行计算,程序如下:function dx = rocket2( t,x )dx=x(2);-0.4*x(2)2/320-9.8;endts2=60:1:80;x1=12189.78,267.26;t2,x2=ode45(rocket2,ts2,x1);h2=x2(
9、:,1);v2=x2(:,2);a2=-0.4*v2.2./320-9.8;t2,h2,v2,a2;数据如下: t2 h2 v2 a2 60.00 12189.78 267.26 -99.08 61.00 12416.32 192.70 -56.22 62.00 12584.73 147.43 -36.97 63.00 12715.67 116.11 -26.65 64.00 12819.59 92.79 -20.56 65.00 12902.81 74.29 -16.70 66.00 12969.22 58.96 -14.15 67.00 13021.42 45.73 -12.41 68.0
10、0 13061.17 33.95 -11.24 69.00 13089.63 23.12 -10.47 70.00 13107.61 12.91 -10.01 71.00 13115.55 3.02 -9.81 72.00 13113.66 -6.80 -9.86 73.00 13101.90 -16.78 -10.15 74.00 13079.96 -27.19 -10.72 75.00 13047.26 -38.34 -11.64 76.00 13002.90 -50.62 -13.00 77.00 12945.47 -64.56 -15.01 78.00 12872.96 -80.96
11、-17.99 79.00 12782.31 -101.08 -22.57 80.00 12668.88 -127.03 -29.97 可以看到:在第60秒瞬间,加速度发生了突变,从0.92m/s2突变为-99.08m/s2;而在第71秒至第72秒之间,速度从正变为负,即速度反向,说明在第71秒中的某个时刻速度为零,火箭达到了最高点。因此需要对这个时间段进行分析,并且找到速度减小到0的时刻和此时的高度。以0.1为步长,在71s到72s中重新求解微分方程的数值解。71.00 13115.16 2.98 -9.8171.10 13115.41 2.00 -9.8171.20 13115.56 1.0
12、2 -9.8071.30 13115.61 0.04 -9.8071.40 13115.57 -0.94 -9.8071.50 13115.42 -1.92 -9.80可见在t=71.3时,速度为0.04,可视为速度为零点,此时最大高度为13115.61,加速度-9.80。综合(1),(2),可以绘出高度,速度和加速度随时间的变化曲线。plot(t,h,t2,h2),xlabel(t/s),ylabel(h/m),title(高度随时间变化曲线);plot(t,v,t2,v2),xlabel(t/s),ylabel(v/(m/s),title(速度随时间变化曲线);aa=a,a2;tt=t,t
13、2;plot(tt,aa),xlabel(t/s),ylabel(a/(m/s2),title(加速度随时间的变化曲线);题6一只小船渡过宽为d的河流,目标是起点A正对着的另一岸B点。已知河水流速v1与船在静水中的速度v2之比为k。(1)建立描述小船航线的数学模型,求其解析解;(2)设d=100m,v1=1m/s,v2=2m/s,用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线,作图,并与解析解比较;(3)若流速v1=0,0.5,1.5,2m/s,结果将如何。模型及其求解(1).假设在航行过程中,人们不知道水流的速度,小船的方向始终指向目标B,因此若以B为原点,我们可以得到如下方程组:解
14、初值为(x, y)=( 0, -d) 的常微分方程组,得到解析解为: 其中 c=1/d=0.01,故有事实上,若用matlab中计算微分方程的语句:x,y=dsolve(Dx=v1-v2*x/sqrt(x2+y2),Dy=-v2*y/sqrt(x2+y2),x(0)=0,y(0)=-d,t);则显示“Warning: Explicit solution could not be found.”即无法得到x,y关于t的分立解。(2)d=100m,v1=1m/s,v2=2m/s。求数值解时,由解析解可以看出,此为刚性方程。为避免用ode45s求解时间过长。采用ode23s进行求解。假定100s可以
15、到达对岸。ts=0:0.1:100;x0=0,-100;t,x=ode23s(boat,ts,x0);t,x;plot(t,x);gtext(x);gtext(y);xlabel(时间t/s);ylabel(距离/m);title(x,y与时间t的关系);可以得到数据如下(部分):t/sx/my/m10.008.93-80.0420.0015.41-60.3630.0018.87-41.4940.0018.71-24.3150.0014.48-10.3450.1014.42-10.2366.500.19-0.0066.600.09-0.0066.700.000可知在t=66.7s时,船到达对岸
16、B点。做 x,y与时间t的关系图:从曲线上可以看出,0到30s这段时间内,y的增长几乎呈线性关系,即小船几乎研直线匀速前进。现在求解析解并将之与数值解对比:function x = jiexi(y,k)x=-0.5*(-0.01)k*y.(k+1)+0.5*(-0.01)(-k).*y.(1-k);endy=-100:1:0;x2=jiexi(y,0.5);plot(x2,y,ro,x(:,1),x(:,2);legend(解析解,数值解,-1);从轨迹曲线中也可以看到,用数值解得到结果与解析解几乎重合,可信度很高。(3)当改变v1的值时,解析式中的k值将发生变化,此时将对结果产生影响。利用M
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 教程 文件
1、咨信平台为文档C2C交易模式,即用户上传的文档直接被用户下载,收益归上传人(含作者)所有;本站仅是提供信息存储空间和展示预览,仅对用户上传内容的表现方式做保护处理,对上载内容不做任何修改或编辑。所展示的作品文档包括内容和图片全部来源于网络用户和作者上传投稿,我们不确定上传用户享有完全著作权,根据《信息网络传播权保护条例》,如果侵犯了您的版权、权益或隐私,请联系我们,核实后会尽快下架及时删除,并可随时和客服了解处理情况,尊重保护知识产权我们共同努力。
2、文档的总页数、文档格式和文档大小以系统显示为准(内容中显示的页数不一定正确),网站客服只以系统显示的页数、文件格式、文档大小作为仲裁依据,平台无法对文档的真实性、完整性、权威性、准确性、专业性及其观点立场做任何保证或承诺,下载前须认真查看,确认无误后再购买,务必慎重购买;若有违法违纪将进行移交司法处理,若涉侵权平台将进行基本处罚并下架。
3、本站所有内容均由用户上传,付费前请自行鉴别,如您付费,意味着您已接受本站规则且自行承担风险,本站不进行额外附加服务,虚拟产品一经售出概不退款(未进行购买下载可退充值款),文档一经付费(服务费)、不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
4、如你看到网页展示的文档有www.zixin.com.cn水印,是因预览和防盗链等技术需要对页面进行转换压缩成图而已,我们并不对上传的文档进行任何编辑或修改,文档下载后都不会有水印标识(原文档上传前个别存留的除外),下载后原文更清晰;试题试卷类文档,如果标题没有明确说明有答案则都视为没有答案,请知晓;PPT和DOC文档可被视为“模板”,允许上传人保留章节、目录结构的情况下删减部份的内容;PDF文档不管是原文档转换或图片扫描而得,本站不作要求视为允许,下载前自行私信或留言给上传者【丰****】。
5、本文档所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用;网站提供的党政主题相关内容(国旗、国徽、党徽--等)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。
6、文档遇到问题,请及时私信或留言给本站上传会员【丰****】,需本站解决可联系【 微信客服】、【 QQ客服】,若有其他问题请点击或扫码反馈【 服务填表】;文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“【 版权申诉】”(推荐),意见反馈和侵权处理邮箱:1219186828@qq.com;也可以拔打客服电话:4008-655-100;投诉/维权电话:4009-655-100。