单步龙格库塔比例导引弹道计算matlab源程序.doc

单步龙格库塔比例导引弹道计算matlab源程序本程序特别适合于弹道计算等方面使用。比ODE45函数的速度加快了很多，且程序充分展示了Matlab向量运算的强大功能，以及编程的简单快捷。但愿本程序，可以为大家提高工作效率提高帮助。一、rightF.m 文件 （调用的时候，要根据自己的应用，改变右函数值）%右函数定义的，用于单步龙格库塔法计算。%如果需要使用全局变量，要使用global定义。% 这种格式的，也完全可以被ode45函数调用。function dy = rightF(t,y) % right functionsdy = zeros(2,1); % a column 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_test.m）%龙格库塔法测试程序，与欧阳联渊书一致 t=1.; a=1 0; h=0.1;for i=1:10 b=stepRK(t,a,h) t=t+h; a=b; arrayT(i)=i; array1(i)=b(1);endplot(arrayT,array1,*)四、比例导引弹道计算程序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*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+RTM22); %弹目距离 XLAM=atan2(RTM2,RTM1);%根据弹目线的连线，利用几何关系来计算 VTM1=-VM1; %相对速度，Zarchan book P14 VTM2=-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*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+gravity_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(也就是导引头的框架角)，把攻角因素也考虑了。 T=T+T_step; n=n+1;endfigureaxis equal plot(ArrayRM1,ArrayRM2),grid%title(Engagement Geometry)xlabel(Downrange (m) )ylabel(Altitude