# JN8: Practice Activity 7#

This notebook use the example of a particle \(P\) sliding on a double pendulum (see figure above); this system was first presented to you in the file “7 Particle Kinematics”. Your objective is to compute the :

\(^N{\bf v}^P\), velocity of \(P\) in \(N\); and

\(^B{\bf v}^P\), velocity of \(P\) in \(A\).

You are to complete this task using:

the defintions of velocity \(^N{\bf v}^P \triangleq \frac{^N d}{dt}{\bf r}_{OP}\) and \(^A{\bf v}^P \triangleq \frac{^A d}{dt}{^N{\bf r}^{QP}}\)

the one-point theorem for velocity.

# Problem set-up#

## Create constants using `symbols`

and time-varying scalars using `dynamicsymbols`

#

Below we create the scalar constant \(l\), which is the length of the links \(A\) and \(B\); this is stored as the variable `l_scalar.`

:

```
from sympy import symbols, sin, cos
from sympy.physics.mechanics import dynamicsymbols, ReferenceFrame, init_vprinting
init_vprinting()
l_scalar = symbols('l')
```

A similar convention for naming variables is used to create `theta_1_scalar`

, `theta_2_scalar`

, and `x_scalar`

to represent \(\theta_1\), \(\theta_2\), and \(x\). These are created as `dynamicsymbols`

as they signifiy the time-varying variables in the figure above. This is done by:

```
theta_1_scalar, theta_2_scalar, x_scalar = dynamicsymbols('theta_1, theta_2, x')
x_scalar_dot = dynamicsymbols('x', 1)
```

## Create and orient reference frames \(N\), \(A\), and \(B\)#

Next, we create the three reference frames \(N\), \(A\), and \(B\), which are stored in the variable names `N`

, `A`

, and `B`

, respectively.

```
N = ReferenceFrame('N')
A = ReferenceFrame('A')
B = ReferenceFrame('B')
```

Since we are given the angles between the three reference frames, we can use the `.orient`

method to define the relative orientations using `sympy`

. We can define the between `N`

and `A`

(and also between `A`

and `B`

) using the `'Axis'`

rotation type within the `.orient`

method, as shown below:

```
A.orient(N, 'Axis', (theta_1_scalar, N.z))
A.ang_vel_in(N)
```

```
B.orient(A, 'Axis', (theta_2_scalar, A.z))
B.ang_vel_in(N)
```

## Create the points \(O\), \(Q\), and \(P\)#

Next, we need to create three points:

\(O\) is the point of contact between link \(A\) and ceiling \(N\); the corresponding variable name we use is

`O`

\(Q\) is the point of contact between link \(A\) and link \(B\); the corresponding variable name we use is

`Q`

\(P\) is also modeled as a point (massless particle) that slides on the link \(B\); the corresponding variable name is

`P`

.

Each of these points is made using `Point`

, which we import from `sympy.physics.mechanics`

, as below:

```
from sympy.physics.mechanics import Point
O = Point('O')
Q = Point('Q')
P = Point('P')
```

You can see that this command is quite similar in its syntax to how we create a reference frame using `sympy`

’s `ReferenceFrame`

feature. As you have seen in this and preceding activities, once we have created a `ReferenceFrame`

, we can exploit methods such as `orient`

, `set_ang_vel`

, and `set_ang_acc`

to describe the orientation kinematics. Similarly, we can use the `Point`

class to define the translational kinematics; in other words, we can exploit methods that are bundled with the `Point`

feature of `sympy`

to:

describe position vectors between points using

`.set_pos`

;define the velocities of points in different frames using

`.set_vel`

; andaccelerations of points in different frames using

`.set_acc`

.

In the remaining discussion of this chapter, we see how to use these methods in computing the kinematics of points by implementing computational versions of the two-point theorem abd one-point theorem. Let’s get started……

From, the figure, we know that \({\bf r}^{OQ}\) is the position vector from \(O\) to \(Q\). We can use `sympy`

to define the location of \(Q\) using the `.set_pos()`

method on that point. The `.set_pos()`

method requires that you (the user) provide two additional details to it:

The point from which \(Q\)’s location is defined relative; in this case, it is \(O\).

The component form of the vector in an appropriate reference frame. In this case, \(l \hat{\bf a}_y\) is the position vector \({\bf r}^{OQ}\), so we provide this as

`l_scalar*A.y`

to the`.set_pos()`

. This is shown below:

```
Q.set_pos(O, l_scalar*A.y)
```

As you can see, this method does not print anything in an output cell. However, the computer stores this information internally. To access this information, we have another method that we can use for any point called `.pos_from`

. This can be used to examine the position vector as shown below. Note that, I choose to actually store the output of `.pos_from`

in a variable name `r_OQ_vector`

as I know I will need it later for my computations of the two-point theorem.

```
r_OQ_vector = Q.pos_from(O)
r_OQ_vector
```

Similarlyl, we can configure the \(P\) by making use of the two methods described above, `set_pos`

and `pos_from`

, on P. Then, we can store the position vector information about P relative to Q (i.e., \({\bf r}^{QP}\)) in the variable name `r_QP_vector`

. This is shown below:

```
P.set_pos(Q, x_scalar*B.y)
r_QP_vector = P.pos_from(Q)
r_QP_vector
```

Now, from out lecture notes, we know that \({\bf r}^{OP}\),the position vector from \(O\) to \(P\), is given by:

So, we can use `sympy`

and our knowledge of vector addition to compute this position vector using the `r_OQ_vector`

and `r_QP_vector`

. I store this information in a variable `r_OP_vector`

, as shown below:

```
r_OP_vector = r_OQ_vector + r_QP_vector
r_OP_vector
```

Next, I describe how we can exploit `sympy`

and our knowledge on vector calculus to compute velocities of \(P\) from diffrerent reference frames.

# Approach 1: Using definitions#

## Compute \(^N{\bf v}^P\)#

The classical definition of velocity is well known (see figure below):

We use that and our knowledge of taking time derivatives in `sympy`

to compute \(^N{\bf v}^P\), which is stored in the variable `N_v_P_approach1`

.

```
N_v_P_approach1 = r_OP_vector.dt(N)
N_v_P_approach1.express(A)
```

## Compute \(^A{\bf v}^P\)#

\(^A{\bf v}^P\) is also computed in a similar fashion using classical definition and sympy; this information is stored in the variable name `A_v_P_approach1`

. See below:

```
A_v_P_approach1 = r_QP_vector.dt(A)
A_v_P_approach1.express(A)
```

# Approach 2: Implementing one point theorem#

In this section, we learn to compute the two velocities \(^N{\bf v}^P\) and \(^A{\bf v}^P\) using the one-point theorem.

## Compute \(^N{\bf v}^P\)#

To begin, we must first identify correctly how to utilize the one-point theorem. This has been explained in the video lecture so reference that to referesh your memory. In the image below, the resulting equation from correctly applying the one point theorem to compute \(^N{\bf v}^P\) is repeated from the handwritten notes :

The right hand side of the above equation tells us what we need to compute. Specifically, it tells us we will need to compute the following to determine \(^N{\bf v}^P\) in `sympy`

:

\(^B{\bf v}^P\), which we can easily evaluate from the figure;

\(^N{\bf v}^Q\), which we compute using the two-point theorem. This is discussed in greater detail below (and also covered in the lecture notes);

\(^N{\bf \omega}^B\); this information is actually already available from approach 1 above, where we used the

`.orient`

method. However, we haven’t stored this in a variable, which we will do in the steps that follow below; and\({\bf r}^{QP}\), which we have already computed above for Approach 1; this is already stored in the variable name

`r_QP_vector`

.

### Step 1: Define \(^B{\bf v}^P\) using `.set_vel`

#

From the figure, we can easily compute the \(^B{\bf v}^P = \dot x \hat {\bf b}_y\). We can then use the `.set_vel`

to set the velocity of `P`

in the `B`

-frame. This is then stored in the variable `B_v_P_vector`

using `.vel`

method.

_**The difference between .set_vel and .vel: the former lets us define the velocity and store it within the computer without an output; the latter method allows a human to access the velocity vector for future computation**. In other words, their behaviour is similar to what we see for creating position vectors using

`.set_pos`

and accessing it using `.pos_from`

.This is shown below:

```
P.set_vel(B, x_scalar_dot*B.y)
B_v_P_vector = P.vel(B)
B_v_P_vector
```

### Step 2: Compute \(^N{\bf v}^Q\)#

To compute this velocity, \(^N{\bf v}^Q\), we make use of the two point theorem. This is shown below from the lecture notes:

From the figure it is clear that point \(O\) is fixed in the two frames \(N\) and \(A\). So, for all our `sympy`

-bsed computations, we need to ensure that \(O\) is fixed in frames \(N\) and \(A\) as a point of zero velocity. This is done using our tried and tested `.set_vel`

method, as shown below:

```
O.set_vel(N, 0)
O.set_vel(A, 0)
```

Further, as we know that the two-point theorem for \(^N{\bf v}^Q\) requires \(^N{\bf v}^O\), we store this latter velocity in the variable named `N_v_O_vector`

using the `.vel`

method.

```
N_v_O_vector = O.vel(N)
```

In approach 1, we defined the orientaion of \(A\) in \(N\). So, we can easily use the `.ang_vel_in`

method to store the angular velocity of A in N (\(^N\omega^A\)) in the variable `N_w_A_vector`

:

```
N_w_A_vector = A.ang_vel_in(N)
```

Finally, we can compute the \(^N{\bf v}^Q\) by defining the the two-point formula, as below:

```
N_v_Q_vector = N_v_O_vector + N_w_A_vector.cross(r_OQ_vector)
N_v_Q_vector
```

### Step 3: Store \(^N\omega^B\) for use in one-point formula#

We can store \(^N\omega^B\) in the variable `N_w_B_vector`

using the `.ang_vel_in`

method on B, as shown below:

```
N_w_B_vector = B.ang_vel_in(N)
```

### Step 4: Finally, compute \(^N{\bf v}^P\)#

Now, we can finally compute the \(^N{\bf v}^P\) by implementing the one point theorem (repeating the hand-written equation so you can avoid scrolling up):

Indeed, we will store \(^N{\bf v}^P\) in a variable named `N_v_P_vector_approach2`

to be clear that we are using our second approach here. In `sympy`

, we implement the handwritten equation to compute \(^N{\bf v}^P\) by making use of vector addition (which you have already learned) as follows:

```
N_v_P_vector_approach2 = B_v_P_vector + N_v_Q_vector + N_w_B_vector.cross(r_QP_vector)
N_v_P_vector_approach2.express(A)
```

## Compute \(^A{\bf v}^P\)#

By one-point theorem, we know that \(^A{\bf v}^P\) can be evaluated using:

So, we can do this entire computation in `sympy`

through the few lines of code written below (This is very similar to how we computed \(^N{\bf v}^P\)). Note that:

we have already computed and stored \(^B{\bf v}^P\) in

`B_v_P_vector`

; andwe have already computed and stored \({\bf r}^{QP}\) in

`r_QP_vector`

.

We discuss the evaluation of the other terms in the below steps.

### Step 1: Set velocity of \(Q\) to zero in \(A\) and \(B\)#

Since Q is a point shared by and also fixed in the bodies \(A\) and \(B\), we use the `.set_vel`

method to set its velocity to be zero. We then store this in the `A_v_Q_vector`

for accessing in later computations:

```
Q.set_vel(A, 0)
Q.set_vel(B, 0)
A_v_Q_vector = Q.vel(A)
A_v_Q_vector
```

### Step 2: Store \(^A\omega^B\) to use in one-point theorem#

In approach 1, we already used the `orient`

method on \(B\) to define its orientaiton relative to \(A\). So, we can readily access its angular velocity using the `.ang_vel_in`

method on B. For future use, we store \(^A\omega^B\) in the variable `A_w_B_vector`

, as shown below:

```
A_w_B_vector = B.ang_vel_in(A)
```

### Step 3: Finally, compute \(^N{\bf v}^P\)#

Now, we can finally compute the \(^A{\bf v}^P\) by implementing the appropriate version of the one-point theorem in `sympy`

. We store \(^A{\bf v}^P\) in the variable `A_v_P_vector_approach2`

.

```
A_v_P_vector_approach2 = B_v_P_vector + A_v_Q_vector + A_w_B_vector.cross(r_QP_vector)
A_v_P_vector_approach2.express(A)
```

We can also compare that the resulting velocities are the same using either approach by comparing the `A_v_P_vector_approach2`

is identical to `A_v_P_vector_approach1`

as below:

```
A_v_P_vector_approach2.express(A) == A_v_P_approach1.express(A)
```

```
True
```

In English, we can say that the above input line is essentially us asking the computer: is the `A_v_P_vector_approach2`

the same as `A_v_P_approach1`

? Or, in other words, are the two approaches to computing velocities identical?

The `==`

sign is one type of a comparison operator in computer programming. For those interested in learning more, they are discussed here.