<jupyter_start><jupyter_text>Defining the Problem Here we will derive the equations of motion for the classic mass-spring-damper system under the influence of gravity. The following figure gives a pictorial description of the problem.<jupyter_code>from IPython.display import SVG
SVG(filename='mass_spring_damper.svg')<jupyter_output><empty_output><jupyter_text>Start by loading in the core functionality of both SymPy and Mechanics.<jupyter_code>import sympy as sym
import sympy.physics.mechanics as me<jupyter_output><empty_output><jupyter_text>We can make use of the pretty printing of our results by loading SymPy's printing extension, in particular we will use the vector printing which is nice for mechanics objects.<jupyter_code>from sympy.physics.vector import init_vprinting
init_vprinting()<jupyter_output><empty_output><jupyter_text>We'll start by defining the variables we will need for this problem:- $x(t)$: distance of the particle from the ceiling- $v(t)$: speed of the particle- $m$: mass of the particle- $c$: damping coefficient of the damper- $k$: stiffness of the spring- $g$: acceleration due to gravity- $t$: time<jupyter_code>x, v = me.dynamicsymbols('x v')
m, c, k, g, t = sym.symbols('m c k g t')<jupyter_output><empty_output><jupyter_text>Now, we define a Newtonian reference frame that represents the ceiling which the particle is attached to, $C$.<jupyter_code>ceiling = me.ReferenceFrame('C')<jupyter_output><empty_output><jupyter_text>We will need two points, one to represent the original position of the particle which stays fixed in the ceiling frame, $O$, and the second one, $P$ which is aligned with the particle as it moves.<jupyter_code>O = me.Point('O')
P = me.Point('P')<jupyter_output><empty_output><jupyter_text>The velocity of point $O$ in the ceiling is zero.<jupyter_code>O.set_vel(ceiling, 0)<jupyter_output><empty_output><jupyter_text>Point $P$ can move downward in the $y$ direction and its velocity is specified as $v$ in the downward direction.<jupyter_code>P.set_pos(O, x * ceiling.x)
P.set_vel(ceiling, v * ceiling.x)
P.vel(ceiling)<jupyter_output><empty_output><jupyter_text>There are three forces acting on the particle. Those due to the acceleration of gravity, the damper, and the spring.<jupyter_code>damping = -c * P.vel(ceiling)
stiffness = -k * P.pos_from(O)
gravity = m * g * ceiling.x
forces = damping + stiffness + gravity
forces<jupyter_output><empty_output><jupyter_text>Now we can use Newton's second law, $0=F-ma$, to form the equation of motion of the system.<jupyter_code>zero = me.dot(forces - m * P.acc(ceiling), ceiling.x)
zero<jupyter_output><empty_output><jupyter_text>We can then form the first order equations of motion by solving for $\frac{dv}{dt}$ and introducing the kinematical differential equation, $v=\frac{dx}{dt}$.<jupyter_code>dv_by_dt = sym.solve(zero, v.diff(t))[0]
dx_by_dt = v
dv_by_dt, dx_by_dt<jupyter_output><empty_output><jupyter_text>Forming the equations of motion can also be done with the automated methods available in the Mechanics package: `LagrangesMethod` and `KanesMethod`. Here we will make use of Kane's method to find the same equations of motion that we found manually above. First, define a particle that represents the mass attached to the damper and spring.<jupyter_code>mass = me.Particle('mass', P, m)<jupyter_output><empty_output><jupyter_text>Now we can construct a `KanesMethod` object by passing in the generalized coordinate, $x$, the generalized speed, $v$, and the kinematical differential equation which relates the two, $0=v-\frac{dx}{dt}$.<jupyter_code>kane = me.KanesMethod(ceiling, q_ind=[x], u_ind=[v], kd_eqs=[v - x.diff(t)])<jupyter_output><empty_output><jupyter_text>Now Kane's equations can be computed, and we can obtain $F_r$ and $F_r^*$.<jupyter_code>fr, frstar = kane.kanes_equations([(P, forces)], [mass])
fr, frstar<jupyter_output><empty_output><jupyter_text>The equations are also available in the form $M\frac{d}{dt}[q,u]^T=f(q, u)$ and we can extract the mass matrix, $M$, and the forcing functions, $f$.<jupyter_code>M = kane.mass_matrix_full
f = kane.forcing_full
M, f<jupyter_output><empty_output><jupyter_text>Finally, we can form the first order differential equations of motion $\frac{d}{dt}[q,u]^T=M^{-1}f(\dot{u}, u, q)$, which is the same as previously found.<jupyter_code>M.inv() * f<jupyter_output><empty_output><jupyter_text>Simulating the system Now that we have defined the mass-spring-damper system, we are going to simulate it.PyDy's `System` is a wrapper that holds the Kanes object to integrate the equations of motion using numerical values of constants.<jupyter_code>from pydy.system import System
sys = System(kane)<jupyter_output>/home/moorepants/src/pydy/pydy/system.py:102: PyDyFutureWarning: PyDy System is experimental and may change in the future.
  warnings.warn(msg, PyDyFutureWarning)<jupyter_text>Now, we specify the numerical values of the constants and the initial values of states in the form of a dict.<jupyter_code>sys.constants = {m:10.0, g:9.8, c:5.0, k:10.0}
sys.initial_conditions = {x:0.0, v:0.0}<jupyter_output><empty_output><jupyter_text>We must generate a time vector over which the integration will be carried out. NumPy's `linspace` is often useful for this.<jupyter_code>from numpy import linspace
sys.times = linspace(0.0, 10.0, 100)<jupyter_output><empty_output><jupyter_text>The trajectory of the states over time can be found by calling the `.integrate()` method.<jupyter_code>x_trajectory = sys.integrate()<jupyter_output><empty_output><jupyter_text>Visualizing the System PyDy has a native module `pydy.viz` which is used to visualize a System in an interactive 3D GUI.<jupyter_code>from pydy.viz import *<jupyter_output>:0: FutureWarning: IPython widgets are experimental and may change in the future.<jupyter_text>For visualizing the system, we need to create shapes for the objects we wish to visualize, and map each of them to a `VisualizationFrame`, which holds the position and orientation of the object. First create a sphere to represent the bob and attach it to the point $P$ and the ceiling reference frame (the sphere does not rotate with respect to the ceiling).<jupyter_code>bob = Sphere(2.0, color="red", material="metal")
bob_vframe = VisualizationFrame(ceiling, P, bob)<jupyter_output><empty_output><jupyter_text>Now create a circular disc that represents the ceiling and fix it to the ceiling reference frame. The circle's default axis is aligned with its local $y$ axis, so we need to attach it to a rotated ceiling reference frame if we want the circle's axis to align with the $\hat{c}_x$ unit vector.<jupyter_code>ceiling_circle = Circle(radius=10, color="white", material="metal")
from numpy import pi
rotated = ceiling.orientnew("C_R", 'Axis', [pi / 2, ceiling.z])
ceiling_vframe = VisualizationFrame(rotated, O, ceiling_circle)<jupyter_output><empty_output><jupyter_text>Now we initialize a Scene. A Scene contains all the information required to visualize a `System` onto a canvas.It takes a ReferenceFrame and Point as arguments.<jupyter_code>scene = Scene(ceiling, O, system=sys)<jupyter_output>/home/moorepants/src/pydy/pydy/viz/camera.py:73: PyDyUserWarning: Rotation of Perspective Camera does not work properly in the visualiser.
  warnings.warn(msg, PyDyUserWarning)<jupyter_text>We provide the VisualizationFrames, which we want to visualize as a list to scene.<jupyter_code>scene.visualization_frames = [bob_vframe, ceiling_vframe]<jupyter_output><empty_output><jupyter_text>The default camera of Scene has the z axis of the base frame pointing out of the screen, and the y axis pointing up. We want the x axis to point downwards, so we supply a new camera that will achieve this.<jupyter_code>camera_frame = ceiling.orientnew('Camera Frame','Axis', [pi / 2, ceiling.z])
camera_point = O.locatenew('Camera Location', 100 * camera_frame.z)
primary_camera = PerspectiveCamera(camera_frame, camera_point)
scene.cameras = [primary_camera]<jupyter_output><empty_output><jupyter_text>Now, we call the display method.<jupyter_code>scene.display_ipython()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.097347316471314, 'fraction_numerical': 0.005305367057371993, 'mean_word_length': 6.512511584800741, 'pattern_counts': {'":': 0, '<': 106, '<?xml version=': 0, '>': 106, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '158660', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2637, 'n_tokens_neox': 2553, 'n_words': 885, 'path': 'examples/mass_spring_damper/mass_spring_damper.ipynb', 'repo_name': 'Sina-Mehdiz/PyDy', 'size': '8105'}