Let an LTI system be described (from previous chapters) by a transfer function \( G(s)=\dfrac{Y(s)}{U(s)} \) under zero initial conditions. This lesson formalizes the mapping from a sinusoidal input at frequency \( \omega \) to its steady-state output and uses that mapping to define the frequency response.

Intuitively, if the system is stable (all poles strictly in the left half-plane), any transient response decays. What remains is a forced response at the input frequency. The entire frequency-response toolbox (Bode, Nyquist, resonance measures) builds on this statement.

flowchart TD
  A["Model (ODE / state-space)"] --> B["Transfer function G(s)"]
  B --> C["Set s = j*w (evaluate on imaginary axis)"]
  C --> D["Complex gain: G(j*w)"]
  D --> E["Input sinusoid u(t)=Um*sin(w*t+phi)"]
  E --> F["Steady-state output y_ss(t)=|G|*Um*sin(w*t+phi+angle(G))"]
  F --> G["Frequency response: w -> G(j*w)"]
        

Two technical details matter and will appear in our proofs:

• Existence of \( G(j\omega) \): \( j\omega \) must not be a pole of \( G(s) \).
• Decay of transients: system poles must satisfy \( \operatorname{Re}(p_i) < 0 \).

A sinusoid can be represented as the real part of a complex exponential. Define a complex phasor \( U \in \mathbb{C} \) such that

\[ u(t) = U_m \cos(\omega t + \varphi) = \operatorname{Re}\!\left\{ U e^{j\omega t} \right\}, \quad U := U_m e^{j\varphi}. \]

This representation is powerful because differentiation corresponds to multiplication by \( j\omega \):

\[ \frac{d}{dt}\left(e^{j\omega t}\right) = j\omega e^{j\omega t}, \quad \frac{d^k}{dt^k}\left(e^{j\omega t}\right) = (j\omega)^k e^{j\omega t}. \]

Therefore, when we analyze an LTI system’s response to \( e^{j\omega t} \), the steady-state response becomes algebraic in \( j\omega \). We will exploit this in both transfer-function and state-space proofs.

Theorem 1 (Sinusoidal Steady-State Response). Consider a proper rational transfer function \( G(s) \) with all poles satisfying \( \operatorname{Re}(p_i) < 0 \) (BIBO stability). Let the input be the complex exponential \( u_c(t)=U e^{j\omega t}\mathbf{1}(t) \), where \( \mathbf{1}(t) \) is the unit step. If \( j\omega \) is not a pole of \( G(s) \), then the output satisfies

\[ y_c(t) = U\,G(j\omega)e^{j\omega t} + y_{\text{tr}}(t), \quad \text{with} \quad \lim_{t→\infty} y_{\text{tr}}(t)=0. \]

For a real sinusoidal input \( u(t)=\operatorname{Re}\{u_c(t)\} \), the real output is \( y(t)=\operatorname{Re}\{y_c(t)\} \), hence the steady-state output is

\[ y_{\text{ss}}(t)=|G(j\omega)|\,U_m \cos\!\big(\omega t+\varphi+\angle G(j\omega)\big). \]

Proof. Take the unilateral Laplace transform (zero initial conditions are built into the transfer function framework):

\[ U_c(s) = \mathcal{L}\{U e^{j\omega t}\mathbf{1}(t)\} = \frac{U}{s-j\omega}, \quad Y_c(s)=G(s)U_c(s)=\frac{U\,G(s)}{s-j\omega}. \]

Because \( j\omega \) is not a pole of \( G(s) \), \( Y_c(s) \) has a simple pole at \( s=j\omega \) whose residue is \( U\,G(j\omega) \). Let the poles of \( G(s) \) be \( \{p_i\}_{i=1}^n \) (possibly repeated). Then \( Y_c(s) \) has poles at \( s=j\omega \) and at \( s=p_i \). The inverse Laplace transform yields a decomposition of the form

\[ y_c(t) = U\,G(j\omega)e^{j\omega t} + \sum_{i=1}^n \sum_{k=1}^{m_i} c_{ik}\, t^{k-1} e^{p_i t}, \]

where \( m_i \) is the multiplicity of pole \( p_i \) and the coefficients \( c_{ik} \) are constants determined by partial fractions (residues). Since \( \operatorname{Re}(p_i) < 0 \) for all poles, each term \( t^{k-1}e^{p_i t} → 0 \) as \( t → \infty \). Hence the transient sum decays to zero, leaving

\[ \lim_{t→\infty}\left(y_c(t)-U\,G(j\omega)e^{j\omega t}\right)=0. \]

Finally, because the system is real and LTI, real and imaginary parts propagate linearly; for \( u(t)=\operatorname{Re}\{U e^{j\omega t}\} \) we have \( y(t)=\operatorname{Re}\{U\,G(j\omega)e^{j\omega t}\} \), which converts directly to the magnitude/phase expression. \(\square\)

Important edge cases. If \( j\omega \) is a pole of \( G(s) \), the forced response may grow without bound (a resonance-like phenomenon in the marginal case), and sinusoidal steady state is not guaranteed.

Definition (Frequency Response). For an LTI system with transfer function \( G(s) \), the frequency response is the complex-valued function of real frequency \( \omega \) defined by

\[ G(j\omega) := G(s)\big|_{s=j\omega}, \quad \omega \in \mathbb{R}, \quad \text{provided } j\omega \text{ is not a pole of } G(s). \]

If a sinusoidal input has phasor \( U \) (so \( u(t)=\operatorname{Re}\{Ue^{j\omega t}\} \)), then the steady-state output phasor is \( Y = G(j\omega)\,U \). This motivates interpreting \( G(j\omega) \) as a complex gain.

The magnitude and phase are defined as:

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

For real-coefficient transfer functions, complex conjugation implies a symmetry:

\[ G(-j\omega) = \overline{G(j\omega)}, \quad |G(-j\omega)| = |G(j\omega)|, \quad \angle G(-j\omega) = -\angle G(j\omega). \]

First-order example (interpretation). For a standard first-order model \( G(s)=\dfrac{1}{1+\tau s} \), the frequency response is

\[ 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). \]

Thus higher frequencies are attenuated and phase-lagged more strongly, matching the time-domain intuition that a slow first-order system cannot track rapid changes.

Consider a continuous-time LTI state-space model (introduced earlier): \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u,\; y=\mathbf{C}\mathbf{x}+\mathbf{D}u \). For complex sinusoidal forcing \( u_c(t)=Ue^{j\omega t} \), seek a particular solution of the form \( \mathbf{x}_p(t)=\mathbf{X}e^{j\omega t} \). Substitution gives:

\[ j\omega \mathbf{X}e^{j\omega t} = \mathbf{A}\mathbf{X}e^{j\omega t} + \mathbf{B}Ue^{j\omega t} \;\;→\;\; (j\omega \mathbf{I}-\mathbf{A})\mathbf{X} = \mathbf{B}U. \]

If \( \det(j\omega \mathbf{I}-\mathbf{A}) \neq 0 \) (i.e., \( j\omega \) is not an eigenvalue of \( \mathbf{A} \)), then

\[ \mathbf{X} = (j\omega \mathbf{I}-\mathbf{A})^{-1}\mathbf{B}U, \quad Y = \mathbf{C}\mathbf{X} + \mathbf{D}U \;\;→\;\; \frac{Y}{U} = \mathbf{C}(j\omega \mathbf{I}-\mathbf{A})^{-1}\mathbf{B} + \mathbf{D}. \]

Hence the state-space frequency response is:

\[ G(j\omega)=\mathbf{C}(j\omega \mathbf{I}-\mathbf{A})^{-1}\mathbf{B} + \mathbf{D}. \]

This expression is numerically important: it avoids explicitly forming transfer functions for high-order systems and directly links resonance behavior to near-singularity of \( j\omega \mathbf{I}-\mathbf{A} \).

In engineering practice, you typically compute \( G(j\omega) \) from the model and then verify it using time simulation (or experimental data). A robust numerical approach is to estimate the steady-state amplitude and phase by fitting the final segment of the simulated output to \( a\sin(\omega t)+b\cos(\omega t) \).

If \( y(t) \approx a\sin(\omega t)+b\cos(\omega t) \), then it can be rewritten as \( y(t)\approx R\sin(\omega t+\phi) \) with

\[ R=\sqrt{a^2+b^2}, \quad \phi=\operatorname{atan2}(b,a). \]

flowchart TD
  S0["Choose model: G(s) or (A,B,C,D)"] --> S1["Pick test frequency w and input amplitude Um"]
  S1 --> S2["Compute Gjw = G(j*w)"]
  S2 --> S3["Predict: Ym = |Gjw|*Um, phase = angle(Gjw)"]
  S3 --> S4["Simulate time response y(t) with u(t)=Um*sin(w*t)"]
  S4 --> S5["Discard transient (late-time window)"]
  S5 --> S6["Fit y ~= a*sin(w*t) + b*cos(w*t)"]
  S6 --> S7["Compute R=sqrt(a^2+b^2), phi=atan2(b,a)"]
  S7 --> S8["Compare (R,phi) with prediction (Ym,phase)"]
        

Libraries: numpy, scipy.signal (and optionally python-control). The script below evaluates \( G(j\omega) \) for a second-order model and validates it by time simulation and sinusoid fitting.

File: Chapter12_Lesson1.py

"""
Chapter12_Lesson1.py
System Dynamics (Control Engineering) — Chapter 12, Lesson 1
Sinusoidal Steady-State Response and Frequency Response Definition

Requires: numpy, scipy (scipy.signal), matplotlib (optional for plots)
Optional: control (python-control) for alternative frequency-response utilities.
"""

import numpy as np
from numpy import pi
from scipy import signal

def tf_eval(num, den, s):
    """
    Evaluate a rational transfer function G(s) = N(s)/D(s) at complex s.
    num, den: 1D arrays of polynomial coefficients in descending powers.
    """
    num = np.asarray(num, dtype=np.complex128)
    den = np.asarray(den, dtype=np.complex128)
    return np.polyval(num, s) / np.polyval(den, s)

def fit_sinusoid(t, y, omega):
    """
    Fit y(t) ≈ a*sin(omega*t) + b*cos(omega*t) by least squares.
    Returns amplitude R and phase phi such that y ≈ R*sin(omega*t + phi).
    """
    A = np.column_stack([np.sin(omega * t), np.cos(omega * t)])
    coeff, *_ = np.linalg.lstsq(A, y, rcond=None)
    a, b = coeff
    R = np.sqrt(a*a + b*b)
    phi = np.arctan2(b, a)  # because R*sin(wt+phi)=R*cos(phi)*sin(wt)+R*sin(phi)*cos(wt)
    return R, phi, a, b

def main():
    # Example: standard 2nd-order low-pass
    wn = 5.0        # rad/s
    zeta = 0.2
    num = [wn**2]
    den = [1.0, 2.0*zeta*wn, wn**2]
    G = signal.TransferFunction(num, den)

    # Sinusoidal input
    Um = 1.0
    omega = 4.0  # rad/s
    phi_u = 0.0  # rad

    # Frequency response prediction
    Gjw = tf_eval(num, den, 1j*omega)
    Ym_pred = abs(Gjw) * Um
    phi_y_pred = np.angle(Gjw) + phi_u

    print("G(jw) =", Gjw)
    print("|G(jw)| =", abs(Gjw), "  angle(G(jw)) =", np.angle(Gjw))
    print("Predicted steady-state amplitude Ym =", Ym_pred)
    print("Predicted steady-state phase phi_y (rad) =", phi_y_pred)

    # Time-domain simulation (zero initial conditions)
    t = np.linspace(0.0, 40.0, 40001)  # long horizon to let transients decay
    u = Um * np.sin(omega*t + phi_u)
    tout, y, _ = signal.lsim(G, U=u, T=t)

    # Estimate steady state from final window (e.g., last 10 seconds)
    mask = tout >= (tout[-1] - 10.0)
    R_hat, phi_hat, a_hat, b_hat = fit_sinusoid(tout[mask], y[mask], omega)

    # Normalize phase to be comparable
    def wrap_to_pi(x):
        return (x + np.pi) % (2*np.pi) - np.pi

    print("\nEstimated from simulation (last 10 s):")
    print("Ym_hat =", R_hat)
    print("phi_y_hat (rad) =", wrap_to_pi(phi_hat))

    print("\nErrors:")
    print("Amplitude error:", R_hat - Ym_pred)
    print("Phase error (rad):", wrap_to_pi(phi_hat - phi_y_pred))

    # Optional plotting (uncomment if desired)
    # import matplotlib.pyplot as plt
    # plt.figure()
    # plt.plot(tout, u, label="u(t)")
    # plt.plot(tout, y, label="y(t)")
    # plt.xlim((tout[-1]-5.0, tout[-1]))
    # plt.legend()
    # plt.xlabel("t (s)")
    # plt.grid(True)
    # plt.show()

if __name__ == "__main__":
    main()

Libraries: standard C++17 only (we implement complex arithmetic via <complex> and RK4 from scratch). For large projects you may use Eigen for linear algebra, but it is not required here.

File: Chapter12_Lesson1.cpp

/*
Chapter12_Lesson1.cpp
System Dynamics (Control Engineering) — Chapter 12, Lesson 1
Sinusoidal Steady-State Response and Frequency Response Definition

This program:
1) Evaluates a 2nd-order transfer function at s = j*omega to obtain G(jw).
2) Simulates the equivalent ODE with RK4 under sinusoidal forcing.
3) Fits the steady-state output to extract amplitude and phase.

Build (example):
  g++ -O2 -std=c++17 Chapter12_Lesson1.cpp -o Chapter12_Lesson1

No external dependencies.
*/

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

static constexpr double PI = 3.14159265358979323846;

std::complex<double> polyval(const std::vector<double>& c, std::complex<double> s) {
    // c in descending powers
    std::complex<double> y = 0.0;
    for (double a : c) y = y * s + a;
    return y;
}

std::complex<double> tf_eval(const std::vector<double>& num,
                             const std::vector<double>& den,
                             std::complex<double> s) {
    return polyval(num, s) / polyval(den, s);
}

struct FitResult {
    double R;   // amplitude
    double phi; // phase (rad) in y ≈ R sin(omega t + phi)
    double a;   // sin coeff
    double b;   // cos coeff
};

// Least squares fit y ≈ a sin(wt) + b cos(wt)
FitResult fit_sinusoid(const std::vector<double>& t,
                       const std::vector<double>& y,
                       double omega) {
    double s11 = 0.0, s12 = 0.0, s22 = 0.0;
    double r1  = 0.0, r2  = 0.0;

    for (size_t i = 0; i < t.size(); ++i) {
        double si = std::sin(omega * t[i]);
        double ci = std::cos(omega * t[i]);
        s11 += si * si;
        s12 += si * ci;
        s22 += ci * ci;
        r1  += si * y[i];
        r2  += ci * y[i];
    }

    // Solve 2x2 normal equations:
    // [s11 s12; s12 s22] [a; b] = [r1; r2]
    double det = s11 * s22 - s12 * s12;
    double a = ( r1 * s22 - r2 * s12) / det;
    double b = (-r1 * s12 + r2 * s11) / det;

    double R = std::sqrt(a*a + b*b);
    double phi = std::atan2(b, a); // because R sin(wt+phi)=R cos(phi) sin(wt)+R sin(phi) cos(wt)

    return {R, phi, a, b};
}

double wrap_to_pi(double x) {
    while (x >  PI) x -= 2.0*PI;
    while (x < -PI) x += 2.0*PI;
    return x;
}

int main() {
    // Second-order low-pass: G(s)=wn^2/(s^2+2*zeta*wn*s+wn^2)
    const double wn   = 5.0;
    const double zeta = 0.2;

    std::vector<double> num{wn*wn};
    std::vector<double> den{1.0, 2.0*zeta*wn, wn*wn};

    // Input: u(t) = Um sin(omega t + phi_u)
    const double Um   = 1.0;
    const double omega = 4.0;
    const double phi_u = 0.0;

    std::complex<double> Gjw = tf_eval(num, den, std::complex<double>(0.0, omega));
    double Ym_pred = std::abs(Gjw) * Um;
    double phi_y_pred = std::arg(Gjw) + phi_u;

    std::cout << "G(jw) = " << Gjw << "\n";
    std::cout << "|G(jw)| = " << std::abs(Gjw) << "  angle(G(jw)) = " << std::arg(Gjw) << "\n";
    std::cout << "Predicted steady-state amplitude Ym = " << Ym_pred << "\n";
    std::cout << "Predicted steady-state phase phi_y (rad) = " << phi_y_pred << "\n\n";

    // ODE: y'' + 2*zeta*wn*y' + wn^2*y = wn^2*u(t)
    // State: x1=y, x2=y'
    auto f = [&](double t, double x1, double x2) {
        double u = Um * std::sin(omega*t + phi_u);
        double dx1 = x2;
        double dx2 = -2.0*zeta*wn*x2 - wn*wn*x1 + wn*wn*u;
        return std::pair<double,double>(dx1, dx2);
    };

    // RK4 simulation
    double t0 = 0.0, tf = 40.0;
    double dt = 0.001;
    size_t N = static_cast<size_t>((tf - t0)/dt) + 1;

    std::vector<double> tvec;
    std::vector<double> yvec;
    tvec.reserve(N);
    yvec.reserve(N);

    double x1 = 0.0, x2 = 0.0;
    double t = t0;
    for (size_t k = 0; k < N; ++k) {
        tvec.push_back(t);
        yvec.push_back(x1);

        auto [k1_1, k1_2] = f(t, x1, x2);
        auto [k2_1, k2_2] = f(t + 0.5*dt, x1 + 0.5*dt*k1_1, x2 + 0.5*dt*k1_2);
        auto [k3_1, k3_2] = f(t + 0.5*dt, x1 + 0.5*dt*k2_1, x2 + 0.5*dt*k2_2);
        auto [k4_1, k4_2] = f(t + dt, x1 + dt*k3_1, x2 + dt*k3_2);

        x1 += (dt/6.0) * (k1_1 + 2.0*k2_1 + 2.0*k3_1 + k4_1);
        x2 += (dt/6.0) * (k1_2 + 2.0*k2_2 + 2.0*k3_2 + k4_2);

        t += dt;
    }

    // Use last 10 seconds for steady-state fitting
    std::vector<double> tss, yss;
    for (size_t i = 0; i < tvec.size(); ++i) {
        if (tvec[i] >= tf - 10.0) {
            tss.push_back(tvec[i]);
            yss.push_back(yvec[i]);
        }
    }

    FitResult fr = fit_sinusoid(tss, yss, omega);

    std::cout << "Estimated from simulation (last 10 s):\n";
    std::cout << "Ym_hat = " << fr.R << "\n";
    std::cout << "phi_y_hat (rad) = " << wrap_to_pi(fr.phi) << "\n\n";

    std::cout << "Errors:\n";
    std::cout << "Amplitude error: " << (fr.R - Ym_pred) << "\n";
    std::cout << "Phase error (rad): " << wrap_to_pi(fr.phi - phi_y_pred) << "\n";

    return 0;
}

Libraries: standard Java only. In larger projects, you can use Apache Commons Math for complex and linear algebra, but this lesson keeps everything self-contained.

File: Chapter12_Lesson1.java

/*
Chapter12_Lesson1.java
System Dynamics (Control Engineering) — Chapter 12, Lesson 1
Sinusoidal Steady-State Response and Frequency Response Definition

This program:
1) Evaluates a 2nd-order transfer function at s = j*omega to obtain G(jw).
2) Simulates the equivalent ODE with RK4 under sinusoidal forcing.
3) Fits the steady-state output to extract amplitude and phase.

Compile and run:
  javac Chapter12_Lesson1.java
  java Chapter12_Lesson1

No external dependencies.
*/

import java.util.ArrayList;
import java.util.List;

public class Chapter12_Lesson1 {

    static final double PI = 3.14159265358979323846;

    // Minimal complex class
    static class Complex {
        final double re, im;

        Complex(double re, double im) { this.re = re; this.im = im; }

        Complex add(Complex z) { return new Complex(this.re + z.re, this.im + z.im); }
        Complex sub(Complex z) { return new Complex(this.re - z.re, this.im - z.im); }
        Complex mul(Complex z) {
            return new Complex(this.re*z.re - this.im*z.im, this.re*z.im + this.im*z.re);
        }
        Complex div(Complex z) {
            double d = z.re*z.re + z.im*z.im;
            return new Complex((this.re*z.re + this.im*z.im)/d, (this.im*z.re - this.re*z.im)/d);
        }
        double abs() { return Math.hypot(this.re, this.im); }
        double arg() { return Math.atan2(this.im, this.re); }

        public String toString() {
            return String.format("(%.6f %+.6f j)", re, im);
        }
    }

    static Complex polyval(double[] c, Complex s) {
        // c in descending powers
        Complex y = new Complex(0.0, 0.0);
        for (double a : c) {
            y = y.mul(s).add(new Complex(a, 0.0));
        }
        return y;
    }

    static Complex tfEval(double[] num, double[] den, Complex s) {
        return polyval(num, s).div(polyval(den, s));
    }

    static class FitResult {
        final double R, phi, a, b;
        FitResult(double R, double phi, double a, double b) {
            this.R = R; this.phi = phi; this.a = a; this.b = b;
        }
    }

    // Least squares fit y ≈ a sin(wt) + b cos(wt)
    static FitResult fitSinusoid(List<Double> t, List<Double> y, double omega) {
        double s11 = 0.0, s12 = 0.0, s22 = 0.0;
        double r1  = 0.0, r2  = 0.0;

        for (int i = 0; i < t.size(); i++) {
            double si = Math.sin(omega * t.get(i));
            double ci = Math.cos(omega * t.get(i));
            double yi = y.get(i);

            s11 += si * si;
            s12 += si * ci;
            s22 += ci * ci;
            r1  += si * yi;
            r2  += ci * yi;
        }

        double det = s11 * s22 - s12 * s12;
        double a = ( r1 * s22 - r2 * s12) / det;
        double b = (-r1 * s12 + r2 * s11) / det;

        double R = Math.sqrt(a*a + b*b);
        double phi = Math.atan2(b, a);

        return new FitResult(R, phi, a, b);
    }

    static double wrapToPi(double x) {
        while (x >  PI) x -= 2.0*PI;
        while (x < -PI) x += 2.0*PI;
        return x;
    }

    public static void main(String[] args) {
        // Second-order low-pass: G(s)=wn^2/(s^2+2*zeta*wn*s+wn^2)
        double wn   = 5.0;
        double zeta = 0.2;

        double[] num = new double[]{wn*wn};
        double[] den = new double[]{1.0, 2.0*zeta*wn, wn*wn};

        // Input: u(t) = Um sin(omega t + phi_u)
        double Um = 1.0;
        double omega = 4.0;
        double phi_u = 0.0;

        Complex Gjw = tfEval(num, den, new Complex(0.0, omega));
        double Ym_pred = Gjw.abs() * Um;
        double phi_y_pred = Gjw.arg() + phi_u;

        System.out.println("G(jw) = " + Gjw);
        System.out.println("|G(jw)| = " + Gjw.abs() + "  angle(G(jw)) = " + Gjw.arg());
        System.out.println("Predicted steady-state amplitude Ym = " + Ym_pred);
        System.out.println("Predicted steady-state phase phi_y (rad) = " + phi_y_pred);
        System.out.println();

        // ODE: y'' + 2*zeta*wn*y' + wn^2*y = wn^2*u(t)
        // State: x1=y, x2=y'
        double t0 = 0.0, tf = 40.0, dt = 0.001;
        int N = (int)Math.round((tf - t0)/dt) + 1;

        List<Double> tvec = new ArrayList<>(N);
        List<Double> yvec = new ArrayList<>(N);

        double x1 = 0.0, x2 = 0.0;
        double t = t0;

        for (int k = 0; k < N; k++) {
            tvec.add(t);
            yvec.add(x1);

            // RK4
            double[] k1 = f(t, x1, x2, Um, omega, phi_u, wn, zeta);
            double[] k2 = f(t + 0.5*dt, x1 + 0.5*dt*k1[0], x2 + 0.5*dt*k1[1], Um, omega, phi_u, wn, zeta);
            double[] k3 = f(t + 0.5*dt, x1 + 0.5*dt*k2[0], x2 + 0.5*dt*k2[1], Um, omega, phi_u, wn, zeta);
            double[] k4 = f(t + dt, x1 + dt*k3[0], x2 + dt*k3[1], Um, omega, phi_u, wn, zeta);

            x1 += (dt/6.0) * (k1[0] + 2.0*k2[0] + 2.0*k3[0] + k4[0]);
            x2 += (dt/6.0) * (k1[1] + 2.0*k2[1] + 2.0*k3[1] + k4[1]);

            t += dt;
        }

        // Fit last 10 seconds
        List<Double> tss = new ArrayList<>();
        List<Double> yss = new ArrayList<>();
        for (int i = 0; i < tvec.size(); i++) {
            if (tvec.get(i) >= tf - 10.0) {
                tss.add(tvec.get(i));
                yss.add(yvec.get(i));
            }
        }

        FitResult fr = fitSinusoid(tss, yss, omega);

        System.out.println("Estimated from simulation (last 10 s):");
        System.out.println("Ym_hat = " + fr.R);
        System.out.println("phi_y_hat (rad) = " + wrapToPi(fr.phi));
        System.out.println();

        System.out.println("Errors:");
        System.out.println("Amplitude error: " + (fr.R - Ym_pred));
        System.out.println("Phase error (rad): " + wrapToPi(fr.phi - phi_y_pred));
    }

    static double[] f(double t, double x1, double x2,
                      double Um, double omega, double phi_u,
                      double wn, double zeta) {
        double u = Um * Math.sin(omega*t + phi_u);
        double dx1 = x2;
        double dx2 = -2.0*zeta*wn*x2 - wn*wn*x1 + wn*wn*u;
        return new double[]{dx1, dx2};
    }
}

Toolboxes: Control System Toolbox (required for tf, freqresp, lsim), Simulink (optional). The script computes \( G(j\omega) \), predicts steady state, simulates the time response, and estimates amplitude/phase by least squares.

File: Chapter12_Lesson1.m

% Chapter12_Lesson1.m
% System Dynamics (Control Engineering) — Chapter 12, Lesson 1
% Sinusoidal Steady-State Response and Frequency Response Definition
%
% Requires: Control System Toolbox (tf, bode, freqresp, lsim)
% Simulink section is optional and requires Simulink.
%
% This script:
% 1) Computes G(jw) and predicts steady-state amplitude/phase.
% 2) Simulates y(t) via lsim and estimates amplitude/phase by least squares.
% 3) (Optional) builds and simulates a Simulink model programmatically.

clear; clc;

wn   = 5.0;          % rad/s
zeta = 0.2;
G = tf([wn^2],[1 2*zeta*wn wn^2]);

Um    = 1.0;
omega = 4.0;         % rad/s
phi_u = 0.0;

% Frequency response at jw
Gjw = squeeze(freqresp(G, omega));   % complex scalar
Ym_pred = abs(Gjw) * Um;
phi_y_pred = angle(Gjw) + phi_u;

fprintf('G(jw) = %.6f%+.6fj\n', real(Gjw), imag(Gjw));
fprintf('|G(jw)| = %.6f, angle(G(jw)) = %.6f rad\n', abs(Gjw), angle(Gjw));
fprintf('Predicted steady-state amplitude Ym = %.6f\n', Ym_pred);
fprintf('Predicted steady-state phase phi_y = %.6f rad\n\n', phi_y_pred);

% Time simulation
t = linspace(0,40,40001).';
u = Um * sin(omega*t + phi_u);
y = lsim(G, u, t);

% Fit last 10 seconds: y ≈ a sin(wt) + b cos(wt)
mask = t >= (t(end)-10);
ts = t(mask); ys = y(mask);
A = [sin(omega*ts), cos(omega*ts)];
coef = A \ ys;
a = coef(1); b = coef(2);
Ym_hat = sqrt(a^2 + b^2);
phi_hat = atan2(b,a);

wrapToPi = @(x) mod(x + pi, 2*pi) - pi;

fprintf('Estimated from simulation (last 10 s):\n');
fprintf('Ym_hat = %.6f\n', Ym_hat);
fprintf('phi_y_hat = %.6f rad\n\n', wrapToPi(phi_hat));

fprintf('Errors:\n');
fprintf('Amplitude error: %.6f\n', Ym_hat - Ym_pred);
fprintf('Phase error (rad): %.6f\n', wrapToPi(phi_hat - phi_y_pred));

% Optional plots
% figure; plot(t,u,'LineWidth',1); hold on; plot(t,y,'LineWidth',1);
% grid on; xlabel('t (s)'); legend('u(t)','y(t)');
% xlim([t(end)-5 t(end)]);

% Optional: Build a Simulink model programmatically
%{
mdl = 'Chapter12_Lesson1_Simulink';
if bdIsLoaded(mdl); close_system(mdl,0); end
new_system(mdl); open_system(mdl);

add_block('simulink/Sources/Sine Wave',[mdl '/Sine']);
set_param([mdl '/Sine'], 'Amplitude', num2str(Um), 'Frequency', num2str(omega), 'Phase', num2str(phi_u));

add_block('simulink/Continuous/Transfer Fcn',[mdl '/G(s)']);
set_param([mdl '/G(s)'], 'Numerator', mat2str([wn^2]), 'Denominator', mat2str([1 2*zeta*wn wn^2]));

add_block('simulink/Sinks/Scope',[mdl '/Scope']);

add_line(mdl,'Sine/1','G(s)/1');
add_line(mdl,'G(s)/1','Scope/1');

set_param(mdl,'StopTime','40');
sim(mdl);

% The Scope shows the time response; to compare with theory, focus on late-time behavior.
%}

Mathematica has native symbolic/numeric frequency-response support via TransferFunctionModel, FrequencyResponse, and plotting via BodePlot.

File: Chapter12_Lesson1.nb

(* Chapter12_Lesson1.nb
   System Dynamics (Control Engineering) — Chapter 12, Lesson 1
   Sinusoidal Steady-State Response and Frequency Response Definition

   This is a Mathematica notebook stored as a plain-text Notebook[] expression.
*)

Notebook[{
  Cell["System Dynamics — Chapter 12, Lesson 1", "Title"],
  Cell["Sinusoidal Steady-State Response and Frequency Response", "Section"],

  Cell[BoxData[
    RowBox[{
      RowBox[{"wn", "=", "5.0"}], ";",
      RowBox[{"zeta", "=", "0.2"}], ";",
      RowBox[{"G", "=", 
        RowBox[{"TransferFunctionModel", "[", 
          RowBox[{
            FractionBox[RowBox[{"wn", "^", "2"}], 
              RowBox[{
                RowBox[{"s", "^", "2"}], "+", RowBox[{"2", " ", "zeta", " ", "wn", " ", "s"}], "+", RowBox[{"wn", "^", "2"}]
              }]
            ], ",", "s"}], "]"}]}], ";"
    }]], "Input"],

  Cell["Frequency response and Bode plot:", "Text"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"Um", "=", "1.0"}], ";",
      RowBox[{"omega", "=", "4.0"}], ";",
      RowBox[{"Gjw", "=", 
        RowBox[{"FrequencyResponse", "[", 
          RowBox[{"G", ",", "omega"}], "]"}]}], ";",
      RowBox[{"YmPred", "=", 
        RowBox[{"Abs", "[", "Gjw", "]"}], ";",
        }], ";",
      RowBox[{"phiPred", "=", 
        RowBox[{"Arg", "[", "Gjw", "]"}]}], ";",
      RowBox[{"{Gjw, YmPred, phiPred}"}]
    }]], "Input"],

  Cell[BoxData[RowBox[{"BodePlot", "[", RowBox[{"G", ",", RowBox[{"{", RowBox[{"0.1", ",", "100"}], "}"}]}], "]"}]], "Input"],

  Cell["Time-domain simulation of y'' + 2 zeta wn y' + wn^2 y = wn^2 u(t):", "Text"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"u", "[", "t_", "]"}], ":=", RowBox[{"Um", " ", RowBox[{"Sin", "[", RowBox[{"omega", " ", "t"}], "]"}]}], ";",
      RowBox[{"sol", "=", RowBox[{"NDSolve", "[", 
        RowBox[{
          RowBox[{"{", 
            RowBox[{
              RowBox[{
                RowBox[{
                  RowBox[{"y", "''"}], "[", "t", "]"}], "+", 
                RowBox[{"2", " ", "zeta", " ", "wn", " ", RowBox[{
                  RowBox[{"y", "'"}], "[", "t", "]"}]}], "+", 
                RowBox[{
                  RowBox[{"wn", "^", "2"}], " ", RowBox[{"y", "[", "t", "]"}]}], "==", 
                RowBox[{
                  RowBox[{"wn", "^", "2"}], " ", RowBox[{"u", "[", "t", "]"}]}
              }], ",",
              RowBox[{"y", "[", "0", "]"}], "==", "0", ",",
              RowBox[{
                RowBox[{"y", "'"}], "[", "0", "]"}], "==", "0"
            }], "}", ",", "y", ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "40"}], "}"}]}], "]"}]}], ";"
    }]], "Input"],

  Cell[BoxData[
    RowBox[{
      RowBox[{"Plot", "[", 
        RowBox[{
          RowBox[{"Evaluate", "[", RowBox[{"y", "[", "t", "]"}], "/.", "sol"}], ",",
          RowBox[{"{", RowBox[{"t", ",", "35", ",", "40"}], "}"}], ",",
          RowBox[{"GridLines", "->", "Automatic"}], ",",
          RowBox[{"PlotRange", "->", "All"}]}], "]"}]
    }]], "Input"]
},
WindowSize -> {1000, 700},
StyleDefinitions -> "Default.nb"
]

Problem 1 (First-order steady-state by direct ODE). Consider \( \tau \dot{y}(t) + y(t) = u(t) \) with \( \tau > 0 \) and input \( u(t)=U_m\cos(\omega t) \). Derive the steady-state output amplitude and phase.

Solution. Seek a particular solution \( y_p(t)=A\cos(\omega t)+B\sin(\omega t) \). Then \( \dot{y}_p(t)=-A\omega\sin(\omega t)+B\omega\cos(\omega t) \). Substitute into \( \tau \dot{y}+y=u \):

\[ \tau\big(-A\omega\sin(\omega t)+B\omega\cos(\omega t)\big) + A\cos(\omega t)+B\sin(\omega t) = U_m\cos(\omega t). \]

Match coefficients of \( \cos(\omega t) \) and \( \sin(\omega t) \):

\[ \begin{aligned} (\tau \omega)B + A &= U_m,\\ -(\tau \omega)A + B &= 0. \end{aligned} \]

Solve: from the second equation \( B=(\tau \omega)A \). Substitute into the first: \( (\tau \omega)^2 A + A = U_m \) so \( A = \dfrac{U_m}{1+(\omega \tau)^2} \), \( B = \dfrac{U_m(\omega \tau)}{1+(\omega \tau)^2} \). Convert to magnitude/phase:

\[ y_{\text{ss}}(t)=R\cos(\omega t - \phi),\quad R=\frac{U_m}{\sqrt{1+(\omega \tau)^2}}, \quad \phi=\arctan(\omega \tau). \]

This equals the frequency-response result for \( G(s)=\dfrac{1}{1+\tau s} \). \(\square\)

Problem 2 (Closed form for \( G(j\omega) \) of a 2nd-order model). For \( G(s)=\dfrac{\omega_n^2}{s^2+2\zeta\omega_n s + \omega_n^2} \), derive \( |G(j\omega)| \) and \( \angle G(j\omega) \).

Solution. Evaluate at \( s=j\omega \):

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

Hence

\[ |G(j\omega)|=\frac{\omega_n^2}{\sqrt{(\omega_n^2-\omega^2)^2 + (2\zeta\omega_n\omega)^2}}, \quad \angle G(j\omega)=-\operatorname{atan2}\!\big(2\zeta\omega_n\omega,\;\omega_n^2-\omega^2\big). \]

The magnitude becomes large when the denominator is small (this motivates resonance analysis in later lessons). \(\square\)

Problem 3 (Impulse response relation). Let \( g(t) \) be the impulse response of a causal LTI system with transfer function \( G(s)=\mathcal{L}\{g(t)\} \). Assume \( g(t) \) is absolutely integrable. Show that \( G(j\omega)=\int_{0}^{\infty} g(t)e^{-j\omega t}\,dt \).

Solution. By definition of the unilateral Laplace transform,

\[ G(s)=\int_{0}^{\infty} g(t)e^{-st}\,dt. \]

Under absolute integrability, the integral converges for \( s=j\omega \) and we can set \( s=j\omega \):

\[ G(j\omega)=\int_{0}^{\infty} g(t)e^{-j\omega t}\,dt, \]

which is exactly the (one-sided) Fourier transform of \( g(t) \) for causal systems. \(\square\)

Problem 4 (State-space frequency response). Prove that for \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u,\; y=\mathbf{C}\mathbf{x}+\mathbf{D}u \) and forcing \( u_c(t)=Ue^{j\omega t} \), a particular solution exists iff \( \det(j\omega \mathbf{I}-\mathbf{A})\neq 0 \), and then \( \dfrac{Y}{U}=\mathbf{C}(j\omega \mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D} \).

Solution. Use the ansatz \( \mathbf{x}_p(t)=\mathbf{X}e^{j\omega t} \) and substitute:

\[ j\omega \mathbf{X}e^{j\omega t}=\mathbf{A}\mathbf{X}e^{j\omega t}+\mathbf{B}Ue^{j\omega t} \;\;→\;\; (j\omega \mathbf{I}-\mathbf{A})\mathbf{X}=\mathbf{B}U. \]

This linear algebraic equation has a unique solution iff \( j\omega \mathbf{I}-\mathbf{A} \) is invertible, i.e. \( \det(j\omega \mathbf{I}-\mathbf{A})\neq 0 \). Then \( \mathbf{X}=(j\omega \mathbf{I}-\mathbf{A})^{-1}\mathbf{B}U \) and \( Y=\mathbf{C}\mathbf{X}+\mathbf{D}U \), giving the claim. \(\square\)

Problem 5 (Why stability is essential for “steady state”). Consider \( G(s)=\dfrac{1}{s} \) (an integrator) and input \( u(t)=\sin(\omega t) \). Determine \( y(t) \) and explain why sinusoidal steady state (in the strict sense) fails.

Solution. In time domain, \( \dot{y}(t)=u(t)=\sin(\omega t) \) so

\[ y(t)=\int_0^t \sin(\omega \tau)\,d\tau =\frac{1-\cos(\omega t)}{\omega}. \]

The output contains a constant (DC) offset \( 1/\omega \) in addition to a cosine term. This violates the “same-frequency-only” steady-state form that holds for strictly stable systems. In Laplace terms, the integrator has a pole at the origin (not satisfying \( \operatorname{Re}(p) < 0 \)), so transients do not decay in the required way. \(\square\)

We proved that strictly stable LTI systems map a sinusoidal input at frequency \( \omega \) to a sinusoidal steady-state output at the same frequency, scaled and phase-shifted by \( G(j\omega) \). We derived \( G(j\omega) \) from both transfer functions and state-space models, clarified existence conditions (no poles on \( j\omega \)), and implemented a predict-and-verify workflow in Python, C++, Java, MATLAB/Simulink, and Mathematica.