In previous chapters, we studied time responses (step, impulse, arbitrary inputs) and frequency responses (Bode, Nyquist, Nichols) of single-input single-output (SISO) LTI systems mainly from an analytic perspective. In practical control engineering, especially in robotics, these responses are almost always computed numerically by software tools from a model \( G(s) \) (transfer function) or directly from differential equations.

For a stable, strictly proper SISO LTI system with transfer function \( G(s) = \dfrac{N(s)}{D(s)} \) and input \( u(t) \), the output \( y(t) \) satisfies

\[ Y(s) = G(s)U(s), \quad y(t) = \mathcal{L}^{-1}\{G(s)U(s)\}(t). \]

In the time domain, this is equivalently expressed as convolution of the impulse response \( g(t) \) with the input:

\[ g(t) = \mathcal{L}^{-1}\{G(s)\}(t), \quad y(t) = (g * u)(t) = \int_0^t g(\tau)\,u(t - \tau)\,\mathrm{d}\tau. \]

In the frequency domain, for sinusoidal steady state, the frequency response \( G(j\omega) \) is obtained by evaluating the rational function along the imaginary axis:

\[ G(j\omega) = \left. G(s)\right|_{s = j\omega}, \quad \omega \ge 0. \]

Automated analysis replaces hand calculations by:

• Algorithmic evaluation of \( y(t) \) on a time grid.
• Algorithmic evaluation of \( G(j\omega) \) on a frequency grid.
• Automatic extraction of performance metrics (rise time, overshoot, bandwidth, margins).

Understanding the underlying formulas is crucial: software tools are essentially implementing these mathematical operations with finite precision and discretization.

Consider a strictly proper rational transfer function with real coefficients and distinct poles:

\[ G(s) = \frac{N(s)}{D(s)} = \frac{b_m s^m + \cdots + b_0}{s^n + a_{n-1}s^{n-1} + \cdots + a_0}, \quad n > m. \]

Suppose the denominator factors as \( D(s) = \prod_{i=1}^n (s - p_i) \) with pairwise distinct poles \( p_i \) satisfying \( \Re(p_i) < 0 \) (asymptotic stability). Partial fractions give

\[ G(s) = \sum_{i=1}^n \frac{k_i}{s - p_i}, \quad k_i = \lim_{s \to p_i} (s - p_i)G(s). \]

Taking the inverse Laplace transform yields the impulse response:

\[ g(t) = \mathcal{L}^{-1}\{G(s)\}(t) = \sum_{i=1}^n k_i e^{p_i t} u(t), \]

where \( u(t) \) is the unit step function. For a unit-step input \( u(t) = u(t) \) with Laplace transform \( U(s) = \dfrac{1}{s} \), we obtain

\[ Y(s) = \frac{G(s)}{s}, \quad y(t) = \mathcal{L}^{-1}\left\{\frac{G(s)}{s}\right\}(t). \]

For the canonical underdamped second-order system

\[ G(s) = \frac{\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2}, \quad \omega_n > 0,\; 0 < \zeta < 1, \]

the impulse and step responses are

\[ g(t) = \frac{\omega_n}{\sqrt{1 - \zeta^2}} e^{-\zeta \omega_n t} \sin\!\big(\omega_d t\big) u(t), \quad \omega_d = \omega_n \sqrt{1 - \zeta^2}, \]

\[ y_{\text{step}}(t) = 1 - \frac{1}{\sqrt{1 - \zeta^2}} e^{-\zeta \omega_n t} \sin\!\big(\omega_d t + \phi\big) u(t), \quad \phi = \arctan\!\left(\frac{\sqrt{1 - \zeta^2}}{\zeta}\right). \]

Automated time-response tools (e.g., step in MATLAB or step_response in Python) are effectively evaluating these expressions (or their numerical analogues) for a wide range of models without requiring the user to perform partial-fraction decompositions manually.

For a general input \( u(t) \), tools compute \( y(t) \) via numerical convolution or numerical integration of the associated differential equation, which we now formalize.

An \( n \)th-order linear differential equation with input \( u(t) \) can be written as

\[ y^{(n)}(t) + a_{n-1} y^{(n-1)}(t) + \cdots + a_0 y(t) = b_m u^{(m)}(t) + \cdots + b_0 u(t), \]

where derivatives of \( u(t) \) may or may not appear (for strictly proper systems we typically have \( m \le n-1 \)). We collect the output and its derivatives into a state-like vector

\[ \mathbf{x}(t) = \begin{bmatrix} y(t) \\ \dot{y}(t) \\ \vdots \\ y^{(n-1)}(t) \end{bmatrix} \in \mathbb{R}^n. \]

This leads to a first-order vector differential equation

\[ \dot{\mathbf{x}}(t) = f\big(t,\mathbf{x}(t),u(t)\big), \quad y(t) = e_1^{\top}\mathbf{x}(t), \]

where \( e_1 \) is the first standard basis vector. Numerical time simulation uses a time grid \( t_k = k h \) with step size \( h > 0 \), approximating the solution by a sequence \( \mathbf{x}_k \approx \mathbf{x}(t_k) \).

The simplest method is explicit (forward) Euler:

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

Using Taylor expansion,

\[ \mathbf{x}(t_k + h) = \mathbf{x}(t_k) + h \dot{\mathbf{x}}(t_k) + \frac{h^2}{2}\ddot{\mathbf{x}}(\xi_k), \quad t_k < \xi_k < t_k + h, \]

and substituting \( \dot{\mathbf{x}}(t_k) = f(t_k,\mathbf{x}(t_k),u(t_k)) \), we see that Euler's method has local truncation error \( \mathcal{O}(h^2) \) and global error \( \mathcal{O}(h) \) under standard Lipschitz conditions.

For higher accuracy, tools use higher-order schemes such as classical Runge–Kutta of order four (RK4):

\[ \begin{aligned} \mathbf{k}_1 &= f(t_k,\mathbf{x}_k,u_k),\\ \mathbf{k}_2 &= f\!\big(t_k + \tfrac{h}{2}, \mathbf{x}_k + \tfrac{h}{2}\mathbf{k}_1, u_k\big),\\ \mathbf{k}_3 &= f\!\big(t_k + \tfrac{h}{2}, \mathbf{x}_k + \tfrac{h}{2}\mathbf{k}_2, u_k\big),\\ \mathbf{k}_4 &= f(t_k + h, \mathbf{x}_k + h\mathbf{k}_3, u_k),\\[4pt] \mathbf{x}_{k+1} &= \mathbf{x}_k + \frac{h}{6}\left( \mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4 \right). \end{aligned} \]

Automated time-response functions like lsim (MATLAB) or forced_response (Python's control library) encapsulate this process:

• convert transfer functions to an internal ODE representation,
• choose an integration method and step size,
• integrate from zero initial conditions,
• return \( y(t_k) \) samples for further analysis.

Time-domain performance metrics are then computed from the sampled trajectory. For example, letting \( y_k = y(t_k) \) and denoting the approximate steady-state value by \( y_{\text{ss}} \approx y_{N} \) for sufficiently large \( N \):

\[ M_p = \frac{\max_k y_k - y_{\text{ss}}}{y_{\text{ss}}}, \quad t_p = t_{k^\star} \text{ where } y_{k^\star} = \max_k y_k, \]

\[ t_s(\delta) = \min\left\{t_k \,\middle|\, \big|y_\ell - y_{\text{ss}}\big| \le \delta\, y_{\text{ss}} \text{ for all } \ell \ge k \right\}, \]

where \( \delta = 0.02 \) or \( 0.05 \) for 2% or 5% settling time. These definitions are implemented internally by many control toolboxes.

flowchart TD
  A["Specify transfer function G(s)"] --> B["Choose input: step / impulse / custom u(t)"]
  B --> C["Create time grid (t0..tN) and integrator options"]
  C --> D["Convert model to ODE form xdot = f(t,x,u)"]
  D --> E["Integrate numerically to obtain x_k and y_k"]
  E --> F["Compute metrics: Mp, tr, ts, steady-state error"]
  F --> G["Generate plots and export data for robotics controller tuning"]
        

For a stable LTI system with impulse response \( g(t) \), the frequency response is the Fourier transform of \( g(t) \):

\[ G(j\omega) = \int_0^{\infty} g(t) e^{-j\omega t}\,\mathrm{d}t. \]

When \( G(s) \) is rational and stable, it is numerically preferable to evaluate the rational function directly at \( s = j\omega \) rather than perform a numerical integral. Let

\[ G(j\omega) = A(\omega) + j B(\omega). \]

The magnitude and phase are then

\[ |G(j\omega)| = \sqrt{A(\omega)^2 + B(\omega)^2}, \quad \angle G(j\omega) = \operatorname{atan2}\big(B(\omega), A(\omega)\big). \]

Bode plots use logarithmic amplitude and log-spaced frequency:

\[ M_{\text{dB}}(\omega) = 20 \log_{10}|G(j\omega)|, \quad \omega_k = 10^{\alpha + k\Delta}, \]

\[ \Delta = \frac{\beta - \alpha}{N-1}, \quad \omega_{\min} = 10^{\alpha},\; \omega_{\max} = 10^{\beta},\; k = 0,\dots,N-1. \]

In practice, \( \omega_{\min} \) and \( \omega_{\max} \) are chosen relative to the real parts of poles and zeros. A common heuristic is

\[ \omega_{\min} \approx 0.1 \min_i \big|\Re(p_i)\big|, \quad \omega_{\max} \approx 10 \max_i \big|\Re(p_i)\big|, \]

ensuring that the frequency range covers all dynamic phenomena relevant for control design (bandwidth, resonant peaks, crossover frequencies).

For the first-order system

\[ G(s) = \frac{1}{\tau s + 1}, \quad \tau > 0, \]

we obtain

\[ G(j\omega) = \frac{1}{1 + j\omega \tau}, \quad |G(j\omega)| = \frac{1}{\sqrt{1 + (\omega \tau)^2}}, \quad \angle G(j\omega) = -\arctan(\omega \tau). \]

At the break frequency \( \omega = \dfrac{1}{\tau} \), the exact magnitude is

\[ M_{\text{dB}}\!\left(\frac{1}{\tau}\right) = 20\log_{10}\!\left(\frac{1}{\sqrt{2}}\right) \approx -3.01 \,\text{dB}, \]

giving the characteristic \( -3\,\text{dB} \) point of a first-order low-pass system used extensively in loop-shaping.

Automated functions such as bode (MATLAB) or bode_plot (Python control libraries) internally:

• build a logarithmic frequency grid,
• evaluate \( G(j\omega_k) \) for each grid point,
• return magnitude and phase arrays for plotting and margin extraction.

flowchart TD
  F0["Specify LTI model G(s)"] --> F1["Choose log-spaced frequency grid w_k"]
  F1 --> F2["Evaluate G(j w_k) for all k"]
  F2 --> F3["Compute magnitude |G| and phase angle"]
  F3 --> F4["Generate Bode, Nyquist, Nichols plots"]
  F4 --> F5["Extract bandwidth, gain and phase margins for controller tuning"]
        

Python, together with numpy, scipy, and the control library, is widely used for prototyping control algorithms in robotics (often inside ROS or ROS 2 nodes via rospy / rclpy). We illustrate automated time- and frequency-response computation for the canonical second-order system:

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

import numpy as np
import matplotlib.pyplot as plt

# Core control libraries
import control  # python-control: pip install control
from control import TransferFunction, step_response, impulse_response, forced_response, bode_plot

# Parameters of a typical servo joint (robotics-like)
omega_n = 20.0   # rad/s
zeta    = 0.4    # damping ratio

# Define transfer function G(s) = omega_n^2 / (s^2 + 2*zeta*omega_n*s + omega_n^2)
num = [omega_n**2]
den = [1.0, 2.0*zeta*omega_n, omega_n**2]
G = TransferFunction(num, den)

# 1) Automated step and impulse responses
t_step, y_step = step_response(G)
t_imp,  y_imp  = impulse_response(G)

# 2) Forced response to a sinusoidal input u(t) = sin(omega_in * t)
omega_in = 5.0  # rad/s
t = np.linspace(0.0, 4.0, 2000)
u = np.sin(omega_in * t)
t_forced, y_forced, x_forced = forced_response(G, T=t, U=u)

# 3) Frequency response (Bode data) on a log-spaced grid
w = np.logspace(-1, 2, 1000)  # 0.1 .. 100 rad/s
mag, phase, omega = control.freqresp(G, w)

# Convert magnitude to dB
mag_db = 20.0 * np.log10(np.squeeze(mag))
phase_deg = np.squeeze(np.degrees(phase))
omega = np.squeeze(omega)

# 4) Plotting (time domain)
plt.figure()
plt.plot(t_step, y_step)
plt.title("Step response")
plt.xlabel("t [s]")
plt.ylabel("y(t)")
plt.grid(True)

plt.figure()
plt.plot(t_imp, y_imp)
plt.title("Impulse response")
plt.xlabel("t [s]")
plt.ylabel("g(t)")
plt.grid(True)

plt.figure()
plt.plot(t_forced, u, label="u(t) = sin(omega_in * t)")
plt.plot(t_forced, y_forced, label="y(t)")
plt.title("Forced response to sinusoidal input")
plt.xlabel("t [s]")
plt.legend()
plt.grid(True)

# 5) Plotting (Bode-like plots)
plt.figure()
plt.semilogx(omega, mag_db)
plt.title("Magnitude response |G(jw)| in dB")
plt.xlabel("omega [rad/s]")
plt.ylabel("Magnitude [dB]")
plt.grid(True, which="both")

plt.figure()
plt.semilogx(omega, phase_deg)
plt.title("Phase response arg G(jw)")
plt.xlabel("omega [rad/s]")
plt.ylabel("Phase [deg]")
plt.grid(True, which="both")

plt.show()
      

In a robotics context, this code would often be wrapped inside a ROS 2 node using rclpy, where the time response of a joint servo loop (with transfer function identified from experiments or from a simplified rigid-body model) is evaluated for various gains. The same model can be used in simulation (e.g., Gazebo) and in the real controller to ensure consistency.

In real-time robotics controllers (e.g., implemented within ros_control or Orocos), C++ is typically used. Time-response simulation can be implemented using explicit numerical integration for the underlying differential equations. For the same second-order system

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

define the state vector \( x_1(t) = y(t) \), \( x_2(t) = \dot{y}(t) \). Then

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

The following C++ code uses RK4 to compute the unit-step response (suitable for offline analysis or for generating reference data for controller tuning). Note the escaping of angle brackets for HTML:

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

struct State {
    double x1; // position y
    double x2; // velocity dy/dt
};

State deriv(const State& x, double t, double omega_n, double zeta, double u) {
    (void)t; // t not used in this time-invariant example
    State dx;
    dx.x1 = x.x2;
    dx.x2 = -2.0 * zeta * omega_n * x.x2
            - omega_n * omega_n * x.x1
            + omega_n * omega_n * u;
    return dx;
}

int main() {
    double omega_n = 20.0;
    double zeta    = 0.4;
    double h       = 0.001;   // step size [s]
    double t_final = 4.0;     // final time [s]

    int n_steps = static_cast<int>(t_final / h);
    std::vector<double> t(n_steps + 1);
    std::vector<double> y(n_steps + 1);

    State x{0.0, 0.0}; // initial conditions y(0) = 0, dy/dt(0) = 0

    for (int k = 0; k <= n_steps; ++k) {
        t[k] = k * h;
        y[k] = x.x1;

        double tk = t[k];
        double u  = 1.0; // unit-step input

        State k1 = deriv(x, tk,          omega_n, zeta, u);
        State k2 = deriv({x.x1 + 0.5*h*k1.x1, x.x2 + 0.5*h*k1.x2},
                         tk + 0.5*h, omega_n, zeta, u);
        State k3 = deriv({x.x1 + 0.5*h*k2.x1, x.x2 + 0.5*h*k2.x2},
                         tk + 0.5*h, omega_n, zeta, u);
        State k4 = deriv({x.x1 + h*k3.x1, x.x2 + h*k3.x2},
                         tk + h, omega_n, zeta, u);

        x.x1 += (h / 6.0) * (k1.x1 + 2.0*k2.x1 + 2.0*k3.x1 + k4.x1);
        x.x2 += (h / 6.0) * (k1.x2 + 2.0*k2.x2 + 2.0*k3.x2 + k4.x2);
    }

    // Export y[k] and t[k] to file or plot (e.g., using gnuplot or Python)
    std::cout << "# t  y\n";
    for (int k = 0; k <= n_steps; ++k) {
        std::cout << t[k] << " " << y[k] << "\n";
    }
    return 0;
}
      

For frequency-response computation in C++, one can use std::complex<double> to evaluate \( G(j\omega) \) on a frequency grid (e.g. for pre-computing Bode data used in an offline robotics tuning tool).

Java is used less frequently for low-level control, but it appears in robotics middleware (e.g., ROSJava) and simulation environments. The following class implements a simple explicit Euler integrator for the same second-order system. While Euler is less accurate than RK4, it suffices for illustrating automation of time-response computation.

public class SecondOrderSystem {

    private final double omegaN;
    private final double zeta;

    public SecondOrderSystem(double omegaN, double zeta) {
        this.omegaN = omegaN;
        this.zeta = zeta;
    }

    public static class Response {
        public final double[] t;
        public final double[] y;

        public Response(double[] t, double[] y) {
            this.t = t;
            this.y = y;
        }
    }

    // Simulate unit-step response y(t) for t in [0, tFinal] with step h
    public Response simulateStep(double tFinal, double h) {
        int nSteps = (int) Math.round(tFinal / h);
        double[] t = new double[nSteps + 1];
        double[] y = new double[nSteps + 1];

        double x1 = 0.0; // y
        double x2 = 0.0; // dy/dt

        for (int k = 0; k <= nSteps; ++k) {
            t[k] = k * h;
            y[k] = x1;

            double u = 1.0; // unit step
            double dx1 = x2;
            double dx2 = -2.0 * zeta * omegaN * x2
                         - omegaN * omegaN * x1
                         + omegaN * omegaN * u;

            x1 += h * dx1;
            x2 += h * dx2;
        }
        return new Response(t, y);
    }

    public static void main(String[] args) {
        SecondOrderSystem sys = new SecondOrderSystem(20.0, 0.4);
        Response resp = sys.simulateStep(4.0, 0.001);

        // Print first few samples (can be plotted in MATLAB, Python, etc.)
        for (int k = 0; k < resp.t.length; ++k) {
            System.out.println(resp.t[k] + " " + resp.y[k]);
        }
    }
}
      

In a robotics setting, such a time-response simulation may be embedded into higher-level Java tools that assist in tuning controller gains for ROSJava-based mobile robots or manipulators. More sophisticated projects would use numerical libraries such as EJML for linear algebra and ODE solvers, but the underlying principle remains the same.

8.1 MATLAB/Simulink

MATLAB and Simulink are standard tools in control and robotics. The basic continuous-time workflow for the second-order system is:

% Parameters
omega_n = 20;    % rad/s
zeta    = 0.4;

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

% 1) Time-domain responses
figure; step(G);
title('Step response');

figure; impulse(G);
title('Impulse response');

t = 0:0.001:4;
u = sin(5*t);   % sinusoidal input (e.g., joint torque disturbance)
[y, t_out] = lsim(G, u, t);

figure;
plot(t_out, u, 'k--', t_out, y, 'b');
legend('u(t)', 'y(t)');
xlabel('t [s]');
ylabel('Amplitude');
title('Forced response to sinusoidal input');
grid on;

% 2) Frequency response (Bode and Nyquist)
figure; bode(G);
grid on;

figure; nyquist(G);
grid on;

% 3) Robotics-oriented example:
%    Combine with Robotics System Toolbox to analyze joint servo dynamics
%    for a rigid-body model robot = robotics.RigidBodyTree;
%    Closed-loop joint model can be mapped to a transfer function and
%    analyzed with the same commands.
      

In Simulink, the same transfer function is implemented using a Transfer Fcn block, excited by a Step or Sine Wave source, and analyzed with a Scope block. Simulation is invoked from MATLAB with:

simOut = sim('SecondOrderServoModel');  % name of the Simulink model file
t_sim  = simOut.tout;
y_sim  = simOut.yout.signals.values;
      

8.2 Wolfram Mathematica

Wolfram Mathematica provides symbolic and numeric capabilities for LTI systems, which is useful for deriving closed-form expressions and then comparing them to numerical results.

(* Parameters *)
omegaN = 20.0;
zeta   = 0.4;

(* Define transfer function G(s) = omegaN^2 / (s^2 + 2*zeta*omegaN*s + omegaN^2) *)
s =. ;
G = TransferFunctionModel[
      omegaN^2 / (s^2 + 2*zeta*omegaN*s + omegaN^2),
      s
    ];

(* 1) Step response y_step(t) *)
yStep[t_] = OutputResponse[G, UnitStep[t], {t, 0, 4}];

(* 2) Impulse response g(t) *)
g[t_] = OutputResponse[G, DiracDelta[t], {t, 0, 4}];

(* 3) Bode plot *)
bodePlot = BodePlot[G, {10^-1, 10^2},
  PlotLayout -> "MagnitudePhase",
  GridLines -> Automatic
];

(* 4) Symbolic check of step response formula *)
assumptions = {omegaN > 0, 0 < zeta < 1};
yStepSymbolic = InverseLaplaceTransform[
  (omegaN^2 / (s^2 + 2*zeta*omegaN*s + omegaN^2)) / s,
  s, t, Assumptions -> assumptions
];
Simplify[yStepSymbolic, Assumptions -> assumptions]
      

Mathematica can thus verify analytic formulas for time responses and provide high-precision numeric data for benchmarking implementations in embedded robotics controllers.

Problem 1 (Analytic and Numeric Step Response): Consider the underdamped second-order system

\[ G(s) = \frac{\omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2}, \quad \omega_n > 0,\; 0 < \zeta < 1. \]

(a) Derive the unit-step response \( y_{\text{step}}(t) \).
(b) Show that the percent overshoot is \( M_p = \exp\!\left(-\dfrac{\pi \zeta}{\sqrt{1-\zeta^2}}\right) \).
(c) Propose an algorithm based on sampled data \( y_k = y(t_k) \) that estimates \( M_p \) numerically.

Solution:

(a) From Section 2, we already know the standard result

\[ y_{\text{step}}(t) = 1 - \frac{1}{\sqrt{1 - \zeta^2}} e^{-\zeta \omega_n t} \sin\!\big(\omega_d t + \phi\big) u(t), \]

\[ \omega_d = \omega_n \sqrt{1 - \zeta^2},\quad \phi = \arctan\!\left(\frac{\sqrt{1 - \zeta^2}}{\zeta}\right). \]

(b) The first peak occurs when \( \dfrac{\mathrm{d}y_{\text{step}}(t)}{\mathrm{d} t} = 0 \), which yields \( t_p = \dfrac{\pi}{\omega_d} \). Evaluating \( y_{\text{step}}(t_p) \) gives

\[ y_{\text{step}}(t_p) = 1 + M_p, \quad M_p = \exp\!\left( -\frac{\pi \zeta}{\sqrt{1-\zeta^2}} \right). \]

(c) Given sampled data, compute the approximate steady-state value \( y_{\text{ss}} \approx y_{N} \) for large \( N \), and the maximum sample \( y_{\max} = \max_k y_k \). Then

\[ \hat{M}_p = \frac{y_{\max} - y_{\text{ss}}}{y_{\text{ss}}}. \]

In software, this is implemented by scanning the vector of samples, which is exactly how, for example, step-response analysis tools extract percent overshoot.

Problem 2 (Frequency Response of a First-Order System): For the first-order system

\[ G(s) = \frac{1}{\tau s + 1}, \quad \tau > 0, \]

(a) Derive expressions for \( |G(j\omega)| \) and \( \angle G(j\omega) \).
(b) Show that at \( \omega = \dfrac{1}{\tau} \), \( |G(j\omega)| = \dfrac{1}{\sqrt{2}} \).
(c) Compute the corresponding magnitude in dB and interpret this in the context of Bode plots.

Solution:

(a) Substituting \( s = j\omega \) gives

\[ G(j\omega) = \frac{1}{1 + j\omega \tau} = \frac{1 - j\omega \tau}{1 + (\omega \tau)^2}. \]

Hence

\[ |G(j\omega)| = \frac{1}{\sqrt{1 + (\omega \tau)^2}}, \quad \angle G(j\omega) = -\arctan(\omega \tau). \]

(b) At \( \omega = \dfrac{1}{\tau} \),

\[ |G(j/\tau)| = \frac{1}{\sqrt{1 + (1)^2}} = \frac{1}{\sqrt{2}}. \]

(c) The dB magnitude is

\[ M_{\text{dB}} = 20\log_{10}\!\left(\frac{1}{\sqrt{2}}\right) = -10\log_{10}(2) \approx -3.01 \,\text{dB}. \]

This explains the \( -3\,\text{dB} \) corner convention for first-order low-pass systems (and for bandwidth definitions in loop shaping).

Problem 3 (Discrete Convolution Approximation): Let the impulse response be \( g(t) = e^{-t} u(t) \) and the input \( u(t) = u(t) \) (unit step). The exact output is

\[ y(t) = \int_0^t e^{-(t-\tau)}\,\mathrm{d}\tau = 1 - e^{-t}. \]

(a) Derive the above expression.
(b) Consider a discrete convolution with step size \( h \):

\[ y_k^{(h)} = h \sum_{i=0}^k g(t_i) u(t_k - t_i), \quad t_i = i h. \]

Show that \( y_k^{(h)} \to y(t_k) \) as \( h \to 0 \), and discuss the order of convergence.

Solution:

(a) Since \( u(t) = 1 \) for \( t \ge 0 \),

\[ y(t) = \int_0^t e^{-(t-\tau)}\,\mathrm{d}\tau = e^{-t}\int_0^t e^{\tau}\,\mathrm{d}\tau = e^{-t}\big(e^t - 1\big) = 1 - e^{-t}. \]

(b) The discrete sum is a Riemann approximation of the integral. Under continuity of \( g(t)u(t_k - t) \), the trapezoidal or rectangular rule yields \( y_k^{(h)} \to y(t_k) \) as \( h \to 0 \). For the rectangular rule used above, the global error is \( \mathcal{O}(h) \), analogous to the error of Euler integration. This is the principle behind discrete convolution algorithms in control libraries.

Problem 4 (Choosing a Frequency Grid from Poles and Zeros): A transfer function has poles at \( p_1 = -1 \), \( p_2 = -100 \) and a zero at \( z_1 = -10 \). Propose reasonable values for \( \omega_{\min} \) and \( \omega_{\max} \) for a Bode plot, and define a log-spaced frequency grid \( \{\omega_k\}_{k=0}^{N-1} \). Justify your choice.

Solution:

The dominant dynamics occur around the magnitudes of poles and zeros: \( 1 \), \( 10 \), and \( 100 \) rad/s. A common heuristic is to extend one decade below the smallest and one decade above the largest magnitude:

\[ \omega_{\min} = 0.1,\quad \omega_{\max} = 1000. \]

For \( N \) frequency points, define

\[ \omega_k = 10^{\alpha + k\Delta}, \quad \alpha = \log_{10}(0.1) = -1,\quad \beta = \log_{10}(1000) = 3, \]

\[ \Delta = \frac{\beta - \alpha}{N - 1} = \frac{4}{N - 1}, \quad k = 0,\dots,N-1. \]

This ensures that the break frequencies at \( \omega = 1, 10, 100 \) are well covered, allowing automated routines to resolve bandwidth and stability margins accurately.

Problem 5 (Estimating Damping Ratio from a Simulated Step Response): Assume you have simulated the step response of an underdamped second-order system (as in Problems 1 and 2) and obtained successive peak values \( y_{p,1} \) and \( y_{p,2} \), with steady-state value \( y_{\text{ss}} \).

(a) Derive the logarithmic decrement \( \delta = \ln\left( \dfrac{y_{p,1} - y_{\text{ss}}}{y_{p,2} - y_{\text{ss}}} \right) \).
(b) Show that

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

(c) Explain how a control tool could automatically estimate \( \zeta \) from simulated data.

Solution:

(a) For an underdamped second-order system, successive peak deviations from steady state satisfy \( y_{p,k} - y_{\text{ss}} = (y_{p,1} - y_{\text{ss}}) e^{-(k-1)\delta} \), where \( \delta \) is the logarithmic decrement. For \( k = 2 \), we get

\[ y_{p,2} - y_{\text{ss}} = (y_{p,1} - y_{\text{ss}}) e^{-\delta} \quad \Rightarrow \quad \delta = \ln\left( \frac{y_{p,1} - y_{\text{ss}}}{y_{p,2} - y_{\text{ss}}} \right). \]

(b) Classical second-order system theory shows \( \delta = \dfrac{2\pi \zeta}{\sqrt{1 - \zeta^2}} \), which can be solved for \( \zeta \):

\[ \delta^2 = \frac{4\pi^2 \zeta^2}{1 - \zeta^2} \quad \Rightarrow \quad \delta^2 (1 - \zeta^2) = 4\pi^2 \zeta^2, \]

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

(c) A software tool can automatically locate consecutive peaks from sampled step-response data, compute \( \delta \) and then estimate \( \zeta \) using the formula above. This is commonly used in automated system identification and tuning for servo drives and robotic joint controllers.

In this lesson we connected the analytic theory of time and frequency responses of LTI systems with their automated computation in modern software tools. We expressed time responses via inverse Laplace transforms and convolution, and we discussed numerical ODE integration (Euler, RK4) underlying time-domain simulations. We showed how frequency responses are computed by evaluating \( G(j\omega) \) on a log-spaced grid to generate Bode and Nyquist plots.

We then illustrated concrete implementations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, with comments on how these workflows appear in robotics controller design. Finally, we solved problems that practice analytic and numeric reasoning about time and frequency responses, including estimation of damping ratio from simulated data and frequency-grid design from pole/zero locations. These concepts form the basis for the software-assisted design techniques developed in subsequent lessons on root locus and frequency-domain design.