import sympy as sm
import sympy.physics.mechanics as me
I recommend doing this problem by inspection of the graphic. It will take time but it will help you visualize the contributions to the vectors. Take $\mathbf{u}=\mathbf{a}_1 + 2\mathbf{a}_2+3\mathbf{a}_3$ for example. This vector is fixed in reference frame $A$, so the vector is not a function of any of the variables in $A$. If you then view the vector from $B$ it is apparent that changing $q_1$ can affect the direction of $u$ in $B$. Viewing from $C$ and $D$ will each add $q_2$ and $q_3$ as dependent variables, respectively.
For v it is important to note that $b_1$ and $c_2$ are always orthogonal, i.e. $q_2$ doesn't change them. Same with $d_3$ and $c_2$, orthogononal.
You can use SymPy to check your work though.
Create the three variables for the positive angles.
q1, q2, q3 = sm.symbols('q1:4')
Create the four reference frames by chaining simple rotations through each angle about the axis of rotation.
A = me.ReferenceFrame('A')
B = A.orientnew('B', 'Axis', (q1, A.x))
C = B.orientnew('C', 'Axis', (q2, B.y))
D = C.orientnew('D', 'Axis', (q3, C.z))
D.dcm(A)
D.dcm?
Signature: D.dcm(otherframe) Docstring: The direction cosine matrix between frames. This gives the DCM between this frame and the otherframe. The format is N.xyz = N.dcm(B) * B.xyz A SymPy Matrix is returned. Parameters ========== otherframe : ReferenceFrame The otherframe which the DCM is generated to. Examples ======== >>> from sympy.physics.vector import ReferenceFrame, Vector >>> from sympy import symbols >>> q1 = symbols('q1') >>> N = ReferenceFrame('N') >>> A = N.orientnew('A', 'Axis', [q1, N.x]) >>> N.dcm(A) Matrix([ [1, 0, 0], [0, cos(q1), -sin(q1)], [0, sin(q1), cos(q1)]]) File: ~/miniconda3/lib/python3.6/site-packages/sympy/physics/vector/frame.py Type: method
Create the three vectors.
u = A.x + 2*A.y + 3*A.z
v = B.x + C.y + D.z
w = D.x + 2*D.y + q3*D.z
u
v
w
You can determine if the variables are present in each vector by expressing each vector in the desired reference frame and inspecting the expression for the variables of interest. For example, for vector $\mathbf{w}$ in reference frame $A$.
w.express(A)
The free_symbols()
method of the vector can show all unique variables in the expression. You have to supply the reference frame.
w.free_symbols(A)
{q1, q2, q3}
Try $\mathbf{v}$ in $B$:
v.express(B)
v.free_symbols(B)
{q2}
Another way to check is to take the derivative of the vector with respect to the variable of interest in the reference frame of interest. If the derivative is zero then the variable is not present. The following nested loop prints N
or Y
if the derivative is zero or not. The derivative of the vector must be simplified to compare to zero. Note that simplify is a slow operation and is not guaranteed to simplify completely.
u.diff(q1, A)
v.diff(q1, A).simplify()
for frame in [A, B, C, D]:
for variable in [q1, q2, q3]:
for vector, vec_name in zip([u, v, w], ['u', 'v', 'w']):
if vector.diff(variable, frame).simplify() != sm.S(0):
ans = 'Y'
else:
ans = 'N'
print('{}, {}, {}: {}'.format(frame, variable, vec_name, ans))
A, q1, u: N A, q1, v: Y A, q1, w: Y A, q2, u: N A, q2, v: Y A, q2, w: Y A, q3, u: N A, q3, v: N A, q3, w: Y B, q1, u: Y B, q1, v: N B, q1, w: N B, q2, u: N B, q2, v: Y B, q2, w: Y B, q3, u: N B, q3, v: N B, q3, w: Y C, q1, u: Y C, q1, v: N C, q1, w: N C, q2, u: Y C, q2, v: Y C, q2, w: N C, q3, u: N C, q3, v: N C, q3, w: Y D, q1, u: Y D, q1, v: N D, q1, w: N D, q2, u: Y D, q2, v: Y D, q2, w: N D, q3, u: Y D, q3, v: Y D, q3, w: Y
Make $q_1$ a function of time and determine which reference frames $\mathbf{v}$ is a function of $t$ in.
q1 = me.dynamicsymbols('q1')
q2, q3 = sm.symbols('q2, q3')
A = me.ReferenceFrame('A')
B = A.orientnew('B', 'Axis', (q1, A.x))
C = B.orientnew('C', 'Axis', (q2, B.y))
D = C.orientnew('D', 'Axis', (q3, C.z))
You can access the default time variable that dynamicsymbols()
uses like so:
t = me.dynamicsymbols._t
v = B.x + C.y + D.z
v.diff(t, A).simplify()
v.diff(t, B).simplify()
v.diff(t, C).simplify()
v.diff(t, D).simplify()
Only in $A$ does $\mathbf{v}$ have a time derivative, thus $A$ is the answer.
This problem can be solved in the same way as 1.2 using SymPy. You can also think about whether each component of $\mathbf{v}$ will change if $q_2$ changes from each frame. From $A$ $\mathbf{d}_3$ will change. From $B$ $\mathbf{d}_3$ will change. From $C$ $\mathbf{b}_1$ will change. From $D$ $\mathbf{b}_1$ will change, thus it is a function of time in all frames.
q2 = me.dynamicsymbols('q2')
q1, q3 = sm.symbols('q1, q3')
A = me.ReferenceFrame('A')
B = A.orientnew('B', 'Axis', (q1, A.x))
C = B.orientnew('C', 'Axis', (q2, B.y))
D = C.orientnew('D', 'Axis', (q3, C.z))
t = me.dynamicsymbols._t
v = B.x + C.y + D.z
v.diff(t, A).simplify()
v.diff(t, B).simplify()
v.diff(t, C).simplify()
v.diff(t, D).simplify()
The derivatives have non-zero values in all frames.
A = me.ReferenceFrame('A')
Create some arbitrary scalars.
a1, a2, a3 = sm.symbols('alpha1:4')
b1, b2, b3 = sm.symbols('beta1:4')
g1, g2, g3 = sm.symbols('gamma1:4')
Create the three vectors using the arbitrary scalars.
alpha = a1 * A.x + a2 * A.y + a3 * A.z
alpha
beta = b1 * A.x + b2 * A.y + b3 * A.z
beta
gamma = g1 * A.x + g2 * A.y + g3 * A.z
gamma
The problem implies that the scalar triple product is commutative and associative. First check if crossing $\alpha$ with $\beta$ before dotting with gamma is equivalent to crossing beta with gamma and dotting with alpha. These can be shown to equivalent by doing the operations and expanding the resulting expressions.
me.dot(alpha, me.cross(beta, gamma)).expand()
me.dot(me.cross(alpha, beta), gamma).expand()
Similarly the gamma alpha beta and beta gamma alpha roder can also been shown to be equivalent.
me.dot(gamma, me.cross(alpha, beta)).expand()
me.dot(beta, me.cross(gamma, alpha)).expand()
We've already caculated these derivatives above, so we do it again but take the magnitudes of the vectors. Some need simplification to get the book result.
q1, q2, q3 = sm.symbols('q1:4')
A = me.ReferenceFrame('A')
B = A.orientnew('B', 'Axis', (q1, A.x))
C = B.orientnew('C', 'Axis', (q2, B.y))
D = C.orientnew('D', 'Axis', (q3, C.z))
v = B.x + C.y + D.z
v
v.express(A)
v.diff(q1, A).simplify().express(C)
v.diff(q1, A).magnitude().simplify()
v.magnitude()
type(v.magnitude())
sympy.core.power.Pow
type(q1 + q2)
sympy.core.add.Add
type(q1*q2)
sympy.core.mul.Mul
This should be the same as differentiating the measure numbers!
sm.sqrt(v.dot(A.x).diff(q1)**2 + v.dot(A.y).diff(q1)**2 + v.dot(A.z).diff(q1)**2).simplify()
v.diff(q1, B).magnitude()
v.diff(q2, C).magnitude().simplify()
v.diff(q3, C).magnitude()
v.diff(q2, D).magnitude().simplify()
v.diff(q1, D).magnitude()
q1, q2, q3 = sm.symbols('q1:4')
A = me.ReferenceFrame('A')
B = A.orientnew('B', 'Axis', (q1, A.x))
C = B.orientnew('C', 'Axis', (q2, B.y))
D = C.orientnew('D', 'Axis', (q3, C.z))
u = A.x + 2*A.y + 3*A.z
The measure numbers can be found by expressing a vector in the desired reference frame and inspecting the expression for the scalars of each of the three unit vectors.
u.diff(q1, B).express(A).simplify()
Each measure number can be extracted by using the dot product with each unit vector, i.e. finding the magnitude of the projection of the vector on that axis.
u.diff(q1, B).dot(A.x)
u.diff(q1, B).dot(A.y).simplify()
u.diff(q1, B).dot(A.z).simplify()
Similiarly for the other two derivatives.
u.diff(q1, B).express(B).simplify()
u.diff(q1, A).express(A).simplify()
There is convenience method for getting measure numbers in matrix form too:
u.diff(q1, B).to_matrix(A).simplify()
q1, q2, q3 = sm.symbols('q1:4')
A = me.ReferenceFrame('A')
B = A.orientnew('B', 'Axis', (q1, A.x))
C = B.orientnew('C', 'Axis', (q2, B.y))
D = C.orientnew('D', 'Axis', (q3, C.z))
alpha = q1 + q2 + q3
beta = q1**2 + q2**2 + q3**2
gamma = q1**3 + q2**3 + q3**3
r_PQ = alpha * A.x + beta * B.y + gamma * C.z
r_PQ
dr_dq3_D = r_PQ.diff(q3, D).simplify()
dr_dq3_D
dr_dq3_D.subs({q1: sm.pi/2, q2: 0, q3: 0}).magnitude()
A = me.ReferenceFrame('A')
v1, v2, v3 = me.dynamicsymbols('v1:4')
v = v1 * A.x + v2 * A.y + v3 * A.z
v.magnitude() * v.magnitude().diff(t)
me.dot(v, v.dt(A))
A unit vector aligned with $\mathbf{v}$ for all time can be found by dividing by the magnitude:
v / v.magnitude()
The normalize()
method of vectors does this operation:
u = v.normalize()
u
u.dt(A)
To show that the expression in the book is the same:
v.dt(A) / v.magnitude() - v.outer(v).dot(v.dt(A))/v.magnitude()**3
q1, q2 = sm.symbols('q1, q2')
f = sm.Function('f')(q1, q2)
f
g = sm.Function('g')(q1, q2)
g
h = sm.Function('h')
type(h)
h(q1, q2).diff(q1)
theta = sm.Function('theta')(q1, q2)
A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient(A, 'Axis', (theta, A.z))
v = f * A.x + g * A.y
v
Left side of the equation simplifies to the expected expression:
(v.diff(q2, B).diff(q1, A) - v.diff(q1, A).diff(q2, B)).simplify()
Make sure you can get these result when taking the derivatives (it is worth doing some by hand):
v.diff(q2, B).simplify()
v.diff(q2, B).diff(q1, A).simplify()
L = me.ReferenceFrame('L') # The laboratory frame
C = me.ReferenceFrame('C') # The disc frame
R = me.ReferenceFrame('R') # The rod frame
We need the q's to be functions of time so you can use dynamicsymbols()
to create these more easily that using Function()
explicitly.
q1, q2 = me.dynamicsymbols('q1, q2')
q1, q2
(q1(t), q2(t))
r = sm.symbols('r', real=True)
I'll assume that when $q_1=0$ the $L$ and $C$ unit vectors are equivalent, so a simple rotation about the 1 axis will orient $C$ with respect to $L$.
C.orient(L, 'Axis', (q1, L.x))
Similarly, I'll assume that when $q_2=0$ the rod points in the $-\mathbf{c}_1$ direction.
R.orient(C, 'Axis', (q2, C.y))
Verify the first two derivatives.
C.y.diff(q1, L).simplify()
C.z.diff(q1, L).simplify()
Form the vector from $O$ to $P$, making use of the two reference frames.
p = r*C.y - 3*r*R.x
p
p.express(C)
p.dt(L)
The following expression matches the result in the book:
p.dt(L).express(C)
Start by creating a vector that has measure numbers that are functions of time.
N = me.ReferenceFrame('N')
a1, a2, a3 = me.dynamicsymbols('a1:4')
v = a1*N.x + a2*N.y + a3*N.z
v
Now create a unit vector that is always aligned with $\mathbf{v}$.
u = v.normalize()
u
We are told that at $t^*$ vector $\mathbf{v}$ has specific properities. We can then infer that:
So twice differentiate $\mathbf{u}$ to get the general expression.
d2udt2 = u.dt(N).dt(N)
d2udt2
Now subsitute the values at time $t^*$ and compute the magnitude. Be careful that you substitue the highest derivatives first.
temp1 = d2udt2.subs({a1.diff(t).diff(t): 0, a2.diff(t).diff(t): 0, a3.diff(t).diff(t): 1})
temp2 = temp1.subs({a1.diff(t): 0, a2.diff(t): 1, a3.diff(t): 0})
temp2.subs({a1: 1, a2:0, a3: 0}).magnitude()