第二章--常微分方程的数值解法

微分方程的数值解法,四阶龙格库塔法（TheFourth-OrderRungeKuttaMethod）,常微分方程（Ordinarydifferentialequations,ODE),初值问题-给出初始值边值问题-给出边界条件,与初值常微分方程解算有关的指令ode23ode45ode113ode23tode15sode23sode23tb,一.解ODE的基本机理:,2.把高阶方程转换成一阶微分方程组,1.列出微分方程,初始条件,令,(2.1),(2.2),(2.3),例：著名的VanderPol方程,令,降为一阶,初始条件,3.根据式(2.2)编写计算导数的M函数文件-ODE文件,把t,Y作为输入宗量，把作为输出宗量,%Mfunctionfilename:dYdt.mfunctionYd=f(t,Y)Yd=f(t,Y)的展开式,例VanderPol方程,%Mfunctionfilename:dYdt.mfunctionYd=f(t,Y)Yd=zeros(size(Y);,4.使编写好的ODE函数文件和初值供微分方程解算指令（solver）调用,Solver解算指令的使用格式,输出宗量形式,说明：t0：初始时刻；tN：终点时刻Y0：初值；tol：计算精度,例题1：著名的VanderPol方程,%主程序（程序名：VanderPol_ex1.m）t0=0;tN=20;tol=1e-6;Y0=0.25;0.0;t,Y=ode45(dYdt,t0,tN,Y0,tol);subplot(121),plot(t,Y)subplot(122),plot(Y(:,1),Y(:,2),解法1：采用ODE命令,VanderPol方程,%子程序（程序名：dYdt.m）functionYdot=dYdt(t,Y),Ydot=Y(2);-Y(2)*(Y(1)2-1)-Y(1);,或写为,functionYdot=dYdt(t,Y)Ydot=zeros(size(Y);Ydot(1)=Y(2);Ydot(2)=-Y(2)*(Y(1).2-1)-Y(1);,各种solver解算指令的特点,二.四阶Runge-Kutta法,对I=a,b作分割,步长,单步法-Runge-Kutta方法多步法-Admas方法,计算的近似值时只用到，是自开始方法,Runge-Kutta法是常微分方程的一种经典解法MATLAB对应命令：ode45,四阶Runge-Kutta公式,四阶Runge-Kutta法计算流程图,开始,Plot,初始条件：；积分步长：迭代次数：,输出结果,子程序计算,End,三.Runge-Kutta法解VanderPol方程的Matlab程序结构主程序：RK_vanderpol.m子程序：RK_sub.m（函数文件）,解法2：采用Runge_Kutta法编程计算,主程序：RK_vanderpol.mt0=0;tN=20;y0=0.25;0;h=0.001;t=t0:h:tN;N=length(t);j=1;fori=1:Nt1=t0+h;K1=RK_sub(t0,y0);K2=RK_sub(t0+h/2,y0+h*K1/2);K3=RK_sub(t0+h/2,y0+h*K2/2);K4=RK_sub(t0+h,y0+h*K3);y1=y0+(h/6)*(K1+2*K2+2*K3+K4);yy1(j)=y1(1);yy2(j)=y1(2);t0=t1;y0=y1;j=j+1;endsubplot(121),plot(t,yy1,t,yy2);gridsubplot(122),plot(yy2,yy1);grid,子程序：RK_sub.mfunctionydot=vdpol(t,y)ydot=zeros(size(y);ydot(1)=y(2);ydot(2)=-y(2)*(y(1)2-1)-y(1);或写为：ydot=y(1);-y(2)*(y(1)2-1)-y(1);,四.Matlab对应命令：ode23,ode45,调用格式：t,y=ode23(函数文件名,t0,tN,y0,tol)t,y=ode45(函数文件名,t0,tN,y0,tol),默认精度：ode231e-3ode451e-6,说明：t0：初始时刻；tN：终点时刻y0：初值；tol：计算精度,3月15日作业:1.VanderPol方程的两种解法:1)采用ode45命令2)Runge-Kutta方法2.Duffing方程的求解（Runge-Kutta方法，计算步长h=0.001，计算时间t0=0.0,tN=100）要求：写出程序体，打印所绘图形，图形标题用个人的名字。,Duffing方程,五.动力学系统的求解,1.动力学方程,2.二阶方程转成一阶方程,(1),令：,(2),其中：,即：,(2),3.Matlab程序(主程序：ZCX),t0;Y0;h;N;P(t0);%输入初始值、步长、迭代次数、初始激励力；fori=1:Nt1=t0+hP(t0)%输入t0时刻的外部输入K1=ZCX_sub(t0,Y0,P)P(t0+h/2)%输入(t0+h/2)时刻的外部输入K2=ZCX_sub(t0+h/2,Y0+hK1/2,P)K3=ZCX_sub(t0+h/2,Y0+hK2/2,P)P(t0+h)K4=ZCX_sub(t0+h,y0+hK3,P)Y1=y0+(h/6)(K1+2K2+2K3+K4)t1,Y1(输出t1,y1)nexti输出数据或图形,Matlab程序(子程序：ZCX_sub.m),functionydot=f(t,Y,P)M=,K=,C=%输入结构参数A=0,eye(n,n);-M-1K,-M-1Cydot=AY+P,例题2：三自由度质量弹簧系统,矩阵表示,其中：,动力学方程：,解析解：,已知参数：m1=m2=m3=1,k1=2,k2=2,K3=1,K4=2,P0=1,要求：采用四阶龙格库塔法编程计算三个质量的响应时程.计算时间050,例如：,4阶龙格库塔法的结果,ode45的结果,第一个质量的位移响应时程,结果完全一致,MATLAB程序（1）4阶RK方法：（2）采用ode45：m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m,例题3：蹦极跳系统的动态仿真,蹦极者系着一根弹性绳从高处的桥梁（或山崖等）向下跳。在下落的过程中，蹦极者几乎处于失重状态。按照牛顿运动规律，自由下落的物体由下式确定：,其中，m为人体的质量，g为重力加速度，x为物体的位置，第二项和第三项表示空气的阻力。其中位置x的基准为桥梁的基准面（即选择桥梁作为位置的起点x0），低于桥梁的位置为正值，高于桥梁的位置为负值。如果人体系在一个弹性常数为k的弹性绳索上，定义绳索下端的初始位置为0，则其对落体位置的影响为：,空气的阻力,整个蹦极系统的数学模型为：,设桥梁距离地面为50m，即h2=50，蹦极者的起始位置为绳索的长度30m，即h1=30，蹦极者起始速度为0，其余的参数分别为k20，a2a11；m70kg，g10m/s2。,初始条件：,已知参数：,初始条件变为：,cleary0=-30;0;%初始位移和初始速度t,y=ode23tb(m_chap2_ex3_sub,0:0.01:100,y0);x1=50.-y(:,1);%x1代表蹦极者与地面之间的距离plot(t,x1);gridplot(t,y(:,1);grid%y(:,1)代表位移,主程序(程序名：m_chap2_ex3),Matlab程序,functionydot=f(t,y)m=70;k=20;a1=1;a2=1;g=10;x=y(1);%x代表蹦极者的位移x_dot=y(2);%x_dot代表x的速度ifx0ydot=0,1;-k/m,-a1/m-(a2/m)*abs(x_dot)*y+0;g;elseydot=0,1;0,-a1/m-(a2/m)*abs(x_dot)*y+0;g;end,子程序（程序名：m_cha