In Chapter 8 you modeled dynamics in state-space form. In Lesson 1 of this chapter you introduced operating points and perturbation variables. Here, we formalize how a nonlinear model becomes a local linear model via Jacobians.

Learning outcomes. By the end of this lesson, you should be able to:

• Write a nonlinear state-space model and define operating-point perturbations \( \delta \mathbf{x} \) and \( \delta \mathbf{u} \).
• Derive the linearized model matrices \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) using Jacobians.
• Prove (under smoothness assumptions) that the linearization error is second-order in perturbation magnitude.
• Compute Jacobians symbolically, by automatic differentiation, and by finite differences; understand numerical error mechanisms.

flowchart TD
  S["Start: nonlinear model xdot = f(x,u), y = g(x,u)"] --> OP["Choose operating point (x*, u*) or nominal trajectory (xbar(t), ubar(t))"]
  OP --> P["Define perturbations: dx = x - x*, du = u - u*"]
  P --> J["Compute Jacobians at operating point: A = df/dx, B = df/du, C = dg/dx, D = dg/du"]
  J --> L["Linearized model: dxdot = A dx + B du, dy = C dx + D du"]
  L --> E["Validate locally: compare nonlinear vs linear for small perturbations"]
        

Consider a general (continuous-time) nonlinear model with state \( \mathbf{x}(t)\in\mathbb{R}^n \) and input \( \mathbf{u}(t)\in\mathbb{R}^m \):

\[ \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t),\mathbf{u}(t)), \qquad \mathbf{y}(t) = \mathbf{g}(\mathbf{x}(t),\mathbf{u}(t)), \]

Let \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \) be an operating point. In Lesson 1, you defined perturbations: \( \delta\mathbf{x}(t) = \mathbf{x}(t)-\mathbf{x}^\ast \) and \( \delta\mathbf{u}(t) = \mathbf{u}(t)-\mathbf{u}^\ast \).

Two common cases:

• Equilibrium operating point: \( \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast)=\mathbf{0} \). Then \( \mathbf{x}(t)=\mathbf{x}^\ast \) is constant when \( \mathbf{u}(t)=\mathbf{u}^\ast \).
• Nominal trajectory (trajectory linearization): you specify a time-varying pair \( (\bar{\mathbf{x}}(t),\bar{\mathbf{u}}(t)) \) that satisfies the dynamics: \( \dot{\bar{\mathbf{x}}}(t)=\mathbf{f}(\bar{\mathbf{x}}(t),\bar{\mathbf{u}}(t)) \). Then you define perturbations \( \delta\mathbf{x}(t)=\mathbf{x}(t)-\bar{\mathbf{x}}(t) \), \( \delta\mathbf{u}(t)=\mathbf{u}(t)-\bar{\mathbf{u}}(t) \).

In this lesson, we derive the linearized model for both cases (equilibrium is the most common and simplest special case).

Assume \( \mathbf{f} \) is continuously differentiable in a neighborhood of the operating point. Define the Jacobians of \( \mathbf{f} \) with respect to state and input:

\[ \mathbf{A}(\mathbf{x},\mathbf{u}) \equiv \frac{\partial \mathbf{f}}{\partial \mathbf{x}}(\mathbf{x},\mathbf{u}) \in \mathbb{R}^{n\times n}, \qquad \mathbf{B}(\mathbf{x},\mathbf{u}) \equiv \frac{\partial \mathbf{f}}{\partial \mathbf{u}}(\mathbf{x},\mathbf{u}) \in \mathbb{R}^{n\times m}. \]

Component-wise, for \( i=1,\dots,n \) and \( j=1,\dots,n \): \( [\mathbf{A}]_{ij} = \frac{\partial f_i}{\partial x_j} \), and for \( k=1,\dots,m \): \( [\mathbf{B}]_{ik} = \frac{\partial f_i}{\partial u_k} \).

Similarly for the output mapping \( \mathbf{g} \), define: \( \mathbf{C} = \frac{\partial \mathbf{g}}{\partial \mathbf{x}} \in \mathbb{R}^{p\times n} \) and \( \mathbf{D} = \frac{\partial \mathbf{g}}{\partial \mathbf{u}} \in \mathbb{R}^{p\times m} \).

3.1 First-order Taylor approximation (vector form)

Let \( \mathbf{z} = \begin{bmatrix}\mathbf{x}\\ \mathbf{u}\end{bmatrix}\in\mathbb{R}^{n+m} \), and define \( \mathbf{F}(\mathbf{z}) = \mathbf{f}(\mathbf{x},\mathbf{u}) \). For small perturbation \( \delta\mathbf{z} = \begin{bmatrix}\delta\mathbf{x}\\ \delta\mathbf{u}\end{bmatrix} \), the first-order Taylor expansion gives:

\[ \mathbf{f}(\mathbf{x}^\ast+\delta\mathbf{x},\mathbf{u}^\ast+\delta\mathbf{u}) = \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) + \mathbf{A}^\ast\,\delta\mathbf{x} + \mathbf{B}^\ast\,\delta\mathbf{u} + \mathbf{r}_f(\delta\mathbf{x},\delta\mathbf{u}), \]

where \( \mathbf{A}^\ast = \mathbf{A}(\mathbf{x}^\ast,\mathbf{u}^\ast) \) and \( \mathbf{B}^\ast = \mathbf{B}(\mathbf{x}^\ast,\mathbf{u}^\ast) \). The remainder term \( \mathbf{r}_f \) is “higher order” (made precise next).

3.2 Remainder term is second-order (a precise bound)

If \( \mathbf{f} \) is twice continuously differentiable in a neighborhood of \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \), then the remainder satisfies a quadratic bound. One convenient statement uses the integral form of Taylor’s theorem.

Theorem (Second-order remainder bound). Suppose \( \mathbf{F}:\mathbb{R}^{n+m}\to\mathbb{R}^n \) is \( C^2 \) on a convex neighborhood containing the line segment \( \{\mathbf{z}^\ast + s\,\delta\mathbf{z}\mid 0\le s\le 1\} \). Then there exists a constant \( M \) such that for sufficiently small \( \delta\mathbf{z} \):

\[ \left\| \mathbf{r}_f(\delta\mathbf{x},\delta\mathbf{u}) \right\| \le \frac{M}{2}\,\left\|\delta\mathbf{z}\right\|^2. \]

Proof (sketch with full mathematical steps). Define \( \phi(s)=\mathbf{F}(\mathbf{z}^\ast+s\delta\mathbf{z}) \). By the chain rule:

\[ \phi'(s) = \mathbf{J}_F(\mathbf{z}^\ast+s\delta\mathbf{z})\,\delta\mathbf{z}, \quad \phi''(s) = \left[\nabla^2 \mathbf{F}(\mathbf{z}^\ast+s\delta\mathbf{z})\right]\,[\delta\mathbf{z},\delta\mathbf{z}], \]

where \( \mathbf{J}_F \) is the Jacobian of \( \mathbf{F} \) (an \( n\times(n+m) \) matrix), and \( \nabla^2\mathbf{F} \) denotes the vector of Hessians (a bilinear form applied to \( \delta\mathbf{z} \) twice). Apply the scalar Taylor theorem to each component of \( \phi \) around \( s=0 \):

\[ \phi(1) = \phi(0) + \phi'(0)\cdot 1 + \int_0^1 (1-s)\,\phi''(s)\,ds. \]

Noting \( \phi(1)=\mathbf{F}(\mathbf{z}^\ast+\delta\mathbf{z}) \), \( \phi(0)=\mathbf{F}(\mathbf{z}^\ast) \), and \( \phi'(0)=\mathbf{J}_F(\mathbf{z}^\ast)\delta\mathbf{z} \), we identify the remainder as:

\[ \mathbf{r}_f(\delta\mathbf{x},\delta\mathbf{u}) = \int_0^1 (1-s)\,\left[\nabla^2 \mathbf{F}(\mathbf{z}^\ast+s\delta\mathbf{z})\right]\,[\delta\mathbf{z},\delta\mathbf{z}]\,ds. \]

If the Hessians are bounded on the segment, i.e., \( \left\|\left[\nabla^2 \mathbf{F}(\mathbf{z})\right]\right\|\le M \), then using \( \int_0^1 (1-s)\,ds = \tfrac{1}{2} \) yields:

\[ \left\|\mathbf{r}_f(\delta\mathbf{x},\delta\mathbf{u})\right\| \le \int_0^1 (1-s)\,M\,\|\delta\mathbf{z}\|^2\,ds = \frac{M}{2}\,\|\delta\mathbf{z}\|^2. \]

Thus the linearization error is quadratic in perturbation magnitude. ■

4.1 Equilibrium operating point

Assume \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \) is an equilibrium: \( \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast)=\mathbf{0} \). Substitute \( \mathbf{x}=\mathbf{x}^\ast+\delta\mathbf{x} \), \( \mathbf{u}=\mathbf{u}^\ast+\delta\mathbf{u} \):

\[ \dot{\mathbf{x}} = \dot{\delta\mathbf{x}} = \mathbf{f}(\mathbf{x}^\ast+\delta\mathbf{x},\mathbf{u}^\ast+\delta\mathbf{u}) = \mathbf{A}^\ast\delta\mathbf{x} + \mathbf{B}^\ast\delta\mathbf{u} + \mathbf{r}_f(\delta\mathbf{x},\delta\mathbf{u}). \]

Neglecting the higher-order remainder yields the linearized (small-signal) model:

\[ \dot{\delta\mathbf{x}} = \mathbf{A}^\ast\,\delta\mathbf{x} + \mathbf{B}^\ast\,\delta\mathbf{u}. \]

For outputs:

\[ \delta\mathbf{y} = \mathbf{y}-\mathbf{y}^\ast = \mathbf{g}(\mathbf{x}^\ast+\delta\mathbf{x},\mathbf{u}^\ast+\delta\mathbf{u}) - \mathbf{g}(\mathbf{x}^\ast,\mathbf{u}^\ast) \approx \mathbf{C}^\ast\delta\mathbf{x} + \mathbf{D}^\ast\delta\mathbf{u}. \]

4.2 Nominal trajectory (time-varying linearization)

If \( \bar{\mathbf{x}}(t) \) and \( \bar{\mathbf{u}}(t) \) satisfy \( \dot{\bar{\mathbf{x}}}(t)=\mathbf{f}(\bar{\mathbf{x}}(t),\bar{\mathbf{u}}(t)) \), then for \( \delta\mathbf{x}(t)=\mathbf{x}(t)-\bar{\mathbf{x}}(t) \), \( \delta\mathbf{u}(t)=\mathbf{u}(t)-\bar{\mathbf{u}}(t) \):

\[ \dot{\delta\mathbf{x}}(t) = \mathbf{A}(t)\,\delta\mathbf{x}(t) + \mathbf{B}(t)\,\delta\mathbf{u}(t) + \mathbf{r}_f(t), \quad \mathbf{A}(t)=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right|_{(\bar{\mathbf{x}}(t),\bar{\mathbf{u}}(t))}, \quad \mathbf{B}(t)=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\right|_{(\bar{\mathbf{x}}(t),\bar{\mathbf{u}}(t))}. \]

Ignoring \( \mathbf{r}_f(t) \) yields a linear time-varying perturbation model. (Lesson 3 will formalize when this approximation is valid over time intervals.)

4.3 Interpretation: what the Jacobian “means” dynamically

The Jacobian \( \mathbf{A}^\ast \) is the best linear map (first-order) that predicts how the vector field changes with state: for “small” perturbations, \( \mathbf{f}(\mathbf{x}^\ast+\delta\mathbf{x},\mathbf{u}^\ast) \approx \mathbf{A}^\ast\delta\mathbf{x} \). Similarly \( \mathbf{B}^\ast \) predicts how the vector field changes with input.

There are three mainstream computational routes to obtain \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \):

1. Symbolic differentiation: exact expressions; can become large for complex models.
2. Automatic differentiation (AD): exact to machine precision by differentiating the program evaluation graph.
3. Finite differences: easy and model-agnostic, but introduces truncation and round-off error.

flowchart TD
  Q["Need Jacobians A,B,C,D"] --> S1["Do you have symbolic model?"]
  S1 -->|yes| SYM["Use symbolic differentiation \n(exact expressions)"]
  S1 -->|no| S2["Is AD available in your toolchain?"]
  S2 -->|yes| AD["Use automatic differentiation \n(machine-accurate derivatives)"]
  S2 -->|no| FD["Use finite differences \n(choose step size carefully)"]
  SYM --> EV["Evaluate at (x*,u*) or (xbar(t), ubar(t))"]
  AD --> EV
  FD --> EV
        

5.1 Finite-difference Jacobian and error order

For a function \( \mathbf{f}:\mathbb{R}^n\to\mathbb{R}^n \) (fixing \( \mathbf{u} \)), a forward-difference approximation for column \( j \) is:

\[ \mathbf{A}_{:,j} \approx \frac{\mathbf{f}(\mathbf{x}^\ast + h\,\mathbf{e}_j,\mathbf{u}^\ast)-\mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast)}{h}, \qquad h > 0, \]

where \( \mathbf{e}_j \) is the \( j \)-th standard basis vector. If \( \mathbf{f} \) is smooth, the truncation error is \( O(h) \). A centered difference improves to \( O(h^2) \):

\[ \mathbf{A}_{:,j} \approx \frac{\mathbf{f}(\mathbf{x}^\ast + h\,\mathbf{e}_j,\mathbf{u}^\ast)-\mathbf{f}(\mathbf{x}^\ast - h\,\mathbf{e}_j,\mathbf{u}^\ast)}{2h}. \]

However, too small \( h \) amplifies floating-point round-off. In practice, \( h \) is chosen by balancing truncation and round-off (Problem 4 formalizes this idea).

Consider a two-state system with a smooth nonlinearity (a sine term is a standard prototype for many mechanical models): state \( \mathbf{x} = \begin{bmatrix}x_1 \\ x_2\end{bmatrix} \) and scalar input \( u \).

\[ \dot{x}_1 = x_2,\qquad \dot{x}_2 = -a\sin(x_1) - b x_2 + c\,u, \qquad y = x_1, \]

Define \( \mathbf{f}(\mathbf{x},u)=\begin{bmatrix} f_1 \\ f_2 \end{bmatrix} \) with \( f_1=x_2 \), \( f_2=-a\sin(x_1)-bx_2+cu \). The Jacobians are:

\[ \mathbf{A}(\mathbf{x},u)=\frac{\partial\mathbf{f}}{\partial\mathbf{x}}= \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2}\\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} \end{bmatrix} = \begin{bmatrix} 0 & 1\\ -a\cos(x_1) & -b \end{bmatrix}, \qquad \mathbf{B}(\mathbf{x},u)=\frac{\partial\mathbf{f}}{\partial u}= \begin{bmatrix} 0\\ c \end{bmatrix}. \]

Linearize about the equilibrium \( x_1^\ast=0,\; x_2^\ast=0,\; u^\ast=0 \). Since \( \cos(0)=1 \):

\[ \mathbf{A}^\ast = \begin{bmatrix} 0 & 1\\ -a & -b \end{bmatrix}, \qquad \mathbf{B}^\ast = \begin{bmatrix} 0\\ c \end{bmatrix}. \]

For the output \( y=x_1 \), we have \( \mathbf{g}(\mathbf{x},u)=x_1 \) so:

\[ \mathbf{C}^\ast = \frac{\partial y}{\partial \mathbf{x}} = \begin{bmatrix} 1 & 0 \end{bmatrix}, \qquad \mathbf{D}^\ast = \frac{\partial y}{\partial u} = \begin{bmatrix} 0 \end{bmatrix}. \]

Therefore the small-signal model is:

\[ \dot{\delta\mathbf{x}} = \begin{bmatrix} 0 & 1\\ -a & -b \end{bmatrix}\delta\mathbf{x} + \begin{bmatrix} 0\\ c \end{bmatrix}\delta u, \qquad \delta y = \begin{bmatrix} 1 & 0 \end{bmatrix}\delta\mathbf{x}. \]

The goal in each implementation is consistent: (i) define \( \mathbf{f}(\mathbf{x},\mathbf{u}) \) and \( \mathbf{g}(\mathbf{x},\mathbf{u}) \), (ii) compute Jacobians \( \mathbf{A}^\ast,\mathbf{B}^\ast,\mathbf{C}^\ast,\mathbf{D}^\ast \) at an operating point, and (iii) use the linearized model for small-signal analysis or as a local surrogate in simulation.

7.1 Python (Symbolic Jacobian with SymPy + optional AD with JAX)

Recommended libraries in Python for system dynamics workflows: \( \texttt{sympy} \) (symbolic Jacobians), \( \texttt{numpy} \) (numerics), \( \texttt{scipy} \) (ODE solvers in later chapters), and \( \texttt{jax} \) (automatic differentiation).

import sympy as sp

# Symbols
x1, x2, u = sp.symbols('x1 x2 u', real=True)
a, b, c = sp.symbols('a b c', real=True)

# Nonlinear dynamics f(x,u) and output g(x,u)
f = sp.Matrix([x2, -a*sp.sin(x1) - b*x2 + c*u])
g = sp.Matrix([x1])

x = sp.Matrix([x1, x2])

# Jacobians
A = f.jacobian(x)     # df/dx
B = f.jacobian(sp.Matrix([u]))  # df/du
C = g.jacobian(x)     # dg/dx
D = g.jacobian(sp.Matrix([u]))  # dg/du

# Evaluate at operating point x*=[0,0], u*=0
subs_op = {x1: 0, x2: 0, u: 0}
A_star = sp.simplify(A.subs(subs_op))
B_star = sp.simplify(B.subs(subs_op))
C_star = sp.simplify(C.subs(subs_op))
D_star = sp.simplify(D.subs(subs_op))

print("A* =", A_star)
print("B* =", B_star)
print("C* =", C_star)
print("D* =", D_star)

# Optional: numerical verification via centered finite differences (model-agnostic)
import numpy as np

def f_num(x_vec, u_val, a_val=2.0, b_val=0.4, c_val=1.5):
    x1v, x2v = x_vec
    return np.array([x2v, -a_val*np.sin(x1v) - b_val*x2v + c_val*u_val], dtype=float)

def jacobian_fd_centered(f_handle, x_star, u_star, h=1e-6):
    n = x_star.size
    Ahat = np.zeros((n, n))
    for j in range(n):
        ej = np.zeros(n); ej[j] = 1.0
        fp = f_handle(x_star + h*ej, u_star)
        fm = f_handle(x_star - h*ej, u_star)
        Ahat[:, j] = (fp - fm) / (2.0*h)
    # For B with scalar u:
    fp_u = f_handle(x_star, u_star + h)
    fm_u = f_handle(x_star, u_star - h)
    Bhat = ((fp_u - fm_u) / (2.0*h)).reshape(n, 1)
    return Ahat, Bhat

A_fd, B_fd = jacobian_fd_centered(f_num, np.array([0.0, 0.0]), 0.0)
print("A_fd =", A_fd)
print("B_fd =", B_fd)
      

If you use JAX, you can obtain Jacobians by differentiating the code of \( \mathbf{f} \) directly (automatic differentiation), which avoids step-size tuning inherent in finite differences.

7.2 C++ (Finite Differences with Eigen; AD note with CppAD/autodiff)

Recommended libraries in C++: \( \texttt{Eigen} \) (matrix computations), and for automatic differentiation, common choices include \( \texttt{CppAD} \) or modern header-only AD libraries (toolchain-dependent).

#include <iostream>
#include <cmath>
#include <Eigen/Dense>

using Eigen::VectorXd;
using Eigen::MatrixXd;

// Example nonlinear f(x,u) for x in R^2, u in R^1
VectorXd f_num(const VectorXd& x, double u, double a, double b, double c) {
    VectorXd fx(2);
    const double x1 = x(0);
    const double x2 = x(1);
    fx(0) = x2;
    fx(1) = -a*std::sin(x1) - b*x2 + c*u;
    return fx;
}

// Centered finite-difference Jacobians A=df/dx and B=df/du at (x*,u*)
void jacobian_centered(
    const VectorXd& x_star, double u_star,
    double a, double b, double c,
    double h, MatrixXd& A, MatrixXd& B
) {
    const int n = (int)x_star.size();
    A = MatrixXd::Zero(n, n);
    B = MatrixXd::Zero(n, 1);

    // A columns
    for (int j = 0; j < n; ++j) {
        VectorXd ej = VectorXd::Zero(n);
        ej(j) = 1.0;

        VectorXd fp = f_num(x_star + h*ej, u_star, a, b, c);
        VectorXd fm = f_num(x_star - h*ej, u_star, a, b, c);
        A.col(j) = (fp - fm) / (2.0*h);
    }

    // B for scalar u
    VectorXd fp_u = f_num(x_star, u_star + h, a, b, c);
    VectorXd fm_u = f_num(x_star, u_star - h, a, b, c);
    B.col(0) = (fp_u - fm_u) / (2.0*h);
}

int main() {
    VectorXd x_star(2);
    x_star << 0.0, 0.0;
    double u_star = 0.0;

    double a = 2.0, b = 0.4, c = 1.5;
    double h = 1e-6;

    MatrixXd A, B;
    jacobian_centered(x_star, u_star, a, b, c, h, A, B);

    std::cout << "A* (FD) =\n" << A << "\n";
    std::cout << "B* (FD) =\n" << B << "\n";

    return 0;
}
      

AD alternative (conceptual): if your model is implemented as code and you require high accuracy without step-size selection, use an AD library to compute \( \mathbf{A}^\ast,\mathbf{B}^\ast \) by differentiating the evaluation trace of \( \mathbf{f} \). This is especially useful for large-scale models where symbolic expressions are unwieldy.

7.3 Java (Finite Differences with EJML matrices)

Recommended libraries in Java: \( \texttt{EJML} \) (Efficient Java Matrix Library) for linear algebra, and \( \texttt{Apache Commons Math} \) for numerical utilities. AD is less standardized in Java; finite differences are common.

import org.ejml.simple.SimpleMatrix;

public class JacobianLinearization {

    // f(x,u) for x in R^2, u scalar
    static SimpleMatrix fNum(SimpleMatrix x, double u, double a, double b, double c) {
        double x1 = x.get(0,0);
        double x2 = x.get(1,0);
        SimpleMatrix fx = new SimpleMatrix(2,1);
        fx.set(0,0, x2);
        fx.set(1,0, -a*Math.sin(x1) - b*x2 + c*u);
        return fx;
    }

    static class Jacobians {
        SimpleMatrix A;
        SimpleMatrix B;
        Jacobians(SimpleMatrix A, SimpleMatrix B) { this.A = A; this.B = B; }
    }

    // Centered finite differences
    static Jacobians jacobianCentered(SimpleMatrix xStar, double uStar,
                                      double a, double b, double c, double h) {
        int n = xStar.numRows();
        SimpleMatrix A = new SimpleMatrix(n, n);
        SimpleMatrix B = new SimpleMatrix(n, 1);

        // A columns
        for (int j = 0; j < n; j++) {
            SimpleMatrix ej = new SimpleMatrix(n,1);
            ej.set(j,0, 1.0);

            SimpleMatrix fp = fNum(xStar.plus(ej.scale(h)), uStar, a, b, c);
            SimpleMatrix fm = fNum(xStar.minus(ej.scale(h)), uStar, a, b, c);
            SimpleMatrix col = fp.minus(fm).divide(2.0*h);
            for (int i = 0; i < n; i++) A.set(i,j, col.get(i,0));
        }

        // B for scalar u
        SimpleMatrix fpU = fNum(xStar, uStar + h, a, b, c);
        SimpleMatrix fmU = fNum(xStar, uStar - h, a, b, c);
        SimpleMatrix bhat = fpU.minus(fmU).divide(2.0*h);
        B.set(0,0, bhat.get(0,0));
        B.set(1,0, bhat.get(1,0));

        return new Jacobians(A, B);
    }

    public static void main(String[] args) {
        SimpleMatrix xStar = new SimpleMatrix(2,1);
        xStar.set(0,0, 0.0);
        xStar.set(1,0, 0.0);
        double uStar = 0.0;

        double a = 2.0, b = 0.4, c = 1.5;
        double h = 1e-6;

        Jacobians J = jacobianCentered(xStar, uStar, a, b, c, h);
        System.out.println("A* (FD) =\n" + J.A);
        System.out.println("B* (FD) =\n" + J.B);
    }
}
      

7.4 MATLAB (Symbolic Jacobian) and Simulink (Linear Analysis)

In MATLAB, symbolic Jacobians use the Symbolic Math Toolbox. For local linear models in state-space form, you can evaluate \( \mathbf{A}^\ast,\mathbf{B}^\ast,\mathbf{C}^\ast,\mathbf{D}^\ast \) at an operating point and then form \( \texttt{ss(A,B,C,D)} \) (Control System Toolbox).

% Symbolic Jacobian linearization for the worked example
syms x1 x2 u a b c real

f = [ x2;
     -a*sin(x1) - b*x2 + c*u ];

g = x1;

x = [x1; x2];

A = jacobian(f, x);      % df/dx
B = jacobian(f, u);      % df/du
C = jacobian(g, x);      % dg/dx
D = jacobian(g, u);      % dg/du

% Operating point
x_star = [0; 0];
u_star = 0;

A_star = subs(A, [x1 x2 u], [x_star.' u_star]);
B_star = subs(B, [x1 x2 u], [x_star.' u_star]);
C_star = subs(C, [x1 x2 u], [x_star.' u_star]);
D_star = subs(D, [x1 x2 u], [x_star.' u_star]);

disp('A* ='); disp(simplify(A_star));
disp('B* ='); disp(simplify(B_star));
disp('C* ='); disp(simplify(C_star));
disp('D* ='); disp(simplify(D_star));

% Optional: form state-space model (requires numeric parameter values)
aN = 2.0; bN = 0.4; cN = 1.5;
A_num = double(subs(A_star, [a b c], [aN bN cN]));
B_num = double(subs(B_star, [a b c], [aN bN cN]));
C_num = double(C_star);
D_num = double(D_star);

sys_lin = ss(A_num, B_num, C_num, D_num);
      

In Simulink, the typical workflow is: build the nonlinear model, define an operating point (equilibrium or snapshot along a trajectory), and then use linearization tools (e.g., \( \texttt{linearize} \) or the Linear Analysis Tool) to obtain \( \mathbf{A}^\ast,\mathbf{B}^\ast,\mathbf{C}^\ast,\mathbf{D}^\ast \). This produces the same Jacobian linearization, but computed from the model structure.

7.5 Wolfram Mathematica (Exact Jacobians with D[])

(* Symbolic Jacobian linearization for the worked example *)
ClearAll[x1, x2, u, a, b, c];

f = {x2, -a Sin[x1] - b x2 + c u};
g = {x1};
x = {x1, x2};

A = D[f, {x}];  (* df/dx *)
B = D[f, u];    (* df/du *)
C = D[g, {x}];  (* dg/dx *)
Dmat = D[g, u]; (* dg/du *)

op = {x1 -> 0, x2 -> 0, u -> 0};

Astar = Simplify[A /. op]
Bstar = Simplify[B /. op]
Cstar = Simplify[C /. op]
Dstar = Simplify[Dmat /. op]
      

The problems below focus on deriving Jacobian linearizations and quantifying approximation and numerical differentiation error. Work symbolically where possible; evaluate at operating points only after obtaining general expressions.

Problem 1 (Derive the small-signal model at an equilibrium): Let \( \dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u}) \), with an equilibrium \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \) satisfying \( \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast)=\mathbf{0} \). Using the perturbations \( \delta\mathbf{x}=\mathbf{x}-\mathbf{x}^\ast \), \( \delta\mathbf{u}=\mathbf{u}-\mathbf{u}^\ast \), show that the linearized model is \( \dot{\delta\mathbf{x}}=\mathbf{A}^\ast\delta\mathbf{x}+\mathbf{B}^\ast\delta\mathbf{u} \), and identify \( \mathbf{A}^\ast,\mathbf{B}^\ast \).

Solution: Expand \( \mathbf{f} \) to first order around \( (\mathbf{x}^\ast,\mathbf{u}^\ast) \):

\[ \mathbf{f}(\mathbf{x}^\ast+\delta\mathbf{x},\mathbf{u}^\ast+\delta\mathbf{u}) = \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast) + \left.\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right|_{(\mathbf{x}^\ast,\mathbf{u}^\ast)}\delta\mathbf{x} + \left.\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\right|_{(\mathbf{x}^\ast,\mathbf{u}^\ast)}\delta\mathbf{u} + \mathbf{r}_f. \]

Since \( \mathbf{f}(\mathbf{x}^\ast,\mathbf{u}^\ast)=\mathbf{0} \) and \( \dot{\mathbf{x}}=\dot{\delta\mathbf{x}} \), dropping the higher-order remainder gives: \( \dot{\delta\mathbf{x}}=\mathbf{A}^\ast\delta\mathbf{x}+\mathbf{B}^\ast\delta\mathbf{u} \), with \( \mathbf{A}^\ast=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right|_{(\mathbf{x}^\ast,\mathbf{u}^\ast)} \), \( \mathbf{B}^\ast=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\right|_{(\mathbf{x}^\ast,\mathbf{u}^\ast)} \).

Problem 2 (Compute Jacobians for a 2-state nonlinear system): For \( \dot{x}_1 = x_2 \), \( \dot{x}_2 = -a\sin(x_1) - b x_2 + c u \), \( y=x_1 \), compute \( \mathbf{A}(\mathbf{x},u),\mathbf{B}(\mathbf{x},u),\mathbf{C}(\mathbf{x},u),\mathbf{D}(\mathbf{x},u) \) and then evaluate at \( (x_1^\ast,x_2^\ast,u^\ast)=(0,0,0) \).

Solution:

\[ \mathbf{A}(\mathbf{x},u)= \begin{bmatrix} 0 & 1\\ -a\cos(x_1) & -b \end{bmatrix}, \quad \mathbf{B}(\mathbf{x},u)= \begin{bmatrix} 0\\ c \end{bmatrix}, \quad \mathbf{C}(\mathbf{x},u)=\begin{bmatrix}1 & 0\end{bmatrix}, \quad \mathbf{D}(\mathbf{x},u)=\begin{bmatrix}0\end{bmatrix}. \]

At \( x_1^\ast=0 \), \( \cos(0)=1 \), hence:

\[ \mathbf{A}^\ast= \begin{bmatrix} 0 & 1\\ -a & -b \end{bmatrix}, \quad \mathbf{B}^\ast= \begin{bmatrix} 0\\ c \end{bmatrix}, \quad \mathbf{C}^\ast=\begin{bmatrix}1 & 0\end{bmatrix}, \quad \mathbf{D}^\ast=\begin{bmatrix}0\end{bmatrix}. \]

Problem 3 (Trajectory linearization / variational equation): Let \( \bar{\mathbf{x}}(t) \), \( \bar{\mathbf{u}}(t) \) satisfy the nonlinear dynamics: \( \dot{\bar{\mathbf{x}}}(t)=\mathbf{f}(\bar{\mathbf{x}}(t),\bar{\mathbf{u}}(t)) \). Define perturbations \( \delta\mathbf{x}(t)=\mathbf{x}(t)-\bar{\mathbf{x}}(t) \), \( \delta\mathbf{u}(t)=\mathbf{u}(t)-\bar{\mathbf{u}}(t) \). Derive the first-order perturbation dynamics.

Solution: Substitute \( \mathbf{x}=\bar{\mathbf{x}}+\delta\mathbf{x} \), \( \mathbf{u}=\bar{\mathbf{u}}+\delta\mathbf{u} \):

\[ \dot{\bar{\mathbf{x}}}+\dot{\delta\mathbf{x}} = \mathbf{f}(\bar{\mathbf{x}}+\delta\mathbf{x},\bar{\mathbf{u}}+\delta\mathbf{u}) = \mathbf{f}(\bar{\mathbf{x}},\bar{\mathbf{u}}) + \mathbf{A}(t)\delta\mathbf{x} + \mathbf{B}(t)\delta\mathbf{u} + \mathbf{r}_f(t), \]

where \( \mathbf{A}(t)=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{x}}\right|_{(\bar{\mathbf{x}}(t),\bar{\mathbf{u}}(t))} \) and \( \mathbf{B}(t)=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\right|_{(\bar{\mathbf{x}}(t),\bar{\mathbf{u}}(t))} \). Using \( \dot{\bar{\mathbf{x}}}=\mathbf{f}(\bar{\mathbf{x}},\bar{\mathbf{u}}) \) cancels the nominal terms, giving:

\[ \dot{\delta\mathbf{x}}(t) = \mathbf{A}(t)\delta\mathbf{x}(t) + \mathbf{B}(t)\delta\mathbf{u}(t) + \mathbf{r}_f(t). \]

Dropping \( \mathbf{r}_f(t) \) yields the time-varying linear perturbation model \( \dot{\delta\mathbf{x}}(t) \approx \mathbf{A}(t)\delta\mathbf{x}(t) + \mathbf{B}(t)\delta\mathbf{u}(t) \).

Problem 4 (Finite-difference step size trade-off): For a scalar smooth function \( \varphi(x) \), the centered difference \( \frac{\varphi(x+h)-\varphi(x-h)}{2h} \) has truncation error \( O(h^2) \). Explain why, in floating-point arithmetic, choosing \( h \) “too small” can worsen accuracy, and give a qualitative rule for choosing \( h \).

Solution: Centered differences reduce truncation error, but when \( h \) is very small, the subtraction \( \varphi(x+h)-\varphi(x-h) \) can lose significant digits (catastrophic cancellation), making round-off dominate. A qualitative balance is: pick \( h \) small enough that the local curvature is well-approximated, but not so small that subtraction errors dominate. In practice, a common starting point is \( h \approx 10^{-6} \) to \( 10^{-8} \) times a characteristic scale of \( x \), then refine by checking sensitivity of the estimated Jacobian to moderate changes in \( h \).

Jacobian linearization replaces a nonlinear state-space model by a local linear approximation around an operating point (or nominal trajectory). The linearized dynamics are determined by the Jacobians \( \mathbf{A}^\ast=\frac{\partial\mathbf{f}}{\partial\mathbf{x}} \), \( \mathbf{B}^\ast=\frac{\partial\mathbf{f}}{\partial\mathbf{u}} \), and similarly \( \mathbf{C}^\ast,\mathbf{D}^\ast \) for outputs. Under smoothness assumptions, the neglected terms are second-order in perturbation size, explaining why the approximation is valid only locally. Computationally, Jacobians can be obtained symbolically, via automatic differentiation, or by finite differences with careful step-size selection.