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

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

本程序特别适合于弹道计算等方面使用。比ODE45函数的速度加快了很多，且程序充分展示了Matlab向量运算的强大功能，以及编程的简单快捷。

但愿本程序，可以为大家提高工作效率提高帮助。

一、rightF.m文件（调用的时候，要根据自己的应用，改变右函数值）

%右函数定义的，用于单步龙格库塔法计算。

%如果需要使用全局变量，要使用global定义。

%这种格式的，也完全可以被ode45函数调用。

function dy=rightF(t,y)%right functions

dy=zeros(2,1);%a column vector

dy(1)=y(2);

dy(2)=3*y(2)/t-4*y(1)/(t^2);

二、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=[10];

h=0.1;

for i=1:10

b=stepRK(t,a,h)

t=t+h;

a=b;

arrayT(i)=i;

array1(i)=b(1);

end

plot(arrayT,array1,'*')

四、比例导引弹道计算程序

clc

clear

disp'***Simulation Start**please wait'

global 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)&(RTM>3)

RTM1=RT1-RM1;

RTM2=RT2-RM2;

RTM=sqrt(RTM1^2+RTM2^2);%弹目距离

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 XNC>Max_G

XNC=Max_G;

end

if XNC<-Max_G

XNC=-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;

end

figure

axis equal

plot(ArrayRM1,ArrayRM2),grid

%title('Engagement Geometry') xlabel('Downrange(m)')

ylabel('Altitude(m)')

figure

plot(ArrayT,ArrayXNCG),grid

%title('Commanded Acceleration') xlabel('t(s)')

ylabel('Commanded Acceleration(g)')