Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % 质点离开竖直光滑半圆轨道的轨迹
- % 如运行有问题,建议使用R2020a或较新版本
- % R2016a或较旧版本请将微分方程函数、事件函数单独存为脚本函数
- % 如需转载此脚本请注明出处:
- % @YBH、ssrs@bilibili - https://space.bilibili.com/9943102
- clc;clear
- close all
- path=char('.\pic_seqs\');% 生成的图片序列所在文件夹
- if exist(path,'dir')==7
- rmdir(path,'s')% 删除旧文件
- end
- mkdir(path)
- dpi=1.25;% 系统缩放125%
- scr=[1920 1080];% 屏幕分辨率
- res=[1280 720];% 图像分辨率
- R=1;%轨道半径
- L=4*R;%水平轨长
- g=10;%重力加速度
- vu=sqrt(g*R);%单位速度
- tu=sqrt(R/g);%单位时间
- T=100*tu;%模拟时间上限
- %半圆轨道
- tht=-0.5*pi:0.01:0.5*pi;
- xc=R*cos(tht);yc=R+R*sin(tht);
- %圆轨道
- tht=-pi:0.01:pi;
- xo=R*cos(tht);yo=R+R*sin(tht);
- % 此脚本目的为生成演示动画,
- % 常见视频分辨率为1920x1080或1280x720等等
- % 一个像素点以下的误差几乎不能察觉,故精度在1e-3即可
- % 考虑到误差积累,保守起见至少设为1e-4 /
- % 实测至少x_0=0 --> 1e-5;x_0=-L --> 1e-6
- odest=odeset('RelTol',1e-7,'AbsTol',1e-7,'NormControl','on',...
- 'MaxStep',tu*1e-2,'InitialStep',tu*1e-2,...
- 'Stats','off','Events',@(t,y) luodi_evns(t,y,R));
- v1=sqrt(2);%恰好到达圆轨道半高的速度
- v2=sqrt(5);%恰好经过圆轨道顶点的速度
- v3=sqrt(7/2);%恰好落在圆轨道底部的速度
- % 恰好通过最高点轨迹
- v=v2*vu;
- [~,yd] = ode45(@(t,y) half_c(t,y,g,R),[0 T],[0 0 -0.5*pi v v 0]',odest);
- % 恰好落至圆轨道底部轨迹
- v=v3*vu;
- [~,yb] = ode45(@(t,y) half_c(t,y,g,R),...
- [0 T],[0 0 -0.5*pi v v 0]',odest);
- sizem=0;tmax=0;
- N1=1;N=16;%小球个数
- % 微分方程数值计算
- for ii=N1:N
- v0=sqrt(ii/2)*vu;%初速度
- eval(['[t,y',num2str(ii),']',...
- ' = ode45(@(t,y) half_c(t,y,g,R),',...
- '[0 T],[0 0 -0.5*pi v0 v0 0]'',odest);']);
- eval(['t',num2str(ii), ' = t;']);
- tmax=max(tmax,max(t));%终止模拟时间
- sizem=max(sizem,size(t,1));%最大数组长度
- end
- % 最大数组不一定对应最大时间!!!
- % 数组填充对齐/与最大的数组对齐
- tmaxp=tmax*1.1;sizemp=fix(1.1*sizem);
- for ii=N1:N
- % 计算每个size(y_ii,1)=sizei
- eval(['sizei=size(y',num2str(ii),',1);']);
- % 计算每个max(t_ii)=tim
- eval(['tim=max(t',num2str(ii),');']);
- % t对齐t-->tf / 线性填充
- eval(['tf',num2str(ii),'=[t',num2str(ii),...
- '(1:sizei-1);linspace(tim,tmaxp,sizemp-sizei+1)''];']);
- % y对齐y-->yf / 拖尾填充
- eval(['yf',num2str(ii),'=[y',num2str(ii),...
- ';zeros(sizemp-sizei,6)+y',num2str(ii),'(sizei,:)];']);
- end
- % 按基准时间线插值
- fps=120;%帧率
- frames=fix(fps*tmaxp);%总帧数
- ti=linspace(0,tmaxp,frames)';% 等间隔的基准时间线
- % yi=interp1(tf,yf,ti)
- for ii=N1:N
- eval(['yi',num2str(ii),'=interp1(tf',num2str(ii),...
- ',yf',num2str(ii),',ti,''linear'');']);
- end
- % 模型绘图
- for n=1:frames
- %半圆轨道
- plot(xc,yc,'Color','k','LineWidth',4)
- hold on
- axis equal
- axis tight
- % axis off
- % 圆轨道
- plot(xo,yo,'-.','Color','k','LineWidth',0.5)
- %水平轨道
- line([-L 0],[0 0],'linestyle','-','Color','k','LineWidth',4)
- % 小球位置及运动轨迹
- for kk=N:-1:N1
- eval(['plot(yi',num2str(kk),'(n,1),yi',num2str(kk),...
- '(n,2),''.'',''markersize'',40)'])% 小球位置
- eval(['plot(yi',num2str(kk),'(1:n,1),yi',num2str(kk),...
- '(1:n,2),''-'',''LineWidth'',2)'])% 小球轨迹
- end
- % 恰好通过最高点轨迹
- pd=plot(yd(:,1),yd(:,2),...
- '--','Color','#444444','LineWidth',1.5);
- % 恰好落至圆轨道底部轨迹
- pb=plot(yb(:,1),yb(:,2),...
- '-','Color','#444444','LineWidth',1);
- legend([pd,pb],{'恰好通过最高点轨迹','恰好落至圆轨道底部轨迹'},...
- 'FontSize',16,'FontName','黑体','Location','northwest');
- legend('boxoff')
- hold off
- set(gcf,'Color','w');% 设置白色背景
- fig=gcf;
- fig.Units = 'pixel';
- fig.Position = [(scr-res)/2,res]/dpi;% 居中显示
- m(n)=getframe(fig);% 获取当前图像
- imwrite(m(n).cdata,[path,num2str(n),'.png']);% 保存图像为png
- end
- % 微分方程
- function dydt = half_c(~,y,g,R)
- a=y(4)^2/R;% 向心加速度
- if y(2)==0 && y(1)<0
- dx=y(4);
- dy=0;
- dtheta=0;
- dv=0;
- dvx=0;
- dvy=0;
- else
- if y(1)>=0 && g*sin(y(3))<a
- dx=-y(4)*sin(y(3));
- dy=y(4)*cos(y(3));
- dtheta=y(4)/R;
- dv=-g*cos(y(3));
- dvx=-a*cos(y(3))-dv*sin(y(3));
- dvy=-a*sin(y(3))+dv*cos(y(3));
- else
- dx=y(5);
- dy=y(6);
- dtheta=0;%记录脱离时角度,同dv=0
- dv=0;%不再计算v,v保持脱离时的速度,同dtheta=0
- dvx=0;
- dvy=-g;
- end
- end
- dydt=[dx dy dtheta dv dvx dvy]';% 需为列向量
- end
- % 事件函数
- function [position,isterminal,direction] = luodi_evns(~,y,R)
- if y(1)>0 && y(2)<R && y(3)>0
- position = y(2)+sqrt(R^2-y(1)^2)-R;% 落在圆轨道
- else
- position = y(2);% 落回水平面
- end
- isterminal = 1;
- direction = -1;
- end
Add Comment
Please, Sign In to add comment