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

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

Answer 1

[Alan Stevens]

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.

[Rachel Ong]

I tried your method and changed to only output rho and speedsound as they will be used in L, D,m and thrust. However, I got the error:

• Output argument "rho" (and maybe others) not assigned during call to "atmos_ode45".

• Error in Untitled_7>dynameqn (line 186)

[rho, speedsound] = atmos_ode45(height);

• Error in Untitled_7>@(t,x)dynameqn(t,x,eledefl(i),thrustvec(j)) (line 51)

[time, state_values] = ode45(@(t,x)dynameqn(t,x,eledefl(i),thrustvec(j)), tspan, initial_conditions);

• Error in ode45 (line 311)

f5 = odeFcn_main(t5, y5);

• Error in Untitled_7 (line 51)

[time, state_values] = ode45(@(t,x)dynameqn(t,x,eledefl(i),thrustvec(j)), tspan, initial_conditions);

[Stephen23]

"Output argument "rho" (and maybe others) not assigned during call to "atmos_ode45"."

Rho and others are undefined when h>h1 and h>h2.

[Rachel Ong]

But the atmos_ode45 function allows for all h values up to 20km though. What does that mean?

[Torsten]

Then most probably the height becomes larger than h2 in the course of the integration.

Output h in "atmos_ode45", and you'll see.

[Rachel Ong]

Do you happen to have any idea how I can fix this?

[Torsten]

If this is unphysical, you should check your equations for errors.

[Stephen23]

"Do you happen to have any idea how I can fix this? "

Make sure that continuous values are defined for all h.