微分方程的数值解法matlab(四阶龙格—库塔法)教学内容.ppt

《微分方程的数值解法matlab(四阶龙格—库塔法)教学内容.ppt》由会员分享，可在线阅读，更多相关《微分方程的数值解法matlab(四阶龙格—库塔法)教学内容.ppt（36页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 微分方程微分方程(wi fn fn chn)的数值解法的数值解法四阶龙格四阶龙格库塔法库塔法（The Fourth-Order RungeKutta Method）第一页，共36页。常微分方程常微分方程(wi fn fn chn)（Ordinary differential equations,ODE)n初值问题初值问题-给出初始值给出初始值n边值问题边值问题-给出边界条件给出边界条件与初值常微分方程与初值常微分方程(wi fn fn chn)解算有关的指令解算有关的指令ode23 ode45 ode113 ode23tode15s ode23sode23tb第二页，共36页。一一.解解OD

2、E的基本的基本(jbn)机理机理:2.把把高高阶阶方方程程(fngchng)转转换换成成一一阶阶微微分方程分方程(fngchng)组组1.列出微分方程列出微分方程(wi fn fn chn)初始条件初始条件令令(2.1)(2.2)(2.3)第三页，共36页。例：著名例：著名(zhmng)的的Van der Pol方程方程令令 降为一阶降为一阶初始条件初始条件第四页，共36页。3.根据根据(gnj)式式(2.2)编写计算导数的编写计算导数的M函数函数文件文件-ODE文件文件把把t,Y作为作为(zuwi)输入宗量，把输入宗量，把 作为作为(zuwi)输出宗量输出宗量%M function file

3、 name:dYdt.m function Yd=f(t,Y)Yd=f(t,Y)的展开式的展开式例例Van der Pol方程方程(fngchng)%M function file name:dYdt.m function Yd=f(t,Y)Yd=zeros(size(Y);第五页，共36页。4.使编写使编写(binxi)好的好的ODE函数文件和初值函数文件和初值 供微分方程解算指令（供微分方程解算指令（solver）调用）调用Solver解算指令解算指令(zhlng)的的使用格式使用格式t,Y=solver(ODE函数文件名函数文件名,t0,tN,Y0,tol);ode45输出输出(shch

4、)宗宗量形式量形式说明：说明：t0：初始时刻；：初始时刻；tN：终点时刻：终点时刻Y0：初值；：初值；tol：计算精度：计算精度第六页，共36页。例题例题1：著名：著名(zhmng)的的Van der Pol方程方程%主程序主程序(chngx)（程序（程序(chngx)名：名：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：采用：采用(ciyng)ODE命令命令第七页，共

5、36页。Van der Pol方程方程(fngchng)%子程序子程序(chngx)（程序（程序(chngx)名：名：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);第八页，共36页。第九页，共36页。解法指令解题类型特 点适合场合ode45非刚性非刚性采用采用4、5阶阶RungeKutta法法大多数场合的首选算法大多数场合的首选算法ode2

6、3非刚性非刚性采用采用Adams算法算法较低精度（较低精度（103）场合）场合ode113非刚性非刚性多步法；采用多步法；采用Adams算法；高算法；高低精度均可（低精度均可（103106）ode45计算时间太长时计算时间太长时取代取代ode45ode23t适度刚适度刚性性采用梯形法则算法采用梯形法则算法适度刚性适度刚性ode15s刚性刚性多步法；采用多步法；采用2阶阶Rosenbrock算式，精度中等算式，精度中等当当ode45失败时使用；失败时使用；或存在质量矩阵时或存在质量矩阵时ode23s刚性刚性一步法；采用一步法；采用2阶阶Rosenbrock算式，低精度算式，低精度低精度时，比低精

7、度时，比ode15s有有效；或存在质量矩阵时效；或存在质量矩阵时ode23tb刚性刚性采用梯形法则反向数值微分两采用梯形法则反向数值微分两阶段算法，低精度阶段算法，低精度低精度时，比低精度时，比ode15s有有效；或存在质量矩阵时效；或存在质量矩阵时各种各种solver 解算指令解算指令(zhlng)的特点的特点第十页，共36页。二二.四四 阶阶 Runge-Kutta 法法对对 I=a,b作分割作分割(fng)步长第十一页，共36页。初值问题的数值初值问题的数值解法分为两大类解法分为两大类单步法-Runge-Kutta 方法(fngf)多步法-Admas方法(fngf)计算 的近似值 时只用

8、到 ，是自开始方法 uRunge-Kutta法是常微分方程的一种(y zhn)经典解法uMATLAB 对应命令：ode45第十二页，共36页。四阶四阶Runge-Kutta公式公式(gngsh)第十三页，共36页。四四 阶阶 Runge-Kutta 法计算法计算(j sun)流程图流程图开始开始(kish)Next ifor i=1:N Plot初始条件：初始条件：；积分步长：积分步长：迭代次数：迭代次数：输出(shch)结果子程序计算End第十四页，共36页。三三.Runge-Kutta 法解法解Van der Pol 方程的方程的Matlab 程序结构程序结构主程序：主程序：RK_vand

9、erpol.m 子程序：子程序：RK_sub.m（函数（函数(hnsh)文件）文件）第十五页，共36页。解法解法2：采用：采用(ciyng)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=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

10、*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第十六页，共36页。第十七页，共36页。子程序：子程序：RK_sub.m function ydot=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);第十八页，共36页。四四

11、.Matlab对应对应(duyng)命令：命令：ode23,ode45调用调用(dioyng)格式：格式：t,y=ode23(函数文件名函数文件名,t0,tN,y0,tol)t,y=ode45(函数文件名函数文件名,t0,tN,y0,tol)默认默认(mrn)精度：精度：ode231e-3 ode451e-6说明：说明：t0：初始时刻；：初始时刻；tN：终点时刻：终点时刻y0：初值；：初值；tol：计算精度：计算精度第十九页，共36页。3月月15日作业日作业:1.Van der Pol 方程的两种解法方程的两种解法:1)采用采用ode45命令命令 2)Runge-Kutta方法方法2.Duff

12、ing 方程的求解（方程的求解（Runge-Kutta方法，计算步长方法，计算步长h=0.005，计算时间，计算时间t0=0.0,tN=100）要求要求(yoqi)：写出程序体，打印所绘图形，图形标题用：写出程序体，打印所绘图形，图形标题用个人的名字。个人的名字。Duffing 方程方程(fngchng)第二十页，共36页。第二十一页，共36页。五五.动力学系统动力学系统(xtng)的求解的求解1.动力学方程动力学方程(fngchng)2.二阶方程二阶方程(fngchng)转转成一阶方程成一阶方程(fngchng)(1)令：令：(2)第二十二页，共36页。其中其中(qzhng)：即：即：(2)

13、第二十三页，共36页。3.Matlab 程序程序(chngx)(主程序主程序(chngx)：ZCX)t0;Y0;h;N;P0,w;%输入初始输入初始(ch sh)值、步长、迭代次数、初始值、步长、迭代次数、初始(ch sh)激励力；激励力；for i=1:N t1=t0+h P=P0*sin(w*t0);0.0;0.0%输入输入t0时刻的外部激励力时刻的外部激励力 K1=ZCX_sub(t0,Y0,P)P=%输入输入(t0+h/2)时刻的外部激励力时刻的外部激励力 K2=ZCX_sub(t0+h/2,Y0+hK1/2,P)K3=ZCX_sub(t0+h/2,Y0+hK2/2,P)P=%输入输入

14、(t0+h)时刻的外部激励力时刻的外部激励力 K4=ZCX_sub(t0+h,y0+hK3,P)Y1=y0+(h/6)(K1+2K2+2K3 +K4)t1,Y1 (输出输出 t1,y1)next i输出数据或图形输出数据或图形第二十四页，共36页。Matlab 程序程序(chngx)(子程序子程序(chngx)：ZCX_sub.m)function ydot=f(t,Y,P)M=,K=,C=%输入结构输入结构(jigu)参数参数 P1=zeros(3,1);inv（M）*P；A=zeros(0,0),eye(n,n);-M-1K,-M-1C ydot=AY+P1第二十五页，共36页。例题例题2

15、：三自由：三自由度质量度质量(zhling)弹簧弹簧系统系统m1m2m3k1k2k3x1x2x3k4P0sin(wt)第二十六页，共36页。矩阵表示矩阵表示其中其中(qzhng)：第二十七页，共36页。动力学方程动力学方程(fngchng)：解析解：已知参数：已知参数：m1=m2=m3=1,k1=2,k2=2,K3=1,K4=2,P0=1,要求：采用四阶龙格库塔法编程计算三个质量的响应要求：采用四阶龙格库塔法编程计算三个质量的响应(xingyng)时时程程.计算时间计算时间 0 50例如例如(lr)：第二十八页，共36页。4阶龙格库塔法的结果阶龙格库塔法的结果(ji gu)ode45 的结果的

16、结果(ji gu)第一个质量的位移第一个质量的位移(wiy)响应时程响应时程结果完全一致结果完全一致MATLAB程序程序（1）4阶阶RK方法：方法：（2）采用）采用ode45：m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m 第二十九页，共36页。例题例题(lt)3：蹦极跳系统的动态仿真蹦极跳系统的动态仿真蹦蹦极极者者系系着着一一根根弹弹性性绳绳从从高高处处的的桥桥梁梁（或或山山崖崖等等）向向下下跳跳。在在下下落落的的过过程程中中，蹦蹦极极者者几几乎乎(jh)处处于于失失重重状状态态。按按照照牛牛顿顿运运动动规规律，自由下落的物体由下式确定：律，自由下落的物体由下式确定：

17、其其中中，m 为为人人体体(rnt)的的质质量量，g 为为重重力力加加速速度度，x 为为物物体体的的位位置置，第第二二项项和和第第三三项项表表示示空空气气的的阻阻力力。其其中中位位置置 x 的的基基准准为为桥桥梁梁的的基基准准面面（即即选选择择桥桥梁梁作作为为位位置置的的起起点点 x0），低低于于桥桥梁梁的的位位置置为为正正值值，高高于于桥桥梁梁的的位位置置为为负负值值。如如果果人人体体(rnt)系系在在一一个个弹弹性性常常数数为为 k 的的弹弹性性绳绳索索上上，定定义义绳绳索索下下端端的的初初始始位位置置为为 0，则则其其对对落落体体位位置置的影响为：的影响为：地面地面x桥梁基准面桥梁基准面

18、0梯梯子子h2h1空气的阻力空气的阻力第三十页，共36页。整个蹦极整个蹦极(bn j)系统的数学模型系统的数学模型为：为：设设桥桥梁梁距距离离地地面面为为 50 m，即即 h2=50，蹦蹦极极者者的的起起始始位位置置为为绳绳索索(shn su)的的长长度度 30 m，即即 h1=30，蹦蹦极极者者起起始始速速度度为为 0，其其余的参数分别为余的参数分别为 k20，a2a11；m70 kg，g10 m/s2。地面地面x桥梁基准面桥梁基准面0梯梯子子h2h1初始条件：初始条件：已知参数已知参数(cnsh)：第三十一页，共36页。令令：初始条件变为：初始条件变为：第三十二页，共36页。y0=-30;

19、0;%初始位移和初始速度初始位移和初始速度t,y=ode45(bengji_sub,0:0.01:100,y0);x1=50.-y(:,1);%x1代表蹦极代表蹦极(bn j)者与地面之间的距离者与地面之间的距离plot(t,x1);gridplot(t,y(:,1);grid%y(:,1)代表位移代表位移主程序主程序(chngx)(程序程序(chngx)名：名：bengji.m)Matlab程序程序(chngx)第三十三页，共36页。function ydot=f(t,y)m=70;k=20;a1=1;a2=1;g=10;x=y(1);%x代表蹦极者的位移代表蹦极者的位移(wiy)x_dot

20、=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子程序子程序(chngx)（程序（程序(chngx)名：名：bengji_sub.m）第三十四页，共36页。y(:,1)x1结果分析：结果分析：右右上上图图为为蹦蹦极极者者与与地地面面之之间间的的距距离离。从从结结果果可可看看出出，对对于于(duy)体体重重为为 70 kg 的的蹦蹦极极者者，此此系系统统是是不不安安全全的的，因因为为蹦蹦极极者者与与地地面面之之间间的的距距离离出出现现了了负负值值。因因此此，必必须须使使用用弹弹性性系系数数较较大大的的弹弹性性绳绳索索，才才能能保保证蹦极者的安全。证蹦极者的安全。第三十五页，共36页。作业作业(zuy)（书面作业（书面作业(zuy)，写出程序体）：，写出程序体）：（1）三自由度模型仿真）三自由度模型仿真 (自编自编Runge-Kutta 法）法）（2）蹦极跳模型仿真）蹦极跳模型仿真 (解算指令解算指令ode45)第三十六页，共36页。