Using Ode45 to solve dynamics problem (ISA model)

How can I express a variable density in diferentialfuntion2? I have modeled athmospheric density in complex_atm_model in terms of height (which is position in the one dimensioinal motion model) but I am having trouble when trying to solve using ode45. Any ideas?

z0 = const.h0;
v0 = const.v0;
t0 = 0;
tf = 800;
N = 60000;
t = linspace(t0, tf, N);                        
X = [z0;v0];
[t,X] = ode45('diferentialfunction2', t, X);

------------------------------------------------------- Main code

function dXdt = diferentialfunction2(t, X)
c=2;
s=0.3;
g=9.81;
m = 100;
G = 6.67E-11;
R = 6371E3; 
M = 5.9722E24;
h = 39045:-1:0;
h(39046) = 0;
[rho, T] = complex_atm_model(h,X);
dXdt(1,1) =  X(2);     
dXdt(2,1) = ((rho*c*s*((X(2))^2))/(2*m))-(G*(M/(R+X(1))^2));     
end

---------------------------------------------------------------------------------

function [density, T] = complex_atm_model (h,X)
% Parameters to be used:
T_SL   = 288.16;      % Sea level temperature [K]
rho_SL = 1.225;       % Sea level density [kg/m3]
Rg     = 287.0;       % Gas constant for the air [J/kg].
dT_dh  = -0.0065;     % Gradient of temperature in the troposphere [K/m]
g0     = 9.81;        % Gravity acceleration at sea level [m/s2]
%Temperature model for all altitudes:
T = (T_SL * h ./ h) + (dT_dh*h);                    %Temperature vector modeled constant
inds = find(h) > 11000;                             %Find indexes for altitudes > 11000 m
T(inds) = T(11000);                                 %Set constant T = T(11000) for h > 11000
%Density model for all altitudes:
inds_low = find(h < 11000);                                %Define inds_low as indexes for altitudes < 11000
rho = rho_SL * h./h;                                       %Prelocate a density vector with uniform density
rho(inds_low) = (rho_SL * (T(inds_low)/T_SL).^4.26);       %Density vector modeled constant
exp_athmos = (Rg * T(11000))/g0;                           %Measure how density decays with altitude
for i=11000:39045
    rho(i) = rho(11000)*exp(-(h(i)-11000)/exp_athmos);      %Set density values for h > 11000
end
density = rho(X(1));
end

[Alan Stevens]

Upload the code, not an image (no-one is going to re-write it all). Also upload the code for complex-atm_model to allow the model to be tested.

[JXT119]

Thanks for your comment, just did it.

[Walter Roberson]

What problem are you encountering?

z0 = const.h0;
v0 = const.v0;

we will need to know those in order to test.

Answer 1

[Alan Stevens]

Like this

z0 = 39045; %const.h0;
v0 = 0;     %const.v0;  
t0 = 0;
tf = 800;
N = 60000;

tspan = linspace(t0, tf, N);
X = [z0;v0];
[t,X] = ode45(@diferentialfunction2, tspan, X);

plot(t,X(:,1)),grid
xlabel('time'), ylabel('height')

function dXdt = diferentialfunction2(~, X)
c=2;
s=0.3;
% g=9.81;
m = 100;
G = 6.67E-11;
R = 6371E3; 
M = 5.9722E24;
h = 0:39045;

rho = complex_atm_model(h,X);
dXdt(1,1) =  X(2);     
dXdt(2,1) = ((rho*c*s*((X(2))^2))/(2*m))-(G*(M/(R+X(1))^2));     
end

function density = complex_atm_model (h,X)
% Parameters to be used:
T_SL   = 288.16;      % Sea level temperature [K]
rho_SL = 1.225;       % Sea level density [kg/m3]
Rg     = 287.0;       % Gas constant for the air [J/kg].
dT_dh  = -0.0065;     % Gradient of temperature in the troposphere [K/m]
g0     = 9.81;        % Gravity acceleration at sea level [m/s2]
%Temperature model for all altitudes:
T = T_SL + dT_dh*h;                    %Temperature vector modeled constant
inds = 11000:max(h);                             %Find indexes for altitudes > 11000 m
T(inds) = T(11000);                                 %Set constant T = T(11000) for h > 11000
%Density model for all altitudes:
inds_low = 1:11000;                                %Define inds_low as indexes for altitudes < 11000
%rho = rho_SL;                                       %Prelocate a density vector with uniform density
rho(inds_low) = (rho_SL * (T(inds_low)/T_SL).^4.26);       %Density vector modeled constant
exp_athmos = (Rg * T(11000))/g0;                           %Measure how density decays with altitude
for i=11000:max(h)+1
rho(i) = rho(11000)*exp(-(h(i)-11000)/exp_athmos);      %Set density values for h > 11000
end
density = interp1(h,rho,X(1));
end

[JXT119]

Thank you so much!! It worked perfectly. I see my main problem were the indexes in the complex_atm_model, thank you for the clarification!