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

例题3： 蹦极跳系统的动态(dòngtài)仿真 蹦极者系着一根弹性绳从高处的桥梁（或山崖等）向下跳。在下落的过程中，蹦极者几乎处于失重状态(zhuàngtài)。按照牛顿运动规律，自由下落的物体由下式确定： 其中，m 为人体的质量，g 为重力加速度，x 为物体的位置，第二项和第三项表示空气的阻力。其中位置 x 的基准为桥梁的基准面（即选择(xuǎnzé)桥梁作为位置的起点 x＝0），低于桥梁的位置为正值，高于桥梁的位置为负值。如果人体系在一个弹性常数为 k 的弹性绳索上，定义绳索下端的初始位置为 0，则其对落体位置的影响为： 地面 x 桥梁基准面 0 梯子 h2 h1 空气的阻力 第三十页，共36页。 整个蹦极(bènɡ jí)系统的数学模型为： 设桥梁距离地面为 50 m，即 h2=50，蹦极者的起始位置为绳索(shénɡ suǒ)的长度 30 m，即 h1=30，蹦极者起始速度为 0，其余的参数分别为 k＝20， a2＝a1＝1；m＝70 kg，g＝10 m/s2。 地面 x 桥梁基准面 0 梯子 h2 h1 初始条件： 已知参数(cānshù)： 第三十一页，共36页。 令： 初始条件变为： 第三十二页，共36页。 y0=[-30; 0]; % 初始位移和初始速度 [t,y]=ode45(‘bengji_sub’, [0:0.01:100], y0); x1=50. - y(:,1); % x1代表(dàibiǎo)蹦极者与地面之间的距离 plot(t,x1); grid plot(t,y(:,1)); grid % y(:,1)代表(dàibiǎo)位移 主程序(chéngxù) (程序(chéngxù)名：bengji.m) Matlab程序(chéngxù) 第三十三页，共36页。 function ydot=f(t,y) m=70; k=20; a1=1; a2=1; g=10; x=y(1); % x代表蹦极(bènɡ jí)者的位移 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 子程序(chéngxù) （程序(chéngxù)名：bengji_sub.m ） 第三十四页，共36页。 微分方程的数值(shùzí)解法 四阶龙格—库塔法 （The Fourth-Order Runge－Kutta Method） 第一页，共36页。 常微分方程(wēi fēn fānɡ chénɡ)（Ordinary differential equations, ODE) 初值问题---给出初始值 边值问题---给出边界条件 与初值常微分方程解算有关(yǒuguān)的指令 ode23 ode45 ode113 ode23t ode15s ode23s ode23tb 第二页，共36页。 一.解ODE的基本(jīběn)机理: 2. 把高阶方程(fāngchéng)转换成一阶微分方程(fāngchéng)组 1. 列出微分方程(wēi fēn fānɡ chénɡ) 初始条件 令 (2.1) (2.2) (2.3) 第三页，共36页。 例：著名(zhùmíng)的Van der Pol方程 令 降为一阶 初始条件 第四页，共36页。 3. 根据式(2.2)编写(biānxiě)计算导数的M函数文件-ODE文件 把t,Y作为输入(shūrù)宗量，把 作为输出宗量 %M function file name: dYdt.m  function Yd = f (t, Y) Yd = f (t,Y) 的展开式 例Van der Pol方程(fāngchéng) %M function file name: dYdt.m  function Yd = f (t, Y) Yd=zeros(size(Y)); 第五页，共36页。 4. 使编写好的ODE函数文件和初值 供微分方程(wēi fēn fānɡ chénɡ)解算指令（solver）调用 Solver解算指令(zhǐlìng)的使用格式 [t