Equation of motion with sympy

It seems that sympy, a computer algebra library in python, can handle equations of motion that appear in physics, so I would like to try it. This time, I will write a script to derive the equation of motion of the mass-spring-damper system, which is also in pydy_example on the reference site. First, suppose you have a mass spring damper system as shown in the figure. If you formulate the equation of motion normally, it will be as follows. Equation of motion

{\bf M}\dot{{\bf x}} = {\bf f}({\bf x})

here

{\bf x} = [{\it x}, \dot{{\it x}}]^T \\
{\bf M} = \begin{bmatrix} 1 & 0 \\ 0 & m \end{bmatrix} \\
{\bf f}({\bf x}) = \begin{bmatrix} \dot{x} \\ -k x - c \dot{x} + m g + f \end{bmatrix}

Suppose that Next is the script. We use a module called sympy.physics.mechanics. I'm using the Kane's method to formulate the equation of motion, but I'm not sure about this. There seems to be a way to solve it with the Lagrange method.

`mass_spring_dumper.py`


#!/usr/bin/python
#coding:utf-8
import sympy as sym
import sympy.physics.mechanics as me

x, v = me.dynamicsymbols('x v')
m, c, k, g, t, f = sym.symbols('m c k g t f')

#Creating a coordinate system
ceiling = me.ReferenceFrame('C')

o = me.Point('o') #Ceiling point
p = me.Point('p') #Mass point
o.set_vel(ceiling, 0)
p.set_pos(o, x * ceiling.x)
p.set_vel(ceiling, v * ceiling.x)

#Calculate the external force applied to the mass point
damping = -c * p.vel(ceiling) #Dumping
stiffness = -k * p.pos_from(o) #Spring
gravity = m * g * ceiling.x #gravity
exf = f * ceiling.x #Other external forces

forces = damping + stiffness + gravity + exf
print forces

mass = me.Particle('mass', p, m)
kane = me.KanesMethod(ceiling, q_ind=[x], u_ind=[v], kd_eqs=[v - x.diff(t)])
kane.kanes_equations([(p, forces)], [mass])
M = kane.mass_matrix_full
f = kane.forcing_full
print M
print f
print M.inv() * f

The output looks like this:

(-c*v + f + g*m - k*x)*C.x
Matrix([[1, 0], [0, m]])
Matrix([[v(t)], [-c*v(t) + f + g*m - k*x(t)]])
Matrix([[v(t)], [(-c*v(t) + f + g*m - k*x(t))/m]])

If you learn how to use it, it will be a very powerful tool. There are maple and mathematica for a fee, but it is expensive to use such software for free, so I'm grateful that such software can be used for free.

Reference site

sympy documentation http://docs.sympy.org/latest/index.html pydy examples https://github.com/PythonDynamics/pydy_examples