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!

3 次查看（过去 30 天）
Rachel Ong
Rachel Ong 2022-3-3
评论： Rachel Ong 2022-3-3
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

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2022-3-3
编辑：Jan 2022-3-3
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 个评论
显示 1更早的评论
Jan
Jan 2022-3-3
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

请先登录，再进行评论。

更多回答（1 个）

Alan Stevens
Alan Stevens 2022-3-3
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by