Is this how I find the area under my graph?

4 次查看（过去 30 天）
Rachel Ong
Rachel Ong 2022-3-10
评论： Rachel Ong 2022-3-11
I have this code which takes too long to run so I would save you the trouble and just attach a snippet of how I plotted the graph and called the trapz for finding the area. However, I did this but my answer seems quite far off from my caluclations so I am wondering if this was the right method to find the area under the graph. I'd like to know if I did something wrong or this answer is correct.
So this is the graph I generated and I'll attach the code below.
%% Maximum ROC vs V by taking the maximum ROC at each altitude
% ------ Functions to solve for FSOLVE (steady climbing flight) ------ %
% func1 = v' = 1/m *(T-D-mgsin(gamma)) = 0
% func2 = gamma' = 1/mu * (L-mgcos(gamma)) = 0
% func3 = Cm = 0
% ------ FSOLVE variables ------ %
% x(1) = velocity
% x(2) = deltae
% x(3) = gamma
% ------ Altitude range from flight envelope ------ %
hlist_2_5 = linspace(0,9000,300);
hlist_3 = linspace(0,17000,300);
% ------ Array for storing ------ %
% 1: V 2:deltae 3:gamma 4:aoa 5:roc 6. altitude 7. thrust 8. 1/ROC
final_list_2_5 = zeros(length(hlist_2_5),8);
final_list_3 = zeros(length(hlist_3),8);
%% PS = 2.5
%% --REMOVED-- the calculations
for j = 1:length(hlist_2_5)
[max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
final_list_2_5(j,5) = max_roc; % ROC
final_list_2_5(j,1) = store_v_de_gam_aoa_roc(index,1); % V
final_list_2_5(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
final_list_2_5(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
final_list_2_5(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
final_list_2_5(j,6) = hlist_2_5(j) % h
final_list_2_5(j,7) = store_v_de_gam_aoa_roc(index,6) % thrust
final_list_2_5(j,8) = 1/max_roc; % 1/ROC
end
%% PS = 3
%% --REMOVED-- the calculations
for j = 1:length(hlist_3)
[max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
final_list_3(j,5) = max_roc; % ROC
final_list_3(j,1) = store_v_de_gam_aoa_roc(index,1); % V
final_list_3(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
final_list_3(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
final_list_3(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
final_list_3(j,6) = hlist_3(j) % h
final_list_3(j,7) = store_v_de_gam_aoa_roc(index,6) % thrust
final_list_3(j,8) = 1/max_roc; % 1/ROC
end
% ------ Finding time to climb to each altitude ------ %
time_2_5_climb = trapz(final_list_2_5(:,8),final_list_2_5(:,6))
% 92475.4221575 s = 25.7 hours <-- this is the result I got
time_3_climb = trapz(final_list_3(:,8),final_list_3(:,6))
% 164693.805554 s = 45.7 hours <-- this is the result I got
% ------ Plotting graph (X vs Y)------ %
figure % 1/ROC vs Altitude
plot(final_list_2_5(:,6),final_list_2_5(:,8))
hold on
plot(final_list_3(:,6),final_list_3(:,8))
xlabel('Altitude [m]')
ylabel('1/Rate of climb [s/m]')
title('1/Rate of climb vs Altitude')
legend('PS = 2.5','PS = 3')
grid on
  1 个评论
Rachel Ong
Rachel Ong 2022-3-10
编辑：Rachel Ong 2022-3-10
Incase the snippet of script didnt make sense, here is the full code. However, it takes almost 600s to run.
%% Maximum ROC vs V by taking the maximum ROC at each altitude
clc
clear all
% ------ Constants ------ %
m = 1248.5;
g = 9.81;
W = m*g;
S = 17.1;
% ------ Functions to solve for FSOLVE (steady climbing flight) ------ %
% func1 = v' = 1/m *(T-D-mgsin(gamma)) = 0
% func2 = gamma' = 1/mu * (L-mgcos(gamma)) = 0
% func3 = Cm = 0
% ------ FSOLVE variables ------ %
% x(1) = velocity
% x(2) = deltae
% x(3) = gamma
% ------ Altitude range from flight envelope ------ %
hlist_2_5 = linspace(0,9000,300);
hlist_3 = linspace(0,17000,300);
% ------ Array for storing ------ %
% 1: V 2:deltae 3:gamma 4:aoa 5:roc 6. altitude 7. thrust 8. 1/ROC
final_list_2_5 = zeros(length(hlist_2_5),8);
final_list_3 = zeros(length(hlist_3),8);
% ------ Initial guess ------ %
x = [40 0 0];
%% PS = 2.5
for j = 1:length(hlist_2_5)
PS = 2.5
h = hlist_2_5(j)
[rho, speedsound] = atmos(h)
[aoa_start,aoa_end] = aoa_guess(PS)
aoalist = linspace(aoa_start,aoa_end,150);
store_v_de_gam_aoa_roc = zeros(length(aoalist),6);
% 1: V 2:deltae 3:gamma 4:aoa 5:roc 6. thrust
for i = 1:length(aoalist)
qbar = @(x) 0.5*rho*x(1)^2*S;
CL = @(x) 6.44*aoalist(i) + 0.355*x(2);
L = @(x) CL(x)*qbar(x);
CD = @(x) 0.03 + 0.05*(CL(x))^2;
D = @(x) CD(x)*qbar(x);
Cm = @(x) 0.05 - 0.683*aoalist(i) - 0.923*x(2);
T = @(x) PS *( (7+x(1)/speedsound)*200/3 + h/1000*(2*x(1)/speedsound -11));
f1 = @(x) (1/m)*(T(x)-D(x)-m*g*sin(x(3)));
f2 = @(x) (1/(m*x(1)))*(L(x)-m*g*cos(x(3)));
f3 = @(x) Cm(x);
func = @(x) [f1(x);f2(x);f3(x)];
x = fsolve(@(x) func(x), x);
store_v_de_gam_aoa_roc(i,1) = x(1); % V
store_v_de_gam_aoa_roc(i,2) = x(2); % delta e
store_v_de_gam_aoa_roc(i,3) = x(3); % gamma
store_v_de_gam_aoa_roc(i,4) = aoalist(i); % aoa
thrust = PS * ( 200/3*(7 + x(1)/speedsound) + h/1000*(2*x(1)/speedsound - 11) );
liftcoeff = 6.44*aoalist(i) + 0.355*x(2);
drag = 0.5*rho*x(1)^2*S*(0.03 + 0.05*liftcoeff^2);
ROC = (thrust-drag)*x(1)/W;
store_v_de_gam_aoa_roc(i,5)= ROC; % ROC
store_v_de_gam_aoa_roc(i,6)= thrust; % ROC
end
[max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
final_list_2_5(j,5) = max_roc; % ROC
final_list_2_5(j,1) = store_v_de_gam_aoa_roc(index,1); % V
final_list_2_5(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
final_list_2_5(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
final_list_2_5(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
final_list_2_5(j,6) = hlist_2_5(j) % h
final_list_2_5(j,7) = store_v_de_gam_aoa_roc(index,6) % thrust
final_list_2_5(j,8) = 1/max_roc; % 1/ROC
end
% ------ Initial guess ------ %
x = [40 0 0];
%% PS = 3
for j = 1:length(hlist_3)
PS = 3;
h = hlist_3(j);
[rho, speedsound] = atmos(h);
[aoa_start,aoa_end] = aoa_guess(PS);
aoalist = linspace(aoa_start,aoa_end,150);
store_v_de_gam_aoa_roc = zeros(length(aoalist),5);
% 1: V 2:deltae 3:gamma 4:aoa 5:roc 6. thrust
for i = 1:length(aoalist)
qbar = @(x) 0.5*rho*x(1)^2*S;
CL = @(x) 6.44*aoalist(i) + 0.355*x(2);
L = @(x) CL(x)*qbar(x);
CD = @(x) 0.03 + 0.05*(CL(x))^2;
D = @(x) CD(x)*qbar(x);
Cm = @(x) 0.05 - 0.683*aoalist(i) - 0.923*x(2);
T = @(x) PS *( (7+x(1)/speedsound)*200/3 + h/1000*(2*x(1)/speedsound -11));
f1 = @(x) (1/m)*(T(x)-D(x)-m*g*sin(x(3)));
f2 = @(x) (1/(m*x(1)))*(L(x)-m*g*cos(x(3)));
f3 = @(x) Cm(x);
func = @(x) [f1(x);f2(x);f3(x)];
x = fsolve(@(x) func(x), x);
store_v_de_gam_aoa_roc(i,1) = x(1); % V
store_v_de_gam_aoa_roc(i,2) = x(2); % delta e
store_v_de_gam_aoa_roc(i,3) = x(3); % gamma
store_v_de_gam_aoa_roc(i,4) = aoalist(i); % aoa
thrust = PS * ( 200/3*(7 + x(1)/speedsound) + h/1000*(2*x(1)/speedsound - 11) );
liftcoeff = 6.44*aoalist(i) + 0.355*x(2);
drag = 0.5*rho*x(1)^2*S*(0.03 + 0.05*liftcoeff^2);
ROC = (thrust-drag)*x(1)/W;
store_v_de_gam_aoa_roc(i,5)= ROC; % ROC
store_v_de_gam_aoa_roc(i,6)= thrust; % thrust
end
[max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
final_list_3(j,5) = max_roc; % ROC
final_list_3(j,1) = store_v_de_gam_aoa_roc(index,1); % V
final_list_3(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
final_list_3(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
final_list_3(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
final_list_3(j,6) = hlist_3(j) % h
final_list_3(j,7) = store_v_de_gam_aoa_roc(index,6) % thrust
final_list_3(j,8) = 1/max_roc; % 1/ROC
end
% ------ Finding time to climb to each altitude ------ %
time_2_5_climb = trapz(final_list_2_5(:,8),final_list_2_5(:,6))
% 92475.4221575 s = 25.7 hours
time_3_climb = trapz(final_list_3(:,8),final_list_3(:,6))
% 164693.805554 s = 45.7 hours
% ------ Plotting graph (X vs Y)------ %
figure % 1/ROC vs Altitude
plot(final_list_2_5(:,6),final_list_2_5(:,8))
hold on
plot(final_list_3(:,6),final_list_3(:,8))
xlabel('Altitude [m]')
ylabel('1/Rate of climb [s/m]')
title('1/Rate of climb vs Altitude')
legend('PS = 2.5','PS = 3')
grid on
figure % Velocity vs Altitude
plot(final_list_2_5(:,6),final_list_2_5(:,1))
hold on
plot(final_list_3(:,6),final_list_3(:,1))
xlabel('Altitude [m]')
ylabel('Velocity [m/s]')
title('Velocity vs Altitude')
legend('PS = 2.5','PS = 3')
grid on
%% ------ Atmos function ------ %%
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [rho, speedsound] = atmos(h)
h1 = 11000; % Height of tropopause
h2 = 20000; % End height of table
g = 9.81;
R = 287;
c = 6.51e-3; % temperature lapse dt/dh = - c = -6.51 degcelcius/km
T0 = 15+273.15; % Temperature sea level
p0 = 101325; % pressure sealevel
rho0 = 101325/R/T0; % density sealevel = pressure / R*T, R=287, T = 15 degcelcius
T1 = T0 - c*h1; % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506; % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506; % Density at 11km
T2 = T1; % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1)); % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1)); % Density at 20km
if h <= h1
% disp('Troposphere');
T = T0 - c*h;
p = p0 * (T/T0)^5.2506;
rho = rho0 * (T/T0)^4.2506;
speedsound = (1.4*R*T)^0.5;
elseif h <= h2
% disp('Tropopause');
T = T1;
p = p1 * exp(-g/(R*T)*(h-h1));
rho = rho1 * exp(-g/(R*T)*(h-h1));
speedsound = (1.4*R*T)^0.5;
end
return
end
%% ------ AOA guess function ------ %%
% To guess the range of AOA at each altitude
function [aoa_start,aoa_end] = aoa_guess(PS)
h = 0;
W = 1248.5*9.81; % Weight of UAV
S = 17.1; % Wing area
% CL = (6.44*aoa + 0.355*eledefl);
% CD = (0.03 + 0.05*(6.44*aoa + 0.355*eledefl)^2);
% Cm = (0.05 - 0.683*aoa - 0.923*eledefl);
% thrust = (PS * ( (7 + V/speedsound )*200/3 + h*(2*(V/speedsound) - 11) ));
%%
% V = x(1);
% aoa = x(2);
% eledefl = x(3);
[rho, speedsound] = atmos(h);
qbar = @(x) 0.5*rho*x(1)^2*S;
CL = @(x) (6.44*x(2) + 0.355*x(3));
CD = @(x) (0.03 + 0.05*(CL(x))^2);
Cm = @(x) 0.05 - 0.683*x(2) - 0.923*x(3);
thrust = @(x) ( PS * ( 200/3*( 7 + x(1)/speedsound ) + h/1000*( 2*x(1)/speedsound - 11) ) );
func_1 = @(x) qbar(x)*CL(x) + thrust(x)*sin(x(2)) - W;
func_2 = @(x) qbar(x)*CD(x) - thrust(x)*cos(x(2));
func_3 = @(x) Cm(x);
FUNC = @(x) [func_1(x); func_2(x); func_3(x)];
% ------ Guess --> Range of AOA ------ %
Xlow = fsolve(@(x) FUNC(x), [10 0 -30/180*pi]); % high AOA
aoa_end = Xlow(1,2);
Xhigh = fsolve(@(x) FUNC(x), [100 25/180*pi -30/180*pi]); % low AOA
aoa_start = Xhigh(1,2);
end

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-3-10
编辑：Torsten 2022-3-10
% ------ Finding time to climb to each altitude ------ %
time_2_5_climb = trapz(final_list_2_5(:,6),final_list_2_5(:,8))
% 92475.4221575 s = 25.7 hours <-- this is the result I got
time_3_climb = trapz(final_list_3(:,6),final_list_3(:,8))
% 164693.805554 s = 45.7 hours <-- this is the result I got
instead of
% ------ Finding time to climb to each altitude ------ %
time_2_5_climb = trapz(final_list_2_5(:,8),final_list_2_5(:,6))
% 92475.4221575 s = 25.7 hours <-- this is the result I got
time_3_climb = trapz(final_list_3(:,8),final_list_3(:,6))
% 164693.805554 s = 45.7 hours <-- this is the result I got
And one further hint:
fsolve will perform much faster in your loops if you choose the last solution as initial guess for the next call:
x = fsolve(@(x) func(x), x);
instead of
x = fsolve(@(x) func(x), [40 0 0]);
Dont forget to initialize
x = [40 0 0]
in front of the two loops.
  8 个评论
显示 6更早的评论
Torsten
Torsten 2022-3-11
You mean
f = @(t)interp1(final_list_2_5(:,6),final_list_2_5(:,1),t);
val = integral(f,final_list_2_5(1,6),x_1000)
?
To see how the area under the curves develop with altitude, you may also use "cumtrapz":
S = cumtrapz(final_list_2_5(:,6),final_list_2_5(:,8));
plot(final_list_2_5(:,6),S)
Rachel Ong
Rachel Ong 2022-3-11
I see, I'll try that and check out cumtrapz as well. Thank you!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by