Beginner Tutorial for SymPy Mechanics

Once you are at the python command line the first step is to import basic functionality from SymPy and the Mechanics module, otherwise you will only have basic python commands available to work with. We will import the symbols function from SymPy core and with the * method bring in all functionality from the mechanics package.

In [1]:
from sympy import symbols
from sympy.physics.mechanics import *

You can now see what functions and variables that are available to you with::

In [2]:
dir()
Out[2]:
['Dyadic',
 'In',
 'KanesMethod',
 'LagrangesMethod',
 'Lagrangian',
 'MechanicsLatexPrinter',
 'MechanicsPrettyPrinter',
 'MechanicsStrPrinter',
 'Out',
 'Particle',
 'Point',
 'ReferenceFrame',
 'RigidBody',
 'Vector',
 '_',
 '__',
 '___',
 '__builtin__',
 '__builtins__',
 '__name__',
 '__package__',
 '_dh',
 '_i',
 '_i1',
 '_i2',
 '_ih',
 '_ii',
 '_iii',
 '_oh',
 '_sh',
 'angular_momentum',
 'cross',
 'dot',
 'dynamicsymbols',
 'exit',
 'express',
 'get_ipython',
 'help',
 'inertia',
 'inertia_of_point_mass',
 'kinematic_equations',
 'kinetic_energy',
 'linear_momentum',
 'mechanics_printing',
 'mlatex',
 'mpprint',
 'mprint',
 'outer',
 'partial_velocity',
 'potential_energy',
 'quit',
 'symbols']

This is a long list of available functions. Read about the python import statement to learn about better ways to import only what you need. One good explanation is http://effbot.org/zone/import-confusion.htm.

To get started working with vectors we will need to create a reference frame, as all vectors need to be defined with respect to a reference frame in the mechanics package. If you know the name of the command that you want to use simply use the builtin help function to bring up the documentation for the function. In our case we need to use the ReferenceFrame class.

In [3]:
help(ReferenceFrame)
Help on class ReferenceFrame in module sympy.physics.mechanics.essential:

class ReferenceFrame(__builtin__.object)
 |  A reference frame in classical mechanics.
 |  
 |  ReferenceFrame is a class used to represent a reference frame in classical
 |  mechanics. It has a standard basis of three unit vectors in the frame's
 |  x, y, and z directions.
 |  
 |  It also can have a rotation relative to a parent frame; this rotation is
 |  defined by a direction cosine matrix relating this frame's basis vectors to
 |  the parent frame's basis vectors.  It can also have an angular velocity
 |  vector, defined in another frame.
 |  
 |  Methods defined here:
 |  
 |  __getitem__(self, ind)
 |      Returns basis vector for the provided index (index being an str)
 |  
 |  __init__(self, name, indices=None, latexs=None)
 |      ReferenceFrame initialization method.
 |      
 |      A ReferenceFrame has a set of orthonormal basis vectors, along with
 |      orientations relative to other ReferenceFrames and angular velocities
 |      relative to other ReferenceFrames.
 |      
 |      Parameters
 |      ==========
 |      
 |      indices : list (of strings)
 |          If custom indices are desired for console, pretty, and LaTeX
 |          printing, supply three as a list. The basis vectors can then be
 |          accessed with the get_item method.
 |      latexs : list (of strings)
 |          If custom names are desired for LaTeX printing of each basis
 |          vector, supply the names here in a list.
 |      
 |      Examples
 |      ========
 |      
 |      >>> from sympy.physics.mechanics import ReferenceFrame, mlatex
 |      >>> N = ReferenceFrame('N')
 |      >>> N.x
 |      N.x
 |      >>> O = ReferenceFrame('O', ('1', '2', '3'))
 |      >>> O.x
 |      O['1']
 |      >>> O['1']
 |      O['1']
 |      >>> P = ReferenceFrame('P', latexs=('A1', 'A2', 'A3'))
 |      >>> mlatex(P.x)
 |      'A1'
 |  
 |  __iter__(self)
 |  
 |  __repr__ = __str__(self)
 |  
 |  __str__(self)
 |      Returns the name of the frame.
 |  
 |  ang_acc_in(self, otherframe)
 |      Returns the angular acceleration Vector of the ReferenceFrame.
 |      
 |      Effectively returns the Vector:
 |      ^N alpha ^B
 |      which represent the angular acceleration of B in N, where B is self, and
 |      N is otherframe.
 |      
 |      Parameters
 |      ==========
 |      
 |      otherframe : ReferenceFrame
 |          The ReferenceFrame which the angular acceleration is returned in.
 |      
 |      Examples
 |      ========
 |      
 |      >>> from sympy.physics.mechanics import ReferenceFrame, Vector
 |      >>> N = ReferenceFrame('N')
 |      >>> A = ReferenceFrame('A')
 |      >>> V = 10 * N.x
 |      >>> A.set_ang_acc(N, V)
 |      >>> A.ang_acc_in(N)
 |      10*N.x
 |  
 |  ang_vel_in(self, otherframe)
 |      Returns the angular velocity Vector of the ReferenceFrame.
 |      
 |      Effectively returns the Vector:
 |      ^N omega ^B
 |      which represent the angular velocity of B in N, where B is self, and
 |      N is otherframe.
 |      
 |      Parameters
 |      ==========
 |      
 |      otherframe : ReferenceFrame
 |          The ReferenceFrame which the angular velocity is returned in.
 |      
 |      Examples
 |      ========
 |      
 |      >>> from sympy.physics.mechanics import ReferenceFrame, Vector
 |      >>> N = ReferenceFrame('N')
 |      >>> A = ReferenceFrame('A')
 |      >>> V = 10 * N.x
 |      >>> A.set_ang_vel(N, V)
 |      >>> A.ang_vel_in(N)
 |      10*N.x
 |  
 |  dcm(self, otherframe)
 |      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.mechanics import ReferenceFrame, Vector
 |      >>> from sympy import symbols
 |      >>> q1 = symbols('q1')
 |      >>> N = ReferenceFrame('N')
 |      >>> A = N.orientnew('A', 'Axis', [q1, N.x])
 |      >>> N.dcm(A)
 |      [1,       0,        0]
 |      [0, cos(q1), -sin(q1)]
 |      [0, sin(q1),  cos(q1)]
 |  
 |  orient(self, parent, rot_type, amounts, rot_order='')
 |      Defines the orientation of this frame relative to a parent frame.
 |      
 |      Supported orientation types are Body, Space, Quaternion, Axis.
 |      Examples show correct usage.
 |      
 |      Parameters
 |      ==========
 |      
 |      parent : ReferenceFrame
 |          The frame that this ReferenceFrame will have its orientation matrix
 |          defined in relation to.
 |      rot_type : str
 |          The type of orientation matrix that is being created.
 |      amounts : list OR value
 |          The quantities that the orientation matrix will be defined by.
 |      rot_order : str
 |          If applicable, the order of a series of rotations.
 |      
 |      Examples
 |      ========
 |      
 |      >>> from sympy.physics.mechanics import ReferenceFrame, Vector
 |      >>> from sympy import symbols
 |      >>> q0, q1, q2, q3, q4 = symbols('q0 q1 q2 q3 q4')
 |      >>> N = ReferenceFrame('N')
 |      >>> B = ReferenceFrame('B')
 |      
 |      Now we have a choice of how to implement the orientation. First is
 |      Body. Body orientation takes this reference frame through three
 |      successive simple rotations. Acceptable rotation orders are of length
 |      3, expressed in XYZ or 123, and cannot have a rotation about about an
 |      axis twice in a row.
 |      
 |      >>> B.orient(N, 'Body', [q1, q2, q3], '123')
 |      >>> B.orient(N, 'Body', [q1, q2, 0], 'ZXZ')
 |      >>> B.orient(N, 'Body', [0, 0, 0], 'XYX')
 |      
 |      Next is Space. Space is like Body, but the rotations are applied in the
 |      opposite order.
 |      
 |      >>> B.orient(N, 'Space', [q1, q2, q3], '312')
 |      
 |      Next is Quaternion. This orients the new ReferenceFrame with
 |      Quaternions, defined as a finite rotation about lambda, a unit vector,
 |      by some amount theta.
 |      This orientation is described by four parameters:
 |      q0 = cos(theta/2)
 |      q1 = lambda_x sin(theta/2)
 |      q2 = lambda_y sin(theta/2)
 |      q3 = lambda_z sin(theta/2)
 |      Quaternion does not take in a rotation order.
 |      
 |      >>> B.orient(N, 'Quaternion', [q0, q1, q2, q3])
 |      
 |      Last is Axis. This is a rotation about an arbitrary, non-time-varying
 |      axis by some angle. The axis is supplied as a Vector. This is how
 |      simple rotations are defined.
 |      
 |      >>> B.orient(N, 'Axis', [q1, N.x + 2 * N.y])
 |  
 |  orientnew(self, newname, rot_type, amounts, rot_order='', indices=None, latexs=None)
 |      Creates a new ReferenceFrame oriented with respect to this Frame.
 |      
 |      See ReferenceFrame.orient() for acceptable rotation types, amounts,
 |      and orders. Parent is going to be self.
 |      
 |      Parameters
 |      ==========
 |      
 |      newname : str
 |          The name for the new ReferenceFrame
 |      rot_type : str
 |          The type of orientation matrix that is being created.
 |      amounts : list OR value
 |          The quantities that the orientation matrix will be defined by.
 |      rot_order : str
 |          If applicable, the order of a series of rotations.
 |      
 |      
 |      Examples
 |      ========
 |      
 |      >>> from sympy.physics.mechanics import ReferenceFrame, Vector
 |      >>> from sympy import symbols
 |      >>> q1 = symbols('q1')
 |      >>> N = ReferenceFrame('N')
 |      >>> A = N.orientnew('A', 'Axis', [q1, N.x])
 |      
 |      
 |      .orient() documentation:
 |      
 |      ========================
 |      
 |      Defines the orientation of this frame relative to a parent frame.
 |      
 |      Supported orientation types are Body, Space, Quaternion, Axis.
 |      Examples show correct usage.
 |      
 |      Parameters
 |      ==========
 |      
 |      parent : ReferenceFrame
 |          The frame that this ReferenceFrame will have its orientation matrix
 |          defined in relation to.
 |      rot_type : str
 |          The type of orientation matrix that is being created.
 |      amounts : list OR value
 |          The quantities that the orientation matrix will be defined by.
 |      rot_order : str
 |          If applicable, the order of a series of rotations.
 |      
 |      Examples
 |      ========
 |      
 |      >>> from sympy.physics.mechanics import ReferenceFrame, Vector
 |      >>> from sympy import symbols
 |      >>> q0, q1, q2, q3, q4 = symbols('q0 q1 q2 q3 q4')
 |      >>> N = ReferenceFrame('N')
 |      >>> B = ReferenceFrame('B')
 |      
 |      Now we have a choice of how to implement the orientation. First is
 |      Body. Body orientation takes this reference frame through three
 |      successive simple rotations. Acceptable rotation orders are of length
 |      3, expressed in XYZ or 123, and cannot have a rotation about about an
 |      axis twice in a row.
 |      
 |      >>> B.orient(N, 'Body', [q1, q2, q3], '123')
 |      >>> B.orient(N, 'Body', [q1, q2, 0], 'ZXZ')
 |      >>> B.orient(N, 'Body', [0, 0, 0], 'XYX')
 |      
 |      Next is Space. Space is like Body, but the rotations are applied in the
 |      opposite order.
 |      
 |      >>> B.orient(N, 'Space', [q1, q2, q3], '312')
 |      
 |      Next is Quaternion. This orients the new ReferenceFrame with
 |      Quaternions, defined as a finite rotation about lambda, a unit vector,
 |      by some amount theta.
 |      This orientation is described by four parameters:
 |      q0 = cos(theta/2)
 |      q1 = lambda_x sin(theta/2)
 |      q2 = lambda_y sin(theta/2)
 |      q3 = lambda_z sin(theta/2)
 |      Quaternion does not take in a rotation order.
 |      
 |      >>> B.orient(N, 'Quaternion', [q0, q1, q2, q3])
 |      
 |      Last is Axis. This is a rotation about an arbitrary, non-time-varying
 |      axis by some angle. The axis is supplied as a Vector. This is how
 |      simple rotations are defined.
 |      
 |      >>> B.orient(N, 'Axis', [q1, N.x + 2 * N.y])
 |  
 |  set_ang_acc(self, otherframe, value)
 |      Define the angular acceleration Vector in a ReferenceFrame.
 |      
 |      Defines the angular acceleration of this ReferenceFrame, in another.
 |      Angular acceleration can be defined with respect to multiple different
 |      ReferenceFrames. Care must be taken to not create loops which are
 |      inconsistent.
 |      
 |      Parameters
 |      ==========
 |      
 |      otherframe : ReferenceFrame
 |          A ReferenceFrame to define the angular acceleration in
 |      value : Vector
 |          The Vector representing angular acceleration
 |      
 |      Examples
 |      ========
 |      
 |      >>> from sympy.physics.mechanics import ReferenceFrame, Vector
 |      >>> N = ReferenceFrame('N')
 |      >>> A = ReferenceFrame('A')
 |      >>> V = 10 * N.x
 |      >>> A.set_ang_acc(N, V)
 |      >>> A.ang_acc_in(N)
 |      10*N.x
 |  
 |  set_ang_vel(self, otherframe, value)
 |      Define the angular velocity vector in a ReferenceFrame.
 |      
 |      Defines the angular velocity of this ReferenceFrame, in another.
 |      Angular velocity can be defined with respect to multiple different
 |      ReferenceFrames. Care must be taken to not create loops which are
 |      inconsistent.
 |      
 |      Parameters
 |      ==========
 |      
 |      otherframe : ReferenceFrame
 |          A ReferenceFrame to define the angular velocity in
 |      value : Vector
 |          The Vector representing angular velocity
 |      
 |      Examples
 |      ========
 |      
 |      >>> from sympy.physics.mechanics import ReferenceFrame, Vector
 |      >>> N = ReferenceFrame('N')
 |      >>> A = ReferenceFrame('A')
 |      >>> V = 10 * N.x
 |      >>> A.set_ang_vel(N, V)
 |      >>> A.ang_vel_in(N)
 |      10*N.x
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |  
 |  __dict__
 |      dictionary for instance variables (if defined)
 |  
 |  __weakref__
 |      list of weak references to the object (if defined)
 |  
 |  x
 |      The basis Vector for the ReferenceFrame, in the x direction.
 |  
 |  y
 |      The basis Vector for the ReferenceFrame, in the y direction.
 |  
 |  z
 |      The basis Vector for the ReferenceFrame, in the z direction.

The ReferenceFrame class manages everything about rotations, angular velocities, and angular accelerations with respect to other reference frames. Now create an inertial reference frame called N for "Newtonian" as was described in the the docstring:

In [4]:
N = ReferenceFrame('N')

Keep in mind that N is the variable name of which the reference frame named "N" is stored. It is important to note that N is an object and it has properties and functions associated with it. To see a list of them type:

In [5]:
dir(N)
Out[5]:
['__class__',
 '__delattr__',
 '__dict__',
 '__doc__',
 '__format__',
 '__getattribute__',
 '__getitem__',
 '__hash__',
 '__init__',
 '__iter__',
 '__module__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_ang_acc_dict',
 '_ang_vel_dict',
 '_cur',
 '_dcm_dict',
 '_dict_list',
 '_dlist',
 '_w_diff_dcm',
 '_x',
 '_y',
 '_z',
 'ang_acc_in',
 'ang_vel_in',
 'dcm',
 'indices',
 'latex_vecs',
 'name',
 'orient',
 'orientnew',
 'pretty_vecs',
 'set_ang_acc',
 'set_ang_vel',
 'str_vecs',
 'x',
 'y',
 'z']

Notice that three of the properties are x, y, and z. These are the orthonormal unit vectors associated with the reference frame and are the building blocks for creating more general vectors. We can create a vector by simply building a linear combination of the unit vectors:

In [6]:
v = 1 * N.x + 2 * N.y + 3 * N.z

Now a vector expressed in the N reference frame is stored in the variable v. We can print v to the screen by typing:

In [7]:
print(v)
N.x + 2*N.y + 3*N.z

The vector v can be manipulated as expected. You can multiply and divide them by scalars:

In [8]:
2 * v
Out[8]:
2*N.x + 4*N.y + 6*N.z
In [9]:
v / 3.0
Out[9]:
0.333333333333333*N.x + 0.666666666666667*N.y + N.z

Note that three is expressed as 3.0 instead of 3. The python language does integer division by default. There are ways around this, but for now simply remember to always declare numbers as floats (i.e. include a decimal).

You can also add and subtract vectors::

In [10]:
v + v
Out[10]:
2*N.x + 4*N.y + 6*N.z
In [11]:
w = 5 * N.x + 7 * N.y
In [12]:
v - w
Out[12]:
- 4*N.x - 5*N.y + 3*N.z

Vectors also have some useful properties:

In [13]:
dir(v)
Out[13]:
['__add__',
 '__and__',
 '__class__',
 '__delattr__',
 '__dict__',
 '__div__',
 '__doc__',
 '__eq__',
 '__format__',
 '__getattribute__',
 '__hash__',
 '__init__',
 '__module__',
 '__mul__',
 '__ne__',
 '__neg__',
 '__new__',
 '__or__',
 '__radd__',
 '__rand__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__rmul__',
 '__ror__',
 '__rsub__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__sub__',
 '__subclasshook__',
 '__truediv__',
 '__weakref__',
 '__xor__',
 '_latex',
 '_pretty',
 '_sympyrepr',
 '_sympystr',
 'args',
 'cross',
 'diff',
 'doit',
 'dot',
 'dt',
 'express',
 'magnitude',
 'normalize',
 'outer',
 'simp',
 'simplify',
 'subs']

You can find the magnitude of a vector by typing:

In [14]:
v.magnitude()
Out[14]:
sqrt(14)

You can compute a unit vector in the direction of v:

In [15]:
v.normalize()
Out[15]:
sqrt(14)/14*N.x + sqrt(14)/7*N.y + 3*sqrt(14)/14*N.z

You can find the measure numbers and the reference frame the vector was defined in with::

In [16]:
v.args
Out[16]:
[([1]
[2]
[3], N)]

Dot and cross products of vectors can also be computed::

In [17]:
dot(v, w)
Out[17]:
19
In [18]:
cross(v, w)
Out[18]:
- 21*N.x + 15*N.y - 3*N.z

We've only used numbers as our measure numbers so far, but it is just as easy to use symbols. We will introduce six symbols for our measure numbers with the SymPy symbols [help(symbols) for the documentation] function:

In [19]:
a1, a2, a3 = symbols('a1 a2 a3')
b1, b2, b3 = symbols('b1 b2 b3')

And create two new vectors that are completely symbolic:

In [20]:
x = a1 * N.x + a2 * N.y + a3 * N.z
y = b1 * N.x + b2 * N.y + b3 * N.z
In [21]:
dot(x, y)
Out[21]:
a1*b1 + a2*b2 + a3*b3
In [22]:
z = cross(x, y)
z
Out[22]:
(a2*b3 - a3*b2)*N.x + (-a1*b3 + a3*b1)*N.y + (a1*b2 - a2*b1)*N.z

Numbers and symbols work together seamlessly:

In [23]:
dot(v, x)
Out[23]:
a1 + 2*a2 + 3*a3

You can also differentiate a vector with respect to a variable in a reference frame:

In [24]:
z.diff(a1, N)
Out[24]:
- b3*N.y + b2*N.z

Vectors don't have be defined with respect to just one reference frame. We can create a new reference frame and orient it with respect to the N frame that has already been created. We will use the orient method of the new frame to do a simple rotation through alpha about the N.x axis:

In [25]:
A = ReferenceFrame('A')
alpha = symbols('alpha')
A.orient(N, 'Axis', [alpha, N.x])

Now the direction cosine matrix with of A with respect to N can be computed:

In [26]:
A.dcm(N)
Out[26]:
[1,           0,          0]
[0,  cos(alpha), sin(alpha)]
[0, -sin(alpha), cos(alpha)]

Now that SymPy knows that A and N are oriented with respect to each other we can express the vectors that we originally wrote in the A frame:

In [27]:
v.express(A)
Out[27]:
A.x + (3*sin(alpha) + 2*cos(alpha))*A.y + (-2*sin(alpha) + 3*cos(alpha))*A.z
In [28]:
z.express(A)
Out[28]:
(a2*b3 - a3*b2)*A.x + ((a1*b2 - a2*b1)*sin(alpha) + (-a1*b3 + a3*b1)*cos(alpha))*A.y + ((a1*b2 - a2*b1)*cos(alpha) - (-a1*b3 + a3*b1)*sin(alpha))*A.z

In dynamics systems, at least some of the relative orientation of reference frames and vectors are time varying. The mechanics module provides a way to specify quantities as time varying. Let's define two variables beta and d as variables which are functions of time:

In [29]:
beta, d = dynamicsymbols('beta d')
beta, d
Out[29]:
(beta(t), d(t))

Now we can create a new reference frame that is oriented with respect to the A frame by beta and create a vector in that new frame that is a function of d. This time we will use the orientnew method of the A frame to create the new reference frame B:

In [30]:
B = A.orientnew('B', 'Axis', (beta, A.y))
In [31]:
vec = d * B.z

We can now compute the angular velocity of the reference frame B with respect to other reference frames:

In [32]:
B.ang_vel_in(N)
Out[32]:
beta'*A.y

This allows us to now differentiate the vector, vec, with respect to time and a reference frame::

In [33]:
vecdot = vec.dt(N)
vecdot
Out[33]:
d*beta'*B.x + d'*B.z
In [34]:
vecdot.express(N)
Out[34]:
(d*cos(beta)*beta' + sin(beta)*d')*N.x + (d*sin(alpha)*sin(beta)*beta' - sin(alpha)*cos(beta)*d')*N.y + (-d*sin(beta)*cos(alpha)*beta' + cos(alpha)*cos(beta)*d')*N.z

The dynamicsymbols function also allows you to easily create the derivatives of time varying variables and store them in a variable.

In [35]:
theta = dynamicsymbols('theta')
thetad = dynamicsymbols('theta', 1)
theta, thetad
Out[35]:
(theta(t), Derivative(theta(t), t))

At this point we now have all the tools needed to setup the kinematics for a dynamic system. Let's start with the classic mass spring damper system under the influence of gravity. Go to the next notebook to follow along: https://github.com/PythonDynamics/pydy_examples/tree/master/mass_spring_damper.