微分方程的数值解法matlab四阶龙格-库塔法课件

微分方程的数值解法四阶龙格—库塔法（TheFourth-OrderRunge－KuttaMethod）1PPT课件常微分方程（Ordinarydifferentialequations,ODE)初值问题---给出初始值边值问题---给出边界条件与初值常微分方程解算有关的指令ode23ode45

ode113ode23tode15sode23sode23tb2PPT课件一.解ODE的基本机理:2.把高阶方程转换成一阶微分方程组1.列出微分方程初始条件令(2.1)

(2.2)(2.3)3PPT课件例：著名的VanderPol方程

令降为一阶初始条件4PPT课件3.根据式(2.2)编写计算导数的M函数文件-ODE文件把t,Y作为输入宗量，把作为输出宗量

%Mfunctionfilename:dYdt.m

function

Yd=f(t,Y)

Yd=f(t,Y)

的展开式例VanderPol方程

%Mfunctionfilename:dYdt.m

function

Yd=f(t,Y)

Yd=zeros(size(Y));5PPT课件4.使编写好的ODE函数文件和初值供微分方程解算指令（solver）调用Solver解算指令的使用格式[t,Y]=solver(‘ODE函数文件名’,t0,tN,Y0,tol);ode45输出宗量形式说明：t0：初始时刻；tN：终点时刻Y0：初值；tol：计算精度6PPT课件例题1：著名的VanderPol方程

%主程序

（程序名：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命令7PPT课件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)];8PPT课件9PPT课件解法指令解题类型特点适合场合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解算指令的特点10PPT课件二.四阶Runge-Kutta

法对I=[a,b]作分割步长11PPT课件初值问题的数值解法分为两大类单步法-Runge-Kutta

方法多步法-Admas方法计算的近似值时只用到，是自开始方法Runge-Kutta法是常微分方程的一种经典解法MATLAB对应命令：ode4512PPT课件四阶Runge-Kutta公式13PPT课件四阶Runge-Kutta

法计算流程图开始Nextifori=1:N

Plot初始条件：；

积分步长：迭代次数：输出结果子程序计算End14PPT课件三.Runge-Kutta

法解VanderPol

方程的Matlab

程序结构

主程序：RK_vanderpol.m

子程序：RK_sub.m（函数文件）15PPT课件解法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);grid16PPT课件17PPT课件子程序：RK_sub.m

functionydot=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)];18PPT课件四.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：计算精度19PPT课件3月15日作业:1.VanderPol

方程的两种解法:1)采用ode45命令2)Runge-Kutta方法2.Duffing方程的求解（Runge-Kutta方法，计算步长h=0.005，计算时间t0=0.0,tN=100）要求：写出程序体，打印所绘图形，图形标题用个人的名字。Duffing

方程20PPT课件21PPT课件五.动力学系统的求解1.动力学方程2.二阶方程转成一阶方程(1)令：(2)22PPT课件其中：即：(2)23PPT课件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输出数据或图形24PPT课件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+P125PPT课件例题2：三自由度质量弹簧系统m1m2m3k1k2k3x1x2x3k4P0sin(wt)26PPT课件矩阵表示其中：27PPT课件动力学方程：解析解：已知参数：m1=m2=m3=1,k1=2,k2=2,K3=1,K4=2,P0=1,要求：采用四阶龙格－库塔法编程计算三个质量的响应时程.计算时间0~50例如：28PPT课件4阶龙格－库塔法的结果ode45的结果第一个质量的位移响应时程结果完全一致MATLAB程序（1）4阶RK方法：

（2）采用ode45：m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m29PPT课件例题3：蹦极跳系统的动态仿真蹦极者系着一根弹性绳从高处的桥梁（或山崖等）向下跳。在下落的过程中，蹦极者几乎处于失重状态。按照牛顿运动规律，自由下落的物体由下式确定：其中，m

为人体的质量，g为重力加速度，x

为物体的位置，第二项和第三项表示空气的阻力。其中位置x

的基准为桥梁的基准面（即选择桥梁作为位置的起点x＝0），低于桥梁的位置为正值，高于桥梁的位置为负值。如果人体系在一个弹性常数为k

的弹性绳索上，定义绳索下端的初始位置为0，则其对落体位置的影响为：地面x桥梁基准面0梯子h2h1空气的阻力30PPT课件整个蹦极系统的数学模型为：设桥梁距离地面为50m，即h2=50，蹦极者的起始位置为绳索的长度30m，即h1=30，蹦极者起始速度为0，其余的参数分别为k＝20，a2＝a1＝1；m＝70kg，g＝10m/s2。地面x桥梁基准面0梯子h2h1初始条件：已知参数：31PPT课件令：初始条件变为：32PPT课件y0=[-30;0];%初始位移和初始速度[t,y]=ode45(‘bengji_sub’,[0:0.01:100],y0);x1=50.-y(:,1);%x1代表蹦极者与地面之间的距离plot(t,x1);gridplot