Integrate a ODE with anonymous functions

11 views (last 30 days)
Giuseppe
Giuseppe on 10 May 2021
Commented: Giuseppe on 10 May 2021
Hi guys, I'm trying to solve a ODE problem by using the ode45 solver that receives in input a matrix that is composed by anonymous functions. I suppose that the problem is in the operations related to the functions inside the matrix.
The part of code of ODE solution is:
zero=zeros(3); %3x3 Matrix of zeros
zero34=zeros(3,4); %3x4 Matrix of zeros
zero43=zeros(4,3); %4x3 Matrix of zeros
zero44=zeros(4,4); %4x4 Matrix of zeros
% Right-hand-sides of the ODE system;
% - rows 1-6 are the equations of traslational and rotational motion of a rigid
%body;
% - rows 7-12 are the navigation and gimbal equations, to calculate CG coordinates
% and Euler angles in the Earth frame (ECEF)
dstate_dt= @(t,state) [-omegatilde(state),zero, zero,zero34;
zero, -(I)\(omegatilde(state)*I), zero,zero34;
T_DCM(state), zero, zero,zero34;
zero43, zero43, zero43,T_q2qdot(state)]*...
[u(state);v(state);w(state);
p(state) ; q(state); r(state);
xEG(state) ;yEG(state);zEG(state);
q_0(state);q_x(state);q_y(state); q_z(state)] +...
[g/W*[X(t,state);Y(t,state);Z(t,state)]; %% <------ Supposed problem HERE!
(I)\[L_roll(t,state);M_pitch(t,state);N_yaw(t,state)];zeros(7,1)];
%Ode45 options
options = odeset('RelTol',1e-9,'AbsTol',1e-9*ones(1,12));
%ODE solver
[Time_vect, State_vect] = ode45(dstate_dt, [0 t_fin], state0, options);
The code part relative to "suspect "functions is:
%Thrust
T = @(t) delta_t(t)*T_max;
%Lift
L = @(t,state) 0.5*rho(state)*(V(state)^2)*S*C_L(alpha(state),delta_e(t));
%Drag
D = @(t,state) 0.5*rho(state)*(V(state)^2)*S*C_D(alpha(state));
%Lateral Force
Y_a = @(t,state) 0.5*rho(state)*(V(state)^2)*S*C_Y(alpha(state),beta(state),delta_a(t),delta_r(t));
% Forces in body axes see Eqns. (7.21)
T_BE = @(state) quat2dcm([q_0(state) q_x(state) q_y(state) q_z(state)]);
F_G = @(state) T_BE(state) * [0; 0; W];
%Extract components of F_G by using anonymous function (ugly syntax!)
X_G = @(state) subsref(F_G(state), struct('type', '()', 'subs', {{1}})); %will return 1st element of F_G(state)
Y_G = @(state) subsref(F_G(state), struct('type', '()', 'subs', {{2}})); %will return 2nd element of F_G(state)
Z_G = @(state) subsref(F_G(state), struct('type', '()', 'subs', {{3}})); %will return 3rd element of F_G(state)
X = @(t,state) T(t)-D(t,state)*cos(alpha(state))+L(t,state)*sin(alpha(state)) + X_G(state);
Y = @(t,state) Y_a(t,state) + Y_G(state) ; %The yb and zb componens of thrust are zero
Z = @(t,state) -D(t,state)*sin(alpha(state))-L(t,state)*cos(alpha(state)) + Z_G(state);
% Aerodynamic Torques in body axes see Eqns. (7.42b)
L_roll = @(t,state) 0.5*rho(state)*V(state)^2*S*b*...
C_roll(alpha(state),beta(state),delta_a(t),delta_r(t),p(state),r(state));
M_pitch = @(t,state) 0.5*rho(state)*V(state)^2*S*c*...
C_pitch(alpha(state),delta_e(t),q(state));
N_yaw = @(t,state) 0.5*rho(state)*V(state)^2*S*b*...
C_yaw(alpha(state),beta(state),delta_a(t),delta_r(t),r(state));
I get the following error:
"Error using odearguments (line 144)
Solving @(T,STATE)[-OMEGATILDE(STATE),ZERO,ZERO,ZERO34;ZERO,-(I)\(OMEGATILDE(STATE)*I),ZERO,ZERO34;T_DCM(STATE),ZERO,ZERO,ZERO34;ZERO43,ZERO43,ZERO43,T_Q2QDOT(STATE)]*[U(STATE);V(STATE);W(STATE);P(STATE);Q(STATE);R(STATE);XEG(STATE);YEG(STATE);ZEG(STATE);Q_0(STATE);Q_X(STATE);Q_Y(STATE);Q_Z(STATE)]+[G/W*[X(T,STATE);Y(T,STATE);Z(T,STATE)];(I)\[L_ROLL(T,STATE);M_PITCH(T,STATE);N_YAW(T,STATE)];ZEROS(7,1)] requires a scalar AbsTol, or a vector AbsTol of length 13
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);"
The complete code and relative matlab functions here:

Accepted Answer

Walter Roberson
Walter Roberson on 10 May 2021
The message implies that your state0 is a vector of length 13, but your 'AbsTol',1e-9*ones(1,12) is only creating a vector of length 12.
  1 Comment
Giuseppe
Giuseppe on 10 May 2021
Oh thanks! A very simple mistakes with respect to what I expected.

Sign in to comment.

More Answers (0)

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!