资源描述：
《一道数学建模问题的matlab求解方法》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、2006年2月安庆师范学院学报(自然科学版)Feb.2006第12卷第1期JournalofAnqingTeachersCollege(NaturalScienceEdition)Vol.12No.1一道数学建模问题的Matlab求解方法朱夜明1,乔宗敏2(1.安徽大学学报编辑部,安徽合肥230039;2.安徽教育学院数学系,安徽合肥230061)摘要:借助于Matlab工具箱中四阶龙格!库塔法求解微分方程的工具,对一类数学建模问题进行了研究。通过对渡河问题的解析解和数值解进行图像和理论分析,说明了四
2、阶龙格!库塔法对求解此类数学建模问题的有效性和实用性。关键词:数学建模;Matlab;龙格!库塔法中图分类号:O241文献标识码:A文章编号:1007-4260(2006)01-0042-04Matlab是英文MatrixLaboratory(矩阵实验室)的缩写[1,2],由美国Mathworks公司开发。Matlab语言是当今国际上科学界(尤其是自动控制领域)最具影响力、也是最有活力的科学计算软件。它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、便捷的与其他程序和语言接口的
3、功能。Matlab语言在各国高校与研究单位正扮演着重要的角色。在美国的一些大学里,Matlab正在成为对数值线性代数以及其他一些高等应用数学课程进行辅助教学的有益工具。在工程技术界,Matlab也被用来解决一些实际课题和数学模型问题。而将Matlab应用于数学教学,不仅可以提高学生的理解和解题能力,而且可以激发学生的学习兴趣[3]。Matlab数学工具箱中集成的龙格-库塔法是求解微分方程的一种有力工具,下面将通过一个具体数学模型来探讨其应用。1问题的提出一只小船渡过宽为d的河流(见图1),目标是起点A
4、正对着的另一岸B点。已知河水流速v1与船在静水中的速度v2之比为k。a)建立小船航线的方程,求其解析解;b)设d=100m,v1=1m/s,v2=2m/s,用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线,作图并与解析解比较;c)若流速v1为0,0.5,1.5,2(m/s),试分析结果将如何?2模型建立如图1,以B为原点,沿河岸向右为x轴正向,垂直河岸向下为y轴正向,建立坐标系。设在t时刻,船在x方向上的位移是x(t),在y方向上的位移是y(t),则:在t时刻,船在x方向上的速度是x′(t),在
5、y方向上的速度是y′(t),将船的速度v和静水速度v1在x,y轴方向上分解,可得ìvx=v1-v2sinqíîvy=-v2cosq图1模型分析图收稿日期:2005-11-23基金项目:安徽省教育厅青年科研基金项目(2005jql1153)作者简介:朱夜明(1970-),女,安徽庐江人,安徽大学学报编辑部副编审。第1期朱夜明1,乔宗敏2:一道数学建模问题的Matlab求解方法·43·x又因为船头始终指向B点,所以tanθ=,所以yìxìx¢(t)=v-vxïvx=v1-v2ï1222x2+y2ïx+yï
6、即ííyïv=-vyïy¢(t)=-v2ïy222ïx2+y2îx+yî这就是本题的微分方程。初始条件为:x(0)=0,y(0)=-d.3模型求解3.1解析解ìx=r×cosq令í,将直角坐标化成极坐标,则原微分方程化成îy=r×sinqìr¢×cosq-r×q¢×sinq=v1-v2×cosqìr¢=v1×cosq-v2í,化简后为ír¢×sinq+r×q¢×cosq=-v×sinqr×q¢=-v×sinqî2î1r¢v1两式相除得=-+2××¢,(cotq)qrvsinq1Rqtan()v于是解得
7、2,其中R=2r=d×vsinq13.2数值解法根据上述微分方程,用龙格-库塔方法求解,源程x0=[0,-d]';%初始条件序如下:[t,x]=ode45(@fun,[0,1000],x0,[],V1,V2);%fun.mplot(x(:,1),x(:,2),'r');holdon;%此函数是微分方程组%作图functionXdot=fun(t,x,V1,V2)[t,x(:,1),x(:,2)]%x(1)代表x;x(2)代表y%打印数据%加入限制条件,防止无限循环%下面是作出精确解的图像d=100,V
8、1=1,V2=2;seta=linspace(-pi/2,0,100);if(norm(x)>1e-5)d=100,V1=1,V2=2;Xdot=[V1-V2*x(1)/sqrt(x(1)^2+x(2)^2),-rou=d*(abs(tan((seta/2))).^(V2/V1))./sin(seta);V2*x(2)/sqrt(x(1)^2+x(2)^2)]';xp=-rou.*cos(seta);yp=-rou.*sin(seta);elseplot
