Mechpy Tutorials

a mechanical engineering toolbox

source code - https://github.com/nagordon/mechpy
documentation - https://nagordon.github.io/mechpy/web/

Neal Gordon
2017-02-20

Solving a second order differential eqaution symbolically with python

In [17]:
# setup 
from sympy import Eq, pi
import sympy as sp
import matplotlib.pyplot as plt
from numpy import linspace
sp.init_printing(use_latex='mathjax')
get_ipython().magic('matplotlib inline') # inline plotting
In [18]:
t,k,m,c = sp.symbols('t,k,m,c')
x = sp.Function('x') # (t)
In [19]:
k_ = 1e3  # spring constant, kN/m
m_ = 50 # mass, Kg
c_ = 3  # damping coefficient
In [20]:
ode = k*x(t) + m*c*x(t).diff(t,1) + m*x(t).diff(t,2)
Eq(ode)
Out[20]:
$$c m \frac{d}{d t} x{\left (t \right )} + k x{\left (t \right )} + m \frac{d^{2}}{d t^{2}} x{\left (t \right )} = 0$$
In [21]:
ode_sol = sp.dsolve(ode)
ode_sol
Out[21]:
$$x{\left (t \right )} = C_{1} e^{\frac{t}{2} \left(- c - \frac{1}{m} \sqrt{m \left(c^{2} m - 4 k\right)}\right)} + C_{2} e^{\frac{t}{2} \left(- c + \frac{1}{m} \sqrt{m \left(c^{2} m - 4 k\right)}\right)}$$
In [22]:
x0 = 1.0
v0 = 0

bc1 = Eq(ode_sol.lhs.subs(x(t),x0), ode_sol.rhs.subs(t,0))
bc2 = Eq(ode_sol.lhs.subs(x(t),v0), ode_sol.rhs.diff(t).subs(t,0))

C_eq = {bc1,bc2}
C_eq
Out[22]:
$$\left\{0 = C_{1} \left(- \frac{c}{2} - \frac{1}{2 m} \sqrt{m \left(c^{2} m - 4 k\right)}\right) + C_{2} \left(- \frac{c}{2} + \frac{1}{2 m} \sqrt{m \left(c^{2} m - 4 k\right)}\right), 1.0 = C_{1} + C_{2}\right\}$$
In [23]:
known_params = {m,c,k,t}
const = ode_sol.free_symbols - known_params
const
Out[23]:
$$\left\{C_{1}, C_{2}\right\}$$
In [24]:
Csol = sp.solve(C_eq,const)
Csol
Out[24]:
$$\left \{ C_{1} : - \frac{0.5 c m}{\sqrt{c^{2} m^{2} - 4.0 k m}} + 0.5, \quad C_{2} : \frac{0.5 c m}{\sqrt{c^{2} m^{2} - 4.0 k m}} + 0.5\right \}$$
In [25]:
ode_sol = ode_sol.subs(Csol)
ode_sol
Out[25]:
$$x{\left (t \right )} = \left(- \frac{0.5 c m}{\sqrt{c^{2} m^{2} - 4.0 k m}} + 0.5\right) e^{\frac{t}{2} \left(- c - \frac{1}{m} \sqrt{m \left(c^{2} m - 4 k\right)}\right)} + \left(\frac{0.5 c m}{\sqrt{c^{2} m^{2} - 4.0 k m}} + 0.5\right) e^{\frac{t}{2} \left(- c + \frac{1}{m} \sqrt{m \left(c^{2} m - 4 k\right)}\right)}$$
In [26]:
ode_sol = ode_sol.subs({m:m_, c:c_, k:k_})
ode_sol
Out[26]:
$$x{\left (t \right )} = \left(0.5 + 0.178017248729078 i\right) e^{\frac{t}{2} \left(-3 - 8.42614977317636 i\right)} + \left(0.5 - 0.178017248729078 i\right) e^{\frac{t}{2} \left(-3 + 8.42614977317636 i\right)}$$
In [27]:
#sp.plot(ode_sol.rhs, (t,0,5)) ;

xfun = sp.lambdify(t,ode_sol.rhs, "numpy")
vfun = sp.lambdify(t,sp.diff(ode_sol.rhs), "numpy")

t = linspace(0,5,1000)

fig, ax1 = plt.subplots(figsize=(12,8))
ax2 = ax1.twinx()
ax1.plot(t,xfun(t),'b',label = r'$x (mm)$', linewidth=2.0)
ax2.plot(t,vfun(t),'g--',label = r'$\dot{x} (m/sec)$', linewidth=2.0)
ax2.legend(loc='lower right')
ax1.legend()
ax1.set_xlabel('time , sec')
ax1.set_ylabel('disp (mm)',color='b')
ax2.set_ylabel('velocity (m/s)',color='g')
plt.title('Mass-Spring System with $v_0=0.1%f' % (v0))
plt.grid()
plt.show()

C:\Users\neal\Anaconda3\lib\site-packages\numpy\core\numeric.py:474: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Solving a second order differential equation numerically

In [28]:
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.pyplot import plot,xlabel,ylabel,title,legend,figure,subplots
from numpy import cos, pi, arange, sqrt, pi, array
get_ipython().magic('matplotlib inline') # inline plotting
In [29]:
def MassSpringDamper(state,t):
    '''
    k=spring constant, Newtons per metre
    m=mass, Kilograms
    c=dampign coefficient, Newton*second / meter    

    for a mass, spring, damper 
        xdd = -k*x/m -c*xd
    '''
  
    k=1e3  # spring constant, kN/m
    m=50 # mass, Kg
    c=3  # damping coefficient 
    # unpack the state vector
    x,xd = state # displacement,x and velocity x'
    # compute acceleration xdd = x''
    xdd = -k*x/m -c*xd
    return [xd, xdd]
In [30]:
x0 = 1.0
v0 = 0
state0 = [x0, v0]  #initial conditions [x0 , v0]  [m, m/sec] 
ti = 0.0  # initial time
tf = 5.0  # final time
step = 0.001  # step
t = arange(ti, tf, step)
state = odeint(MassSpringDamper, state0, t)
x = array(state[:,[0]])
xd = array(state[:,[1]])
In [31]:
# Plotting displacement and velocity
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 14

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(t,x,'b',label = r'$x (mm)$', linewidth=2.0)
ax2.plot(t,xd,'g--',label = r'$\dot{x} (m/sec)$', linewidth=2.0)
ax2.legend(loc='lower right')
ax1.legend()
ax1.set_xlabel('time , sec')
ax1.set_ylabel('disp (mm)',color='b')
ax2.set_ylabel('velocity (m/s)',color='g')
plt.title('Mass-Spring System with $v_0=%0.1f$' % (v0))
plt.grid()
plt.show()
In [ ]: