How to Use Symbolic Python to Solve and Visualize Dynamical Systems

Written by Neil Sawhney

If you've worked with some advanced dynamics problems, you'll know that it's often next to impossible to derive the equations of motion by hand.

While most courses will teach you how to solve these kinds of systems using matrices and numerical methods, I'll show you how I use the Sympy.physics.mechanics module in python to solve for the equations of motion symbolically. These can then be simulated by numerical integration.

Below is a simulation I made of a basketball spinning on a finger, which we will use as an example. Feel free to orbit around and change the initial conditions.

Pros of using this method

When you're dealing with vectors, typically you have to be very concious about what frame of reference the vector is expressed in. When using multiple vectors in an equation, you must first make sure to apply the approriate rotation matrices to each vector to ensure that they are expressed in the same frame. However, when using the sympy.physics.mechanics module, vectors are defined in any frame of reference you want. The module will take care of the rotation matrices in the background, allowing you to multiply vectors without worrying about which frame they are expressed in. When you want to print out the result of your equation, you can use the express() funtion, which will print it out in the frame of reference you specify.
Because the sympy.physics.mechanics module uses symbolic python (sympy) you obtain exact solutions instead of numerical approximations.
The sympy.physics.mechanics module includes many useful functions for solving dynamics problems that bypass some of the more tedious steps. For example, the angular velocity of a frame e_b with respect to frame N whose position has been defined with respect to any number of other fully defined frames using e_b.ang_vel_in(N). The direction cosine matrix of frame e_b with respect to N can be easily found by e_b.dcm(N).

Cons of using this method

The sympy.physics.mechanics module has minimal documentation and it can be difficult to find examples of how to use it.
If you're taking a class on advanced dynamics, you should probably learn how to solve these problems by hand first, as it will help you understand the concepts better and catch mistakes if you notice that an equation doesn't make sense.

Some prerequisites/supplies to follow this guide:

A basic understanding of dynamics

A basic understanding of python

The Jupyter Notebook code for this project can be foundhere.
The javascript code for the simulation shown above can be found in the source code (using the developer console of your browser) here.

Step 1: Setup reference frames using a 3-1-3 Euler angle sequence.

\begin{equation}\label{omegabn}\omega_{B/N}^{(b)} = (\sin{\left(\psi{\left(t \right)} \right)} \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} + \cos{\left(\psi{\left(t \right)} \right)} \frac{d}{d t} \theta{\left(t \right)})\mathbf{\hat{e_b}_x} + (- \sin{\left(\psi{\left(t \right)} \right)} \frac{d}{d t} \theta{\left(t \right)} + \sin{\left(\theta{\left(t \right)} \right)} \cos{\left(\psi{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)})\mathbf{\hat{e_b}_y} + (\cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} + \frac{d}{d t} \psi{\left(t \right)})\mathbf{\hat{e_b}_z}\end{equation}

On line 4, we define three dynamic symbols, $\theta$, $\psi$ and $\phi$, this tells sympy that these variables are functions of time.

The code cell above sets up the reference frames that define a 313 euler angle sequence. The interactive window below helps to visualize the sequence.

Change the sequence to 313 using the dropdown to match this tutorial, but feel free to explore how other sequences affect the outcome.

Frame A represents a rotation of the inertial frame about the 3 axis (N.z) by an angle $\phi$,
Frame B represents a rotation of the A frame about the A frame's 1 axis (A.x) by an angle $\theta$,
Frame e_b (the body fixed frame) represents a rotation of the B frame about the B frame's 3 axis (B.z) by an angle $\psi$.

Now that we're set up, notice that sympy allows us to immediately find the angular velocity of any of our defined frames with respect to any other in terms of the three euler angles we defined in a singular line of code. To get the angular velocity of frame e_b with respect to frame N ( $\omega_{B/N}$ ), we simply use omega = e_b.ang_vel_in(N)

If we want to print $\omega_{B/N}$ we can choose to express it in any frame of reference we want. For example, in equation \ref{omegabn}, $\omega_{B/N}$ is expressed in the b frame (hence the (b) that looks like an exponent) using omega.express(e_b). If we want to express it in the inertial frame N instead, we can similarly use omega.express(N).

Step 2: Define the kinematics.

\begin{equation}\label{velocity}\vec{v}_{CM} = r \left(- \sin{\left(\psi{\left(t \right)} \right)} \frac{d}{d t} \theta{\left(t \right)} + \sin{\left(\theta{\left(t \right)} \right)} \cos{\left(\psi{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)}\right)\mathbf{\hat{e_b}_x} - r \left(\sin{\left(\psi{\left(t \right)} \right)} \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} + \cos{\left(\psi{\left(t \right)} \right)} \frac{d}{d t} \theta{\left(t \right)}\right)\mathbf{\hat{e_b}_y}\end{equation}
\begin{equation}\label{acceleration}\vec{a}_{CM} = r \left(\frac{\left(\cos{\left(\psi{\left(t \right)} - 2 \theta{\left(t \right)} \right)} - \cos{\left(\psi{\left(t \right)} + 2 \theta{\left(t \right)} \right)}\right) \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2}}{4} - \sin{\left(\psi{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} \theta{\left(t \right)} + \sin{\left(\theta{\left(t \right)} \right)} \cos{\left(\psi{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} \phi{\left(t \right)} + 2 \cos{\left(\psi{\left(t \right)} \right)} \cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} \frac{d}{d t} \theta{\left(t \right)}\right)\mathbf{\hat{e_b}_x} + r \left(\frac{\left(- \sin{\left(\psi{\left(t \right)} - 2 \theta{\left(t \right)} \right)} + \sin{\left(\psi{\left(t \right)} + 2 \theta{\left(t \right)} \right)}\right) \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2}}{4} - \sin{\left(\psi{\left(t \right)} \right)} \sin{\left(\theta{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} \phi{\left(t \right)} - 2 \sin{\left(\psi{\left(t \right)} \right)} \cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} \frac{d}{d t} \theta{\left(t \right)} - \cos{\left(\psi{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} \theta{\left(t \right)}\right)\mathbf{\hat{e_b}_y} - r \left(\sin^{2}{\left(\theta{\left(t \right)} \right)} \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2} + \left(\frac{d}{d t} \theta{\left(t \right)}\right)^{2}\right)\mathbf{\hat{e_b}_z}\end{equation}

In this example, the kinematics are defined by enforcing the rolling constraint. The velocity of the center of mass of the ball about the point of contact between the ball and the finger is defined to be equal to $\omega_{B/N} \times r_{cm}$. The vector pointing from the point of contact between the finger and the ball to the center of mass of the ball is $r_{cm}$

On lines 4 and 5, We define the position of our points P and CM, we tell sympy that the center of mass, (CM) is located r units away from point of contact (P) in the z direction with respect to the ball's frame of reference. The location of the point P doesn't actually matter since everything else is defined with respect to it.

Line 6 simply defines our $r_{cm}$ vector by stating that it points to point CM from point P.

Lines 7 through 9 is the simplest way to define points P and CM with respect to the previously defined frames of reference.

The point of contact is not moving in the inertial frame because the finger is perfectly still.
The point of contact is also not moving in the body fixed frame, because the ball is always rotating about its bottom most point.
The center of mass is not moving in the body fixed frame... almost by definition.

Line 10 finally enforces the rolling constraint by stating that the velocity of point CM in the intertial frame is equal to $\omega_{B/N} \times r_{cm}$

Once this is done, we can easily obtain the velocity and acceleration of the center of mass with respect to any frame we want and expressed in any frame we want. For example, the velocity and acceleration of the center of mass with respect to the inertial frame as expressed in the body fixed frame is shown in equations \ref{velocity} and \ref{acceleration}. These expression can easily be obtained with the following 2 lines of code. CM.vel(N).express(e_b) and CM.acc(N).express(e_b).

Step 3: Define the kinetics.

Now that the kinematics are set up, to solve the equations of motion, we need to apply Euler's first and second laws.

Lines 3 through 5 simply define the inertia of the ball. The inertia of a sphere about its center of mass is $I_{sphere} = \frac{2}{5} m r^2$. Due to the symmetry of the sphere, the inertia about all three axes (in the body fixed frame e_b) is the same.

Euler's first law is $\sum F = \frac{d}{dt}p = ma$. The force acting on the ball due to gravity is easy to define in sympy. It is found by multiplying the mass of the ball with the gravitational acceleration vector which acts in the negative z direction of the inertial frame.

Once we have the force due to gravity vector defined, it is easy to define the reaction force acting on the ball from the finger. We use the equation $\sum F = Fr + Fg$, since $F = ma$, $F_r = ma - F_g$. The acceleration of the center of mass with respect to the inertial frame can be written using sympy with CM.acc(N). Therefore, the reaction force acting on the ball from the finger can be written as Fr = m*CM.acc(N) - Fg as written on line 10 of the code block.

Euler's second law is $M = \frac{d}{dt}H_c = I \alpha$, therefore $H_c = I \omega$. Since $M = r_p \times F$, we can solve for the equations of motion by setting $\dot{H_c} = M$ and letting sympy figure out the rest.

Lines 13 through 18 are pretty straight forward and define $\dot{H_c}$ and $M$ using the equations above. On line 14, $\dot{H_c}$ is found by telling sympy to differentiate with respect to time taken in the inertial frame.

Step 4: Let sympy solve the problem.

\begin{equation}\label{psi}\ddot{\psi} = - \frac{\left(5 \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} + \frac{2 \frac{d}{d t} \psi{\left(t \right)}}{\tan{\left(\theta{\left(t \right)} \right)}} - \frac{12 \frac{d}{d t} \phi{\left(t \right)}}{\sin{\left(\theta{\left(t \right)} \right)}}\right) \frac{d}{d t} \theta{\left(t \right)}}{7}\end{equation}
\begin{equation}\label{theta}\ddot{\theta} = \frac{5 g \sin{\left(\theta{\left(t \right)} \right)}}{7 r} - \frac{5 \sin{\left(\psi{\left(t \right)} - 3 \theta{\left(t \right)} \right)} \sin{\left(\psi{\left(t \right)} \right)} \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2}}{56 \sin{\left(\theta{\left(t \right)} \right)}} + \frac{5 \sin{\left(\psi{\left(t \right)} + \theta{\left(t \right)} \right)} \sin{\left(\psi{\left(t \right)} \right)} \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2}}{56 \sin{\left(\theta{\left(t \right)} \right)}} - \frac{5 \sin{\left(\psi{\left(t \right)} \right)} \cos{\left(\psi{\left(t \right)} + \theta{\left(t \right)} \right)} \cos{\left(\theta{\left(t \right)} \right)} \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2}}{14} + \frac{5 \sin{\left(\theta{\left(t \right)} \right)} \cos^{2}{\left(\psi{\left(t \right)} \right)} \cos{\left(\theta{\left(t \right)} \right)} \left(\frac{d}{d t} \phi{\left(t \right)}\right)^{2}}{7} - \frac{2 \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} \frac{d}{d t} \psi{\left(t \right)}}{7}\end{equation}
\begin{equation}\label{phi}\ddot{\phi} = - \frac{2 \cdot \left(6 \cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \phi{\left(t \right)} - \frac{d}{d t} \psi{\left(t \right)}\right) \frac{d}{d t} \theta{\left(t \right)}}{7 \sin{\left(\theta{\left(t \right)} \right)}}\end{equation}

Well this part is easy, we hand off the dirty work to sympy and let it figure it out. Sympy solves the problem and the result is easily found using `sol[psi.diff(t,2)]`, `sol[theta.diff(t,2)]`, and `sol[phi.diff(t,2)]` as shown in equations \ref{psi}, \ref{theta}, and \ref{phi}.

Sympy solve doesn't understand the vector notation of the physics.mechanics module (at the time of writing this), so on lines 1-3, I use the dot product to split the vector equation into 3 regular equations.

Step 5: Numerically integrate and simulate.

I wont go as deep into the low level details here because honestly the hard part is done. There are a million well documented ways to numerically integrate. Simulation involves creating an object in any visualization software and using its built in rotate tools in tandem with the numerically integrated solution.

In the jupyter notebook code I provided, I used the odeint function in the scipy python module for numerical integration, and manim to render and display a simulation with preset initial conditions. However, a major flaw in my implementation is that the simulation accuracy relies heavily on the frame rate of the renderer, which makes for long render times and creates some chaotic behavior at low fps.

For the interactive javascript simulation I previewed at the beginning of this article, I used the dopri funtion in the numeric library for the numerical integration. For the visualization I used the three.js library.

However you choose to do it, the process is the same:

Put your numerically integrated solution in an array.
Create the basketball.
Get the interpolated psi, theta, and phi values for the current time in the simulation.
Rotate about the bottom point of the basketball first by phi about the z axis, then by theta about the x axis, then by psi about the z axis again.
(This procedure is shown in the following code snippet using three.js. However, the camera is rotated, so to keep things vertical I'm technically doing a 2-3-2 euler angle sequence instead of a 3-1-3.)

Well, that's basically it! I hope you enjoyed!