Many engineered subsystems behave as (or are well-approximated by) a second-order model: mass–spring–damper motion, RLC electrical networks, actuator–load pairs, and lightly damped structural modes. A key reason second-order systems are foundational is that a large class of responses can be described by just two dimensionless parameters: \( \omega_n \) (time scale) and \( \zeta \) (dissipation relative to the critical case).

The lesson’s workflow is: start from a second-order model (ODE or transfer function), normalize it, extract \( \omega_n \) and \( \zeta \), and then interpret or simulate the response.

flowchart TD
  A["Start: 2nd-order model (ODE or G(s))"] --> B["Normalize coefficients (make leading term = 1)"]
  B --> C["Match: s^2 + 2*zeta*wn*s + wn^2"]
  C --> D["Extract parameters: wn, zeta, gain K"]
  D --> E["Choose standard form: unity DC gain or general K"]
  E --> F["Simulate time response (step/impulse)"]
  F --> G["Interpret time scale (wn) and damping (zeta)"]
        

Consider the (strictly proper) second-order LTI input–output model with constant coefficients:

\[ a_2 \frac{d^2 y(t)}{dt^2} + a_1 \frac{dy(t)}{dt} + a_0 y(t) = b_0 u(t), \quad a_2 \neq 0. \]

Assuming zero initial conditions (as in the transfer-function derivations of Chapter 6), applying the Laplace transform yields:

\[ a_2 s^2 Y(s) + a_1 s Y(s) + a_0 Y(s) = b_0 U(s). \]

Hence the transfer function is:

\[ G(s) \equiv \frac{Y(s)}{U(s)} = \frac{b_0}{a_2 s^2 + a_1 s + a_0}. \]

This form is not yet “interpretable” in a universal way. The objective is to rewrite the denominator so that the coefficients directly encode the natural time scale and the relative damping.

Divide the ODE by \( a_2 \) to obtain the monic form:

\[ \frac{d^2 y}{dt^2} + \alpha \frac{dy}{dt} + \beta y = \gamma u, \quad \alpha \equiv \frac{a_1}{a_2},\; \beta \equiv \frac{a_0}{a_2},\; \gamma \equiv \frac{b_0}{a_2}. \]

The corresponding transfer function becomes:

\[ G(s) = \frac{\gamma}{s^2 + \alpha s + \beta}. \]

The canonical (control-engineering) second-order denominator is written as:

\[ s^2 + 2\zeta \omega_n s + \omega_n^2. \]

Coefficient-to-parameter mapping (proof). Match coefficients term-by-term:

\[ \omega_n^2 = \beta, \quad 2\zeta \omega_n = \alpha. \]

Therefore:

\[ \omega_n = \sqrt{\beta}, \quad \zeta = \frac{\alpha}{2\sqrt{\beta}} = \frac{a_1}{2\sqrt{a_0 a_2}} \quad (\text{when } a_0, a_2 > 0). \]

This mapping is algebraic and unique for \( \beta > 0 \). In practice, physical models often ensure \( a_0, a_2 > 0 \), yielding \( \omega_n > 0 \). The sign and magnitude of \( \zeta \) then encode the damping relative to the critical case.

Standard transfer function with gain. Define the static (DC) gain: \( K \equiv G(0) \). From \( G(s)=\gamma/(s^2+\alpha s+\beta) \):

\[ K = G(0) = \frac{\gamma}{\beta}. \]

Substitute \( \gamma = K\beta = K\omega_n^2 \) to obtain the widely used standard form:

\[ G(s) = \frac{K\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2}. \]

This parameterization separates: (i) overall gain \( K \), (ii) time scaling \( \omega_n \), and (iii) damping shape parameter \( \zeta \).

The natural frequency is defined by the undamped homogeneous dynamics. Consider the homogeneous equation with no damping and no forcing:

\[ \frac{d^2 y}{dt^2} + \omega_n^2 y = 0. \]

Claim. The solutions are sinusoidal with angular frequency \( \omega_n \).

Proof. Substitute a trial solution \( y(t) = e^{st} \). Then:

\[ s^2 e^{st} + \omega_n^2 e^{st} = 0 \;\;\Rightarrow\;\; s^2 + \omega_n^2 = 0 \;\;\Rightarrow\;\; s = \pm j\omega_n. \]

Therefore, \( y(t) = c_1 \cos(\omega_n t) + c_2 \sin(\omega_n t) \), which oscillates at angular frequency \( \omega_n \).

In physical models, \( \omega_n \) is determined by inertia and stiffness. For example, a mass–spring system:

\[ m \ddot{x} + k x = 0 \quad \Rightarrow \quad \omega_n = \sqrt{\frac{k}{m}}. \]

Damping ratio is defined by comparing actual damping to the critical damping that makes the characteristic polynomial have a repeated real root. Start from the mass–spring–damper free dynamics:

\[ m \ddot{x} + c \dot{x} + k x = 0. \]

Substituting \( x(t)=e^{st} \) gives the characteristic equation:

\[ m s^2 + c s + k = 0. \]

The discriminant is \( \Delta = c^2 - 4mk \). The boundary between “oscillatory-like” and “non-oscillatory-like” behavior occurs at \( \Delta = 0 \), i.e., when the two roots coincide.

Critical damping coefficient (proof).

\[ \Delta = 0 \;\;\Rightarrow\;\; c^2 - 4mk = 0 \;\;\Rightarrow\;\; c_c = 2\sqrt{mk}. \]

The damping ratio is the normalized quantity:

\[ \zeta \equiv \frac{c}{c_c} = \frac{c}{2\sqrt{mk}}. \]

Also, dividing the ODE by \( m \) yields:

\[ \ddot{x} + \frac{c}{m}\dot{x} + \frac{k}{m}x = 0 \quad \Rightarrow \quad \ddot{x} + 2\zeta\omega_n \dot{x} + \omega_n^2 x = 0, \]

where \( \omega_n = \sqrt{k/m} \) and \( 2\zeta\omega_n = c/m \).

A major advantage of the standard form is that it cleanly separates time scaling from response “shape.” Consider the normalized second-order ODE with input:

\[ \ddot{y} + 2\zeta\omega_n \dot{y} + \omega_n^2 y = K\omega_n^2 u(t). \]

Define a dimensionless time variable \( \theta \equiv \omega_n t \). Then: \( \frac{d}{dt} = \omega_n \frac{d}{d\theta} \) and \( \frac{d^2}{dt^2} = \omega_n^2 \frac{d^2}{d\theta^2} \).

Derivation (proof). Substitute into the ODE:

\[ \omega_n^2 \frac{d^2 y}{d\theta^2} + 2\zeta\omega_n \left(\omega_n \frac{dy}{d\theta}\right) + \omega_n^2 y = K\omega_n^2 u\!\left(\frac{\theta}{\omega_n}\right). \]

Divide both sides by \( \omega_n^2 \):

\[ \frac{d^2 y}{d\theta^2} + 2\zeta \frac{dy}{d\theta} + y = K\, u\!\left(\frac{\theta}{\omega_n}\right). \]

For a step input (constant after application), the right side becomes constant after the step, and the response shape in \( \theta \)-time depends primarily on \( \zeta \) (with \( K \) setting the final value). Thus \( \omega_n \) scales time, while \( \zeta \) governs the qualitative transient form.

Starting from the standard second-order ODE:

\[ \ddot{y} + 2\zeta\omega_n \dot{y} + \omega_n^2 y = K\omega_n^2 u, \]

define states \( x_1 \equiv y \) and \( x_2 \equiv \dot{y} \). Then:

\[ \dot{x}_1 = x_2, \quad \dot{x}_2 = -\omega_n^2 x_1 - 2\zeta\omega_n x_2 + K\omega_n^2 u, \quad y = x_1. \]

In matrix form:

\[ \dot{\mathbf{x}} = \begin{bmatrix} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{bmatrix}\mathbf{x} + \begin{bmatrix} 0 \\ K\omega_n^2 \end{bmatrix}u, \quad y = \begin{bmatrix} 1 & 0 \end{bmatrix}\mathbf{x}. \]

Transfer function consistency (proof sketch). Using the standard state-space to transfer-function identity: \( G(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} \), one recovers:

\[ G(s) = \frac{K\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \]

matching the standard form derived earlier from the ODE.

flowchart TD
  ODE["ODE: y'' + 2*zeta*wn*y' + wn^2*y = K*wn^2*u"] <--> TF["G(s)=K*wn^2 / (s^2 + 2*zeta*wn*s + wn^2)"]
  TF <--> SS["State space: x1'=x2; x2'=-wn^2*x1-2*zeta*wn*x2+K*wn^2*u"]
  SS <--> ODE
        

We adopt a common example parameter set for all languages: \( \omega_n = 5 \) rad/s, \( \zeta = 0.3 \), and \( K = 2 \). The standard transfer function is:

\[ G(s) = \frac{K\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2} = \frac{50}{s^2 + 3s + 25}. \]

8.1 Python (python-control, SciPy) + From-Scratch RK4

Recommended libraries: control (python-control) for classical LTI tools, scipy.signal for signal processing LTI utilities, and numpy.

import numpy as np

# Parameters
wn = 5.0
zeta = 0.3
K = 2.0

# --- Library-based (python-control) ---
# pip install control
import control as ct

num = [K * wn**2]
den = [1.0, 2.0 * zeta * wn, wn**2]
G = ct.tf(num, den)

t = np.linspace(0.0, 5.0, 2000)
t_out, y_out = ct.step_response(G, T=t)

print("G(s) =", G)
print("Approx final value:", y_out[-1])

# --- From-scratch RK4 simulation of state space ---
# x1 = y, x2 = ydot
A = np.array([[0.0, 1.0],
              [-wn**2, -2.0 * zeta * wn]])
B = np.array([[0.0],
              [K * wn**2]])
C = np.array([[1.0, 0.0]])

def rk4_step(x, u, dt):
    def f(xv):
        return (A @ xv.reshape(2,1) + B * u).reshape(2,)
    k1 = f(x)
    k2 = f(x + 0.5 * dt * k1)
    k3 = f(x + 0.5 * dt * k2)
    k4 = f(x + dt * k3)
    return x + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4)

dt = 0.001
N = int(5.0 / dt)
x = np.array([0.0, 0.0])
y = np.zeros(N)
tt = np.arange(N) * dt

# Unit step input u(t)=1 for t>=0
for i in range(N):
    x = rk4_step(x, 1.0, dt)
    y[i] = (C @ x.reshape(2,1)).item()

print("RK4 final value:", y[-1])
      

8.2 C++ (Eigen) + From-Scratch RK4

Recommended library: Eigen for linear algebra. The simulation uses RK4 on the second-order state-space form. Note: in HTML code display, angle brackets are written as < and >.

#include <iostream>
#include <vector>
#include <Eigen/Dense>

using Eigen::Matrix2d;
using Eigen::Vector2d;

struct SecondOrder {
  double wn;
  double zeta;
  double K;
  Matrix2d A;
  Vector2d B;
  Eigen::RowVector2d C;

  SecondOrder(double wn_, double zeta_, double K_) : wn(wn_), zeta(zeta_), K(K_) {
    A << 0.0, 1.0,
         -wn*wn, -2.0*zeta*wn;
    B << 0.0, K*wn*wn;
    C << 1.0, 0.0;
  }

  Vector2d f(const Vector2d& x, double u) const {
    return A * x + B * u;
  }

  Vector2d rk4_step(const Vector2d& x, double u, double dt) const {
    Vector2d k1 = f(x, u);
    Vector2d k2 = f(x + 0.5*dt*k1, u);
    Vector2d k3 = f(x + 0.5*dt*k2, u);
    Vector2d k4 = f(x + dt*k3, u);
    return x + (dt/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);
  }
};

int main() {
  double wn = 5.0, zeta = 0.3, K = 2.0;
  SecondOrder sys(wn, zeta, K);

  double dt = 1e-3;
  int N = static_cast<int>(5.0 / dt);

  Vector2d x(0.0, 0.0);
  std::vector<double> y(N);

  for (int i = 0; i < N; ++i) {
    double u = 1.0; // unit step for t>=0
    x = sys.rk4_step(x, u, dt);
    y[i] = (sys.C * x)(0);
  }

  std::cout << "Approx final value: " << y.back() << std::endl;
  return 0;
}
      

8.3 Java (Apache Commons Math optional) + From-Scratch RK4

Recommended library: Apache Commons Math provides ODE solvers, but below is a minimal from-scratch RK4 to keep the numerical method transparent.

public class SecondOrderRK4 {
    // x1 = y, x2 = ydot
    static double wn = 5.0;
    static double zeta = 0.3;
    static double K = 2.0;

    static double[] f(double[] x, double u) {
        double x1dot = x[1];
        double x2dot = -wn*wn*x[0] - 2.0*zeta*wn*x[1] + K*wn*wn*u;
        return new double[]{x1dot, x2dot};
    }

    static double[] rk4Step(double[] x, double u, double dt) {
        double[] k1 = f(x, u);

        double[] x2 = new double[]{x[0] + 0.5*dt*k1[0], x[1] + 0.5*dt*k1[1]};
        double[] k2 = f(x2, u);

        double[] x3 = new double[]{x[0] + 0.5*dt*k2[0], x[1] + 0.5*dt*k2[1]};
        double[] k3 = f(x3, u);

        double[] x4 = new double[]{x[0] + dt*k3[0], x[1] + dt*k3[1]};
        double[] k4 = f(x4, u);

        return new double[]{
            x[0] + (dt/6.0)*(k1[0] + 2.0*k2[0] + 2.0*k3[0] + k4[0]),
            x[1] + (dt/6.0)*(k1[1] + 2.0*k2[1] + 2.0*k3[1] + k4[1])
        };
    }

    public static void main(String[] args) {
        double dt = 1e-3;
        int N = (int)(5.0 / dt);

        double[] x = new double[]{0.0, 0.0};
        double y = 0.0;

        for (int i = 0; i < N; i++) {
            double u = 1.0; // unit step
            x = rk4Step(x, u, dt);
            y = x[0];
        }
        System.out.println("Approx final value: " + y);
    }
}
      

8.4 MATLAB and Simulink (Control System Toolbox)

MATLAB’s Control System Toolbox directly supports transfer functions and step responses. Simulink provides block-level realization of the same standard form.

% Parameters
wn = 5;
zeta = 0.3;
K = 2;

% Transfer function G(s) = K*wn^2 / (s^2 + 2*zeta*wn*s + wn^2)
num = [K*wn^2];
den = [1 2*zeta*wn wn^2];
G = tf(num, den);

figure; step(G);
grid on;
title('Second-Order Step Response');

% State-space form with x1 = y, x2 = ydot
A = [0 1; -wn^2 -2*zeta*wn];
B = [0; K*wn^2];
C = [1 0];
D = 0;
sys_ss = ss(A,B,C,D);

figure; step(sys_ss);
grid on;
title('State-Space Step Response');
      

Simulink realization (conceptual checklist). Build a model with:

• Step block (Amplitude = 1)
• Transfer Fcn block with Numerator = [K*wn^2] and Denominator = [1 2*zeta*wn wn^2]
• Scope block to view output

8.5 Wolfram Mathematica (Symbolic + Numeric Response)

Mathematica can represent transfer functions symbolically and compute time responses.

wn = 5;
zeta = 0.3;
K = 2;

(* Transfer function model *)
G = TransferFunctionModel[{K*wn^2}, {s^2 + 2*zeta*wn*s + wn^2}, s];

(* Step response (unit step) *)
y[t_] = OutputResponse[G, UnitStep[t], t];

(* Plot over 0..5 seconds *)
Plot[y[t], {t, 0, 5}, PlotRange -> All, AxesLabel -> {"t", "y(t)"}]

(* Verify DC gain *)
dcGain = TransferFunctionExpand[G] /. s -> 0
      

Problem 1 (Coefficient Matching): Given \( 0.5\,\ddot{y} + 2\,\dot{y} + 8\,y = 16\,u(t) \), compute \( \omega_n \), \( \zeta \), and express \( G(s) \) in the standard form \( \frac{K\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2} \).

Solution: Divide by 0.5:

\[ \ddot{y} + 4\dot{y} + 16y = 32u(t). \]

Match coefficients:

\[ \omega_n^2 = 16 \Rightarrow \omega_n = 4, \quad 2\zeta\omega_n = 4 \Rightarrow \zeta = \frac{4}{2\cdot 4} = 0.5. \]

Compute \(K\) using \(K\omega_n^2 = 32\):

\[ K = \frac{32}{16} = 2. \]

Therefore:

\[ G(s) = \frac{32}{s^2 + 4s + 16} = \frac{2\cdot 4^2}{s^2 + 2(0.5)(4)s + 4^2}. \]

Problem 2 (Mass–Spring–Damper Parameters): For a mass–spring–damper system \( m\ddot{x} + c\dot{x} + kx = f(t) \), show that \( \omega_n = \sqrt{k/m} \) and \( \zeta = c/(2\sqrt{mk}) \).

Solution: Divide by \(m\):

\[ \ddot{x} + \frac{c}{m}\dot{x} + \frac{k}{m}x = \frac{1}{m}f(t). \]

Match with \( \ddot{x} + 2\zeta\omega_n \dot{x} + \omega_n^2 x = (\cdot)f(t) \):

\[ \omega_n^2 = \frac{k}{m} \Rightarrow \omega_n = \sqrt{\frac{k}{m}}, \quad 2\zeta\omega_n = \frac{c}{m} \Rightarrow \zeta = \frac{c}{2m\omega_n} = \frac{c}{2\sqrt{mk}}. \]

Problem 3 (Critical Damping Condition): For \( ms^2 + cs + k = 0 \), derive the critical damping value \( c_c = 2\sqrt{mk} \) by the repeated-root condition.

Solution: Repeated roots occur when the discriminant is zero:

\[ \Delta = c^2 - 4mk = 0 \Rightarrow c_c = 2\sqrt{mk}. \]

Problem 4 (Time-Scaling): Starting from \( \ddot{y} + 2\zeta\omega_n \dot{y} + \omega_n^2 y = K\omega_n^2 u(t) \), prove that with \( \theta=\omega_n t \), the equation becomes \( y''(\theta) + 2\zeta y'(\theta) + y(\theta) = K\,u(\theta/\omega_n) \), where primes denote derivatives with respect to \( \theta \).

Solution: Use \( \frac{d}{dt}=\omega_n\frac{d}{d\theta} \) and \( \frac{d^2}{dt^2}=\omega_n^2\frac{d^2}{d\theta^2} \). Substitute and divide by \( \omega_n^2 \) (see Section 6).

Problem 5 (State-Space Construction): For the standard second-order model, define \( x_1=y \), \( x_2=\dot{y} \). Derive the matrices \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) in \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \), \( y=\mathbf{C}\mathbf{x}+\mathbf{D}u \).

Solution: From Section 7:

\[ \mathbf{A}= \begin{bmatrix} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{bmatrix}, \quad \mathbf{B}= \begin{bmatrix} 0\\ K\omega_n^2 \end{bmatrix}, \quad \mathbf{C}= \begin{bmatrix} 1 & 0 \end{bmatrix}, \quad \mathbf{D}=0. \]

We established the second-order standard form by normalizing a general second-order LTI ODE and matching coefficients to define \( \omega_n \) and \( \zeta \). We proved the coefficient mapping, derived critical damping for the mass–spring–damper model, and showed how time-scaling separates response shape (largely set by \( \zeta \)) from time scale (set by \( \omega_n \)). We also connected ODE, transfer-function, and state-space forms and provided multi-language implementations. The next lesson will use these standard forms to derive and classify step responses.