How do I solve a symbolic system for a certain vector?

3 views (last 30 days)
Vítor Narciso
Vítor Narciso on 31 May 2022
Edited: Paul on 1 Jun 2022
Hello everyone!
I've run into the rigid-body dynamics equation below:
which, reduced for a constant ω (i.e., ), results in:
In theory, for any inertia tensor J, there is a certain angular velocity ω that spins the rigid-body along its major axis of inertia, thus the spin is stable. Therefore, in that case, and
I've been able to prove it in a simulation environment using Simulink, but I want to find the theoretical relation between ω and τ, i.e., an equation system that takes equation (2) and writes it in terms of:
I've tried to expand it manually but I believe it is easier to solve using MATLAB Symbolic Toolbox. So far, I've been ably to write it for the thrust τ, i.e., simply using equation (2) and the following code snippet:
syms(sym('w',[1 3],'real'));
w = [w1; w2; w3];
syms(sym('T',[1 3],'real'));
T = [T1; T2; T3];
syms(sym('J',[3 3],'real'));
J = [J1_1 J1_2 J1_3; J2_1 J2_2 J2_3; J3_1 J3_2 J3_3];
T = cross(w,J*w);
My question is, is there any way to redefine the result from T = cross(w,J*w); and rewrite the equation in termos of ?
Thank you for your time!
Best regards,
Vítor Narciso
Vítor Narciso on 1 Jun 2022
Yes, that's correct! In fact, I have solved by hand those equations for a diagonal inertia tensor (aka the Principal Moments of Inertia, as you stated). I wanted to use the full inertia tensor to solve it symbolically as, in my problem, I am dealing with Products of Inertia and it proves to be non-trivial to solve by hand.

Sign in to comment.

Answers (2)

Sam Chak
Sam Chak on 1 Jun 2022
Edited: Sam Chak on 1 Jun 2022
It's an interesting problem as mentioned by @Paul. The following does not solve your problem symbolically in MATLAB, but an attempt to investigate your rigid body problem in two scenarios by the numerical approach.
To begin, run the following code accordingly:
tspan = [0 10];
init = [1; 2; 3];
opt = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t, y] = ode45(@odefcn1, tspan, init, opt); % Scenario 1
% [t, y] = ode45(@odefcn2, tspan, init, opt); % Scenario 2
plot(t, y, 'linewidth', 1.5)
ylim([0 4]);
Scenario 1:
If the torque vector is designed as , then . Let's see what happens, but I guess you have already simulated this event in Simulink.
function dydt = odefcn1(t, y)
dydt = zeros(3,1);
J1 = 0.005; % principal moment of inertia about x-axis
J2 = 0.012; % principal moment of inertia about y-axis
J3 = 0.013; % principal moment of inertia about z-axis
dydt(1) = 0;
dydt(2) = 0;
dydt(3) = 0;
The angular velocities are stable (not diverging) but also not 'attractive' (not converging to the equilibrium point at ). It means that the rigid body has a stable rotation but tumbling. Can watch a video here:
However, if the initial angular velocity is large, then this is generally undesirable, because the rigid body is at the risk of spinning and breaking up structurally weak components in the real world due to the centrifugal force.
More importantly, it is obeserved that when is applied to rigid body, the rotational motions are successfully decoupled and they become independent.
Scenario 2:
Thus, if the torque vector is designed as with , then and the angular velocities are expected to exponentially decay to zero from their non-zero initial angular velocities.
function dydt = odefcn2(t, y)
dydt = zeros(3,1);
J1 = 0.005; % principal moment of inertia about x-axis
J2 = 0.012; % principal moment of inertia about y-axis
J3 = 0.013; % principal moment of inertia about z-axis
k = 1;
tau1 = - (J2 - J3)*y(2)*y(3) - J1*k*y(1); % torque about x-axis
tau2 = - (J3 - J1)*y(3)*y(1) - J2*k*y(2); % torque about y-axis
tau3 = - (J1 - J2)*y(1)*y(2) - J3*k*y(3); % torque about z-axis
dydt(1) = ((J2 - J3)*y(2)*y(3) + tau1)/J1;
dydt(2) = ((J3 - J1)*y(3)*y(1) + tau2)/J2;
dydt(3) = ((J1 - J2)*y(1)*y(2) + tau3)/J3;
As theoretically predicted above, the response of the rigid body is exponentially stable and stop rotating at steady state.
Non-trivial Solutions:
This problem requires rigorous mathematical thinking to solve, at least to me. Need pencil and paper. Solving the 3 equations simultaneously yields two sets solutions (due to having 3 sets of quadratics). The solutions are valid for three non-zero torque components () and three different principal moments of inertia ().
Solution Set #1
Solution Set #2
subject to
J1 = 0.005; % principal moment of inertia about x-axis
J2 = 0.012; % principal moment of inertia about y-axis
J3 = 0.013; % principal moment of inertia about z-axis
omega10 = 1; % initial omega1
omega20 = 2; % initial omega2
omega30 = 3; % initial omega3
tau1 = - (J2 - J3)*omega20*omega30 % torque about x-axis
tau2 = - (J3 - J1)*omega30*omega10 % torque about y-axis
tau3 = - (J1 - J2)*omega10*omega20 % torque about z-axis
Verifying the solution (set #1):
omega1 = (sqrt((J2 - J3)*tau2*tau3))/(sqrt((J1 - J2)*(J1 - J3)*tau1))
omega2 = (tau3*sqrt((J1 - J2)*(J1 - J3)*tau1))/((J2 - J1)*sqrt((J2 - J3)*tau2*tau3))
omega3 = (tau2*sqrt((J1 - J2)*(J1 - J3)*tau1))/((J1 - J3)*sqrt((J2 - J3)*tau2*tau3))
omega1 =
omega2 =
omega3 =

Paul on 1 Jun 2022
Edited: Paul on 1 Jun 2022
Hi Vitor,
To make sure I understand the question, it assumes that a rigid body can have a constant, non-zero angular velocity vector resolved in body-fixed coordinates while being subjected to a non-zero torque, also resolved in body-fixed coordinates.
Let's try to solve this problem for the simple case where the body-fixed frame is defined by the principal axes
syms w [3 1] real
syms T [3 1] real
syms J1_1 J2_2 J3_3 positive
eq = T == cross(w,diag([J1_1 J2_2 J3_3])*w);
sol = solve(eq,w,'ReturnConditions',true,'Real',true)
sol = struct with fields:
w1: [15×1 sym] w2: [15×1 sym] w3: [15×1 sym] parameters: [u x y z] conditions: [15×1 sym]
So Matlab found 15 solutions under various conditions. Let's look at the first three conditions
ans = 
So the first three conditions all have the three torque components equal to zero, i..e., torque free motion. But those solutions are not of interest and so can be discarded.
Let's look at the last solution
ans = 
[sol.w1(end) sol.w2(end) sol.w3(end)]
ans = 
So that might be a solution of interest for some constant z for the case where T3 = 0 and J1_1 = J2_2 and J2_2 ~= J3_3.
Presumably you can go through all of the sol.conditions and find the ones of interest.
If you try to expand this approach to the general case, you'll need to enforce all of constraints on the inertia tensor, like symmetric and positive definite to begin with. There may be others that I'm forgetting.
Having said that, it seems to me that you don't need to solve the general case explicitly. Because the w and T vectors are constant in the body-fixed, principal axes frame, they are constant in any body-fixed frame. So start with whatever body-fixed frame you have, derive the direction cosine matrix from that frame to the principal axes, tranform the intertia tensor to principal axes, solve the problem in the principal axes frame, then tansform the solution back to the original body-fixed frame.
Interesting problem.

Community Treasure Hunt

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

Start Hunting!

Translated by