导弹运动方程求解数值解离不开MATLAB，本文提出一些基本的方法与技巧。
1.微分方程的求解：ode45。给出例子：
function dy = lateral( t,y )
%求解偶然干扰作用下侧向扰动运动过渡过程
%y1=wx,y2=wy,y3=bt,y4=gm
b11=1.86;b12=0.66;b14=8.8;b21=0.02;b22=0.2;bp=0;%bp=b24'
b24=1.34;b36=-1;b34=0.06;a33=0;b35=-0.028;b17=-0.78;
b56=0.0012;b18=0.98;b37=0.018;b27=-0.9;
af=15/57.3;
My=0;
Jy=5.3;
dy=zeros(4,1);
dy(1)=-b11*y(1)-b12*y(2)-b14*y(3);
dy(2)=-(b21+bp*af)*y(1)-(b22-bp*b36+bp*af*b56)*y(2)-(b24-bp*b34-bp*a33)*y(3)+bp*b35*y(4)+My/Jy;
dy(3)=af*y(1)-(b36-af*b56)*y(2)-(b34+a33)*y(3)-b35*y(4);
dy(4)=y(1)+b56*y(2);
end
函数结束，命令行代码
[t,y]=ode45(@lateral,[0 10],[0 02/57.3 0]);%0~10s的过渡过程，有初值2°
plot(t,y(:,1),'r-',t,y(:,2),'b-',t,y(:,3),'y-',t,y(:,4),'m-');
title('自由扰动');
legend('Wx','Wy','贝塔','伽马')
plot基本画图函数，legend是给每条线标注的
2.求取特征值与特征向量
状态矩阵A，定义符号变量s，使用det函数，例子：
A=[-b11 -b12 -b14 0;
-(b21+bp*af) -(b22-bp*b36+bp*af*b56) -(b24-bp*b34-bp*a33) bp*b35;
af -(b36-af*b56) -(b34+a33) -b35;
1 b56 0 0];%状态矩阵
symss;%定义符号变量s
B=s*eye(4);%eye(4):四阶单位矩阵
G=det(B-A);%矩阵求值，即特征方程式，使用MATLAB内置函数
x=eig(A);%求状态矩阵A的特征值，使用MATLAB内置函数
3.画稳定边界图
求解出直线方程后，一般式直线方程用ezplot画线，而ezplot不能设置颜色，默认薄荷绿色，可以使用set函数设置颜色，例子：
h=ezplot(A2,[-3,6,-3,2]);%画A2=0的图像，x轴范围[-3,6],y轴范围[-3,2]
set(h,'color','r')%设置图像颜色为红色
可用text函数标注稳定边界对应的线条，例子
text(1.5,-1,'A2=0\rightarrow')%在坐标（1.5，-1）处画右箭头，写A2=0
4.使用终止条件
想要知道导弹恰好飞多高或速度恰好多少时的各解，难道我们要一次次试吗？使用ode事件属性events即可。
假设我们要让导弹飞到9000m时终止，首先设置events事件
function [value,isterminal,direction] = judge1(t,y)
%终止条件
value = y(3)-9000;%y(3)即高度h
isterminal= 1;%isterminal表示检测到指定条件时是否终止ode45函数，1终止
direction = 1;%direction表示过零点检测的方向，-1表示由负到正，+1表示由正到负
end
这时，在求解微分方程时，要在ode45中加入options，
例子：op=odeset('Events',@judge1)
[t,y]=ode45(@function,tspan,y0,op).
5.隐函数求导
隐函数用ode15i求导，给出导弹纵向运动的例子：
f=[y(5)*dy(1)-2000*cos(0.24*(k1*(y(2)-k*(y(7)-atan(-0.5083)))+k2*(dy(2)-k*dy(7))))+(0.2+0.005*((0.24*(k1*(y(2)-k*(y(7)-atan(-0.5083)))+k2*(dy(2)-k*dy(7)))*57.3)^2))*0.45*0.62475*(y(1))^2*(1-0.0065/288.15*y(4))^4.25588+y(5)*9.8*sin(y(2))
y(5)*y(1)*dy(2)-2000*sin(0.24*(k1*(y(2)-k*(y(7)-atan(-0.5083)))+k2*(dy(2)-k*dy(7))))-(0.11*(k1*(y(2)-k*(y(7)-atan(-0.5083)))+k2*(dy(2)-k*dy(7))*57.3))*0.45*0.62475*(y(1))^2*(1-0.0065/288.15*y(4))^4.25588+y(5)*9.8*cos(y(2))
dy(3)-y(1)*cos(y(2))
dy(4)-y(1)*sin(y(2))
dy(5)+0.46
dy(6)+y(1)*cos(y(7)-y(2))
y(6)*dy(7)-y(1)*sin(y(7)-y(2))];
6.传递函数与波特图
求解出导弹运动方程后，即可求出动力系数，代入传函表达式即可，注意要定义符号变量s。画波特图，例子：
logspace（-2，3，500），从10^-2到10^3,共500个点
[mag,phase],幅值mag,相位phase,角频率w
magdb,对数幅值，画对数曲线
subplot(211),共两行一列，图是第一个图
semilogx,半对数坐标函数
7.插值
如密度插值如下：
h6=[0 500 1000 1500 2000 2500 3000 3500 4000 4500 ...
5000 5500 6000 6500 7000 7500 8000 8500 9000 ...
9500 10000 11000 12000 13000 14000 15000 16000 17000 ...
18000 19000 20000 25000 30000 35000 40000 45000 50000 ...
55000 60000 65000 70000 75000 80000];
a1=[340.3 338.4 336.4 334.5 332.5 330.6 328.6 326.6 324.6 322.6 ...
320.5 318.5 316.5 314.4 312.3 310.2 308.1 306.0 303.8 301.7 ...
299.5 295.1 295.1 295.1 295.1 295.1 295.1 295.1 295.1 295.1 ...
295.1 298.4 301.7 308.3 317.2 325.8 329.8 326.7 320.6 310.1 ...
297.1 283.6 269.4];
if h>80000
a=269.4;
else
a=interp1(h6,a1,h,'spline');
end
8.求解过渡过程函数
原理：
部分代码如下：
Hwx=det([C',G(:,2),G(:,3),G(:,4)]);
syms s;
P=diff(g,s);
syms t;
wx=subs(Hwx/g,s,0);
subs的用法：将0代为表达式H/g中的s值
for i=1:4
f1=subs(Hwx/(P*s)*exp(s*t),s,x(i));%将x(i)的值赋给s
wx=wx+f1;%循环相加
end
time=0:0.1:10;
plot(time,subs(wx,t,time),'r-',time,subs(wy,t,time),'b-',time,subs(bt,t,time),'y-',time,s ubs(gm,t,time),'m-');