Informally, a system is a part of the physical or abstract world that we isolate for study, and whose behavior in time we want to predict or influence. Examples include electrical circuits, mechanical structures, vehicles, and even population models. \( \text{System dynamics} \) is the study of how such systems evolve over time in response to inputs and initial conditions.

In control engineering, we are usually interested in:

• Mapping \( u(\cdot) \) (input signals) and initial state \( x_0 \) to trajectories of the state \( x(\cdot) \) and output \( y(\cdot) \).
• Understanding qualitative properties such as stability, transient response, and steady-state behavior.
• Designing controllers that modify the input to achieve desired behavior.

Conceptually, we can think of the workflow from a real system to a mathematical model and then to analysis and design as:

flowchart TD
  P["Physical system"] --> M["Modeling: assumptions and laws"]
  M --> EQ["Mathematical model (typically ODEs)"]
  EQ --> AN["Analysis (stability, response)"]
  EQ --> SI["Simulation (numerical solution)"]
  AN --> DE["Design of controllers"]
  SI --> DE
        

The rest of this lesson turns this intuition into precise mathematical concepts, restricted to settings that rely only on your existing knowledge of calculus, linear algebra, and differential equations.

Fix a time axis \( \mathcal{T} \subset \mathbb{R} \), typically \( \mathcal{T} = [t_0,\infty) \). Let \( \mathcal{U} \) be a space of admissible input functions \( u : \mathcal{T} \to \mathbb{R}^m \), and \( \mathcal{Y} \) a space of output functions \( y : \mathcal{T} \to \mathbb{R}^p \). A (deterministic) dynamic system can be seen abstractly as an operator

\[ \mathcal{S} : \mathbb{R}^n \times \mathcal{U} \to \mathcal{Y}, \quad (x_0, u(\cdot)) \mapsto y(\cdot), \]

where \( x_0 \in \mathbb{R}^n \) is the initial condition (the initial state). The key point is that for dynamic systems, \( y(t) \) generally depends on the entire past history of \( u(\cdot) \) and on \( x_0 \), not only on the instantaneous value \( u(t) \).

In most of this course, we will represent such operators via state-space models of the form:

\[ \dot{x}(t) = f\bigl(t, x(t), u(t)\bigr), \quad y(t) = g\bigl(t, x(t), u(t)\bigr), \]

where:

• \( x(t) \in \mathbb{R}^n \) is the state vector,
• \( u(t) \in \mathbb{R}^m \) is the input vector,
• \( y(t) \in \mathbb{R}^p \) is the output vector,
• \( f \) and \( g \) are given functions describing the system's dynamics and output.

For the moment, we treat \( f \) and \( g \) as sufficiently smooth functions so that the initial value problem has a unique solution; this will be made more precise later in the lesson.

The notion of state is central to system dynamics. Intuitively, the state at time \( t \) contains exactly the information about the past that is needed to predict the future evolution, together with the future inputs.

Formally, consider the state-space model above. For fixed \( f \), input \( u(\cdot) \), and initial time \( t_0 \), there exists a (local) flow map \( \Phi \) such that

\[ x(t) = \Phi\bigl(t, t_0, x_0, u(\cdot)\bigr), \quad t \in \mathcal{T}. \]

The Markov property of the state dynamics can be expressed as

\[ x(t_2) = \Phi\bigl(t_2, t_1, x(t_1), u(\cdot)\bigr) \quad \text{for all } t_2 \ge t_1 \ge t_0, \]

meaning that once \( x(t_1) \) and the future input are known, the specific history \( x(\tau), u(\tau) \) for \( t_0 \le \tau < t_1 \) no longer needs to be stored.

In this course, we will almost always work with \( x(t) \) as a finite-dimensional vector in \( \mathbb{R}^n \) (lumped-parameter models). Distributed-parameter systems with infinite-dimensional state (e.g., PDEs) appear in later chapters.

A very common class of models in control engineering is the linear time-invariant (LTI) system. In its simplest scalar form, this can be written as

\[ a_1 \dot{x}(t) + a_0 x(t) = b_0 u(t), \]

where \( a_1 \neq 0 \). Dividing both sides by \( a_1 \) gives

\[ \dot{x}(t) = -\alpha x(t) + \beta u(t), \quad \alpha = \frac{a_0}{a_1}, \quad \beta = \frac{b_0}{a_1}. \]

This already fits the state-space form with scalar state: \( x(t) \in \mathbb{R} \), \( u(t) \in \mathbb{R} \), and \( f(t,x,u) = -\alpha x + \beta u \).

For the homogeneous case (no input), \( u(t) \equiv 0 \), we obtain

\[ \dot{x}(t) = -\alpha x(t), \quad x(t_0) = x_0. \]

This ODE can be solved using separation of variables:

\[ \frac{\dot{x}(t)}{x(t)} = -\alpha \quad \Rightarrow \quad \int_{t_0}^{t} \frac{\dot{x}(\xi)}{x(\xi)}\,\mathrm{d}\xi = -\alpha (t - t_0). \]

By the chain rule, the left-hand side is \( \ln |x(t)| - \ln |x(t_0)| \), so

\[ \ln |x(t)| - \ln |x_0| = -\alpha (t - t_0) \quad \Rightarrow \quad x(t) = x_0 \exp\bigl(-\alpha (t - t_0)\bigr). \]

This exponential behavior is fundamental to the study of system dynamics. In later chapters, we will relate the sign and magnitude of \( \alpha \) to stability and speed of response.

The state-space model defines an initial value problem (IVP). A fundamental theoretical question is: under what conditions does a solution \( x(t) \) exist and is unique for given input \( u(\cdot) \) and initial state \( x_0 \)?

Consider the IVP

\[ \dot{x}(t) = f\bigl(t, x(t)\bigr), \quad x(t_0) = x_0. \]

(For simplicity, we treat the input \( u \) as fixed and absorb it into the function \( f \).) A standard result from the theory of ODEs is the Picard–Lindelöf theorem:

\[ \text{If } f(t,x) \text{ is continuous in } t \text{ and locally Lipschitz in } x, \text{ then there exists } \delta > 0 \\ \text{ such that the IVP has a unique solution on } [t_0-\delta,t_0+\delta]. \]

Here, locally Lipschitz in \( x \) means that for each compact set \( K \subset \mathbb{R}^n \), there exists \( L_K \ge 0 \) such that

\[ \bigl\| f(t,x_1) - f(t,x_2) \bigr\| \le L_K \bigl\| x_1 - x_2 \bigr\| \quad \text{for all } x_1,x_2 \in K, \; t \in \mathcal{T}. \]

For linear systems with continuous coefficients, these conditions are satisfied, guaranteeing a unique trajectory for each initial state and input. This is the mathematical backbone behind the idea that a well-posed dynamic model produces deterministic, reproducible behavior.

In practice, we rarely obtain closed-form expressions like the exponential solution above for complex systems. Instead, we approximate trajectories numerically. A simple way is the (explicit) Euler method applied to \( \dot{x}(t) = f(t,x(t),u(t)) \) with step size \( h > 0 \):

\[ x_{k+1} = x_k + h \, f\bigl(t_k, x_k, u(t_k)\bigr), \quad t_k = t_0 + k h. \]

The sequence \( \{x_k\} \) approximates the continuous trajectory \( x(t_k) \). This mapping from equations to algorithms is what software tools like MATLAB, SciPy, and others implement in much more sophisticated forms (Runge–Kutta methods, implicit schemes for stiff systems, etc.).

The remainder of this lesson illustrates how to implement a simple first-order system

\[ \dot{x}(t) = a x(t) + b u(t), \quad x(0) = x_0, \quad u(t) \equiv 1, \]

using widely used scientific computing environments.

For Python, typical libraries used in system dynamics and control are:

• numpy for vector and matrix operations,
• scipy.integrate (e.g. solve_ivp) for ODE solvers,
• python-control (control package) for state-space and transfer-function models (used extensively in later chapters).

The following script simulates \( \dot{x} = a x + b \) with \( a = -1 \), \( b = 1 \), \( x(0) = 0 \) using solve_ivp and compares it with the analytical solution \( x(t) = 1 - \exp(-t) \).

import numpy as np
from scipy.integrate import solve_ivp

# Parameters
a = -1.0
b = 1.0

def rhs(t, x):
    """
    Right-hand side of the ODE: dx/dt = a*x + b*u, with u(t) = 1.
    x is a 1D array of length 1 here.
    """
    u = 1.0  # unit-step input
    return a * x + b * u

# Time horizon
t0 = 0.0
tf = 5.0
t_eval = np.linspace(t0, tf, 201)

# Initial condition
x0 = np.array([0.0])

sol = solve_ivp(rhs, (t0, tf), x0, t_eval=t_eval)

# Analytical solution for comparison
x_analytical = 1.0 - np.exp(a * t_eval)  # since a = -1

# Print a few sample points
for t, x_num, x_ana in zip(t_eval[::50], sol.y[0, ::50], x_analytical[::50]):
    print(f"t = {t:.2f}, numerical x = {x_num:.5f}, analytical x = {x_ana:.5f}")
      

Later, we will introduce the control library to construct transfer functions and state-space objects for the same system and to analyze its stability and frequency response.

In C++, a common choice for ODE integration is boost::numeric::odeint. It provides generic Runge–Kutta and other schemes. For our scalar example, we can write:

#include <iostream>
#include <boost/numeric/odeint.hpp>

using state_type = double;

// dx/dt = a * x + b * u, with u(t) = 1
void rhs(const state_type &x, state_type &dxdt, const double t)
{
    const double a = -1.0;
    const double b = 1.0;
    const double u = 1.0;
    dxdt = a * x + b * u;
}

int main()
{
    state_type x = 0.0;                 // initial condition
    double t0 = 0.0;
    double tf = 5.0;
    double dt = 0.01;

    // integrate using a simple stepper
    boost::numeric::odeint::integrate(rhs, x, t0, tf, dt);

    std::cout << "x(" << tf << ") = " << x << std::endl;
    return 0;
}
      

For multi-dimensional systems, one would typically use std::vector<double> or Eigen vectors as the state type, and define the right-hand side accordingly. More advanced usage (adaptive step sizes, implicit methods) will be relevant for stiff system dynamics in later chapters.

In Java, the org.apache.commons.math3.ode package (Apache Commons Math) provides ODE solvers such as Runge–Kutta schemes. Below is a minimal implementation of the same scalar system:

import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.FirstOrderIntegrator;
import org.apache.commons.math3.ode.nonstiff.ClassicalRungeKuttaIntegrator;

public class SimpleSystem implements FirstOrderDifferentialEquations {

    private final double a;
    private final double b;

    public SimpleSystem(double a, double b) {
        this.a = a;
        this.b = b;
    }

    @Override
    public int getDimension() {
        return 1; // scalar state
    }

    @Override
    public void computeDerivatives(double t, double[] x, double[] dxdt) {
        double u = 1.0; // unit-step input
        dxdt[0] = a * x[0] + b * u;
    }

    public static void main(String[] args) {
        double a = -1.0;
        double b = 1.0;
        SimpleSystem system = new SimpleSystem(a, b);

        FirstOrderIntegrator integrator =
                new ClassicalRungeKuttaIntegrator(0.01);

        double t0 = 0.0;
        double tf = 5.0;
        double[] x0 = new double[]{0.0};

        integrator.integrate(system, t0, x0, tf, x0);
        System.out.println("x(" + tf + ") = " + x0[0]);
    }
}
      

This object-oriented view emphasizes that the system dynamics are encoded in a model class, which is then passed to a generic integrator. Conceptually, this mirrors the abstract operator view introduced earlier in the lesson.

MATLAB provides high-level ODE solvers such as ode45. The same scalar system can be simulated as:

% Parameters
a = -1.0;
b = 1.0;

% ODE right-hand side
f = @(t, x) a * x + b * 1.0;   % u(t) = 1

tspan = [0 5];
x0 = 0.0;

[t, x] = ode45(f, tspan, x0);

% Analytical solution for comparison
x_analytical = 1 - exp(a * t);

% Plot
figure;
plot(t, x, 'o', t, x_analytical, '-');
xlabel('t');
ylabel('x(t)');
legend('Numerical', 'Analytical');
grid on;
      

In Simulink, the same dynamics are typically represented by an integrator block implementing \( \dot{x} = a x + b u \), with a gain block a on the feedback path, a gain block b on the input u, and a summing junction. This block-diagram view is a graphical counterpart of the mathematical and programmatic descriptions above, and will be formalized in later chapters on block diagrams and state-space modeling.

In Wolfram Mathematica, NDSolve is the standard tool for solving IVPs. The scalar system is expressed almost identically to its mathematical form:

a = -1.0;
b = 1.0;

(* u[t] = 1 *)
sol = NDSolve[
  {
    x'[t] == a x[t] + b*1,
    x[0] == 0
  },
  x,
  {t, 0, 5}
];

Plot[Evaluate[x[t] /. sol], {t, 0, 5},
  AxesLabel -> {"t", "x(t)"},
  PlotLegends -> {"x(t)"}
]
      

Mathematica is particularly convenient for symbolic manipulations (e.g., Laplace transforms, eigenvalue computations) and for visualizing state trajectories and phase portraits, which we will leverage in later lessons dealing with nonlinear and multi-dimensional system dynamics.

Problem 1 (Exponential Solution of a First-Order System). Consider the homogeneous system \( \dot{x}(t) = a x(t) \) with \( a \in \mathbb{R} \), \( x(0) = x_0 \). Derive the solution and discuss qualitatively what happens for \( a < 0 \), \( a = 0 \), and \( a > 0 \).

Solution. Separating variables and integrating, as in Section 4, we obtain

\[ \frac{\dot{x}(t)}{x(t)} = a \quad \Rightarrow \quad \int_{0}^{t} \frac{\dot{x}(\xi)}{x(\xi)} \,\mathrm{d}\xi = a t \quad \Rightarrow \quad \ln|x(t)| - \ln|x_0| = a t. \]

Exponentiating both sides yields

\[ x(t) = x_0 \exp(a t). \]

If \( a < 0 \), then \( \exp(a t) \to 0 \) as \( t \to \infty \), so \( x(t) \) decays to zero (asymptotic stability of the equilibrium). If \( a = 0 \), the solution is constant \( x(t) = x_0 \). If \( a > 0 \), \( \exp(a t) \) grows without bound and \( x(t) \) diverges unless \( x_0 = 0 \). These distinctions foreshadow the stability theory developed in later chapters.

Problem 2 (Static vs. Dynamic Input–Output Relations). Let \( u(t) \in \mathbb{R} \) be an input and \( y(t) \in \mathbb{R} \) an output. Consider the relations:

1. \( y(t) = 3 u(t) \).
2. \( \dot{y}(t) = -2 y(t) + u(t) \), \( y(0) = y_0 \).

Decide which relation is static and which is dynamic, and justify your answer.

Solution. In (1), for each time \( t \), \( y(t) \) depends only on the current value \( u(t) \). There is no dependence on past values or on an initial condition, so this is a static relation.

In (2), the differential equation together with the initial condition implies that \( y(t) \) at time \( t \) depends on the entire past input \( u(\tau) \) for \( 0 \le \tau \le t \) and on \( y_0 \). Hence (2) represents a dynamic system. This illustrates at a simple level how system dynamics extends static input–output relations.

Problem 3 (State Representation of a Second-Order System). Consider the second-order ODE

\[ \ddot{q}(t) + c \dot{q}(t) + k q(t) = u(t), \]

with given constants \( c, k > 0 \). Define an appropriate state vector \( x(t) \) and write the system in first-order state-space form \( \dot{x} = f(x,u) \).

Solution. A natural choice of state is \( x_1(t) = q(t) \), \( x_2(t) = \dot{q}(t) \). Then

\[ \dot{x}_1(t) = x_2(t), \quad \dot{x}_2(t) = -c x_2(t) - k x_1(t) + u(t). \]

In vector form, with \( x(t) = \begin{bmatrix} x_1(t) \\ x_2(t) \end{bmatrix} \),

\[ \dot{x}(t) = \begin{bmatrix} 0 & 1 \\ -k & -c \end{bmatrix} x(t) + \begin{bmatrix} 0 \\ 1 \end{bmatrix} u(t). \]

This example shows how higher-order ODEs are systematically converted to first-order systems, a standard step in system dynamics and control.

Problem 4 (Uniqueness via Difference of Solutions). Suppose \( x_1(t) \) and \( x_2(t) \) are two solutions of the homogeneous linear system

\[ \dot{x}(t) = A x(t), \quad x(0) = x_0, \]

where \( A \in \mathbb{R}^{n \times n} \) is constant. Show that \( x_1(t) \equiv x_2(t) \) for all \( t \), using only basic linear algebra and the properties of derivatives.

Solution. Consider the difference \( z(t) = x_1(t) - x_2(t) \). Then

\[ \dot{z}(t) = \dot{x}_1(t) - \dot{x}_2(t) = A x_1(t) - A x_2(t) = A z(t), \quad z(0) = x_1(0) - x_2(0) = 0. \]

By linearity, the zero function \( z(t) \equiv 0 \) satisfies the same IVP. By the uniqueness theorem for linear systems (a special case of the Picard–Lindelöf result), the zero solution is the only solution of this IVP. Hence \( z(t) \equiv 0 \) and \( x_1(t) \equiv x_2(t) \). This argument is typical in system dynamics for proving that a model is well-posed.

Problem 5 (Modeling Workflow Overview). Sketch a conceptual workflow from a real-world engineering problem to a dynamic model suitable for analysis and control design.

Solution (conceptual flow chart). One possible workflow is summarized below:

flowchart TD
  ST["Problem description"] --> VAR["Identify variables and signals"]
  VAR --> LAW["Apply physical laws / conservation principles"]
  LAW --> EQN["Derive governing equations (ODEs)"]
  EQN --> MOD["Choose state variables and write state-space model"]
  MOD --> SIM["Select numerical methods / software for simulation"]
  SIM --> ANA["Analyze trajectories and system properties"]
  ANA --> REF["Refine model or design controller"]
        

In this lesson, we introduced system dynamics as the study of how systems evolve in time under the influence of inputs and initial conditions. We formalized a dynamic system as an operator mapping signals and initial states to trajectories, and adopted state-space models as our primary representation. We emphasized the role of the state in capturing the Markov property, reviewed existence and uniqueness results for ODEs, and used a simple first-order system to illustrate exponential behavior.

On the computational side, we implemented the same basic dynamic system in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, highlighting that the core mathematical structure is independent of the programming language or software. In the next lesson, we will refine the classification of systems into dynamic vs. static, linear vs. nonlinear, and time-invariant vs. time-varying, building directly on the foundational concepts introduced here.