In [1]:
import time
print(time.ctime())
%matplotlib inline

Wed Oct 26 15:33:09 2022

This notebook contains some additional material on the topic "Simulation of Dynamic Systems" and addresses the following points:

  • What is a dynamical system?
  • How can a dynamical system be described by a differential equation?
  • What is a state space representation?
  • How can differential equations be solved numerically?
  • How can one use Python to simulate a dynamical system?

The document tries to make an appropriate compromise between generality and comprehensibility. In case of doubt it favors comprehensibility so not every technical detail might be 100% accurate.

What is a dynamic system?

A system whose "state" can change over time. Examples:

  • Bathtub with inlet and outlet (state variable: filling level)
  • Cup with tea (state variables: temperature of contents, temperature of cup)
  • Pendulum (state variables: angle of deflection and angular velocity)

How can a dynamic system be described by a differential equation?

If a quantity can change over time, it makes sense to model it mathematically as a function of time. In the bathtub example: $h(t)$.

The velocity or the rate of change of a quantity over time is described by the first derivative. In the example: $\dot h(t)$. This is again a time function. In the process of modeling, equations are derived, e.g. by applying balances and physical laws, which relate the derivatives of the system quantities to the system quantities themselves. In the bathtub example $$\dot h(t) = \frac{1}{A}(q_\mathrm{to}(t) - q_\mathrm{ab}(t)), \tag{1}$$ i.e., the rate of change in water surface elevation is proportional to the difference in inflow and outflow. For simplicity, a cross-sectional area $A$ constant over height is assumed.

The balance equation (1) can be regarded as an odinary differential equation (ODE), because a derivative ($\dot h$) occurs, but a "typical" ODE is only obtained, if for $q_\mathrm{ab}(t)$ we add the pyhsical relation $$q_\mathrm{ab}(t) = c \sqrt{h(t)}$$ (with $c$ depending on the cross section of the outflow opening and the acceleration due to gravity $g$). If one then still assumes the inflow to be constant, i.e., sets $q_\mathrm{to}(t)\equiv q_0$, one obtains $$\dot h(t) = \frac{1}{A}\Big(q_0 - c \sqrt{h(t)}\Big). \tag{2}$$ This ODE states that the rate of change of the level depends on the level itself. Just as the solution of a (scalar) algebraic equation is a number, the solution of a ODE is a function, in our case $h(t)$, the height as a function of time.

Another example is a (mathematical) pendulum with the deflection $\varphi$. From a torque balance around the axis we get (modeling) $$ml^2 \ddot \varphi = - m g l\sin \varphi. \tag{3}$$ This is a 2nd order scalar ODE (because the highest derivative occurring is 2nd order).

When considering ODEs, one can often show mathematically the existence and the uniqueness of the solution - but the analytical calculation of the solution is only possible in special cases, e.g. linear ODEs. Therefore, one usually has to be satisfied with a numerical approximate solution. If the ODE describes a dynamic system, one speaks then also of a simulation of the system.

What is a state space representation?

In order to avoid having to derive a separate numerical solution method for each ODE, it makes sense to start from a representation that covers as many relevant cases as possible. The state space representation $$\dot {\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t))\qquad \text{with the state vector} ~ \mathbf{x}(t) = \left(\begin{matrix} x_1(t) \\ \vdots \\ x_n(t) \end{matrix} \right) \in \mathbb{R}^n \tag{4} $$ is therefore very useful and common.

It is a system of $n$ coupled first order ODEs. Obviously, equation (2) is already in state space representation. It holds $n=1$ and $\mathbf{x} = h$.

The scalar second order ODE (3) can be transformed into a system of two first order ODEs by using a so-called definitional equation (in German: "definitorische Gleichung"): \begin{align} x_1:= \varphi \\ x_2:= \dot \varphi \end{align} From this, solving (3) for $\ddot \varphi = \dot x_2$ yields the state space representation of the pendulum:

\begin{align} \dot{\mathbf{x}}(t) = \left(\begin{matrix} \dot x_1(t) \\ \dot x_2(t) \end{matrix} \right) = \left(\begin{matrix} x_2(t) \\ -\frac{g}{l} \sin x_1 \end{matrix} \right). \tag{5} \end{align}

How to solve a ODE system numerically?

When solving a state space model (5) numerically, one typically assumes that the initial state $\mathbf{x}(0)=:\mathbf{x}_0$ is known and the course of $\mathbf{x}(\cdot)$ on the time interval $[0, T]$ is sought. Therefore one speaks of an initial value problem (IVP) concerning the computation of the solution. Because one calculates the quantities $\mathbf{x}(t)$ from the time derivative of $\mathbf{x}(t)$, one also speaks of numerical integration of a ODE. The simplest numerical integration or solution method is the explicit Euler method. Thereby the rate of change of the state $\mathbf{f}(\mathbf{x})$ is assumed to be constant during the short time interval $\delta t$. One obtains in the first step $$\mathbf{x}(\Delta t) \approx \mathbf{x}_0 + \Delta t \cdot \mathbf{f}(\mathbf{x}_0),$$ i.e., the change of state (i.e., $\mathbf{x}(\Delta t) - \mathbf{x}_0$ ) is equal to $\Delta t$ times (constant) rate of change. This procedure can be continued iteratively: $$\mathbf{x}\big(\Delta t \cdot (k+1)\big) \approx \mathbf{x}(\Delta t \cdot k) + \Delta t \cdot \mathbf{f}\big(\mathbf{x}(\Delta t \cdot k)\big)$$ or in short notation $$\mathbf{x}_{k+1} = \mathbf{x}_{k} + \delta t \cdot \mathbf{f}(\mathbf{x}_k).$$ This does not give the complete solution process $t \mapsto \mathbf{x}(t)$, but only pointwise the approximated function values of the state vector $\mathbf{x_1} := \mathbf{x}(\Delta t), \mathbf{x_2} := \mathbf{x}(2 \Delta t), \ldots, \mathbf{x_N} := \mathbf{x}(N \Delta t)$ with $N \Delta t =T$, but one can show that for sufficiently small values of $\Delta t$ this approximate solution converges to the exact solution - which is known to exist, but generally cannot be written down as a formula.

This explicit Euler method is easy to explain, but has acceptable appoximation accuracy only for small step sizes $\Delta t$. In practice, more complicated methods like Runge-Kutta methods with a step size control are often used. However, from an application point of view, this does not matter much. For most applications one can fall back on already prefabricated numerical integration methods.

Wie kann man mit Python ein dynamisches System simulieren?

In the Python course we use the ODE solver solve_ivp from the package scipy.integrate.

In [2]:
from scipy.integrate import solve_ivp
print(solve_ivp.__doc__[:520]) # print the beginning of the docstring of the solver

Solve an initial value problem for a system of ODEs.

    This function numerically integrates a system of ordinary differential
    equations given an initial value::

        dy / dt = f(t, y)
        y(t0) = y0

    Here t is a 1-D independent variable (time), y(t) is an
    N-D vector-valued function (state), and an N-D
    vector-valued function f(t, y) determines the differential equations.
    The goal is to find y(t) approximately satisfying the differential
    equations, given an initial value y(t0)=y0.

You can get more information with solve_ivp?.

Note: in SciPy the state vector is called $\mathbf{y}$ (instead of $\mathbf{x}$ as above).

The function solve_ivp(rhs_fun, t_span, xx0) expects three obligatory arguments:

  1. a function of the right-hand side of the ODE i.e. $\mathbf{f}(t, \mathbf{x})$.
  2. a sequence t_span = (t_start, t_end) (tuple, list, ...) with start and end time.
  3. the initial state $\mathbf{x}_0$.

Many other optional arguments are allowed, see docs.

The return value is a container object (e.g. called res) which contains the following important attributes:

  • res.success: Binary flag (True (normal case) or False(simulation not successful))
  • res.t: Evaluation times of the solution
  • res.y: Time course of the state vector. Rows: State components; columns: Time points (corresponding to res.t).

For other attributes, see docs.

Notes

  1. the function $\mathbf{f}$ is often called rhs(t, x) in the Python code, because it describes the "right hand side" of the ODE. In equation (4), the function $\mathbf{f}$ depends only on $\mathbf{x}$ and not on $t$. However, the function solve_ivp always expects a Python function accepting two parameters as the first argument. Thus, time-variant ODEs can also be simulated (modified bathtub example: inflow is an explicit function of time). If an ODE is to be simulated which, like the pendulum example, does not depend explicitly on time, the 1st argument inside rhs can be ignored.
  2. it does not matter what the arguments are called. In this text the state vector is called $\mathbf{x}$ or xx in the docstring of solve_ivp it is called y.
  3. besides the mandatory argument t_span, optionally the keyword argument t_eval can be passed as array evaluation times. If this is not done, the integration algorithm decides alone (e.g. based on its step size control).

Example

In the following, the pendulum from equation (5) is simulated for 25 seconds for different initial states ($\varphi(0)= 20^\circ$ as well as $\varphi(0)= 170^\circ$ without and with additional angular velocity).

In [3]:
from numpy import sqrt, sin, array, linspace, pi
import matplotlib.pyplot as plt

g = 9.81
l = 1

# Definition of the right hand side:

def rhs_pendulum(t, xx):
    """
    Right hand side of pendulum ode.
    
    :param t:   time instant (ignored here)
    :param xx:  state vector (1D-array or list or tuple)
    
    :returns:   time derivative of the state vector (as 1D-array)
    """
    x1, x2 = xx
    
    x1dot = x2
    x2dot = -sqrt(g/l)*sin(x1)
    
    return array([x1dot, x2dot])
In [4]:
# Various initial values

xx0_20 = [20*pi/180, 0]
xx0_170a = [170*pi/180, 0]
xx0_170b = [170*pi/180, 0.31]

# evaluation time steps (25s)
tt = linspace(0, 25, 1000)
In [5]:
# Execution of the simulation (with explicit specification of the evaluation times (t_eval)).
# In addition, specification of the maximum step size (max_step) for higher accuracy.
res_20 = solve_ivp(rhs_pendulum, (tt[0], tt[-1]), xx0_20, t_eval=tt, max_step=0.1)
res_170a = solve_ivp(rhs_pendulum, (tt[0], tt[-1]), xx0_170a, t_eval=tt, max_step=0.1)
res_170b = solve_ivp(rhs_pendulum, (tt[0], tt[-1]), xx0_170b, t_eval=tt, max_step=0.1)
In [6]:
# Graphical representation (3 angle curves in one diagram)

plt.figure(figsize=(15, 10))

# Reminder: res.y is 2D array of simulation result (rows: state components; columns: time instants)
plt.plot(tt, res_20.y[0, :], label=r"$\varphi(0)= 20^\circ$")
plt.plot(tt, res_170a.y[0, :], label=r"$\varphi(0)= 170^\circ$")
plt.plot(tt, res_170b.y[0, :], label=r"$\varphi(0)= 170^\circ$, $\dot \varphi > 0$")

plt.legend()
Out[6]:
<matplotlib.legend.Legend at 0x7f9428ddfc10>

Evaluation of the three cases:

At $20°$ initial deflection the pendulum swings approximately harmonically (i.e. $\varphi(t) \approx \varphi(0) \cdot \cos\Big(\sqrt{\frac{g}{l}} \cdot t\Big)$). This solution could also be calculated analytically via a linearization of (3) or (5), because for small deflections $\sin \varphi \approx \varphi$ holds.

At $170°$ initial deflection, the pendulum swings significantly non-harmonically (orange curve looks different from sin or cos function). In addition, the period has increased significantly. (At $180°$ initial deflection it would be infinitely long, because according to the mathematical model the pendulum would remain in the upper rest position).

If a suitable initial angular velocity is added to the $170°$ initial deflection, then the pendulum no longer swings but performs overshoots (green curve).

Final remark: For the nonlinear pendulum ODE (3) an exact analytical solution does exist, see e.g. Wikipedia#Arbitrary-amplitude_period). Nevertheless, the statement that typically and especially for practically relevant systems one cannot find an analytic solution remains valid. That is why numerical integration of ODE systems is a central technique in science and engineering.

In [7]:
import ipydex
ipydex.save_current_nb_as_html()

`simulation_of_dynamical_systems.ipynb` written.