Exposing Noncontributing Forces

Note

You can download this example as a Python script: noncontributing.py or Jupyter Notebook: noncontributing.ipynb.

import numpy as np
import sympy as sm
import sympy.physics.mechanics as me
me.init_vprinting(use_latex='mathjax')

class ReferenceFrame(me.ReferenceFrame):

    def __init__(self, *args, **kwargs):

        kwargs.pop('latexs', None)

        lab = args[0].lower()
        tex = r'\hat{{{}}}_{}'

        super(ReferenceFrame, self).__init__(*args,
                                             latexs=(tex.format(lab, 'x'),
                                                     tex.format(lab, 'y'),
                                                     tex.format(lab, 'z')),
                                             **kwargs)
me.ReferenceFrame = ReferenceFrame

Learning Objectives

After completing this chapter readers will be able to:

• apply Newton’s Second Law to write equations of motion with maximal coordinates, which naturally expose noncontributing forces

• use the auxiliary generalized speed method to expose noncontributing forces in Kane’s minimal coordinate formulation

Introduction

Kane’s formulation relieves us from having to consider noncontributing forces (See Sec. Contributing and Noncontributing Forces), but often we are interested in one or more of these noncontributing forces. In this chapter, I will show how you can find the equation for a noncontributing force by introducing auxiliary generalized speeds. But first, let’s solve the equations of motion for a system by directly applying Newton’s Second Law of motion, which requires us to explicitly define all contributing and noncontributing forces.

Double Pendulum Example

Fig. 50 shows a schematic of a simple planar double pendulum described by two generalized coordinates \(q_1\) and \(q_2\). The particles \(P_1\) and \(P_2\) have masses \(m_1\) and \(m_2\), respectively. The lengths of the first and second pendulum arms are \(l_1\) and \(l_2\), respectively. On the right, the free body diagrams depict the two tension forces \(T_1\) and \(T_2\) that act on each particle to keep them at their respective radial locations.

Fig. 50 On the left, a kinematic diagram of a simple double planar pendulum with two links \(A\) of length \(l_1\) and \(B\) of length \(l_2\). On the right are free body diagrams of each particle showing all of the contributing and noncontributing forces acting on them. Gravity acts in the \(-\hat{n}_y\) direction.

Start by creating all of the necessary variables. The tension forces are time varying quantities.

m1, m2, l1, l2, g = sm.symbols('m1, m2, l1, l2, g')
q1, q2, u1, u2, T1, T2 = me.dynamicsymbols('q1, q2, u1, u2, T1, T2')
t = me.dynamicsymbols._t

p = sm.Matrix([m1, m2, l1, l2, g])
q = sm.Matrix([q1, q2])
u = sm.Matrix([u1, u2])
r = sm.Matrix([T1, T2])

ud = u.diff(t)

p, q, u, r, ud

\[\begin{split}\displaystyle \left( \left[\begin{matrix}m_{1}\\m_{2}\\l_{1}\\l_{2}\\g\end{matrix}\right], \ \left[\begin{matrix}q_{1}\\q_{2}\end{matrix}\right], \ \left[\begin{matrix}u_{1}\\u_{2}\end{matrix}\right], \ \left[\begin{matrix}T_{1}\\T_{2}\end{matrix}\right], \ \left[\begin{matrix}\dot{u}_{1}\\\dot{u}_{2}\end{matrix}\right]\right)\end{split}\]

Both pendulums’ configuration are described by angles relative to the vertical direction. We will choose the generalized speeds to be \(\bar{u} = \dot{\bar{q}}\) and set the angular velocities to be in terms of them.

N = me.ReferenceFrame('N')
A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')

A.orient_axis(N, q1, N.z)
B.orient_axis(N, q2, N.z)

A.set_ang_vel(N, u1*N.z)
B.set_ang_vel(N, u2*N.z)

Now the positions, velocities, and accelerations of each particle can be formed.

O = me.Point('O')
P1 = O.locatenew('P1', -l1*A.y)
P2 = P1.locatenew('P2', -l2*B.y)

O.set_vel(N, 0)
P1.v2pt_theory(O, N, A)

\[\displaystyle l_{1} u_{1}\hat{a}_x\]

P2.v2pt_theory(P1, N, B)

\[\displaystyle l_{1} u_{1}\hat{a}_x + l_{2} u_{2}\hat{b}_x\]

P1.a2pt_theory(O, N, A)

\[\displaystyle l_{1} \dot{u}_{1}\hat{a}_x + l_{1} u_{1}^{2}\hat{a}_y\]

P2.a2pt_theory(P1, N, B)

\[\displaystyle l_{1} \dot{u}_{1}\hat{a}_x + l_{1} u_{1}^{2}\hat{a}_y + l_{2} \dot{u}_{2}\hat{b}_x + l_{2} u_{2}^{2}\hat{b}_y\]

All of the kinematics are strictly in terms of the generalized coordinates and the generalized speeds.

Apply Newton’s Second Law Directly

Direction application of Newton’s Second Law can be done if all of the forces (noncontributing and contributing) are described for each of the two particles. Vector equations representing the law for each particle are:

(213)\[\begin{split}\sum\bar{F}^{P_1} = m_1 {}^N\bar{a}^{P_1} \\ \sum\bar{F}^{P_2} = m_2 {}^N\bar{a}^{P_2}\end{split}\]

From the free body diagram (Fig. 50) we see that all of the forces acting on \(P_1\) are:

F_P1 = T1*A.y - T2*B.y - m1*g*N.y
F_P1.express(N)

\[\displaystyle (- T_{1} \sin{\left(q_{1} \right)} + T_{2} \sin{\left(q_{2} \right)})\hat{n}_x + (- g m_{1} + T_{1} \cos{\left(q_{1} \right)} - T_{2} \cos{\left(q_{2} \right)})\hat{n}_y\]

and all of the forces acting on \(P_2\) are:

F_P2 = T2*B.y - m2*g*N.y
F_P2.express(N)

\[\displaystyle - T_{2} \sin{\left(q_{2} \right)}\hat{n}_x + (- g m_{2} + T_{2} \cos{\left(q_{2} \right)})\hat{n}_y\]

Now we can form the two vector expressions of Newton’s Second Law for each particle. Moving everything to the right hand side gives:

(214)\[\begin{split}\bar{0} = \sum\bar{F}^{P_1} - m_1 {}^N\bar{a}^{P_1} \\ \bar{0} = \sum\bar{F}^{P_2} - m_2 {}^N\bar{a}^{P_2}\end{split}\]

zero_P1 = F_P1 - m1*P1.acc(N)
zero_P2 = F_P2 - m2*P2.acc(N)

These two planar vector equations can then be written as four scalar equations by extracting the \(\hat{n}_x\) and \(\hat{n}_y\) measure numbers.

fd = sm.Matrix([
    zero_P1.dot(N.x),
    zero_P1.dot(N.y),
    zero_P2.dot(N.x),
    zero_P2.dot(N.y),
])
fd

\[\begin{split}\displaystyle \left[\begin{matrix}- l_{1} m_{1} \cos{\left(q_{1} \right)} \dot{u}_{1} - \left(- l_{1} m_{1} u_{1}^{2} + T_{1}\right) \sin{\left(q_{1} \right)} + T_{2} \sin{\left(q_{2} \right)}\\- g m_{1} - l_{1} m_{1} \sin{\left(q_{1} \right)} \dot{u}_{1} + \left(- l_{1} m_{1} u_{1}^{2} + T_{1}\right) \cos{\left(q_{1} \right)} - T_{2} \cos{\left(q_{2} \right)}\\l_{1} m_{2} u_{1}^{2} \sin{\left(q_{1} \right)} - l_{1} m_{2} \cos{\left(q_{1} \right)} \dot{u}_{1} - l_{2} m_{2} \cos{\left(q_{2} \right)} \dot{u}_{2} - \left(- l_{2} m_{2} u_{2}^{2} + T_{2}\right) \sin{\left(q_{2} \right)}\\- g m_{2} - l_{1} m_{2} u_{1}^{2} \cos{\left(q_{1} \right)} - l_{1} m_{2} \sin{\left(q_{1} \right)} \dot{u}_{1} - l_{2} m_{2} \sin{\left(q_{2} \right)} \dot{u}_{2} + \left(- l_{2} m_{2} u_{2}^{2} + T_{2}\right) \cos{\left(q_{2} \right)}\end{matrix}\right]\end{split}\]

It is important to note that these scalar equations are linear in both the time derivatives of the generalized speeds \(\dot{u}_1,\dot{u}_2\) as well as the two noncontributing force magnitudes \(T_1,T_2\) and that all four equations are coupled in these four variables.

(me.find_dynamicsymbols(fd[0]), me.find_dynamicsymbols(fd[1]),
 me.find_dynamicsymbols(fd[2]), me.find_dynamicsymbols(fd[3]))

\[\displaystyle \left( \left\{T_{1}, T_{2}, q_{1}, q_{2}, u_{1}, \dot{u}_{1}\right\}, \ \left\{T_{1}, T_{2}, q_{1}, q_{2}, u_{1}, \dot{u}_{1}\right\}, \ \left\{T_{2}, q_{1}, q_{2}, u_{1}, u_{2}, \dot{u}_{1}, \dot{u}_{2}\right\}, \ \left\{T_{2}, q_{1}, q_{2}, u_{1}, u_{2}, \dot{u}_{1}, \dot{u}_{2}\right\}\right)\]

That means we can write the equations as:

(215)\[\begin{split}\bar{f}_d(\dot{\bar{u}}, \bar{q}, \bar{r}, t) = \mathbf{M}_d \begin{bmatrix} \dot{\bar{u}} \\ \bar{r} \end{bmatrix} + \bar{g}_d\end{split}\]

where \(\bar{r} = \left[T_1 \ T_2 \right]^T\). The linear coefficient matrix and the remainder can be extracted as usual:

ud, r

\[\begin{split}\displaystyle \left( \left[\begin{matrix}\dot{u}_{1}\\\dot{u}_{2}\end{matrix}\right], \ \left[\begin{matrix}T_{1}\\T_{2}\end{matrix}\right]\right)\end{split}\]

udr = ud.col_join(r)
udr_zero = {v: 0 for v in udr}

Md = fd.jacobian(udr)
gd = fd.xreplace(udr_zero)

Md, udr, gd

\[\begin{split}\displaystyle \left( \left[\begin{matrix}- l_{1} m_{1} \cos{\left(q_{1} \right)} & 0 & - \sin{\left(q_{1} \right)} & \sin{\left(q_{2} \right)}\\- l_{1} m_{1} \sin{\left(q_{1} \right)} & 0 & \cos{\left(q_{1} \right)} & - \cos{\left(q_{2} \right)}\\- l_{1} m_{2} \cos{\left(q_{1} \right)} & - l_{2} m_{2} \cos{\left(q_{2} \right)} & 0 & - \sin{\left(q_{2} \right)}\\- l_{1} m_{2} \sin{\left(q_{1} \right)} & - l_{2} m_{2} \sin{\left(q_{2} \right)} & 0 & \cos{\left(q_{2} \right)}\end{matrix}\right], \ \left[\begin{matrix}\dot{u}_{1}\\\dot{u}_{2}\\T_{1}\\T_{2}\end{matrix}\right], \ \left[\begin{matrix}l_{1} m_{1} u_{1}^{2} \sin{\left(q_{1} \right)}\\- g m_{1} - l_{1} m_{1} u_{1}^{2} \cos{\left(q_{1} \right)}\\l_{1} m_{2} u_{1}^{2} \sin{\left(q_{1} \right)} + l_{2} m_{2} u_{2}^{2} \sin{\left(q_{2} \right)}\\- g m_{2} - l_{1} m_{2} u_{1}^{2} \cos{\left(q_{1} \right)} - l_{2} m_{2} u_{2}^{2} \cos{\left(q_{2} \right)}\end{matrix}\right]\right)\end{split}\]

The four equations are fully coupled, so we must solve for the four variables simultaneously. When applying Newton’s Second Law directly, additional coupled equations for each noncontributing force are necessary to solve the dynamical differential equations. When formulating the equations with Kane’s method, similar equations for the noncontributing forces can be generated, but the noncontributing forces will remain absent from the dynamical differential equations.

Auxiliary Generalized Speeds

When we form Kane’s equations, noncontributing forces will not be present in the equations of motion as they are above in the classical Newton formulation, but it is possible to expose select noncontributing forces by taking advantage of the role of the partial velocities. Forces and torques that are not normal to the partial velocity will contribute to the equations of motion. It is then possible to introduce fictitious partial velocities via an auxiliary generalized speed, along with a force or torque that acts in the same direction of the fictitious motion to generate extra equations for the noncontributing forces or torques. See [Kane1985] pg. 114 for more explanation of this idea.

As an example , here I introduce two fictitious generalized speeds, \(u_3\) and \(u_4\) that lets each particle have motion relative to its fixed location on the pendulum arm in the direction of the two noncontributing forces that we desire to know. Fig. 51 shows the two additional speeds and the associated forces. We introduce these speeds without introducing any related generalized coordinates.

Fig. 51 Kinematic diagram of the double pendulum showing the fictitious auxiliarly generalized speeds \(u_3\) and \(u_4\) and the associated contributing forces.

First find the velocity of \(P_1\) with the additional velocity component and store this separately in N_v_P1a to indicate it is affected by this auxiliary generalized speed.

u3, u4 = me.dynamicsymbols('u3, u4')

N_v_P1a = P1.vel(N) - u3*A.y
N_v_P1a

\[\displaystyle l_{1} u_{1}\hat{a}_x - u_{3}\hat{a}_y\]

Similarly, write the velocity of \(P_2\) using the velocity two point theorem and adding the auxiliary component. Note that the pendulum arm does not change in length because we have not added any generalized coordinates, so the two auxiliary velocities can be simply added in each step.

N_v_P2a = N_v_P1a + me.cross(B.ang_vel_in(N), P2.pos_from(P1)) - u4*B.y
N_v_P2a

\[\displaystyle l_{1} u_{1}\hat{a}_x - u_{3}\hat{a}_y + l_{2} u_{2}\hat{b}_x - u_{4}\hat{b}_y\]

These two velocities will be used to generate the partial velocities for two additional generalized active forces and generalized inertia forces, one for each of the auxiliary generalized speeds \(u_3\) and \(u_4\).

Auxiliary Generalized Active Forces

We now have four generalized speeds, two of which are auxiliary generalized speeds. With these speeds we will formulate four generalized active forces. The generalized active forces associated with \(u_1\) and \(u_2\) are no different than if we were not exposing the noncontributing forces, so we follow the usual procedure.

R_P1 = -m1*g*N.y
R_P2 = -m2*g*N.y

F1 = P1.vel(N).diff(u1, N).dot(R_P1) + P2.vel(N).diff(u1, N).dot(R_P2)
F1

\[\displaystyle - g l_{1} m_{1} \sin{\left(q_{1} \right)} - g l_{1} m_{2} \sin{\left(q_{1} \right)}\]

F2 = P1.vel(N).diff(u2, N).dot(R_P1) + P2.vel(N).diff(u2, N).dot(R_P2)
F2

\[\displaystyle - g l_{2} m_{2} \sin{\left(q_{2} \right)}\]

For \(F_3\) and \(F_4\), the contributing forces we wish to know that are associated with the auxiliary generalized speeds are added to the resultant acting on the two particles.

R_P1_aux = R_P1 + T1*A.y - T2*B.y
R_P2_aux = R_P2 + T2*B.y

Now the velocities of the particles that include the auxiliary generalized speeds are used to calculate the partial velocities and the auxiliary generalized active forces are formed.

F3 = N_v_P1a.diff(u3, N).dot(R_P1_aux) + N_v_P2a.diff(u3, N).dot(R_P2_aux)
F3

\[\displaystyle g m_{1} \cos{\left(q_{1} \right)} + g m_{2} \cos{\left(q_{1} \right)} - T_{1}\]

F4 = N_v_P1a.diff(u4, N).dot(R_P1_aux) + N_v_P2a.diff(u4, N).dot(R_P2_aux)
F4

\[\displaystyle g m_{2} \cos{\left(q_{2} \right)} - T_{2}\]

Finally, we form \(\bar{F}_r\) that consists of the two normal generalized active forces and the two auxiliary generalized active forces, the later two containing the unknown force magnitudes \(T_1\) and \(T_2\).

Fr = sm.Matrix([F1, F2, F3, F4])
Fr

\[\begin{split}\displaystyle \left[\begin{matrix}- g l_{1} m_{1} \sin{\left(q_{1} \right)} - g l_{1} m_{2} \sin{\left(q_{1} \right)}\\- g l_{2} m_{2} \sin{\left(q_{2} \right)}\\g m_{1} \cos{\left(q_{1} \right)} + g m_{2} \cos{\left(q_{1} \right)} - T_{1}\\g m_{2} \cos{\left(q_{2} \right)} - T_{2}\end{matrix}\right]\end{split}\]

Auxiliary Generalized Inertia Forces

Similar to the generalized active forces, the generalized inertia forces for \(u_1\) and \(u_2\) are computed as usual. See [Kane1985] pg. 169 and pg. 217 for more explanation.

Rs_P1 = -m1*P1.acc(N)
Rs_P2 = -m2*P2.acc(N)

F1s = P1.vel(N).diff(u1, N).dot(Rs_P1) + P2.vel(N).diff(u1, N).dot(Rs_P2)
F1s

\[\displaystyle - l_{1}^{2} m_{1} \dot{u}_{1} - l_{1}^{2} m_{2} \dot{u}_{1} + l_{1} \left(- l_{2} m_{2} \left(\sin{\left(q_{1} \right)} \sin{\left(q_{2} \right)} + \cos{\left(q_{1} \right)} \cos{\left(q_{2} \right)}\right) \dot{u}_{2} - l_{2} m_{2} \left(\sin{\left(q_{1} \right)} \cos{\left(q_{2} \right)} - \sin{\left(q_{2} \right)} \cos{\left(q_{1} \right)}\right) u_{2}^{2}\right)\]

F2s = P1.vel(N).diff(u2, N).dot(Rs_P1) + P2.vel(N).diff(u2, N).dot(Rs_P2)
F2s

\[\displaystyle - l_{2}^{2} m_{2} \dot{u}_{2} + l_{2} \left(- l_{1} m_{2} \left(\sin{\left(q_{1} \right)} \sin{\left(q_{2} \right)} + \cos{\left(q_{1} \right)} \cos{\left(q_{2} \right)}\right) \dot{u}_{1} - l_{1} m_{2} \left(- \sin{\left(q_{1} \right)} \cos{\left(q_{2} \right)} + \sin{\left(q_{2} \right)} \cos{\left(q_{1} \right)}\right) u_{1}^{2}\right)\]

The auxiliary generalized inertia forces are found using the velocities where \(u_3\) and \(u_4\) are present, but the acceleration of the particles need not include \(u_3\) and \(u_4\), because they are equal to zero because \(u_3\) and \(u_4\) are actually equal to zero.

F3s = N_v_P1a.diff(u3, N).dot(Rs_P1) + N_v_P2a.diff(u3, N).dot(Rs_P2)
F3s

\[\displaystyle l_{1} m_{1} u_{1}^{2} + l_{1} m_{2} u_{1}^{2} + l_{2} m_{2} \left(\sin{\left(q_{1} \right)} \sin{\left(q_{2} \right)} + \cos{\left(q_{1} \right)} \cos{\left(q_{2} \right)}\right) u_{2}^{2} + l_{2} m_{2} \left(- \sin{\left(q_{1} \right)} \cos{\left(q_{2} \right)} + \sin{\left(q_{2} \right)} \cos{\left(q_{1} \right)}\right) \dot{u}_{2}\]

F4s = N_v_P1a.diff(u4, N).dot(Rs_P1) + N_v_P2a.diff(u4, N).dot(Rs_P2)
F4s

\[\displaystyle l_{1} m_{2} \left(\sin{\left(q_{1} \right)} \sin{\left(q_{2} \right)} + \cos{\left(q_{1} \right)} \cos{\left(q_{2} \right)}\right) u_{1}^{2} + l_{1} m_{2} \left(\sin{\left(q_{1} \right)} \cos{\left(q_{2} \right)} - \sin{\left(q_{2} \right)} \cos{\left(q_{1} \right)}\right) \dot{u}_{1} + l_{2} m_{2} u_{2}^{2}\]

And finally, \(\bar{F}_r^*\) is formed for all four generalized speeds:

Frs = sm.Matrix([F1s, F2s, F3s, F4s])
Frs = sm.trigsimp(Frs)
Frs

\[\begin{split}\displaystyle \left[\begin{matrix}- l_{1} \left(l_{1} m_{1} \dot{u}_{1} + l_{1} m_{2} \dot{u}_{1} + l_{2} m_{2} u_{2}^{2} \sin{\left(q_{1} - q_{2} \right)} + l_{2} m_{2} \cos{\left(q_{1} - q_{2} \right)} \dot{u}_{2}\right)\\l_{2} m_{2} \left(l_{1} u_{1}^{2} \sin{\left(q_{1} - q_{2} \right)} - l_{1} \cos{\left(q_{1} - q_{2} \right)} \dot{u}_{1} - l_{2} \dot{u}_{2}\right)\\l_{1} m_{1} u_{1}^{2} + l_{1} m_{2} u_{1}^{2} + l_{2} m_{2} u_{2}^{2} \cos{\left(q_{1} - q_{2} \right)} - l_{2} m_{2} \sin{\left(q_{1} - q_{2} \right)} \dot{u}_{2}\\m_{2} \left(l_{1} u_{1}^{2} \cos{\left(q_{1} - q_{2} \right)} + l_{1} \sin{\left(q_{1} - q_{2} \right)} \dot{u}_{1} + l_{2} u_{2}^{2}\right)\end{matrix}\right]\end{split}\]

Warning

In this example, \(u_3,u_4,\dot{u}_3,\dot{u}_4\) are not present in the auxiliary generalized inertia forces but you may end up with auxiliary speeds and their derivatives in your auxiliary generalized inertia forces. If you do, you need to set them all to zero to arrive at the desired equations.

Augmented Dynamical Differential Equations

We can now form Kane’s dynamical differential equations which I will name \(\bar{f}_a\) to indicate they include the auxiliary equations. These equations are linear in \(\dot{u}_1,\dot{u}_2,T_1\) and \(T_2\).

fa = Frs + Fr
me.find_dynamicsymbols(fa)

\[\displaystyle \left\{T_{1}, T_{2}, q_{1}, q_{2}, u_{1}, u_{2}, \dot{u}_{1}, \dot{u}_{2}\right\}\]

Now when we extract the linear coefficients, we see that the dynamical differential equations (the first two rows) are independent of the unknown force magnitudes, allowing us to use the equations for \(\dot{\bar{u}}\) independently.

Ma = fa.jacobian(udr)
ga = fa.xreplace(udr_zero)

Ma, udr, ga

\[\begin{split}\displaystyle \left( \left[\begin{matrix}- l_{1} \left(l_{1} m_{1} + l_{1} m_{2}\right) & - l_{1} l_{2} m_{2} \cos{\left(q_{1} - q_{2} \right)} & 0 & 0\\- l_{1} l_{2} m_{2} \cos{\left(q_{1} - q_{2} \right)} & - l_{2}^{2} m_{2} & 0 & 0\\0 & - l_{2} m_{2} \sin{\left(q_{1} - q_{2} \right)} & -1 & 0\\l_{1} m_{2} \sin{\left(q_{1} - q_{2} \right)} & 0 & 0 & -1\end{matrix}\right], \ \left[\begin{matrix}\dot{u}_{1}\\\dot{u}_{2}\\T_{1}\\T_{2}\end{matrix}\right], \ \left[\begin{matrix}- g l_{1} m_{1} \sin{\left(q_{1} \right)} - g l_{1} m_{2} \sin{\left(q_{1} \right)} - l_{1} l_{2} m_{2} u_{2}^{2} \sin{\left(q_{1} - q_{2} \right)}\\- g l_{2} m_{2} \sin{\left(q_{2} \right)} + l_{1} l_{2} m_{2} u_{1}^{2} \sin{\left(q_{1} - q_{2} \right)}\\g m_{1} \cos{\left(q_{1} \right)} + g m_{2} \cos{\left(q_{1} \right)} + l_{1} m_{1} u_{1}^{2} + l_{1} m_{2} u_{1}^{2} + l_{2} m_{2} u_{2}^{2} \cos{\left(q_{1} - q_{2} \right)}\\g m_{2} \cos{\left(q_{2} \right)} + m_{2} \left(l_{1} u_{1}^{2} \cos{\left(q_{1} - q_{2} \right)} + l_{2} u_{2}^{2}\right)\end{matrix}\right]\right)\end{split}\]

We can solve the system to find functions for \(T_1\) and \(T_2\), if desired.

udr_sol = -Ma.LUsolve(ga)

T1_sol = sm.trigsimp(udr_sol[2])
T1_sol

\[\displaystyle g m_{1} \cos{\left(q_{1} \right)} + g m_{2} \cos{\left(q_{1} \right)} + l_{1} m_{1} u_{1}^{2} + l_{1} m_{2} u_{1}^{2} + l_{2} m_{2} u_{2}^{2} \cos{\left(q_{1} - q_{2} \right)} + \frac{l_{2} m_{2} \left(- g l_{2} m_{2} \sin{\left(q_{2} \right)} + l_{1} l_{2} m_{2} u_{1}^{2} \sin{\left(q_{1} - q_{2} \right)} + \frac{l_{1} l_{2} m_{2} \left(g m_{1} \sin{\left(q_{1} \right)} + g m_{2} \sin{\left(q_{1} \right)} + l_{2} m_{2} u_{2}^{2} \sin{\left(q_{1} - q_{2} \right)}\right) \cos{\left(q_{1} - q_{2} \right)}}{l_{1} m_{1} + l_{1} m_{2}}\right) \sin{\left(q_{1} - q_{2} \right)}}{\frac{l_{1} l_{2}^{2} m_{2}^{2} \cos^{2}{\left(q_{1} - q_{2} \right)}}{l_{1} m_{1} + l_{1} m_{2}} - l_{2}^{2} m_{2}}\]

T2_sol = sm.trigsimp(udr_sol[3])
T2_sol

\[\displaystyle g m_{2} \cos{\left(q_{2} \right)} + \frac{l_{1} l_{2} m_{2}^{2} \left(- g l_{2} m_{2} \sin{\left(q_{2} \right)} + l_{1} l_{2} m_{2} u_{1}^{2} \sin{\left(q_{1} - q_{2} \right)} + \frac{l_{1} l_{2} m_{2} \left(g m_{1} \sin{\left(q_{1} \right)} + g m_{2} \sin{\left(q_{1} \right)} + l_{2} m_{2} u_{2}^{2} \sin{\left(q_{1} - q_{2} \right)}\right) \cos{\left(q_{1} - q_{2} \right)}}{l_{1} m_{1} + l_{1} m_{2}}\right) \sin{\left(q_{1} - q_{2} \right)} \cos{\left(q_{1} - q_{2} \right)}}{\left(l_{1} m_{1} + l_{1} m_{2}\right) \left(\frac{l_{1} l_{2}^{2} m_{2}^{2} \cos^{2}{\left(q_{1} - q_{2} \right)}}{l_{1} m_{1} + l_{1} m_{2}} - l_{2}^{2} m_{2}\right)} - \frac{l_{1} m_{2} \left(g m_{1} \sin{\left(q_{1} \right)} + g m_{2} \sin{\left(q_{1} \right)} + l_{2} m_{2} u_{2}^{2} \sin{\left(q_{1} - q_{2} \right)}\right) \sin{\left(q_{1} - q_{2} \right)}}{l_{1} m_{1} + l_{1} m_{2}} + m_{2} \left(l_{1} u_{1}^{2} \cos{\left(q_{1} - q_{2} \right)} + l_{2} u_{2}^{2}\right)\]

Compare Newton and Kane Results

To ensure that the Newton approach and the Kane approach do produce equivalent results, we can numerically evaluate the equations with the same inputs and see if the results are the same. Here are some arbitrary numerical values for the states and constants.

q0 = np.array([
    np.deg2rad(15.0),  # q1 [rad]
    np.deg2rad(25.0),  # q2 [rad]
])

u0 = np.array([
    np.deg2rad(123.0),  # u1 [rad/s]
    np.deg2rad(-41.0),  # u2 [rad/s]
])

p_vals = np.array([
    1.2,  # m1 [kg]
    5.6,  # m2 [kg]
    1.34,  # l1 [m]
    6.7,  # l2 [m]
    9.81,  # g [m/2^2]
])

Create numeric functions to evaluate the two sets of matrices and execute both functions with the same numerical inputs from above.

eval_d = sm.lambdify((q, u, p), (Md, gd))
eval_a = sm.lambdify((q, u, p), (Ma, ga))

Md_vals, gd_vals = eval_d(q0, u0, p_vals)
Ma_vals, ga_vals = eval_a(q0, u0, p_vals)

Now compare the solutions for \(\begin{bmatrix}\dot{\bar{u}} & \bar{r} \end{bmatrix}\).

-np.linalg.solve(Md_vals, np.squeeze(gd_vals))

array([  8.09538007,  -2.37332094, 109.88598116,  92.50997719])

-np.linalg.solve(Ma_vals, np.squeeze(ga_vals))

array([  8.09538007,  -2.37332094, 109.88598116,  92.50997719])

For this set of inputs, the outputs are the same showing that using the auxiliary speed approach gives the same results, with the slight advantage that the dynamical differential equations are not coupled to the equations for the noncontributing forces in Kane’s method.

The forces can also be evaluated directly from the symbolic solutions, which is useful for post simulation application.

eval_forces = sm.lambdify((q, u, p), (T1_sol, T2_sol))
eval_forces(q0, u0, p_vals)

\[\displaystyle \left( 109.885981161619, \ 92.5099771909879\right)\]