【气动学】基于龙格库塔算法实现外弹道仿真含Matlab源码

 1 简介

文着重介绍弹丸外弹道运动轨迹仿真分析,得出弹丸质点运动方程和榴弹刚体弹道方程;运用Matlab软件代入初始值用龙格-库塔进行仿真分析,仿真结果与弹道表进行对比,误差在5%以内。表明弹丸质点弹道仿真结果与弹道表基本吻合,仿真具有一定的参考意义。

2 部分代码

function comet3m(varargin)%COMET3 3-D Comet-like trajectories.% COMET3(Z) displays an animated three dimensional plot of the vector Z.% COMET3(X,Y,Z) displays an animated comet plot of the curve through the% points [X(i),Y(i),Z(i)].% COMET3(X,Y,Z,p) uses a comet of length p*length(Z). Default is p = 0.1.%% COMET3(AX,...) plots into AX instead of GCA.%% Example:% t = -pi:pi/500:pi;% comet3(sin(5*t),cos(3*t),t)%% See also COMET.% Charles R. Denham, MathWorks, 1989.% Revised 2-9-92, LS and DTP; 8-18-92, 11-30-92 CBM.% Copyright 1984-2006 The MathWorks, Inc. % $Revision: 5.11.4.4 $ $Date: 2011/03/09 07:03:36 $% Parse possible Axes input[ax,args,nargs] = axescheck(varargin{:});error(nargchk(1,10,nargs,'struct'));% Parse the rest of the inputsif nargs < 2, x = args{1}; endif nargs == 2, y = args{2}; endif nargs < 3, z = x; x = 1:length(z); y = 1:length(z); endif nargs == 3, [x,y,z] = deal(args{:}); endif nargs < 4, p = 0.10; endif nargs == 4, [x,y,z,p] = deal(args{:}); endif nargs ==10, [x,y,z,vx,vy,vz,ph,fy,hg,p] = deal(args{:}); end if ~isscalar(p) || ~isreal(p) || p < 0 || p >= 1 error(message('MATLAB:comet3:InvalidP'));endax = newplot(ax);if ~ishold(ax), [minx,maxx] = minmax(x); [miny,maxy] = minmax(y); [minz,maxz] = minmax(z); axis(ax,[minx maxx miny maxy minz maxz])endco = get(ax,'colororder');grid onif size(co,1)>=3, % Choose first three colors for head, body, and tail head = line('parent',ax,'color',co(1,:),'marker','.','MarkerSize',20,'erase','xor', ... 'xdata',x(1),'ydata',y(1),'zdata',z(1)); body = line('parent',ax,'color',co(2,:),'linestyle','-','erase','none', ... 'xdata',[],'ydata',[],'zdata',[]); tail = line('parent',ax,'color',co(3,:),'linestyle','-','erase','none', ... 'xdata',[],'ydata',[],'zdata',[]);else % Choose first three colors for head, body, and tail head = line('parent',ax,'color',co(1,:),'marker','.','MarkerSize',20,'erase','xor', ... 'xdata',x(1),'ydata',y(1),'zdata',z(1)); body = line('parent',ax,'color',co(1,:),'linestyle','--','erase','none', ... 'xdata',[],'ydata',[],'zdata',[]); tail = line('parent',ax,'color',co(1,:),'linestyle','-','erase','none', ... 'xdata',[],'ydata',[],'zdata',[]);endmm = length(z);k = round(p*mm);TIME=0;try% Grow the body% for i = 2:k+1% j = i-1:i;% set(head,'xdata',x(i),'ydata',y(i),'zdata',z(i))% set(body,'xdata',x(j),'ydata',y(j),'zdata',z(j))% drawnow% % end% Primary loop% h =waitbar(0,'Please wait...');% set(h,'name','仿真进度');tic;m = length(x);for i = 2:m j = i-1:i; set(head,'xdata',x(i),'ydata',y(i),'zdata',z(i)) set(body,'xdata',x(j),'ydata',y(j),'zdata',z(j)) %set(tail,'xdata',x(j-k),'ydata',y(j-k),'zdata',z(j-k)) drawnow t = toc; str = format_time(TIME + t); set(findobj('tag','txt_shijian'), 'String', str); %toc 停止计数，设置显示格式，在display 显示结果 set(findobj('tag','txt_feixinggaodu'), 'String', z(i)); set(findobj('tag','txt_vx'), 'String', vx(i)); set(findobj('tag','txt_vy'), 'String', vy(i)); set(findobj('tag','txt_vz'), 'String', vz(i)); set(findobj('tag','txt_hesudu'), 'String', sqrt(vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i))); set(findobj('tag','txt_pianhang'), 'String', ph(i)); set(findobj('tag','txt_fuyang'), 'String', fy(i)); set(findobj('tag','txt_henggun'), 'String', hg(i)); set(findobj('tag','txt_weizhizuobiao'), 'String', ['(' num2str(x(i)) ',' num2str(z(i)) ',' num2str(y(i)) ')']);% waitbar(i/m,h,[num2str(i*100/m) '%']);end% Clean up the tailfor i = m+1:m+k j = i-1:i; set(tail,'xdata',x(j-k),'ydata',y(j-k),'zdata',z(j-k)) drawnowendcatch E if ~strcmp(E.identifier, 'MATLAB:class:InvalidHandle') rethrow(E); endend % same subfunction as in cometfunction [minx,maxx] = minmax(x)minx = min(x(isfinite(x)));maxx = max(x(isfinite(x)));if minx == maxx minx = maxx-1; maxx = maxx+1;endfunction str = format_time(t) hrs = floor(t/3600); min = floor(t/60 - 60*hrs); sec = t - 60*(min + 60*hrs); if hrs < 10 h = sprintf('0%1.0f:', hrs); else h = sprintf('%1.0f:', hrs); end if min < 10 m = sprintf('0%1.0f:', min); else m = sprintf('%1.0f:', min); end if sec < 9.9995 s = sprintf('0%1.3f', sec); else s = sprintf('%1.3f', sec); end str = [h m s];

3 仿真结果

​

编辑

4 参考文献

[1]董理赢;王志军;焦志刚;王少宏;. 基于Matlab对弹丸外弹道运动轨迹仿真分析[C]// OSEC首届兵器工程大会论文集. 2017.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。

​

编辑