Advertisement
Y-BH

half_cycle

Mar 26th, 2022
1,129
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. function half_cycle
  2. % 质点离开竖直光滑半圆轨道的轨迹
  3. % 模型参数
  4. R=1;L=4*R;g=10; % 轨道半径、水平轨长、重力加速度
  5. vu=sqrt(g*R);tu=sqrt(R/g); % 单位速度、单位时间
  6. T=100*tu;%模拟时间上限
  7. %半圆轨道
  8. tht=-0.5*pi:0.01:0.5*pi;
  9. xC=R*cos(tht);yC=R+R*sin(tht);
  10. %圆轨道
  11. tht=-pi:0.01:pi;
  12. xO=R*cos(tht);yO=R+R*sin(tht);
  13. % ode设置
  14. opts=odeset('RelTol',1e-6,'AbsTol',1e-6, ...
  15.     'MaxStep',tu*1e-2,'InitialStep',tu*1e-2, ...
  16.     'Events',@(t,u) luodi(t,u,R));
  17. v1=sqrt(2)*vu;%恰好到达圆轨道半高的速度
  18. v2=sqrt(5)*vu;%恰好经过圆轨道顶点的速度
  19. v3=sqrt(7/2)*vu;%恰好落在圆轨道底部的速度
  20. % 恰好通过最高点轨迹
  21. [~,yd] = ode45(@(t,u) half_c(t,u,g,R), ...
  22.     [0 T],[0 0 -0.5*pi v2 v2 0]',opts);
  23. % 恰好落至圆轨道底部轨迹
  24. [~,yb] = ode45(@(t,u) half_c(t,u,g,R), ...
  25.     [0 T],[0 0 -0.5*pi v3 v3 0]',opts);
  26. NN=16;%小球个数
  27. % 微分方程数值计算
  28. for ii=1:NN
  29.     v0=sqrt(ii/2)*vu;%初速度
  30.     [tN{ii},uN{ii}]=ode45(@(t,u)half_c(t,u,g,R), ...
  31.         0:0.005:T,[0 0 -0.5*pi v0 v0 0]',opts);
  32.     LT(ii)=length(tN{ii});
  33. end
  34. % 数组填充对齐
  35. for ii=1:NN
  36.     dL=max(LT)-LT(ii);
  37.     if dL>0
  38.         u=uN{ii};
  39.         uN{ii}=[u;repmat(u(end,:),dL,1)];
  40.     end
  41.     u=uN{ii};
  42.     X(:,ii)=u(:,1);Y(:,ii)=u(:,2);
  43. end
  44. %% 模型绘图
  45. pos1 = [0 0 1 1];pos2 = [0.05 0.05 0.9 0.9];
  46. subplot('Position',pos2)
  47. dpi=1.25;scr=[1920 1080];res=[1280 720];
  48. fig=gcf;fig.Color='w';
  49. fig.Position = [(scr-res)/2,res]/dpi;% 居中显示
  50. %半圆轨道
  51. plot(xC,yC,'Color','k','LineWidth',4);hold on
  52. % 圆轨道
  53. plot(xO,yO,'-.','Color','k','LineWidth',0.5)
  54. %水平轨道
  55. line([-L 0],[0 0],'linestyle','-','Color','k','LineWidth',4)
  56. % 恰好通过最高点轨迹
  57. pd=plot(yd(:,1),yd(:,2),'--','Color','#444444','LineWidth',1.5);
  58. % 恰好落至圆轨道底部轨迹
  59. pb=plot(yb(:,1),yb(:,2),'-','Color','#444444','LineWidth',1);
  60. % 动画线条
  61. for k=NN:-1:1
  62.     p(k)=animatedline('Color','#D95319','LineWidth',1.5);
  63. end
  64. % 小球
  65. pp=plot(0,0,'.','markersize',30,'Color','#D95319');
  66. % 图例
  67. legend([pd pb],{'恰好通过最高点轨迹','恰好落至圆轨道底部轨迹'},...
  68.     'FontSize',16,'FontName','黑体','Location','northwest');
  69. legend('boxoff')
  70. axis([-4 1 0 2]);axis equal;axis tight% axis off
  71. %% 动画部分
  72. mv = 0;% 是否录制视频,1录制,其它不录制
  73. if mv == 1
  74. v = VideoWriter('D:\Programs\MATLAB\小球脱轨.mp4','MPEG-4');
  75. v.Quality = 10;v.FrameRate = 30;
  76. open(v);
  77. end
  78. for ii=1:max(LT)
  79.     for jj=NN:-1:1
  80.         u=uN{jj};
  81.         addpoints(p(jj),u(ii,1),u(ii,2));
  82.     end
  83.     pp.XData=X(ii,:);pp.YData=Y(ii,:);
  84.     drawnow% limitrate
  85.     if mv == 1;F = getframe(fig);writeVideo(v,F);end
  86. end
  87. if mv == 1;close(v);end
  88. %% 微分方程
  89. function du = half_c(~,u,g,R)
  90. x=u(1);y=u(2);t=u(3);v=u(4);vx=u(5);vy=u(6);
  91. a=v^2/R;
  92. if y==0 && x<0
  93.     dx=v;dy=0;dtheta=0;
  94.     dv=0;dvx=0;dvy=0;
  95. elseif x>=0 && g*sin(t)<a
  96.     dx=-v*sin(t);dy=v*cos(t);
  97.     dtheta=v/R;dv=-g*cos(t);
  98.     dvx=-a*cos(t)-dv*sin(t);
  99.     dvy=-a*sin(t)+dv*cos(t);
  100. else
  101.     dx=vx;dy=vy;dtheta=0;
  102.     dv=0;dvx=0;dvy=-g;
  103. end
  104. du=[dx dy dtheta dv dvx dvy]';% 应返回列向量
  105. %% 事件函数
  106. function [position,isterminal,direction] = luodi(~,u,R)
  107. if u(1)>0 && u(2)<R && u(3)>0
  108.     position = u(2)+sqrt(R^2-u(1)^2)-R;% 落在圆轨道
  109. else
  110.     position = u(2);% 落回水平面
  111. end
  112. isterminal = 1;
  113. direction = -1;
Advertisement
RAW Paste Data Copied
Advertisement