《常微分方程的数值解法.pptx》由会员分享，可在线阅读，更多相关《常微分方程的数值解法.pptx（35页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、一.解ODE的基本机理:2.把高阶方程转换成一阶微分方程组1.列出微分方程初始条件令(2.1)(2.2)(2.3)第1页/共35页例：著名的Van der Pol方程令 降为一阶初始条件第2页/共35页3.根据式(2.2)编写计算导数的M函数文件-ODE文件把t,Y作为输入宗量，把 作为输出宗量%M function file name:dYdt.m function Yd=f(t,Y)Yd=f(t,Y)的展开式例Van der Pol方程%M function file name:dYdt.m function Yd=f(t,Y)Yd=zeros(size(Y);第3页/共35页4.使编写好
2、的ODE函数文件和初值 供微分方程解算指令（solver）调用Solver解算指令的使用格式t,Y=solver(ODE函数文件名,t0,tN,Y0,tol);ode45输出宗量形式说明：t0：初始时刻；tN：终点时刻Y0：初值；tol：计算精度第4页/共35页例题1：著名的Van der Pol方程%主程序 （程序名：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：采用OD
3、E命令第5页/共35页Van der Pol方程%子程序 （程序名：dYdt.m）function Ydot=dYdt(t,Y)Ydot=Y(2);-Y(2)*(Y(1)2-1)-Y(1);或写为function Ydot=dYdt(t,Y)Ydot=zeros(size(Y);Ydot(1)=Y(2);Ydot(2)=-Y(2)*(Y(1).2-1)-Y(1);第6页/共35页第7页/共35页解法指令解题类型特 点适合场合ode45非刚性非刚性采用采用4、5阶阶RungeKutta法法大多数场合的首选算法大多数场合的首选算法ode23非刚性非刚性采用采用Adams算法算法较低精度（较低精度（
4、103）场合）场合ode113非刚性非刚性多步法；采用多步法；采用Adams算法；高算法；高低精度均可（低精度均可（103106）ode45计算时间太长时计算时间太长时取代取代ode45ode23t适度刚适度刚性性采用梯形法则算法采用梯形法则算法适度刚性适度刚性ode15s刚性刚性多步法；采用多步法；采用2阶阶Rosenbrock算式，精度中等算式，精度中等当当ode45失败时使用；失败时使用；或存在质量矩阵时或存在质量矩阵时ode23s刚性刚性一步法；采用一步法；采用2阶阶Rosenbrock算式，低精度算式，低精度低精度时，比低精度时，比ode15s有有效；或存在质量矩阵时效；或存在质量矩
5、阵时ode23tb刚性刚性采用梯形法则反向数值微分采用梯形法则反向数值微分两阶段算法，低精度两阶段算法，低精度低精度时，比低精度时，比ode15s有有效；或存在质量矩阵时效；或存在质量矩阵时各种solver 解算指令的特点第8页/共35页二.四 阶 Runge-Kutta 法对 I=a,b作分割步长第9页/共35页初值问题的数值解法分为两大类单步法-Runge-Kutta 方法多步法-Admas方法计算 的近似值 时只用到 ，是自开始方法 uRunge-Kutta法是常微分方程的一种经典解法uMATLAB 对应命令：ode45第10页/共35页四阶Runge-Kutta公式第11页/共35页四
6、 阶 Runge-Kutta 法计算流程图开始Next ifor i=1:N Plot初始条件：；积分步长：迭代次数：输出结果子程序计算End第12页/共35页三.Runge-Kutta 法解Van der Pol 方程的Matlab 程序结构主程序：RK_vanderpol.m 子程序：RK_sub.m（函数文件）第13页/共35页解法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;for i=1:N t1=t0+h;K1=RK_sub(t0,y0);K2
7、=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第14页/共35页第15页/共35页子程序：RK_sub.m function ydot=vdpol(t,y)ydot=zeros(size(y);ydot
8、(1)=y(2);ydot(2)=-y(2)*(y(1)2-1)-y(1);或写为：ydot=y(1);-y(2)*(y(1)2-1)-y(1);第16页/共35页四.Matlab对应命令：ode23,ode45调用格式：t,y=ode23(函数文件名,t0,tN,y0,tol)t,y=ode45(函数文件名,t0,tN,y0,tol)默认精度：ode231e-3 ode451e-6说明：t0：初始时刻；tN：终点时刻y0：初值；tol：计算精度第17页/共35页3月15日作业:1.Van der Pol 方程的两种解法:1)采用ode45命令 2)Runge-Kutta方法2.Duffing
9、 方程的求解（Runge-Kutta方法，计算步长h=0.005，计算时间t0=0.0,tN=100）要求：写出程序体，打印所绘图形，图形标题用个人的名字。Duffing 方程第18页/共35页第19页/共35页五.动力学系统的求解1.动力学方程2.二阶方程转成一阶方程(1)令：(2)第20页/共35页其中：即：(2)第21页/共35页3.Matlab 程序(主程序：ZCX)t0;Y0;h;N;P0,w;%输入初始值、步长、迭代次数、初始激励力；for i=1:N t1=t0+h P=P0*sin(w*t0);0.0;0.0%输入t0时刻的外部激励力 K1=ZCX_sub(t0,Y0,P)P=
10、%输入(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)next i输出数据或图形第22页/共35页Matlab 程序(子程序：ZCX_sub.m)function ydot=f(t,Y,P)M=,K=,C=%输入结构参数 P1=zeros(3,1);inv（M）*P；A=zeros(0,0),eye(n,n);-M-1K
11、,-M-1C ydot=AY+P1第23页/共35页例题2：三自由度质量弹簧系统m1m2m3k1k2k3x1x2x3k4P0sin(wt)第24页/共35页矩阵表示其中：第25页/共35页动力学方程：解析解：已知参数：m1=m2=m3=1,k1=2,k2=2,K3=1,K4=2,P0=1,要求：采用四阶龙格库塔法编程计算三个质量的响应时程.计算时间 0 50例如：第26页/共35页4阶龙格库塔法的结果ode45 的结果第一个质量的位移响应时程结果完全一致MATLAB程序（1）4阶RK方法：（2）采用ode45：m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m 第27页/共
12、35页例题3：蹦极跳系统的动态仿真蹦极者系着一根弹性绳从高处的桥梁（或山崖等）向下跳。在下落的过程中，蹦极者几乎处于失重状态。按照牛顿运动规律，自由下落的物体由下式确定：其中，m 为人体的质量，g 为重力加速度，x 为物体的位置，第二项和第三项表示空气的阻力。其中位置 x 的基准为桥梁的基准面（即选择桥梁作为位置的起点 x0），低于桥梁的位置为正值，高于桥梁的位置为负值。如果人体系在一个弹性常数为 k 的弹性绳索上，定义绳索下端的初始位置为 0，则其对落体位置的影响为：地面x桥梁基准面0梯子h2h1空气的阻力第28页/共35页整个蹦极系统的数学模型为：设桥梁距离地面为 50 m，即 h2=50
13、，蹦极者的起始位置为绳索的长度 30 m，即 h1=30，蹦极者起始速度为 0，其余的参数分别为 k20，a2a11；m70 kg，g10 m/s2。地面x桥梁基准面0梯子h2h1初始条件：已知参数：第29页/共35页令：初始条件变为：第30页/共35页y0=-30;0;%初始位移和初始速度t,y=ode45(bengji_sub,0:0.01:100,y0);x1=50.-y(:,1);%x1代表蹦极者与地面之间的距离plot(t,x1);gridplot(t,y(:,1);grid%y(:,1)代表位移主程序 (程序名：bengji.m)Matlab程序第31页/共35页function
14、ydot=f(t,y)m=70;k=20;a1=1;a2=1;g=10;x=y(1);%x代表蹦极者的位移x_dot=y(2);%x_dot 代表 x 的速度if x0 ydot=0,1;-k/m,-a1/m-(a2/m)*abs(x_dot)*y+0;g;else ydot=0,1;0,-a1/m-(a2/m)*abs(x_dot)*y+0;g;end子程序 （程序名：bengji_sub.m）第32页/共35页y(:,1)x1结果分析：右上图为蹦极者与地面之间的距离。从结果可看出，对于体重为 70 kg 的蹦极者，此系统是不安全的，因为蹦极者与地面之间的距离出现了负值。因此，必须使用弹性系数较大的弹性绳索，才能保证蹦极者的安全。第33页/共35页作业（书面作业，写出程序体）：（1）三自由度模型仿真 (自编Runge-Kutta 法）（2）蹦极跳模型仿真 (解算指令ode45)第34页/共35页感谢您的观看！第35页/共35页