• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    迪恩网络公众号

四阶龙格库塔法matlab解微分方程组

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

这是我在学习飞行器制导与控制时的课程作业。用四阶龙格库塔法解微分方程组。我一开始的想法是分别利用龙格库塔法解每一个微分方程,但变量很多,算法会比较复杂。后来明白可以把多变量看作是一个变量,利用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

 


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
Delphi格式化字符串函数发布时间:2022-07-18
下一篇:
Delphi–TDataSet确定它是否在插入/编辑状态时被修改发布时间:2022-07-18
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap