ODE45 Output size

I am trying to subtract w_ode from w. I used two methods to get w and I want to see the difference between them. But, the ode45 returns a w_ode(1:3.:) of size 3x645 and the other w has a size 3x101. Shouldn't ode45 return a 3x101 because the t = 0:100. So, the inputed time is a matrix of size 1x101.

I_T = 2;

J = 1.5;

I = [I_T; I_T; J];

w = [-0.1; 0.3; -1.1];

L = [0;0;0];

t = 0:100;

% Using ode45 to integrate w dot

options = odeset('RelTol',1e-10,'AbsTol',1e-10);

result = ode45(@dwdt_Bframe, t, [w; I; L], options);

% Extracting information from ode solver

t = result.x;

w_ode = result.y;

A1 = w(1);

A2 = w(2);

B1 = w(2);

B2 = -w(1);

K = (I_T - J)/I_T;

w3 = w(3);

t = 0:100;

w1 = A1*cos(K*w3*t) + B1*sin(K*w3*t);

w2 = A2*cos(K*w3*t) + B2*sin(K*w3*t);

w3 = w3;

w = [w1; w2; w3];

function dwdt = dwdt_Bframe(t, wIL)

w = wIL(1:3);

I = wIL(4:6);

L = wIL(7:9);

% The w dot equations

dwdt = zeros(9,1);

dwdt(1) = (-(I(3) - I(2))*w(2)*w(3) + L(1)) / I(1);

dwdt(2) = (-(I(1) - I(3))*w(3)*w(1) + L(2)) / I(2);

dwdt(3) = (-(I(2) - I(1))*w(1)*w(2) + L(3)) / I(3);

% Keep the values of I and L constant

dwdt(4:6) = 0;

dwdt(7:9) = 0;

end

Answer 1

[Sam Chak]

Hi @Sneh

I think you want to compare the numerical solutions with the analytical solutions.

%% Method 1: Numerical Solutions
% Settings for ode45 solver
tspan  = 0:100;                 % integration time
w0     = [-0.1; 0.3; -1.1];     % initial omega values
opts   = odeset('RelTol', 1e-10, 'AbsTol', 1e-10);

% Using ode45 to integrate w dot
[t, wn] = ode45(@dwdt_Bframe, tspan, w0, opts);  % numerical solutions are stored in wn

% Plot the numerical solutions, {ω₁, ω₂, ω₃}
subplot(211)
plot(t, wn), grid on, title('Numerical Solutions (via ode45)')
ylim([-1.5 0.5])

%% Method 2: Analytical Solutions (by formulas)
I_T = 2;
J   = 1.5;

A1  =  w0(1);   % ω₁'s magnitude of cos component
A2  =  w0(2);   % ω₂'s magnitude of cos component
B1  =  w0(2);   % ω₁'s magnitude of sin component
B2  = -w0(1);   % ω₂'s magnitude of sin component
K   = (I_T - J)/I_T;
w30 =  w0(3);   % initial value of ω₃

% the formulas
w1  = A1*cos(K*w30*t) + B1*sin(K*w30*t);
w2  = A2*cos(K*w30*t) + B2*sin(K*w30*t);
w3  = w30*ones(numel(t), 1);
wa  = [w1, w2, w3];

% Plot the analytical solutions, {ω₁, ω₂, ω₃}
subplot(212)
plot(t, wa), grid on, title('Analytical Solutions (by formulas)')
ylim([-1.5 0.5]), xlabel('Time / seconds')

% Rotational dynamics
function dwdt = dwdt_Bframe(t, w)
    I_T     = 2;
    J       = 1.5;
    I       = [I_T; I_T; J];
    L       = [0; 0; 0];
    
    % The w dot equations
    dwdt    = zeros(3, 1);
    dwdt(1) = (- (I(3) - I(2))*w(2)*w(3) + L(1)) / I(1);
    dwdt(2) = (- (I(1) - I(3))*w(3)*w(1) + L(2)) / I(2);
    dwdt(3) = (- (I(2) - I(1))*w(1)*w(2) + L(3)) / I(3);
end