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!

1 view (last 30 days)
Rachel Ong
Rachel Ong on 3 Mar 2022
Commented: Rachel Ong on 3 Mar 2022
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

Accepted Answer

Jan
Jan on 3 Mar 2022
Edited: Jan on 3 Mar 2022
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;
  3 Comments

Sign in to comment.

More Answers (1)

Alan Stevens
Alan Stevens on 3 Mar 2022
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

Tags

Community Treasure Hunt

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

Start Hunting!