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!

7 visualizaciones (últimos 30 días)
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

Respuesta aceptada

Jan
Jan el 3 de Mzo. de 2022
Editada: Jan el 3 de Mzo. de 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 comentarios
Jan
Jan el 3 de Mzo. de 2022
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

Iniciar sesión para comentar.

Más respuestas (1)

Alan Stevens
Alan Stevens el 3 de Mzo. de 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

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by