Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Project: Czwartek
Views: 271
Image: ubuntu1804
Kernel: Python 3 (system-wide)

Sympy initiation and dynamic system definition

In the next block, necessary modules are called and dynamic system is defined by its kinetic and potential energy.

from sympy import * from sympy.physics.mechanics import * mechanics_printing(pretty_print=True) t=Symbol('t') m_1,m_2,l_m,l_M,k,l,c,u_0,Omega = symbols('m_1,m_2,l_m,l_M,k,l,c,u_0,Omega',positive=True) g=Symbol('g',positive=True) x, phi = dynamicsymbols('x_1,phi') q=[x,phi] u=[x.diff(t),phi.diff(t)] Y=q+u x_2=x+l*sin(phi) y_2=l*cos(phi) v_2=sqrt((x_2.diff(t)**2+y_2.diff(t)**2)) T=(m_1*x.diff(t)**2)/2+(m_2*v_2**2)/2 V=-m_2*g*y_2+(k*(x-u_0*sin(Omega*t))**2)/2
Eq(Symbol('T'),T),Eq(Symbol('V'),V)

(T=m1x˙122+m2(l2sin2(ϕ)ϕ˙2+(lcos(ϕ)ϕ˙+x˙1)2)2, V=glm2cos(ϕ)+k(u0sin(Ωt)+x1)22)\displaystyle \left( T = \frac{m_{1} \dot{x}_{1}^{2}}{2} + \frac{m_{2} \left(l^{2} \operatorname{sin}^{2}\left(\phi\right) \dot{\phi}^{2} + \left(l \operatorname{cos}\left(\phi\right) \dot{\phi} + \dot{x}_{1}\right)^{2}\right)}{2}, \ V = - g l m_{2} \operatorname{cos}\left(\phi\right) + \frac{k \left(- u_{0} \operatorname{sin}\left(\Omega t\right) + x_{1}\right)^{2}}{2}\right)

System Lagrangian

Lagrange function is given by the difference of kinetic and potential energy

L=T-V #Lagrange function Eq(Symbol('L'),L.expand().simplify())

L=glm2cos(ϕ)ku02sin2(Ωt)2+ku0x1sin(Ωt)kx122+l2m2ϕ˙22+lm2cos(ϕ)ϕ˙x˙1+m1x˙122+m2x˙122\displaystyle L = g l m_{2} \operatorname{cos}\left(\phi\right) - \frac{k u_{0}^{2} \operatorname{sin}^{2}\left(\Omega t\right)}{2} + k u_{0} x_{1} \operatorname{sin}\left(\Omega t\right) - \frac{k x^{2}_{1}}{2} + \frac{l^{2} m_{2} \dot{\phi}^{2}}{2} + l m_{2} \operatorname{cos}\left(\phi\right) \dot{\phi} \dot{x}_{1} + \frac{m_{1} \dot{x}_{1}^{2}}{2} + \frac{m_{2} \dot{x}_{1}^{2}}{2}

Determining of the governing equations

In order to obtain equations of motion, the Lagrangian Mechanics was applied. Reduction of the differential equations system allows to obtain the state space equations of the dynamic system. The right hand sides of system under consideration were found in this way.

First method

Direct differentiation of the Lagrange function (L)

L_mat=Matrix([L.expand().doit().simplify()]) # matrix form of the scalar function L which enables .jacobian() method - (1x1) matrix eoms=(L_mat.jacobian(u).diff(t) - L_mat.jacobian(q)).T eoms_rhs=solve(eoms,Matrix(q).diff(t,t)) state_space=Matrix(u+[eoms_rhs[comp] for comp in Matrix(q).diff(t,t)]) ode_rhs=state_space.expand().applyfunc(simplify) #Determine the right hand sides of the govergning equations in order to use it for solve_ivp() integrating procedure. This function solves the initial value problem with utilization of numerical algorithms (as RK45 ord LSODA) implemented within scipy module ode_rhs

[x˙1ϕ˙gm2sin(2ϕ)2+ku0sin(Ωt)kx1+lm2sin(ϕ)ϕ˙2m1+m2sin2(ϕ)gm1sin(ϕ)gm2sin(ϕ)ku0sin(Ωt)cos(ϕ)+kx1cos(ϕ)lm2sin(2ϕ)ϕ˙22l(m1+m2sin2(ϕ))]\displaystyle \left[\begin{matrix}\dot{x}_{1}\\\dot{\phi}\\\frac{\frac{g m_{2} \operatorname{sin}\left(2 \phi\right)}{2} + k u_{0} \operatorname{sin}\left(\Omega t\right) - k x_{1} + l m_{2} \operatorname{sin}\left(\phi\right) \dot{\phi}^{2}}{m_{1} + m_{2} \operatorname{sin}^{2}\left(\phi\right)}\\\frac{- g m_{1} \operatorname{sin}\left(\phi\right) - g m_{2} \operatorname{sin}\left(\phi\right) - k u_{0} \operatorname{sin}\left(\Omega t\right) \operatorname{cos}\left(\phi\right) + k x_{1} \operatorname{cos}\left(\phi\right) - \frac{l m_{2} \operatorname{sin}\left(2 \phi\right) \dot{\phi}^{2}}{2}}{l \left(m_{1} + m_{2} \operatorname{sin}^{2}\left(\phi\right)\right)}\end{matrix}\right]

Second method

Application of the LagrangesMethod class.

LM=LagrangesMethod(L.expand().simplify().doit(),q) eqns=LM.form_lagranges_equations() LM.rhs().expand().applyfunc(simplify) #do nothing, just for the testing purposes ode_rhs=LM.rhs().expand().applyfunc(simplify) #Determine the right hand sides of the govergning equations in order to use it for solve_ivp() integrating procedure. This function solves the initial value problem with utilization numerical algorithms of (as RK45 ord LSODA) implemented within scipy module ode_rhs

[x˙1ϕ˙gm2sin(2ϕ)2+ku0sin(Ωt)kx1+lm2sin(ϕ)ϕ˙2m1+m2sin2(ϕ)gm1sin(ϕ)gm2sin(ϕ)ku0sin(Ωt)cos(ϕ)+kx1cos(ϕ)lm2sin(2ϕ)ϕ˙22l(m1+m2sin2(ϕ))]\displaystyle \left[\begin{matrix}\dot{x}_{1}\\\dot{\phi}\\\frac{\frac{g m_{2} \operatorname{sin}\left(2 \phi\right)}{2} + k u_{0} \operatorname{sin}\left(\Omega t\right) - k x_{1} + l m_{2} \operatorname{sin}\left(\phi\right) \dot{\phi}^{2}}{m_{1} + m_{2} \operatorname{sin}^{2}\left(\phi\right)}\\\frac{- g m_{1} \operatorname{sin}\left(\phi\right) - g m_{2} \operatorname{sin}\left(\phi\right) - k u_{0} \operatorname{sin}\left(\Omega t\right) \operatorname{cos}\left(\phi\right) + k x_{1} \operatorname{cos}\left(\phi\right) - \frac{l m_{2} \operatorname{sin}\left(2 \phi\right) \dot{\phi}^{2}}{2}}{l \left(m_{1} + m_{2} \operatorname{sin}^{2}\left(\phi\right)\right)}\end{matrix}\right]

Numeriacal form of right hand sides of the equations of motion

In order to determine the numerical model, regular Python function has to be created.

import numpy as np #a library for the numerical computations #dictionary (something like a list) with the given system parameters,there won't be possibility to change when the integrator (solving procedure) will be called; all of the system invariants might be declared here case_data={} # .subs() substitute the data to the equations ode_rhs_num=ode_rhs.subs(case_data) odes_numpy=lambdify( (t,Y,m_1,m_2,l,k,g,u_0,Omega),ode_rhs_num) #numerical form of rhs - creates an input for solver - sympy lambdify. It is regular Python function defined by the conversion of Sympy expresion to Python code test_result=odes_numpy(60,range(4),2,4,8,10,10,20,10) #exemplary result computed in order to verify generated expression. Just for the testing purposes. # initial conditons y_0=[0,0,0,0] t_span=np.linspace(0,7,350) #creates time span and assigns it to t_span variable. It is the list of the evalutation points where the solution will be calculated.

Numerical simulations

Application of the appropriate integrating procedure enables to obtain a numerical solution of the sytem under consideration. The solve_ivp code from scipy.integrate module was used in this case.

#importing the necessary libraries for integrating, ploting and data handling from scipy import integrate import matplotlib.pyplot as plt import pandas as pd #default values of the system parameters m_1_num=8 m_2_num=12 l_num=30 k_num=2000 g_num=10 u_0_num=0.1 Omega_num=40 for m_1_num in [4,8,16,32,64,128]: y=integrate.solve_ivp(odes_numpy, [t_span[0],t_span[-1]],y_0,method='LSODA',vectorized=False,t_eval=t_span,args=(m_1_num,m_2_num,l_num,k_num,g_num,u_0_num,Omega_num)) #integrates the given system of the ODEs. plt.plot(y.t,y.y[3,:], y.t,y.y[2,:],y.t,y.y[1,:],y.t,y.y[0,:]) #plots the calculated solution plt.grid(True) #enables grid plt.xlabel("t") #shows the x axis label plt.ylabel("A") #shows the y axis label plt.title("Time history") #shows the plot name plt.show() #shows the scene with plots df_sol=pd.DataFrame(data=y.y.T,index=np.round(y.t,2)).round(3)#saves the results, The variable y.y is a matrix of computed values. The matrix y.y.T is transposition of t #version with labels df_sol=pd.DataFrame(data={'Y'+str(no):y.y[no,:] for no,key in enumerate(Y)},index=np.round(y.t,2)).round(3)#saves the results, The variable y.y is a matrix of computed values. The matrix y.y.T is transposition of t df_sol.to_csv('./DATA/results_'+str(m_1_num) +'.csv')
Image in a Jupyter notebookImage in a Jupyter notebookImage in a Jupyter notebookImage in a Jupyter notebookImage in a Jupyter notebookImage in a Jupyter notebook
df_sol
Y0 Y1 Y2 Y3
0.00 0.000 0.0 0.000 0.000
0.02 0.000 -0.0 0.012 -0.000
0.04 0.001 -0.0 0.040 -0.001
0.06 0.002 -0.0 0.068 -0.002
0.08 0.003 -0.0 0.077 -0.003
... ... ... ... ...
6.92 0.007 -0.0 -0.062 0.002
6.94 0.006 -0.0 -0.043 0.001
6.96 0.006 -0.0 -0.015 0.000
6.98 0.005 -0.0 0.006 -0.000
7.00 0.006 -0.0 0.004 -0.000

350 rows × 4 columns