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

微分方程的数值解法四阶龙格—库塔法（TheFourth-OrderRunge－KuttaMethod）常微分方程（Ordinarydifferentialequations,ODE)初值问题---给出初始值边值问题---给出边界条件与初值常微分方程解算有关的指令ode23ode45ode113ode23tode15sode23sode23tb一.解ODE的基本机理:2.把高阶方程转换成一阶微分方程组1.列出微分方程初始条件令)0()0(210yyYyyyy21,),(),(),(2121YtfYtfYtfyyY(2.1)(2.2)(2.3)),,(tyyfy00)0(,)0(yyyy例：著名的VanderPol方程0)1(2yyyy令yyyy21,降为一阶1221221)1(yyyyyyY21yyY初始条件2010210)0()0(yyyyY3.根据式(2.2)编写计算导数的M函数文件-ODE文件把t,Y作为输入宗量，把作为输出宗量Y%Mfunctionfilename:dYdt.mfunctionYd=f(t,Y)Yd=f(t,Y)的展开式例VanderPol方程%Mfunctionfilename:dYdt.mfunctionYd=f(t,Y)Yd=zeros(size(Y)););1()2(*)12).^1(()2();2()1(YYYYdYYd4.使编写好的ODE函数文件和初值供微分方程解算指令（solver）调用0YSolver解算指令的使用格式[t,Y]=solver(‘ODE函数文件名’,t0,tN,Y0,tol);ode45)()()()()()(222112110201tytytytytytyY输出宗量形式)()1(:,1tyY)()2(:,2tyY说明：t0：初始时刻；tN：终点时刻Y0：初值；tol：计算精度例题1：著名的VanderPol方程0)1(2yyyy%主程序（程序名：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)];解法指令解题类型特点适合场合ode45非刚性采用4、5阶Runge－Kutta法大多数场合的首选算法ode23非刚性采用Adams算法较低精度（10－3）场合ode113非刚性多步法；采用Adams算法；高低精度均可（10－3～10－6）ode45计算时间太长时取代ode45ode23t适度刚性采用梯形法则算法适度刚性ode15s刚性多步法；采用2阶Rosenbrock算式，精度中等当ode45失败时使用；或存在质量矩阵时ode23s刚性一步法；采用2阶Rosenbrock算式，低精度低精度时，比ode15s有效；或存在质量矩阵时ode23tb刚性采用梯形法则－反向数值微分两阶段算法，低精度低精度时，比ode15s有效；或存在质量矩阵时各种solver解算指令的特点二.四阶Runge-Kutta法00)(),(ytybtaytfdtdybtttaN110对I=[a,b]作分割1,,1,0,1Nihhhiii步长初值问题的数值解法分为两大类单步法-Runge-Kutta方法多步法-Admas方法计算的近似值时只用到，是自开始方法)(1ntynnyt,1nyRunge-Kutta法是常微分方程的一种经典解法MATLAB对应命令：ode45四阶Runge-Kutta公式),()2,21()2,21(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn四阶Runge-Kutta法计算流程图开始Nextifori=1:NPlot)22(643211KKKKhyynn11,nnnnyytt初始条件：；积分步长：迭代次数：0t0yhN输出结果httnn10ttn子程序计算iKEnd三.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)默认精度：ode23——1e-3ode45——1e-6说明：t0：初始时刻；tN：终点时刻y0：初值；tol：计算精度3月15日作业:1.VanderPol方程的两种解法:1)采用ode45命令2)Runge-Kutta方法2.Duffing方程的求解（Runge-Kutta方法，计算步长h=0.005，计算时间t0=0.0,tN=100）要求：写出程序体，打印所绘图形，图形标题用个人的名字。)2.1cos(28.03.03tyyyyDuffing方程0.0)0(,01.0)0(yy五.动力学系统的求解1.动力学方程2.二阶方程转成一阶方程00)0(,)0()(),()()()(XXXXXPKXXCXMnRttttt(1))0()0()0(,)(,)()()(2XXYYXXYnRtttt令：)()()(tttPAYY(2)12112211)()(,nnnnnnnnRttRPM0PCMKMI0A其中：)()()()()(1111tttttnnnnnPM0XXCMKMI0XX即：00)0(,)0()(),()()()(XXXXXPKXXCXMnRttttt)()()(tttPAYY(2)3.Matlab程序(主程序：ZCX)t0;Y0;h;N;P0,w;%输入初始值、步长、迭代次数、初始激励力；fori=1:Nt1=t0+hP=[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=%输入(t0+h)时刻的外部激励力K4=ZCX_sub(t0+h,y0+hK3,P)Y1=y0+(h/6)(K1+2K2+2K3+K4)t1,Y1(输出t1,y1)nexti输出数据或图形)()()(tttPAYYMatlab程序(子程序：ZCX_sub.m)functionydot=f(t,Y,P)M=,K=,C=%输入结构参数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()(02121111tPxxkxkxm0)()(32321222xxkxxkxm0)(3432333xkxxkxm00)sin(0000000003214333322221321321tPxxxkkkkkkkkkkxxxmmm矩阵表示其中：321xxxX321000000mmmM433332222100kkkkkkkkkkK00)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305101520253035404550-0.4-0.3-0.2-0.100.10.20.30.405101520253035404550-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：蹦极跳系统的动态仿真蹦极者系着一根弹性绳从高处的桥梁（或山崖等）向下跳。在下落的过程中，蹦极者几乎处于失重状态。按照牛顿运动规律，自由下落的物体由下式确定：xxaxamgxm210,00,)(xxkxxb其中，m为人体的质量，g为重力加速度，x为物体的位置，第二项和第三项表示空气的阻力。其中位置x的基准为桥梁的基准面（即选择桥梁作为位置的起点x＝0），低于桥梁的位置为正值，高于桥梁的位置为负值。如果人体系在一个弹性常数为k的弹性绳索上，定义绳索下端的初始位置为0，则其对落体位置的影响为：地面x桥梁基准面0梯子h2h1空气的阻力xxaxaxbmgxm21)(整个蹦极系统的数学模型为：设桥梁距离地面为50m，即h2=50，蹦极者的起始位置为绳索的长度30m，即h1=30，蹦极者起始速度为0，其余的参数分别为k＝20，a2＝a1＝1；m＝70kg，g＝10m/s2。地面x桥梁基准面0梯子h2h1初始