《单步龙格库塔比例导引弹道计算matlab源程序.doc》由会员分享，可在线阅读，更多相关《单步龙格库塔比例导引弹道计算matlab源程序.doc（4页珍藏版）》请在金锄头文库上搜索。
1、1单步龙格库塔比例导引弹道计算 matlab 源程序本程序特别适合于弹道计算等方面使用。比 ODE45 函数的速度加快了很多，且程序充分展示了 Matlab 向量运算的强大功能，以及编程的简单快捷。但愿本程序，可以为大家提高工作效率提高帮助。一、rightF.m 文件 （调用的时候，要根据自己的应用，改变右函数值）%右函数定义的，用于单步龙格库塔法计算。%如果需要使用全局变量，要使用 global 定义。% 这种格式的，也完全可以被 ode45 函数调用。function dy = rightF(t,y) % right functionsdy = zeros(2,1); % a column
2、 vectordy(1) = y(2);dy(2) = 3*y(2)/t-4*y(1)/(t2);二、stepRK.m 文件 （本函数任何情况都不用变的）%函数说明，t 表示自变量（一般为时间） ，y 表示变量的初值（是一个行向量） ，h 表示步长%运行的结果，也是一个行向量function dy = stepRK(t,y,h)K1=rightF(t,y);K2=rightF(t+h/2,y+h*K1/2);K3=rightF(t+h/2,y+h*K2/2);K4=rightF(t+h,y+h*K3);dy=y+(K1+2*K2+2*K3+K4)*h/6;dy=dy;三、示例程序（rk_tes
3、t.m）%龙格库塔法测试程序，与欧阳联渊书一致t=1.;a=1 0;h=0.1;for i=1:10b=stepRK(t,a,h)t=t+h;a=b;arrayT(i)=i;array1(i)=b(1);endplot(arrayT,array1,*)2四、比例导引弹道计算程序clccleardisp * Simulation Start * please waitglobal XNC VM;RM1=0.; %导弹初始位置RM2=1500.;RT1=3000.;%目标（或者假想目标点）初始位置RT2=0.;RTM1=RT1-RM1;RTM2=RT2-RM2;VM=260.;Max_G=9.8*
4、15;XLAMF=-90/57.3; %命中目标时候导弹的攻击角（impact angle）要求theta=0; %导弹初始速度方向(角度)VM1=VM*cos(theta);VM2=VM*sin(theta);y(1)=theta; y(2)=VM1;y(3)=VM2;y(4)=RM1;y(5)=RM2; T=0.;T_step=.01; %仿真步长n=1; %这个序号用于记录数据用VC=100; RTM=100; % 只是为了不为 0,好编写循环运算的代码%-while (VC = 0)&(RTM3) RTM1=RT1-RM1; RTM2=RT2-RM2;RTM=sqrt(RTM12+RT
5、M22); %弹目距离XLAM=atan2(RTM2,RTM1);%根据弹目线的连线，利用几何关系来计算VTM1=-VM1; %相对速度， Zarchan book P143VTM2=-VM2;VC=-(RTM1*VTM1+RTM2*VTM2)/RTM;XLAM_dot=(RTM1*VTM2-RTM2*VTM1)/(RTM*RTM);XNC=4.*VC*XLAM_dot;%-过载限制-if XNCMax_GXNC=Max_G;endif XNC-Max_GXNC=-Max_G;end%-过载限制代码结束 - %-垂直于 LOS 使用这个代码，同时右函数计算的，也需要更改-% AM1=-XNC*
6、sin(XLAM);% AM2=XNC*cos(XLAM);%-%调用子程序,进行一步龙格库塔积分运算,解出弹道参数y=stepRK(T,y,T_step);theta=y(1);VM1=y(2);VM2=y(3);RM1=y(4);RM2=y(5); ArrayT(n)=T;ArrayRM1(n)=RM1;ArrayRM2(n)=RM2;ArrayRT1(n)=RT1;ArrayRT2(n)=RT2;Array_theta_degree(n)=theta*57.3;XNCG=XNC/9.8; gravity_compensate=cos(theta);ArrayXNCG(n)=XNCG+gr
7、avity_compensate;%“重力补偿”因素，也考虑进导弹的过载里面了。XLAM_degree=XLAM*57.3;ArrayXLAM_degree(n)=XLAM_degree; %这个值是负的。ArraySeeker_degree(n)=(theta-XLAM)*57.3+XNCG*3; % the angle between LOS and missile body(也就是导引头的框架角)，把攻角因素也考虑了。 4T=T+T_step;n=n+1;endfigureaxis equal plot(ArrayRM1,ArrayRM2),grid%title(Engagement Geometry)xlabel(Downrange (m) )ylabel(Altitude (m)figureplot(ArrayT,ArrayXNCG),grid%title(Commanded Acceleration)xlabel(t(s) )ylabel(Commanded Acceleration(g)
