Consider an autonomous nonlinear system \( \dot{\mathbf{x} } = \mathbf{f}(\mathbf{x}) \), with \( \mathbf{x} \in \mathbb{R}^n \) and \( \mathbf{f} \) continuously differentiable. An equilibrium (operating point) \( \mathbf{x}^\star \) satisfies \( \mathbf{f}(\mathbf{x}^\star)=\mathbf{0} \). Linearization replaces the nonlinear vector field by its first-order Taylor model near \( \mathbf{x}^\star \).

The key practical question is: When does the local linear model correctly predict qualitative behavior (stability, trajectories, phase portrait) of the nonlinear system? The rigorous answer depends on whether the equilibrium is hyperbolic (no eigenvalues on the imaginary axis).

flowchart TD
  S["Start: nonlinear model xdot = f(x)"] --> E["Find equilibrium x_star, f(x_star)=0"]
  E --> J["Compute Jacobian A = df/dx at x_star"]
  J --> H["Check eigenvalues of A"]
  H --> HG["Case 1: hyperbolic, Re(lambda) not 0"]
  H --> NH["Case 2: nonhyperbolic, some Re(lambda)=0"]
  HG --> R["Estimate valid neighborhood \nand region of attraction"]
  NH --> TOOLS["Use nonlinear tools: \nLyapunov, invariants, normal forms"]
  R --> END["Use linear model only inside verified neighborhood"]
  TOOLS --> END
  

We will make each box above mathematically explicit and then show failure modes: (i) inconclusive linearization (e.g., eigenvalue 0); (ii) qualitative mismatch (e.g., linear center but nonlinear spiral); (iii) local vs global mismatch (stable locally, unstable globally due to other equilibria).

Let \( \mathbf{f} \) be \( C^1 \) near \( \mathbf{x}^\star \). Define perturbation \( \boldsymbol{\delta} = \mathbf{x} - \mathbf{x}^\star \). By Taylor's theorem,

\[ \mathbf{f}(\mathbf{x}^\star + \boldsymbol{\delta}) = \mathbf{f}(\mathbf{x}^\star) + \mathbf{A}\,\boldsymbol{\delta} + \mathbf{r}(\boldsymbol{\delta}), \quad \mathbf{A} = \left.\frac{\partial \mathbf{f} }{\partial \mathbf{x} }\right|_{\mathbf{x}^\star}, \quad \lim_{\|\boldsymbol{\delta}\|\to 0}\frac{\|\mathbf{r}(\boldsymbol{\delta})\|}{\|\boldsymbol{\delta}\|}=0. \]

Since \( \mathbf{f}(\mathbf{x}^\star)=\mathbf{0} \), the exact perturbation dynamics are \( \dot{\boldsymbol{\delta} } = \mathbf{A}\boldsymbol{\delta} + \mathbf{r}(\boldsymbol{\delta}) \). The linearized model is \( \dot{\boldsymbol{\delta} } = \mathbf{A}\boldsymbol{\delta} \).

A common (and dangerous) implicit assumption is that “small remainder” means “correct conclusions”. In fact:

• If the linearized equilibrium is hyperbolic, the nonlinear flow is locally topologically equivalent to the linear flow (Hartman–Grobman). Qualitative features (stable vs unstable manifold dimensions) are preserved.
• If the linearized equilibrium is nonhyperbolic, the remainder term can dominate the qualitative behavior no matter how small, because the linear part is too weak along neutral directions.

A useful analytic bound comes from a local Lipschitz estimate on the Jacobian: if \( \mathbf{f} \in C^2 \), then for small \( \boldsymbol{\delta} \), \( \|\mathbf{r}(\boldsymbol{\delta})\| \le c\,\|\boldsymbol{\delta}\|^2 \) for some constant \( c \). This quadratic smallness is powerful only when the linear part has exponential contraction/expansion.

An equilibrium \( \mathbf{x}^\star \) is hyperbolic if the Jacobian \( \mathbf{A} \) has no eigenvalues with zero real part: \( \Re(\lambda_i) \neq 0 \) for all eigenvalues. This creates an exponential separation between stable and unstable directions.

3.1 Hartman–Grobman theorem (statement)

For \( \dot{\mathbf{x} } = \mathbf{f}(\mathbf{x}) \) with \( \mathbf{f}\in C^1 \) and a hyperbolic equilibrium \( \mathbf{x}^\star \), there exists a neighborhood \( \mathcal{U} \) of \( \mathbf{x}^\star \) and a homeomorphism \( \mathbf{h}:\mathcal{U}\to \mathbb{R}^n \) such that trajectories of the nonlinear system in \( \mathcal{U} \) map to trajectories of the linear system \( \dot{\boldsymbol{\delta} }=\mathbf{A}\boldsymbol{\delta} \). In other words, the two flows are topologically conjugate locally.

Consequence: local qualitative properties that depend only on the phase portrait topology are preserved: dimensions of stable/unstable manifolds, existence of exponential convergence/divergence, and saddle structure. However, metric quantities (exact decay rates far from the equilibrium, overshoot, nonlinear distortions) need not match even in the hyperbolic case unless the neighborhood is verified.

3.2 Lyapunov’s indirect method (hyperbolic case)

If all eigenvalues satisfy \( \Re(\lambda_i) < 0 \), the equilibrium is locally asymptotically stable. If any eigenvalue satisfies \( \Re(\lambda_i) > 0 \), it is unstable. This is the practical “indirect method” conclusion, but note the hidden assumption: hyperbolicity.

\[ \text{If } \max_i \Re(\lambda_i) < 0 \text{ then } \exists\,\epsilon > 0\text{ s.t. }\|\mathbf{x}(0)-\mathbf{x}^\star\| < \epsilon \Rightarrow \mathbf{x}(t) \to \mathbf{x}^\star \text{ as } t\to\infty. \]

The proof idea (sketch): write perturbation dynamics \( \dot{\boldsymbol{\delta} } = \mathbf{A}\boldsymbol{\delta}+\mathbf{r}(\boldsymbol{\delta}) \). Use the exponential stability of \( \dot{\boldsymbol{\delta} }=\mathbf{A}\boldsymbol{\delta} \) and the smallness of \( \mathbf{r} \) to show that trajectories remain in a neighborhood and decay (via variation of constants and Grönwall-type arguments).

If \( \mathbf{A} \) has an eigenvalue with \( \Re(\lambda)=0 \) (e.g., \( \lambda=0 \) or purely imaginary), the equilibrium is nonhyperbolic. The linear model then contains neutral directions, and the higher-order terms decide stability.

4.1 Example: stable nonlinear system with zero Jacobian

Consider the scalar nonlinear system \( \dot{x} = -x^3 \). The equilibrium \( x^\star=0 \) has Jacobian \( A = f'(0)=0 \). The linearized system is \( \dot{\delta}=0 \), which is neither asymptotically stable nor unstable. Yet the true nonlinear system is asymptotically stable.

\[ \dot{x} = -x^3,\quad x^\star=0,\quad f'(0)=0. \]

A direct proof uses the Lyapunov candidate \( V(x)=\tfrac{1}{2}x^2 \), which is positive definite. Along trajectories:

\[ \dot{V}(x)=\frac{\mathrm{d} }{\mathrm{d}t}\left(\tfrac{1}{2}x^2\right)=x\dot{x}=x(-x^3)=-x^4 \le 0. \]

Moreover, \( \dot{V}=0 \) only at \( x=0 \), and the dynamics strictly decrease \( |x| \) for any nonzero state. One can integrate explicitly: \( x(t)=\frac{x(0)}{\sqrt{1+2x(0)^2 t} } \to 0 \). Thus the nonlinear equilibrium is asymptotically stable even though the linearization is inconclusive.

4.2 Example: linear center vs nonlinear spiral

Consider the planar system \( \dot{x}=-y + x(x^2+y^2),\; \dot{y}=x + y(x^2+y^2) \). The Jacobian at the origin is \( \mathbf{A}=\begin{bmatrix}0 & -1 \\ 1 & 0\end{bmatrix} \), whose eigenvalues are purely imaginary \( \pm j \). The linearized system is a center (closed orbits), but the nonlinear system spirals outward.

\[ \dot{x}=-y + x(x^2+y^2),\qquad \dot{y}=x + y(x^2+y^2). \]

Convert to polar coordinates \( x=r\cos\theta,\; y=r\sin\theta \). Using identities \( \dot{r} = \frac{x\dot{x}+y\dot{y} }{r} \) and \( \dot{\theta} = \frac{x\dot{y}-y\dot{x} }{r^2} \), we obtain:

\[ \dot{r} = r^3,\qquad \dot{\theta}=1. \]

Therefore \( r(t) \) strictly increases for any \( r(0)\neq 0 \) and the origin is unstable. The linearized center is qualitatively wrong: it predicts neutral closed orbits, while the nonlinear system escapes.

Even when the equilibrium is hyperbolic and locally stable, linearization says nothing about what happens far away. A common failure is assuming global stability from local eigenvalues.

5.1 Example: locally stable equilibrium with finite region of attraction

Consider the scalar system \( \dot{x}=-x+x^3 \). Equilibria are \( x^\star\in\{-1,0,1\} \). The Jacobian is \( f'(x)=-1+3x^2 \), so \( f'(0)=-1 \) and linearization predicts local asymptotic stability at \( 0 \). However, \( \pm 1 \) are unstable equilibria, and trajectories with \( |x(0)| > 1 \) diverge in magnitude.

\[ \dot{x}=-x+x^3 = x(x^2-1),\qquad f'(x)=-1+3x^2. \]

The phase line analysis is elementary: \( x(x^2-1) < 0 \) for \( 0 < |x| < 1 \) so trajectories move toward \( 0 \), while \( x(x^2-1) > 0 \) for \( |x| > 1 \) so trajectories move away. Hence the region of attraction of \( 0 \) is exactly \( (-1,1) \).

flowchart TD
  L["Local conclusion from A: stable near x_star"] --> Q1["Other equilibria exist"]
  Q1 --> Q2["Yes: analyze stability \nand basin boundaries"]
  Q1 --> Q3["No or unknown: check invariants \nor energy-like bounds"]
  Q2 --> Q4["Estimate region of attraction near x_star"]
  Q3 --> Q4
  Q4 --> SIM["Simulate nonlinear model for many initial states"]
  SIM --> USE["Use linear model only in verified neighborhood"]
  

5.2 Practical meaning for control engineering

• Linear controllers (e.g., state feedback, observer design) computed at an operating point guarantee local properties for the linear model. They must be validated on the nonlinear plant for initial conditions and disturbances within the intended envelope.
• A “good” linear model can still fail due to actuator limits, saturations, friction, and switching (addressed in Lesson 4), but even without those, nonlinear geometry can create multiple basins of attraction.

The following reference implementations reproduce three canonical failure modes: (A) nonhyperbolic but stable, (B) linear center but nonlinear spiral-out, (C) locally stable but not globally stable.

6.1 Python lab (from-scratch RK4 + Jacobian finite-differences)

Chapter14_Lesson3.py

"""
Chapter 14 - Nonlinear System Dynamics
Lesson 3 - Linearization vs. True Nonlinear Behavior: When Linear Models Fail

File: Chapter14_Lesson3.py
Dependencies: numpy, matplotlib (optional), scipy (optional)
This script:
  1) Computes Jacobian linearizations at equilibria
  2) Simulates nonlinear and linearized models (RK4 from scratch)
  3) Demonstrates cases where linearization is conclusive / inconclusive / misleading globally
"""

from __future__ import annotations
import math
from dataclasses import dataclass
from typing import Callable, List, Tuple

import numpy as np

try:
    import matplotlib.pyplot as plt
    HAS_MPL = True
except Exception:
    HAS_MPL = False


Vector = np.ndarray
Func = Callable[[float, Vector], Vector]


def rk4(f: Func, t0: float, tf: float, x0: Vector, h: float) -> Tuple[np.ndarray, np.ndarray]:
    """Classic fixed-step RK4 integrator (from scratch)."""
    n_steps = int(math.ceil((tf - t0) / h))
    ts = np.zeros(n_steps + 1, dtype=float)
    xs = np.zeros((n_steps + 1, x0.size), dtype=float)
    ts[0] = t0
    xs[0] = x0.copy()
    t = t0
    x = x0.copy()
    for k in range(n_steps):
        if t + h > tf:
            h = tf - t
        k1 = f(t, x)
        k2 = f(t + 0.5 * h, x + 0.5 * h * k1)
        k3 = f(t + 0.5 * h, x + 0.5 * h * k2)
        k4 = f(t + h, x + h * k3)
        x = x + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        t = t + h
        ts[k + 1] = t
        xs[k + 1] = x
    return ts, xs


def jacobian_fd(g: Callable[[Vector], Vector], x: Vector, eps: float = 1e-6) -> np.ndarray:
    """Finite-difference Jacobian J = dg/dx at x (central differences)."""
    x = np.asarray(x, dtype=float)
    n = x.size
    J = np.zeros((n, n), dtype=float)
    for i in range(n):
        dx = np.zeros(n, dtype=float)
        dx[i] = eps
        gp = g(x + dx)
        gm = g(x - dx)
        J[:, i] = (gp - gm) / (2.0 * eps)
    return J


# ---------------------------
# Example A: Nonhyperbolic equilibrium (linearization inconclusive)
# xdot = -x^3  (equilibrium at 0 is asymptotically stable)
# ---------------------------

def fA(t: float, x: Vector) -> Vector:
    return np.array([-x[0] ** 3], dtype=float)

def linA_at0(t: float, dx: Vector) -> Vector:
    # Jacobian A = 0 at x*=0
    return np.array([0.0], dtype=float)


# ---------------------------
# Example B: Linear center vs nonlinear spiral (linearization qualitatively wrong)
# xdot = -y + x (x^2 + y^2)
# ydot =  x + y (x^2 + y^2)
# At origin: A = [[0,-1],[1,0]] (pure rotation, center)
# Nonlinear term makes rdot = r^3 (spiral out)
# ---------------------------

def fB(t: float, z: Vector) -> Vector:
    x, y = z[0], z[1]
    r2 = x * x + y * y
    return np.array([-y + x * r2, x + y * r2], dtype=float)

A_B = np.array([[0.0, -1.0],
                [1.0,  0.0]], dtype=float)

def linB_at0(t: float, dz: Vector) -> Vector:
    return A_B @ dz


# ---------------------------
# Example C: Locally stable but not globally stable
# xdot = -x + x^3  (equilibria at -1,0,1; 0 locally stable, ROA is |x0|<1)
# ---------------------------

def fC(t: float, x: Vector) -> Vector:
    return np.array([-x[0] + x[0] ** 3], dtype=float)

def linC_at0(t: float, dx: Vector) -> Vector:
    # Jacobian A = -1 at x*=0
    return np.array([-dx[0]], dtype=float)


@dataclass
class SimResult:
    t: np.ndarray
    x_nl: np.ndarray
    x_lin: np.ndarray


def simulate_compare(f_nl: Func, f_lin: Func, xeq: Vector, x0: Vector,
                     t0: float = 0.0, tf: float = 20.0, h: float = 1e-3) -> SimResult:
    # Linear model is in perturbation coordinates: d = x - xeq
    d0 = x0 - xeq
    t_nl, x_nl = rk4(f_nl, t0, tf, x0, h)
    t_lin, d_lin = rk4(f_lin, t0, tf, d0, h)
    x_lin = d_lin + xeq.reshape(1, -1)
    # Ensure same time grid (RK4 uses same grid by construction)
    return SimResult(t=t_nl, x_nl=x_nl, x_lin=x_lin)


def main() -> None:
    np.set_printoptions(precision=4, suppress=True)

    # ----- Example A
    print("\nExample A: xdot = -x^3 (nonhyperbolic at 0, but stable)")
    xeq = np.array([0.0])
    x0 = np.array([0.8])
    resA = simulate_compare(fA, linA_at0, xeq, x0, tf=10.0, h=1e-3)
    print("x(T) nonlinear:", resA.x_nl[-1, 0], "  |  x(T) linear:", resA.x_lin[-1, 0])

    # ----- Example B
    print("\nExample B: linear center vs nonlinear spiral-out")
    z0 = np.array([0.2, 0.0])
    zeq = np.array([0.0, 0.0])
    resB = simulate_compare(fB, linB_at0, zeq, z0, tf=25.0, h=1e-3)
    r_nl = np.sqrt(resB.x_nl[:, 0] ** 2 + resB.x_nl[:, 1] ** 2)
    r_lin = np.sqrt(resB.x_lin[:, 0] ** 2 + resB.x_lin[:, 1] ** 2)
    print("r(0) nonlinear:", r_nl[0], " r(T) nonlinear:", r_nl[-1])
    print("r(0) linear   :", r_lin[0], " r(T) linear   :", r_lin[-1])

    # ----- Example C
    print("\nExample C: local vs global mismatch (region of attraction issue)")
    xeq = np.array([0.0])
    for x_init in [0.2, 0.9, 1.1, 1.5]:
        resC = simulate_compare(fC, linC_at0, xeq, np.array([x_init]), tf=10.0, h=1e-3)
        print(f"x0={x_init:.2f} -> x(T) nonlinear={resC.x_nl[-1,0]:.4f} | linear={resC.x_lin[-1,0]:.4f}")

    # ----- Jacobian verification (finite difference)
    print("\nFinite-difference Jacobian checks:")
    JB = jacobian_fd(lambda z: fB(0.0, z), np.array([0.0, 0.0]))
    print("J_B(0) =\n", JB)
    eigB = np.linalg.eigvals(JB)
    print("eig(J_B(0)) =", eigB)

    JC0 = jacobian_fd(lambda x: fC(0.0, x), np.array([0.0]))
    print("J_C(0) =", JC0[0,0], " expected -1")

    # ----- Optional plots
    if HAS_MPL:
        # Phase-plane comparison for Example B
        plt.figure()
        plt.plot(resB.x_lin[:, 0], resB.x_lin[:, 1], label="Linearized (center)")
        plt.plot(resB.x_nl[:, 0], resB.x_nl[:, 1], label="Nonlinear (spiral out)")
        plt.xlabel("x")
        plt.ylabel("y")
        plt.title("Example B: Phase portrait (linear vs nonlinear)")
        plt.axis("equal")
        plt.grid(True)
        plt.legend()
        plt.savefig("Chapter14_Lesson3_ExampleB_Phase.png", dpi=150)

        # Radius vs time for Example B
        plt.figure()
        plt.plot(resB.t, r_lin, label="r(t) linear")
        plt.plot(resB.t, r_nl, label="r(t) nonlinear")
        plt.xlabel("t")
        plt.ylabel("r")
        plt.title("Example B: Radius growth (linear vs nonlinear)")
        plt.grid(True)
        plt.legend()
        plt.savefig("Chapter14_Lesson3_ExampleB_Radius.png", dpi=150)

        print("\nSaved plots: Chapter14_Lesson3_ExampleB_Phase.png, Chapter14_Lesson3_ExampleB_Radius.png")
    else:
        print("\nmatplotlib not available; skipping plots.")

    # Save CSV for Example B
    out = np.column_stack([resB.t, resB.x_lin[:,0], resB.x_lin[:,1], resB.x_nl[:,0], resB.x_nl[:,1]])
    np.savetxt("Chapter14_Lesson3_ExampleB.csv", out, delimiter=",",
               header="t,x_lin,y_lin,x_nl,y_nl", comments="")
    print("Saved CSV: Chapter14_Lesson3_ExampleB.csv")


if __name__ == "__main__":
    main()

A companion exercise estimates the region of attraction numerically for \( \dot{x}=-x+x^3 \) by sampling initial conditions and classifying convergence to \( 0 \).

Chapter14_Lesson3_Ex1.py

"""
Chapter 14 - Nonlinear System Dynamics
Lesson 3 - Exercise 1

File: Chapter14_Lesson3_Ex1.py
Goal: Empirically estimate a region of attraction for xdot = -x + x^3 around x*=0
by simulating many initial conditions and classifying convergence to 0.

Dependencies: numpy
Optional: matplotlib
"""

from __future__ import annotations
import numpy as np

try:
    import matplotlib.pyplot as plt
    HAS_MPL = True
except Exception:
    HAS_MPL = False


def rk4_scalar(f, t0, tf, x0, h):
    n_steps = int(np.ceil((tf - t0) / h))
    t = t0
    x = float(x0)
    xs = np.zeros(n_steps + 1)
    ts = np.zeros(n_steps + 1)
    xs[0] = x
    ts[0] = t
    for k in range(n_steps):
        if t + h > tf:
            h = tf - t
        k1 = f(t, x)
        k2 = f(t + 0.5*h, x + 0.5*h*k1)
        k3 = f(t + 0.5*h, x + 0.5*h*k2)
        k4 = f(t + h, x + h*k3)
        x = x + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)
        t = t + h
        xs[k+1] = x
        ts[k+1] = t
    return ts, xs


def f(t, x):
    return -x + x**3


def main():
    rng = np.random.default_rng(7)
    N = 400
    t0, tf, h = 0.0, 10.0, 1e-3
    eps = 1e-2

    x0s = rng.uniform(-1.6, 1.6, size=N)
    xT = np.zeros(N)
    converged = np.zeros(N, dtype=int)

    for i, x0 in enumerate(x0s):
        _, xs = rk4_scalar(f, t0, tf, x0, h)
        xT[i] = xs[-1]
        converged[i] = 1 if abs(xs[-1]) < eps else 0

    # Save results
    out = np.column_stack([x0s, xT, converged])
    np.savetxt("Chapter14_Lesson3_Ex1_ROA.csv", out, delimiter=",",
               header="x0,xT,converged_to_0", comments="")

    # Quick summary
    conv_rate = converged.mean()
    print("Convergence fraction (to |x(T)| < eps):", conv_rate)

    # Optional visualization
    if HAS_MPL:
        plt.figure()
        plt.scatter(x0s, xT, s=10)
        plt.xlabel("x0")
        plt.ylabel("x(T)")
        plt.title("Empirical convergence for xdot = -x + x^3")
        plt.grid(True)
        plt.savefig("Chapter14_Lesson3_Ex1_ROA.png", dpi=150)
        print("Saved: Chapter14_Lesson3_Ex1_ROA.png")

    print("Saved CSV: Chapter14_Lesson3_Ex1_ROA.csv")


if __name__ == "__main__":
    main()

6.2 C++ lab (CSV export for linear vs nonlinear comparison)

Chapter14_Lesson3.cpp

/*
Chapter 14 - Nonlinear System Dynamics
Lesson 3 - Linearization vs. True Nonlinear Behavior: When Linear Models Fail

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

This program simulates Example B:
  xdot = -y + x (x^2 + y^2)
  ydot =  x + y (x^2 + y^2)
and its linearization at the origin:
  zdot = A z,  A = [[0,-1],[1,0]]
Outputs CSV: Chapter14_Lesson3_ExampleB_cpp.csv
*/

#include <array>
#include <cmath>
#include <fstream>
#include <iostream>

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double, 2>, 2>;

static Vec2 add(const Vec2& a, const Vec2& b) { return {a[0] + b[0], a[1] + b[1]}; }
static Vec2 sub(const Vec2& a, const Vec2& b) { return {a[0] - b[0], a[1] - b[1]}; }
static Vec2 mul(double s, const Vec2& a) { return {s * a[0], s * a[1]}; }
static double norm2(const Vec2& a) { return std::sqrt(a[0] * a[0] + a[1] * a[1]); }

static Vec2 matvec(const Mat2& A, const Vec2& x) {
    return {A[0][0] * x[0] + A[0][1] * x[1],
            A[1][0] * x[0] + A[1][1] * x[1]};
}

static Vec2 f_nl(double /*t*/, const Vec2& z) {
    const double x = z[0], y = z[1];
    const double r2 = x * x + y * y;
    return {-y + x * r2, x + y * r2};
}

static Vec2 f_lin(double /*t*/, const Vec2& z) {
    static const Mat2 A = { { {0.0, -1.0}, {1.0, 0.0} } };
    return matvec(A, z);
}

template <typename F>
static Vec2 rk4_step(F f, double t, const Vec2& x, double h) {
    const Vec2 k1 = f(t, x);
    const Vec2 k2 = f(t + 0.5 * h, add(x, mul(0.5 * h, k1)));
    const Vec2 k3 = f(t + 0.5 * h, add(x, mul(0.5 * h, k2)));
    const Vec2 k4 = f(t + h, add(x, mul(h, k3)));
    Vec2 sum = add(k1, add(mul(2.0, k2), add(mul(2.0, k3), k4)));
    return add(x, mul(h / 6.0, sum));
}

int main() {
    const double t0 = 0.0, tf = 25.0, h = 1e-3;
    const int nSteps = static_cast<int>(std::ceil((tf - t0) / h));

    Vec2 z0 = {0.2, 0.0};

    double t = t0;
    Vec2 z_nl = z0;
    Vec2 z_lin = z0; // perturbation coords = state here (equilibrium at origin)

    std::ofstream ofs("Chapter14_Lesson3_ExampleB_cpp.csv");
    ofs << "t,x_lin,y_lin,r_lin,x_nl,y_nl,r_nl\n";

    for (int k = 0; k <= nSteps; ++k) {
        const double rlin = norm2(z_lin);
        const double rnl  = norm2(z_nl);
        ofs << t << "," << z_lin[0] << "," << z_lin[1] << "," << rlin << ","
            << z_nl[0] << "," << z_nl[1] << "," << rnl << "\n";

        if (k == nSteps) break;
        z_lin = rk4_step(f_lin, t, z_lin, h);
        z_nl  = rk4_step(f_nl,  t, z_nl,  h);
        t += h;
    }

    ofs.close();
    std::cout << "Saved: Chapter14_Lesson3_ExampleB_cpp.csv\n";
    std::cout << "Final r_lin=" << norm2(z_lin) << "  r_nl=" << norm2(z_nl) << "\n";
    return 0;
}

6.3 Java lab (CSV export for linear vs nonlinear comparison)

Chapter14_Lesson3.java

/*
Chapter 14 - Nonlinear System Dynamics
Lesson 3 - Linearization vs. True Nonlinear Behavior: When Linear Models Fail

File: Chapter14_Lesson3.java
Build & run:
  javac Chapter14_Lesson3.java
  java Chapter14_Lesson3

Outputs CSV: Chapter14_Lesson3_ExampleB_java.csv
*/

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

public class Chapter14_Lesson3 {

    // Example B nonlinear dynamics
    static double[] fNL(double t, double[] z) {
        double x = z[0], y = z[1];
        double r2 = x*x + y*y;
        return new double[] { -y + x*r2, x + y*r2 };
    }

    // Linearization at origin: zdot = A z, A=[[0,-1],[1,0]]
    static double[] fLIN(double t, double[] z) {
        double x = z[0], y = z[1];
        return new double[] { -y, x };
    }

    static double[] add(double[] a, double[] b) {
        return new double[] { a[0] + b[0], a[1] + b[1] };
    }

    static double[] mul(double s, double[] a) {
        return new double[] { s * a[0], s * a[1] };
    }

    static double norm2(double[] a) {
        return Math.sqrt(a[0]*a[0] + a[1]*a[1]);
    }

    interface F {
        double[] eval(double t, double[] x);
    }

    static double[] rk4Step(F f, double t, double[] x, double h) {
        double[] k1 = f.eval(t, x);
        double[] k2 = f.eval(t + 0.5*h, add(x, mul(0.5*h, k1)));
        double[] k3 = f.eval(t + 0.5*h, add(x, mul(0.5*h, k2)));
        double[] k4 = f.eval(t + h, add(x, mul(h, k3)));
        double[] sum = add(k1, add(mul(2.0, k2), add(mul(2.0, k3), k4)));
        return add(x, mul(h/6.0, sum));
    }

    public static void main(String[] args) throws IOException {
        double t0 = 0.0, tf = 25.0, h = 1e-3;
        int nSteps = (int)Math.ceil((tf - t0) / h);

        double[] z0 = new double[] { 0.2, 0.0 };
        double[] zNL = new double[] { z0[0], z0[1] };
        double[] zLIN = new double[] { z0[0], z0[1] };

        PrintWriter out = new PrintWriter(new FileWriter("Chapter14_Lesson3_ExampleB_java.csv"));
        out.println("t,x_lin,y_lin,r_lin,x_nl,y_nl,r_nl");

        double t = t0;
        for (int k = 0; k <= nSteps; k++) {
            double rlin = norm2(zLIN);
            double rnl  = norm2(zNL);
            out.printf(java.util.Locale.US, "%.6f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n",
                    t, zLIN[0], zLIN[1], rlin, zNL[0], zNL[1], rnl);

            if (k == nSteps) break;
            zLIN = rk4Step(Chapter14_Lesson3::fLIN, t, zLIN, h);
            zNL  = rk4Step(Chapter14_Lesson3::fNL,  t, zNL,  h);
            t += h;
        }

        out.close();
        System.out.println("Saved: Chapter14_Lesson3_ExampleB_java.csv");
        System.out.println("Final r_lin=" + norm2(zLIN) + "  r_nl=" + norm2(zNL));
    }
}

6.4 MATLAB/Simulink lab (ODE45 + optional programmatic Simulink build)

Chapter14_Lesson3.m

% Chapter 14 - Nonlinear System Dynamics
% Lesson 3 - Linearization vs. True Nonlinear Behavior: When Linear Models Fail
%
% File: Chapter14_Lesson3.m
%
% This script demonstrates Example B (2D) and Example C (1D):
%   Example B:
%     xdot = -y + x (x^2 + y^2)
%     ydot =  x + y (x^2 + y^2)
%     Linearization at origin: zdot = A z, A = [0 -1; 1 0]
%   Example C:
%     xdot = -x + x^3 (0 locally stable, but not globally)
%
% Optional: programmatically generate a Simulink model for Example B.

clear; clc;

%% Example B: Compare nonlinear vs linearized
A = [0 -1; 1 0];

f_nl = @(t,z) [ -z(2) + z(1)*(z(1)^2 + z(2)^2);
                 z(1) + z(2)*(z(1)^2 + z(2)^2) ];

f_lin = @(t,z) A*z;

tspan = [0 25];
z0 = [0.2; 0.0];

opts = odeset('RelTol',1e-9,'AbsTol',1e-11);

[t_nl, z_nl] = ode45(f_nl, tspan, z0, opts);
[t_li, z_li] = ode45(f_lin, tspan, z0, opts);

r_nl = sqrt(z_nl(:,1).^2 + z_nl(:,2).^2);
r_li = sqrt(z_li(:,1).^2 + z_li(:,2).^2);

fprintf('Example B: r(0)=%.4f, r(T) nonlinear=%.4f, r(T) linear=%.4f\n', r_nl(1), r_nl(end), r_li(end));

figure; plot(z_li(:,1), z_li(:,2), 'LineWidth', 1.5); hold on;
plot(z_nl(:,1), z_nl(:,2), 'LineWidth', 1.5);
axis equal; grid on; xlabel('x'); ylabel('y');
title('Example B: Phase portrait (linear vs nonlinear)');
legend('Linearized','Nonlinear');
saveas(gcf, 'Chapter14_Lesson3_ExampleB_Phase_matlab.png');

figure; plot(t_li, r_li, 'LineWidth', 1.5); hold on;
plot(t_nl, r_nl, 'LineWidth', 1.5);
grid on; xlabel('t'); ylabel('r');
title('Example B: Radius growth (linear vs nonlinear)');
legend('r(t) linear','r(t) nonlinear');
saveas(gcf, 'Chapter14_Lesson3_ExampleB_Radius_matlab.png');

T = min(length(t_li), length(t_nl));
csvwrite('Chapter14_Lesson3_ExampleB_matlab.csv', [t_li(1:T), z_li(1:T,:), z_nl(1:T,:)]);

%% Example C: Region of attraction issue
fC = @(t,x) -x + x.^3;
fC_lin = @(t,dx) -dx;

x0s = [0.2, 0.9, 1.1, 1.5];
for k = 1:length(x0s)
    x0 = x0s(k);
    [tC, xC] = ode45(fC, [0 10], x0, opts);
    [tL, xL] = ode45(fC_lin, [0 10], x0, opts);
    fprintf('Example C: x0=%.2f -> x(T) nonlinear=%.4f | linear=%.4f\n', x0, xC(end), xL(end));
end

%% Optional Simulink model generation
buildSimulink = false;   % set true if you want to generate a Simulink model
if buildSimulink
    mdl = 'Chapter14_Lesson3_Simulink_ExampleB';
    if bdIsLoaded(mdl)
        close_system(mdl, 0);
    end
    new_system(mdl); open_system(mdl);

    % Integrators for x and y
    add_block('simulink/Continuous/Integrator', [mdl '/Int_x'], 'Position', [140 90 170 120]);
    add_block('simulink/Continuous/Integrator', [mdl '/Int_y'], 'Position', [140 170 170 200]);

    % Blocks to compute r2 = x^2 + y^2
    add_block('simulink/Math Operations/Product', [mdl '/x_sq'], 'Position', [260 90 290 120]);
    add_block('simulink/Math Operations/Product', [mdl '/y_sq'], 'Position', [260 170 290 200]);
    add_block('simulink/Math Operations/Sum', [mdl '/sum_r2'], 'Inputs', '++', 'Position', [330 130 360 160]);

    % Products for x*r2 and y*r2
    add_block('simulink/Math Operations/Product', [mdl '/x_r2'], 'Position', [410 90 440 120]);
    add_block('simulink/Math Operations/Product', [mdl '/y_r2'], 'Position', [410 170 440 200]);

    % Sums for xdot = -y + x*r2 and ydot = x + y*r2
    add_block('simulink/Math Operations/Gain', [mdl '/neg_y'], 'Gain', '-1', 'Position', [500 170 530 200]);
    add_block('simulink/Math Operations/Sum', [mdl '/sum_xdot'], 'Inputs', '++', 'Position', [560 110 590 140]);
    add_block('simulink/Math Operations/Sum', [mdl '/sum_ydot'], 'Inputs', '++', 'Position', [560 190 590 220]);

    % Scope
    add_block('simulink/Sinks/Scope', [mdl '/Scope'], 'Position', [700 120 730 190]);

    % Wiring
    add_line(mdl, 'Int_x/1', 'x_sq/1'); add_line(mdl, 'Int_x/1', 'x_sq/2');
    add_line(mdl, 'Int_y/1', 'y_sq/1'); add_line(mdl, 'Int_y/1', 'y_sq/2');
    add_line(mdl, 'x_sq/1', 'sum_r2/1'); add_line(mdl, 'y_sq/1', 'sum_r2/2');

    add_line(mdl, 'Int_x/1', 'x_r2/1'); add_line(mdl, 'sum_r2/1', 'x_r2/2');
    add_line(mdl, 'Int_y/1', 'y_r2/1'); add_line(mdl, 'sum_r2/1', 'y_r2/2');

    add_line(mdl, 'Int_y/1', 'neg_y/1');
    add_line(mdl, 'neg_y/1', 'sum_xdot/1');
    add_line(mdl, 'x_r2/1', 'sum_xdot/2');

    add_line(mdl, 'Int_x/1', 'sum_ydot/1');
    add_line(mdl, 'y_r2/1', 'sum_ydot/2');

    add_line(mdl, 'sum_xdot/1', 'Int_x/1');
    add_line(mdl, 'sum_ydot/1', 'Int_y/1');

    add_line(mdl, 'Int_x/1', 'Scope/1');
    add_line(mdl, 'Int_y/1', 'Scope/2');

    set_param(mdl, 'StopTime', '25', 'Solver', 'ode45');
    save_system(mdl);
    fprintf('Saved Simulink model: %s.slx\n', mdl);
end

6.5 Wolfram Mathematica lab (NDSolve + phase portrait)

Chapter14_Lesson3.nb

(*
Chapter 14 - Nonlinear System Dynamics
Lesson 3 - Linearization vs. True Nonlinear Behavior: When Linear Models Fail

File: Chapter14_Lesson3.nb

This is a plain-text Wolfram Notebook expression. Open it in Mathematica.
It compares nonlinear vs linearized dynamics for Example B and illustrates
a nonhyperbolic case in Example A.
*)

Notebook[{
  Cell["Chapter 14 - Nonlinear System Dynamics\nLesson 3 - Linearization vs. True Nonlinear Behavior: When Linear Models Fail", "Title"],

  Cell["Example B: linear center vs nonlinear spiral-out", "Section"],

  Cell[BoxData[
    ToBoxes[
      Module[{mu = 1, fNL, A, fLIN, solNL, solLIN, tmax = 25, z0 = {0.2, 0.0} },
        fNL[{x_, y_}] := {-y + x (x^2 + y^2), x + y (x^2 + y^2)};
        A = { {0, -1}, {1, 0} };
        fLIN[z_] := A.z;

        solNL = NDSolveValue[{z'[t] == fNL[z[t]], z[0] == z0}, z, {t, 0, tmax}];
        solLIN = NDSolveValue[{w'[t] == fLIN[w[t]], w[0] == z0}, w, {t, 0, tmax}];

        Print["r(T) nonlinear = ", Norm[solNL[tmax]]];
        Print["r(T) linear    = ", Norm[solLIN[tmax]]];

        Show[
          ParametricPlot[{solLIN[t][[1]], solLIN[t][[2]]}, {t, 0, tmax},
            PlotLegends -> {"Linearized"}, AxesLabel -> {"x", "y"}],
          ParametricPlot[{solNL[t][[1]], solNL[t][[2]]}, {t, 0, tmax},
            PlotLegends -> {"Nonlinear"}, AxesLabel -> {"x", "y"}]
        ]
      ]
    ]
  ], "Input"],

  Cell["Example A: nonhyperbolic equilibrium (linearization inconclusive)", "Section"],

  Cell[BoxData[
    ToBoxes[
      Module[{f, sol, x0 = 0.8, tmax = 10},
        f[x_] := -x^3;
        sol = NDSolveValue[{x'[t] == f[x[t]], x[0] == x0}, x, {t, 0, tmax}];
        Plot[sol[t], {t, 0, tmax}, AxesLabel -> {"t", "x(t)"}]
      ]
    ]
  ], "Input"]
}]

6.6 Notes on libraries in System Dynamics workflows

• Python: scipy.integrate.solve_ivp (IVP solvers), control (linear control tools), sympy (symbolic Jacobians), numpy/matplotlib (numerics/plots).
• MATLAB: ode45/ode15s and Simulink integrator blocks for state-space modeling.
• C++/Java: RK4 implemented explicitly for transparency; professional workflows may use established ODE libraries.

Problem 1 (Hyperbolicity and exponential stability): Let \( \dot{\boldsymbol{\delta} }=\mathbf{A}\boldsymbol{\delta} \) with all eigenvalues satisfying \( \Re(\lambda_i) < 0 \). Prove there exist constants \( M>0 \), \( \alpha>0 \) such that \( \|\boldsymbol{\delta}(t)\| \le M e^{-\alpha t}\|\boldsymbol{\delta}(0)\| \).

Solution: Since \( \mathbf{A} \) is Hurwitz, for any \( \mathbf{Q}\succ 0 \) there exists \( \mathbf{P}\succ 0 \) solving \( \mathbf{A}^\top \mathbf{P} + \mathbf{P}\mathbf{A} = -\mathbf{Q} \). Using \( V(\boldsymbol{\delta})=\boldsymbol{\delta}^\top \mathbf{P}\boldsymbol{\delta} \), we get \( \dot{V}=-\boldsymbol{\delta}^\top \mathbf{Q}\boldsymbol{\delta} \le -\lambda_{\min}(\mathbf{Q})\|\boldsymbol{\delta}\|^2 \). Combine with \( \lambda_{\min}(\mathbf{P})\|\boldsymbol{\delta}\|^2 \le V \le \lambda_{\max}(\mathbf{P})\|\boldsymbol{\delta}\|^2 \) to obtain \( V(t)\le e^{-ct}V(0) \) with \( c=\lambda_{\min}(\mathbf{Q})/\lambda_{\max}(\mathbf{P}) \), implying exponential decay.

Problem 2 (Nonhyperbolic but asymptotically stable): Prove that \( x^\star=0 \) of \( \dot{x}=-x^3 \) is globally asymptotically stable.

Solution: Separate variables: \( \mathrm{d}x/x^3 = -\mathrm{d}t \).

\[ -\frac{1}{2x(t)^2} + \frac{1}{2x(0)^2} = -t \quad\Rightarrow\quad x(t)=\frac{x(0)}{\sqrt{1+2x(0)^2 t} } \to 0. \]

Problem 3 (Polar reduction): For \( \dot{x}=-y + x(x^2+y^2) \), \( \dot{y}=x + y(x^2+y^2) \), derive \( \dot{r}=r^3 \) and \( \dot{\theta}=1 \).

Solution: Compute \( x\dot{x}+y\dot{y} = r^4 \), hence \( \dot{r}=(x\dot{x}+y\dot{y})/r=r^3 \). Similarly \( x\dot{y}-y\dot{x} = r^2 \), hence \( \dot{\theta}=(x\dot{y}-y\dot{x})/r^2=1 \).

Problem 4 (Region of attraction): For \( \dot{x}=-x+x^3 \), prove the region of attraction of \( 0 \) is \( (-1,1) \).

Solution: Equilibria are \( -1,0,1 \). For \( 0<x<1 \), \( \dot{x}=x(x^2-1)<0 \) so \( x(t) \) decreases toward 0. For \( -1<x<0 \), \( \dot{x}>0 \) so \( x(t) \) increases toward 0. For \( |x|>1 \), \( \dot{x} \) pushes states away to infinity. Thus the basin is exactly \( (-1,1) \).

Problem 5 (Nonlinear perturbation bound): Suppose \( \mathbf{A} \) is Hurwitz and \( \|\mathbf{r}(\boldsymbol{\delta})\|\le c\|\boldsymbol{\delta}\|^2 \) locally. Show that sufficiently small \( \|\boldsymbol{\delta}(0)\| \) implies exponential decay of \( \|\boldsymbol{\delta}(t)\| \).

Solution (sketch): Use \( \boldsymbol{\delta}(t)=e^{\mathbf{A}t}\boldsymbol{\delta}(0)+\int_0^t e^{\mathbf{A}(t-s)}\mathbf{r}(\boldsymbol{\delta}(s))\,\mathrm{d}s \), the bound \( \|e^{\mathbf{A}t}\|\le M e^{-\alpha t} \), and bootstrap: assume \( \|\boldsymbol{\delta}(t)\|\le \rho \), then \( \|\boldsymbol{\delta}(t)\|\le M e^{-\alpha t}\|\boldsymbol{\delta}(0)\| + Mc\int_0^t e^{-\alpha(t-s)}\|\boldsymbol{\delta}(s)\|^2\,\mathrm{d}s \). For small enough initial condition, the quadratic integral term remains dominated and the inequality closes.

• Linearization is a first-order local model: \( \dot{\boldsymbol{\delta} }=\mathbf{A}\boldsymbol{\delta} \). It is reliable for qualitative local behavior only at hyperbolic equilibria.
• If any eigenvalue has \( \Re(\lambda)=0 \), linearization can be silent or misleading; higher-order terms decide the outcome.
• Local stability does not imply global stability; the region of attraction must be analyzed or estimated.
• In control engineering practice, validate linear designs on the nonlinear plant inside the operating envelope.