Consider a single-input single-output (SISO) LTI system described (in the time domain) by a constant-coefficient linear differential equation or, equivalently, by a transfer function \( G(s) \). We apply a sinusoidal input of the form \( u(t) = U \cos(\omega t + \varphi_u) \), where \( U > 0 \) is the amplitude and \( \omega \) is the angular frequency.

For an asymptotically stable LTI system, the total output decomposes into \( y(t) = y_{\text{tr}}(t) + y_{\text{ss}}(t) \), where \( y_{\text{tr}}(t) \) is a transient term that decays as \( t → \infty \), and \( y_{\text{ss}}(t) \) is the sinusoidal steady-state response. The central theorem of sinusoidal steady state is that:

\[ \text{If the LTI system is asymptotically stable and } u(t) = U \cos(\omega t + \varphi_u), \text{ then} \]

\[ y_{\text{ss}}(t) = \left| G(j\omega) \right| U \cos\!\big( \omega t + \varphi_u + \arg G(j\omega) \big). \]

Thus, sinusoidal inputs are eigenfunctions of LTI systems: the output is a sinusoid at the same frequency, scaled in amplitude and shifted in phase by the complex number \( G(j\omega) \).

flowchart TD
  U["Input: u(t) = U cos(omega t + phi_u)"] --> G["LTI system with G(s)"]
  G --> Ytr["Transient response y_tr(t) (decays if stable)"]
  G --> Yss["Steady-state y_ss(t)"]
  Yss --> FORM["Same frequency, amplitude = |G(jw)| U, phase shift = arg G(jw)"]
        

In this lesson, we derive this result using two complementary viewpoints: (i) complex exponentials and phasors, and (ii) Laplace-transform analysis.

A causal LTI system can be described by an nth-order ordinary differential equation (ODE)

\[ a_n \frac{d^n y(t)}{dt^n} + a_{n-1} \frac{d^{n-1} y(t)}{dt^{n-1}} + \dots + a_1 \frac{dy(t)}{dt} + a_0 y(t) = b_m \frac{d^m u(t)}{dt^m} + \dots + b_1 \frac{du(t)}{dt} + b_0 u(t), \]

with constant real coefficients \( a_i, b_j \). Taking Laplace transforms (assuming zero initial conditions), we obtain the transfer function representation

\[ G(s) = \frac{Y(s)}{U(s)} = \frac{b_m s^m + \dots + b_1 s + b_0}{a_n s^n + \dots + a_1 s + a_0}. \]

A fundamental property of LTI systems is that complex exponentials \( e^{st} \) are eigenfunctions. Consider an input \( u(t) = U e^{st} \) for complex \( s \). We seek a particular solution of the form \( y_p(t) = Y e^{st} \), where \( Y \) is a complex constant to be determined.

Substituting into the ODE and dividing by \( e^{st} \) yields

\[ \big(a_n s^n + a_{n-1} s^{n-1} + \dots + a_0\big) Y = \big(b_m s^m + \dots + b_0\big) U. \]

Hence,

\[ Y = G(s)\,U, \quad \text{so that} \quad y_p(t) = G(s)\,U\,e^{st}. \]

Therefore, for any exponential input \( e^{st} \), the forced response of the LTI system has the same exponential factor multiplied by the complex gain \( G(s) \) (plus homogeneous transients corresponding to the system poles).

Real sinusoids can be represented as the real part of complex exponentials. For example,

\[ \cos(\omega t + \varphi_u) = \Re\!\big\{ e^{j(\omega t + \varphi_u)} \big\} = \Re\!\big\{ e^{j\varphi_u} e^{j\omega t} \big\}. \]

Consider the complex input \( u_c(t) = U e^{j(\omega t + \varphi_u)} \). By the exponential eigenfunction property, a particular response is

\[ y_c(t) = G(j\omega)\,U\,e^{j(\omega t + \varphi_u)}. \]

Let us write \( G(j\omega) \) in polar form \( G(j\omega) = \left| G(j\omega) \right| e^{j\varphi_G(\omega)} \). Then

\[ y_c(t) = \left| G(j\omega) \right| U e^{j\big( \omega t + \varphi_u + \varphi_G(\omega) \big)}. \]

The physical output is the real part of \( y_c(t) \). Therefore,

\[ y_{\text{ss}}(t) = \Re\{ y_c(t) \} = \left| G(j\omega) \right| U \cos\!\big( \omega t + \varphi_u + \varphi_G(\omega) \big), \]

which proves the sinusoidal steady-state formula of Section 1 under the assumption that the system is asymptotically stable (so transients decay).

This calculation shows that \( G(j\omega) \) acts as a complex gain: its magnitude scales the amplitude, and its argument shifts the phase.

In the phasor representation, a sinusoid \( u(t) = U \cos(\omega t + \varphi_u) \) at fixed frequency \( \omega \) is encoded as a complex number (phasor) \( \underline{U} = U e^{j\varphi_u} \). The system, for that fixed frequency, is represented by the complex gain \( G(j\omega) \). The output phasor is then

\[ \underline{Y} = G(j\omega)\,\underline{U}. \]

Writing \( \underline{Y} = Y_{\text{amp}} e^{j\varphi_y} \), we have

\[ Y_{\text{amp}} = \left| G(j\omega) \right| U, \quad \varphi_y = \varphi_u + \arg G(j\omega). \]

Converting back to the time domain:

\[ y_{\text{ss}}(t) = Y_{\text{amp}} \cos\!\big( \omega t + \varphi_y \big) = \left| G(j\omega) \right| U \cos\!\big( \omega t + \varphi_u + \arg G(j\omega) \big). \]

This viewpoint is heavily used in electrical engineering and robotics, for instance when analyzing the effect of sinusoidal disturbances (e.g., periodic torque ripple, sensor vibrations) on robot joints modeled as LTI systems.

flowchart TD
  IN["Input phasor: U angle phi_u"] --> GAIN["Complex gain: G(jw) = |G| angle phi_G"]
  GAIN --> OUT["Output phasor: |G| U angle (phi_u + phi_G)"]
  OUT --> TIME["Time signal y_ss(t) = |G(jw)| U cos(omega t + phi_u + phi_G)"]
        

Consider the standard first-order transfer function (studied earlier in the time domain)

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

Evaluating at \( s = j\omega \) gives

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

The magnitude and phase are

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

For an input \( u(t) = U \cos(\omega t) \) (zero input phase for simplicity), the steady-state output is

\[ y_{\text{ss}}(t) = \frac{K U}{\sqrt{1 + (\omega \tau)^2}} \cos\!\big( \omega t - \arctan(\omega \tau) \big). \]

As \( \omega \tau \rightarrow 0 \) (very low frequency), the amplitude ratio \( \left| G(j\omega) \right| \rightarrow K \) and the phase shift tends to zero. As \( \omega \tau \rightarrow \infty \) (very high frequency), the amplitude ratio tends to \( 0 \) and the phase lag approaches \( -\tfrac{\pi}{2} \).

This frequency-dependent behavior is the starting point for understanding filtering, bandwidth, and disturbance rejection in robot actuators and servo drives.

We now illustrate sinusoidal input response using Python. We use the python-control library (commonly used in robotics and mechatronics) to simulate a first-order actuator model, which could represent a simplified robot joint velocity loop.

import numpy as np
import control  # pip install control
import matplotlib.pyplot as plt

# Physical parameters (example: velocity actuator in a robot joint)
K = 5.0      # static gain
tau = 0.1    # time constant [s]

# Transfer function G(s) = K / (tau s + 1)
num = [K]
den = [tau, 1.0]
G = control.tf(num, den)

# Sinusoidal input parameters
U = 1.0          # input amplitude
omega = 10.0     # rad/s
t_final = 2.0
n_points = 2000
t = np.linspace(0.0, t_final, n_points)
u = U * np.sin(omega * t)

# Simulate time response to sinusoidal input
t_out, y_out, x_out = control.forced_response(G, T=t, U=u)

# Theoretical steady-state amplitude and phase
Gjw = control.evalfr(G, 1j * omega)
mag = abs(Gjw)
phase = np.angle(Gjw)  # radians

print("Complex gain G(jw) =", Gjw)
print("Amplitude ratio |G(jw)| =", mag)
print("Phase shift arg(G(jw)) [rad] =", phase)

# Optional plotting (remove or adapt for headless servers)
plt.figure()
plt.plot(t_out, u, label="input u(t)")
plt.plot(t_out, y_out, label="output y(t)")
plt.xlabel("t [s]")
plt.legend()
plt.grid(True)
plt.show()
      

In a robotics context, the same methodology applies to more complex LTI models extracted from rigid-body dynamics (e.g., after linearization around a trajectory using a robotics toolbox), where G represents the transfer function from torque commands to joint position or velocity.

For embedded control in robotic systems (e.g., with ROS control loops), we often simulate or implement LTI dynamics numerically. Below is a minimal C++ example that integrates a first-order model \( \tau \dot{y}(t) + y(t) = K u(t) \) driven by a sinusoidal input \( u(t) = U \sin(\omega t) \) using forward Euler integration.

#include <iostream>
#include <cmath>

int main() {
    double K = 5.0;
    double tau = 0.1;
    double U = 1.0;
    double omega = 10.0;

    double dt = 1e-4;
    double t_final = 2.0;
    int steps = static_cast<int>(t_final / dt);

    double y = 0.0;   // output state
    double t = 0.0;

    for (int k = 0; k < steps; ++k) {
        double u = U * std::sin(omega * t);
        double dy = (-y + K * u) / tau;  // tau * dy/dt + y = K u

        // Forward Euler update
        y += dt * dy;
        t += dt;

        if (k % 1000 == 0) {
            std::cout << t << " " << u << " " << y << std::endl;
        }
    }

    return 0;
}
      

In more advanced robotic software, one would typically represent the LTI system in state-space form using a linear algebra library such as Eigen, and integrate it inside a real-time control loop (e.g., within a ROS controller) while still exploiting the theoretical insight that sinusoidal disturbances are scaled and phase-shifted according to \( G(j\omega) \).

Java is sometimes used in robotics (e.g., educational robots, FIRST robotics with WPILib). Below is a basic Java implementation of the same first-order plant, using discrete-time Euler integration. More sophisticated implementations could use matrix libraries such as Apache Commons Math for multivariable models.

public class SinusoidalLtiSimulation {

    public static void main(String[] args) {
        double K = 5.0;
        double tau = 0.1;
        double U = 1.0;
        double omega = 10.0;

        double dt = 1e-4;
        double tFinal = 2.0;
        int steps = (int) (tFinal / dt);

        double y = 0.0;
        double t = 0.0;

        for (int k = 0; k < steps; k++) {
            double u = U * Math.sin(omega * t);
            double dy = (-y + K * u) / tau;

            y += dt * dy;
            t += dt;

            if (k % 1000 == 0) {
                System.out.println(t + " " + u + " " + y);
            }
        }
    }
}
      

In robot control frameworks, this type of simulation supports design and verification of control algorithms that must tolerate sinusoidal disturbances (e.g., periodic load variations).

MATLAB and Simulink are standard tools in control and robotics. Using the Control System Toolbox, one can directly build \( G(s) \) and simulate its response to sinusoidal inputs.

% Parameters for first-order LTI system
K = 5;
tau = 0.1;

s = tf('s');
G = K / (tau * s + 1);

% Sinusoidal input
U = 1;
omega = 10;          % rad/s
t_final = 2;
t = linspace(0, t_final, 2000);
u = U * sin(omega * t);

% Time response to sinusoidal input
[y, t_out] = lsim(G, u, t);

% Theoretical complex gain at frequency omega
Gjw = freqresp(G, omega);
mag = abs(Gjw);
phase = angle(Gjw);  % radians

fprintf('Amplitude ratio |G(jw)| = %g\n', mag);
fprintf('Phase shift arg(G(jw)) = %g rad\n', phase);

figure;
plot(t_out, u, 'LineWidth', 1.2); hold on;
plot(t_out, y, 'LineWidth', 1.2);
xlabel('t [s]');
ylabel('Signals');
legend('u(t)', 'y(t)');
grid on;
title('Sinusoidal Input Response of First-Order System');
      

In Simulink, the same system can be built using a Sine Wave block feeding a Transfer Fcn block with numerator [K] and denominator [tau 1]. For robotics applications, Simulink models derived from the Robotics System Toolbox (e.g., linearized joint dynamics) can be driven by sinusoidal test inputs to study frequency-dependent behavior.

Wolfram Mathematica provides symbolic and numeric tools for analyzing sinusoidal responses. Below we solve a first-order ODE and extract its steady-state behavior.

(* Parameters *)
K = 5;
tau = 1/10;
U = 1;
omega = 10;

(* Differential equation: tau y'(t) + y(t) == K u(t) *)
u[t_] := U*Sin[omega*t];
eq = tau*y'[t] + y[t] == K*u[t];
ic = y[0] == 0;

(* Solve ODE *)
sol = DSolve[{eq, ic}, y, t][[1]];
ySol[t_] = y[t] /. sol // Simplify

(* Extract long-time behavior (steady-state) by assuming t is large *)
ySS[t_] = ComplexExpand[
  Assuming[t > 0, TrigToExp[ySol[t]]]
] // Simplify

(* Alternatively: use transfer function approach *)
sys = TransferFunctionModel[K/(tau*s + 1), s];
Gjw = TransferFunctionExpand[sys /. s -> I*omega]
mag = Abs[Gjw]
phase = Arg[Gjw]

(* Numeric plot *)
Plot[{u[t], ySol[t]}, {t, 0, 2},
 PlotLegends -> {"u(t)", "y(t)"},
 AxesLabel -> {"t", "Signal"}]
      

The symbolic approach makes it easy to see the transient term (typically containing exponentials) and the steady-state sinusoidal term, confirming the general theory derived earlier.

Problem 1 (First-Order Sinusoidal Response): Consider the first-order system \( G(s) = \dfrac{K}{\tau s + 1} \) with \( K > 0 \), \( \tau > 0 \). For input \( u(t) = U \cos(\omega t) \), derive the steady-state output \( y_{\text{ss}}(t) \).

Solution:

We use the complex exponential method. Consider the complex input \( u_c(t) = U e^{j\omega t} \). The corresponding particular response is

\[ y_c(t) = G(j\omega)\,U\,e^{j\omega t} = \frac{K}{1 + j\omega \tau} U e^{j\omega t}. \]

Write \( 1 + j\omega \tau = \sqrt{1 + (\omega \tau)^2} \, e^{j\arctan(\omega \tau)} \), so

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

Thus

\[ y_c(t) = \frac{K U}{\sqrt{1 + (\omega \tau)^2}} e^{j\big( \omega t - \arctan(\omega \tau) \big)}. \]

Taking the real part (physical signal) gives

\[ y_{\text{ss}}(t) = \frac{K U}{\sqrt{1 + (\omega \tau)^2}} \cos\!\big( \omega t - \arctan(\omega \tau) \big). \]

This is the sinusoidal steady-state response of the first-order system.

Problem 2 (Second-Order Standard Form): Consider the second-order system

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

For input \( u(t) = U \cos(\omega t) \), derive the magnitude and phase of the steady-state output.

Solution:

Evaluate \( G(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^2 + j(2\zeta \omega_n \omega) + \omega_n^2}. \]

Group real and imaginary parts in the denominator:

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

Therefore, the magnitude is

\[ \left| G(j\omega) \right| = \frac{\omega_n^2}{ \sqrt{ \big(\omega_n^2 - \omega^2\big)^2 + \big(2\zeta \omega_n \omega\big)^2 } }, \]

and the phase is

\[ \arg G(j\omega) = -\arctan\!\left( \frac{2\zeta \omega_n \omega}{\omega_n^2 - \omega^2} \right). \]

Thus the steady-state output is

\[ y_{\text{ss}}(t) = \left| G(j\omega) \right| U \cos\!\big( \omega t + \arg G(j\omega) \big). \]

The magnitude expression reveals a resonant peak when \( \omega \approx \omega_n \) and \( \zeta \) is small, a phenomenon important in vibration and robot structure design.

Problem 3 (Resonance When Input Matches an Imaginary-Axis Pole): Consider the undamped oscillator

\[ G(s) = \frac{1}{s^2 + \omega_0^2}, \]

which has poles at \( s = \pm j\omega_0 \). Let \( u(t) = U \sin(\omega_0 t) \). Show that the output is unbounded and grows linearly with time.

Solution:

The time-domain equation is

\[ \ddot{y}(t) + \omega_0^2 y(t) = U \sin(\omega_0 t). \]

The homogeneous solution is \( y_h(t) = C_1 \cos(\omega_0 t) + C_2 \sin(\omega_0 t) \). For the particular solution we try \( y_p(t) = t\big( A \cos(\omega_0 t) + B \sin(\omega_0 t) \big) \), because the forcing frequency coincides with the natural frequency.

Differentiating and substituting into the ODE (details omitted for brevity) yields \( A = 0 \) and \( B = \dfrac{U}{2\omega_0} \), so

\[ y_p(t) = \frac{U}{2\omega_0} t \sin(\omega_0 t). \]

Therefore,

\[ y(t) = y_h(t) + y_p(t) = C_1 \cos(\omega_0 t) + C_2 \sin(\omega_0 t) + \frac{U}{2\omega_0} t \sin(\omega_0 t), \]

which is unbounded as \( t → \infty \). This shows that the sinusoidal steady-state framework requires an asymptotically stable system (no poles on the imaginary axis).

Problem 4 (Numerical Amplitude Ratio for a Robotic Joint Actuator): A simplified velocity actuator of a robot joint has first-order model \( G(s) = \dfrac{100}{0.02 s + 1} \). A periodic disturbance torque is modeled as \( u(t) = 0.5 \cos(50 t) \). Compute the steady-state velocity oscillation amplitude \( Y_{\text{amp}} \).

Solution:

We have \( K = 100 \), \( \tau = 0.02 \), \( U = 0.5 \), \( \omega = 50 \) rad/s. The magnitude is

\[ \left| G(j\omega) \right| = \frac{K}{\sqrt{1 + (\omega \tau)^2}} = \frac{100}{\sqrt{1 + (50 \cdot 0.02)^2}} = \frac{100}{\sqrt{1 + 1^2}} = \frac{100}{\sqrt{2}}. \]

Therefore the output amplitude is

\[ Y_{\text{amp}} = \left| G(j\omega) \right| U = \frac{100}{\sqrt{2}} \cdot 0.5 = \frac{50}{\sqrt{2}} \approx 35.36. \]

The steady-state velocity oscillation is approximately \( 35.36 \cos(50 t + \varphi) \) with \( \varphi = -\arctan(\omega \tau) = -\arctan(1) = -\tfrac{\pi}{4} \).

In this lesson we showed that for an asymptotically stable LTI system, sinusoids are eigenfunctions: a sinusoidal input at frequency \( \omega \) yields, in steady state, a sinusoidal output at the same frequency, with amplitude and phase determined by \( G(j\omega) \). Using exponential and Laplace-transform arguments, we derived explicit formulas for first- and second-order systems and discussed resonance phenomena when inputs align with natural frequencies. Practical implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica illustrated how to simulate and verify sinusoidal responses, providing a bridge from abstract theory to robotic control applications. Subsequent lessons will generalize these ideas into full frequency-response analysis and Bode/Nyquist design tools.