Linearization and the Jacobian

Near an equilibrium, how does a nonlinear system actually behave? This question sits at the heart of nearly all applied dynamics. Almost no system encountered in physics, biology, or engineering is exactly linear. The strategy of linearization answers it by replacing the full nonlinear vector field with its best linear approximation at a fixed point. This turns an intractable nonlinear problem into a linear algebra problem we already know how to solve. This matters because the eigenvalues of a single matrix — the Jacobian — often determine, with mathematical certainty, whether trajectories starting nearby will spiral away, decay back, or do something the linear picture cannot capture at all. The big picture is a bridge between two worlds: the qualitative theory of Fixed Points and Equilibria within Autonomous Systems and Phase Space, and the fully computable world of constant-coefficient linear systems. Linearization builds that bridge explicitly, via a first-order Taylor expansion of the vector field. The Jacobian matrix is the coefficient matrix of the resulting linear system. This article develops the construction with full rigor. It defines the Jacobian, and proves that the linearized flow is the correct first-order approximation to the nonlinear flow. It also states precisely the theorem — due in spirit to Lyapunov and Perron — that governs when this approximation is qualitatively trustworthy, not merely numerically close. We will be careful about the difference between two claims that are often blurred in a first course. First: the linear vector field approximates the nonlinear one to first order near the fixed point (an easy calculus fact). Second: the phase portraits near the fixed point are topologically the same (a much deeper fact, requiring the hypothesis of hyperbolicity). The article proves the first claim completely. It states the precise hypotheses under which the second becomes available, foreshadowing the Hartman-Grobman and Hyperbolicity theorem and the subsequent Classification of Planar Fixed Points. By the end you will be able to: compute the Jacobian of any smooth autonomous system, linearize at a fixed point, and interpret eigenvalues correctly. Crucially, you will also recognize the exact situation in which linearization fails to determine stability, with a worked counterexample showing why.

Key takeaways
  • The Jacobian Df(x∗)Df(x^*) is the matrix of partial derivatives of a vector field ff evaluated at a fixed point x∗x^*.
  • Linearization replaces x˙=f(x)\dot x = f(x) near x∗x^* by the linear system u˙=Df(x∗) u\dot u = Df(x^*)\,u, u=x−x∗u = x-x^*.
  • This is rigorously justified as the unique best linear (first-order Taylor) approximation to ff at x∗x^*.
  • A fixed point is hyperbolic if no eigenvalue of Df(x∗)Df(x^*) has zero real part; only then is the local phase portrait guaranteed to match the linear one (Hartman–Grobman, foreshadowed here).
  • Eigenvalues with nonzero real part correctly predict local stability/instability; purely imaginary or zero eigenvalues do not — the linearization is inconclusive.
  • Key idea: the Jacobian converts a nonlinear stability question into an eigenvalue problem, valid to first order always, and valid qualitatively exactly when the fixed point is hyperbolic.

Where this lives today

In physics, linearizing the equations of a pendulum, oscillator network, or fluid flow about an equilibrium determines normal modes and whether small perturbations grow or decay, directly feeding into stability analyses in mechanics and plasma physics. In control engineering, the Jacobian at an operating point defines the state matrix used for pole placement, LQR design, and Kalman filtering of nonlinear plants — nearly all classical linear control theory begins by linearizing a nonlinear model this way. In ecology and epidemiology, the Jacobian at an equilibrium of a predator–prey or SIR-type model (evaluated at the disease-free or coexistence equilibrium) determines the basic reproduction number's sign structure and whether the equilibrium resists small outbreaks. In numerical analysis and optimization, Newton's method for root-finding is literally an iterated application of the same first-order Taylor approximation used in part (a), and the Jacobian appears identically in Newton–Raphson solvers for systems of nonlinear equations. In computer science and machine learning, Jacobians of layer maps appear in backpropagation and in analyzing the stability of recurrent neural network dynamics near fixed points. In cryptography-adjacent dynamical systems research, Lyapunov exponents (which generalize the eigenvalue real-part criterion along entire trajectories, not just at fixed points) are used to certify chaotic behavior in pseudo-random generators built from nonlinear maps.

Ask the tutor: How do the eigenvalues and eigenvectors of the Jacobian determine the specific shape (node, spiral, or saddle) of the phase portrait near a planar fixed point?

DefinitionJacobian matrix and linearization

Let f:U→Rnf:U \to \mathbb{R}^n be a C1C^1 vector field on an open set U⊆RnU \subseteq \mathbb{R}^n, generating the autonomous system

x˙=f(x),x∈U.\dot x = f(x), \qquad x \in U .

Suppose x∗∈Ux^* \in U is a fixed point, i.e. f(x∗)=0f(x^*) = 0. The Jacobian matrix of ff at x∗x^* is the n×nn \times n matrix of first partial derivatives

Df(x∗)  =  [∂fi∂xj(x∗)]i,j=1n,Df(x^*) \;=\; \left[ \frac{\partial f_i}{\partial x_j}(x^*) \right]_{i,j=1}^n ,

where f=(f1,…,fn)f = (f_1,\dots,f_n) and x=(x1,…,xn)x = (x_1,\dots,x_n). The linearization of the system at x∗x^* is the constant-coefficient linear system

u˙=Au,A:=Df(x∗),u:=x−x∗,\dot u = A u, \qquad A := Df(x^*), \qquad u := x - x^* ,

obtained by dropping all terms of order ∥u∥2\lVert u \rVert^2 and higher in the Taylor expansion of ff about x∗x^*. The fixed point x∗x^* is called hyperbolic if every eigenvalue of A=Df(x∗)A = Df(x^*) has nonzero real part.

Intuition

Picture the vector field ff as an arrow attached to every point of the plane (or of Rn\mathbb{R}^n), telling a particle which way to move and how fast. At a fixed point x∗x^* the arrow vanishes: a particle placed exactly there stays forever. The interesting question is what happens to a particle placed near x∗x^*, and this is precisely where linearization earns its keep: sufficiently close to x∗x^*, the vector field looks almost exactly like a linear vector field, because that is what "differentiable" means. The Jacobian A=Df(x∗)A = Df(x^*) is the matrix that encodes this local linear vector field: it tells you, for a small displacement uu from x∗x^*, the (leading-order) velocity AuAu that the true nonlinear field produces. Two pictures help fix the idea. First, think of the graph of a single nonlinear function y=f(x)y = f(x) near a root x∗x^*: the tangent line y≈f′(x∗)(x−x∗)y \approx f'(x^*)(x - x^*) is the one-dimensional shadow of exactly this construction, with the "Jacobian" reducing to the scalar derivative f′(x∗)f'(x^*). Second, in two dimensions, the Jacobian's eigenvalues determine whether nearby trajectories are pulled in (stable node/spiral), pushed out (unstable node/spiral), or pass through in a hyperbolic saddle pattern — the classification developed in Classification of Planar Fixed Points. The figures below show a concrete nonlinear vector field, its linearization at a fixed point, and how closely the two agree in a small neighborhood while diverging further out.

Figure
A saddle-type fixed point: nonlinear flow (left) vs. its linearization (right), both showing the same local eigen-directions
x y x* = (0,0) Nonlinear f(x,y) u₁ u₂ u = 0 Linear A·u, A = Df(x*)
Figure
Scalar picture: tangent line as the n=1 shadow of the Jacobian, f(x) = x - x³ near x* = 1
x f(x) x* = 1, f(x*) = 0 tangent: f'(1)(x−1) = −2(x−1)
TheoremLinearization as the first-order flow approximation

Let f:U→Rnf: U \to \mathbb{R}^n be C2C^2 on an open set U⊆RnU \subseteq \mathbb{R}^n containing a fixed point x∗x^* (so f(x∗)=0f(x^*) = 0), and let A=Df(x∗)A = Df(x^*). Write x˙=f(x)\dot x = f(x) and let u=x−x∗u = x - x^*. Then:

(a) (Local Taylor approximation.) There is a function R:U−x∗→RnR: U - x^* \to \mathbb{R}^n such that

f(x∗+u)=Au+R(u),with ∥R(u)∥=O(∥u∥2) as u→0.f(x^*+u) = Au + R(u), \qquad \text{with } \lVert R(u)\rVert = O(\lVert u \rVert^2) \text{ as } u \to 0.

Equivalently R(u)/∥u∥→0R(u)/\lVert u\rVert \to 0 as u→0u \to 0, so AuAu is the unique linear map approximating f(x∗+u)f(x^*+u) to first order.

(b) (Hyperbolicity controls local stability.) If every eigenvalue of AA has negative real part, then x∗x^* is asymptotically stable for the full nonlinear system x˙=f(x)\dot x = f(x). If at least one eigenvalue of AA has positive real part, then x∗x^* is unstable for the full nonlinear system. (Part (b) is the Hartman–Grobman / Lyapunov linearization principle; hypothesis "hyperbolic," i.e. no eigenvalue with zero real part, is essential and is discussed further below and in Hartman-Grobman and Hyperbolicity.)

Proof

We prove (a) directly from multivariable Taylor's theorem, then prove (b) using a Lyapunov-function argument built from the linear theory.

Proof of (a). Since f∈C2(U)f \in C^2(U), Taylor's theorem with the Lagrange (or Peano) remainder form applied to each component fif_i at the point x∗x^* gives, for uu small enough that the segment from x∗x^* to x∗+ux^*+u lies in UU,

fi(x∗+u)=fi(x∗)+∇fi(x∗)⋅u+12uTHi(ξi)u,f_i(x^*+u) = f_i(x^*) + \nabla f_i(x^*)\cdot u + \tfrac12 u^{\mathsf T} H_i(\xi_i) u,

where HiH_i is the Hessian matrix of fif_i and ξi\xi_i is a point on the segment [x∗,x∗+u][x^*, x^*+u] depending on ii and uu (this is the multivariable mean value form of the remainder, valid since fi∈C2f_i \in C^2). Since x∗x^* is a fixed point, fi(x∗)=0f_i(x^*) = 0. Stacking the nn equations,

f(x∗+u)=Au+R(u),Ri(u)=12uTHi(ξi)u,f(x^*+u) = Au + R(u), \qquad R_i(u) = \tfrac12 u^{\mathsf T} H_i(\xi_i) u,

where A=Df(x∗)A = Df(x^*) is exactly the matrix with rows ∇fi(x∗)T\nabla f_i(x^*)^{\mathsf T}, i.e. the Jacobian as defined above. To bound RR, note f∈C2(U)f \in C^2(U) means each Hessian HiH_i is continuous, hence bounded on a compact neighborhood Bρ(x∗)‾⊂U\overline{B_\rho(x^*)} \subset U for some ρ>0\rho>0: there is M>0M>0 with ∥Hi(y)∥op≤M\lVert H_i(y)\rVert_{\mathrm{op}} \le M for all y∈Bρ(x∗)‾y \in \overline{B_\rho(x^*)} and all i=1,…,ni=1,\dots,n. Then for ∥u∥≤ρ\lVert u \rVert \le \rho,

∣Ri(u)∣=∣12uTHi(ξi)u∣≤12M∥u∥2,|R_i(u)| = \left|\tfrac12 u^{\mathsf T} H_i(\xi_i) u\right| \le \tfrac12 M \lVert u \rVert^2,

by the definition of operator norm. Summing in quadrature over i=1,…,ni = 1,\dots,n,

∥R(u)∥≤n2M∥u∥2=O(∥u∥2).\lVert R(u) \rVert \le \tfrac{\sqrt n}{2} M \lVert u \rVert^2 = O(\lVert u\rVert^2).

Dividing by ∥u∥\lVert u \rVert and letting u→0u \to 0 shows R(u)/∥u∥→0R(u)/\lVert u \rVert \to 0, which is precisely the statement that u↦Auu \mapsto Au is the (unique) best linear approximation to u↦f(x∗+u)u \mapsto f(x^*+u) at u=0u=0: uniqueness follows because if BuBu also satisfies f(x∗+u)−Bu=o(∥u∥)f(x^*+u) - Bu = o(\lVert u\rVert), subtracting gives (A−B)u=o(∥u∥)(A-B)u = o(\lVert u \rVert), and evaluating along u=teju = t e_j for fixed standard basis vector eje_j and t→0t \to 0 forces column jj of A−BA - B to be zero for every jj, i.e. A=BA = B. This proves (a).

Proof of (b), stable case. Assume every eigenvalue of AA has real part ≤−α\le -\alpha for some α>0\alpha > 0. We invoke the standard fact from linear system theory (a consequence of the Jordan form / matrix exponential estimates, not dependent on the theorem being proved here) that there exists a symmetric positive definite matrix PP solving the Lyapunov equation

ATP+PA=−I,A^{\mathsf T} P + P A = -I,

which exists precisely because AA has no eigenvalues on the imaginary axis with the correct sign structure (all eigenvalues in the open left half-plane); this is the standard continuous Lyapunov equation existence result, proved by taking P=∫0∞eATteAt dtP = \int_0^\infty e^{A^{\mathsf T} t} e^{At}\,dt, an integral that converges absolutely because ∥eAt∥≤Ce−αt\lVert e^{At}\rVert \le Ce^{-\alpha t} for some C>0C>0 when all eigenvalues of AA have real part ≤−α<0\le -\alpha<0. Define V(u)=uTPuV(u) = u^{\mathsf T} P u on a neighborhood of 00. Since PP is symmetric positive definite, V(u)>0V(u) > 0 for u≠0u \ne 0 and V(0)=0V(0)=0, and there are constants 0<λmin⁡≤λmax⁡0 < \lambda_{\min} \le \lambda_{\max} (the extreme eigenvalues of PP) with

λmin⁡∥u∥2≤V(u)≤λmax⁡∥u∥2.\lambda_{\min}\lVert u\rVert^2 \le V(u) \le \lambda_{\max}\lVert u\rVert^2 .

Compute the derivative of VV along trajectories of the full nonlinear system u˙=f(x∗+u)=Au+R(u)\dot u = f(x^*+u) = Au + R(u):

V˙(u)=u˙TPu+uTPu˙=uT(ATP+PA)u+2uTPR(u)=−∥u∥2+2uTPR(u).\dot V(u) = \dot u^{\mathsf T} P u + u^{\mathsf T} P \dot u = u^{\mathsf T}(A^{\mathsf T} P + PA)u + 2u^{\mathsf T} P R(u) = -\lVert u \rVert^2 + 2 u^{\mathsf T} P R(u).

By part (a), given ε>0\varepsilon>0 there is δ>0\delta>0 such that ∥R(u)∥≤ε∥u∥\lVert R(u)\rVert \le \varepsilon \lVert u \rVert whenever ∥u∥≤δ\lVert u \rVert \le \delta. Using Cauchy–Schwarz and the operator norm bound ∥P∥op=λmax⁡\lVert P \rVert_{\mathrm{op}} = \lambda_{\max},

∣2uTPR(u)∣≤2λmax⁡∥u∥ ∥R(u)∥≤2λmax⁡ε∥u∥2.|2u^{\mathsf T} P R(u)| \le 2\lambda_{\max}\lVert u\rVert\,\lVert R(u)\rVert \le 2\lambda_{\max}\varepsilon \lVert u\rVert^2 .

Choose ε=1/(4λmax⁡)\varepsilon = 1/(4\lambda_{\max}) and let δ\delta be the corresponding radius; then for 0<∥u∥≤δ0<\lVert u\rVert \le \delta,

V˙(u)≤−∥u∥2+12∥u∥2=−12∥u∥2<0.\dot V(u) \le -\lVert u\rVert^2 + \tfrac12 \lVert u \rVert^2 = -\tfrac12 \lVert u \rVert^2 < 0.

Thus VV is a strict Lyapunov function on the ball Bδ(0)B_\delta(0): V>0V>0 away from 00, V(0)=0V(0)=0, and V˙<0\dot V < 0 on trajectories away from 00. By the standard Lyapunov stability theorem (asymptotic stability from a strict Lyapunov function, itself independent of the present theorem), u=0u=0, i.e. x=x∗x = x^*, is an asymptotically stable fixed point of x˙=f(x)\dot x = f(x). This proves the stable case of (b).

Proof of (b), unstable case. Suppose AA has an eigenvalue with positive real part. Decompose Rn=Eu⊕Ecs\mathbb{R}^n = E^u \oplus E^{cs} where EuE^u is the (nontrivial) sum of generalized eigenspaces for eigenvalues with positive real part and EcsE^{cs} the sum for the rest. A standard instability argument (Chetaev's instability theorem, again a linear-algebra-based fact not circular with the present statement) constructs, from the spectral decomposition of AA restricted to EuE^u, a quadratic function WW that is positive on a cone containing directions in EuE^u and satisfies W˙>0\dot W > 0 near 00 along nonlinear trajectories in that cone, using the same first-order remainder bound from part (a) to control the nonlinear correction exactly as above. The existence of points arbitrarily close to 00 inside this cone with W>0W>0, together with W˙>0\dot W>0 forcing WW to increase along the trajectory, shows the trajectory must leave any small ball around 00, which is precisely instability. This establishes the unstable case of (b), completing the proof. ■\blacksquare

∎

Worked examples

Example 1 (asymmetric planar system — asymptotically stable spiral). Consider

x˙=−x+2y−x3,y˙=−3x−y+xy.\dot x = -x + 2y - x^3, \qquad \dot y = -3x - y + xy.

Check (0,0)(0,0) is a fixed point: plugging in x=y=0x=y=0 gives x˙=0,y˙=0\dot x = 0, \dot y=0. Good. Compute the Jacobian. With f1(x,y)=−x+2y−x3f_1(x,y) = -x+2y-x^3 and f2(x,y)=−3x−y+xyf_2(x,y) = -3x-y+xy,

∂f1∂x=−1−3x2,∂f1∂y=2,∂f2∂x=−3+y,∂f2∂y=−1+x.\frac{\partial f_1}{\partial x} = -1 - 3x^2,\quad \frac{\partial f_1}{\partial y} = 2,\quad \frac{\partial f_2}{\partial x} = -3+y,\quad \frac{\partial f_2}{\partial y} = -1+x.

At (x,y)=(0,0)(x,y)=(0,0):

A=Df(0,0)=(−12−3−1).A = Df(0,0) = \begin{pmatrix} -1 & 2 \\ -3 & -1 \end{pmatrix}.

Find eigenvalues from det⁡(A−λI)=0\det(A - \lambda I) = 0:

(−1−λ)2−(2)(−3)=0  ⟹  (λ+1)2+6=0  ⟹  λ=−1±i6.(-1-\lambda)^2 - (2)(-3) = 0 \;\Longrightarrow\; (\lambda+1)^2 + 6 = 0 \;\Longrightarrow\; \lambda = -1 \pm i\sqrt6.

Both eigenvalues have real part −1<0-1 < 0. By Theorem part (b), the origin is asymptotically stable for the nonlinear system, and moreover the negative real part with nonzero imaginary part predicts (foreshadowing Classification of Planar Fixed Points) that trajectories spiral into the origin rather than approaching along fixed directions. Note this conclusion required only the linear term; the cubic and quadratic terms −x3-x^3 and xyxy affect the picture far from the origin but are invisible to the Jacobian at (0,0)(0,0), exactly as guaranteed by the O(∥u∥2)O(\lVert u\rVert^2) remainder bound in part (a).

Example 2 (three-dimensional system with an explicit non-origin fixed point). Consider the system

x˙=y−x2,y˙=x+z−4,z˙=xy−z−1.\dot x = y - x^2, \qquad \dot y = x + z - 4, \qquad \dot z = xy - z - 1.

First, locate a fixed point other than the origin (the origin is not a fixed point here, since y˙∣0=−4≠0\dot y|_{0} = -4 \ne 0). Try to solve

y=x2,x+z=4,xy−z=1.y = x^2, \qquad x + z = 4, \qquad xy - z = 1.

From the second equation, z=4−xz = 4-x. Substitute into the third: x⋅x2−(4−x)=1x\cdot x^2 - (4-x) = 1, i.e. x3+x−5=0x^3 + x - 5 = 0. Testing x=1x=1: 1+1−5=−3≠01+1-5=-3\ne0. Testing x=32x=\tfrac{3}{2}: 27/8+3/2−5=3.375+1.5−5=−0.12527/8 + 3/2 - 5 = 3.375+1.5-5 = -0.125, close. Rather than hunt numerically, pick a cleaner route: choose the fixed point to be exactly x∗=1.5x^*=1.5 is inconvenient; instead verify x=5−x3x=\sqrt[3]{5-x} has a root near 1.5161.516, but to keep the Jacobian computation exact we instead exhibit the linearization at an arbitrary fixed point (x∗,y∗,z∗)(x^*,y^*,z^*) satisfying the three equations above, which suffices to illustrate the method without needing a closed-form root. Compute the Jacobian symbolically at a fixed point (x∗,y∗,z∗)(x^*,y^*,z^*) with y∗=(x∗)2y^*=(x^*)^2, z∗=4−x∗z^*=4-x^*:

Df(x,y,z)=(−2x10101yx−1),A=Df(x∗,y∗,z∗)=(−2x∗10101(x∗)2x∗−1).Df(x,y,z) = \begin{pmatrix} -2x & 1 & 0 \\ 1 & 0 & 1 \\ y & x & -1 \end{pmatrix}, \qquad A = Df(x^*,y^*,z^*) = \begin{pmatrix} -2x^* & 1 & 0 \\ 1 & 0 & 1 \\ (x^*)^2 & x^* & -1 \end{pmatrix}.

To make this concrete, take the numerically located root of x3+x−5=0x^3+x-5=0: since the function g(x)=x3+x−5g(x)=x^3+x-5 satisfies g(1.5)=−0.125g(1.5) = -0.125 and g(1.55)=1.553+1.55−5=3.724+1.55−5=0.274g(1.55) = 1.55^3+1.55-5 = 3.724+1.55-5=0.274, there is a root x∗∈(1.5,1.55)x^*\in(1.5,1.55) by the Intermediate Value Theorem (continuity of gg); refine: g(1.52)=3.512+1.52−5=0.032g(1.52)=3.512+1.52-5=0.032, g(1.51)=3.443+1.51−5=−0.047g(1.51)=3.443+1.51-5=-0.047, so x∗≈1.516x^* \approx 1.516 (to three decimals). Then y∗=(x∗)2≈2.298y^* = (x^*)^2 \approx 2.298 and z∗=4−x∗≈2.484z^* = 4-x^*\approx 2.484. Numerically,

A≈(−3.032101012.2981.516−1).A \approx \begin{pmatrix} -3.032 & 1 & 0 \\ 1 & 0 & 1 \\ 2.298 & 1.516 & -1 \end{pmatrix}.

The characteristic polynomial is det⁡(A−λI)=0\det(A-\lambda I)=0; expanding along the first row,

det⁡(A−λI)=(−3.032−λ)[(−λ)(−1−λ)−1⋅1.516]−1⋅[1⋅(−1−λ)−1⋅2.298].\det(A-\lambda I) = (-3.032-\lambda)\big[(-\lambda)(-1-\lambda) - 1\cdot 1.516\big] - 1\cdot\big[1\cdot(-1-\lambda) - 1\cdot 2.298\big].

Compute the bracket: (−λ)(−1−λ)−1.516=λ+λ2−1.516(-\lambda)(-1-\lambda) - 1.516 = \lambda + \lambda^2 - 1.516. And the second bracket: −1−λ−2.298=−3.298−λ-1-\lambda-2.298 = -3.298-\lambda. So

det⁡(A−λI)=(−3.032−λ)(λ2+λ−1.516)−(−3.298−λ)=(−3.032−λ)(λ2+λ−1.516)+3.298+λ.\det(A-\lambda I) = (-3.032-\lambda)(\lambda^2+\lambda-1.516) - (-3.298-\lambda) = (-3.032-\lambda)(\lambda^2+\lambda-1.516) + 3.298+\lambda.

Expanding the product:

(−3.032)(λ2+λ−1.516)−λ(λ2+λ−1.516)=−3.032λ2−3.032λ+4.597−λ3−λ2+1.516λ,(-3.032)(\lambda^2+\lambda-1.516) - \lambda(\lambda^2+\lambda-1.516) = -3.032\lambda^2-3.032\lambda+4.597 -\lambda^3-\lambda^2+1.516\lambda,

so the product equals −λ3−4.032λ2−1.516λ+4.597-\lambda^3 -4.032\lambda^2 -1.516\lambda + 4.597. Adding 3.298+λ3.298+\lambda:

det⁡(A−λI)=−λ3−4.032λ2−0.516λ+7.895.\det(A-\lambda I) = -\lambda^3 - 4.032\lambda^2 - 0.516\lambda + 7.895.

Setting this to zero and multiplying by −1-1: λ3+4.032λ2+0.516λ−7.895=0\lambda^3 + 4.032\lambda^2 + 0.516\lambda - 7.895 = 0. Evaluate at λ=1\lambda=1: 1+4.032+0.516−7.895=−2.3471+4.032+0.516-7.895 = -2.347. At λ=1.3\lambda=1.3: 2.197+6.814+0.671−7.895=1.7872.197+6.814+0.671-7.895 = 1.787. So there is a real root between 11 and 1.31.3 — a positive real eigenvalue. By Theorem part (b), since at least one eigenvalue has positive real part, this fixed point is unstable. This illustrates that linearization at a non-origin, numerically located fixed point works exactly the same way, and again the instability conclusion follows purely from the sign pattern of the Jacobian's spectrum.

Counterexamplewhy hyperbolicity is needed

Drop the hypothesis that A=Df(x∗)A=Df(x^*) has no eigenvalue with zero real part, and consider

x˙=−y+μx(x2+y2),y˙=x+μy(x2+y2),\dot x = -y + \mu x(x^2+y^2), \qquad \dot y = x + \mu y(x^2+y^2),

for a real parameter μ\mu. The origin is a fixed point for every μ\mu. The Jacobian at the origin is

A=(0−110),eigenvalues λ=±i,A = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}, \qquad \text{eigenvalues } \lambda = \pm i,

independent of μ\mu: purely imaginary, so the origin is not hyperbolic, regardless of μ\mu. Convert to polar coordinates x=rcos⁡θx=r\cos\theta, y=rsin⁡θy=r\sin\theta. A direct computation using rr˙=xx˙+yy˙r\dot r = x\dot x+y\dot y gives

rr˙=x(−y+μxr2)+y(x+μyr2)=μr2(x2+y2)=μr4,r\dot r = x(-y+\mu x r^2) + y(x+\mu y r^2) = \mu r^2 (x^2+y^2) = \mu r^4,

so r˙=μr3\dot r = \mu r^3, and similarly θ˙=1\dot\theta = 1 (uniform rotation). Thus:

  • If μ<0\mu < 0: r˙=μr3<0\dot r = \mu r^3 < 0 for r>0r>0, so r(t)→0r(t)\to 0 — the origin is asymptotically stable.
  • If μ>0\mu > 0: r˙=μr3>0\dot r = \mu r^3 > 0 for r>0r>0, so r(t)r(t) grows — the origin is unstable.
  • If μ=0\mu = 0: r˙=0\dot r = 0, every circle r=constr=\text{const} is an invariant periodic orbit — the origin is a stable but not asymptotically stable center.

All three of these genuinely different nonlinear behaviors occur while the linearization AA is exactly the same rotation matrix with eigenvalues ±i\pm i in every case. The linear system u˙=Au\dot u = Au alone (a center, with all orbits closed circles) cannot distinguish μ<0\mu<0 from μ>0\mu>0 from μ=0\mu=0: it predicts a center regardless of μ\mu, which is correct only when μ=0\mu=0 and wrong (as a stability classification) otherwise. This shows precisely why Theorem part (b) requires hyperbolicity: when eigenvalues lie exactly on the imaginary axis, the O(∥u∥2)O(\lVert u \rVert^2) remainder term — however small — determines the true qualitative behavior, and the linearization is uninformative about stability.

Scope, edge cases & common traps

SituationWhat linearization tells youWhat it does NOT tell you
All eigenvalues of AA have negative real partAsymptotic stability of x∗x^* for the nonlinear system (proved above)Exact trajectories, basin of attraction size
At least one eigenvalue has positive real partInstability of x∗x^* for the nonlinear system (proved above)Which directions trajectories leave along, globally
Some eigenvalue has zero real part (non-hyperbolic)Nothing about nonlinear stabilityStability is genuinely undetermined by AA alone; see counterexample
ff only C1C^1, not C2C^2Part (a) with a weaker o(∥u∥)o(\lVert u\rVert) remainder still holds (differentiability alone gives this)The clean quadratic remainder bound used in part (b)'s proof requires C2C^2 (or Lipschitz DfDf)
Fixed point at x∗≠0x^*\ne 0Same theorem, with u=x−x∗u = x-x^*The Jacobian must be evaluated at x∗x^*, not at the origin

The mathematician behind it

Origins · 1 minAleksandr Mikhailovich Lyapunov and how this was discovered
AM
Aleksandr Mikhailovich Lyapunov
1857–1918 · Imperial Russian mathematics, stability theory of motion

Lyapunov's 1892 doctoral thesis, The General Problem of the Stability of Motion, established both the linearization principle proved above (his "first method," valid precisely when the linearization is hyperbolic) and the direct, energy-like method of Lyapunov functions used in the proof of Theorem part (b). He worked in the tradition of Chebyshev's school in St. Petersburg, motivated by problems in celestial mechanics and the stability of rotating fluid masses, and his methods remain the standard toolkit for stability analysis in every modern textbook on nonlinear dynamics and control theory.

History
- 1892 — Lyapunov's thesis proves the linearization ("first method") and direct ("second method," Lyapunov functions) stability theorems.
- 1901 — Oskar Perron begins developing what becomes Perron–Frobenius and related spectral techniques underlying rigorous linearization estimates.
- 1929 — Oskar Perron formalizes conditions on nonlinear perturbations of linear systems (Perron's theorem on stability under perturbation), sharpening Lyapunov's linearization result.
- 1959–1963 — Philip Hartman and David Grobman independently prove the topological conjugacy theorem for hyperbolic fixed points, showing the full local phase portrait (not just stability) matches the linearization.
- 1970s–present — Linearization and the Jacobian become foundational computational tools across bifurcation theory, control engineering, and numerical continuation software.
Is the Jacobian the same thing as the derivative?

Yes, in the multivariable sense: for a vector-valued function f:Rn→Rnf:\mathbb{R}^n \to \mathbb{R}^n, the Jacobian matrix is the total derivative Df(x∗)Df(x^*), the unique linear map satisfying f(x∗+u)=f(x∗)+Df(x∗)u+o(∥u∥)f(x^*+u) = f(x^*)+Df(x^*)u + o(\lVert u\rVert). When n=1n=1 it reduces to the familiar scalar derivative f′(x∗)f'(x^*). The word "Jacobian" is often reserved for the matrix representation (the array of partial derivatives ∂fi/∂xj\partial f_i/\partial x_j), while "derivative" or "differential" refers to the linear map itself; they carry identical information once a basis (coordinates) is fixed, which is standard in Rn\mathbb{R}^n.

Why does linearization only work near the fixed point, and how near is "near"?

The remainder bound ∥R(u)∥=O(∥u∥2)\lVert R(u)\rVert = O(\lVert u\rVert^2) is a statement about the limit as u→0u\to0; it guarantees the ratio R(u)/∥u∥R(u)/\lVert u\rVert vanishes in the limit but gives no explicit radius of validity without more information (a bound MM on the Hessian, as used in the proof). In practice the linearization is trustworthy in a neighborhood whose size depends on how large the second-derivative terms are and how far the eigenvalues of AA are from the imaginary axis (the "spectral gap"); the closer the fixed point is to non-hyperbolicity, the smaller this neighborhood tends to be, since the linear restoring force weakens relative to the nonlinear terms.

Does a stable linearization make the nonlinear system stable everywhere, or just locally?

Only locally. Theorem part (b) is inherently a statement about a (possibly small) neighborhood of x∗x^*; the proof produces a specific radius δ\delta inside which the Lyapunov function argument works, but says nothing about behavior far from x∗x^*. A system can have a locally asymptotically stable fixed point while trajectories starting far away diverge to infinity or approach a different attractor; establishing global stability requires separate techniques, such as constructing a Lyapunov function valid on the entire domain.

What if the fixed point is not at the origin — do I need a different formula?

No new formula is needed, only the substitution u=x−x∗u = x - x^*. You still compute the same partial-derivative matrix DfDf, but you evaluate it at x∗x^* rather than at 00, and the linear system is u˙=Df(x∗)u\dot u = Df(x^*)u in the shifted coordinate uu. Example 2 above illustrates this: the fixed point was located by solving three equations simultaneously, and the Jacobian was evaluated at that specific point, not at the origin, since the origin was not even a fixed point of that system.

Can linearization fail even when hyperbolicity holds, due to higher-dimensional subtleties?

No — this is exactly the content of the Hartman–Grobman theorem (foreshadowed here, proved in full in Hartman-Grobman and Hyperbolicity): whenever the fixed point is hyperbolic, there is a continuous, invertible change of coordinates (a homeomorphism) taking the nonlinear flow near x∗x^* to the linear flow of u˙=Au\dot u = Au near 00, valid in any finite dimension nn. The stability conclusion proved in this article (part (b)) is a consequence of that stronger topological statement, though we proved stability directly via Lyapunov functions without needing the full conjugacy.

Why did the proof need C² instead of just C¹?

Differentiability alone (C1C^1) gives R(u)=o(∥u∥)R(u) = o(\lVert u \rVert), which is enough to define the linearization and prove part (a) in its basic form. However the rate estimate R(u)=O(∥u∥2)R(u) = O(\lVert u\rVert^2) used in the stability proof (part (b)) requires a uniform bound on second derivatives (or at least a Lipschitz bound on DfDf) so that the quadratic remainder in Taylor's theorem is controlled uniformly near x∗x^*; without this, the Lyapunov function argument's cancellation −∥u∥2+2uTPR(u)<0-\lVert u\rVert^2 + 2u^{\mathsf T}PR(u) < 0 cannot be guaranteed to hold on a fixed-size neighborhood.

How does this relate to Newton's method?

Newton's method for solving f(x)=0f(x)=0 iterates xk+1=xk−[Df(xk)]−1f(xk)x_{k+1} = x_k - [Df(x_k)]^{-1}f(x_k), which is precisely the linear approximation from part (a) solved for the root of the linearized equation f(xk)+Df(xk)(x−xk)≈0f(x_k) + Df(x_k)(x-x_k) \approx 0. The same Taylor expansion and remainder estimate used to justify linearization near a dynamical fixed point also justifies the quadratic convergence rate of Newton's method near a simple root, making this article's Theorem part (a) a shared foundation for both stability theory and root-finding numerics.

Before this

Where this leads