常微分方程的数值解法.pptx
《常微分方程的数值解法.pptx》由会员分享,可在线阅读,更多相关《常微分方程的数值解法.pptx(35页珍藏版)》请在得力文库 - 分享文档赚钱的网站上搜索。
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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 解法
限制150内