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
