Hello, I am trying to solve the ODE but I have an error I do not understand how to solve. Would appreciate some help. Thank you!

I have this set of ode and I would like to solve it with ode45. I have this one variable deltae and this depends on the solution on x(6) but I am not sure if this is how I should be calling it. How can I save the values of deltae at each step of ode45 to plot at the end later?

The error comes from the [t,x] = ode45 line. I have pasted the error as shown below.

Error using odearguments (line 95)
@(T,X)ODEFCN(T,X) returns a vector of length 1, but the length of initial conditions
vector is 7. The vector returned by @(T,X)ODEFCN(T,X) and the initial conditions
vector must have the same number of elements.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in ODE_max_speed (line 25)
[t,x] = ode45(@(t,x) odefcn(t,x), t_span, init_cond)

My script is attached as below:

%% ODE45_estimation

% ------ Constants ------ %
m = 1248.5;
g = 9.81;
W = m*g;
S = 17.1;

% ------ Time interval for simulation ------ %
t_span = [0,20];         % [start time, end time]

% ------ Initial conditions ------ %
alpha0 = 0.0584957579231716;      % Using trim conditions at sealevel for PS = 2.5
u0 = 55.2682238485092;
q0 = 0;
theta0 = 0.0584957579231716;
mass0 = 1248.5;
h0 = 0;
displacement0 = 0;

init_cond = [alpha0;u0;q0;theta0;mass0;h0;displacement0];

% ------ Solution ------ %
% eledefl_table = zeros(length(t),1)
[t,x] = ode45(@(t,x) odefcn(t,x), t_span, init_cond)


%% ODE equations function
function dxdt = odefcn(t,x)
    g = 9.81;
    b = 10.2;
    c = 1.74;
    S = 17.1;
    initial_mass = 1248.5;
    Ix = 1421;
    Iy = 4067.5;
    Iz = 4786;
    Ixz = 200;
    csfc = 200/60/60/1000;
    
    PS = 2.5;
    
    alpha = x(1);
    u = x(2); 
    q = x(3);
    theta = x(4);
    mass = x(5);
    height = x(6);
    displacement = x(7);
    
    [rho,speedsound] = atmos(x(6));
    
    deltae = 1.387e-15*height^3 - 1.314e-11*height^2 - 9.405e-07*height - 0.02487    % INSERT THE EQUATION
    
    epsilon = 0;    % CHANGE EPSILON LATER 
    
    qbar = 0.5*rho*u^2*S;
    CL = 6.44*alpha + 3.8*c*q/(2*u) + 0.355*deltae;
    L = qbar*CL;
    CD = 0.03 + 0.05*CL^2;
    D = qbar*CD;
    Cm = 0.05 - 0.683*alpha - 9.96*c*q/(2*u) - 0.923*deltae;
    m = qbar*c*Cm;
    
    thrust = PS * ((7+u/speedsound)*200/3 + height/1000*(2*u/speedsound - 11) );
    
    % dxdt(1) = alpha, dxdt(2) = u, dxdt(3) = q, dxdt(4) = theta, 
    % dxdt(5) = mass, dxdt(6) = height, dxdt(7) = displacement 
    
    dxdt = @(x) [q - (thrust*sin(alpha+epsilon) + L)/(mass*u) + g/u*cos(theta-alpha);
                 (thrust*cos(alpha-epsilon) - D)/mass - g*sin(theta-alpha);
                 m/Iy;
                 q;
                 -csfc*thrust;
                 u*sin(theta-alpha);
                 u*cos(theta-alpha)];
             
end

%% 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

Answer 1

[Jan]

Omit the "@(x)" in:

    dxdt = @(x) [q - (thrust*sin(alpha+epsilon) + L)/(mass*u) + g/u*cos(theta-alpha);
    ... %  ^^^^
                 (thrust*cos(alpha-epsilon) - D)/mass - g*sin(theta-alpha);
                 m/Iy;
                 q;
                 -csfc*thrust;
                 u*sin(theta-alpha);
                 u*cos(theta-alpha)]; 

You want to reply a vector, not a scalar function handle, which produces a vector.

Some hints:

• x^0.5 is much slower than sqrt(x).
• The exp() function is very expensive. Avoid to call it twice with the same argument:

    p = p1 * exp(-g/(R*T)*(h-h1));
    rho = rho1 * exp(-g/(R*T)*(h-h1));
    % Faster:
    tmp = exp(-g/(R*T)*(h-h1));
    p   = p1 * tmp;
    rho = rho1 * tmp;
    

[Rachel Ong]

Ah thank you! I didnt see that.

Do you perhaps know how I can possibly store my values of deltae somewhere so that I can plot it against each time step later?

Thanks for your help on the error once again.

[Jan]

It is not useful to store data from inside the function to be integrated. Remember that the integrator can reject steps, if they exceed the tolerances. In addition the function is called multiple times for one step of the integrator.

Better so not store the values, but let ODE45 run at first. Afterwards use the output of the integeration as input of the function to be integrated or to calculate the wanted value directly:

[t, x] = ode45(...);
height = x(:, 6);
deltae = 1.387e-15 * height.^3 - 1.314e-11 * height.^2 - 9.405e-07.*height - 0.02487

[Rachel Ong]

I see, I will try that as suggested. Thank you!

Answer 2

[Alan Stevens]

Like this?

%% ODE45_estimation

% ------ Constants ------ %
m = 1248.5;
g = 9.81;
W = m*g;
S = 17.1;

% ------ Time interval for simulation ------ %
t_span = [0,20];         % [start time, end time]

% ------ Initial conditions ------ %
alpha0 = 0.0584957579231716;      % Using trim conditions at sealevel for PS = 2.5
u0 = 55.2682238485092;
q0 = 0;
theta0 = 0.0584957579231716;
mass0 = 1248.5;
h0 = 0;
displacement0 = 0;

init_cond = [alpha0;u0;q0;theta0;mass0;h0;displacement0];

% ------ Solution ------ %
% eledefl_table = zeros(length(t),1)
[t,x] = ode45(@(t,x) odefcn(t,x), t_span, init_cond);

plot(t,x(:,7)),grid
xlabel('t'), ylabel('displacement?')

%% ODE equations function
function dxdt = odefcn(~,x)
    g = 9.81;
   % b = 10.2;
    c = 1.74;
    S = 17.1;
    %initial_mass = 1248.5;
    %Ix = 1421;
    Iy = 4067.5;
    %Iz = 4786;
    %Ixz = 200;
    csfc = 200/60/60/1000;
    
    PS = 2.5;
    
    alpha = x(1);
    u = x(2); 
    q = x(3);
    theta = x(4);
    mass = x(5);
    height = x(6);
   % displacement = x(7);
    
    [rho,speedsound] = atmos(x(6));
    
    deltae = 1.387e-15*height^3 - 1.314e-11*height^2 - 9.405e-07*height - 0.02487;    % INSERT THE EQUATION
    
    epsilon = 0;    % CHANGE EPSILON LATER 
    
    qbar = 0.5*rho*u^2*S;
    CL = 6.44*alpha + 3.8*c*q/(2*u) + 0.355*deltae;
    L = qbar*CL;
    CD = 0.03 + 0.05*CL^2;
    D = qbar*CD;
    Cm = 0.05 - 0.683*alpha - 9.96*c*q/(2*u) - 0.923*deltae;
    m = qbar*c*Cm;
    
    thrust = PS * ((7+u/speedsound)*200/3 + height/1000*(2*u/speedsound - 11) );
    
    % dxdt(1) = alpha, dxdt(2) = u, dxdt(3) = q, dxdt(4) = theta, 
    % dxdt(5) = mass, dxdt(6) = height, dxdt(7) = displacement 
    
    dxdt =  [q - (thrust*sin(alpha+epsilon) + L)/(mass*u) + g/u*cos(theta-alpha);
                 (thrust*cos(alpha-epsilon) - D)/mass - g*sin(theta-alpha);
                 m/Iy;
                 q;
                 -csfc*thrust;
                 u*sin(theta-alpha);
                 u*cos(theta-alpha)];
             
end

%% 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