MATLAB的答案

指数在位置2超出数组边界(不得超过1)。

18视图(30天)
奥Ataseven”class=
奥Ataseven 2022年1月4日
编辑: Torsten 2022年1月5日
嗨。首先很抱歉问同一个问题一次又一次但我真的接近解决方案(可能)。
请不要介意排长队,因为大部分都是常量和循环。
这段代码试图解决6常微分方程与6状态变量(水平位置(x1和x2),高度(x3),真正的空速(x4),航向角(x5)和飞机的质量(x6)]和3控制输入(发动机推力(u1)、银行角度(u2)和飞行路线角(u3)]用欧拉方法。
不同的飞行动作执行指定的时间间隔。
速度。米,Cruise_Vel。米,Des_Vel。米,Thr_cl。米,Thr_cr。米,Thr_des。米,fuel_cl。米,fuel_cr。m, fuel_des.m、den.m den_cr.m、den_des.m drag.m, drag_cr.m, drag_des.m, lift.m, lift_cr.m lift_des。m是功能分离选项卡。
主要代码:
%从h1 = 1100 [m] h2与α= 5 = 1600 [m]飞行路线角。
%为t = 60分钟执行巡航飞行。
%与β= 30银行角度,直到转航向由η= 270◦改变。
从h2 = 1600 [m] %下降h1和ζ= 4 = 1100 [m]◦飞行路线角。
%完成360◦在水平飞行(徘徊)转。
%下降和κh3 = 800 [m] = 4.5◦飞行路线角。
%飞机属性
W = .44225E + 06;% .44225E + 3吨= .44225E + 6公斤
S = .51097E + 03;%的表面积(m ^ 2)
g0 = 9.80665;%重力加速度(m / s2)
%使用数值方法求解一阶的颂歌
t0 = 0;
往往= 3960;
h = 0.05;
N = (tend-t0) / h;
t = t0: h:一般;
%预先配置
x = 0(6、长度(t));
x1 = 0(1、长度(t));
x2 = 0(1、长度(t));
x3 = 0(1、长度(t));
x4 = 0(1、长度(t));
x5 = 0(1、长度(t));
x6 = 0(1、长度(t));
u1 = 0(1、长度(t));
u2 = 0(1、长度(t));
u3 = 0(1、长度(t));
重金属镉= 0(1、长度(t));
p = 0(1、长度(t));
Cl = 0(1、长度(t));
f = 0(1、长度(t));
dx1dt = 0(1、长度(t));
dx2dt = 0(1、长度(t));
dx3dt = 0(1、长度(t));
dx4dt = 0(1、长度(t));
dx5dt = 0(1、长度(t));
dx6dt = 0(1、长度(t));
%初始条件
x (: 1) = (0, 0, 3608.92; 1.0 e + 2 * 1.161544478045788; 0; W);
为我= 2:长度(t)
如果和(t(1张)> = 0,t(1张)< 60)%从h1 = 1100 [m] h2与α= 5 = 1600 [m]飞行路线角。
x3 = linspace (3608.92, 5249.3, 79201);
x4 =速度(x3);%改变速度(米/秒)
x5 = 0;%改变头部角(度)
f = fuel_cl (x3);%改变燃料流量(公斤/分钟)
u1 = Thr_cl (x3);%改变推力[N]
u2 = 0;%改变银行角度(度)
u3 = 5;%改变飞行路线角(度)
V_ver = x4 * sin (u3);%改变垂直速度(米/秒)
重金属镉=阻力(x3, x4);%改变阻力系数
Cl =电梯(x3, x4);%改变升力系数
p =窝(x3);%改变密度(公斤/立方米)
elseif和(t(1张)> = 60 t(1张)< 3660)%为t = 60分钟执行巡航飞行。
x3 = 5249.3;
x4 = Cruise_Vel (x3);%改变速度(米/秒)
x5 = 0;%改变头部角(度)
f = fuel_cr (x3);%改变燃料流量(公斤/分钟)
u1 = Thr_cr (x3);%改变推力[N]
u2 = 0;%改变银行角度(度)
u3 = 0;%改变飞行路线角(度)
V_ver = x4 * sin (u3);%改变垂直速度(米/秒)
重金属镉= drag_cr (x3, x4);%改变阻力系数
Cl = lift_cr (x3, x4);%改变升力系数
p = den_cr (x3);%改变密度(公斤/立方米)
elseif和(t(1张)> = 3660,t(1张)< 3720)%与β= 30银行角度,直到转航向由η= 270◦改变。
x3 = 5249.3;
x4 = Cruise_Vel (x3);%改变速度(米/秒)
x5 = 0:30:270;%改变头部角(度)
f = fuel_cr (x3);%改变燃料流量(公斤/分钟)
u1 = Thr_cr (x3);%改变推力[N]
u2 = 30;%改变银行角度(度)
u3 = 0;%改变飞行路线角(度)
V_ver = x4 * sin (u3);%改变垂直速度(米/秒)
重金属镉= drag_cr (x3, x4);%改变阻力系数
Cl = lift_cr (x3, x4);%改变升力系数
p = den_cr (x3);%改变密度(公斤/立方米)
elseif和(t(1张)> = 3720,t(1张)< 3780)从h2 = 1600 [m] %下降h1和ζ= 4 = 1100 [m]◦飞行路线角。
x3 = linspace (5249.3, 3608.92, 79201);
x4 = Des_Vel (x3);%改变速度(米/秒)
x5 = 270;%改变头部角(度)
f = fuel_des (x3);%改变燃料流量(公斤/分钟)
u1 = Thr_des (x3);%改变推力[N]
u2 = 0;%改变银行角度(度)
u3 = 4;%改变飞行路线角(度)
V_ver = x4 * sin (u3);%改变垂直速度(米/秒)
重金属镉= drag_des (x3, x4);%改变阻力系数
Cl = lift_des (x3, x4);%改变升力系数
p = den_des (x3);%改变密度(公斤/立方米)
elseif和(t(1张)> = 3780,t(1张)< 3900)%完成360◦在水平飞行(徘徊)转。
x3 = 3608.9;
x4 = Cruise_Vel (x3);%改变速度(米/秒)
朗= (60 120 180 240 270 270 300 360);
x5 = wrapTo360(朗);%改变头部角(度)
f = fuel_cr (x3);%改变燃料流量(公斤/分钟)
u1 = Thr_cr (x3);%改变推力[N]
u2 = 0;%改变银行角度(度)
u3 = 0;%改变飞行路线角(度)
V_ver = x4 * sin (u3);%改变垂直速度(米/秒)
重金属镉= drag_cr (x3, x4);%改变阻力系数
Cl = lift_cr (x3, x4);%改变升力系数
p = den_cr (x3);%改变密度(公斤/立方米)
elseif和(t(1张)> = 3900,t(1张)< 3960)%下降和κh3 = 800 [m] = 4.5◦飞行路线角。
x3 = linspace (3608.92, 2624.67, 79201);
x4 = Des_Vel (x3);%改变速度(米/秒)
x5 = 270;%改变头部角(度)
f = fuel_des (x3);%改变燃料流量(公斤/分钟)
u1 = Thr_des (x3);%改变推力[N]
u2 = 0;%改变银行角度(度)
u3 = 4.5;%改变飞行路线角(度)
V_ver = x4 * sin (u3);%改变垂直速度(米/秒)
重金属镉= drag_des (x3, x4);%改变阻力系数
Cl = lift_des (x3, x4);%改变升力系数
p = den_des (x3);%改变密度(公斤/立方米)
其他的
流(“一个问题发生。”);
结束
dx1dt = x4。* cos (x5)。* cos (u3);
dx2dt = x4。* sin (x5)。* cos (u3);
dx3dt = x4。* sin (u3);
dx4dt =重金属镉。* * p。* (x4 ^ 2)。/ (2。* x6) g0。* sin (u3) + u1. / x6;
dx5dt = cl。* * p。* x4. / (2 * x6)。*罪(u2);
dx6dt = - f;
(我)= x(1张)+ h * dx1dt(1张);% % % % % % % % %线138% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
(我)= x(2张)+ h * dx2dt(1张);
(我)= x(3张)+ h * dx3dt(1张);
(我)= x(4张)+ h * dx4dt(1张);
(我)= x(5张)+ h * dx5dt(1张);
我x(6日)= x(6张)+ h * dx6dt(1张);
结束
合计= cell2mat (f);%的总燃料消耗量在任务(公斤/分钟)
Tot_fuel =金额(合计);
图(1)
plot3 (x1 (:), x2 (:), x3 (:));% 3 d图形位置
图(2)
情节(t, x4 (:));% vta−时间图
图(3)
情节(t, V_ver (:));% V_vertical−时间图
图(4)
情节(t, x5 (:));%标题−时间图
图(5)
情节(t, x6 (:));%质量−时间图
图(6)
情节(t, u1 (:));%推力−时间图
图(7)
情节(t, u2 (:));%的倾斜角−时间图
图(8)
情节(t, u3 (:));%飞行路线角−时间图
流(在任务是%的总燃料消耗量。2 f(公斤),Tot_fuel * / 60);
为什么我用79201大小的数组长度(t) = 79201。
当我运行:
指数在位置2超出数组边界(不得超过1)。
错误在论坛(第138行)
(我)= x(1张)+ h * dx1dt(1张);
我应该做什么?谢谢。
下面的函数是独立的标签之一,休息是类似的:
函数[Vtas_cl] =速度(x3)
% % % % % % % % % % % % % % % % % % % %常量% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
Vcl_1 = 335;%的标准校准空速(kt)
Vcl_2 = 172.3;%的标准校准空速(kt) - > (m / s)(找到马赫过渡高度)
Vcl_2_in_knots = 335;%标准校准空速(kt)(在海里找到结果,如果海拔10000英尺和马赫数之间的过渡高度)
M_cl = 0.86;%的标准校准空速(kt)
K = 1.4;%的空气绝热指数
R = 287.05287;%真实气体常数空气(平方米/ (K·s2))
Bt = - 0.0065;% ISA温度梯度与对流层顶高度低于(K / m)
T0 = 288.15;%标准大气温度与韩剧[K]
g0 = 9.80665;%重力加速度(m / s2)
a0 = 340.294;%声速(米/秒)
Vd_CL1 = 5;%上升速度增量低于1500英尺(飞机)
Vd_CL2 = 10;%上升速度增量低于3000英尺(飞机)
Vd_CL3 = 30;%上升速度增量低于4000英尺(飞机)
Vd_CL4 = 60;%上升速度增量低于5000英尺(飞机)
Vd_CL5 = 80;%上升速度增量低于6000英尺(飞机)
CV_min = 1.3;%的最低速度系数
Vstall_TO = .14200E + 03;%失速速度起飞(王者文化)
CAS_climb = Vcl_2;
Mach_climb = M_cl;
delta_trans = (((1 + ((K - 1) / 2) * (CAS_climb / a0) ^ 2) ^ (K / (K - 1))) 1) / (((1 + (K - 1) / 2 * Mach_climb ^ 2) ^ (K / (K - 1))) 1);%的压力比过渡高度
teta_trans = delta_trans ^ (bt * R / g0);%在转变温度比高度
H_p_trans_climb = (1000/0.348/6.5) * (T0 * (1-teta_trans));%转换高度攀登(英尺)
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %的常数
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
H_climb = x3;输入% % % % % % % % % % % % % % % % % % % % % % % % %
Vnom_climb_jet = 0(1、长度(H_climb));
为k1 = 1:长度(H_climb)
如果(0 < = H_climb (k1) &&H_climb (k1) < 1500)
Vnom_climb_jet (k1) = CV_min * Vstall_TO + Vd_CL1;
elseif(1500 < = H_climb (k1) &&H_climb (k1) < 3000)
Vnom_climb_jet (k1) = CV_min * Vstall_TO + Vd_CL2;
elseif(3000 < = H_climb (k1) &&H_climb (k1) < 4000)
Vnom_climb_jet (k1) = CV_min * Vstall_TO + Vd_CL3;
elseif(4000 < = H_climb (k1) &&H_climb (k1) < 5000)
Vnom_climb_jet (k1) = CV_min * Vstall_TO + Vd_CL4;
elseif(5000 < = H_climb (k1) &&H_climb (k1) < 6000)
Vnom_climb_jet (k1) = CV_min * Vstall_TO + Vd_CL5;
elseif(6000 < = H_climb (k1) &&H_climb (k1) < 10000)
Vnom_climb_jet (k1) = min (Vcl_1,250);
elseif(10000 < = H_climb (k1) &&H_climb (k1) < = H_p_trans_climb)
Vnom_climb_jet (k1) = Vcl_2_in_knots;
elseif(H_p_trans_climb < H_climb (k1))
Vnom_climb_jet (k1) = M_cl;
结束
Vcas_cl (k1) = Vnom_climb_jet (k1) * 0.514;% (kn) - >(米/秒)
H_climb (k1) = H_climb (k1) * 0.3048;%[脚]- > [m]
K = 1.4;%的空气绝热指数
R = 287.05287;%真实气体常数空气(平方米/ (K·s2))
Bt = - 0.0065;% ISA温度梯度与对流层顶高度低于(K / m)
deltaT = 0;%的价值真正的温度T ISA条件[K]
T0 = 288.15;%标准大气温度与韩剧[K]
P0 = 101325;%火星科学实验室标准大气压力(Pa)
g0 = 9.80665;%重力加速度(m / s2)
p0 = 1.225;%的标准大气密度与韩剧(公斤/立方米)
visc = (K - 1) / K;
T (k1) = T0 + deltaT + Bt * H_climb (k1);%温度[K]
P (k1) = P0 * ((T (k1) -deltaT) / T0)。^ ((g0) / (Bt * R));% (Pa)的压力
p (k1) = p (k1)。/ (R * T (k1));%密度(公斤/ m ^ 3)
Vtas_cl (k1) = (2 * P (k1) / visc / P (k1) * ((1 + P0 / P (k1) * ((1 + visc * P0 * Vcas_cl (k1) * Vcas_cl (k1) / 2 / P0)。^ (1 / visc) 1))。^ (visc) 1))。^ (1/2);%真实的空气速度(米/秒)
结束
%输出
结束
8的评论
显示 7年长的评论
Torsten”class=
Torsten 2022年1月4日
编辑:Torsten 2022年1月5日
我不明白你的问题。
你不需要拯救dx1dt,…, dx6dt。所以不需要把它们写在数组dx1dt(1张),…,dx6dt(1张)。他们可以是标量。

登录置评。

在回答这个问题。

接受的答案

1月”class=
1月 2022年1月4日
编辑:1月 2022年1月4日
dx1dt = x4。* cos (x5)。* cos (u3);%现在dx1dt是一个标量
…
(我)= x(1张)+ h * dx1dt(1张);%这将dx1dt视为一个向量
我猜,你想要的:
dx1dt (i) = x4。* cos (x5)。* cos (u3);
% ^ ^ ^ ^ ^
顺便说一下:
如果和(t(1张)> = 0,t(1张)< 60)
…
elseif和(t(1张)> = 60 t(1张)< 3660)
可以simpilified。所有t > = 0和第一个测试后,t < 60已经排除:
如果t(1张)< 60
…
elseift(1张)< 3660
有更多的问题,例如:
plot3 (x1 (:), x2 (:), x3 (:));
x1是预先分配,但现在值被应用。x3多次重新定义。
1评论
奥Ataseven”class=
奥Ataseven 2022年1月5日
谢谢你的回应。
” x3重新定义了很多次了。” 打开了我的眼睛,我能说的。

登录置评。

更多的答案(0)

类别

找到更多的在制导、导航和控制(GNC)在帮助中心和文件交换

标签

s manbetx 845


释放

R2021a

社区寻宝

找到宝藏在MATLAB中央,发现社区如何帮助你!

开始狩猎!