matlab第二章常微分方程的数值解法课件.ppt

《matlab第二章常微分方程的数值解法课件.ppt》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、微分方程的数值解法四阶龙格—库塔法（TheFourth-OrderRunge－KuttaMethod）常微分方程（Ordinarydifferentialequations,ODE)初值问题---给出初始值边值问题---给出边界条件与初值常微分方程解算有关的指令ode23ode45ode113ode23tode15sode23sode23tb一.解ODE的基本机理:2.把高阶方程转换成一阶微分方程组1.列出微分方程初始条件令(2.1)(2.2)(2.3)例：著名的VanderPol方程令降为一阶初始条件3.根据式(2.2)编写计算导数的M函数文件-ODE文件把t,Y作为

2、输入宗量，把作为输出宗量%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解算指令的使用格式[t,Y]=solver(‘ODE函数文件名’,t0,tN,Y0,tol);ode45输出宗量形式说明：t0：初始时刻；tN：终点时刻Y0：初值；tol：计算精度例题1：著名的VanderPol方

3、程%主程序（程序名：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));Yd

4、ot(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算式，低精度

5、低精度时，比ode15s有效；或存在质量矩阵时ode23tb刚性采用梯形法则－反向数值微分两阶段算法，低精度低精度时，比ode15s有效；或存在质量矩阵时各种solver解算指令的特点二.四阶Runge-Kutta法对I=[a,b]作分割步长初值问题的数值解法分为两大类单步法-Runge-Kutta方法多步法-Admas方法计算的近似值时只用到，是自开始方法Runge-Kutta法是常微分方程的一种经典解法MATLAB对应命令：ode45四阶Runge-Kutta公式四阶Runge-Kutta法计算流程图开始Nextifori=1:NPlot初始条件：；积分步长：迭代次

6、数：输出结果子程序计算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(

7、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);