In earlier chapters, many systems were represented accurately by linear ordinary differential equations because (i) motions were small, (ii) components operated within designed ranges, and (iii) material and geometric relations were approximately linear. In practice, nonlinearities are not “rare exceptions”; they are the default outside narrow operating regions.

A useful engineering view is: nonlinearities arise when one (or more) of the following is true:

• Geometry is nonlinear (e.g., exact pendulum kinematics uses \( \sin(\theta) \)).
• Constitutive laws are nonlinear (stress–strain, magnetic hysteresis, aerodynamic drag).
• Constraints/switching exist (contact, impacts, diodes, relays, backlash).
• Components have limits (actuator saturation, rate limits, dead-zones, quantization).
• Dissipation is nonlinear (Coulomb friction, stick–slip, hysteresis with memory).

In this chapter, we represent nonlinear systems primarily in the state-space form \( \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},\mathbf{u},t) \), \( \mathbf{y} = \mathbf{g}(\mathbf{x},\mathbf{u},t) \). A system is linear only if these maps are linear in the precise sense given next.

Let \( \mathcal{T} \) denote the input–output operator mapping an input signal \( u(t) \) to an output signal \( y(t)=\mathcal{T}[u](t) \). The operator is linear if and only if it satisfies superposition:

\[ \mathcal{T}[a u_1 + b u_2] = a\,\mathcal{T}[u_1] + b\,\mathcal{T}[u_2] \quad \text{for all scalars } a,b \text{ and all admissible inputs } u_1,u_2. \]

This single axiom implies both additivity and homogeneity. Two immediate consequences are often used as quick checks:

\[ \mathcal{T}[0] = 0, \qquad \mathcal{T}[\alpha u] = \alpha\,\mathcal{T}[u]. \]

Proof (that linearity implies \( \mathcal{T}[0]=0 \)): Take \( a=b=0 \) and any inputs \( u_1,u_2 \). Then \( \mathcal{T}[0] = 0\cdot\mathcal{T}[u_1] + 0\cdot\mathcal{T}[u_2]=0 \).

For memoryless nonlinearities \( y = \varphi(u) \), linearity reduces to: \( \varphi(a u_1 + b u_2)=a\varphi(u_1)+b\varphi(u_2) \). For example, \( \varphi(u)=u^2 \) is nonlinear because:

\[ (u_1+u_2)^2 = u_1^2 + 2u_1u_2 + u_2^2 \neq u_1^2 + u_2^2. \]

For state-space models, linearity means:

\[ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u},\qquad \mathbf{y} = \mathbf{C}\mathbf{x} + \mathbf{D}\mathbf{u}, \]

where \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) are constant matrices (time-invariant case). Any dependence like \( x_1 x_2 \), \( \sin(x_1) \), \( |x_2| \), or piecewise switching breaks linearity.

Two complementary classifications are always useful: (i) where the nonlinearity comes from physically, and (ii) how it appears mathematically. The diagram below summarizes this mapping.

flowchart LR
  A["Observed nonlinearity \n(data or physics)"] --> B["Locate source"]
  B --> C1["Structure/geometry"]
  B --> C2["Materials/constitutive law"]
  B --> C3["Components \n(actuator/sensor)"]
  B --> C4["Interconnections/constraints"]
  B --> C5["Dissipation \n(friction/hysteresis)"]

  C1 --> T1["Examples: \nsin(theta), products like x*v"]
  C2 --> T2["Examples: \nsigma = E*eps + alpha*eps^3"]
  C3 --> T3["Examples: saturation, \ndead-zone, relay, quantization"]
  C4 --> T4["Examples: impacts, \nunilateral contact, backlash"]
  C5 --> T5["Examples: Coulomb friction, \nPreisach-type hysteresis"]

  A --> D["Classify mathematically"]
  D --> M1["Static (memoryless) vs \nDynamic (state-dependent)"]
  D --> M2["Smooth (C1) vs Nonsmooth \n(kinks/discontinuities)"]
  D --> M3["Polynomial vs transcendental \nvs rational"]
  D --> M4["Time-invariant vs time-varying \nparameter dependence"]
  D --> M5["Local vs global nonlinearity \n(operating region)"]
        

A third classification is by regularity. If \( \mathbf{f} \) is continuously differentiable (smooth), local linearization tools (Chapter 9) are well-defined. If the model has discontinuities (relay, ideal Coulomb friction, impacts), careful treatment with event handling or nonsmooth analysis becomes important (we will revisit this in numerical simulation chapters).

Many mechanical systems are linear only after small-angle/small-deflection approximation. The classic example is the pendulum:

\[ \ddot{\theta} + \frac{g}{\ell}\sin(\theta) = 0. \]

The nonlinearity is the \( \sin(\theta) \) term (geometry of the arc). Around \( \theta \approx 0 \), the Taylor approximation \( \sin(\theta) \approx \theta \) yields the linear model:

\[ \ddot{\theta} + \frac{g}{\ell}\theta = 0. \]

Another common geometric source is when forces depend on products of variables, e.g., centripetal terms \( \omega^2 r \) or Coriolis-like couplings. In state-space form, such couplings appear as bilinear terms like \( x_1 x_2 \).

Many nonlinearities are static (memoryless): the output at time \( t \) depends only on the input at the same time. These are ubiquitous in sensors/actuators and signal conditioning.

Actuator saturation. A symmetric saturation nonlinearity is:

\[ \operatorname{sat}(u;U) = \begin{cases} -U & \text{if } u \; < \; -U \\ u & \text{if } |u| \le U \\ U & \text{if } u \; > \; U \end{cases} \]

Proof of nonlinearity (homogeneity fails): choose \( u=U/2 \) and \( \alpha=3 \). Then \( \operatorname{sat}(\alpha u;U)=\operatorname{sat}(3U/2;U)=U \), while \( \alpha\operatorname{sat}(u;U)=3(U/2)=3U/2 \). Since \( U \neq 3U/2 \), homogeneity does not hold.

Dead-zone. A symmetric dead-zone (no response within a band) is:

\[ \operatorname{dz}(x;d) = \begin{cases} 0 & \text{if } |x| \; < \; d \\ x - d\,\operatorname{sgn}(x) & \text{if } |x| \ge d \end{cases} \]

Relay (on–off). A simple relay with output levels \( \pm M \):

\[ \operatorname{relay}(u;M) = M\,\operatorname{sgn}(u). \]

Quantization. For step size \( \Delta \): \( q(u)=\Delta\,\operatorname{round}(u/\Delta) \). Quantization is static but discontinuous, often requiring careful simulation.

Dissipation mechanisms often produce nonlinear force–velocity relations. A standard idealization is Coulomb + viscous friction:

\[ F_{\text{fric}}(\dot{x}) = c\dot{x} + F_c\,\operatorname{sgn}(\dot{x}). \]

The sign function is discontinuous at \( \dot{x}=0 \), which can cause numerical difficulties. A common regularization used in simulation replaces \( \operatorname{sgn}(v) \) by \( \tanh(v/v_s) \):

\[ F_{\text{fric}}(v) \approx c v + F_c\tanh\!\left(\frac{v}{v_s}\right), \quad v_s > 0. \]

Hysteresis introduces memory: the output depends on the history of the input, not only its instantaneous value. A mathematically clean viewpoint is to model hysteresis as an operator \( y(t)=\mathcal{H}[u](t) \), where \( \mathcal{H} \) is not memoryless. In later lessons, hysteresis can produce multiple effective equilibria and loop-like input–output curves.

In System Dynamics, we almost always want a state model \( \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u},t) \). The figure below is a practical, repeatable workflow.

flowchart TD
  S["Start from physical principles"] --> I["Choose generalized coordinates & states"]
  I --> E["Write energy/force balances"]
  E --> N["Identify nonlinear terms (geometry, material, components)"]
  N --> R["Select representation: xdot=f(x,u,t), y=g(x,u,t)"]
  R --> SM["Check smoothness: differentiable? kinks/switching?"]
  SM --> SIM["Pick simulation approach: adaptive ODE / event handling"]
  SIM --> ID["Estimate/validate parameters with data"]
  ID --> V["Verify model: units, limits, special cases"]
  V --> O["Use model for analysis (equilibria, trajectories)"]
        

A useful special form is input-affine dynamics: \( \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})+\mathbf{G}(\mathbf{x})\mathbf{u} \). Many mechanical systems with forces that add linearly with inputs fit this form, even when their internal physics is nonlinear.

Consider a single-DOF mass–spring–damper driven by an actuator. Three nonlinearities are included: (i) cubic stiffness (\( k_3 x^3 \)), (ii) Coulomb friction (regularized), and (iii) actuator saturation on the input.

\[ m\ddot{x} + c\dot{x} + kx + k_3 x^3 + F_c\tanh\!\left(\frac{\dot{x}}{v_s}\right) = b\,\operatorname{sat}(u(t);U). \]

Define states \( x_1=x \), \( x_2=\dot{x} \). Then:

\[ \dot{x}_1 = x_2, \qquad \dot{x}_2 = \frac{ b\operatorname{sat}(u(t);U) - c x_2 - k x_1 - k_3 x_1^3 - F_c\tanh\left(\frac{x_2}{v_s}\right) }{m}. \]

This model is nonlinear due to \( x_1^3 \), \( \tanh(\cdot) \), and saturation. In Lesson 2 we will visualize its trajectories in the phase plane; here we focus on identifying and encoding the nonlinear terms correctly.

The following implementations simulate the model from Section 8 and export trajectories to CSV. Recommended libraries:

• Python: numpy, scipy.integrate, matplotlib.
• C++: standard library only (RK4 from scratch; optional plotting via external tools).
• Java: standard library only (RK4 from scratch; CSV output).
• MATLAB: ode45 and a script that builds a Simulink model.
• Wolfram Mathematica: NDSolve and CSV export.

File: Chapter14_Lesson1.py

"""
System Dynamics — Chapter 14 (Nonlinear System Dynamics)
Lesson 1: Sources and Types of Nonlinearities in Engineering Systems

This script demonstrates a single-DOF nonlinear mass–spring–damper with:
- Cubic stiffness (geometric/material nonlinearity proxy)
- Coulomb + viscous friction (dissipative nonlinearity, smoothed)
- Actuator saturation and optional sensor dead-zone (static, memoryless nonlinearities)

Outputs: plots and a CSV trajectory for cross-language comparison.
"""

from __future__ import annotations

import math
import csv
from dataclasses import dataclass

import numpy as np

# Optional SciPy; the script falls back to an explicit RK4 if SciPy is unavailable.
try:
    from scipy.integrate import solve_ivp
    _HAVE_SCIPY = True
except Exception:
    _HAVE_SCIPY = False

import matplotlib.pyplot as plt


def sat(u: float, umax: float) -> float:
    """Symmetric saturation."""
    return float(np.clip(u, -umax, umax))


def deadzone(x: float, d: float) -> float:
    """
    Symmetric dead-zone:
        y = 0                 if |x| < d
        y = x - d*sign(x)     otherwise
    """
    ax = abs(x)
    if ax < d:
        return 0.0
    return x - d * math.copysign(1.0, x)


def coulomb_friction(v: float, Fc: float, vs: float) -> float:
    """
    Smoothed Coulomb friction ~ Fc*sign(v), using tanh(v/vs) to avoid discontinuity at v=0.
    Smaller vs -> sharper transition.
    """
    return Fc * math.tanh(v / vs)


@dataclass
class Params:
    m: float = 1.0       # kg
    c: float = 0.4       # N*s/m
    k: float = 4.0       # N/m
    k3: float = 8.0      # N/m^3 (cubic stiffness)
    Fc: float = 0.8      # N (Coulomb level)
    vs: float = 0.02     # m/s (smoothing velocity)
    b: float = 1.0       # N per unit input
    umax: float = 1.5    # saturation limit (input units)
    dz: float = 0.0      # sensor dead-zone (meters); set >0 to enable


def input_u(t: float) -> float:
    """
    Excitation signal (can be replaced by any input profile):
    a combination of low-frequency sine + small bias.
    """
    return 1.2 * math.sin(1.0 * t) + 0.3 * math.sin(3.0 * t)


def dynamics(t: float, x: np.ndarray, p: Params) -> np.ndarray:
    """
    State x = [position, velocity].
    Model:
        m xdd + c xd + k x + k3 x^3 + Fc*tanh(xd/vs) = b*sat(u(t), umax)
    """
    pos, vel = float(x[0]), float(x[1])

    # Optional sensor dead-zone (used only to illustrate a measurement nonlinearity)
    pos_meas = deadzone(pos, p.dz) if p.dz > 0 else pos

    u = input_u(t)
    u_sat = sat(u, p.umax)

    spring = p.k * pos + p.k3 * (pos ** 3)
    fric = p.c * vel + coulomb_friction(vel, p.Fc, p.vs)

    acc = (p.b * u_sat - spring - fric) / p.m
    return np.array([vel, acc], dtype=float)


def rk4(f, t0: float, tf: float, x0: np.ndarray, h: float, p: Params):
    """Fixed-step RK4 integrator (educational baseline)."""
    ts = [t0]
    xs = [x0.astype(float).copy()]
    t = t0
    x = x0.astype(float).copy()

    while t < tf - 1e-12:
        if t + h > tf:
            h = tf - t
        k1 = f(t, x, p)
        k2 = f(t + 0.5 * h, x + 0.5 * h * k1, p)
        k3 = f(t + 0.5 * h, x + 0.5 * h * k2, p)
        k4 = f(t + h, x + h * k3, p)
        x = x + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        t = t + h
        ts.append(t)
        xs.append(x.copy())

    return np.array(ts), np.vstack(xs)


def simulate(p: Params, t0: float = 0.0, tf: float = 20.0):
    x0 = np.array([0.4, 0.0], dtype=float)

    if _HAVE_SCIPY:
        sol = solve_ivp(
            fun=lambda t, x: dynamics(t, x, p),
            t_span=(t0, tf),
            y0=x0,
            max_step=0.01,
            rtol=1e-7,
            atol=1e-9,
        )
        t = sol.t
        x = sol.y.T
    else:
        t, x = rk4(dynamics, t0, tf, x0, h=0.005, p=p)

    return t, x


def export_csv(path: str, t: np.ndarray, x: np.ndarray):
    with open(path, "w", newline="", encoding="utf-8") as f:
        w = csv.writer(f)
        w.writerow(["t", "x", "xdot"])
        for ti, (xi, vi) in zip(t, x):
            w.writerow([float(ti), float(xi), float(vi)])


def main():
    p_nl = Params()
    p_lin = Params(k3=0.0, Fc=0.0, umax=1e9)  # "linear-ish" comparison (no cubic, no Coulomb, no sat)

    t1, x1 = simulate(p_nl)
    t2, x2 = simulate(p_lin)

    export_csv("Chapter14_Lesson1_trace_nonlinear.csv", t1, x1)
    export_csv("Chapter14_Lesson1_trace_linear.csv", t2, x2)

    plt.figure()
    plt.plot(t1, x1[:, 0], label="x(t) nonlinear")
    plt.plot(t2, x2[:, 0], label="x(t) linear-ish", linestyle="--")
    plt.xlabel("t [s]")
    plt.ylabel("position x [m]")
    plt.legend()
    plt.grid(True)

    plt.figure()
    plt.plot(t1, x1[:, 1], label="xdot(t) nonlinear")
    plt.plot(t2, x2[:, 1], label="xdot(t) linear-ish", linestyle="--")
    plt.xlabel("t [s]")
    plt.ylabel("velocity xdot [m/s]")
    plt.legend()
    plt.grid(True)

    plt.show()


if __name__ == "__main__":
    main()
      

File: Chapter14_Lesson1_Ex1.py (superposition test for static maps)

"""
Exercise — Chapter 14 Lesson 1 (Nonlinearities)

A quick computational check of linearity properties for memoryless maps f(x).
We test additivity and homogeneity:
    f(x1 + x2) ?= f(x1) + f(x2)
    f(a x) ?= a f(x)
"""

from __future__ import annotations
import numpy as np


def check_linear(f, xs, a=1.7, tol=1e-10):
    x1, x2 = xs
    add = np.max(np.abs(f(x1 + x2) - (f(x1) + f(x2))))
    hom = np.max(np.abs(f(a * x1) - a * f(x1)))
    return add < tol and hom < tol, add, hom


def main():
    rng = np.random.default_rng(7)
    x1 = rng.normal(size=1000)
    x2 = rng.normal(size=1000)

    funcs = {
        "f(x)=x": lambda x: x,
        "f(x)=2x": lambda x: 2.0 * x,
        "f(x)=x+1": lambda x: x + 1.0,
        "f(x)=x^2": lambda x: x**2,
        "f(x)=tanh(x)": lambda x: np.tanh(x),
        "f(x)=sat(x,1)": lambda x: np.clip(x, -1.0, 1.0),
    }

    for name, f in funcs.items():
        ok, add_err, hom_err = check_linear(f, (x1, x2))
        print(f"{name:14s}  linear? {ok}  add_err={add_err:.3e}  hom_err={hom_err:.3e}")


if __name__ == "__main__":
    main()
      

File: Chapter14_Lesson1.cpp

// System Dynamics — Chapter 14 (Nonlinear System Dynamics)
// Lesson 1: Sources and Types of Nonlinearities in Engineering Systems
//
// Nonlinear mass–spring–damper demonstration with RK4 integration.
// Writes trajectory to CSV for plotting in any tool.
//
// Build (GCC/Clang):
//   g++ -O2 -std=c++17 Chapter14_Lesson1.cpp -o Chapter14_Lesson1
//
// Run:
//   ./Chapter14_Lesson1
//
// Output:
//   Chapter14_Lesson1_trace_cpp.csv

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

struct Params {
    double m   = 1.0;
    double c   = 0.4;
    double k   = 4.0;
    double k3  = 8.0;
    double Fc  = 0.8;
    double vs  = 0.02;
    double b   = 1.0;
    double umax= 1.5;
};

static inline double sat(double u, double umax) {
    if (u >  umax) return  umax;
    if (u < -umax) return -umax;
    return u;
}

static inline double coulomb_smooth(double v, double Fc, double vs) {
    return Fc * std::tanh(v / vs);
}

static inline double input_u(double t) {
    return 1.2 * std::sin(1.0 * t) + 0.3 * std::sin(3.0 * t);
}

struct State {
    double x;   // position
    double v;   // velocity
};

static inline State f(double t, const State& s, const Params& p) {
    double u = sat(input_u(t), p.umax);
    double spring = p.k * s.x + p.k3 * s.x * s.x * s.x;
    double fric = p.c * s.v + coulomb_smooth(s.v, p.Fc, p.vs);
    double a = (p.b * u - spring - fric) / p.m;
    return State{ s.v, a };
}

static inline State add(const State& a, const State& b) {
    return State{ a.x + b.x, a.v + b.v };
}

static inline State mul(double h, const State& a) {
    return State{ h * a.x, h * a.v };
}

int main() {
    Params p;
    const double t0 = 0.0, tf = 20.0;
    double h = 0.005;

    State s{0.4, 0.0};

    std::ofstream out("Chapter14_Lesson1_trace_cpp.csv");
    out << "t,x,xdot\n";

    double t = t0;
    while (t < tf - 1e-12) {
        out << t << "," << s.x << "," << s.v << "\n";

        if (t + h > tf) h = tf - t;

        State k1 = f(t, s, p);
        State k2 = f(t + 0.5*h, add(s, mul(0.5*h, k1)), p);
        State k3 = f(t + 0.5*h, add(s, mul(0.5*h, k2)), p);
        State k4 = f(t + h,     add(s, mul(h,     k3)), p);

        s = add(s, mul(h/6.0, add(add(k1, mul(2.0, k2)), add(mul(2.0, k3), k4))));
        t += h;
    }
    out << tf << "," << s.x << "," << s.v << "\n";

    std::cout << "Wrote: Chapter14_Lesson1_trace_cpp.csv\n";
    return 0;
}
      

File: Chapter14_Lesson1.java

// System Dynamics — Chapter 14 (Nonlinear System Dynamics)
// Lesson 1: Sources and Types of Nonlinearities in Engineering Systems
//
// Nonlinear mass–spring–damper with RK4 integration.
// Writes a CSV trajectory.
//
// Compile:
//   javac Chapter14_Lesson1.java
//
// Run:
//   java Chapter14_Lesson1
//
// Output:
//   Chapter14_Lesson1_trace_java.csv

import java.io.FileWriter;
import java.io.IOException;

public class Chapter14_Lesson1 {

    static class Params {
        double m = 1.0;
        double c = 0.4;
        double k = 4.0;
        double k3 = 8.0;
        double Fc = 0.8;
        double vs = 0.02;
        double b = 1.0;
        double umax = 1.5;
    }

    static class State {
        double x;
        double v;
        State(double x, double v) { this.x = x; this.v = v; }
    }

    static double sat(double u, double umax) {
        if (u >  umax) return  umax;
        if (u < -umax) return -umax;
        return u;
    }

    static double coulombSmooth(double v, double Fc, double vs) {
        return Fc * Math.tanh(v / vs);
    }

    static double inputU(double t) {
        return 1.2 * Math.sin(1.0 * t) + 0.3 * Math.sin(3.0 * t);
    }

    static State f(double t, State s, Params p) {
        double u = sat(inputU(t), p.umax);
        double spring = p.k * s.x + p.k3 * s.x * s.x * s.x;
        double fric = p.c * s.v + coulombSmooth(s.v, p.Fc, p.vs);
        double a = (p.b * u - spring - fric) / p.m;
        return new State(s.v, a);
    }

    static State add(State a, State b) {
        return new State(a.x + b.x, a.v + b.v);
    }

    static State mul(double h, State a) {
        return new State(h * a.x, h * a.v);
    }

    public static void main(String[] args) throws IOException {
        Params p = new Params();

        double t0 = 0.0, tf = 20.0;
        double h = 0.005;

        State s = new State(0.4, 0.0);

        try (FileWriter out = new FileWriter("Chapter14_Lesson1_trace_java.csv")) {
            out.write("t,x,xdot\n");

            double t = t0;
            while (t < tf - 1e-12) {
                out.write(t + "," + s.x + "," + s.v + "\n");

                double hh = h;
                if (t + hh > tf) hh = tf - t;

                State k1 = f(t, s, p);
                State k2 = f(t + 0.5*hh, add(s, mul(0.5*hh, k1)), p);
                State k3 = f(t + 0.5*hh, add(s, mul(0.5*hh, k2)), p);
                State k4 = f(t + hh,     add(s, mul(hh,     k3)), p);

                State incr = add(add(k1, mul(2.0, k2)), add(mul(2.0, k3), k4));
                s = add(s, mul(hh/6.0, incr));
                t += hh;
            }
            out.write(tf + "," + s.x + "," + s.v + "\n");
        }

        System.out.println("Wrote: Chapter14_Lesson1_trace_java.csv");
    }
}
      

File: Chapter14_Lesson1.m

% System Dynamics — Chapter 14 (Nonlinear System Dynamics)
% Lesson 1: Sources and Types of Nonlinearities in Engineering Systems
%
% Nonlinear mass–spring–damper simulation using ode45.
% Saves trajectory to CSV for cross-language comparison.

clear; clc;

p.m = 1.0;
p.c = 0.4;
p.k = 4.0;
p.k3 = 8.0;
p.Fc = 0.8;
p.vs = 0.02;
p.b = 1.0;
p.umax = 1.5;

tspan = [0 20];
x0 = [0.4; 0.0];  % [position; velocity]

opts = odeset('RelTol',1e-7,'AbsTol',1e-9,'MaxStep',0.01);
[t, X] = ode45(@(t,x) dyn(t,x,p), tspan, x0, opts);

% Export to CSV
T = table(t, X(:,1), X(:,2), 'VariableNames', {'t','x','xdot'});
writetable(T, 'Chapter14_Lesson1_trace_matlab.csv');

% Plot
figure; plot(t, X(:,1)); grid on;
xlabel('t [s]'); ylabel('x [m]'); title('Nonlinear position');

figure; plot(t, X(:,2)); grid on;
xlabel('t [s]'); ylabel('xdot [m/s]'); title('Nonlinear velocity');

disp('Wrote: Chapter14_Lesson1_trace_matlab.csv');

function dx = dyn(t, x, p)
    pos = x(1);
    vel = x(2);

    u = input_u(t);
    u_sat = sat(u, p.umax);

    spring = p.k*pos + p.k3*pos^3;
    fric = p.c*vel + p.Fc*tanh(vel/p.vs);

    acc = (p.b*u_sat - spring - fric) / p.m;

    dx = [vel; acc];
end

function u = input_u(t)
    u = 1.2*sin(1.0*t) + 0.3*sin(3.0*t);
end

function y = sat(u, umax)
    y = min(max(u, -umax), umax);
end
      

File: Chapter14_Lesson1_Simulink.m (builds Chapter14_Lesson1_Simulink.slx)

% System Dynamics — Chapter 14 (Nonlinear System Dynamics)
% Lesson 1: Sources and Types of Nonlinearities in Engineering Systems
%
% Programmatically builds a Simulink model that implements:
%   m xdd + c xd + k x + k3 x^3 + Fc*tanh(xd/vs) = b*sat(u)
%
% The model uses standard Simulink blocks: Integrator, Sum, Gain, Math Function,
% Saturation, Trigonometric Function, and a MATLAB Function block for tanh friction.
%
% Run this script in MATLAB. It creates and opens 'Chapter14_Lesson1_Simulink.slx'.

clear; clc;

model = 'Chapter14_Lesson1_Simulink';
if bdIsLoaded(model), close_system(model, 0); end
new_system(model); open_system(model);

% Parameters (in base workspace so blocks can use them)
m = 1.0; c = 0.4; k = 4.0; k3 = 8.0; Fc = 0.8; vs = 0.02; b = 1.0; umax = 1.5; %#ok<NASGU>
assignin('base','m',m); assignin('base','c',c); assignin('base','k',k);
assignin('base','k3',k3); assignin('base','Fc',Fc); assignin('base','vs',vs);
assignin('base','b',b); assignin('base','umax',umax);

% Layout helpers
x0 = 40; y0 = 50; dx = 120; dy = 70;

% Input u(t) = 1.2*sin(t) + 0.3*sin(3t)
add_block('simulink/Sources/Clock', [model '/t'], 'Position', [x0 y0 x0+30 y0+30]);
add_block('simulink/Math Operations/Gain', [model '/w1'], 'Gain', '1', 'Position', [x0+dx y0 x0+dx+40 y0+30]);
add_block('simulink/Math Operations/Gain', [model '/w3'], 'Gain', '3', 'Position', [x0+dx y0+dy x0+dx+40 y0+dy+30]);
add_block('simulink/Math Operations/Trigonometric Function', [model '/sin1'], 'Operator', 'sin', 'Position', [x0+2*dx y0 x0+2*dx+50 y0+30]);
add_block('simulink/Math Operations/Trigonometric Function', [model '/sin3'], 'Operator', 'sin', 'Position', [x0+2*dx y0+dy x0+2*dx+50 y0+dy+30]);
add_block('simulink/Math Operations/Gain', [model '/a1'], 'Gain', '1.2', 'Position', [x0+3*dx y0 x0+3*dx+40 y0+30]);
add_block('simulink/Math Operations/Gain', [model '/a3'], 'Gain', '0.3', 'Position', [x0+3*dx y0+dy x0+3*dx+40 y0+dy+30]);
add_block('simulink/Math Operations/Sum', [model '/sum_u'], 'Inputs', '++', 'Position', [x0+4*dx y0+dy/2 x0+4*dx+40 y0+dy/2+40]);

% Saturation block
add_block('simulink/Discontinuities/Saturation', [model '/sat_u'], ...
    'UpperLimit', 'umax', 'LowerLimit', '-umax', ...
    'Position', [x0+5*dx y0+dy/2 x0+5*dx+60 y0+dy/2+40]);

% Dynamics: two integrators xdd -> xd -> x
add_block('simulink/Continuous/Integrator', [model '/int_v'], 'Position', [x0+7*dx y0 x0+7*dx+40 y0+40]);
add_block('simulink/Continuous/Integrator', [model '/int_x'], 'Position', [x0+8*dx y0 x0+8*dx+40 y0+40]);

% Compute spring: k*x + k3*x^3
add_block('simulink/Math Operations/Gain', [model '/k*x'], 'Gain', 'k', 'Position', [x0+8*dx y0+dy x0+8*dx+60 y0+dy+30]);
add_block('simulink/Math Operations/Math Function', [model '/x^3'], 'Operator', 'pow', 'Exponent', '3', ...
    'Position', [x0+8*dx y0+2*dy x0+8*dx+60 y0+2*dy+30]);
add_block('simulink/Math Operations/Gain', [model '/k3*x^3'], 'Gain', 'k3', 'Position', [x0+9*dx y0+2*dy x0+9*dx+70 y0+2*dy+30]);

add_block('simulink/Math Operations/Sum', [model '/sum_spring'], 'Inputs', '++', 'Position', [x0+10*dx y0+dy+25 x0+10*dx+40 y0+dy+65]);

% Friction: c*v + Fc*tanh(v/vs) via MATLAB Function block
add_block('simulink/Math Operations/Gain', [model '/c*v'], 'Gain', 'c', 'Position', [x0+8*dx y0+3*dy x0+8*dx+60 y0+3*dy+30]);
add_block('simulink/User-Defined Functions/MATLAB Function', [model '/tanh_fric'], ...
    'Position', [x0+9*dx y0+3*dy x0+9*dx+130 y0+3*dy+60]);
set_param([model '/tanh_fric'], 'Script', sprintf([ ...
    'function y = f(v)\n' ...
    '%% Smoothed Coulomb friction: Fc*tanh(v/vs)\n' ...
    'Fc = evalin(''base'',''Fc'');\n' ...
    'vs = evalin(''base'',''vs'');\n' ...
    'y = Fc*tanh(v/vs);\n' ...
    'end\n' ...
]));

add_block('simulink/Math Operations/Sum', [model '/sum_fric'], 'Inputs', '++', 'Position', [x0+10*dx y0+3*dy x0+10*dx+40 y0+3*dy+40]);

% Force balance: (b*sat(u) - spring - fric)/m
add_block('simulink/Math Operations/Gain', [model '/b*u'], 'Gain', 'b', 'Position', [x0+6*dx y0+dy/2 x0+6*dx+50 y0+dy/2+40]);

add_block('simulink/Math Operations/Sum', [model '/sum_force'], 'Inputs', '+--', 'Position', [x0+11*dx y0+2*dy x0+11*dx+50 y0+2*dy+60]);
add_block('simulink/Math Operations/Gain', [model '/1/m'], 'Gain', '1/m', 'Position', [x0+12*dx y0+2*dy x0+12*dx+50 y0+2*dy+30]);

% Scopes
add_block('simulink/Sinks/Scope', [model '/Scope_x'], 'Position', [x0+9*dx y0-10 x0+9*dx+70 y0+30]);
add_block('simulink/Sinks/Scope', [model '/Scope_v'], 'Position', [x0+8*dx y0-10 x0+8*dx+70 y0+30]);

% Wiring
add_line(model,'t/1','w1/1');        add_line(model,'t/1','w3/1');
add_line(model,'w1/1','sin1/1');     add_line(model,'w3/1','sin3/1');
add_line(model,'sin1/1','a1/1');     add_line(model,'sin3/1','a3/1');
add_line(model,'a1/1','sum_u/1');    add_line(model,'a3/1','sum_u/2');
add_line(model,'sum_u/1','sat_u/1'); add_line(model,'sat_u/1','b*u/1');

add_line(model,'b*u/1','sum_force/1');
add_line(model,'1/m/1','int_v/1');
add_line(model,'int_v/1','int_x/1');

% Feedback signals
add_line(model,'int_x/1','k*x/1');
add_line(model,'int_x/1','x^3/1');
add_line(model,'x^3/1','k3*x^3/1');

add_line(model,'k*x/1','sum_spring/1');
add_line(model,'k3*x^3/1','sum_spring/2');

add_line(model,'int_v/1','c*v/1');
add_line(model,'int_v/1','tanh_fric/1');
add_line(model,'c*v/1','sum_fric/1');
add_line(model,'tanh_fric/1','sum_fric/2');

add_line(model,'sum_spring/1','sum_force/2');
add_line(model,'sum_fric/1','sum_force/3');
add_line(model,'sum_force/1','1/m/1');

% Scopes
add_line(model,'int_x/1','Scope_x/1');
add_line(model,'int_v/1','Scope_v/1');

set_param(model,'StopTime','20');
save_system(model);
open_system(model);
disp('Created and opened Chapter14_Lesson1_Simulink.slx');
      

File: Chapter14_Lesson1.nb

(* ::Package:: *)

(* System Dynamics — Chapter 14 (Nonlinear System Dynamics)
   Lesson 1: Sources and Types of Nonlinearities in Engineering Systems

   Mathematica notebook content (plain Wolfram Language).
   Simulates the same nonlinear mass–spring–damper model using NDSolve.
*)

ClearAll["Global`*"];

params = <|
  "m" -> 1.0,
  "c" -> 0.4,
  "k" -> 4.0,
  "k3" -> 8.0,
  "Fc" -> 0.8,
  "vs" -> 0.02,
  "b" -> 1.0,
  "umax" -> 1.5
|>;

sat[u_, umax_] := Clip[u, {-umax, umax}];
uIn[t_] := 1.2 Sin[t] + 0.3 Sin[3 t];
fric[v_, Fc_, vs_] := Fc Tanh[v/vs];

m = params["m"]; c = params["c"]; k = params["k"]; k3 = params["k3"];
Fc = params["Fc"]; vs = params["vs"]; b = params["b"]; umax = params["umax"];

eqn = m x''[t] + c x'[t] + k x[t] + k3 x[t]^3 + fric[x'[t], Fc, vs] == b sat[uIn[t], umax];

ic = {x[0] == 0.4, x'[0] == 0.0};

sol = NDSolve[{eqn, ic}, x, {t, 0, 20}, MaxStepFraction -> 1/2000][[1]];

p1 = Plot[Evaluate[x[t] /. sol], {t, 0, 20}, GridLines -> Automatic, PlotLabel -> "x(t)"];
p2 = Plot[Evaluate[x'[t] /. sol], {t, 0, 20}, GridLines -> Automatic, PlotLabel -> "x'(t)"];

Print[p1];
Print[p2];

data = Table[{tt, (x[tt] /. sol), (x'[tt] /. sol)}, {tt, 0, 20, 0.005}];
Export["Chapter14_Lesson1_trace_mathematica.csv", Prepend[data, {"t", "x", "xdot"}]];
Print["Wrote: Chapter14_Lesson1_trace_mathematica.csv"];
      

Cross-checking tip: all languages export CSV traces. If you use the same step size and parameters, the curves should closely match (small differences are expected due to solver details).

Problem 1 (Linearity of an affine map): Let \( y = a x + b \) with constants \( a,b \). Determine for which values of \( b \) the mapping is linear.

Solution: A linear map must satisfy \( \varphi(0)=0 \). Here \( \varphi(0)=b \), so linearity requires \( b=0 \). Therefore \( y=ax \) is linear; any nonzero offset makes it nonlinear (affine).

Problem 2 (Pendulum nonlinearity and small-angle model): Starting from \( \ddot{\theta} + \frac{g}{\ell}\sin(\theta)=0 \), show that the system is nonlinear and derive the small-angle linear approximation about \( \theta=0 \).

Solution: The restoring term is \( \sin(\theta) \), which violates superposition (e.g., \( \sin(\theta_1+\theta_2)\neq\sin\theta_1+\sin\theta_2 \)). For small angles, use the Taylor series:

\[ \sin(\theta) = \theta - \frac{\theta^3}{6} + \cdots \approx \theta \quad (|\theta| \text{ small}). \]

Substituting gives the linear model \( \ddot{\theta} + \frac{g}{\ell}\theta=0 \).

Problem 3 (Saturation breaks homogeneity): For \( y=\operatorname{sat}(u;U) \), prove it is nonlinear by showing homogeneity fails.

Solution: Pick any \( u \) with \( 0 < u < U \) and choose \( \alpha \) such that \( \alpha u > U \). Then \( \operatorname{sat}(\alpha u;U)=U \), but \( \alpha\operatorname{sat}(u;U)=\alpha u > U \). Hence \( \operatorname{sat}(\alpha u;U) \neq \alpha\operatorname{sat}(u;U) \).

Problem 4 (Equilibria with a nonlinear spring): Consider \( m\ddot{x} + kx + k_3 x^3 = F_0 \) with constant force \( F_0 \). Find the equilibrium points.

Solution: At equilibrium \( \dot{x}=0 \) and \( \ddot{x}=0 \), so \( kx_e + k_3 x_e^3 = F_0 \). This is a cubic equation:

\[ k_3 x_e^3 + k x_e - F_0 = 0. \]

Depending on parameters, it can have one or three real roots. Even before doing stability analysis, this illustrates a key nonlinear feature: multiple equilibria can exist for the same constant input.

Problem 5 (Classify nonlinearities in a block description): A servo axis has (i) amplifier saturation, (ii) gearbox backlash, and (iii) Coulomb friction. Classify each as static/dynamic and smooth/nonsmooth (conceptually).

Solution:

• Saturation: static (memoryless), nonsmooth (kinks at limits).
• Backlash: dynamic with memory (hysteresis-like), nonsmooth (switching when engaging).
• Coulomb friction: static with respect to velocity in the ideal model, nonsmooth (discontinuity at velocity zero).

Problem 6 (Detect a nonlinearity from data): Suppose you apply inputs \( u_1(t) \) and \( u_2(t) \) separately to a system and measure outputs \( y_1(t), y_2(t) \). You then apply \( u_1(t)+u_2(t) \) and measure \( y_{12}(t) \). Propose a quantitative test statistic for linearity.

Solution: A practical test is the normalized superposition error:

\[ \varepsilon = \frac{\|y_{12} - (y_1+y_2)\|_2}{\|y_{12}\|_2 + \|y_1+y_2\|_2}. \]

For an ideal linear system (with the same initial conditions), \( \varepsilon=0 \). Nonzero values indicate nonlinearity and/or changing operating conditions.

We defined linearity precisely using the superposition axiom, then organized nonlinearities by physical origin (geometry, materials, components, constraints, dissipation) and mathematical form (static vs dynamic, smooth vs nonsmooth, polynomial/transcendental/piecewise). We also encoded a canonical nonlinear mass–spring–damper model and implemented it in multiple languages. In Lesson 2 we will use these models to study equilibria and trajectories in the phase plane.