Hi, is it possible to have variables in an ode45 that varies according to another function?

4 views (last 30 days)
This is my ode45 that I have to solve and, actually, the variables rho and speedsound (i have commented them out in the first script) is supposed to vary with altitude (x(1)). How do I allow this to work if I have another function called atmos and I want it to go back to atmos to retrieve the values of rho and speed sound for ever altitude it is before the next time step of integration?
%% Dynamic equations to do ODE45
function xdot = dynameqn(t,x,eledefl, thrustvec)
xdot = zeros(1,7);
g = 9.81;
c = 1.74;
S = 17.1;
initial_mass = 1248.5;
Ix = 1421;
Iy = 4067.5;
Iz = 4786;
Ixz = 200;
thrust_sfc = 200/60/60/1000; %kg/kN/h --> kg/N/s --> /60^2 & /1000
% rho = 1.225;
% speedsound = 340.26;
powersetting = 4;
height = x(1);
horz_displacement = x(2);
alpha = x(3);
u = x(4);
q = x(5);
theta = x(6);
mass = x(7);
CL = 6.44*alpha + 3.8*c*q/(2*u) + 0.355*eledefl; % lift coefficient
L = 0.5*rho*u^2 * CL * S; % lift force
CD = 0.03 + 0.05*CL^2; % drag coefficient
D = 0.5*rho*u^2 * CD* S; % drag force
Cm = 0.05 - 0.683*alpha - 9.96*c*q/(2*u) - 0.923*eledefl; % pitching moment coefficient
m = 0.5*rho*u^2*S*c*Cm; % pitching moment
thrust = powersetting * ( (7+u/speedsound)*200/3 + height/1000*(2*u/speedsound -11) ); % thrust force
xdot(1) = u*sin(theta-alpha); % h'= rate of vertical climb
xdot(2) = u*cos(theta-alpha); % x'= rate of change of horizontal displacement
xdot(3) = q - (thrust*sin(alpha+thrustvec) + L) / (mass* u) + g/u*cos(theta - alpha); % alpha'
xdot(4) = (thrust * cos(alpha+thrustvec) - D)/ mass - g*sin(theta-alpha); % u'
xdot(5) = m/Iy; % q'
xdot(6) = q; % theta'
xdot(7) = -thrust_sfc * thrust; % mass'
xdot=xdot';
end
atmos script below:
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [rho, speedsound] = atmos_ode45(h)
h1 = 11; % Height of tropopause
h2 = 20; % 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*1000; % 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)*1000); % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1)*1000); % Density at 20km
if h <= h1
% disp('Troposphere');
T = T0 - c*h*1000
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)*1000)
rho = rho1 * exp(-g/(R*T)*(h-h1)*1000)
speedsound = (1.4*R*T)^0.5
end
return
end

Answers (1)

Alan Stevens
Alan Stevens on 12 Jan 2022
Put something like
[T, p, rho, speedsound] = atmos(height);
immediately after
height = x(1);
in function dynameqn.
Not sure why you return T and p from atmos, as you don't seem to use them in dynameqn. Perhaps you use them elsewhere.
  7 Comments
Stephen23
Stephen23 on 12 Jan 2022
"Do you happen to have any idea how I can fix this? "
Make sure that continuous values are defined for all h.

Sign in to comment.

Categories

Find more on Dates and Time in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!