微分方程的数值解法matlab(四阶龙格—库塔法)

1、 微分方程的数值解法微分方程的数值解法四阶龙格四阶龙格库塔法库塔法（The Fourth-Order RungeKutta Method）常微分方程常微分方程（Ordinary differential equations, ODE)n初值问题初值问题-给出初始值给出初始值n边值问题边值问题-给出边界条件给出边界条件与初值常微分方程解算有关的指令与初值常微分方程解算有关的指令ode23 ode45 ode113 ode23tode15s ode23sode23tb一一.解解ODE的基本机理的基本机理:2. 把高阶方程转换成一阶微分方程组把高阶方程转换成一阶微分方程组1. 列出微分方程列出微分方

2、程初始条件初始条件令令)0()0(210yyYyyyy21,),(),(),(2121YtfYtfYtfyyY(2.1) (2.2)(2.3),(tyyfy 00)0(,)0(yyyy例：著名的例：著名的Van der Pol方程方程0)1(2yyyy 令令 yyyy21,降为一阶降为一阶1221221) 1(yyyyyyY21yyY初始条件初始条件2010210)0()0(yyyyY3. 根据式根据式(2.2)编写计算导数的编写计算导数的M函数文件函数文件-ODE文件文件把把t,Y作为输入宗量，把作为输入宗量，把 作为输出宗量作为输出宗量Y %M function file name: dY

3、dt.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););1 () 2 (*) 12).1 () 2 ();2 () 1 (YYYYdYYd4. 使编写好的使编写好的ODE函数文件和初值函数文件和初值 供微分供微分方程解算指令（方程解算指令（solver）调用）调用0YSolver解算指令的使用格式解算指令的使用格式t, Y=solver (ODE函数文件名函数文件名, t0,

4、 tN, Y0, tol);ode45)()()()()()(222112110201tytytytytytyY输出宗量形式输出宗量形式)() 1(:,1tyY)()2(:,2tyY说明：说明：t0：初始时刻；：初始时刻；tN：终点时刻：终点时刻Y0：初值；：初值； tol：计算精度：计算精度例题例题1：著名的著名的Van der Pol方程方程0)1(2yyyy % 主程序主程序 （程序名：程序名：VanderPol _ex1.m） t0 = 0; tN = 20; tol = 1e-6; Y0 = 0.25; 0.0; t, Y=ode45 (dYdt, t0, tN, Y0, tol);

5、 subplot (121), plot (t, Y) subplot (122), plot (Y( :, 1), Y( :, 2)解法解法1：采用：采用ODE命令命令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);解法指令解题类型特 点适合场合ode4

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

7、23s刚性刚性一步法；采用一步法；采用2阶阶Rosenbrock算式，低精度算式，低精度低精度时，比低精度时，比ode15s有有效；或存在质量矩阵时效；或存在质量矩阵时ode23tb刚性刚性采用梯形法则反向数值微分采用梯形法则反向数值微分两阶段算法，低精度两阶段算法，低精度低精度时，比低精度时，比ode15s有有效；或存在质量矩阵时效；或存在质量矩阵时各种各种solver 解算指令的特点解算指令的特点二二. 四四 阶阶 Runge-Kutta 法法00)(),(ytybtaytfdtdybtttaN110对对 I=a,b作分割作分割1, 1 , 0,1Nihhhiii步长初值问题的数值初值问题

8、的数值解法分为两大类解法分为两大类单步法-Runge-Kutta 方法多步法-Admas方法计算 的近似值 时只用到 ，是自开始方法 )(1ntynnyt ,1nyuRunge-Kutta法是常微分方程的一种经典解法常微分方程的一种经典解法uMATLAB 对应命令：对应命令：ode45四阶四阶Runge-Kutta公式公式),()2,21()2,21(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn四四 阶阶 Runge-Kutta 法计算流程图法计算流程图开始开始Next ifor i = 1 : N Plot)22

9、(643211KKKKhyynn11,nnnnyytt初始条件：初始条件： ； 积分步长：积分步长： 迭代次数：迭代次数：0t0yhN输出结果httnn10ttn子程序计算iKEnd三三. Runge-Kutta 法解法解Van der Pol 方程的方程的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;

10、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*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, yy

11、2); gridsubplot (122), plot (yy2, yy1); grid子程序：子程序：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);四四. Matlab对应命令：对应命令：ode23,ode45调用格式调用格式： t, y=ode23 (函数文件名函数文件名, t0, tN, y0, tol) t, y=ode45 (函数文件名

12、函数文件名, t0, tN, y0, tol)默认精度：默认精度： ode231e-3 ode451e-6说明：说明：t0：初始时刻；：初始时刻；tN：终点时刻：终点时刻y0：初值；：初值； tol：计算精度：计算精度3月月15日作业日作业: 1.Van der Pol 方程的两种解法方程的两种解法:1)采用采用ode45命令命令 2)Runge-Kutta方法方法2.Duffing 方程的求解方程的求解（Runge-Kutta方法，计算步长方法，计算步长h=0.005，计算时间，计算时间t0=0.0,tN=100）要求：要求：写出程序体，打印所绘图形，图形标题用个写出程序体，打印所绘图形，图

13、形标题用个人的名字。人的名字。)2.1cos(28.03.03tyyyy Duffing 方程方程0.0)0(,01.0)0(yy五五. 动力学系统的求解动力学系统的求解1. 动力学方程动力学方程2. 二阶方程转成一阶方程二阶方程转成一阶方程00)0(,)0()(),()()()(XXXXXPKXXCXM nRttttt(1) 0 () 0 () 0 (,)(,)()()(2XXYYXXYnRtttt令：令：)()()(tttPAYY(2)12112211)()(,nnnnnnnnRttRPM0PCMKMI0A其中：其中：)()()()()(1111tttttnnnnnPM0XXCMKMI0X

14、X 即：即：00)0(,)0()(),()()()(XXXXXPKXXCXM nRttttt)()()(tttPAYY(2)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= %输入输入 (t0+h/2) 时刻的外部激励力时刻的外部激励力 K2 = ZCX_sub (t0 +

15、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输出数据或图形输出数据或图形)()()(tttPAYYMatlab 程序程序(子程序：子程序：ZCX_sub.m)function ydot = f (t, Y,P) M=, K=, C= %输入结构参数输入结构参

16、数 P1=zeros(3,1);inv（M）*P； A=zeros(0,0), eye(n,n) ; -M-1K, -M-1C ydot =AY+P1)()()(tttPAYY例题例题2：三自由：三自由度质量弹簧系统度质量弹簧系统m1m2m3k1k2k3x1x2x3k4P0sin(wt)0)sin()(02121111tPxxkxkxm 0)()(32321222xxkxxkxm 0)(3432333xkxxkxm 00)sin(0000000003214333322221321321tPxxxkkkkkkkkkkxxxmmm 矩阵表示矩阵表示其中：其中：321xxxX321000000mmm

17、M433332222100kkkkkkkkkkK00)sin()(0tPtP动力学方程：动力学方程：)sin(22221. 0)sin(30363. 2)sin(121088. 0)()()(000321twkPtwkPtwkPtxtxtx解析解：已知参数已知参数：m1=m2=m3=1, k1=2, k2=2, K3=1, K4=2, P0=1, 要求：要求：采用四阶龙格库塔法编程计算三个质量的响应时程.计算时间 0 50例如：例如：)()()(tttPKXXM 305101520253035404550-0.4-0.3-0.2-0.100.10.20.30.40510152025303540

18、4550-0.4-0.3-0.2-0.100.10.20.30.44阶龙格库塔法的结果阶龙格库塔法的结果ode45 的结果的结果第一个质量的位移响应时程第一个质量的位移响应时程结果完全一致结果完全一致MATLAB程序程序（1）4阶阶RK方法：方法： （2）采用）采用ode45： m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m 例题例题3： 蹦极跳系统的动态仿真蹦极跳系统的动态仿真蹦极者系着一根弹性绳从高处的桥梁（或山崖等）向下跳。在蹦极者系着一根弹性绳从高处的桥梁（或山崖等）向下跳。在下落的过程中，蹦极者几乎处于失重状态。按照牛顿运动规律，下落的过程中，蹦极者几乎处于失

19、重状态。按照牛顿运动规律，自由下落的物体由下式确定自由下落的物体由下式确定：xxaxamgxm 210 , 00,)(xxkxxb其中，其中，m 为人体的质量，为人体的质量，g 为重力加速度，为重力加速度，x 为物体的位置，第二项和第三项表示空气的阻为物体的位置，第二项和第三项表示空气的阻力。其中位置力。其中位置 x 的基准为桥梁的基准面（即选的基准为桥梁的基准面（即选择桥梁作为位置的起点择桥梁作为位置的起点 x0），低于桥梁的），低于桥梁的位置为正值，高于桥梁的位置为负值。如果人位置为正值，高于桥梁的位置为负值。如果人体系在一个弹性常数为体系在一个弹性常数为 k 的弹性绳索上，定义的弹性绳索

20、上，定义绳索下端的初始位置为绳索下端的初始位置为 0，则其对落体位置的，则其对落体位置的影响为：影响为：地面地面x桥梁基准面桥梁基准面0梯梯子子h2h1空气的阻力空气的阻力xxaxaxbmgxm 21)(整个蹦极系统的数学模型为：整个蹦极系统的数学模型为：设桥梁距离地面为设桥梁距离地面为 50 m，即，即 h2=50，蹦极者的起始位置为绳索，蹦极者的起始位置为绳索的长度的长度 30 m，即，即 h1=30，蹦极者起始速度为，蹦极者起始速度为 0，其余的参数分，其余的参数分别为别为 k20， a2a11；m70 kg，g10 m/s2。地面地面x桥梁基准面桥梁基准面0梯梯子子h2h1初始条件：初始条件：0)0(;30)0(xx已知参数：已知参数：10,70, 1,20,50,302121gmaakhh0 , 00,)(xxkxxbxxY0,00100,0102121xgYxmamaYxgYxmamamkY令令： 0300)0()0(xxY初始条件变为：初始条件变为：xxaxaxbmgxm 21)(0 , 00,)(xxkxxby0=-30; 0; % 初始位移和初始速度初始位移和初始速度t,y=ode45(bengji_sub, 0:0.01:100, y0);x1=50. - y(:,1); % x1代表蹦极