In control engineering, the dynamics of a single-input single-output (SISO) system are typically described by ordinary differential equations relating an output variable (e.g., position, angle, current) to input signals (e.g., torque, voltage). A general n-th order ODE for a scalar state \( y(t) \) has the form

\[ F\!\big(t, y(t), y'(t), \dots, y^{(n)}(t)\big) = 0, \quad t \in I \subset \mathbb{R}. \]

For mechanical and electrical systems, Newton's and Kirchhoff's laws give ODEs where the highest derivative captures the inertial or storage effects. For instance, the standard mass–spring–damper system with position \( q(t) \), mass \( m > 0 \), viscous damping \( b \ge 0 \), stiffness \( k \ge 0 \), and external force \( u(t) \) satisfies

\[ m\,\ddot q(t) + b\,\dot q(t) + k\,q(t) = u(t). \]

This equation is the starting point for many control problems, including vibration control and robot joint position control. It is a linear, second-order ODE with constant coefficients and will serve as our running example throughout the lesson.

flowchart TD
  P["Physical system (robot link, motor, etc.)"]
  L["Physical laws (Newton, Kirchhoff)"]
  M["Mathematical model: ODE"]
  A["Analytical properties (order, linearity, stability)"]
  N["Numerical simulation (Euler / Runge-Kutta)"]
  C["Later: controller synthesis"]
  P --> L
  L --> M
  M --> A
  M --> N
  A --> C
  N --> C
        

We recall some standard classifications, specialized to dynamical systems encountered in control engineering:

• Order: The highest derivative of \( y(t) \) appearing in the equation. For example, \( m\ddot q + b\dot q + kq = u(t) \) is second order.
• Autonomous vs. non-autonomous: An ODE is autonomous if it does not depend explicitly on time \( t \). For example

\[ \dot x(t) = f\big(x(t)\big) \]

is autonomous, whereas

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

is non-autonomous. Most regulation problems with constant setpoints lead to autonomous models once reference tracking is encoded in the state.

• Linear vs. nonlinear: A scalar ODE is linear in \( y \) if it can be written as

\[ a_n(t)\,y^{(n)}(t) + a_{n-1}(t)\,y^{(n-1)}(t) + \cdots + a_0(t)\,y(t) = b(t), \]

where the coefficients \( a_i(t) \) and forcing \( b(t) \) are known functions of \( t \) only. No products like \( y\,\dot y \) or nonlinearities such as \( \sin(y) \) may appear on the left-hand side.

• Time-invariant vs. time-varying: In control, we are particularly interested in linear time-invariant (LTI) systems, where the coefficients are constants:

\[ a_n\,y^{(n)}(t) + a_{n-1}\,y^{(n-1)}(t) + \cdots + a_0\,y(t) = b(t), \quad a_i \in \mathbb{R}. \]

Time invariance underpins the superposition principle and allows systematic analysis via exponentials and (later) Laplace transforms.

• Homogeneous vs. inhomogeneous (forced): The associated homogeneous equation sets \( b(t) = 0 \):

\[ a_n\,y^{(n)}(t) + \cdots + a_0\,y(t) = 0. \]

For linear systems, the general solution is the sum of a homogeneous (free) response and a particular (forced) response:

\[ y(t) = y_h(t) + y_p(t). \]

Many mechanical and electro-mechanical systems approximate to an LTI ODE around an operating point. Two canonical forms are:

(a) First-order system

\[ \dot x(t) + \frac{1}{\tau}\,x(t) = \frac{1}{\tau}\,u(t), \]

where \( x(t) \) is a state (e.g., capacitor voltage, fluid level), \( u(t) \) is an input, and \( \tau > 0 \) is the time constant determining how quickly the state reacts to changes in the input.

(b) Second-order system

The mass–spring–damper example can be written in the normalized form

\[ \ddot q(t) + 2\zeta \omega_n \dot q(t) + \omega_n^2 q(t) = \frac{\omega_n^2}{k}\,u(t), \]

where \( \omega_n = \sqrt{k/m} \) is the natural frequency and \( \zeta = \tfrac{b}{2\sqrt{mk}} \) is the damping ratio. Later chapters will interpret \( \omega_n \) and \( \zeta \) in terms of time response metrics such as rise time and overshoot; here we only use this normalization to simplify the structure of solutions.

Consider the homogeneous LTI ODE

\[ a_n y^{(n)}(t) + a_{n-1} y^{(n-1)}(t) + \cdots + a_0 y(t) = 0, \quad a_n \ne 0. \]

The standard ansatz is an exponential solution \( y(t) = e^{\lambda t} \). Substituting into the ODE yields

\[ a_n \lambda^n e^{\lambda t} + a_{n-1} \lambda^{n-1} e^{\lambda t} + \cdots + a_0 e^{\lambda t} = 0 \quad \Rightarrow \quad p(\lambda) = a_n \lambda^n + a_{n-1} \lambda^{n-1} + \cdots + a_0 = 0. \]

The polynomial \( p(\lambda) \) is the characteristic polynomial. Its roots \( \lambda_1, \dots, \lambda_n \) (counted with multiplicity) determine the homogeneous solution:

• If all roots are distinct and real, the solution is \( y_h(t) = \sum_{i=1}^n C_i e^{\lambda_i t} \).
• If a root \( \lambda \) has multiplicity \( m \), the contribution is \( e^{\lambda t}(C_0 + C_1 t + \cdots + C_{m-1} t^{m-1}) \).
• Complex conjugate pairs \( \lambda_{1,2} = \alpha \pm j\beta \), with \( \beta \ne 0 \), give oscillatory terms

\[ e^{\alpha t}\big(C_1 \cos(\beta t) + C_2 \sin(\beta t)\big). \]

In control, these exponentials correspond to modes of the system. Modes with negative real part decay, modes with zero real part persist, and modes with positive real part grow unbounded. This connection between characteristic roots and qualitative behavior will be central to stability analysis in later chapters.

For an inhomogeneous linear ODE

\[ a_n y^{(n)}(t) + \cdots + a_0 y(t) = b(t), \]

the general solution has the structure \( y(t) = y_h(t) + y_p(t) \), where \( y_h(t) \) solves the homogeneous equation with \( b(t) = 0 \), and \( y_p(t) \) is any particular solution of the full equation. Linearity implies that

\[ \mathcal{L}\big[y_h\big] = 0,\quad \mathcal{L}\big[y_p\big] = b(t) \quad \Rightarrow \quad \mathcal{L}\big[y_h + y_p\big] = b(t), \]

where \( \mathcal{L} \) denotes the linear differential operator \( \mathcal{L}[y] = a_n y^{(n)} + \dots + a_0 y \).

Classical methods such as undetermined coefficients and variation of parameters provide explicit constructions of \( y_p \) for simple input signals \( b(t) \) (e.g., steps, ramps, sinusoids). In control design, we will frequently use input types (step, ramp, etc.) as standard tests of performance; having explicit forms for \( y_p(t) \) will simplify analysis in those contexts.

High-order ODEs are often converted to an equivalent system of first-order equations. This representation is essential for both numerical integration and later controller design.

Example: mass–spring–damper. Define the state

\[ x_1(t) = q(t), \quad x_2(t) = \dot q(t), \quad x(t) = \begin{bmatrix} x_1(t) \\ x_2(t) \end{bmatrix}. \]

Then the second-order equation \( m\ddot q + b\dot q + kq = u(t) \) becomes the first-order system

\[ \dot x(t) = \begin{bmatrix} 0 & 1 \\ -\tfrac{k}{m} & -\tfrac{b}{m} \end{bmatrix} x(t) + \begin{bmatrix} 0 \\ \tfrac{1}{m} \end{bmatrix} u(t), \quad x(0) = x_0. \]

More generally, a linear time-invariant first-order system with input \( u(t) \) and state \( x(t) \in \mathbb{R}^n \) can be written as

\[ \dot x(t) = A x(t) + B u(t), \quad x(0) = x_0, \]

where \( A \in \mathbb{R}^{n \times n} \) and \( B \in \mathbb{R}^{n \times 1} \) are constant matrices. The solution can be expressed using the matrix exponential

\[ e^{At} = \sum_{k=0}^{\infty} \frac{(At)^k}{k!}, \]

giving

\[ x(t) = e^{At} x_0 + \int_0^t e^{A(t-\sigma)} B u(\sigma)\, d\sigma. \]

This formula generalizes the scalar exponential solution and highlights how eigenvalues of \( A \) control the qualitative behavior of the system. In later chapters, this representation will connect directly to transfer functions and stability criteria.

We recall a central result for first-order systems of ODEs, which justifies treating dynamical systems as deterministic maps from initial conditions and inputs to trajectories.

(Picard–Lindelöf type theorem, informal.)

Consider an initial-value problem (IVP)

\[ \dot x(t) = f(t,x(t)), \quad x(t_0) = x_0, \]

where \( f : D \subset \mathbb{R} \times \mathbb{R}^n \to \mathbb{R}^n \). Assume that in a neighborhood of \( (t_0, x_0) \):

• \( f \) is continuous, and
• there exists \( L > 0 \) such that \( \|f(t,x_1) - f(t,x_2)\| \le L \|x_1 - x_2\| \) for all \( x_1, x_2 \) in that neighborhood (Lipschitz continuity in \( x \)).

Then there exists a unique solution \( x(t) \) defined on some interval containing \( t_0 \). In particular, for linear systems \( \dot x = A x + B u(t) \) with continuous \( u(t) \), these conditions are automatically satisfied, so solutions exist and are unique for all \( t \).

This guarantees that, for well-posed dynamical models used in control, each initial state and input signal leads to a unique trajectory, a key requirement when we later talk about closed-loop behavior.

Exact analytical solutions are not always available or convenient. In practical control engineering, especially for robotics, numerical methods approximate trajectories of the IVP

\[ \dot x(t) = f(t,x(t)), \quad x(t_0) = x_0. \]

Let \( t_k = t_0 + k h \) be a time grid with step size \( h > 0 \). Two classical methods are:

(a) Explicit Euler method)

\[ x_{k+1} = x_k + h\,f(t_k, x_k). \]

Euler's method is first-order accurate (local truncation error \( \mathcal{O}(h^2) \), global error \( \mathcal{O}(h) \)) but may be unstable for stiff systems or large step sizes.

(b) Classical 4th-order Runge–Kutta (RK4)

\[ \begin{aligned} k_1 &= f(t_k, x_k), \\ k_2 &= f\big(t_k + \tfrac{h}{2}, x_k + \tfrac{h}{2} k_1\big), \\ k_3 &= f\big(t_k + \tfrac{h}{2}, x_k + \tfrac{h}{2} k_2\big), \\ k_4 &= f(t_k + h, x_k + h k_3), \\ x_{k+1} &= x_k + \tfrac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4). \end{aligned} \]

RK4 is fourth-order accurate and widely used for offline simulation of robot dynamics. Many robotics libraries wrap such schemes in high-level APIs for simulating manipulators and mobile platforms.

flowchart TD
  S["Initial state x0 at t0"] --> STEP["Choose step h"]
  STEP --> LOOP["For each step: compute slope(s) from f(t,x)"]
  LOOP --> UPDATE["Update state approximation x_{k+1}"]
  UPDATE --> CHECK["Reached final time?"]
  CHECK -->|no| LOOP
  CHECK -->|yes| OUT["Store trajectory x(t_k) for analysis or control design"]
        

We now simulate the mass–spring–damper system \( m\ddot q + b\dot q + kq = u(t) \) under a unit-step input \( u(t) = 1 \) with zero initial conditions. This example mirrors a single revolute joint of a robot arm driven by a torque step.

The state-space form is

\[ \dot x_1 = x_2, \quad \dot x_2 = \frac{1}{m}\big(u(t) - b x_2 - k x_1\big), \quad x_1(0) = 0, \; x_2(0) = 0. \]

9.1 Python — scipy.integrate.solve_ivp and robotics context

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Physical parameters (e.g., simplified robot joint)
m, b, k = 1.0, 0.2, 1.0  # mass, damping, stiffness

def dynamics(t, x):
    q, qdot = x
    u = 1.0  # unit step input
    qddot = (u - b * qdot - k * q) / m
    return [qdot, qddot]

t_span = (0.0, 10.0)
x0 = [0.0, 0.0]

sol = solve_ivp(dynamics, t_span, x0, max_step=0.01, dense_output=True)

t_plot = np.linspace(0.0, 10.0, 1000)
q_plot = sol.sol(t_plot)[0]

plt.plot(t_plot, q_plot)
plt.xlabel("t (s)")
plt.ylabel("q(t)")
plt.title("Mass-spring-damper step response")
plt.grid(True)
plt.show()

# In robotics, python-control and roboticstoolbox-python use similar ODE models
# to simulate joint and manipulator dynamics before controller design.
      

9.2 C++ — Explicit Euler Integration (with Robotics Notes)

#include <iostream>

int main() {
    double m = 1.0, b = 0.2, k = 1.0;
    double dt = 0.001;
    double q = 0.0;      // position
    double qdot = 0.0;   // velocity

    for (int i = 0; i < 10000; ++i) {
        double t = i * dt;
        double u = 1.0;  // unit step
        double qddot = (u - b * qdot - k * q) / m;

        // Explicit Euler step
        q     += dt * qdot;
        qdot  += dt * qddot;

        if (i % 1000 == 0) {
            std::cout << t << " " << q << std::endl;
        }
    }

    return 0;
}

// In robotics C++ stacks, Eigen is often used for vector/matrix operations,
// while ROS control frameworks integrate similar ODEs for joint-level dynamics.
      

9.3 Java — Simple Euler Integrator

public class MassSpringDamperEuler {
    public static void main(String[] args) {
        double m = 1.0, b = 0.2, k = 1.0;
        double dt = 0.001;
        double q = 0.0;
        double qdot = 0.0;

        for (int i = 0; i < 10000; i++) {
            double t = i * dt;
            double u = 1.0;
            double qddot = (u - b * qdot - k * q) / m;

            q    += dt * qdot;
            qdot += dt * qddot;

            if (i % 1000 == 0) {
                System.out.println(t + " " + q);
            }
        }
    }
}

// Java robotics frameworks (e.g., WPILib in mobile robotics) rely on
// numerically integrating ODEs representing drivetrain or arm dynamics.
      

9.4 MATLAB/Simulink — ode45 and Block-Diagram Simulation

m = 1.0; b = 0.2; k = 1.0;

f = @(t, x) [x(2);
             (1.0 - b * x(2) - k * x(1)) / m];

[t, x] = ode45(f, [0 10], [0; 0]);  % [position; velocity]

plot(t, x(:, 1));
xlabel('t (s)');
ylabel('q(t)');
title('Mass-spring-damper step response');
grid on;

% In Simulink, the same dynamics are implemented with two Integrator blocks
% in series (for q and qdot) plus Sum and Gain blocks for the forces.
% Robotics System Toolbox represents manipulator dynamics via similar ODEs.
      

9.5 Wolfram Mathematica — NDSolve for Symbolic-Numeric Workflows

m = 1.0; b = 0.2; k = 1.0;

sol = NDSolve[{
    m * q''[t] + b * q'[t] + k * q[t] == 1.0,
    q[0] == 0, q'[0] == 0
    },
    q, {t, 0, 10}
];

Plot[Evaluate[q[t] /. sol], {t, 0, 10},
     AxesLabel -> {"t", "q(t)"},
     PlotLabel -> "Mass-spring-damper step response"];

(* Mathematica is often used to derive closed-form solutions before exporting
   simplified ODEs for numerical simulation in robotics toolchains. *)
      

For the homogeneous system \( \dot x(t) = A x(t) \), a fundamental matrix solution is a matrix-valued function \( \Phi(t) \in \mathbb{R}^{n \times n} \) satisfying

\[ \dot\Phi(t) = A \Phi(t), \quad \Phi(0) = I. \]

It follows that \( \Phi(t) = e^{At} \) is a fundamental matrix. For any initial state \( x_0 \), the corresponding solution is \( x(t) = \Phi(t) x_0 \). If \( A \) is diagonalizable with \( A = V \Lambda V^{-1} \), then

\[ e^{At} = V e^{\Lambda t} V^{-1}, \quad e^{\Lambda t} = \mathrm{diag}(e^{\lambda_1 t},\dots,e^{\lambda_n t}), \]

and each eigenvalue \( \lambda_i \) contributes a mode \( e^{\lambda_i t} \) to the system response. This modal representation provides the bridge between algebraic properties of \( A \) and time-domain behavior of the dynamical system, which we will exploit in later discussions of stability and performance.

In this lesson we reviewed ordinary differential equations as the core modeling framework for dynamical systems in linear control. We recalled the classification of ODEs (order, linearity, time-variance, homogeneity), analyzed homogeneous and forced responses of linear time-invariant systems, and introduced the first-order vector form \( \dot x = A x + B u \) and its matrix-exponential solution. We also revisited existence and uniqueness conditions and surveyed basic numerical methods such as Euler and Runge–Kutta schemes, illustrating their use through multi-language simulations of a mass–spring–damper system relevant to robotic joint dynamics. These concepts form the mathematical foundation for the Laplace-transform tools and transfer function analysis developed in the remaining lessons of this chapter.