资源预览内容
第1页 / 共36页
第2页 / 共36页
第3页 / 共36页
第4页 / 共36页
第5页 / 共36页
第6页 / 共36页
第7页 / 共36页
第8页 / 共36页
第9页 / 共36页
第10页 / 共36页
亲,该文档总共36页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述
微分方程的数值解法,四阶龙格库塔法 (The Fourth-Order RungeKutta Method),常微分方程(Ordinary differential equations, ODE),初值问题-给出初始值 边值问题-给出边界条件,与初值常微分方程解算有关的指令 ode23 ode45 ode113 ode23t ode15s ode23s ode23tb,一.解ODE的基本机理:,2. 把高阶方程转换成一阶微分方程组,1. 列出微分方程,初始条件,令,(2.1),(2.2),(2.3),例:著名的Van der Pol方程,令,降为一阶,初始条件,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);,4. 使编写好的ODE函数文件和初值 供微分方程解算指令(solver)调用,Solver解算指令的使用格式,输出宗量形式,说明: t0:初始时刻;tN:终点时刻Y0:初值; tol:计算精度,例题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:采用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);,各种solver 解算指令的特点,二. 四 阶 Runge-Kutta 法,对 I=a,b作分割,步长,单步法-Runge-Kutta 方法 多步法-Admas方法,计算 的近似值 时只用到 ,是自开始方法,Runge-Kutta法是常微分方程的一种经典解法 MATLAB 对应命令:ode45,四阶Runge-Kutta公式,四 阶 Runge-Kutta 法计算流程图,开始,Plot,初始条件: ; 积分步长: 迭代次数:,输出结果,子程序计算,End,三. Runge-Kutta 法解Van der Pol 方程的Matlab 程序结构主程序:RK_vanderpol.m 子程序:RK_sub.m(函数文件),解法2:采用Runge_Kutta法编程计算,主程序:RK_vanderpol.m t0=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 = 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; end subplot (121), plot (t, yy1, t, yy2); grid subplot (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 (函数文件名, 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) 要求:写出程序体,打印所绘图形,图形标题用个人的名字。,Duffing 方程,五. 动力学系统的求解,1. 动力学方程,2. 二阶方程转成一阶方程,(1),令:,(2),其中:,即:,(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 + 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 输出数据或图形,Matlab 程序(子程序:ZCX_sub.m),function ydot = 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,例题2:三自由度质量弹簧系统,矩阵表示,其中:,动力学方程:,解析解:,已知参数:m1=m2=m3=1, k1=2, k2=2, K3=1, K4=2, P0=1, 要求:采用四阶龙格库塔法编程计算三个质量的响应时程.计算时间 0 50,例如:,4阶龙格库塔法的结果,ode45 的结果,第一个质量的位移响应时程,结果完全一致,MATLAB程序 (1)4阶RK方法: (2)采用ode45: m_chap2_ex2_1.m,m_chap2_ex2_1_sub.m,例题3: 蹦极跳系统的动态仿真,蹦极者系着一根弹性绳从高处的桥梁(或山崖等)向下跳。在下落的过程中,蹦极者几乎处于失重状态。按照牛顿运动规律,自由下落的物体由下式确定:,其中,m 为人体的质量,g 为重力加速度,x 为物体的位置,第二项和第三项表示空气的阻力。其中位置 x 的基准为桥梁的基准面(即选择桥梁作为位置的起点 x0),低于桥梁的位置为正值,高于桥梁的位置为负值。如果人体系在一个弹性常数为 k 的弹性绳索上,定义绳索下端的初始位置为 0,则其对落体位置的影响为:,空气的阻力,整个蹦极系统的数学模型为:,设桥梁距离地面为 50 m,即 h2=50,蹦极者的起始位置为绳索的长度 30 m,即 h1=30,蹦极者起始速度为 0,其余的参数分别为 k20, a2a11;m70 kg,g10 m/s2。,初始条件:,已知参数:,初始条件变为:,y0=-30; 0; % 初始位移和初始速度 t,y=ode45(bengji_sub, 0:0.01:100, y0); x1=50. - y(:,1); % x1代表蹦极者与地面之间的距离 plot(t,x1); grid plot(t,y(:,1); grid % y(:,1)代表位移,主程序 (程序名:bengji.m),Matlab程序,function ydot=f(t,y) m=70; k=20; a1=1; a2=1; g=10; x=y(1); % x代表蹦极者的位移 x_dot=y(2); % x_dot 代表 x 的速度 if x0 ydot=0,1;-k/m,-a1/m-(a2/m)*abs(x_dot)*y+0;g; else ydot=0,1;0,-a1/m-(a2/m)*abs(x_dot)*y+0;g; end,子程序 (程序名:bengji_sub.m ),y(:,1),x1,结果分析: 右上图为蹦极者与地面之间的距离。从结果可看出,对于体重为 70 kg 的蹦极者,此系统是不安全的,因为蹦极者与地面之间的距离出现了负值 。因此,必须使用弹性系数较大的弹性绳索,才能保证蹦极者的安全。,作业(书面作业,写出程序体): (1)三自由度模型仿真 (自编Runge-Kutta 法) (2)蹦极跳模型仿真 (解算指令ode45),
网站客服QQ:2055934822
金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号