All three example systems can be described by the same second-order matrix form introduced informally in earlier lessons for coupled translational motions:

\[ \mathbf{M}\,\ddot{\mathbf{q}}(t) + \mathbf{C}\,\dot{\mathbf{q}}(t) + \mathbf{K}\,\mathbf{q}(t) = \mathbf{f}(t), \]

where \( \mathbf{q}(t) \) collects the translational coordinates (displacements), \( \mathbf{M} \) is the mass matrix, \( \mathbf{C} \) the viscous damping matrix, and \( \mathbf{K} \) the stiffness matrix. The vector \( \mathbf{f}(t) \) represents external forces or equivalent base-excitation terms. For single-degree-of-freedom (SDOF) systems these reduce to scalar coefficients.

In this lesson we:

• Derive \( \mathbf{M},\mathbf{C},\mathbf{K} \) for each example.
• Express the ODEs with clear physical interpretation of each term.
• Show how to rewrite second-order equations as equivalent first-order systems, preparing the ground for later state-space modeling.
• Implement direct time-domain simulations using basic numerical integration.

flowchart TD
  A["Physical description (vehicle, building, actuator)"]
    --> B["Choose generalized coordinates q(t)"]
  B --> C["Draw free-body diagrams for each mass"]
  C --> D["Sum of forces = mass * acceleration"]
  D --> E["Identify equivalent C (damping) and K (stiffness)"]
  E --> F["Assemble M q''(t) + C q'(t) + K q(t) = f(t)"]
  F --> G["Check units, limiting cases, parameter ranges"]
        

A quarter-car model captures the vertical dynamics of one corner of a vehicle:

• \( m_s \): sprung mass (portion of vehicle body).
• \( m_u \): unsprung mass (wheel, hub, part of suspension).
• \( k_s, c_s \): suspension spring and damper between body and wheel.
• \( k_t \): tire stiffness between wheel and road.
• \( x_s(t) \), \( x_u(t) \): vertical displacements of sprung and unsprung mass (inertial frame).
• \( x_r(t) \): prescribed road displacement input.

We choose upward displacements as positive. The forces in the suspension element are:

\[ F_s(t) = k_s\big(x_s(t) - x_u(t)\big) + c_s\big(\dot{x}_s(t) - \dot{x}_u(t)\big), \]

acting downward on the sprung mass and upward on the unsprung mass when \( x_s(t) > x_u(t) \). The tire force is

\[ F_t(t) = k_t\big(x_u(t) - x_r(t)\big), \]

acting upward on the unsprung mass when the wheel is compressed relative to the road.

Equations of motion. Applying Newton's law to each mass:

\[ m_s \ddot{x}_s(t) = -F_s(t), \quad m_u \ddot{x}_u(t) = F_s(t) - F_t(t). \]

Substituting the expressions for \( F_s(t) \) and \( F_t(t) \) and rearranging gives:

\[ m_s \ddot{x}_s(t) + c_s\big(\dot{x}_s(t) - \dot{x}_u(t)\big) + k_s\big(x_s(t) - x_u(t)\big) = 0, \]

\[ m_u \ddot{x}_u(t) + c_s\big(\dot{x}_u(t) - \dot{x}_s(t)\big) + k_s\big(x_u(t) - x_s(t)\big) + k_t\big(x_u(t) - x_r(t)\big) = 0. \]

These two coupled second-order ODEs fully describe the vertical quarter-car dynamics in the chosen coordinates.

Matrix form. Define

\[ \mathbf{q}(t) = \begin{bmatrix} x_s(t)\\[4pt] x_u(t) \end{bmatrix},\quad \mathbf{M} = \begin{bmatrix} m_s & 0\\[4pt] 0 & m_u \end{bmatrix}, \]

\[ \mathbf{C} = \begin{bmatrix} c_s & -c_s\\[4pt] -c_s & c_s \end{bmatrix},\quad \mathbf{K} = \begin{bmatrix} k_s & -k_s\\[4pt] -k_s & k_s + k_t \end{bmatrix}, \]

\[ \mathbf{f}(t) = \begin{bmatrix} 0\\[4pt] k_t x_r(t) \end{bmatrix}. \]

Then the quarter-car model can be written compactly as

\[ \mathbf{M}\,\ddot{\mathbf{q}}(t) + \mathbf{C}\,\dot{\mathbf{q}}(t) + \mathbf{K}\,\mathbf{q}(t) = \mathbf{f}(t). \]

This structure mirrors the general framework and will reappear in the actuator–load system below.

For some control and comfort studies, one is mainly interested in the motion of the vehicle body \( x_s(t) \) at relatively low frequencies. Under simplifying assumptions, the quarter-car can be approximated by a single-degree-of-freedom base-excited oscillator.

A common simplification is to suppose that the unsprung mass follows the road closely: \( x_u(t) \approx x_r(t) \). Then the sprung mass equation reduces to

\[ m_s \ddot{x}_s(t) + c_s\big(\dot{x}_s(t) - \dot{x}_r(t)\big) + k_s\big(x_s(t) - x_r(t)\big) = 0. \]

Introducing the relative coordinate \( z(t) = x_s(t) - x_r(t) \) gives

\[ \dot{z}(t) = \dot{x}_s(t) - \dot{x}_r(t),\quad \ddot{z}(t) = \ddot{x}_s(t) - \ddot{x}_r(t). \]

Substituting into the equation for \( x_s(t) \) leads to

\[ m_s\big(\ddot{z}(t) + \ddot{x}_r(t)\big) + c_s \dot{z}(t) + k_s z(t) = 0, \]

hence

\[ m_s \ddot{z}(t) + c_s \dot{z}(t) + k_s z(t) = - m_s \ddot{x}_r(t). \]

This is exactly the form of an SDOF mass–spring–damper system with \( z(t) \) as generalized coordinate and an effective inertial excitation \( -m_s \ddot{x}_r(t) \).

Such reduced models are very useful in early control design (e.g., designing a vertical ride controller) and are much cheaper to simulate analytically and numerically than full multi-degree-of-freedom vehicle models.

A single floor of a building can often be idealized as a lumped mass supported by vertical columns that provide effective stiffness and damping. Let \( m \) denote the floor mass, \( k \) the equivalent vertical stiffness, and \( c \) the equivalent viscous damping representing structural damping and energy dissipation mechanisms.

If \( x(t) \) is the vertical displacement of the floor relative to a fixed reference and \( f(t) \) is an externally applied vertical force (for example, from a heavy machine), then the governing equation is

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

Dividing by \( m \) yields the normalized form

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

If the base of the building moves during an earthquake, an analogous base-excited model can be obtained using the same relative-coordinate technique as for the quarter-car SDOF approximation. This example demonstrates that building floor dynamics are structurally identical to vehicle body bounce dynamics after suitable identification of parameters.

For numerical work it is convenient to write this second-order ODE as a first-order system. Introduce

\[ x_1(t) = x(t),\quad x_2(t) = \dot{x}(t), \]

which gives

\[ \dot{x}_1(t) = x_2(t),\quad \dot{x}_2(t) = -\frac{c}{m} x_2(t) - \frac{k}{m} x_1(t) + \frac{1}{m} f(t). \]

This first-order representation is directly implementable in numerical integrators, and the same idea will be used for vehicle and actuator–load models in the code sections.

Translational actuators (e.g., linear electric actuators, hydraulic cylinders) are often required to move a mechanical load. Depending on the mechanical coupling, different dynamic structures arise. We consider a simple but widely used model: two masses coupled by a spring–damper pair, with the actuator force applied to the first mass.

• \( m_a \): equivalent mass of actuator and moving internal parts.
• \( m_L \): load mass.
• \( k_c, c_c \): coupling stiffness and damping between actuator and load.
• \( b_a, b_L \): viscous friction at actuator and load interfaces.
• \( x_a(t), x_L(t) \): actuator and load displacements.
• \( u(t) \): control input, modeled here as an external force applied to the actuator mass.

The coupling force is

\[ F_c(t) = k_c\big(x_a(t) - x_L(t)\big) + c_c\big(\dot{x}_a(t) - \dot{x}_L(t)\big), \]

acting on the actuator in the negative direction of \( x_a(t) \) when \( x_a(t) > x_L(t) \), and on the load in the positive direction.

Equations of motion. Newton's law for each mass gives

\[ m_a \ddot{x}_a(t) = u(t) - F_c(t) - b_a \dot{x}_a(t), \]

\[ m_L \ddot{x}_L(t) = F_c(t) - b_L \dot{x}_L(t). \]

Substituting \( F_c(t) \) and rearranging yields the coupled ODE system

\[ m_a \ddot{x}_a(t) + (b_a + c_c)\dot{x}_a(t) - c_c \dot{x}_L(t) + (k_c) x_a(t) - k_c x_L(t) = u(t), \]

\[ m_L \ddot{x}_L(t) - c_c \dot{x}_a(t) + (b_L + c_c)\dot{x}_L(t) - k_c x_a(t) + k_c x_L(t) = 0. \]

In matrix form with

\[ \mathbf{q}(t) = \begin{bmatrix} x_a(t)\\[4pt] x_L(t) \end{bmatrix},\quad \mathbf{M} = \begin{bmatrix} m_a & 0\\[4pt] 0 & m_L \end{bmatrix}, \]

\[ \mathbf{C} = \begin{bmatrix} b_a + c_c & -c_c\\[4pt] -c_c & b_L + c_c \end{bmatrix},\quad \mathbf{K} = \begin{bmatrix} k_c & -k_c\\[4pt] -k_c & k_c \end{bmatrix}, \]

\[ \mathbf{B} = \begin{bmatrix} 1\\[4pt] 0 \end{bmatrix}, \quad \mathbf{f}(t) = \mathbf{B}\,u(t), \]

we obtain again \( \mathbf{M}\ddot{\mathbf{q}}(t) + \mathbf{C}\dot{\mathbf{q}}(t) + \mathbf{K}\mathbf{q}(t) = \mathbf{f}(t) \), structurally identical to the quarter-car formulation. This structural similarity is important for later controller design.

flowchart LR
  Q["2-DOF mass-spring-damper chain"]
    --> V["Vehicle suspension: (x_s, x_u) with road input x_r(t)"]
  Q --> B["Building SDOF: special case with single mass and base motion"]
  Q --> A["Actuator-load: (x_a, x_L) with input force u(t)"]
        

We now simulate the quarter-car model using Python with numpy, scipy.integrate.solve_ivp, and matplotlib. The second-order system is rewritten as a first-order system in the state \( \mathbf{z}(t) = [x_s(t), \dot{x}_s(t), x_u(t), \dot{x}_u(t)]^{\top} \).

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Quarter-car parameters (typical orders of magnitude)
m_s = 300.0    # sprung mass [kg]
m_u = 40.0     # unsprung mass [kg]
k_s = 15000.0  # suspension stiffness [N/m]
c_s = 1500.0   # suspension damping [N*s/m]
k_t = 200000.0 # tire stiffness [N/m]

def road_profile(t):
    """
    Simple road input: a bump of height h for a short duration.
    """
    h = 0.05  # 5 cm bump
    if 0.0 <= t and t <= 0.2:
        return h
    return 0.0

def quarter_car_ode(t, z):
    """
    State vector z = [x_s, v_s, x_u, v_u].
    """
    x_s, v_s, x_u, v_u = z
    x_r = road_profile(t)

    # Suspension and tire forces
    F_s = k_s * (x_s - x_u) + c_s * (v_s - v_u)
    F_t = k_t * (x_u - x_r)

    # Accelerations from Newton's laws
    a_s = -F_s / m_s
    a_u = (F_s - F_t) / m_u

    return [v_s, a_s, v_u, a_u]

# Initial conditions: all displacements and velocities zero
z0 = [0.0, 0.0, 0.0, 0.0]

t_span = (0.0, 3.0)
t_eval = np.linspace(t_span[0], t_span[1], 2000)

sol = solve_ivp(quarter_car_ode, t_span, z0, t_eval=t_eval)

x_s = sol.y[0, :]
x_u = sol.y[2, :]
x_r = np.array([road_profile(t) for t in t_eval])

plt.figure()
plt.plot(t_eval, x_s, label="sprung mass x_s(t)")
plt.plot(t_eval, x_u, label="unsprung mass x_u(t)")
plt.plot(t_eval, x_r, "--", label="road x_r(t)")
plt.xlabel("time [s]")
plt.ylabel("vertical displacement [m]")
plt.legend()
plt.grid(True)
plt.show()
      

This script demonstrates the direct translation of the modeled equations into code. More advanced Python packages for system dynamics and control (e.g., higher-level modeling interfaces) can be used later, but here we rely only on ODE solvers and basic arrays.

In C++, we implement a simple explicit Euler integrator for the building floor SDOF model \( m \ddot{x}(t) + c \dot{x}(t) + k x(t) = f(t) \). Euler's method for a first-order system \( \dot{\mathbf{x}}(t) = \mathbf{g}(\mathbf{x}(t),t) \) with time step \( h \) reads \( \mathbf{x}_{k+1} = \mathbf{x}_k + h\,\mathbf{g}(\mathbf{x}_k,t_k) \).

#include <iostream>
#include <vector>

int main() {
    // Physical parameters
    double m = 1000.0;     // floor mass [kg]
    double c = 1500.0;     // damping [N*s/m]
    double k = 2.0e5;      // stiffness [N/m]

    // Time discretization
    double t0 = 0.0;
    double tf = 10.0;
    int n_steps = 10000;
    double h = (tf - t0) / static_cast<double>(n_steps);

    std::vector<double> t(n_steps + 1);
    std::vector<double> x(n_steps + 1);  // displacement
    std::vector<double> v(n_steps + 1);  // velocity

    // Initial conditions
    x[0] = 0.0;
    v[0] = 0.0;
    t[0] = t0;

    auto forcing = [](double time) {
        // Example: half-sine pulse load
        double T = 1.0;
        double F0 = 1000.0;
        if (time >= 0.0 && time <= T) {
            return F0 * std::sin(3.141592653589793 * time / T);
        }
        return 0.0;
    };

    // Time stepping (explicit Euler)
    for (int k_step = 0; k_step < n_steps; ++k_step) {
        t[k_step + 1] = t[k_step] + h;
        double f = forcing(t[k_step]);
        double a = (f - c * v[k_step] - k * x[k_step]) / m;

        // Euler update
        v[k_step + 1] = v[k_step] + h * a;
        x[k_step + 1] = x[k_step] + h * v[k_step];
    }

    // Print a subset of samples
    for (int k_step = 0; k_step <= n_steps; k_step += 1000) {
        std::cout << "t = " << t[k_step]
                  << " s, x = " << x[k_step] << " m" << std::endl;
    }

    return 0;
}
      

For more sophisticated and accurate C++ simulations one would typically use dedicated numerical libraries such as Boost.odeint, Eigen-based integrators, or wrappers to high-order Runge–Kutta methods.

The actuator–load 2-DOF system can be simulated similarly using Java. We again use a simple explicit Euler method for illustrative purposes. In practice one might employ numerical ODE libraries such as Apache Commons Math, but the basic structure is the same.

public class ActuatorLoadSimulation {

    // Physical parameters
    static double mA = 5.0;       // actuator mass [kg]
    static double mL = 20.0;      // load mass [kg]
    static double kC = 1.0e4;     // coupling stiffness [N/m]
    static double cC = 200.0;     // coupling damping [N*s/m]
    static double bA = 50.0;      // actuator friction [N*s/m]
    static double bL = 20.0;      // load friction [N*s/m]

    static double inputForce(double t) {
        // Step-like input force after t = 0.1 s
        if (t >= 0.1) {
            return 500.0;
        }
        return 0.0;
    }

    public static void main(String[] args) {
        double t0 = 0.0;
        double tf = 2.0;
        int nSteps = 4000;
        double h = (tf - t0) / (double) nSteps;

        double t = t0;
        double xA = 0.0;  // actuator displacement
        double vA = 0.0;  // actuator velocity
        double xL = 0.0;  // load displacement
        double vL = 0.0;  // load velocity

        for (int k = 0; k <= nSteps; ++k) {
            double u = inputForce(t);

            // Coupling force
            double Fc = kC * (xA - xL) + cC * (vA - vL);

            // Accelerations from Newton's laws
            double aA = (u - Fc - bA * vA) / mA;
            double aL = (Fc - bL * vL) / mL;

            // Explicit Euler update
            xA = xA + h * vA;
            vA = vA + h * aA;
            xL = xL + h * vL;
            vL = vL + h * aL;

            if (k % 400 == 0) {
                System.out.printf(
                    "t = %.3f s, xA = %.4f m, xL = %.4f m%n",
                    t, xA, xL
                );
            }

            t += h;
        }
    }
}
      

The mapping from the continuous-time model \( \mathbf{M}\ddot{\mathbf{q}}(t) + \mathbf{C}\dot{\mathbf{q}}(t) + \mathbf{K}\mathbf{q}(t) = \mathbf{B}u(t) \) to the discrete updates in the code is direct once the system is written in first-order form.

MATLAB offers built-in ODE solvers such as ode45 that are well suited to simulating the models derived above. The following script simulates the building floor model \( m \ddot{x}(t) + c \dot{x}(t) + k x(t) = f(t) \).

function building_floor_demo
    m = 1000;      % mass [kg]
    c = 1500;      % damping [N*s/m]
    k = 2e5;       % stiffness [N/m]

    tspan = [0 10];
    x0 = [0; 0];   % [x(0); xdot(0)]

    [t, x] = ode45(@(t, x) floor_ode(t, x, m, c, k), tspan, x0);

    figure;
    plot(t, x(:, 1), 'LineWidth', 1.5);
    xlabel('time [s]');
    ylabel('displacement x(t) [m]');
    grid on;
end

function dx = floor_ode(t, x, m, c, k)
    % x(1) = displacement, x(2) = velocity
    F = external_force(t);
    dx = zeros(2, 1);
    dx(1) = x(2);
    dx(2) = (F - c * x(2) - k * x(1)) / m;
end

function F = external_force(t)
    % Half-sine pulse, similar to the C++ example
    T = 1.0;
    F0 = 1000.0;
    if t >= 0 && t <= T
        F = F0 * sin(pi * t / T);
    else
        F = 0.0;
    end
end
      

A Simulink diagram implementing the same model would use the following blocks:

1. A Sum block to compute (F - c*xdot - k*x)/m.
2. Two cascaded Integrator blocks to obtain x from the acceleration.
3. A MATLAB Function or Signal Builder block to generate F(t).
4. Scope blocks to visualize displacement and velocity.

The structure of the Simulink model mirrors the first-order ODE form used in floor_ode.

Wolfram Mathematica provides symbolic and numeric capabilities for ODEs. The quarter-car model can be implemented symbolically and then simulated numerically using NDSolve.

(* Quarter-car parameters *)
ms = 300.0;
mu = 40.0;
ks = 15000.0;
cs = 1500.0;
kt = 200000.0;

xr[t_] := Piecewise[{ {0.05, 0 <= t <= 0.2} }, 0]; (* road bump *)

(* Unknown functions xs[t], xu[t] *)
eqs = {
  ms*xs''[t] + cs*(xs'[t] - xu'[t]) + ks*(xs[t] - xu[t]) == 0,
  mu*xu''[t] + cs*(xu'[t] - xs'[t]) + ks*(xu[t] - xs[t]) +
    kt*(xu[t] - xr[t]) == 0,
  xs[0] == 0, xs'[0] == 0,
  xu[0] == 0, xu'[0] == 0
};

sol = NDSolve[eqs, {xs, xu}, {t, 0, 3}];

Plot[
  Evaluate[{xs[t], xu[t], xr[t]} /. sol],
  {t, 0, 3},
  PlotLegends -> {"xs(t)", "xu(t)", "xr(t)"},
  AxesLabel -> {"t [s]", "displacement [m]"}
]
      

Here, the model equations are written almost exactly as in the analytical derivation. Mathematica's symbolic features can later be used to obtain analytic solutions for simplified cases (e.g., undamped systems or constant inputs).

Problem 1 (Matrix Form of Quarter-Car Model). Starting from the scalar equations of motion for the quarter-car model, derive the mass, damping, and stiffness matrices \( \mathbf{M},\mathbf{C},\mathbf{K} \) and the input vector \( \mathbf{f}(t) \) in \( \mathbf{M}\ddot{\mathbf{q}}(t) + \mathbf{C}\dot{\mathbf{q}}(t) + \mathbf{K}\mathbf{q}(t) = \mathbf{f}(t) \).

Solution. From Section 2 we have

\[ m_s \ddot{x}_s + c_s(\dot{x}_s - \dot{x}_u) + k_s(x_s - x_u) = 0, \]

\[ m_u \ddot{x}_u + c_s(\dot{x}_u - \dot{x}_s) + k_s(x_u - x_s) + k_t(x_u - x_r) = 0. \]

Collect coefficients of \( \ddot{x}_s,\ddot{x}_u,\dot{x}_s,\dot{x}_u,x_s,x_u \). Defining \( \mathbf{q}(t) = [x_s(t), x_u(t)]^{\top} \), the mass matrix has diagonal entries \( m_s,m_u \), giving

\[ \mathbf{M} = \begin{bmatrix} m_s & 0\\[4pt] 0 & m_u \end{bmatrix}. \]

The damping terms are \( c_s(\dot{x}_s - \dot{x}_u) \) in the first equation and \( c_s(\dot{x}_u - \dot{x}_s) \) in the second, which yields

\[ \mathbf{C} = \begin{bmatrix} c_s & -c_s\\[4pt] -c_s & c_s \end{bmatrix}. \]

The stiffness terms have the same structure with \( k_s \), and additionally \( k_t x_u \) in the second equation. Thus

\[ \mathbf{K} = \begin{bmatrix} k_s & -k_s\\[4pt] -k_s & k_s + k_t \end{bmatrix}. \]

The term involving \( x_r(t) \) appears only in the second equation with coefficient \( -k_t x_r(t) \), which can be moved to the right-hand side, giving

\[ \mathbf{f}(t) = \begin{bmatrix} 0\\[4pt] k_t x_r(t) \end{bmatrix}. \]

This recovers the matrices in Section 2 and shows that the quarter-car model fits the standard second-order MCK structure.

Problem 2 (Free Vibration of the Building Floor). Consider the undamped building floor model with no external force: \( m \ddot{x}(t) + k x(t) = 0 \). Solve this ODE for the initial conditions \( x(0) = x_0 \) and \( \dot{x}(0) = v_0 \).

Solution. The characteristic equation of \( m \ddot{x}(t) + k x(t) = 0 \) is \( m r^2 + k = 0 \), whose roots are \( r = \pm \mathrm{i} \sqrt{k/m} \). Thus the general solution is

\[ x(t) = A \cos\big(\omega t\big) + B \sin\big(\omega t\big), \quad \omega = \sqrt{\frac{k}{m}}. \]

Imposing \( x(0) = x_0 \) gives \( A = x_0 \). Differentiating, \( \dot{x}(t) = -A \omega \sin(\omega t) + B \omega \cos(\omega t) \), and the initial condition \( \dot{x}(0) = v_0 \) yields \( v_0 = B \omega \), so \( B = v_0 / \omega \). Therefore

\[ x(t) = x_0 \cos(\omega t) + \frac{v_0}{\omega} \sin(\omega t), \quad \omega = \sqrt{\frac{k}{m}}. \]

This solution describes undamped oscillations with angular frequency \( \omega \), which later will be interpreted as the natural frequency of the floor.

Problem 3 (First-Order Form for Actuator–Load System). For the actuator–load system in Section 5, derive a first-order system in the state vector \( \mathbf{z}(t) = [x_a(t), \dot{x}_a(t), x_L(t), \dot{x}_L(t)]^{\top} \).

Solution. Define \( z_1 = x_a, z_2 = \dot{x}_a, z_3 = x_L, z_4 = \dot{x}_L \). Then

\[ \dot{z}_1(t) = z_2(t),\quad \dot{z}_3(t) = z_4(t). \]

The second derivatives follow from the equations of motion. Recall the coupling force

\[ F_c(t) = k_c\big(z_1(t) - z_3(t)\big) + c_c\big(z_2(t) - z_4(t)\big). \]

The accelerations are

\[ \ddot{x}_a(t) = \dot{z}_2(t) = \frac{u(t) - F_c(t) - b_a z_2(t)}{m_a}, \]

\[ \ddot{x}_L(t) = \dot{z}_4(t) = \frac{F_c(t) - b_L z_4(t)}{m_L}. \]

Substituting \( F_c(t) \) yields a closed first-order system \( \dot{\mathbf{z}}(t) = \mathbf{g}(\mathbf{z}(t),u(t)) \) with four state variables, suitable for numerical integration as in the Java code.

Problem 4 (Limiting Case of Quarter-Car Tire Stiffness). For the quarter-car model, consider the limit \( k_t \to \infty \) while keeping \( x_r(t) \) bounded. Show that the unsprung mass displacement \( x_u(t) \) approaches the road displacement \( x_r(t) \) and derive the limiting equation for the sprung mass.

Solution. The unsprung mass equation can be written as

\[ m_u \ddot{x}_u(t) + c_s\big(\dot{x}_u(t) - \dot{x}_s(t)\big) + k_s\big(x_u(t) - x_s(t)\big) + k_t\big(x_u(t) - x_r(t)\big) = 0. \]

If \( k_t \) becomes very large while the other parameters remain fixed, the term \( k_t(x_u(t) - x_r(t)) \) dominates unless \( x_u(t) - x_r(t) \) becomes small. To keep the overall sum finite, the only consistent limit is \( x_u(t) \approx x_r(t) \), i.e., the wheel follows the road closely. In this limit the sprung mass equation reduces to the SDOF approximation in Section 3:

\[ m_s \ddot{x}_s(t) + c_s\big(\dot{x}_s(t) - \dot{x}_r(t)\big) + k_s\big(x_s(t) - x_r(t)\big) = 0, \]

or in relative coordinates \( z(t) = x_s(t) - x_r(t) \)

\[ m_s \ddot{z}(t) + c_s \dot{z}(t) + k_s z(t) = - m_s \ddot{x}_r(t). \]

This formalizes the intuitive statement that an infinitely stiff tire enforces equality between wheel and road displacement.

In this lesson we modeled three important classes of translational mechanical systems: vehicle suspensions, building floors, and actuator–load chains. Despite their different physical contexts, all admit a common representation in terms of a mass matrix, damping matrix, stiffness matrix, and external inputs. We derived the governing ODEs using Newton's laws, assembled the corresponding matrices, and highlighted limiting cases (such as very stiff tires) that lead to reduced-order SDOF approximations.

We then translated these models into numerical simulations in Python, C++, Java, MATLAB/Simulink, and Mathematica using first-order ODE formulations. These concrete implementations form the basis for the input–output, transfer-function, and state-space analyses developed in subsequent chapters.