In robot mechanics, the evolution of joint positions and velocities over time is governed by ordinary differential equations. An ordinary differential equation is a relation between an unknown function and a finite number of its derivatives with respect to a single independent variable, here time \( t \). The most common abstract form for mechanical systems is the first-order initial value problem (IVP)

\[ \dot{\mathbf{x}}(t) = \mathbf{f}(t,\mathbf{x}(t),\mathbf{u}(t)), \quad \mathbf{x}(t_0) = \mathbf{x}_0, \]

where \( \mathbf{x}(t) \in \mathbb{R}^n \) is the state (e.g., joint angles and velocities), \( \mathbf{u}(t) \) is an input (e.g., joint torque or motor command), and \( \mathbf{f} : \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n \) encodes the mechanics.

In simple mechanical systems, we often start from Newton's law or torque balance, which naturally yields second-order ODEs in generalized coordinates. We then convert them into first-order form compatible with state-space modeling (recall linear algebra and eigenvalue tools from previous lessons).

flowchart TD
  A["Physical description (links, masses, joints)"] --> B["Apply Newton / energy laws"]
  B --> C["Derive equations: m*qdd + b*qdot + k*q = u(t)"]
  C --> D["Rewrite as first-order system xdot = f(t,x,u)"]
  D --> E["Integrate over time to obtain trajectories"]
  E --> F["Use trajectories in kinematics / dynamics analysis"]
        

Consider a single revolute joint of a robot arm driven by a motor. Let \( q(t) \) be the joint angle, with inertia \( J > 0 \), viscous friction coefficient \( b \ge 0 \), and torsional stiffness \( k \ge 0 \). Let \( u(t) \) denote the applied motor torque. A standard torque balance gives

\[ J\ddot{q}(t) + b\dot{q}(t) + k q(t) = u(t). \]

This is a linear, time-invariant, second-order ODE in the scalar variable \( q(t) \). It is convenient to normalize the equation by dividing by \( J \) and introducing the natural frequency and damping ratio:

\[ \omega_n = \sqrt{\frac{k}{J}}, \quad \zeta = \frac{b}{2\sqrt{kJ}} \quad (k > 0). \]

Then the equation becomes

\[ \ddot{q}(t) + 2\zeta \omega_n \dot{q}(t) + \omega_n^2 q(t) = \frac{1}{J} u(t). \]

For a robot joint without stiffness (direct drive or negligible elasticity, \( k = 0 \)), we obtain a pure inertial and damping model

\[ J\ddot{q}(t) + b\dot{q}(t) = u(t), \]

which is still a second-order ODE but no longer has a restoring force.

To connect with the state-space viewpoint used later in dynamics, define state components \( x_1(t) = q(t) \), \( x_2(t) = \dot{q}(t) \). Then

\[ \begin{aligned} \dot{x}_1(t) &= x_2(t),\\ \dot{x}_2(t) &= -\frac{k}{J} x_1(t) - \frac{b}{J} x_2(t) + \frac{1}{J}u(t). \end{aligned} \]

In vector form with \( \mathbf{x}(t) = \begin{bmatrix} x_1(t) \\ x_2(t) \end{bmatrix} \) we can write

\[ \dot{\mathbf{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}u(t), \quad \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -\dfrac{k}{J} & -\dfrac{b}{J} \end{bmatrix},\quad \mathbf{B} = \begin{bmatrix} 0 \\ \dfrac{1}{J} \end{bmatrix}. \]

This is precisely the LTI state-space form that will be central in later dynamics and control courses.

Robot and mechanical models exhibit several important ODE classes:

• Order: number of time derivatives of the unknown. Single-DOF joint dynamics yield a second-order ODE in \( q \).
• Scalar vs. vector: for an \( n \)-DOF mechanical system, we obtain vector ODEs in \( \mathbf{q}(t) \in \mathbb{R}^n \) and \( \dot{\mathbf{q}}(t) \). A typical linear approximation is

  \[ \mathbf{M}\ddot{\mathbf{q}}(t) + \mathbf{B}\dot{\mathbf{q}}(t) + \mathbf{K}\mathbf{q}(t) = \mathbf{u}(t), \]

  where \( \mathbf{M} \) is symmetric positive definite, representing inertial coupling between joints.
• Linearity: an ODE is linear in the unknown if it can be written as

  \[ \mathbf{P}_0(t)\mathbf{x}(t) + \mathbf{P}_1(t)\dot{\mathbf{x}}(t) + \dots + \mathbf{P}_k(t)\mathbf{x}^{(k)}(t) = \mathbf{r}(t), \]

  with matrices \( \mathbf{P}_i(t) \) independent of \( \mathbf{x} \). The joint model with constant \( J,b,k \) is linear; the nonlinear pendulum equation

  \[ \ddot{\theta}(t) + \frac{g}{\ell}\sin\theta(t) = 0 \]

  is not, due to \( \sin\theta(t) \).
• Time-invariant vs. time-varying: if coefficients do not depend on \( t \), the system is time-invariant. Our joint model with constant \( J,b,k \) is time-invariant; a joint whose inertia effectively changes with a moving load is approximately time-varying.
• Homogeneous vs. forced: if the right-hand side \( \mathbf{r}(t) \) or input \( \mathbf{u}(t) \) is identically zero, the ODE is homogeneous; otherwise it is forced. For the joint, \( u(t)\equiv 0 \) describes free motion (e.g., passive joint after a disturbance).

Classification guides the choice of analytical solution methods and numerical integration schemes that will be used in later chapters.

A key theoretical question is: given a mechanical IVP, does there exist a trajectory, and is it unique? For state-space systems \( \dot{\mathbf{x}} = \mathbf{f}(t,\mathbf{x}) \), a standard result is the Picard–Lindelöf theorem.

Theorem (Picard–Lindelöf, informal).

Let \( \mathbf{f} : D \subset \mathbb{R}\times\mathbb{R}^n \to \mathbb{R}^n \) be continuous and locally Lipschitz in \( \mathbf{x} \) on a domain \( D \). That is, for every compact \( K \subset D \) there exists \( L \ge 0 \) such that

\[ \|\mathbf{f}(t,\mathbf{x}_1) - \mathbf{f}(t,\mathbf{x}_2)\| \le L \|\mathbf{x}_1 - \mathbf{x}_2\| \quad \text{for all } (t,\mathbf{x}_1),(t,\mathbf{x}_2) \in K. \]

Then for any initial condition \( (t_0,\mathbf{x}_0) \in D \) there exists a time interval \( [t_0 - \delta, t_0 + \delta] \) and a unique solution \( \mathbf{x}(t) \) of \( \dot{\mathbf{x}} = \mathbf{f}(t,\mathbf{x}) \) with \( \mathbf{x}(t_0) = \mathbf{x}_0 \).

For mechanical systems built from smooth force laws (spring, viscous damping, gravity, smooth actuator torques), the function \( \mathbf{f} \) is typically smooth, hence Lipschitz, and solutions are guaranteed to exist and be unique, at least locally. Physically: given the present configuration and velocity of a robot, plus the applied torques, the future trajectory is deterministic.

The classical proof constructs successive approximations via Picard iteration, showing that the integral operator

\[ (\mathcal{T}\mathbf{x})(t) = \mathbf{x}_0 + \int_{t_0}^{t} \mathbf{f}(\tau,\mathbf{x}(\tau))\,d\tau \]

is a contraction on a space of continuous functions, so it has a unique fixed point, which is the unique solution trajectory.

For the homogeneous joint model

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

we seek solutions of the form \( q(t) = e^{\lambda t} \). Substitution yields the characteristic equation

\[ \lambda^2 + 2\zeta \omega_n \lambda + \omega_n^2 = 0. \]

This quadratic has roots

\[ \lambda_{1,2} = -\zeta \omega_n \pm \omega_n \sqrt{\zeta^2 - 1}. \]

Three regimes arise:

• Underdamped: \( 0 < \zeta < 1 \). Roots are complex conjugates \( \lambda_{1,2} = -\zeta\omega_n \pm \mathrm{i}\omega_d \) with \( \omega_d = \omega_n\sqrt{1-\zeta^2} \). The solution is

  \[ q(t) = e^{-\zeta \omega_n t} \left( C_1 \cos(\omega_d t) + C_2 \sin(\omega_d t) \right), \]

  where \( C_1,C_2 \) are determined by \( q(0),\dot{q}(0) \).
• Critically damped: \( \zeta = 1 \). Double real root \( \lambda = -\omega_n \); solution

  \[ q(t) = (C_1 + C_2 t)e^{-\omega_n t}. \]

• Overdamped: \( \zeta > 1 \). Two distinct real roots \( \lambda_1,\lambda_2 < 0 \) and

  \[ q(t) = C_1 e^{\lambda_1 t} + C_2 e^{\lambda_2 t}. \]

For a forced system \( \ddot{q} + 2\zeta \omega_n \dot{q} + \omega_n^2 q = \frac{1}{J} u(t) \), the solution decomposes into a homogeneous part and a particular part:

\[ q(t) = q_h(t) + q_p(t), \]

where \( q_h(t) \) solves the homogeneous equation and \( q_p(t) \) is any specific solution to the full forced equation. For simple inputs (steps, sinusoids), one can find \( q_p(t) \) by the method of undetermined coefficients or via Laplace transforms; in robotics, such formulas are useful as benchmarks for numerical integration.

Using the state definition \( \mathbf{x}(t) = [q(t),\dot{q}(t)]^\top \), we have the linear time-invariant (LTI) state equation

\[ \dot{\mathbf{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}u(t), \]

where \( \mathbf{A},\mathbf{B} \) were derived in Section 2. For the homogeneous case \( u(t)\equiv 0 \), the solution is

\[ \mathbf{x}(t) = e^{\mathbf{A}(t-t_0)}\mathbf{x}_0, \]

where \( e^{\mathbf{A}t} \) is the matrix exponential, defined via its Taylor series

\[ e^{\mathbf{A}t} = \mathbf{I} + \mathbf{A}t + \frac{(\mathbf{A}t)^2}{2!} + \frac{(\mathbf{A}t)^3}{3!} + \dots. \]

If \( \mathbf{A} \) is diagonalizable, \( \mathbf{A} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} \), then

\[ e^{\mathbf{A}t} = \mathbf{V}e^{\boldsymbol{\Lambda} t}\mathbf{V}^{-1}, \quad e^{\boldsymbol{\Lambda} t} = \operatorname{diag}(e^{\lambda_1 t},\dots,e^{\lambda_n t}), \]

which directly connects eigenvalues of \( \mathbf{A} \) (from Lesson 2) to the modes of motion of the mechanical system.

For nonzero input \( u(t) \), the solution of the LTI system is

\[ \mathbf{x}(t) = e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 + \int_{t_0}^{t} e^{\mathbf{A}(t-\tau)} \mathbf{B} u(\tau)\,d\tau, \]

which forms the basis for many simulation and control algorithms in robotics.

Exact formulas are available only for special ODEs. For realistic robot models, we rely on numerical integration to approximate trajectories. A simple explicit Euler step for \( \dot{\mathbf{x}} = \mathbf{f}(t,\mathbf{x}) \) with step size \( h > 0 \) reads

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + h\,\mathbf{f}(t_k,\mathbf{x}_k), \quad t_{k+1} = t_k + h. \]

Higher-order methods (Runge–Kutta, multi-step) will be treated in the next lesson on numerical computation, but it is useful already to view ODE simulation as repeated application of an integration step.

flowchart TD
  I["Initial time t0, state x0, step h"] --> K["k = 0"]
  K --> F["Evaluate xdot = f(tk, xk)"]
  F --> U["Update x_{k+1} = x_k + h * xdot"]
  U --> T["Update tk+1 = tk + h"]
  T --> L["Check time horizon"]
  L -->|continue| F
  L -->|stop| O["Return discrete trajectory {tk, xk}"]
        

We now implement the single-DOF joint model \( J\ddot{q} + b\dot{q} + k q = u(t) \) in several programming environments commonly used in robotics. We use:

• Python with scipy.integrate (and optionally roboticstoolbox-python in later labs).
• C++ with Eigen for linear algebra and Boost.Numeric.Odeint for integration.
• Java with Apache Commons Math ODE solvers.
• MATLAB/Simulink with ode45 and a basic block diagram.
• Wolfram Mathematica with NDSolve.

8.1 Python: scipy.integrate.solve_ivp and Explicit Euler

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

# Joint parameters
J = 0.01   # inertia
b = 0.1    # viscous damping
k = 1.0    # stiffness

def joint_ode(t, x):
    """
    x[0] = q (angle)
    x[1] = qd (angular velocity)
    """
    q, qd = x
    # Constant step torque input
    u = 1.0
    qdd = (u - b * qd - k * q) / J
    return [qd, qdd]

# High-level integrator (Runge-Kutta) via solve_ivp
t0, tf = 0.0, 5.0
x0 = [0.0, 0.0]  # q(0) = 0, qd(0) = 0

sol = solve_ivp(joint_ode, (t0, tf), x0, max_step=0.01, dense_output=True)

t_grid = np.linspace(t0, tf, 500)
x_grid = sol.sol(t_grid)

plt.figure()
plt.plot(t_grid, x_grid[0], label="q(t)")
plt.plot(t_grid, x_grid[1], label="qd(t)")
plt.xlabel("t [s]")
plt.legend()
plt.title("Single-DOF joint response (Python solve_ivp)")
plt.show()

# From-scratch explicit Euler integrator (for concept illustration)
def euler_integrate(f, t0, tf, x0, h):
    t_values = [t0]
    x_values = [np.array(x0, dtype=float)]
    t = t0
    x = np.array(x0, dtype=float)
    while t < tf:
        dx = np.array(f(t, x))
        x = x + h * dx
        t = t + h
        t_values.append(t)
        x_values.append(x.copy())
    return np.array(t_values), np.vstack(x_values)

t_e, x_e = euler_integrate(joint_ode, 0.0, 5.0, x0, h=0.001)
      

8.2 C++: Eigen + Boost.Odeint

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

using state_type = std::array<double, 2>;

// Parameters
const double J = 0.01;
const double b = 0.1;
const double k = 1.0;

// ODE: x[0] = q, x[1] = qd
void joint_ode(const state_type &x, state_type &dxdt, double t)
{
    (void)t; // unused
    double q  = x[0];
    double qd = x[1];
    double u  = 1.0; // constant torque
    double qdd = (u - b * qd - k * q) / J;
    dxdt[0] = qd;
    dxdt[1] = qdd;
}

int main()
{
    using namespace boost::numeric::odeint;

    state_type x = {0.0, 0.0}; // initial state

    auto observer = [](const state_type &x, double t) {
        std::cout << t << " " << x[0] << " " << x[1] << std::endl;
    };

    size_t steps = integrate(joint_ode, x, 0.0, 5.0, 0.01, observer);
    std::cerr << "Total steps: " << steps << std::endl;

    return 0;
}
      

8.3 Java: Apache Commons Math ODE Solvers

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

public class JointODEExample {

    // Parameters
    static final double J = 0.01;
    static final double b = 0.1;
    static final double k = 1.0;

    public static class JointODE implements FirstOrderDifferentialEquations {

        @Override
        public int getDimension() {
            return 2; // q, qd
        }

        @Override
        public void computeDerivatives(double t, double[] x, double[] dxdt) {
            double q  = x[0];
            double qd = x[1];
            double u  = 1.0; // constant torque
            double qdd = (u - b * qd - k * q) / J;
            dxdt[0] = qd;
            dxdt[1] = qdd;
        }
    }

    public static void main(String[] args) {
        FirstOrderDifferentialEquations ode = new JointODE();

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

        double minStep = 1.0e-3;
        double maxStep = 0.05;
        double absTol = 1.0e-8;
        double relTol = 1.0e-8;

        FirstOrderIntegrator integrator =
            new DormandPrince54Integrator(minStep, maxStep, absTol, relTol);

        integrator.addStepHandler((interpolator, isLast) -> {
            double t = interpolator.getCurrentTime();
            double[] x = interpolator.getInterpolatedState();
            System.out.println(t + " " + x[0] + " " + x[1]);
        });

        integrator.integrate(ode, t0, x0, tf, x0);
    }
}
      

8.4 MATLAB / Simulink

% Joint parameters
J = 0.01;
b = 0.1;
k = 1.0;

% ODE as first-order system: x(1) = q, x(2) = qd
joint_ode = @(t,x) [ ...
    x(2); ...
    (1.0 - b * x(2) - k * x(1)) / J ...
];

tspan = [0 5];
x0 = [0; 0];  % q(0) = 0, qd(0) = 0

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

figure;
plot(t, x(:,1), 'LineWidth', 1.5); hold on;
plot(t, x(:,2), 'LineWidth', 1.5);
xlabel('t [s]');
legend('q(t)', 'qd(t)');
title('Single-DOF joint response (MATLAB ode45)');
grid on;

% Minimal Simulink-style script (conceptual)
% 1) create a new model
% new_system('joint_model');
% open_system('joint_model');
%
% 2) add blocks (Integrator, Gain, Sum, Step, Scope) and connect them
% programmatically using add_block / add_line, or use the Simulink GUI.
% 3) parameterize gains with J, b, k and run:
% sim('joint_model');
      

8.5 Wolfram Mathematica

J = 0.01;
b = 0.1;
k = 1.0;

u[t_] := 1.0; (* constant torque *)

eq = J*q''[t] + b*q'[t] + k*q[t] == u[t];
ic = {q[0] == 0, q'[0] == 0};

sol = NDSolve[{eq, ic}, q, {t, 0, 5}];

Plot[Evaluate[q[t] /. sol], {t, 0, 5},
     AxesLabel -> {"t", "q(t)"},
     PlotLabel -> "Single-DOF joint response (Mathematica)"]
      

In this lesson we connected physical reasoning about robot joints and mechanical elements to ordinary differential equations. Starting from Newtonian torque balance, we derived second-order ODEs, normalized them in terms of natural frequency and damping ratio, and converted them into first-order state-space form. We discussed classification of ODEs (linear vs. nonlinear, time-invariant vs. time-varying, homogeneous vs. forced) and stated the existence and uniqueness theorem guaranteeing well-defined mechanical trajectories. Analytical solutions of second-order linear systems revealed how eigenvalues determine decay and oscillation, while numerical integration schemes such as Euler and Runge–Kutta provide practical tools for simulation. Finally, we implemented these concepts across Python, C++, Java, MATLAB/Simulink, and Mathematica, preparing the ground for the next lesson on numerical computation practices for robotics.