pydykit - A Python-based Dynamics toolKIT

Repository	https://github.com/pydykit/pydykit
Copyright	Copyright © 2024 - 2025 The authors of pydykit

Welcome to

What is pydykit?

pydykit is a Python-based dynamics simulation toolkit for dynamical systems. The package is based on time stepping methods, which are discrete versions of the corresponding dynamics equations - either ordinary differential equations (ODEs) or differential-algebraic equations (DAEs).

How to install?

Directly install from the Python Package Index (PyPI) with

pip install pydykit

or clone the repository and install the package locally, following these steps:

1. Create a new virtual environment and activate it. We recommend using venv:

   python3.12 -m venv .venv
   source .venv/bin/activate
   

2. Clone the repository
3. Install in editable- / development-mode:

   pip install --editable .
   

4. Run your first script, e.g.

   python scripts/s*.py
   

How to use?

Check out the getting started to familiarize yourself with pydykit and examples for more details.

Who are we?

• Julian K. Bauer - Code architect - @JulianKarlBauer
• Philipp L. Kinon - Core developer - @plkinon

Getting Started

A pydykit-simulation contains several building blocks, which are at least a system, integrator, simulator and a time_stepper. These building blocks are defined within a config file.

Config File

A pydykit-config file configures a simulation devided into several building blocks. Each building blocks consists of parameters and methods. Example config files can be found in pydykit/example_files and some of them are discussed in more detail in the section examples.

System

The system defines the system of ordinary differential equations (ODEs) or differential algebraic equations (DAEs), to be solved. Pydykit covers three system families, which are:

1. Multibody systems
2. Port-Hamiltonian systems
3. Quasilinear DAEs

Lets briefly introduce them one by one.

1. Multibody system (MBS) dynamics

This family covers DAEs belonging to mechanical systems composed of multiple rigid bodies, governed by

\[ \begin{align} \dot q &= M^{-1} p \\ \dot p &= -\nabla V(q) - D(q) M^{-1} p - \nabla g(q)^{\mathrm{T}} \lambda \\ 0 &= g(q) \end{align} \]

where \(q\) and \(p\) are generalized coordinates and momenta, respectively, \(M\) is the mass matrix, \(V\) is the potential energy, \(D\) is the Rayleigh dissipation matrix and \(g\) are independent, holonomic constraint functions enforced by Lagrange multipliers \(\lambda\). One example can be studied in more detail, the 3D pendulum.

2. Port-Hamiltonian systems (PHS)

Port-Hamiltonian systems are an extension of classical Hamiltonian dyamics with respect to dissipative effects and interconnection terms from the environment. This energy-based approach combines dynamics from various physical disciplines in one formalism and emerges from control systems theory. The dynamics follow

\[ \begin{align} E(x) \dot x &= \left( J(x) -R(x) \right) z(x) + B(x)u \\ E(x)^{\mathrm{T}} z(x) &= \nabla H(x) \\ y &= B(x)^{\mathrm{T}} z(x) \end{align} \]

with state \(x\), co-state function \(z\), possibly singular descriptor matrix \(E\), structure matrix \(J=-J^{\mathrm{T}}\), positive semi-definite dissipatrix matrix \(R=R^{\mathrm{T}}\), input matrix \(B\), control input \(u\) with collocated output \(y\) and Hamiltonian \(H\).

3. Quasilinear DAEs

This last framework is an even more general one. It comprises all DAEs of the general form

\[ E(x) \dot x = f(x,t) \]

where \(E\) is a possibly singular descriptor matrix, \(x\) are the unknowns/states and \(f\) denotes an arbitrary right-hand side function. Examples that can be studied in more detail, are the Lorenz attractor and the chemical reactor.

Integrator

The integrator defines the numerical integration scheme to be used for the system at hand. As the integration scheme defines whether a numerical procedure preserves the underlying structure of the problem, pydykit places particular emphasis on integrator classed. Although a system may fit into even two or three of the abovementioned framework, the user might want to introduce specialized integrators to capture as much of the underlying structure as possible, e.g. exactly preserve the constraint functions of the MBS or the skew-symmetry of the structure matrix of a PHS.

Simulator

The simulator defines the solution procedure to be used. This includes the equation solver to be used, e.g. Newton's method and its characeteristics such as accuracy.

Time Stepper

The time_stepper defines the discrete temporal grid, where approximative solutions are computed. Currently, pydykit supports one-step methods only, i.e. for a known solution at one time-instance pydykit computes the solution at the next time instance on the specified grid.

Examples

Pendulum 3D

Let's use the system class pydykit.systems_multi_body.ParticleSystem to simulate a single particle with concentrated mass within a three-dimensional space.

The particle's initial position is (1.0, 0.0, 0.0) and its initial velocity points into the y-direction. The particle's distance to a fixed support at (0.0, 0.0, 0.0) is constraint to be of length 1.0 and there is a gravitational field into negative z-direction.

In consequence, the particle is called a 3D pendulum and its motion is visualized in belows Result tab. The source code of leading to this viosualization is given within the Source tab.

ResultSource

import plotly.graph_objects as go

from pydykit.configuration import Configuration
from pydykit.examples import ExampleManager
from pydykit.managers import Manager
from pydykit.plotters import Plotter

####### Simulation
file_content = ExampleManager().get_example(name="pendulum_3d")
configuration = Configuration(**file_content)  # Validate config file content

manager = Manager()
manager.configure(configuration=configuration)

result = manager.manage()  # Run the simulation

####### Plotting
df = result.to_df()
fig = go.Figure()
plotter = Plotter(results_df=df)
for index in range(manager.system.nbr_particles):
    plotter.plot_3d_trajectory(
        figure=fig,
        x_components=df[f"position0_particle{index}"],
        y_components=df[f"position1_particle{index}"],
        z_components=df[f"position2_particle{index}"],
        time=df["time"],
    )
    index_time = 0
    plotter.add_3d_annotation(
        figure=fig,
        x=df[f"position0_particle{index}"][index_time],
        y=df[f"position1_particle{index}"][index_time],
        z=df[f"position2_particle{index}"][index_time],
        text=f"start of particle {index}",
    )
plotter.fix_scene_bounds_to_extrema(figure=fig, df=df)

####### Show visualization
script_is_executed_on_local_machine = False  # You might want to change this

if script_is_executed_on_local_machine:
    fig.show()
else:
    print(
        fig.to_html(
            full_html=False,
            include_plotlyjs="cdn",
        )
    )

Config File

Let's have a closer look at the configuration file of this simulation:

name: pendulum_3d
system:
  class_name: ParticleSystem
  nbr_spatial_dimensions: 3
  particles:
    - index: 0
      initial_position: [1.0, 0.0, 0.0]
      initial_momentum: [0.0, 1.0, 0.0]
      mass: 1.0
  supports:
    - index: 0
      type: fixed
      position: [0.0, 0.0, 0.0]
  springs: []
  dampers: []
  constraints:
    - start:
        type: support
        index: 0
      end:
        type: particle
        index: 0
      length: 1.0
  gravity: [0.0, 0.0, -9.81]
integrator:
  class_name: MidpointMultibody
simulator:
  class_name: OneStep
  solver_name: NewtonPlainPython
  newton_epsilon: 1.e-07
  max_iterations: 40
time_stepper:
  class_name: FixedIncrementHittingEnd
  step_size: 0.08
  start: 0.0
  end: 1.3

System

The system-section in the above configuration file defines the scene, aka. system, to be simulated. This includes definition of the particle, the fixed support, the constraint and the gravitation.

The system-sections variable class_name tells pydykit that the scene shall be based on the system class ParticleSystem, which belongs to the family of MBS. This class is known to pydykit as it has been registered within the SystemFactory.

Integrator

The integrator-section in the above configuration file defines the integration scheme to be used. Here, the implicit midpoint rule will be applied to the system which is formulated as a MBS.

Similar to the registration pattern of the system, the variable class_name within section integrator tells pydykit to use the class MidpointMultibody. This class is known to pydykit as it has been registered within the IntegratorFactory. The pattern of referencing a registered Python class in terms of class_name also applies to the sections simulator and time_stepper.

Simulator

The simulator-section defines the solution procedure, aka. simulator, to be used. The simulator uses a Newton method with a specific accuracy for the norm of the residual newton_epsilon and a maximum number of iterations per time step before the simulation procedure is stopped.

Time Stepper

The time_stepper-section defines the time stepping algorithm based on settings for start time, end time and step size.

Lorenz System

Let's use the system class pydykit.systems_dae.Lorenz to simulate the famous Lorenz attractor, which is modelled in terms a quasilinear differential algebraic equations.

ResultSource

import plotly.graph_objects as go

from pydykit.configuration import Configuration
from pydykit.examples import ExampleManager
from pydykit.managers import Manager
from pydykit.plotters import Plotter

####### Simulation
file_content = ExampleManager().get_example(name="lorenz")
configuration = Configuration(**file_content)  # Validate config file content

manager = Manager()
manager.configure(configuration=configuration)

result = manager.manage()  # Run the simulation

####### Plotting
fig = go.Figure()
df = result.to_df()
plotter = Plotter(results_df=df)
plotter.plot_3d_trajectory(
    figure=fig,
    x_components=df["x"],
    y_components=df["y"],
    z_components=df["z"],
    time=df["time"],
)
plotter.fix_scene_bounds_to_extrema(figure=fig, df=df)

####### Show visualization
script_is_executed_on_local_machine = False  # You might want to change this

if script_is_executed_on_local_machine:
    fig.show()
else:
    print(
        fig.to_html(
            full_html=False,
            include_plotlyjs="cdn",
        )
    )

Config File

The configuration file of this simulation reads as

name: lorenz
system:
  class_name: Lorenz
  sigma: 10.0
  rho: 28.0
  beta: 2.6666666667
  state:
    state: [2.0, 1.0, 1.0]
simulator:
  class_name: OneStep
  solver_name: NewtonPlainPython
  newton_epsilon: 1.e-07
  max_iterations: 40
integrator:
  class_name: MidpointDAE
time_stepper:
  class_name: FixedIncrement
  step_size: 0.02
  start: 0.0
  end: 4.0

For details on the usage of class_name-variable, see the example Pendulum 3D.

Chemical Reactor

Let's use the system class pydykit.systems_dae.ChemicalReactor to represent a chemical reactor in terms of a differential-algebraic system and solve it's behavior in time. This system is modelled as quasilinear differential-algebraic equations.

ResultSource

import plotly.graph_objects as go

from pydykit.configuration import Configuration
from pydykit.examples import ExampleManager
from pydykit.managers import Manager
from pydykit.plotters import Plotter

####### Simulation
file_content = ExampleManager().get_example(name="reactor")
configuration = Configuration(**file_content)  # Validate config file content

manager = Manager()
manager.configure(configuration=configuration)

result = manager.manage()  # Run the simulation

####### Plotting
fig = go.Figure()
df = result.to_df()
plotter = Plotter(results_df=df)
plotter.plot_3d_trajectory(
    figure=fig,
    x_components=df["concentration"],
    y_components=df["temperature"],
    z_components=df["reaction_rate"],
    time=df["time"],
)

####### Show visualization
script_is_executed_on_local_machine = False  # You might want to change this

if script_is_executed_on_local_machine:
    fig.show()
else:
    print(
        fig.to_html(
            full_html=False,
            include_plotlyjs="cdn",
        )
    )

Config File

The configuration file of this simulation reads as

name: reactor
system:
  class_name: ChemicalReactor
  constants: [1.0, 1.0, 1.0, -100.0]
  cooling_temperature: 10.0
  reactant_concentration: 0.5
  initial_temperature: 50.0
  state:
    state: [0.0, 50.0, 0.0]
simulator:
  class_name: OneStep
  solver_name: NewtonPlainPython
  newton_epsilon: 1.e-10
  max_iterations: 40
integrator:
  class_name: MidpointDAE
time_stepper:
  class_name: FixedIncrement
  step_size: 0.01
  start: 0.0
  end: 1

For details on the usage of class_name-variable, see the example Pendulum 3D.

List Examples

Examples shipped with pydykit can accessed with pydykit.examples.ExampleManager. This is demonstrated by the list of available examples within the tab Result shown below. The source code generating this list, is shown in tab Source.

ResultSource

• four_particle_system_discrete_gradient_dissipative config file
• four_particle_system_midpoint config file
• four_particle_system_ph_discrete_gradient_dissipative config file
• four_particle_system_ph_midpoint config file
• lorenz config file
• pendulum_2d config file
• pendulum_3d config file
• reactor config file
• rigid_body_rotating_quaternion config file
• two_particle_system config file
• visco_pendulum config file

from markdown import markdown

from pydykit import examples

example_manager = examples.ExampleManager()
keys = example_manager.examples.keys()

base_url = example_manager.BASE_URL_EXAMPLE_FILES
html = markdown(
    "\n".join([f"- {item} [config file]({base_url}{item}.yml)" for item in keys])
)
print(html)

Ended: Examples