资源描述:
《matlab微分方程的求解》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库。
1、Matlab程序设计:微分方程求解主讲人:王佐才1.引言Matlab能够求解的微分方程包括:常微分方程的初值问题,常微分方程的边值问题,时变常微分方程的初值问题,以及偏微分方程。2.常微分方程的初值问题Matlab可以求解的常微分方程包括:显示常微分方程:y'=f(t,y)线性隐式常微分方程:M(t,y)y'=f(t,y),其中,M(t,y)为矩阵。全隐式常微分方程:f(t,y,y')=03.利用Matlab编程时需要用的主要命令ode45:基于显示Runge-Kutta(4,5)方法求解。对于多数方程来讲,ode45是最佳的尝试函数命令。ode23:基于显示Runge-Kutta(2,3)
2、方法求解。ode113:利用变阶Adams-Bashforth-Moulton算法求解。与ode45函数相比,该方程对于精密度步长及方程难于估计时效更好,但是该方法是多步算法,需要用到前面几个节点的信息来求解当前节点,效率较低。4.ode45调用格式[t,x]=ode45(@myode,[t0:dt:t1],initial_condition)从调用格式来看,首先必须生成函数文件“myode”。functiondydt=myode(t,y)dydt=[fun1;fun2;…;funn];end5.应用实例一:求解微分方程:(1+y2)y'-2y=0,初值y(0)=1时,t从0至10时刻的解。
3、functiondydt=myode(t,y)dydt=[y(1)*2./(1+y(1).^2)];end[t,y]=ode45(@myode,[0:0.01:10],[1])6.应用实例二:3求解范德蒙方程y''-1-y2y'+y=0的解。时间0至20,初值y(0)=2,y’(0)=0。functiondydt=myode(t,y)dydt=[y(2);(1-y(1)^2)*y(2)-y(1)];end[t,y]=ode45(@myode,[0:0.01:20],[2;0])1.应用实例三:一质点在空中飞行,所受的空气阻力方向始终与速度方向相反,大小与速度平方成正比。求质点的飞行轨迹和飞行
4、路程。根据题意:我们可以建立如下方程(vx,vy分别为质点沿x、y方向的速度,c为空气阻力系数):dx/dt=vxdy/dt=vymdvx/dt=-cv^2cosa=-cv*vxmdvy/dt=-cv*vy-mgfunctiondydt=myode(t,x)c=0.01;m=1;g=9.8;dydt=[x(3);x(4);-1/m.*(c.*sqrt((x(3).^2+x(4).^2))).*x(3);-1/m.*(c.*sqrt((x(3).^2+x(4).^2))).*x(4)-g];end[t,y]=ode45(@myode,[0:0.01:6],[0;0;50*cos(pi/6);5
5、0*sin(pi/6)];2.应用实例四:3如图所示的弹簧-阻尼-质量系统的质量为m,刚度为k,阻尼为c。其自由振动的微分方程可以表示为:mx+cx+kx=0,其中x为任一时刻质量块离平衡位置的位移,x与x分别为位移对时间的一阶、二阶导数,分别称为速度与加速度。若m=10,k=100,c=4。编写Matlab程序求该质量块在初始位移为1,初始速度为0时,前20秒的位移解。若c=8,50时,从新求解上述方程,并画出c=4,8,50时的位移曲线图(提示,先编写函数文件,再调用matlab命令求解)。functiondydt=myode(t,x)m=10;k=100;c=4;dydt=[x(2);
6、-c/m*x(2)-k/m*x(1)];end[t,y]=ode45(@myode,[0:0.01:20],[1,0]);1.作业:任意找一具有物理背景的二阶或者高阶微分方程或方程组,在给定初始条件下,编写matlab程序求解该微分方程或者方程组。3