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 visualizzazioni (ultimi 30 giorni)
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

Risposta accettata

Jan
Jan il 3 Mar 2022
Modificato: Jan il 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 Commenti
Jan
Jan il 3 Mar 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

Accedi per commentare.

Più risposte (1)

Alan Stevens
Alan Stevens il 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

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by