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

Question

Rachel Ong 2022-1-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1627225-hi-is-it-possible-to-have-variables-in-an-ode45-that-varies-according-to-another-function

评论： Stephen23 2022-1-12

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