Consider an autonomous nonlinear state model (introduced in Chapter 8) with an equilibrium: \( \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}),\;\mathbf{x}\in\mathbb{R}^n \), and an equilibrium point \( \mathbf{x}_e \) satisfying \( \mathbf{f}(\mathbf{x}_e)=\mathbf{0} \). Define the perturbation (from Lesson 1): \( \boldsymbol{\eta} \triangleq \mathbf{x}-\mathbf{x}_e \).

If \( \mathbf{f} \) is continuously differentiable in a neighborhood of \( \mathbf{x}_e \), then the Jacobian linearization (Lesson 2) is \( \dot{\boldsymbol{\eta}} = \mathbf{A}\boldsymbol{\eta} \), where \( \mathbf{A} = \left.\dfrac{\partial \mathbf{f}}{\partial \mathbf{x}}\right|_{\mathbf{x}_e} \). The key question here is: for which initial perturbations does the nonlinear trajectory behave “like” the linear one, and for which initial perturbations does it still converge to the same equilibrium?

1.1 Decomposition into linear part + nonlinear remainder

Translate the equilibrium to the origin by redefining coordinates so that \( \mathbf{x}_e=\mathbf{0} \). (This is purely notational; the results apply around any \( \mathbf{x}_e \).) We can always write

\[ \mathbf{f}(\mathbf{x}) = \mathbf{A}\mathbf{x} + \mathbf{r}(\mathbf{x}), \quad \mathbf{A} = \left.\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\right|_{\mathbf{0}}, \quad \mathbf{r}(\mathbf{0})=\mathbf{0},\; \left.\frac{\partial \mathbf{r}}{\partial \mathbf{x}}\right|_{\mathbf{0}}=\mathbf{0}. \]

The remainder term \( \mathbf{r}(\mathbf{x}) \) contains all higher-order effects. “Validity of linearization” is governed by how small \( \mathbf{r}(\mathbf{x}) \) is relative to \( \mathbf{A}\mathbf{x} \) in the region being considered.

1.2 Second-order bound (core fact behind “local”)

A standard sufficient condition (easy to check in practice) is that the Jacobian is locally Lipschitz. Assume there exists \( L_J > 0 \) and \( \rho > 0 \) such that for all \( \mathbf{z} \) with \( \|\mathbf{z}\| < \rho \),

\[ \left\| \mathbf{J}(\mathbf{z}) - \mathbf{J}(\mathbf{0}) \right\| \le L_J \|\mathbf{z}\|, \quad \text{where } \mathbf{J}(\mathbf{x}) \triangleq \frac{\partial \mathbf{f}}{\partial \mathbf{x}}(\mathbf{x}). \]

Then the nonlinear remainder is quadratically small near the equilibrium:

\[ \|\mathbf{r}(\mathbf{x})\| \le \frac{L_J}{2}\|\mathbf{x}\|^2, \quad \forall\, \|\mathbf{x}\| < \rho. \]

Proof (integral form of Taylor’s theorem):

Since \( \mathbf{f} \) is differentiable, for any \( \mathbf{x} \) we have

\[ \mathbf{f}(\mathbf{x}) - \mathbf{f}(\mathbf{0}) = \int_{0}^{1} \mathbf{J}(s\mathbf{x})\,\mathbf{x}\,ds. \]

With \( \mathbf{f}(\mathbf{0})=\mathbf{0} \) and \( \mathbf{A}=\mathbf{J}(\mathbf{0}) \), define \( \mathbf{r}(\mathbf{x}) \triangleq \mathbf{f}(\mathbf{x})-\mathbf{A}\mathbf{x} \). Then

\[ \mathbf{r}(\mathbf{x}) = \int_{0}^{1} \left(\mathbf{J}(s\mathbf{x})-\mathbf{J}(\mathbf{0})\right)\mathbf{x}\,ds. \]

Taking norms and using the Lipschitz bound,

\[ \|\mathbf{r}(\mathbf{x})\| \le \int_{0}^{1} \left\|\mathbf{J}(s\mathbf{x})-\mathbf{J}(\mathbf{0})\right\|\,\|\mathbf{x}\|\,ds \le \int_{0}^{1} L_J \|s\mathbf{x}\|\,\|\mathbf{x}\|\,ds = \frac{L_J}{2}\|\mathbf{x}\|^2. \]

This inequality is the mathematical reason why linearization is inherently local: the neglected term scales like \( \|\mathbf{x}\|^2 \), so it becomes non-negligible as the state grows.

flowchart TD
  A["Pick equilibrium x_e"] --> B["Shift coordinates: eta = x - x_e"]
  B --> C["Compute Jacobian A = d f / d x at x_e"]
  C --> D["Write nonlinear dynamics: eta_dot = A eta + r(eta)"]
  D --> E["Bound remainder: ||r(eta)|| <= c ||eta||^2 in a ball"]
  E --> F["Use bound to certify: inner region where linear model dominates"]
  F --> G["Simulate nonlinear system to validate (optional)"]
        

For the equilibrium at the origin, define the (maximal) region of attraction as the set of initial conditions that converge to the equilibrium under the nonlinear dynamics.

\[ \mathcal{R} \triangleq \left\{ \mathbf{x}_0 \in \mathbb{R}^n \;:\; \mathbf{x}(0)=\mathbf{x}_0,\; \text{and } \mathbf{x}(t)\ \text{exists for all } t\ge 0,\; \mathbf{x}(t)\; →\; \mathbf{0}\ \text{as } t\; →\; \infty \right\}. \]

The region of attraction is typically hard to compute exactly for nonlinear systems. However, several properties can be stated at the level needed in this chapter:

• Invariance: If \( \mathbf{x}_0 \in \mathcal{R} \), then \( \mathbf{x}(t;\mathbf{x}_0)\in \mathcal{R} \) for all \( t\ge 0 \).
• Openness (under smoothness): If the equilibrium is asymptotically stable and \( \mathbf{f} \) is locally Lipschitz, then \( \mathcal{R} \) is an open set.
• Non-globality: Even if the linearization suggests stability, \( \mathcal{R} \) may be a strict subset of \( \mathbb{R}^n \) because other equilibria or nonlinear phenomena may exist.

2.1 A canonical example where linearization over-promises

Consider the scalar nonlinear system: \( \dot{x} = -x + x^3 \). The equilibrium \( x=0 \) has linearization \( \dot{x} = -x \), which is stable. But the nonlinear system has additional equilibria at \( x=\pm 1 \), and the region of attraction of the origin is only \( |x_0| < 1 \).

Derivation of the region of attraction (exact, 1D):

Factor the dynamics: \( \dot{x} = x(x^2-1) \). Then:

• If \( 0 < x < 1 \), then \( x^2-1 < 0 \) so \( \dot{x} < 0 \), hence \( x(t) \) decreases toward 0.
• If \( -1 < x < 0 \), then \( x < 0 \) and \( x^2-1 < 0 \), so \( \dot{x} > 0 \), hence \( x(t) \) increases toward 0.
• If \( x > 1 \), then \( x > 0 \) and \( x^2-1 > 0 \), so \( \dot{x} > 0 \) and the state moves away from 0.
• If \( x < -1 \), then \( x < 0 \) and \( x^2-1 > 0 \), so \( \dot{x} < 0 \) and the state moves away from 0.

Therefore: \( \mathcal{R} = (-1,1) \) for the equilibrium at 0, while the linear model would (incorrectly) suggest convergence for all \( x_0 \).

We now derive a rigorous (but conservative) method to obtain a guaranteed inner subset of \( \mathcal{R} \) using the linearized model. This is one of the main engineering uses of linearization: it can provide certificates of local convergence.

3.1 Assumptions and setup

Consider the shifted system around the equilibrium at the origin:

\[ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{r}(\mathbf{x}), \quad \|\mathbf{r}(\mathbf{x})\|\le c\|\mathbf{x}\|^2 \ \text{for}\ \|\mathbf{x}\| < \rho, \]

where \( c > 0 \) comes from the remainder bound (Section 1.2). Assume the linearized matrix \( \mathbf{A} \) is Hurwitz, meaning all eigenvalues satisfy \( \operatorname{Re}(\lambda_i(\mathbf{A})) < 0 \). (You already have the linear-algebra background to interpret eigenvalues from Chapter 2.)

3.2 Quadratic Lyapunov function induced by the linearization

Choose any symmetric positive definite matrix \( \mathbf{Q}\succ \mathbf{0} \). Since \( \mathbf{A} \) is Hurwitz, the continuous Lyapunov equation has a unique symmetric solution \( \mathbf{P}\succ \mathbf{0} \):

\[ \mathbf{A}^T\mathbf{P} + \mathbf{P}\mathbf{A} = -\mathbf{Q}. \]

Define \( V(\mathbf{x}) \triangleq \mathbf{x}^T\mathbf{P}\mathbf{x} \). Along trajectories of the nonlinear system:

\[ \dot{V}(\mathbf{x}) = \mathbf{x}^T(\mathbf{A}^T\mathbf{P}+\mathbf{P}\mathbf{A})\mathbf{x} + 2\mathbf{x}^T\mathbf{P}\mathbf{r}(\mathbf{x}) = -\mathbf{x}^T\mathbf{Q}\mathbf{x} + 2\mathbf{x}^T\mathbf{P}\mathbf{r}(\mathbf{x}). \]

Bound the nonlinear term using Cauchy–Schwarz and operator norms:

\[ 2\mathbf{x}^T\mathbf{P}\mathbf{r}(\mathbf{x}) \le 2\|\mathbf{x}\|\,\|\mathbf{P}\|\,\|\mathbf{r}(\mathbf{x})\| \le 2\|\mathbf{x}\|\,\|\mathbf{P}\|\,c\|\mathbf{x}\|^2 = 2c\|\mathbf{P}\|\|\mathbf{x}\|^3. \]

Also, since \( \mathbf{Q}\succ\mathbf{0} \), \( \mathbf{x}^T\mathbf{Q}\mathbf{x} \ge \lambda_{\min}(\mathbf{Q})\|\mathbf{x}\|^2 \). Therefore:

\[ \dot{V}(\mathbf{x}) \le -\lambda_{\min}(\mathbf{Q})\|\mathbf{x}\|^2 + 2c\|\mathbf{P}\|\|\mathbf{x}\|^3 = -\|\mathbf{x}\|^2\left(\lambda_{\min}(\mathbf{Q}) - 2c\|\mathbf{P}\|\|\mathbf{x}\|\right). \]

Hence, if \( \|\mathbf{x}\| < r \) where

\[ r \triangleq \min\!\left(\rho,\; \frac{\lambda_{\min}(\mathbf{Q})}{2c\|\mathbf{P}\|}\right), \]

then the bracket is positive and \( \dot{V}(\mathbf{x}) < 0 \) for all \( \mathbf{x}\neq \mathbf{0} \) in that ball. This yields a certified local convergence region.

3.3 Converting the radius bound into an invariant ellipsoid

Since \( V(\mathbf{x}) = \mathbf{x}^T\mathbf{P}\mathbf{x} \ge \lambda_{\min}(\mathbf{P})\|\mathbf{x}\|^2 \), the set

\[ \Omega_\alpha \triangleq \left\{\mathbf{x}:\; \mathbf{x}^T\mathbf{P}\mathbf{x} \le \alpha \right\} \]

satisfies \( \Omega_\alpha \subset \{\|\mathbf{x}\|\le \sqrt{\alpha/\lambda_{\min}(\mathbf{P})}\} \). Choose \( \alpha \) so that \( \sqrt{\alpha/\lambda_{\min}(\mathbf{P})} < r \), i.e.

\[ \alpha < \lambda_{\min}(\mathbf{P})\,r^2. \]

Then \( \dot{V} < 0 \) on \( \Omega_\alpha \setminus \{\mathbf{0}\} \), implying trajectories starting in \( \Omega_\alpha \) remain in it and converge to the equilibrium. This provides a conservative inner approximation \( \Omega_\alpha \subset \mathcal{R} \).

Conceptually, linearization provides: \( \mathbf{A} \) (local dynamics) → solve Lyapunov equation for \( \mathbf{P} \) → choose a level set of \( V \) that dominates higher-order nonlinearities.

4.1 Non-hyperbolic equilibria (zero/imaginary eigenvalues)

Linearization gives a definitive local stability conclusion only when the equilibrium is hyperbolic (no eigenvalues with zero real part). If \( \mathbf{A} \) has an eigenvalue with \( \operatorname{Re}(\lambda)=0 \), higher-order terms decide stability.

Example (scalar): \( \dot{x} = x^2 \). The Jacobian at 0 is \( A=0 \).

The linearized system \( \dot{x}=0 \) is inconclusive. But the nonlinear system has \( \dot{x} > 0 \) for all \( x \neq 0 \), so any small positive perturbation grows away. Thus the equilibrium is not asymptotically stable.

4.2 Multiple equilibria and finite regions of attraction

As in \( \dot{x}=-x+x^3 \), the presence of other equilibria can “cap” the region of attraction. Linearization around one equilibrium does not encode the global phase-space structure.

4.3 Large-signal behavior and neglected nonlinearities

The bound \( \|\mathbf{r}(\mathbf{x})\|\le c\|\mathbf{x}\|^2 \) is inherently local. When \( \|\mathbf{x}\| \) is not small, the neglected terms can:

• shift the effective stiffness/damping (state-dependent dynamics),
• create or destroy equilibria as parameters/operating points change,
• introduce behaviors not possible in linear models (e.g., self-sustained oscillations; treated later in Chapter 14).

4.4 Input-dependent validity (brief note)

If the system includes inputs \( \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u}) \) and you linearize around \( (\mathbf{x}_e,\mathbf{u}_e) \), then the approximation is valid only for perturbations \( \|\delta\mathbf{x}\| \) and \( \|\delta\mathbf{u}\| \) small enough that higher-order coupling terms remain dominated by the linear model: \( \delta\dot{\mathbf{x}} \approx \mathbf{A}\delta\mathbf{x}+\mathbf{B}\delta\mathbf{u} \). This is particularly relevant when actuators saturate (a nonlinearity not captured by a Jacobian).

flowchart TD
  X0["Initial condition x0"] --> Q1["Near equilibrium?"]
  Q1 -->|"yes"| L["Linear model predicts convergence"]
  Q1 -->|"no"| N["Nonlinear effects dominate"]
  L --> Q2["Other equilibria/attractors nearby?"]
  Q2 -->|"no"| OK1["Large basin likely \n(but still verify)"]
  Q2 -->|"yes"| BASIN["Finite basin: \nboundary/separatrix exists"]
  N --> RISK["Linearization unreliable: \nsimulate or use nonlinear analysis"]
  

We demonstrate a full workflow on a simple 2D example that exhibits a finite region of attraction in one coordinate:

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

Equilibrium at \( (x,y)=(0,0) \). Linearization yields \( \dot{x}=-x,\; \dot{y}=-2y \) (stable), but the true ROA is \( |x_0| < 1 \) with any \( y_0 \), because the \( x \)-subsystem matches the scalar counterexample from Section 2.1 while \( y \) always decays.

5.1 Python (SymPy + SciPy + Control-related tooling)

Suggested libraries: sympy (symbolic Jacobians), scipy.integrate (simulation), scipy.linalg (Lyapunov equation), and optionally control (python-control) for LTI checks.

import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp
from scipy.linalg import solve_continuous_lyapunov, eigvals, norm

# ----------------------------
# Define nonlinear system
# xdot = -x + x^3
# ydot = -2y
# ----------------------------
def f(t, z):
    x, y = z
    return np.array([-x + x**3, -2.0*y])

# Symbolic Jacobian at equilibrium
x, y = sp.symbols('x y', real=True)
f_sym = sp.Matrix([ -x + x**3, -2*y ])
J = f_sym.jacobian([x, y])
A = np.array(J.subs({x: 0, y: 0})).astype(float)  # [[-1,0],[0,-2]]

print("A =\n", A)
print("eig(A) =", eigvals(A))

# Lyapunov equation: A^T P + P A = -Q
Q = np.eye(2)
P = solve_continuous_lyapunov(A.T, -Q)  # solves A^T P + P A = -Q
print("P =\n", P)

# Conservative inner ellipsoid: x^T P x <= alpha
# Here remainder is r(x,y) = [x^3, 0], and one can bound ||r|| <= c ||z||^2 locally.
# For demonstration, pick alpha by inspection (must be small enough).
alpha = 0.5

def in_ellipsoid(z):
    return z.T @ P @ z <= alpha

# Monte-Carlo grid to approximate ROA numerically (simulation-based, not certified)
grid = np.linspace(-1.5, 1.5, 61)
Tfinal = 10.0
tol = 1e-3

def converges(z0):
    sol = solve_ivp(f, [0, Tfinal], z0, rtol=1e-9, atol=1e-12)
    zT = sol.y[:, -1]
    return np.linalg.norm(zT) < tol

count_conv = 0
count_total = 0
for xi in grid:
    for yi in grid:
        z0 = np.array([xi, yi])
        count_total += 1
        if converges(z0):
            count_conv += 1

print("Converged fraction (grid-based heuristic):", count_conv / count_total)

# Note: The true ROA is |x0| < 1 (all y0), but simulations illustrate the basin boundary.

5.2 C++ (Eigen + Boost.odeint; include “from scratch” Lyapunov via Kronecker)

Suggested libraries: Eigen (linear algebra), boost::numeric::odeint (ODE integration). For small matrices, you can solve the Lyapunov equation by vectorization:

\[ \mathbf{A}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{Q} \;\Longleftrightarrow\; \left(\mathbf{I}\otimes \mathbf{A}^T + \mathbf{A}^T\otimes \mathbf{I}\right)\operatorname{vec}(\mathbf{P}) = -\operatorname{vec}(\mathbf{Q}). \]

#include <iostream>
#include <Eigen/Dense>
#include <boost/numeric/odeint.hpp>

using Eigen::MatrixXd;
using Eigen::VectorXd;
using namespace boost::numeric::odeint;

// Nonlinear system: [x' = -x + x^3, y' = -2y]
struct System {
  void operator()(const std::array<double,2> &z,
                  std::array<double,2> &dzdt,
                  const double /*t*/) const {
    const double x = z[0];
    const double y = z[1];
    dzdt[0] = -x + x*x*x;
    dzdt[1] = -2.0*y;
  }
};

// Helper: vec operator for 2x2 (column-stacking)
VectorXd vec2(const MatrixXd &M) {
  VectorXd v(4);
  v << M(0,0), M(1,0), M(0,1), M(1,1);
  return v;
}

// Helper: unvec for 2x2
MatrixXd unvec2(const VectorXd &v) {
  MatrixXd M(2,2);
  M(0,0) = v(0); M(1,0) = v(1);
  M(0,1) = v(2); M(1,1) = v(3);
  return M;
}

// Kronecker product (small utility for 2x2)
MatrixXd kron(const MatrixXd &A, const MatrixXd &B) {
  MatrixXd K(A.rows()*B.rows(), A.cols()*B.cols());
  for(int i=0;i<A.rows();++i){
    for(int j=0;j<A.cols();++j){
      K.block(i*B.rows(), j*B.cols(), B.rows(), B.cols()) = A(i,j)*B;
    }
  }
  return K;
}

int main() {
  // Linearization matrix A
  MatrixXd A(2,2);
  A << -1, 0,
        0, -2;

  MatrixXd Q = MatrixXd::Identity(2,2);

  // Solve Lyapunov: A^T P + P A = -Q via vec
  MatrixXd I = MatrixXd::Identity(2,2);
  MatrixXd K = kron(I, A.transpose()) + kron(A.transpose(), I);
  VectorXd b = -vec2(Q);
  VectorXd pvec = K.fullPivLu().solve(b);
  MatrixXd P = unvec2(pvec);

  std::cout << "P =\n" << P << "\n";

  // Simulate a trajectory
  std::array<double,2> z = {0.8, 2.0};   // try x0 near basin boundary
  runge_kutta4< std::array<double,2> > stepper;

  double t = 0.0, dt = 1e-3, T = 10.0;
  System sys;
  for (int k=0; t < T; ++k) {
    stepper.do_step(sys, z, t, dt);
    t += dt;
  }
  std::cout << "z(T) = [" << z[0] << ", " << z[1] << "]\n";
  return 0;
}

5.3 Java (Apache Commons Math ODE + EJML for linear algebra)

Suggested libraries: org.apache.commons:commons-math3 (ODE integrators) and org.ejml:ejml-all (matrix algebra). For Lyapunov, either use a built-in routine if available, or implement the Kronecker/vectorization solve as in C++ (small dimensions).

import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;

public class LinearizationValidityDemo {

  // x' = -x + x^3, y' = -2y
  static class NonlinearSystem implements FirstOrderDifferentialEquations {
    @Override
    public int getDimension() { return 2; }

    @Override
    public void computeDerivatives(double t, double[] z, double[] dzdt) {
      double x = z[0];
      double y = z[1];
      dzdt[0] = -x + x*x*x;
      dzdt[1] = -2.0*y;
    }
  }

  public static void main(String[] args) {
    FirstOrderDifferentialEquations sys = new NonlinearSystem();

    double minStep = 1e-6;
    double maxStep = 1e-2;
    double absTol = 1e-10;
    double relTol = 1e-10;

    DormandPrince853Integrator integrator =
        new DormandPrince853Integrator(minStep, maxStep, absTol, relTol);

    double t0 = 0.0;
    double t1 = 10.0;
    double[] z = new double[] {0.8, 2.0}; // x0 near boundary, y0 large

    integrator.integrate(sys, t0, z, t1, z);

    System.out.println("State at t=10: x=" + z[0] + ", y=" + z[1]);

    // Tip: For ROA estimation, loop over a grid of initial conditions and test convergence numerically.
    // For certification, compute P from Lyapunov equation A^T P + P A = -Q (via Kronecker solve).
  }
}

5.4 MATLAB + Simulink (lyap, ode45, and linearization tools)

In MATLAB, the Lyapunov equation is solved by lyap. Nonlinear simulation is via ode45. In Simulink, you can build the nonlinear model with Integrator blocks and a MATLAB Function block for \( x^3 \), then use linearize (or the Model Linearizer app) around an operating point.

% Nonlinear system: xdot = -x + x^3, ydot = -2y
f = @(t,z) [-z(1) + z(1)^3; -2*z(2)];

% Linearization at origin: A = diag(-1,-2)
A = [-1 0; 0 -2];

% Lyapunov equation: A'P + P A = -Q
Q = eye(2);
P = lyap(A', Q);

disp('P ='); disp(P);

% Certified ellipsoid candidate: z' P z <= alpha (choose alpha small)
alpha = 0.5;

% Simulation from an initial condition
z0 = [0.8; 2.0];  % try near boundary
[t,z] = ode45(f, [0 10], z0);

zT = z(end,:)';
disp('State at t=10 ='); disp(zT);

% Grid-based basin approximation (simulation heuristic)
grid = linspace(-1.5, 1.5, 61);
tol = 1e-3;
Tfinal = 10;
conv = zeros(length(grid), length(grid));
for i=1:length(grid)
  for j=1:length(grid)
    z0 = [grid(i); grid(j)];
    [~,zz] = ode45(f, [0 Tfinal], z0);
    conv(i,j) = (norm(zz(end,:)') < tol);
  end
end
fprintf('Converged fraction: %g\n', mean(conv(:)));

Simulink implementation notes (no image, procedural):

• Use two Integrator blocks for \( x \) and \( y \).
• Compute \( -x \) and \( x^3 \) (Math Function block set to power or a MATLAB Function block).
• Sum to produce \( \dot{x} \); multiply \( y \) by -2 for \( \dot{y} \).
• For linearization: define operating point at \( x=0, y=0 \), then call linearize to obtain \( A \).

5.5 Wolfram Mathematica (symbolic Jacobian + LyapunovSolve + NDSolve)

(* Define the vector field *)
f[{x_, y_}] := { -x + x^3, -2 y };

(* Jacobian *)
J = D[f[{x, y}], {{x, y}}];
A = J /. {x -> 0, y -> 0} // N
(* A should be { {-1,0},{0,-2} } *)

(* Solve Lyapunov equation: Transpose[A].P + P.A == -Q *)
Q = IdentityMatrix[2];
P = LyapunovSolve[Transpose[A], Q];
P // MatrixForm

(* Simulate nonlinear system *)
sol = NDSolve[
  {x'[t] == -x[t] + x[t]^3, y'[t] == -2 y[t], x[0] == 0.8, y[0] == 2.0},
  {x, y}, {t, 0, 10}
];

{x[10], y[10]} /. sol[[1]]

5.6 Practical checklist (what you should do in engineering practice)

1. Compute \( \mathbf{A} \) and check eigenvalues. If any satisfy \( \operatorname{Re}(\lambda)\ge 0 \), local asymptotic stability is not guaranteed by linearization.
2. Estimate a remainder bound constant \( c \) (analytically for simple models; numerically by sampling Jacobian variation in a small ball).
3. Solve \( \mathbf{A}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{Q} \) for some \( \mathbf{Q}\succ 0 \).
4. Choose a level set \( \mathbf{x}^T\mathbf{P}\mathbf{x}\le \alpha \) small enough to dominate the nonlinear term, yielding a certified inner ROA.
5. Use simulation to validate conservatism and to explore the basin boundary (simulation is not a proof).

Problem 1 (Exact ROA in 1D): Consider \( \dot{x} = -x + x^3 \). (a) Find all equilibria. (b) Determine the region of attraction of the equilibrium at \( x=0 \). (c) Compare with the prediction from linearization at \( x=0 \).

Solution:

(a) Equilibria satisfy \( -x + x^3 = 0 \) so \( x(x^2-1)=0 \), hence \( x\in\{-1,0,1\} \).

(b) For \( 0 < x < 1 \), \( \dot{x} = x(x^2-1) < 0 \) so trajectories decrease to 0. For \( -1 < x < 0 \), \( \dot{x} > 0 \) so trajectories increase to 0. If \( x > 1 \), \( \dot{x} > 0 \) so the trajectory moves away from 0; similarly if \( x < -1 \), then \( \dot{x} < 0 \) moving away from 0. Therefore \( \mathcal{R} = (-1,1) \).

(c) Linearization at 0 yields \( \dot{x}=-x \), predicting convergence for all \( x_0 \). The nonlinear model shows the true basin is finite: linearization is valid only locally.

Problem 2 (Remainder bound via Lipschitz Jacobian): Assume \( \mathbf{f} \) is differentiable and its Jacobian satisfies \( \|\mathbf{J}(\mathbf{z})-\mathbf{J}(\mathbf{0})\|\le L_J\|\mathbf{z}\| \) for all \( \|\mathbf{z}\| < \rho \). Prove that \( \|\mathbf{r}(\mathbf{x})\|\le \dfrac{L_J}{2}\|\mathbf{x}\|^2 \) for all \( \|\mathbf{x}\| < \rho \), where \( \mathbf{r}(\mathbf{x})=\mathbf{f}(\mathbf{x})-\mathbf{J}(\mathbf{0})\mathbf{x} \).

Solution:

Use the integral identity \( \mathbf{f}(\mathbf{x})=\int_0^1 \mathbf{J}(s\mathbf{x})\mathbf{x}\,ds \) (since \( \mathbf{f}(\mathbf{0})=\mathbf{0} \)). Subtract \( \mathbf{J}(\mathbf{0})\mathbf{x}=\int_0^1 \mathbf{J}(\mathbf{0})\mathbf{x}\,ds \) to get

\[ \mathbf{r}(\mathbf{x}) = \int_0^1 \left(\mathbf{J}(s\mathbf{x})-\mathbf{J}(\mathbf{0})\right)\mathbf{x}\,ds. \]

Then \( \|\mathbf{r}(\mathbf{x})\|\le \int_0^1 L_J \|s\mathbf{x}\|\,\|\mathbf{x}\|\,ds = \dfrac{L_J}{2}\|\mathbf{x}\|^2 \).

Problem 3 (Certified inner ROA from a quadratic Lyapunov function): Consider \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{r}(\mathbf{x}) \), and assume \( \|\mathbf{r}(\mathbf{x})\|\le c\|\mathbf{x}\|^2 \) for \( \|\mathbf{x}\| < \rho \). Let \( \mathbf{A} \) be Hurwitz. Show that for any \( \mathbf{Q}\succ 0 \), if \( \mathbf{P}\succ 0 \) solves \( \mathbf{A}^T\mathbf{P}+\mathbf{P}\mathbf{A}=-\mathbf{Q} \), then there exists \( r>0 \) such that \( \dot{V}(\mathbf{x}) < 0 \) for all \( 0 < \|\mathbf{x}\| < r \), where \( V(\mathbf{x})=\mathbf{x}^T\mathbf{P}\mathbf{x} \).

Solution:

From Section 3, compute \( \dot{V} = -\mathbf{x}^T\mathbf{Q}\mathbf{x}+2\mathbf{x}^T\mathbf{P}\mathbf{r}(\mathbf{x}) \). Use bounds: \( \mathbf{x}^T\mathbf{Q}\mathbf{x}\ge \lambda_{\min}(\mathbf{Q})\|\mathbf{x}\|^2 \) and \( 2\mathbf{x}^T\mathbf{P}\mathbf{r}(\mathbf{x})\le 2c\|\mathbf{P}\|\|\mathbf{x}\|^3 \). Thus \( \dot{V}\le -\|\mathbf{x}\|^2(\lambda_{\min}(\mathbf{Q})-2c\|\mathbf{P}\|\|\mathbf{x}\|) \). Choose \( r=\min\!\left(\rho,\dfrac{\lambda_{\min}(\mathbf{Q})}{2c\|\mathbf{P}\|}\right) \) so that the parenthesis is positive for \( \|\mathbf{x}\| < r \), giving \( \dot{V} < 0 \).

Problem 4 (Inconclusive linearization): Analyze the equilibrium at 0 for \( \dot{x} = x^2 \). (a) Compute the linearization. (b) Determine stability of the true nonlinear equilibrium.

Solution:

(a) \( f(x)=x^2 \), so \( f'(0)=0 \) and linearization is \( \dot{x}=0 \), which does not imply attraction or repulsion.

(b) For any \( x_0 > 0 \), \( \dot{x}=x^2 > 0 \), so \( x(t) \) increases away from 0. For any \( x_0 < 0 \), \( \dot{x}=x^2 > 0 \), so the state increases toward 0 initially, crosses 0, and then grows positive. Hence the equilibrium is not asymptotically stable.

Problem 5 (ROA of a decoupled extension): For the 2D system used in Section 5, \( \dot{x}=-x+x^3,\;\dot{y}=-2y \), determine the region of attraction of \( (0,0) \).

Solution:

The system is decoupled. The \( y \)-subsystem satisfies \( y(t)=y_0 e^{-2t} \), so it converges for all \( y_0 \). The \( x \)-subsystem is identical to Problem 1, so convergence to 0 occurs iff \( |x_0| < 1 \). Therefore the region of attraction is

\[ \mathcal{R} = \{(x_0,y_0)\in\mathbb{R}^2:\; |x_0| < 1,\; y_0\in\mathbb{R}\}. \]

Linearization is a local approximation because the neglected nonlinear remainder is typically second-order: \( \|\mathbf{r}(\mathbf{x})\| = O(\|\mathbf{x}\|^2) \). The region of attraction \( \mathcal{R} \) formalizes which initial conditions converge to an equilibrium. Even when the Jacobian is stable, \( \mathcal{R} \) may be finite due to other equilibria or nonlinear effects. When the Jacobian is Hurwitz, a quadratic Lyapunov function derived from the linear model can certify a conservative inner subset \( \Omega_\alpha \subset \mathcal{R} \). Practical workflows combine (i) Jacobian eigenvalue checks, (ii) remainder bounding, (iii) Lyapunov-based certification, and (iv) simulation-based validation.