防灾数学分享 http://blog.sciencenet.cn/u/fzmath 防灾科技学院数学教研室

博文

用MATLAB模拟简单小球跳动

已有 5346 次阅读 2017-9-30 16:06 |系统分类:教学心得

Bouncing Ball 的图像结果      

这里使用了ode求解器的event事件技巧


小球满足的微分方程

\begin{array}{cl} y'_1 &= y_2\\ y'_2 &= -9.8 . \end{array}

写作

function dydt = f(t,y)
dydt = [y(2); -9.8];

小球在接触地面时会反弹,事件函数可写为

function [value,isterminal,direction] = bounceEvents(t,y)
value = y(1);     % Detect height = 0
isterminal = 1;   % Stop the integration
direction = -1;   % Negative direction only

运行主函数为

tstart = 0;
tfinal = 30;
y0 = [0; 20];
refine = 4;
options = odeset('Events',@events,'OutputFcn',@odeplot,'OutputSel',1,...
  'Refine',refine);
fig = figure;
ax = axes;
ax.XLim = [0 30];
ax.YLim = [0 25];
box on
hold on;
tout = tstart;
yout = y0.';
teout = [];
yeout = [];
ieout = [];
for i = 1:10
  % Solve until the first terminal event.
  [t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options);
  if ~ishold
     hold on
  end
  % Accumulate output.  This could be passed out as output arguments.
  nt = length(t);
  tout = [tout; t(2:nt)];
  yout = [yout; y(2:nt,:)];
  teout = [teout; te];          % Events at tstart are never reported.
  yeout = [yeout; ye];
  ieout = [ieout; ie];
 
  ud = fig.UserData;
  if ud.stop
     break;
  end
 
  % Set the new initial conditions, with .9 attenuation.
  y0(1) = 0;
  y0(2) = -.9*y(nt,2);
 
  % A good guess of a valid first timestep is the length of the last valid
  % timestep, so use it for faster computation.  'refine' is 4 by default.
  options = odeset(options,'InitialStep',t(nt)-t(nt-refine),...
     'MaxStep',t(nt)-t(1));
 
  tstart = t(nt);
end
plot(teout,yeout(:,1),'ro')
xlabel('time');
ylabel('height');
title('Ball trajectory and the events');
hold off
odeplot([],[],'done');







https://blog.sciencenet.cn/blog-292361-1078519.html

上一篇:用MATLAB绘制空间曲面曲面及其切平面
下一篇:用MATLAB编写预估-校正法程序求分数阶常微分方程组数值解
收藏 IP: 124.238.135.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-20 05:34

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部