SymPy Mechanics for Autolev Users

Introduction

Autolev (now superseded by MotionGenesis) is a domain specific programming language which is used for symbolic multibody dynamics. The SymPy mechanics module now has enough power and functionality to be a fully featured symbolic dynamics module. The PyDy package extends the SymPy output to the numerical domain for simulation, analyses and visualization. Autolev and SymPy Mechanics have a lot in common but there are also many differences between them. This page shall expand upon their differences. It is meant to be a go-to reference for Autolev users who want to transition to SymPy Mechanics.

It would be nice to have a basic understanding of SymPy and SymPy Mechanics before going over this page. If you are completely new to Python, you can check out the official Python Tutorial. Check out the SymPy Documentation, especially the tutorial to get a feel for SymPy. For an introduction to Multibody dynamics in Python, this lecture is very helpful.

You might also find the Autolev Parser which is a part of SymPy to be helpful.

Some Key Differences

Autolev

SymPy Mechanics
Autolev is a domain specific programming language designed to perform multibody dynamics. Since it is a language of its own, it has a very rigid language specification. It predefines, assumes and computes many things based on the input code. Its code is a lot cleaner and concise as a result of this.
SymPy is a library written in the general purpose language Python. Although Autolev’s code is more compact, SymPy (by virtue of being an add on to Python) is more flexible. The users have more control over what they can do. For example, one can create a class in their code for let’s say a type of rigibodies with common properties. The wide array of scientific Python libraries available is also a big plus.
Autolev generates Matlab, C, or Fortran code from a small set of symbolic mathematics.
SymPy generates numerical Python, C or Octave/Matlab code from a large set of symbolic mathematics created with SymPy. It also builds on the popular scientific Python stack such as NumPy, SciPy, IPython, matplotlib, Cython and Theano.
Autolev uses 1 (one) based indexing. The initial element of a sequence is found using a[1].
Python uses 0 (zero) based indexing. The initial element of a sequence is found using a[0].
Autolev is case insensitive.
SymPy code being Python code is case sensitive.
One can define their own commands in Autolev by making .R and .A files which can be used in their programs.
SymPy code is Python code, so one can define functions in their code. This is a lot more convenient.
Autolev is proprietary.
SymPy is open source.

Rough Autolev-SymPy Equivalents

The tables below give rough equivalents for some common Autolev expressions. These are not exact equivalents, but rather should be taken as hints to get you going in the right direction. For more detail read the built-in documentation on SymPy vectors, SymPy mechanics and PyDy .

In the tables below, it is assumed that you have executed the following commands in Python:

import sympy.physics.mechanics as me
import sympy as sm

Mathematical Equivalents

Autolev

SymPy

Notes
Constants A, B
a, b = sm.symbols(‘a b’, real=True)
Note that the names of the symbols can be different from the names of the variables they are assigned to. We can define a, b = symbols(‘b a’) but its good practice to follow the convention.
Constants C+
c = sm.symbols(‘c’, real=True, nonnegative=True)
Refer to SymPy assumptions for more information.
Constants D-
d = sm.symbols(‘d’, real=True, nonpositive=True)

Constants K{4}
k1, k2, k3, k4 = sm.symbols('k1 k2 k3 k4', real=True)

Constants a{2:4}
a2, a3, a4 = sm.symbols('a2 a3 a4', real=True)

Constants b{1:2, 1:2}
b11, b12, b21, b22 = sm.symbols('b11 b12 b21 b22', real=True)

Specified Phi
phi = me.dynamicsymbols(‘phi ')

Variables q, s
q, s = me.dynamicsymbols(q, s)

Variables x’’

x = me.dynamicsymbols(‘x’ )

xd = me.dynamicsymbols(‘x’ , 1)

xd2 = me.dynamicsymbols(‘x’ , 2)

Variables y{2}’

y1 = me.dynamicsymbols(‘y1’ )

y2 = me.dynamicsymbols(‘y2’ )

y1d = me.dynamicsymbols(‘y1’ , 1)

y2d = me.dynamicsymbols(‘y2' , 1)

MotionVariables u{2}

u1 = me.dynamicsymbols(‘u1’ )

u2 = me.dynamicsymbols('u2' )

SymPy doesn’t differentiate between variables, motionvariables and specifieds during declaration. Instead, it takes different lists of these as parameters in objects like the KanesMethod.

Imaginary j

j = sm.I

I is a sympy object which stands for the imaginary unit. One can define complex numbers using it.

z = x + I*y

where x, y and z are symbols.

Tina = 2*pi

s = u*t + a*t^2/2

tina = 2*sm.pi

tina = tina.evalf()

t = me.dynamicsymbols._t

s = u*t + a*t**2/2

Using .evalf() will result in the numeric value.
abs(x)^3 + sin(x)^2 + acos(x)
sm.abs(x)**3 + sm.sin(x)**2 + sm.acos(x)

E = (x+2*y)^2 + 3*(7+x)*(x+y)

Expand(E)

Factor(E, x)

Coef(y, x)

Replace(y, sin(x)=3)

Exclude(E,x)

Include(E,x)

Arrange(E,2,y)

E = (x+2*y)**2 + 3*(7+x)*(x+y)

sm.expand(E)

sm.horner(E, wrt=x)

y.coeff(x)

y.subs({sm.sin(x): 3})

e.collect(x).coeff( x, 0)

e.collect(x).coeff( x, 1)

e.collect(y)

For more information refer to simplification.

These SymPy functions do not work in place. They just return expressions. If you want to overwrite the original expression you would have to do something like:

y = y.subs({sm.sin(x): 3})

Dy = D(E, y)

Dt = Dt(E)

Dt2 = Dt(V, A) where V is a vector and A is a frame

Dy2 = D(V, y, A)

E.diff(y)

E.diff( me.dynamicsymbols._t )

Works if the expression is made up of dynamicsymbols.

dt2 = v.dt(A)

dy2 = v.diff(y, A)

For more information refer to calculus.

E = COS(X*Y)

TY = Taylor(E, 0:2, x=0, y=0)

e = sm.cos(x*y)

b = e.series(x, 0, 2).removeO().series(y, 0, 2).removeO()

For more information refer to series.

F = Evaluate(E, x=a, y=2)

E.subs([(x, a), (y, 2)])

To get floating point numbers from numerical expressions use .evalf()

E.evalf((a + sm.pi).subs({a: 3}))

P = Polynomial([a, b, c], x)

p = sm.Poly(sm.Matrix([a, b, c]).reshape(1, 3), x)

For more information refer to polys.

Roots(Polynomial( a*x^2 + b*x + c, x, 2)

Roots([1;2;3])

sm.solve( sm.Poly(a*x**2 + b*x + c))

sm.solve(sm.Poly( sm.Matrix([1,2,3]). reshape(3, 1), x), x)

For more information refer to Solvers.

For numerical computation related to polynomials and roots refer to mpmath/calculus.

Solve(A, x1, x2)

where A is an augmented matrix that represents the linear equations and x1, x2 are the variables to solve for.

sm.linsolve(A, (x1, x2))

where A is an augmented matrix

For more information refer to :ref:` solvers/solveset. <solveset>`

For non linear solvers refer to nonlinsolve and nsolve in solvers.

RowMatrix = [1, 2, 3, 4]

ColMatrix = [1; 2; 3; 4]

MO = [a, b; c, 0]

MO[2, 2] := d

A + B*C

Cols(A)

Cols(A, 1)

Rows(A)

Rows(A, 1)

Det(A)

Element(A, 2, 3)

Inv(A)

Trace(A)

Transpose(A)

Diagmat(4, 1)

Eig(A)

Eig(A, EigVal, EigVec)

row_matrix = sm.Matrix([[1],[2], [3],[4]])

col_matrix = sm.Matrix([1, 2, 3, 4])

MO = sm.Matrix([[a, b], [c, 0]])

MO[1, 1] = d

A + B*C

A.cols

A.col(0)

A.rows

A.row(0)

M.det()

M[2, 3]

M**-1

sm.trace(A)

A.T

sm.diag(1,1,1,1)

A.eigenvals()

eigval = A.eigenvals()

eigvec = A.eigenvects()

For more information refer to matrices.

Physical Equivalents

Autolev

SymPy

Notes

Bodies A

Declares A, its masscenter Ao, and orthonormal vectors A1>, A2> and A3> fixed in A.

m =sm.symbols(‘m’)

Ao = sm.symbols(‘Ao’)

Af = me.ReferenceFrame(‘Af’ )

I = me.outer(Af.x,Af.x)

P = me.Point(‘P’)

A =me.RigidBody(‘A’, Ao, Af, m, (I, P))

Af.x, Af.y and Af.z are equivalent to A1>, A2> and A3>.

The 4th and 5th arguments are for the mass and inertia. These are specified after the declaration in Autolev.

One can pass a dummy for the parameters and use setters A.mass = \_ and A.inertia = \_ to set them later.

For more information refer to mechanics/masses .

Frames A

V1> = X1*A1> + X2*A2>

A = me.ReferenceFrame(‘A’ )

v1 = x1*A.x + x2*A.y

For more information refer to physics/vectors.

Newtonian N

N = me.ReferenceFrame(‘N’ )

SymPy doesn’t specify that a frame is inertial during declaration. Many functions such as set_ang_vel() take the inertial reference frame as a parameter.

Particles C

m = sm.symbols(‘m’)

Po = me.Point(‘Po’)

C = me.Particle(‘C’, Po, m)

The 2nd and 3rd arguments are for the point and mass. In Autolev, these are specified after the declaration.

One can pass a dummy and use setters (A.point = \_ and A.mass = \_) to set them later.

Points P, Q

P = me.Point(‘P’)

Q = me.Point(‘Q’)

Mass B=mB

mB = symbols(‘mB’)

B.mass = mB

Inertia B, I1, I2, I3, I12, I23, I31

I = me.inertia(Bf, i1, i2, i3, i12, i23, i31)

B.inertia = (I, P) where B is a rigidbody, Bf is the related frame and P is the center of mass of B.

Inertia dyadics can also be formed using vector outer products.

I = me.outer(N.x, N.x)

For more information refer to the mechanics api.

vec> = P_O_Q>/L

vec> = u1*N1> + u2*N2>

Cross(a>, b>)

Dot(a>, b>)

Mag(v>)

Unitvec(v>)

DYAD>> = 3*A1>*A1> + A2>*A2> + 2*A3>*A3>

vec  = (Qo.pos_from(O))/L

vec = u1*N.x + u2*N.y

cross(a, b)

dot(a, b)

v.magnitude()

v.normalize()

dyad = 3*me.outer(a.x ,a.x) + me.outer(a.y, a.y) + 2*me.outer(a.z ,a.z)

For more information refer to physics/vectors.

P_O_Q> = LA*A1>

P_P_Q> = LA*A1>

Q.point = O.locatenew(‘Qo’, LA*A.x)

where A is a reference frame.

Q.point = P.point.locatenew(‘Qo ’, LA*A.x)

For more information refer to the kinematics api.

All these vector and kinematic functions are to be used on Point objects and not Particle objects so .point must be used for particles.

V_O_N> = u3*N.1> + u4*N.2>

Partials(V_O_N>, u3)

O.set_vel(N, u1*N.x + u2*N.y)

O.partial_velocity(N , u3)

The getter would be O.vel(N).

A_O_N> = 0>

Acceleration of point O in reference frame N.

O.set_acc(N, 0)

The getter would be O.acc(N).

W_B_N> = qB’*B3>

Angular velocity of body B in reference frame F.

B.set_ang_vel(N, qBd*Bf.z)

where Bf is the frame associated with the body B.

The getter would be B.ang_vel_in(N).

ALF_B_N> =Dt(W_B_N>, N)

Angular acceleration of body B in reference frame N.

B.set_ang_acc(N, diff(B.ang_vel_in(N) )

The getter would be B.ang_acc_in(N).

Force_O> = F1*N1> + F2*N2>

Torque_A> = -c*qA’*A3>

In SymPy one should have a list which contains all the forces and torques.

fL.append((O, f1*N.x + f2*N.y))

where fL is the force list.

fl.append((A, -c*qAd*A.z))

A_B = M where M is a matrix and A, B are frames.

D = A_B*2 + 1

B.orient(A, 'DCM', M) where M is a SymPy Matrix.

D = A.dcm(B)*2 + 1

CM(B)

B.masscenter

Mass(A,B,C)

A.mass + B.mass + C.mass

V1pt(A,B,P,Q)

Q.v1pt_theory(P, A, B)

P and Q are assumed to be Point objects here. Remember to use .point for particles.

V2pts(A,B,P,Q)

Q.v2pt_theory(P, A, B)

A1pt(A,B,P,Q)

Q.a1pt_theory(P, A, B)

A2pts(A,B,P,Q)

Q.a2pt_theory(P, A, B)

Angvel(A,B)

B.ang_vel_in(A)

Simprot(A, B, 1, qA)

B.orient(A, ‘Axis’, qA, A.x)

Gravity(G*N1>)

fL.extend(gravity( g*N.x, P1, P2, ...))

In SymPy we must use a forceList (here fL) which contains tuples of the form (point, force_vector). This is passed to the kanes_equations() method of the KanesMethod object.

CM(O,P1,R)

me.functions. center_of_mass(o, p1, r)

Force(P/Q, v>)

fL.append((P, -1*v), (Q, v))

Torque(A/B, v>)

fL.append((A, -1*v), (B, v))

Kindiffs(A, B ...)

KM.kindiffdict()

Momentum(option)

linear_momentum(N, B1, B2 ...)

reference frame followed by one or more bodies

angular_momentum(O, N, B1, B2 ...)

point, reference frame followed by one or more bodies

KE()

kinetic_energy(N, B1, B2 ...)

reference frame followed by one or more bodies

Constrain(...)

velocity_constraints = [...]

u_dependent = [...]

u_auxiliary = [...]

These lists are passed to the KanesMethod object.

For more details refer to mechanics/kane and the kane api.

Fr() FrStar()

KM = KanesMethod(f, q_ind, u_ind, kd_eqs, q_dependent, configura tion_constraints, u_de pendent, velocity_cons traints, acceleration_ constraints, u_auxilia ry)

The KanesMethod object takes a reference frame followed by multiple lists as arguments.

(fr, frstar) = KM.kanes_equations(fL, bL) where fL and bL are lists of forces and bodies respectively.

For more details refer to mechanics/kane and the kane api.

Numerical Evaluation and Visualization

Autolev’s CODE Option() command allows one to generate Matlab, C, or Fortran code for numerical evaluation and visualization. Option can be Dynamics, ODE, Nonlinear or Algebraic.

Numerical evaluation for dynamics can be achieved using PyDy. One can pass in the KanesMethod object to the System class along with the values for the constants, specifieds, initial conditions and time steps. The equations of motion can then be integrated. The plotting is achieved using matlplotlib. Here is an example from the PyDy Documentation on how it is done:

from numpy import array, linspace, sin
from pydy.system import System

sys = System(kane,
             constants = {mass: 1.0, stiffness: 1.0,
                          damping: 0.2, gravity: 9.8},
             specifieds = {force: lambda x, t: sin(t)},
             initial_conditions = {position: 0.1, speed:-1.0},
             times = linspace(0.0, 10.0, 1000))

y = sys.integrate()

import matplotlib.pyplot as plt
plt.plot(sys.times, y)
plt.legend((str(position), str(speed)))
plt.show()

For information on all the things PyDy can accomplish refer to the PyDy Documentation.

The tools in the PyDy workflow are :

  • SymPy: SymPy is a Python library for

    symbolic computation. It provides computer algebra capabilities either as a standalone application, as a library to other applications, or live on the web as SymPy Live or SymPy Gamma.

  • NumPy: NumPy is a library for the

    Python programming language, adding support for large, multi-dimensional arrays and matrices, along with a large collection of high-level mathematical functions to operate on these arrays.

  • SciPy: SciPy is an open source

    Python library used for scientific computing and technical computing. SciPy contains modules for optimization, linear algebra, integration, interpolation, special functions, FFT, signal and image processing, ODE solvers and other tasks common in science and engineering.

  • IPython: IPython is a command shell

    for interactive computing in multiple programming languages, originally developed for the Python programming language, that offers introspection, rich media, shell syntax, tab completion, and history.

  • Aesara: Aesara is

    a numerical computation library for Python. In Aesara, computations are expressed using a NumPy-esque syntax and compiled to run efficiently on either CPU or GPU architectures.

  • Cython: Cython is a superset of the

    Python programming language, designed to give C-like performance with code that is mostly written in Python. Cython is a compiled language that generates CPython extension modules.

  • matplotlib: matplotlib is a

    plotting library for the Python programming language and its numerical mathematics extension NumPy.

One will be able to write code equivalent to the Matlab, C or Fortran code generated by Autolev using these scientific computing tools. It is recommended to go over these modules to gain an understanding of scientific computing with Python.