四階龍格庫塔法matlab解微分方程組

這是我在學習飛行器制導與控制時的課程作業。用四階龍格庫塔法解微分方程組。我一開始的想法是分別利用龍格庫塔法解每一個微分方程，但變量很多，算法會比較復雜。后來明白可以把多變量看作是一個變量，利用matlab的feval函數進行代入變量的函數運算。

matlab中feval函數的作用：feval(f,x,y)；將x，y代入函數f中。

四階龍格-庫塔法：

需要解的四個微分方程組為：

 

算法代碼：

主函數main.m

main.m
clear all;
clc;
%% 定義初始參數
M0=2;            %馬赫數
h0=5000;       %初始高度
theta0=-30*pi/180;    %初始彈道傾角，注意度數的表示
a=340;            %音速
g0=9.81;
%% 運算
v0=M0*a;
y0=[v0;theta0;0;h0];
h=0.01;   %步長
t0=0;
t1=200;
[finaltime,n,y,t] = RK4(@equations,y0,h,t0,t1);
%% 顯示結果
title('彈道曲線');
plot(y(3,1:n),y(4,1:n));
xlabel('x/m');ylabel('h/m');
figure;
title('速度變化曲線');
plot(t(1:n),y(1,1:n));
xlabel('t/s');ylabel('v/(m/s)');
figure;

四階龍格-庫塔法RK4.m

function [finaltime,n,y,t] = RK4(f,y0,h,t0,t1)
%龍格庫塔四階
%   f:微分方程組
%   y0:[v;theta;x;h]初始值
%   h:步長
%   t0-t1:時間區間

t=t0:h:t1;
y=zeros(length(y0),length(t));
y(:,1)=y0;
flag=1;
n=1;
while(flag)           %判斷導彈是否落地
    if (y(4,n)>0)
        k1=feval(f,t(n),y(:,n));
        k2=feval(f,t(n)+h/2,y(:,n)+h/2*k1);
        k3=feval(f,t(n)+h/2,y(:,n)+h/2*k2);
        k4=feval(f,t(n)+h,y(:,n)+h*k3);
        y(:,n+1)=y(:,n)+(h/6)*(k1+2*k2+2*k3+k4);
        n=n+1;
    else
        finaltime=t(n-1);
        flag=0;
     end
end

end

導彈質點運動的動力學微分方程組equations.m

function f= equations(t,y)
% 四個微分方程組
%   t為時間
%   y=[v;theta;x;h]
m=260;         %質量
s=0.24;           %參考面積
a=340;            %音速
alpha=2;       %攻角
g0=9.81;

M=y(1)/a;
G=m*g0;
rho=1.225*exp(-0.00015*y(4));

Cx=[M^2,M,1]*[0.0002, 0.0038, 0.1582;
                           -0.0022, -0.0132, -0.8520;
                           0.0115, -0.0044, 1.9712;]*[alpha^2,alpha,1]';
X=Cx*1/2*rho*(y(1)^2)*s;
Cy=[M^2,M,1]*[-0.026,0.0651,0.4913]'*alpha;
Y=Cy*1/2*rho*(y(1)^2)*s;

f(1)=-X/m-G*sin(y(2))/m;
f(2)=Y/(m*y(1))-G*cos(y(2))/(m*y(1));
f(3)=y(1)*cos(y(2));
f(4)=y(1)*sin(y(2));
f=f(:);
end