Integrate a ODE with anonymous functions

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:

Maltab code and functions

Answer 1

[Walter Roberson]

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.

[Giuseppe]

Oh thanks! A very simple mistakes with respect to what I expected.