導彈運動方程求解數值解離不開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-');
?

總結

以上是生活随笔為你收集整理的MATLAB求解导弹运动的一些基础方法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。