導彈運動方程求解數值解離不開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);%矩陣求值，即特徵方程式，使用

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-');