Consider a point mass that can move along a straight line (translational degree of freedom). The mass is connected to a rigid support by a linear spring and a viscous damper. Let \( x(t) \) denote the displacement of the mass from a reference position, and let \( f(t) \) denote an externally applied force (e.g., from an actuator).

We adopt the following modeling assumptions, standard in linear system dynamics:

• The mass is rigid and its motion is constrained to one dimension.
• The spring is linear (Hookean): force is proportional to displacement.
• The damper is linear and viscous: force is proportional to velocity.
• All parameters are constant in time (time-invariant model).
• Deformations and velocities are small enough that nonlinear effects can be neglected.

The constitutive relations for the spring and damper are

\[ F_s(t) = -k\,x(t), \quad k > 0, \]

\[ F_d(t) = -c\,\dot{x}(t), \quad c \ge 0, \]

where \( F_s(t) \) is the spring force, \( F_d(t) \) is the damping force, \( k \) is the stiffness (units N/m), and \( c \) is the viscous damping coefficient (units N·s/m).

flowchart TD
  W["Rigid support"] --> S["Spring (stiffness k)"]
  S --> D["Damper (damping c)"]
  D --> M["Mass m with displacement x(t)"]
  M --> F["Applied force f(t)"]
        

The free-body diagram of the mass, which you constructed in Lesson 1, shows the forces acting on the mass along the motion direction: the spring force, the damping force, and the external force.

Let the positive direction of \( x(t) \) be to the right. The forces acting on the mass in the direction of motion are:

• Spring force: \( F_s(t) = -k\,x(t) \) (opposes displacement).
• Damping force: \( F_d(t) = -c\,\dot{x}(t) \) (opposes velocity).
• External force: \( f(t) \) (actuator or disturbance).

Newton’s second law for translational motion states that \( \sum F = m\,\ddot{x}(t) \), where \( m \) is the mass. Taking forces positive in the direction of increasing \( x \),

\[ f(t) + F_s(t) + F_d(t) = m\,\ddot{x}(t). \]

Substituting the constitutive relations,

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

Rearranging terms gives the standard form of the mass–spring–damper ODE:

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

This is a linear, constant-coefficient, second-order ODE that relates the displacement \( x(t) \) to the input force \( f(t) \). For many control-engineering applications, \( x(t) \) is the output of interest, and \( f(t) \) is the control input.

If the external force is zero, \( f(t) = 0 \), we obtain the free vibration model:

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

If \( f(t) = F_0 \) is constant, the system has a static equilibrium where the spring force balances the applied load:

\[ k\,x_{\text{eq}} = F_0 \quad \Rightarrow \quad x_{\text{eq}} = \frac{F_0}{k}. \]

It is convenient to express the ODE in a normalized form that separates inertial, stiffness, and damping effects. Dividing the governing equation by \( m \) yields

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

We define the undamped natural frequency and the damping ratio as

\[ \omega_n = \sqrt{\frac{k}{m}}, \quad \zeta = \frac{c}{2\sqrt{k m}} = \frac{c}{2 m \omega_n}, \]

and obtain the standard second-order form

\[ \ddot{x}(t) + 2\,\zeta\,\omega_n\,\dot{x}(t) + \omega_n^2\,x(t) = \frac{1}{m}f(t). \]

The parameters \( \omega_n \) and \( \zeta \) summarize the dynamic character of the system and will be used repeatedly in later chapters.

3.1 Homogeneous Equation and Characteristic Polynomial

For free vibration (\( f(t) = 0 \)), the ODE is

\[ \ddot{x}(t) + 2\,\zeta\,\omega_n\,\dot{x}(t) + \omega_n^2\,x(t) = 0. \]

Seeking exponential solutions of the form \( x(t) = e^{r t} \) leads to the characteristic polynomial

\[ r^2 + 2\,\zeta\,\omega_n\,r + \omega_n^2 = 0, \]

with roots

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

3.2 Damping Regimes

The qualitative behavior depends on the discriminant \( \zeta^2 - 1 \):

• Underdamped: \( 0 \le \zeta < 1 \).
  The roots are complex:

  \[ r_{1,2} = -\zeta\,\omega_n \pm \mathrm{i}\,\omega_d, \quad \omega_d = \omega_n \sqrt{1 - \zeta^2}. \]

  The solution is exponentially decaying oscillation:

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

  where constants \( C_1, C_2 \) are determined by initial conditions.
• Critically damped: \( \zeta = 1 \).
  The characteristic polynomial has a repeated real root \( r = -\omega_n \), and the general solution is

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

• Overdamped: \( \zeta > 1 \).
  The roots are real and distinct, negative:

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

  and the solution is a sum of two decaying exponentials:

  \[ x(t) = C_1 e^{r_1 t} + C_2 e^{r_2 t}. \]

flowchart TD
  A["Start from m x'' + c x' + k x = f(t)"] --> B["Divide by m to get x'' + (c/m) x' + (k/m) x = f(t)/m"]
  B --> C["Define omega_n^2 = k/m, zeta = c/(2 sqrt(k m))"]
  C --> D["Homogeneous: x'' + 2 zeta omega_n x' + omega_n^2 x = 0"]
  D --> E["Characteristic r^2 + 2 zeta omega_n r + omega_n^2 = 0"]
  E --> F["zeta < 1: underdamped \n(complex roots)"]
  E --> G["zeta = 1: critical \n(repeated root)"]
  E --> H["zeta > 1: overdamped \n(real distinct roots)"]
        

A key structural fact is that for \( m > 0 \), \( k > 0 \), and \( c \ge 0 \), the real parts of \( r_{1,2} \) are always non-positive, so the homogeneous response does not grow unbounded in time.

For a general input \( f(t) \), the solution of the ODE can be decomposed as

\[ x(t) = x_h(t) + x_p(t), \]

where \( x_h(t) \) solves the homogeneous equation and \( x_p(t) \) is any particular solution that accounts for the input. For control engineering, an important case is a constant load.

4.1 Constant Force Input

Consider \( f(t) = F_0 \) (constant). Looking for a constant particular solution \( x_p(t) = x_{\text{ss}} \), we have

\[ m \cdot 0 + c \cdot 0 + k\,x_{\text{ss}} = F_0 \quad \Rightarrow \quad x_{\text{ss}} = \frac{F_0}{k}. \]

Thus the total solution is

\[ x(t) = x_h(t) + \frac{F_0}{k}. \]

All transient behavior is contained in \( x_h(t) \), which decays to zero for \( \zeta > 0 \), while the steady displacement is the static deflection \( F_0/k \).

For simulation and later state-space modeling, it is useful to rewrite the second-order ODE as two coupled first-order equations. Introduce auxiliary variables

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

Then

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

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

A simple (though not always accurate) numerical scheme is the explicit Euler method with step size \( h > 0 \) and time grid \( t_k = k h \):

\[ x_1^{k+1} = x_1^{k} + h\,x_2^{k}, \]

\[ x_2^{k+1} = x_2^{k} + h\left( -\frac{c}{m} x_2^{k} - \frac{k}{m} x_1^{k} + \frac{1}{m} f(t_k) \right). \]

This discrete-time approximation will be analyzed more systematically in the numerical methods chapter; here we use it to illustrate simulation code across several languages.

The following Python script simulates free vibration of a mass–spring–damper system using explicit Euler integration. In practice, scientific libraries such as scipy.integrate (e.g., solve_ivp) provide more accurate and robust solvers.

import numpy as np

# Parameters
m = 1.0     # mass [kg]
k = 20.0    # stiffness [N/m]
c = 1.0     # viscous damping [N*s/m]

def f(t):
    # External force (free vibration example)
    return 0.0

# Simulation settings
t_end = 10.0
h = 0.001
N = int(t_end / h)

t = np.linspace(0.0, t_end, N + 1)
x1 = np.zeros(N + 1)  # displacement
x2 = np.zeros(N + 1)  # velocity

# Initial conditions
x1[0] = 0.1  # initial displacement [m]
x2[0] = 0.0  # initial velocity [m/s]

for k_step in range(N):
    # Euler integration
    x1_dot = x2[k_step]
    x2_dot = -(c / m) * x2[k_step] - (k / m) * x1[k_step] + (1.0 / m) * f(t[k_step])

    x1[k_step + 1] = x1[k_step] + h * x1_dot
    x2[k_step + 1] = x2[k_step] + h * x2_dot

# Post-processing (e.g., simple textual output)
print("Final displacement:", x1[-1])
print("Final velocity:", x2[-1])

# For visualization, one would typically use matplotlib:
# import matplotlib.pyplot as plt
# plt.plot(t, x1)
# plt.xlabel("t [s]")
# plt.ylabel("x(t) [m]")
# plt.show()
      

The function f(t) can be modified to represent more general inputs, such as step or sinusoidal forces.

In C++, one can implement the same Euler integrator directly. For larger projects, libraries such as Boost.Odeint provide higher-order and adaptive solvers.

#include <iostream>
#include <vector>

double external_force(double t) {
    // Free vibration example
    return 0.0;
}

int main() {
    // Parameters
    double m = 1.0;
    double k = 20.0;
    double c = 1.0;

    // Simulation settings
    double t_end = 10.0;
    double h = 0.001;
    int N = static_cast<int>(t_end / h);

    std::vector<double> t(N + 1);
    std::vector<double> x1(N + 1);
    std::vector<double> x2(N + 1);

    // Initial conditions
    x1[0] = 0.1;
    x2[0] = 0.0;
    t[0] = 0.0;

    for (int k_step = 0; k_step < N; ++k_step) {
        t[k_step + 1] = t[k_step] + h;

        double f = external_force(t[k_step]);

        double x1_dot = x2[k_step];
        double x2_dot = -(c / m) * x2[k_step]
                        - (k / m) * x1[k_step]
                        + (1.0 / m) * f;

        x1[k_step + 1] = x1[k_step] + h * x1_dot;
        x2[k_step + 1] = x2[k_step] + h * x2_dot;
    }

    std::cout << "Final displacement: " << x1.back() << std::endl;
    std::cout << "Final velocity: " << x2.back() << std::endl;

    return 0;
}
      

In Java, one can implement the same update logic. For more advanced applications, Apache Commons Math provides ODE solvers that are useful for system dynamics simulations.

public class MassSpringDamperEuler {

    static double externalForce(double t) {
        // Free vibration example
        return 0.0;
    }

    public static void main(String[] args) {
        // Parameters
        double m = 1.0;
        double k = 20.0;
        double c = 1.0;

        // Simulation settings
        double tEnd = 10.0;
        double h = 0.001;
        int N = (int) (tEnd / h);

        double[] t = new double[N + 1];
        double[] x1 = new double[N + 1];
        double[] x2 = new double[N + 1];

        // Initial conditions
        x1[0] = 0.1;
        x2[0] = 0.0;
        t[0] = 0.0;

        for (int kStep = 0; kStep < N; ++kStep) {
            t[kStep + 1] = t[kStep] + h;

            double f = externalForce(t[kStep]);

            double x1Dot = x2[kStep];
            double x2Dot = -(c / m) * x2[kStep]
                           - (k / m) * x1[kStep]
                           + (1.0 / m) * f;

            x1[kStep + 1] = x1[kStep] + h * x1Dot;
            x2[kStep + 1] = x2[kStep] + h * x2Dot;
        }

        System.out.println("Final displacement: " + x1[N]);
        System.out.println("Final velocity: " + x2[N]);
    }
}
      

MATLAB's built-in ODE solvers (e.g., ode45) are well suited for simulating mass–spring–damper systems. Below is a simple script using ode45 for free vibration:

function mass_spring_damper_ode45
    % Parameters
    m = 1.0;
    k = 20.0;
    c = 1.0;

    % Time span
    tspan = [0 10];

    % Initial conditions: [x; xdot]
    x0 = [0.1; 0.0];

    % Solve ODE
    [t, x] = ode45(@(t, x) dynamics(t, x, m, c, k), tspan, x0);

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

function dx = dynamics(t, x, m, c, k)
    % External force (free vibration)
    f = 0.0;

    dx = zeros(2, 1);
    dx(1) = x(2);
    dx(2) = -(c / m) * x(2) - (k / m) * x(1) + (1.0 / m) * f;
end
      

In Simulink, the same system can be built graphically by:

• Adding two integrator blocks to represent velocity and displacement.
• Implementing the expression for \( \dot{x}_2(t) \) with gain and sum blocks.
• Feeding an external force block (e.g., step or sine) into the acceleration calculation.

These constructions match the first-order representation derived above.

Wolfram Mathematica provides symbolic and numerical solvers. A direct formulation with NDSolve is:

m = 1.0;
k = 20.0;
c = 1.0;

f[t_] := 0.0; (* free vibration *)

x0 = 0.1;
v0 = 0.0;

sol = NDSolve[
   {
    m * x''[t] + c * x'[t] + k * x[t] == f[t],
    x[0] == x0,
    x'[0] == v0
   },
   x,
   {t, 0, 10}
];

Plot[Evaluate[x[t] /. sol], {t, 0, 10},
  AxesLabel -> {"t [s]", "x(t) [m]"},
  PlotRange -> All
]
      

Mathematica can also obtain closed-form symbolic solutions using DSolve, which is useful for verifying analytical derivations of homogeneous and particular solutions.

Problem 1 (Derivation of the Governing ODE):
A mass \( m \) is attached horizontally to a wall by a linear spring of stiffness \( k \) and a viscous damper with coefficient \( c \) in parallel. An external force \( f(t) \) pulls the mass along a frictionless surface. Using Newton's second law, derive the governing ODE relating the displacement \( x(t) \) to the applied force \( f(t) \).

Solution:
Taking positive displacement to the right, the forces along the direction of motion are:

• Spring: \( F_s(t) = -k\,x(t) \).
• Damper: \( F_d(t) = -c\,\dot{x}(t) \).
• External: \( f(t) \).

Newton's law gives

\[ f(t) + F_s(t) + F_d(t) = m\,\ddot{x}(t). \]

Substituting the constitutive relations,

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

or equivalently

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

Problem 2 (Natural Frequency and Damping Ratio):
A mass–spring–damper system has \( m = 2 \,\text{kg} \), \( k = 18 \,\text{N/m} \), and \( c = 4 \,\text{N·s/m} \). Compute the undamped natural frequency \( \omega_n \), the damping ratio \( \zeta \), and determine whether the system is underdamped, critically damped, or overdamped.

Solution:
The undamped natural frequency is

\[ \omega_n = \sqrt{\frac{k}{m}} = \sqrt{\frac{18}{2}} = \sqrt{9} = 3 \,\text{rad/s}. \]

The damping ratio is

\[ \zeta = \frac{c}{2\sqrt{k m}} = \frac{4}{2\sqrt{18 \cdot 2}} = \frac{4}{2\sqrt{36}} = \frac{4}{2 \cdot 6} = \frac{4}{12} = \frac{1}{3}. \]

Since \( 0 \le \zeta = 1/3 < 1 \), the system is underdamped and exhibits decaying oscillations.

Problem 3 (Free Response with Given Initial Conditions):
Consider the system in Problem 2 with \( f(t) = 0 \), initial displacement \( x(0) = x_0 \) and initial velocity \( \dot{x}(0) = v_0 \). Derive the explicit homogeneous solution \( x_h(t) \) using the underdamped form.

Solution:
As computed, \( \omega_n = 3 \) and \( \zeta = 1/3 \). The damped natural frequency is

\[ \omega_d = \omega_n \sqrt{1 - \zeta^2} = 3 \sqrt{1 - \left(\frac{1}{3}\right)^2} = 3 \sqrt{1 - \frac{1}{9}} = 3 \sqrt{\frac{8}{9}} = 3 \cdot \frac{\sqrt{8}}{3} = \sqrt{8}. \]

For underdamped systems, the homogeneous solution is

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

Imposing initial conditions:

\[ x_h(0) = C_1 = x_0, \]

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

Evaluating at \( t = 0 \):

\[ \dot{x}_h(0) = -\zeta\,\omega_n C_1 + C_2 \omega_d = v_0. \]

With \( C_1 = x_0 \), we solve for \( C_2 \):

\[ C_2 = \frac{v_0 + \zeta\,\omega_n x_0}{\omega_d}. \]

Therefore,

\[ x_h(t) = e^{-\zeta\,\omega_n t} \left( x_0 \cos(\omega_d t) + \frac{v_0 + \zeta\,\omega_n x_0}{\omega_d} \sin(\omega_d t) \right). \]

Problem 4 (Static Deflection under Constant Force):
A mass–spring–damper system with parameters \( m, c, k \) is subjected to a constant force \( f(t) = F_0 \) applied at time \( t = 0 \). Assume the mass starts at rest at the origin: \( x(0) = 0 \), \( \dot{x}(0) = 0 \). Show that the displacement converges to \( F_0/k \) for all damping ratios \( \zeta > 0 \).

Solution:
The governing equation is

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

A constant particular solution \( x_p(t) = x_{\text{ss}} \) satisfies

\[ k\,x_{\text{ss}} = F_0 \quad \Rightarrow \quad x_{\text{ss}} = \frac{F_0}{k}. \]

Let \( y(t) = x(t) - x_{\text{ss}} \). Then

\[ \ddot{y}(t) = \ddot{x}(t), \quad \dot{y}(t) = \dot{x}(t), \]

and substituting into the ODE gives

\[ m\,\ddot{y}(t) + c\,\dot{y}(t) + k\,y(t) = 0. \]

Thus \( y(t) \) satisfies the homogeneous equation. For \( c \ge 0 \) and \( m, k > 0 \), we have \( \zeta \ge 0 \), and the homogeneous solution \( y(t) \) decays to zero as \( t \rightarrow \infty \). Therefore,

\[ \lim_{t \rightarrow \infty} x(t) = \lim_{t \rightarrow \infty} \left( y(t) + x_{\text{ss}} \right) = x_{\text{ss}} = \frac{F_0}{k}. \]

In this lesson, we constructed the canonical mass–spring–damper model for translational motion starting from Newton's laws and linear constitutive relations. The resulting governing equation \( m\,\ddot{x}(t) + c\,\dot{x}(t) + k\,x(t) = f(t) \) is a linear, time-invariant, second-order ODE. By normalizing the equation, we introduced the undamped natural frequency \( \omega_n = \sqrt{k/m} \) and damping ratio \( \zeta = c/(2\sqrt{k m}) \), and showed how these parameters determine the qualitative behavior (underdamped, critically damped, overdamped) of the homogeneous solution. We also examined static deflection under constant forcing and constructed first-order representations suitable for numerical simulation.

Finally, we implemented explicit Euler simulations in Python, C++, Java, MATLAB, and Mathematica, laying the computational foundation for more advanced numerical and control-oriented analyses in later chapters.