% P1
% Direct Dynamics
% bar on the wall
clear all; clc; close;
L=2; m=1; g=10;
syms t
theta=sym('theta(t)');
omega = [0, 0, -diff(theta,t)];
alpha = diff(omega,t);
% diff(X,'t') or diff(X,sym('t'))
% differentiates a symbolic expression
% X with respect to t
% diff(X,'t',n) and diff(X,n,'t')
% differentiates X n times
% n is a positive integer
xC = (L/2)*cos(theta);
yC = (L/2)*sin(theta);
rC = [xC yC 0];
G = [0 -m*g 0];
IC = m*L^2/12;
vC = diff(rC,t);
aC = diff(vC,t);
fprintf('aC =')
pretty(aC)
fprintf('\n\n')

FAx = m*aC(1);
FBy = m*aC(2)+m*g;

fprintf('FAx =')
pretty(FAx)
fprintf('\n')
fprintf('FBy =')
pretty(FBy)
fprintf('\n\n')

rA = [0, L*sin(theta), 0];
rB = [L*cos(theta), 0, 0];
FA = [FAx, 0, 0];
FB = [0, FBy, 0];

MC = cross(rA, FA) + cross(rB, FB);

eq = -IC*alpha + MC;
eqI = simplify(eq(3));
fprintf('eom: ');
pretty(eqI)
fprintf('\n\n')

eqI1 = subs(eqI,diff('theta(t)',t,2),'ddtheta');
eqI2 = subs(eqI1,diff('theta(t)',t), sym('x(2)'));
eqI3 = subs(eqI2,'theta(t)', sym('x(1)'));

% solve a second order ODE using MATLAB
% ODE45 function
% write the second order equation as a
% system of two first order equations,
% by introducing x(2) = x(1)'
%  x(1)' = x(2)  ==>   x' = g(t,x)
%  x(2)' = f
%
% define the vector x = [x(1); x(2)]
% (column-vector).
% theta(t) =  x(1)
% d(theta)/dt = x(2)


% first differential equation d[x(1)]/dt = x(2)
dx1 = sym('x(2)');
% second differential equation d[x(2)]/dt =
dx2 = solve(eqI3,'ddtheta');

digits(3)
fprintf('d(x1)/dt = %s \n',char(vpa(dx1)))
fprintf('d(x2)/dt = %s \n',char(vpa(dx2)))

% define right-hand side vector
dx1dt = char(dx1);
dx2dt = char(dx2);

f = inline(sprintf('[%s; %s]',dx1dt,dx2dt),'t','x')

t0 = 0;  % define initial time
tf = 0.5; % define final time
time = [0 tf];
x0 = [pi/4 0]; % define initial conditions

[ts,xs] = ode45(f, time, x0);

x1 = xs(:,1); % extract x1 & x2 components from xs
x2 = xs(:,2);

subplot(2,1,1),plot(ts,x1,'r')
xlabel('t'),ylabel('\theta'), grid
subplot(2,1,2),plot(ts,x2)
xlabel('t'),ylabel('d \theta'),grid

% [ts,xs] = ode45(f, 0:0.5:10, x0);
% [ts,xs]

aC =
  +-                                 2                                                                                                           2   -+
  | - cos(theta(t)) diff(theta(t), t)  - sin(theta(t)) diff(theta(t), t, t), cos(theta(t)) diff(theta(t), t, t) - sin(theta(t)) diff(theta(t), t) , 0 |
  +-                                                                                                                                                 -+


FAx =
                                 2
- cos(theta(t)) diff(theta(t), t)  -

   sin(theta(t)) diff(theta(t), t, t)

FBy =
                                 2
- sin(theta(t)) diff(theta(t), t)  +

   cos(theta(t)) diff(theta(t), t, t) +

   10


eom: 
20 cos(theta(t)) +

   7 diff(theta(t), t, t)
   ----------------------
             3


d(x1)/dt = x(2) 
d(x2)/dt = (-8.57)*cos(x(1)) 

f =

     Inline function:
     f(t,x) = [x(2); -(60*cos(x(1)))/7]

Published with MATLAB® 7.9