From Lesson 2, we use the standard second-order closed-loop form with unity DC gain: \( G(s)=\dfrac{\omega_n^2}{s^2+2\zeta\omega_n s+\omega_n^2} \), where \( \omega_n \) is the natural frequency and \( \zeta \) is the damping ratio.

For a unit-step input \( u(t)=1(t) \) (so \( U(s)=\dfrac{1}{s} \)), the output transform is \( Y(s)=G(s)U(s) \):

\[ Y(s)=\frac{\omega_n^2}{s\left(s^2+2\zeta\omega_n s+\omega_n^2\right)}. \]

In the time domain, the corresponding differential equation (zero initial conditions) is:

\[ \ddot{y}(t)+2\zeta\omega_n\dot{y}(t)+\omega_n^2 y(t)=\omega_n^2 u(t), \qquad y(0)=0,\;\dot{y}(0)=0. \]

A key simplification for step inputs is to define the tracking error \( e(t)=y(t)-1 \). Since a constant particular solution is \( y_p(t)=1 \), the error satisfies the homogeneous ODE:

\[ \ddot{e}(t)+2\zeta\omega_n\dot{e}(t)+\omega_n^2 e(t)=0, \qquad e(0)=-1,\;\dot{e}(0)=0. \]

Therefore, the classification of the step response is entirely determined by the roots of the characteristic polynomial \( s^2+2\zeta\omega_n s+\omega_n^2=0 \).

The characteristic equation is:

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

Its roots are:

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

The discriminant determines the qualitative response:

\[ \Delta=(2\zeta\omega_n)^2-4\omega_n^2=4\omega_n^2(\zeta^2-1). \]

flowchart TD
  S["Start: s^2 + 2*zeta*wn*s + wn^2 = 0"] --> D["Compute Delta = 4*wn^2*(zeta^2 - 1)"]
  D -->|"zeta < 1 (Delta < 0)"| U["Underdamped: \ncomplex poles"]
  D -->|"zeta = 1 (Delta = 0)"| C["Critically damped: \nrepeated real pole"]
  D -->|"zeta > 1 (Delta > 0)"| O["Overdamped: \ndistinct real poles"]
  U --> UF["y(t)=1 - exp(-zeta*wn*t)*(cos(wd*t) + \n(zeta/sqrt(1-zeta^2))*sin(wd*t))"]
  C --> CF["y(t)=1 - exp(-wn*t)*(1 + wn*t)"]
  O --> OF["y(t)=1 + (s2*exp(s1*t) - \ns1*exp(s2*t))/(s1 - s2)"]
  

In the underdamped case, we also define the damped natural frequency: \( \omega_d=\omega_n\sqrt{1-\zeta^2} \) (valid for \( 0 < \zeta < 1 \)).

For \( 0 < \zeta < 1 \), the roots are complex: \( s_{1,2}=-\zeta\omega_n \pm j\omega_d \), where \( \omega_d=\omega_n\sqrt{1-\zeta^2} \). The general homogeneous error solution is:

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

Apply the step-response initial conditions \( e(0)=-1 \), \( \dot{e}(0)=0 \). First, from \( e(0)=C_1=-1 \). Differentiate:

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

Evaluate at \( t=0 \):

\[ \dot{e}(0)= -\zeta\omega_n C_1 + \omega_d C_2 = 0 \;\; → \;\; C_2=\frac{\zeta\omega_n}{\omega_d}C_1 =-\frac{\zeta\omega_n}{\omega_d}. \]

Substitute \( C_1=-1 \), \( C_2=-\dfrac{\zeta\omega_n}{\omega_d} \), and use \( \dfrac{\omega_n}{\omega_d}=\dfrac{1}{\sqrt{1-\zeta^2}} \):

\[ e(t)=-e^{-\zeta\omega_n t}\left(\cos(\omega_d t)+\frac{\zeta}{\sqrt{1-\zeta^2}}\sin(\omega_d t)\right). \]

Since \( y(t)=1+e(t) \), the underdamped unit-step response is:

\[ y(t)=1-e^{-\zeta\omega_n t}\left(\cos(\omega_d t)+\frac{\zeta}{\sqrt{1-\zeta^2}}\sin(\omega_d t)\right), \qquad \omega_d=\omega_n\sqrt{1-\zeta^2}. \]

Oscillation and overshoot (derived from the closed form). The exponential envelope \( e^{-\zeta\omega_n t} \) guarantees decay, while the sinusoidal terms create oscillation. The first peak occurs when \( \dot{y}(t)=0 \); for the canonical step response this yields:

\[ t_p=\frac{\pi}{\omega_d}, \qquad M_p=\exp\!\left(\frac{-\pi\zeta}{\sqrt{1-\zeta^2}}\right), \qquad (0 < \zeta < 1). \]

These expressions will be used more systematically in Lesson 5 (performance metrics), but the key lesson here is structural: \( 0 < \zeta < 1 \) implies decaying oscillations because the poles are complex conjugates.

For \( \zeta=1 \), the characteristic equation has a repeated root: \( s=-\omega_n \) with multiplicity 2. The homogeneous error solution becomes:

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

Apply the initial conditions. From \( e(0)=C_1=-1 \). Differentiate:

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

At \( t=0 \):

\[ \dot{e}(0)=C_2-\omega_n C_1=0 \;\; → \;\; C_2=\omega_n C_1=-\omega_n. \]

Therefore:

\[ e(t)=(-1-\omega_n t)e^{-\omega_n t}, \qquad y(t)=1-(1+\omega_n t)e^{-\omega_n t}. \]

Proof of monotonicity and no overshoot. For all \( t > 0 \), we have \( (1+\omega_n t)e^{-\omega_n t} > 0 \), hence: \( y(t) < 1 \). Also,

\[ \dot{y}(t)=\frac{d}{dt}\left(1-(1+\omega_n t)e^{-\omega_n t}\right) =\omega_n^2 t e^{-\omega_n t}. \]

Since \( \omega_n^2 t e^{-\omega_n t} > 0 \) for all \( t > 0 \), the response is strictly increasing and approaches 1 from below. This formalizes the statement: \( \zeta=1 \) is the fastest non-oscillatory step response within the second-order family.

For \( \zeta > 1 \), the poles are distinct real numbers:

\[ s_{1,2}=-\zeta\omega_n \pm \omega_n\sqrt{\zeta^2-1}, \qquad s_1 > s_2,\;\; s_1<0,\;\; s_2<0. \]

The homogeneous error solution is:

\[ e(t)=C_1 e^{s_1 t}+C_2 e^{s_2 t}. \]

Apply \( e(0)=-1 \) and \( \dot{e}(0)=0 \):

\[ C_1+C_2=-1, \qquad s_1 C_1+s_2 C_2=0. \]

Solving the linear system yields:

\[ C_1=\frac{s_2}{s_2-s_1}, \qquad C_2=\frac{s_1}{s_1-s_2}. \]

Hence the overdamped unit-step response is:

\[ y(t)=1+e(t) =1+\frac{s_2 e^{s_1 t}-s_1 e^{s_2 t}}{s_1-s_2}, \qquad (\zeta > 1). \]

Proof that the response is monotone increasing. Differentiate:

\[ \dot{y}(t)=\frac{s_1 s_2}{s_1-s_2}\left(e^{s_1 t}-e^{s_2 t}\right). \]

Because \( s_1<0 \) and \( s_2<0 \), we have \( s_1 s_2 > 0 \). Also, since \( s_1 > s_2 \), it follows that \( e^{s_1 t} > e^{s_2 t} \) for all \( t>0 \). Finally, \( s_1-s_2 > 0 \), hence \( \dot{y}(t) > 0 \) for all \( t>0 \). Therefore the overdamped step response increases monotonically to 1 and cannot overshoot.

flowchart TD
  P["Poles from s^2 + 2*zeta*wn*s + wn^2 = 0"] --> U["0 < zeta < 1: poles = -zeta*wn ± j*wd"]
  P --> C["zeta = 1: double pole at -wn"]
  P --> O["zeta > 1: poles = s1, s2 \n(both real, negative)"]
  U --> UE["Time response: decaying oscillation \n(exp envelope + sin/cos)"]
  C --> CE["Time response: fastest \nnon-oscillatory (t*exp term)"]
  O --> OE["Time response: \nsum of two decaying \nexponentials (two modes)"]
        

Below, we implement the unit-step response \( y(t) \) in two ways: (i) via control libraries (transfer function → step response), and (ii) via the derived closed forms (from scratch), using careful branching near \( \zeta=1 \).

6.1 Python (SciPy / python-control) + Closed-Form Implementation

import numpy as np

def step_response_closed_form(t, wn, zeta, tol=1e-8):
    t = np.asarray(t, dtype=float)

    if zeta < 1.0 - tol:
        wd = wn * np.sqrt(1.0 - zeta**2)
        return 1.0 - np.exp(-zeta*wn*t) * (np.cos(wd*t) + (zeta/np.sqrt(1.0 - zeta**2))*np.sin(wd*t))

    elif abs(zeta - 1.0) <= tol:
        return 1.0 - np.exp(-wn*t) * (1.0 + wn*t)

    else:
        a = np.sqrt(zeta**2 - 1.0)
        s1 = -wn*(zeta - a)
        s2 = -wn*(zeta + a)
        return 1.0 + (s2*np.exp(s1*t) - s1*np.exp(s2*t)) / (s1 - s2)

# Library-based validation (SciPy)
if __name__ == "__main__":
    import matplotlib.pyplot as plt
    from scipy import signal

    wn, zeta = 5.0, 0.3
    sys = signal.TransferFunction([wn**2], [1.0, 2.0*zeta*wn, wn**2])

    t = np.linspace(0, 5, 2000)
    tout, y_lib = signal.step(sys, T=t)
    y_cf = step_response_closed_form(t, wn, zeta)

    plt.figure()
    plt.plot(tout, y_lib, label="SciPy step()")
    plt.plot(t, y_cf, "--", label="Closed form")
    plt.xlabel("t (s)")
    plt.ylabel("y(t)")
    plt.grid(True)
    plt.legend()
    plt.show()

Typical Python control workflows use scipy.signal (included with SciPy) or the control package (python-control) for transfer-function modeling and response simulation.

6.2 C++ (From Scratch Closed Form)

#include <cmath>
#include <iostream>
#include <vector>

double stepResponseClosedForm(double t, double wn, double zeta, double tol = 1e-8) {
    if (zeta < 1.0 - tol) {
        double wd = wn * std::sqrt(1.0 - zeta*zeta);
        double expo = std::exp(-zeta*wn*t);
        double term = std::cos(wd*t) + (zeta/std::sqrt(1.0 - zeta*zeta))*std::sin(wd*t);
        return 1.0 - expo * term;
    } else if (std::abs(zeta - 1.0) <= tol) {
        return 1.0 - std::exp(-wn*t) * (1.0 + wn*t);
    } else {
        double a = std::sqrt(zeta*zeta - 1.0);
        double s1 = -wn*(zeta - a);
        double s2 = -wn*(zeta + a);
        return 1.0 + (s2*std::exp(s1*t) - s1*std::exp(s2*t)) / (s1 - s2);
    }
}

int main() {
    double wn = 5.0;
    std::vector<double> zetas = {0.2, 1.0, 2.0};

    for (double zeta : zetas) {
        std::cout << "zeta=" << zeta << "\n";
        for (int k = 0; k <= 10; ++k) {
            double t = 0.1 * k;
            double y = stepResponseClosedForm(t, wn, zeta);
            std::cout << "  t=" << t << "  y=" << y << "\n";
        }
        std::cout << "\n";
    }
    return 0;
}

For larger system-dynamics projects in C++, common supporting libraries include Eigen (linear algebra) and Boost (numerics). Here, the second-order closed form avoids any dependency while remaining numerically stable if you treat the neighborhood of \( \zeta=1 \) with a tolerance.

6.3 Java (Closed Form + Notes on Libraries)

public class SecondOrderStepResponse {

    public static double stepResponseClosedForm(double t, double wn, double zeta, double tol) {
        if (zeta < 1.0 - tol) {
            double wd = wn * Math.sqrt(1.0 - zeta*zeta);
            double expo = Math.exp(-zeta*wn*t);
            double term = Math.cos(wd*t) + (zeta/Math.sqrt(1.0 - zeta*zeta))*Math.sin(wd*t);
            return 1.0 - expo * term;
        } else if (Math.abs(zeta - 1.0) <= tol) {
            return 1.0 - Math.exp(-wn*t) * (1.0 + wn*t);
        } else {
            double a = Math.sqrt(zeta*zeta - 1.0);
            double s1 = -wn*(zeta - a);
            double s2 = -wn*(zeta + a);
            return 1.0 + (s2*Math.exp(s1*t) - s1*Math.exp(s2*t)) / (s1 - s2);
        }
    }

    public static void main(String[] args) {
        double wn = 5.0;
        double[] zetas = {0.2, 1.0, 2.0};
        double tol = 1e-8;

        for (double zeta : zetas) {
            System.out.println("zeta=" + zeta);
            for (int k = 0; k <= 10; k++) {
                double t = 0.1 * k;
                double y = stepResponseClosedForm(t, wn, zeta, tol);
                System.out.println("  t=" + t + "  y=" + y);
            }
            System.out.println();
        }
    }
}

For more general simulation (beyond closed forms), Java users often rely on Apache Commons Math for ODE integrators and numerics. In this lesson, the closed-form approach directly encodes the theory derived in Sections 3–5.

6.4 MATLAB + Simulink (Transfer Function and Step)

% MATLAB: second-order step response using Control System Toolbox
wn = 5; zeta = 0.3;
G = tf(wn^2, [1 2*zeta*wn wn^2]);

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

% Closed-form verification
t = linspace(0, 5, 2000);
wd = wn*sqrt(1 - zeta^2);
y_cf = 1 - exp(-zeta*wn*t).*(cos(wd*t) + (zeta/sqrt(1-zeta^2))*sin(wd*t));

hold on;
plot(t, y_cf, '--');
legend('step(G)', 'closed form');

Simulink guidance (conceptual, no images). Use a Transfer Fcn block with numerator [wn^2] and denominator [1 2*zeta*wn wn^2]. Drive it with a Step block (step time 0, initial 0, final 1) and view output on a Scope. This exactly matches the mathematical model \( G(s) \) used above.

6.5 Wolfram Mathematica (TransferFunctionModel + StepResponse)

(* Mathematica: second-order step response *)
wn = 5; zeta = 0.3;
G = TransferFunctionModel[wn^2/(s^2 + 2 zeta wn s + wn^2), s];

(* Step response *)
StepResponse[G, {t, 0, 5}]

(* Closed form (piecewise) *)
yClosed[t_] := Piecewise[{
  {1 - Exp[-zeta wn t] (Cos[wn Sqrt[1 - zeta^2] t] + (zeta/Sqrt[1 - zeta^2]) Sin[wn Sqrt[1 - zeta^2] t]), 0 < zeta < 1},
  {1 - Exp[-wn t] (1 + wn t), zeta == 1},
  {Module[{a, s1, s2},
     a = Sqrt[zeta^2 - 1];
     s1 = -wn (zeta - a); s2 = -wn (zeta + a);
     1 + (s2 Exp[s1 t] - s1 Exp[s2 t])/(s1 - s2)
   ], zeta > 1}
}];

Plot[{StepResponse[G, t], yClosed[t]}, {t, 0, 5}, PlotLegends -> {"Library", "Closed form"}]

Mathematica’s control-systems functions provide an immediate symbolic/numeric check against the derived expressions, which is particularly useful when studying limiting behavior as \( \zeta \) approaches 1.

Problem 1 (Classify and write the step response): Consider \( G(s)=\dfrac{25}{s^2+4s+25} \). Find \( \omega_n \) and \( \zeta \), classify the response, and write the unit-step response \( y(t) \) in closed form.

Solution: Match coefficients with \( s^2+2\zeta\omega_n s+\omega_n^2 \).

\[ \omega_n^2=25 \;\; → \;\; \omega_n=5, \qquad 2\zeta\omega_n=4 \;\; → \;\; \zeta=\frac{4}{2\cdot 5}=0.4. \]

Since \( 0 < \zeta < 1 \), the system is underdamped. Compute \( \omega_d=\omega_n\sqrt{1-\zeta^2}=5\sqrt{1-0.16}=5\sqrt{0.84} \). The unit-step response is:

\[ y(t)=1-e^{-0.4\cdot 5\, t}\left(\cos(\omega_d t)+\frac{0.4}{\sqrt{1-0.4^2}}\sin(\omega_d t)\right), \quad \omega_d=5\sqrt{0.84}. \]

Problem 2 (Compute peak time and overshoot from the underdamped form): For \( \omega_n=10 \), \( \zeta=0.2 \), compute \( \omega_d \), the first peak time \( t_p \), and the overshoot ratio \( M_p \).

Solution:

\[ \omega_d=\omega_n\sqrt{1-\zeta^2}=10\sqrt{1-0.04}=10\sqrt{0.96}. \]

\[ t_p=\frac{\pi}{\omega_d}=\frac{\pi}{10\sqrt{0.96}}. \]

\[ M_p=\exp\!\left(\frac{-\pi\zeta}{\sqrt{1-\zeta^2}}\right) =\exp\!\left(\frac{-\pi\cdot 0.2}{\sqrt{0.96}}\right). \]

These follow directly from the canonical underdamped step-response structure; they quantify the oscillatory character implied by complex poles.

Problem 3 (Monotonicity of critical damping): Show that for \( \zeta=1 \) the unit-step response \( y(t)=1-(1+\omega_n t)e^{-\omega_n t} \) satisfies \( 0 \le y(t) < 1 \) and is strictly increasing for all \( t>0 \).

Solution: For \( t\ge 0 \), \( (1+\omega_n t)e^{-\omega_n t} \ge 0 \), hence \( y(t)\le 1 \). Also \( y(0)=0 \). Differentiate:

\[ \dot{y}(t)=\omega_n^2 t e^{-\omega_n t}. \]

Since \( \omega_n^2 t e^{-\omega_n t} > 0 \) for \( t>0 \), the response is strictly increasing. Therefore it approaches 1 from below with no overshoot.

Problem 4 (Overdamped derivative sign): Let \( \zeta>1 \), and define poles \( s_1>s_2 \) as in Section 5. Starting from \( y(t)=1+\dfrac{s_2 e^{s_1 t}-s_1 e^{s_2 t}}{s_1-s_2} \), prove \( \dot{y}(t)>0 \) for all \( t>0 \).

Solution: Differentiate explicitly:

\[ \dot{y}(t)=\frac{s_1 s_2}{s_1-s_2}\left(e^{s_1 t}-e^{s_2 t}\right). \]

For \( \zeta>1 \), both poles are real negative, so \( s_1 s_2 > 0 \). Also, since \( s_1>s_2 \), we have \( e^{s_1 t} > e^{s_2 t} \) for all \( t>0 \). Finally, \( s_1-s_2 > 0 \). Hence every factor is positive, so \( \dot{y}(t)>0 \).

Problem 5 (Estimating damping ratio from decay of peaks): In an underdamped step response, suppose two successive peak deviations from steady state satisfy \( |e(t_k)|/|e(t_{k+1})| = R \) with \( R>1 \). Show that the logarithmic decrement \( \delta=\ln(R) \) satisfies: \( \delta=\dfrac{2\pi\zeta}{\sqrt{1-\zeta^2}} \), and solve for \( \zeta \) in terms of \( \delta \).

Solution: For \( 0 < \zeta < 1 \), the envelope of oscillation is proportional to \( e^{-\zeta\omega_n t} \). Successive peaks are separated by one damped period \( T_d=\dfrac{2\pi}{\omega_d} \), where \( \omega_d=\omega_n\sqrt{1-\zeta^2} \). Therefore:

\[ \frac{|e(t_k)|}{|e(t_{k+1})|}=\frac{e^{-\zeta\omega_n t_k}}{e^{-\zeta\omega_n (t_k+T_d)}}=e^{\zeta\omega_n T_d} \;\; → \;\; \delta=\ln(R)=\zeta\omega_n\frac{2\pi}{\omega_d}. \]

\[ \delta=\zeta\omega_n\frac{2\pi}{\omega_n\sqrt{1-\zeta^2}} =\frac{2\pi\zeta}{\sqrt{1-\zeta^2}}. \]

Solve for \( \zeta \) by squaring:

\[ \delta^2=\frac{4\pi^2\zeta^2}{1-\zeta^2} \;\; → \;\; \delta^2(1-\zeta^2)=4\pi^2\zeta^2 \;\; → \;\; \delta^2=\zeta^2(\delta^2+4\pi^2). \]

\[ \zeta=\frac{\delta}{\sqrt{\delta^2+4\pi^2}}. \]

We derived the unit-step responses for the canonical second-order system by reducing the problem to the homogeneous error dynamics and solving the characteristic equation. The discriminant (equivalently, the damping ratio) partitions responses into: underdamped (complex poles and decaying oscillations), critically damped (repeated pole and fastest non-oscillatory response), and overdamped (two real modes and monotone behavior). Finally, we implemented the closed forms and validated them against standard control libraries across multiple programming environments.