Shows that the trace T and the determinant D of a real 2×2 matrix decide the entire orbit picture of x˙=Ax near the origin — saddle, node, spiral, centre, star node, degenerate node, or a line of fixed points — and places each type in the (T,D) plane relative to the two axes and the parabola T2=4D.
A fixed point is a point x∗ with Ax∗=0; the origin is always one, and it is isolated exactly when A is invertible. The two numbers that will decide everything are
T=trA=a+d,D=detA=ad−bc,Δ=T2−4D,
the trace, the determinant and the discriminant. An eigenline is a line Rv spanned by a real eigenvector v of A; it is invariant, because x(0)∈Rv forces x(t)=eλtx(0)∈Rv.
x∗ is Lyapunov stable if every neighbourhood U of x∗ contains a neighbourhood U1 of x∗ such that every solution starting in U1 exists and stays in U for all t>0.
x∗ is attracting if for some δ>0 every solution with ∣x(0)−x∗∣<δ satisfies x(t)→x∗ as t→∞.
x∗ is asymptotically stable if it is both; neutrally stable if it is Lyapunov stable but not attracting; unstable if it is not Lyapunov stable.
The two halves are independent: Strogatz (Strogatz (1994), p. 129) notes that a fixed point may be attracting without being Lyapunov stable, and the centre below is Lyapunov stable without being attracting.
Definition(The six types)AdaptedAdapted from Strogatz (1994), §5.2, pp. 132-136
Suppose A is invertible, so that the origin is an isolated fixed point. It is called a
saddle if the eigenvalues are real with opposite signs;
node if they are real, distinct and of the same sign;
spiral if they are a complex conjugate pair with nonzero real part;
centre if they are a nonzero purely imaginary pair;
star node if A=λI with λ=0;
degenerate node if the eigenvalue is real and repeated but the eigenspace is only one-dimensional.
These six cases are exhaustive and mutually exclusive. The theorem below says which region of the (T,D) plane each occupies, what the orbits look like, and which are stable.
Let A be a real 2×2 matrix with T=trA, D=detA, Δ=T2−4D. The eigenvalues of A are
λ±=21(T±Δ),soλ++λ−=T,λ+λ−=D.
D<0: saddle, unstable. The eigenvalues are real with λ−<0<λ+. The eigenline of λ− (the stable line) consists of the origin and orbits tending to the origin as t→+∞; the eigenline of λ+ (the unstable line) of orbits tending to it as t→−∞. Every other orbit leaves every bounded set in both time directions.
D>0 and Δ>0: node. The eigenvalues are real, distinct and both of the sign of T (in particular T=0). The node is asymptotically stable if T<0 and unstable if T>0. Call the eigenvalue of smaller modulus slow and the other fast. For T<0 every orbit tends to the origin; every orbit off the fast eigenline arrives at the origin tangent to the slow eigendirection, and every orbit off the slow eigenline has its direction tending to the fast eigendirection as t→−∞. The two halves exclude different orbits: an orbit on the slow eigenline satisfies the first and not the second, and an orbit on the fast eigenline the second and not the first. For T>0 the same holds with time reversed.
Δ<0: spiral or centre. Then D>T2/4≥0 automatically. The eigenvalues are α±iω with
α=21T,ω=214D−T2>0.
In a suitable linear coordinate system every orbit winds about the origin with constant angular speed ω, its distance from the origin being multiplied by eαt. The sense of rotation is the same for every orbit, but it is not determined by T and D; the remark below gives two matrices with the same (T,D) that turn opposite ways. Hence the origin is an asymptotically stable stable spiral if T<0, an unstable spiral if T>0, and if T=0 a centre: every orbit other than the origin is then a closed curve of period 2π/ω=2π/D, and the origin is neutrally stable.
D>0 and Δ=0: star or degenerate node. The eigenvalue λ=T/2 is real, repeated and nonzero. If A=λI the origin is a star node and every orbit is an open ray from the origin along a straight line. Otherwise the eigenspace is a single line Rv, the origin is a degenerate node, and every orbit off that line has direction tending to ±v both as t→+∞ and as t→−∞. In both cases the origin is asymptotically stable if T<0 and unstable if T>0.
D=0: the fixed point is not isolated. At least one eigenvalue vanishes. If A=0 the fixed points fill the line kerA; if A=0 every point of the plane is fixed.
In the (T,D) plane this reads off at a glance: the open half-plane D<0 is saddles; above the T-axis the parabola Δ=0, that is D=T2/4, separates nodes (Δ>0, outside the parabola) from spirals (Δ<0, inside it), with star and degenerate nodes on the parabola itself; the sign of T gives stability throughout the region D>0; and the positive D-axis T=0 carries the centres.
Idea of the proofOriginal
Everything comes from one line of algebra: the characteristic polynomial of a 2×2 matrix is λ2−Tλ+D, so the trace and the determinant are the sum and the product of the eigenvalues, and nothing else about A enters the eigenvalues at all. The product tells us whether the eigenvalues have the same sign; the sign of the discriminant tells us whether they are real; and the sum then fixes which sign they share. Each of the resulting eigenvalue configurations has a general solution written down in Linear systems and the matrix exponential — a combination of two real exponentials, or an exponential times a rotation, or an exponential times a polynomial of degree one — and the shape of the orbits is read directly off that formula by asking which term dominates as t→±∞.
Its roots are λ±=21(T±Δ), and comparing with (λ−λ+)(λ−λ−) gives λ++λ−=T and λ+λ−=D. This is the only place the entries of A are used; from here on we argue with T and D.
Step 2 (case 1, D<0). Subtracting a negative number enlarges T2, so Δ=T2−4D>T2≥0. A positive discriminant makes both roots real, and Δ>∣T∣. Hence λ+=21(T+Δ)>0 and λ−=21(T−Δ)<0. The eigenvalues are distinct, so eigenvectors v± for them are independent and, by Linear systems and the matrix exponential, every solution is
x(t)=c+eλ+tv++c−eλ−tv−.
If c+=0 and c−=0 the orbit lies on Rv− and tends to 0 as t→+∞; symmetrically for c−=0. If both coefficients are nonzero, then because v+,v− are independent there is m>0 with ∣αv++βv−∣≥mmax(∣α∣,∣β∣) for all real α,β, so ∣x(t)∣≥m∣c+∣eλ+t→∞ as t→+∞ and ∣x(t)∣≥m∣c−∣eλ−t→∞ as t→−∞. In particular arbitrarily small initial data with c+=0 leave any fixed ball, so the origin is unstable.
Step 3 (case 2, D>0, Δ>0). The eigenvalues are real and distinct. Their product is D>0, so they share a sign; their sum is T, so that sign is the sign of T, and T=0. (Consistently, Δ>0 and D>0 force T2>4D>0.) Take T<0; the case T>0 follows by replacing t by −t, which replaces A by −A and so T by −T while leaving D and Δ unchanged. So let λf<λs<0, with eigenvectors vf,vs; here ∣λs∣<∣λf∣, which is why λs is called slow. Every solution x(t)=cseλstvs+cfeλftvf tends to 0, and the origin is asymptotically stable because both exponentials are bounded by 1 for t≥0, so ∣x(t)∣≤C∣x(0)∣ with C depending only on the eigenbasis. For the tangency, suppose cs=0 and divide by the leading scalar:
since λf−λs<0. Rescaling by a positive number does not change the direction of a vector, so the unit vector x(t)/∣x(t)∣ converges to ±vs/∣vs∣: the orbit reaches the origin tangent to the slow eigendirection. This half of the claim used cs=0, that is, the orbit is off the fast eigenline.
For the backward-time half suppose instead that cf=0, that is, the orbit is off the slow eigenline, and divide by cfeλft:
since λs−λf>0 makes the exponential tend to 0 as t→−∞; so the unit vector x(t)/∣x(t)∣ converges to ±vf/∣vf∣. The two exclusions are genuinely different, and neither can be dropped: on the slow eigenline (cf=0) the direction of the orbit is ±vs/∣vs∣ for every t and never approaches vf, and on the fast eigenline (cs=0) it is ±vf/∣vf∣ for every t and never approaches vs.
Step 4 (case 3, Δ<0). From T2−4D<0 we get D>T2/4≥0, so D>0 is automatic here. The eigenvalues are λ±=α±iω with α=T/2 and ω=21−Δ=214D−T2>0. Let v=p+iq be a complex eigenvector for λ+; by Linear systems and the matrix exponential the real vectors p and q are independent and the general real solution is
x(t)=eαt[c1(cosωtp−sinωtq)+c2(sinωtp+cosωtq)].
Collecting the coefficients of p and of q, the coordinate pair of x(t) in the basis (p,q) is
a rotation through the angle −ωt scaled by eαt. So in these coordinates ξ2+η2=e2αt(c12+c22) and the polar angle decreases at the constant rate ω. If α<0 the orbit spirals into the origin through infinitely many turns, and ∣x(t)∣≤Ceαt∣x(0)∣ for a constant C fixed by the basis (p,q), so the origin is asymptotically stable; if α>0 it spirals out and is unstable. If α=0 then ξ2+η2 is constant: in the (p,q) coordinates every orbit is a circle, hence in the original coordinates an ellipse, closed and of period 2π/ω, and ω=214D=D. Such a fixed point is Lyapunov stable, because ∣x(t)∣≤M∣x(0)∣ with M=max0≤s≤2π/ω∥esA∥ by periodicity, but not attracting, because x(2πk/ω)=x(0) for every integer k: it is neutrally stable.
Step 5 (case 4, D>0, Δ=0). Now λ+=λ−=T/2=:λ, and T2=4D>0 gives λ=0. If the eigenspace of λ is all of R2 then Ax=λx for every x, that is A=λI; then x(t)=eλtx(0), so each orbit is the open ray from the origin through x(0), traversed inward if λ<0 and outward if λ>0. Otherwise the eigenspace is a line Rv, and by Jordan normal form there is a vector w with (A−λI)w=v, and (v,w) is a basis. Solving in that basis gives
x(t)=eλt[(c1+c2t)v+c2w],
as one checks by differentiating and using Av=λv, Aw=λw+v. For c2=0 the bracket has length growing like ∣c2t∣ and direction
because (c1+c2t)/∣c1+c2t∣ equals sgn(c2) once t is large and positive and −sgn(c2) once t is large and negative, while c2w/∣c1+c2t∣→0 in both directions. The two limits are the two opposite directions along Rv, so the orbit is tangent to the eigenline at both ends. Since a polynomial times eλt tends to 0 as t→+∞ exactly when λ<0, the fixed point is asymptotically stable for T<0 and unstable for T>0.
Step 6 (case 5, D=0). Then λ+λ−=0, so A is singular and kerA={0}. Every x∗∈kerA satisfies Ax∗=0 and is therefore a fixed point, so the origin is not isolated. If A=0 then rankA=1 and kerA is a line; if A=0 then kerA=R2. □
D=0 is what makes the fixed point isolated. It is used in Steps 2-5 through λ+λ−=0. Drop it and the classification collapses: for A=(0010) (T=D=0) the whole x-axis consists of fixed points and x(t)=x(0)+tAx(0) moves along horizontal lines at constant speed, so the fixed points are unstable; for A=(−1000) (T=−1, D=0) every point of the y-axis is a fixed point and is neutrally stable. Both are excluded by D=0.
Δ=0 is what gives two independent eigenvectors. It is used in Steps 2-4 to write the solution in an eigenbasis. On the parabola Δ=0 the Jordan form is needed instead (Step 5), and the pair (T,D) no longer determines the portrait: (−100−1) and (−101−1) both have T=−2, D=1, but the first is a star node and the second a degenerate node.
(T,D) never decides the sense of rotation.(0−110) and (01−10) both have T=0, D=1, and both are centres of period 2π, but the first turns clockwise and the second counterclockwise. More generally, Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), p. 63) point out that a single point of the (T,D) plane carries infinitely many different matrices — the entries supply four parameters and (T,D) imposes two conditions — all with the same eigenvalue configuration but with phase portraits that may differ in ways such as the sense of rotation or the number of independent eigenvectors. Those are examples, not a complete list: two matrices with the same (T,D) also have, in general, different eigendirections, so their portraits agree only after a linear change of coordinates.
The two halves of the node tangency exclude different orbits. Tangency to the slow eigendirection as t→+∞ needs cs=0, while convergence of the direction to the fast eigendirection as t→−∞ needs cf=0, and each hypothesis fails on one of the two eigenlines: for A=diag(−1,−2) the orbit x(t)=(e−t,0) lies on the slow eigenline, so its direction is (1,0) for every t and tends to the slow, not the fast, eigendirection in backward time. Strogatz (Strogatz (1994), pp. 133, 136) hedges the forward-time half with "typically" and states the backward-time half for all trajectories; case 2 above instead names the excluded eigenline in each half, which is what makes both statements hold of every orbit they mention.
The centre is the fragile case. It sits on the boundary T=0 between the stable and unstable half of the region D>0, so an arbitrarily small change of A turns it into a spiral of either sign. What happens when nonlinear terms rather than linear ones are added is the subject of Linearization and the Hartman-Grobman theorem.
Method(Classifying a planar linear system)Original
Given a real 2×2 matrix A:
Compute T=trA and D=detA. If D=0, stop: the fixed points are not isolated, and they fill kerA (or the plane, if A=0).
If D<0, the origin is a saddle and unstable. Its stable and unstable lines are the eigenlines of λ∓=21(T∓T2−4D).
Otherwise D>0; compute Δ=T2−4D. The sign of T now gives the stability: asymptotically stable for T<0, unstable for T>0, neutrally stable for T=0.
If Δ>0: a node. Compute λ± and their eigenvectors; the one with the smaller ∣λ∣ is the slow eigendirection, to which the orbits are tangent at the origin.
If Δ<0: a spiral, or a centre when T=0. The decay rate is α=T/2 and the angular speed is ω=214D−T2, so one turn takes time 2π/ω and multiplies the distance to the origin by e2πα/ω.
If Δ=0: a star node when A=(T/2)I, otherwise a degenerate node, with orbits tangent to the single eigenline at both ends.
For a spiral or centre, get the sense of rotation from one vector of the field: evaluate Ae1, where e1=(1,0). Its second component is positive for counterclockwise rotation and negative for clockwise.
Example(A stable spiral, solved and measured)Original
Classify and solve x˙=Ax with
A=(−3−421),x(0)=(10).
Step 1 (invariants).T=−3+1=−2 and D=(−3)(1)−(2)(−4)=−3+8=5. Hence Δ=T2−4D=4−20=−16<0. By case 3 of the theorem the origin is a spiral, and T<0 makes it a stable one: asymptotically stable.
Step 2 (rates).α=T/2=−1 and ω=2116=2, so the eigenvalues are −1±2i. Indeed λ2+2λ+5=(λ+1)2+4 vanishes at λ=−1±2i. One turn takes 2π/ω=π units of time.
Step 3 (eigenvector). For λ=−1+2i,
A−λI=(−2−2i−422−2i),
and the first row gives (−2−2i)v1+2v2=0, that is v2=(1+i)v1. Take v=(1,1+i); the second row is satisfied too, since −4+(2−2i)(1+i)=−4+4=0. So p=Rev=(1,1) and q=Imv=(0,1).
Step 4 (the two real solutions). With α=−1, ω=2 the formula of Step 4 of the proof gives the basis
Step 5 (fit the data).u1(0)=(1,1) and u2(0)=(0,1), so c1(1,1)+c2(0,1)=(1,0) gives c1=1, c2=−1 and
x(t)=e−t(cos2t−sin2t,−2sin2t).
Check: x˙1=e−t(−3cos2t−sin2t), and −3x1+2x2=e−t(−3cos2t+3sin2t−4sin2t) agrees; x˙2=e−t(2sin2t−4cos2t), and −4x1+x2=e−t(−4cos2t+4sin2t−2sin2t) agrees.
Step 6 (sense of rotation and a number).Ae1=(−3,−4) has negative second component, so the orbit turns clockwise. After one full turn, at t=π,
x(π)=e−π(cos2π−sin2π,−2sin2π)=e−π(10),
so the orbit returns to the positive x1-axis with its distance to the origin multiplied by
e−π=0.0432139182…
That is, one revolution leaves 4.32% of the distance: the spiral is a tight one, as ∣α∣/ω=1/2 predicts.
ApplicationOriginal
The classification is the standard first move on any two-dimensional model, and the (T,D) plane is where the familiar regimes of the damped oscillator live. Writing y¨+2γy˙+ω02y=0 of The damped and driven harmonic oscillator as the planar system y˙=v, v˙=−ω02y−2γv gives
A=(0−ω021−2γ),T=−2γ,D=ω02,Δ=4(γ2−ω02).
So D>0 always, and increasing the damping from γ=0 drives the point (T,D)=(−2γ,ω02) horizontally across the (T,D) plane: it starts on the positive D-axis (a centre: the undamped oscillator), moves into the spiral region (underdamped, γ<ω0, that is quality factor Q=ω0/2γ>1/2), touches the parabola Δ=0 exactly at critical damping γ=ω0 (a degenerate node, since A is then not a multiple of the identity), and passes into the node region beyond it (overdamped). The three classical regimes are thus three regions of one picture, and the same statement applies verbatim to a series RLC circuit with 2γ=R/L, ω02=1/(LC), and to any lightly damped mechanical mode.
The same reading is what makes the classification useful for nonlinear models: by Linearization and the Hartman-Grobman theorem the Jacobian at an equilibrium of a nonlinear planar system is a 2×2 matrix, and its trace and determinant decide the local picture whenever the equilibrium is hyperbolic. In an ecological or chemical model, where the Jacobian entries are measurable response coefficients, the sign of D separates saddle-type coexistence states from genuine attractors and repellers, and the sign of T then says whether the state is approached or left; whether the approach oscillates is the sign of Δ.
Defines the Jacobian linearization at an equilibrium and the notion of a hyperbolic equilibrium, proves that an equilibrium whose Jacobian has only eigenvalues of negative real part is asymptotically stable with an explicit exponential rate, states the Hartman-Grobman theorem that the nonlinear flow near any hyperbolic equilibrium is topologically conjugate to its linearization, and shows by two examples what goes wrong when hyperbolicity fails.
An equilibrium x∗ is hyperbolic if no eigenvalue of DF(x∗) has zero real part. In the plane this excludes exactly the two borderline columns of Classification of linear planar systems: a centre (T=0, D>0) and any matrix with D=0. Strogatz (Strogatz (1994), p. 155) calls the name unfortunate, since it suggests a saddle, and records that it has nevertheless become standard.
Let x∗ be an equilibrium of x˙=F(x) with F of class C1 near x∗. If every eigenvalue of DF(x∗) has negative real part, then x∗ is asymptotically stable.
Let x∗ be a hyperbolic equilibrium of x˙=F(x), F of class C1, and let A=DF(x∗). Then on some neighbourhood of x∗ the flow of x˙=F(x) is topologically conjugate to the flow of u˙=Au near 0: there is a homeomorphism h of a neighbourhood of x∗ onto a neighbourhood of 0, with h(x∗)=0, carrying orbits of the nonlinear system onto orbits of the linear one and preserving the direction and the parametrisation of time.
Strogatz (Strogatz (1994), p. 155) states the same result as topological equivalence of the two local phase portraits, the homeomorphism being required only to send trajectories to trajectories and to preserve the sense of time. Either way the consequence is the one we use: at a hyperbolic equilibrium the local picture — saddle, node or spiral, and its stability — is exactly the one predicted by the Jacobian. Strogatz (Strogatz (1994), p. 151) puts it as the rule that a linearized saddle, node or spiral is a genuine saddle, node or spiral of the nonlinear system, while the borderline cases (centres, stars, degenerate nodes, non-isolated fixed points) may be changed by the nonlinear terms. Neither book proves the theorem; Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), p. 168) say the general case needs techniques beyond their scope, and prove it only for sinks and sources.
Idea of the proofOriginal
For a sink the whole proof is one derivative. In coordinates adapted to A the linear field points strictly inwards across every small circle about the equilibrium, which one sees by taking the dot product of Au with u: it equals −λx2−μy2 for a diagonal A, a negative multiple of ∣u∣2. The nonlinear remainder h is o(∣u∣), so its contribution to that dot product is o(∣u∣2): for small enough circles it cannot reverse the sign. Since the dot product of the field with the position vector is exactly 21dtd∣u∣2, we get a differential inequality for the squared distance, and Gronwall's elementary one-line version of it gives exponential decay. Hartman-Grobman says the same comparison succeeds at every hyperbolic equilibrium, not only at sinks, but there one must match an expanding and a contracting direction at once and the conjugating homeomorphism is no longer explicit.
We give the planar case; the argument in Rn is the same once A is put in a real canonical form whose symmetric part is negative definite.
Step 0 (normalisations). Replacing x by u=x−x∗ moves the equilibrium to the origin and does not change A. Replacing u by Pu for an invertible real P replaces A by P−1AP and h by P−1h(P⋅), which is still o(∣u∣) because P and P−1 are bounded. Both changes are homeomorphisms fixing the origin, so they preserve asymptotic stability; we may therefore assume x∗=0 and that A is in whichever real canonical form is convenient.
Step 1 (the three canonical forms and the inward inequality). Write u=(x,y), r2=x2+y2, and h=(h1,h2). The eigenvalues of A have negative real parts, so exactly one of the following holds, and in each case we compute Au⋅u.
Distinct real eigenvalues −λ<−μ<0. In an eigenbasis A=diag(−λ,−μ) and
Au⋅u=−λx2−μy2≤−μ(x2+y2)=−μr2,
since λ>μ.
A complex pair α±iβ with α<0, β=0. In the real canonical basis A=(α−ββα), so the off-diagonal contributions βxy and −βxy cancel and
Au⋅u=αr2.
Put μ=−α>0.
A repeated eigenvalue −λ<0. If the eigenspace is two-dimensional then A=−λI and Au⋅u=−λr2. If it is a line Rv, then by Jordan normal form there is w with (A+λI)w=v; in the basis (v,εw) the matrix is (−λ0ε−λ) for any ε>0 we like, whence
Au⋅u=−λr2+εxy≤(−λ+2ε)r2,
using ∣xy∣≤21r2. Choosing ε=λ gives the bound with μ=λ/2>0.
In every case there is μ>0 and a linear coordinate system in which Au⋅u≤−μr2 for all u.
Step 2 (the nonlinear terms cannot reverse the sign). Since ∣h(u)∣/∣u∣→0, there is δ>0 with Bδ⊂O and
∣h(u)∣≤2μ∣u∣for ∣u∣≤δ.
By the Cauchy-Schwarz inequality, ∣h(u)⋅u∣≤2μr2 there. Hence for a solution inside the ball,
21dtdr(t)2=u⋅u˙=Au⋅u+h(u)⋅u≤−μr2+2μr2=−2μr2.
Here differentiability of F at the equilibrium is what makes h small, and the canonical form of Step 1 is what makes the linear part strictly negative.
Step 3 (integrate the inequality). Let u(0) satisfy ∣u(0)∣≤δ and let u be the solution through it, which exists and is unique by Existence and uniqueness of solutions, F being C1 and hence locally Lipschitz. As long as the solution stays in Bδ, Step 2 gives dtd(eμtr2)=eμt(μr2+dtdr2)≤0, so
r(t)≤r(0)e−μt/2.
In particular r is non-increasing, so the solution never reaches the boundary of Bδ and, remaining in a compact subset of O, continues for all t≥0. The displayed bound therefore holds for all t≥0.
Step 4 (read off stability). Given a neighbourhood U of the origin, choose ρ≤δ with Bρ⊂U; by Step 3 every solution starting in Bρ stays in Bρ⊂U for all t≥0, so the origin is Lyapunov stable. And r(t)≤r(0)e−μt/2→0, so every solution starting in Bδ converges to the origin: the origin is attracting, hence asymptotically stable. □
Remark(What the proof gives beyond the statement)Original
Step 3 proves more than asymptotic stability. It exhibits a radius δ and a rate: in the adapted coordinates the closed ball of radius δ is positively invariant and lies in the basin of attraction, and the distance to the equilibrium decays at least like e−μt/2, where −μ is (for distinct real eigenvalues) the largest eigenvalue, or (for a complex pair) the common real part, or (for a defective repeated eigenvalue −λ) −λ/2. The extra halving in the defective case is the price of the shear: for A=(−λ01−λ) the linear flow itself decays only like te−λt, so no estimate of the form ∣u(t)∣≤C∣u(0)∣e−λt can hold, and a rate strictly smaller than λ must be accepted. Returning to the original coordinates multiplies the constant by the condition number of the change of basis but not the rate.
Remark(Where hyperbolicity is used, and what fails without it)AdaptedAdapted from Strogatz (1994), §6.3, Example 6.3.2, pp. 153-154; Hirsch, Smale & Devaney (2013), §8.1, pp. 164-165
Hyperbolicity enters the proof above as the strict inequality Au⋅u≤−μr2 with μ>0 in Step 1: it is the gap between the spectrum and the imaginary axis that gives the nonlinear remainder no room to act. When the gap closes, the remainder decides, and the linearization tells us nothing.
A centre that is not a centre. For the family studied by Strogatz (Strogatz (1994), Example 6.3.2, p. 153) and by Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), p. 164),
x˙=−y+ax(x2+y2),y˙=x+ay(x2+y2),
the Jacobian at the origin is (01−10) for every value of a: a centre, with T=0 and D=1. The nonlinear system is transparent in polar coordinates. From rr˙=xx˙+yy˙ the cubic terms give rr˙=a(x2+y2)2, so r˙=ar3; and from θ˙=(xy˙−yx˙)/r2 the cubic terms cancel and θ˙=1. Thus the origin is a stable spiral for a<0, an unstable spiral for a>0, and a genuine centre only for a=0. Arbitrarily small nonlinear terms change the answer, and even its sign.
A zero eigenvalue. For x˙=x2, y˙=−y the Jacobian at the origin is diag(0,−1), so the equilibrium is not hyperbolic. The linearized system x˙=0, y˙=−y has an entire line of equilibria, while the nonlinear system has only the origin: x˙=x2 integrates to x(t)=x0/(1−x0t), which blows up in finite time for x0>0 and decays to 0 for x0<0. No homeomorphism can match a line of equilibria to a single one, so the two flows are not conjugate near the origin. This is the example of Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), p. 165).
What survives at a hyperbolic saddle. The straight stable and unstable lines of Classification of linear planar systems do not survive the nonlinear terms, but curves do: the Stable Curve Theorem of Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), p. 169) supplies, near a hyperbolic saddle, a smooth curve of initial conditions whose solutions tend to the equilibrium as t→∞, tangent at the equilibrium to the stable eigendirection, with a mirror statement in backward time; all other nearby solutions leave a fixed small disc.
Method(Linearizing and classifying the equilibria of a planar system)Original
For x˙=f(x,y), y˙=g(x,y) with f,g of class C1:
Solve f=g=0 for all equilibria. Keep them exact where possible; the classification below uses their coordinates.
Differentiate once to get J(x,y)=(fxgxfygy) as a matrix of functions, then evaluate it at each equilibrium in turn.
Check hyperbolicity: the equilibrium is hyperbolic unless D=0, or T=0 with D>0. If it is, Hartman-Grobman licenses the picture of step 3 as the local picture of the nonlinear system, and the eigenvectors give the tangent directions of the stable and unstable curves, or of the slow and fast approach at a node.
If it is not hyperbolic, stop and use another tool: a conserved quantity or a Lyapunov function (Lyapunov functions), a symmetry, or a change to polar coordinates as in the family above. The linearization is then genuinely silent.
Assemble the local pictures into a phase portrait by adding the nullclines f=0 and g=0, on which the field is vertical and horizontal respectively.
Example(Two equilibria of a nonlinear planar system)Original
Classify the equilibria of
x˙=y−x,y˙=x2−4.
Step 1 (equilibria).y˙=0 gives x=±2 and x˙=0 gives y=x, so the equilibria are P+=(2,2) and P−=(−2,−2).
Step 2 (Jacobian).
J(x,y)=(−12x10).
Step 3 (the saddle at P+).J(2,2)=(−1410), so T=−1, D=(−1)(0)−(1)(4)=−4<0. By case 1 of Classification of linear planar systems this is a saddle, and D=0 makes it hyperbolic. Its eigenvalues solve λ2+λ−4=0:
λ±=2−1±17,λ+=1.5615528128…,λ−=−2.5615528128…
For an eigenvector, the second row of J−λI reads 4v1−λv2=0, so v=(λ,4) and the eigendirection through P+ has slope 4/λ. Because λ+λ−=−4, that slope is −λ∓: the unstable curve leaves P+ with slope
λ+4=−λ−=21+17=2.5615528128…,
and the stable curve arrives with slope 4/λ−=−λ+=−1.5615528128…. By Hartman-Grobman these are genuine curves of the nonlinear flow, tangent to those directions at P+.
Step 4 (the spiral at P−).J(−2,−2)=(−1−410), so T=−1 and D=(−1)(0)−(1)(−4)=4>0, hence Δ=1−16=−15<0. By case 3 this is a spiral, stable because T<0, and hyperbolic because T=0. Its data are
α=2T=−21,ω=214D−T2=215=1.9364916731…
Step 5 (two numbers). One turn of the spiral takes
ω2π=154π=3.2446229408…
units of time, and over that turn the distance to P− is multiplied by
e2πα/ω=e−π/ω=e−2π/15=0.1974417904…,
a shrinkage to about 19.7% per revolution. (These describe the linearized motion exactly and the nonlinear motion to leading order near P−; Hartman-Grobman guarantees the qualitative picture, not the rates.)
Step 6 (sense of rotation). At P− take the displacement u=(1,0): then J(−2,−2)u=(−1,−4), whose second component is negative, so the spiral turns clockwise.
So the portrait has a saddle at (2,2) and a clockwise stable spiral at (−2,−2); the divergence ∂x(y−x)+∂y(x2−4)=−1 is negative everywhere, which by the criterion of Limit cycles and the Poincare-Bendixson theorem already rules out any closed orbit in the plane.
ApplicationOriginal
Linearization is how stability is decided in practice wherever a steady state can be found but the full dynamics cannot be solved: the small-signal model of a transistor amplifier, the response of a laser to a perturbation of its steady photon number, the stability of a coexistence state in an ecological model, the onset of oscillation in a feedback controller. In each case one differentiates the right-hand side once at the operating point and looks at two numbers, the trace and the determinant of the Jacobian, as in Classification of linear planar systems.
The same principle is used far beyond the finite-dimensional setting. Stakgold and Holst (Stakgold & Holst (2011), §9.5, pp. 625-626) write an evolution problem as du/dt=F(u) in a Hilbert space, set u=uˉ+h around a steady state uˉ, and obtain dh/dt=Fu(uˉ)h+r with Fu(uˉ) the Fréchet derivative and r of higher order in h — exactly the decomposition of our definition block, with the Jacobian replaced by a linear operator. They then adopt the principle of linearized stability as an explicit assumption: if all solutions of the linearized equation decay exponentially, so do solutions of the full equation for small enough initial disturbance, and if the linearized equation has an exponentially growing solution, some small disturbance grows. Their stability criterion is that the spectrum of Fu(uˉ) lie in the open left half-plane, and they are careful to add that this footing is secure for ordinary differential equations, and for certain families of partial differential equations, while elsewhere the principle remains an assumption they adopt rather than a theorem they prove. The theorem proved above is the ordinary-differential-equation case of that principle, and Hartman-Grobman is the statement that at a hyperbolic steady state the linearization captures not merely the stability but the entire local geometry. Their worked example is a continuously stirred tank reactor in which the flow rate is the control parameter, and the loss of stability of a steady state as the parameter crosses a critical value is what The Hopf bifurcation and Bifurcations in one dimension take up.
Gives Lyapunov's criterion — a positive definite function that does not increase along orbits makes an equilibrium stable, and strictly decreases makes it asymptotically stable — proves it, states LaSalle's invariance principle for the non-strict case, and turns both into a recipe for building such functions and for estimating the basin of attraction with a number.
Let x˙=F(x) with F of class C1 on an open set, let x∗ be an equilibrium, and let L:O→R be differentiable on an open set O containing x∗. The derivative of L along the flow is the function
L˙(x)=DLx(F(x))=∇L(x)⋅F(x),
so that by the chain rule dtdL(x(t))=L˙(x(t)) for every solution. Note that L˙ is computed from L and F alone: no solution of the differential equation is needed.
L is a Lyapunov function for x∗ if
(a) L(x∗)=0 and L(x)>0 for every x∈O with x=x∗ (that is, L is positive definite at x∗), and
(b) L˙≤0 on O∖{x∗}.
It is a strict Lyapunov function if in addition
(c) L˙<0 on O∖{x∗}.
A set P is positively invariant if every solution starting in P is defined and stays in P for all t≥0, and invariant if this holds for all t∈R. An entire solution is an orbit {x(t):t∈R} defined for all time, forward and backward.
Theorem(Lyapunov's stability theorem)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §9.2, the Liapunov Stability Theorem, p. 193
Suppose the equilibrium x∗ of x˙=F(x) admits a Lyapunov function L on some open set O around it. Then:
x∗ is Lyapunov stable;
if moreover L is strict, so that L˙ is negative at every point of O other than x∗, then x∗ is asymptotically stable.
What makes this usable is that both hypotheses are checked by differentiating two known functions. Neither the flow nor any solution formula appears, which is why the criterion survives in situations where the equations cannot be integrated at all.
Theorem(LaSalle's invariance principle)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §9.2, Lasalle's Invariance Principle, p. 199
Take an equilibrium x∗ of x˙=F(x) and suppose three things hold.
There is a Lyapunov function L for x∗ on an open set U around it — conditions (a) and (b) only, with no strictness demanded anywhere.
There is a set P with x∗∈P⊆U which is closed, is a neighbourhood of x∗, and is positively invariant.
Apart from the equilibrium itself, no orbit that stays in P for all time, forwards and backwards, carries a single constant value of L.
Then x∗ is asymptotically stable and every point of P is attracted to it, so that P lies inside the basin of attraction.
Two things are gained over the previous theorem. First, L˙ may vanish away from the equilibrium, as long as no complete orbit can sit inside the set where it does; this is the usual situation for a dissipative mechanical system, where the energy stops falling at every instant when the velocity passes through zero. Second, the conclusion names an explicit region of the state space rather than an unspecified small neighbourhood, and that is what turns the criterion into a numerical estimate of a basin.
Idea of the proofOriginal
A positive definite L has sublevel sets {L<c} that shrink down onto x∗, and condition (b) says that the flow can never move from a smaller sublevel set to a larger one. Stability is exactly that statement: given a ball we must stay inside, take c smaller than the minimum of L on its boundary sphere — a positive number, by positive definiteness and compactness — and the sublevel set {L<c} is a trap. For asymptotic stability the extra condition (c) forces L to keep falling, so L(x(t)) converges to some limit c∞; if the orbit accumulated anywhere but x∗, the flow would carry that accumulation point strictly below c∞, and continuity in the initial condition would carry nearby points of the orbit below c∞ too — which L(x(t))>c∞ forbids. LaSalle's principle replaces "keeps falling" by "cannot stop falling forever", which is all the argument really needs.
Write Bδ, Bδ and Sδ for the open ball, closed ball and sphere of radius δ about x∗.
Step 1 (a trapping sublevel set). Let δ>0 be small enough that Bδ⊂O. The sphere Sδ is compact and misses x∗, so by hypothesis (a) and continuity
α=x∈SδminL(x)>0.
Put Uδ={x∈Bδ:L(x)<α}. It is open, and it contains x∗ because L(x∗)=0<α.
Step 2 (no solution escapes Uδ). Let x(0)∈Uδ and suppose the solution left Bδ. Being continuous, it would first meet Sδ at some finite time t1>0, and L(x(t1))≥α by the definition of α. But t↦L(x(t)) has derivative L˙(x(t))≤0 on [0,t1] by hypothesis (b) — the orbit lies in Bδ⊂O there — so L(x(t1))≤L(x(0))<α, a contradiction. Hence the solution stays in Bδ; in particular it stays in a compact subset of the domain and so is defined for all t≥0, and since L keeps decreasing it stays in Uδ.
Step 3 (stability). Given any neighbourhood U of x∗, choose δ with Bδ⊂U∩O and run Steps 1-2: the open neighbourhood Uδ⊆U has the property that every solution starting in it remains in Bδ⊆U for all t≥0. That is Lyapunov stability.
Step 4 (the limit of L). Now assume (c) and take x(0)∈Uδ∖{x∗}. By uniqueness of solutions (Existence and uniqueness of solutions) the orbit never reaches x∗ in finite time, so L˙(x(t))<0 for all t≥0 and t↦L(x(t)) is strictly decreasing. It is bounded below by 0, so it has a limit
c∞=t→∞limL(x(t))≥0,andL(x(t))>c∞ for every t≥0.
Step 5 (the only accumulation point is x∗). The orbit lies in the compact set Bδ, so it has at least one accumulation point: there are tn→∞ with x(tn)→z0∈Bδ. By continuity L(z0)=c∞. Suppose z0=x∗. Let z(s) be the solution with z(0)=z0 and fix s>0 small enough that z(s) is defined; by (c), L(z(s))<L(z0)=c∞. The map y0↦y(s) is continuous in the initial value by Dependence on initial data and Gronwall's inequality, and L is continuous, so there is a neighbourhood W of z0 such that every solution with y(0)∈W satisfies L(y(s))<c∞. Taking y(0)=x(tn) for n large enough that x(tn)∈W gives
L(x(tn+s))<c∞,
contradicting Step 4. Hence z0=x∗.
Step 6 (convergence). If x(t)→x∗, there would be ε>0 and tn→∞ with ∣x(tn)−x∗∣≥ε; compactness of Bδ would give a subsequence converging to some z0 with ∣z0−x∗∣≥ε, contradicting Step 5. So x(t)→x∗ for every initial point of Uδ, and with Step 3 this is asymptotic stability. □
Positive definiteness (a) is used twice: to make α>0 in Step 1, and to bound L(x(t)) below in Step 4. Without it the sublevel sets need not shrink onto x∗ and nothing is trapped. For L(x,y)=x2−y2 and the flow x˙=−x, y˙=−y one has L˙=−2L, which is negative on part of every neighbourhood and positive on the rest, and L says nothing at all.
L˙≤0 (b) is used only in Step 2, and only on the ball actually reached. This is why the theorem is local: L need not be defined, let alone a Lyapunov function, outside O.
The strict inequality (c) is used in Steps 4 and 5 and nowhere else. It cannot simply be dropped: the harmonic oscillator x˙=y, y˙=−x has L=x2+y2 with L˙≡0, so (a) and (b) hold, the origin is stable — and it is not asymptotically stable, since every orbit is a circle. What LaSalle adds is the observation that this example fails its hypothesis: the circles are entire solutions in P∖{x∗} on which L is constant.
Two sources, two forms.Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), p. 193) state the criterion locally, on an open set O containing the equilibrium, and conclude stability or asymptotic stability. Strogatz (Strogatz (1994), p. 201) states it globally: for a continuously differentiable V with V(x)>0 for all x=x∗, V(x∗)=0 and V˙<0 for all x=x∗, he concludes that x∗ is globally asymptotically stable and that the system has no closed orbits, referring to Jordan and Smith for the proof. We follow the local version and prove it. Where we want a global conclusion (Example 2 below) we obtain it from LaSalle applied to the compact sublevel sets of a Lyapunov function whose sublevel sets are bounded, which keeps the argument inside what is proved here.
No recipe exists. Both books say so: Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), p. 193) describe finding a Lyapunov function as a matter of ingenuity or trial and error, and Strogatz (Strogatz (1994), p. 201) that there is no systematic construction, though sums of squares often serve. The method below collects the cases in which one can be built rather than guessed.
Try the energy. For a mechanical or electrical model, take L to be the total energy measured from the equilibrium. Then L˙ is minus the dissipated power, automatically ≤0, and usually not strict — which is what LaSalle is for.
Quadratic ansatz with undetermined coefficients. Put L=ax2+2bxy+cy2, positive definite exactly when a>0 and ac>b2 (Quadratic forms and the Rayleigh quotient). Compute L˙ and choose the coefficients so that the indefinite cross terms cancel, leaving a negative definite leading part; this is the trick of Strogatz's Example 7.2.3 (Strogatz (1994), p. 201).
Borrow one from the linearization. If A=DF(x∗) has all eigenvalues with negative real part, pick any symmetric positive definite Q and set
P=∫0∞eATtQeAtdt,
which converges because ∥eAt∥ decays exponentially. Differentiating the integrand gives dtd(eATtQeAt)=eATt(ATQ+QA)eAt, so by the fundamental theorem of calculus,
ATP+PA=∫0∞dtd(eATtQeAt)dt=[eATtQeAt]0∞=−Q.P is symmetric and positive definite, and L(u)=uTPu satisfies L˙=−uTQu along the linearized flow, hence L˙<0 near x∗ for the nonlinear flow as well. This always works at a hyperbolic sink; its value is that it produces an explicit L with which to estimate a basin.
Find the region where L˙<0. Bound the higher-order terms of L˙ in terms of r=∣x−x∗∣, and solve for the radius at which the bound changes sign.
Convert it into a basin estimate. Take the largest sublevel set {L≤c} that is contained in the region of step 4 and is bounded. It is closed, positively invariant (because L cannot increase) and, if no entire orbit inside it other than the equilibrium keeps L constant, LaSalle puts all of it in the basin of attraction.
If L˙≤0 but not <0, identify the set Z={L˙=0} and ask which orbits can stay in Z forever. Substituting the constraint L˙≡0 into the equations usually forces the orbit to be the equilibrium, and then LaSalle applies.
Example(A basin of attraction, with a number)Original
Consider
x˙=−2x+y2,y˙=−y+x2.
Step 1 (equilibria and what linearization gives). Setting both right-hand sides to zero gives y=x2 and 2x=y2=x4, so x(x3−2)=0: the equilibria are the origin and
q=(21/3,22/3)=(1.2599210499…,1.5874010520…).
The Jacobian is J(x,y)=(−22x2y−1), and J(0,0)=diag(−2,−1). By Linearization and the Hartman-Grobman theorem the origin is a hyperbolic sink, hence asymptotically stable — but linearization says nothing about how far the basin extends. That is what the Lyapunov function will supply.
Step 2 (the function). Take the simplest positive definite choice, L(x,y)=x2+y2=r2. Then
L˙=2x(−2x+y2)+2y(−y+x2)=−4x2−2y2+2xy(x+y).
Step 3 (bound the two parts). The quadratic part obeys −4x2−2y2≤−2(x2+y2)=−2r2. For the cubic part, 2∣x∣∣y∣≤x2+y2=r2 (arithmetic-geometric mean) and ∣x∣+∣y∣≤2r (Cauchy-Schwarz), so ∣2xy(x+y)∣≤2r3. Therefore
L˙≤−2r2+2r3=r2(2r−2)<0for 0<r<22=2.
Step 4 (conclusion). On the open disc O={r<2}, L is positive definite and L˙<0 off the origin: a strict Lyapunov function, so the origin is asymptotically stable. For the basin, fix ρ<2 and let P={r≤ρ}. It is closed and a neighbourhood of the origin; it is positively invariant, because L=r2 cannot increase along orbits in O; and L is strictly decreasing off the origin, so no entire orbit in P∖{0} keeps it constant. LaSalle therefore puts P inside the basin of attraction. Letting ρ↑2:
{x2+y2<2}⊆basin of attraction of the origin,2=1.4142135624…
Step 5 (how conservative is this?). Two computations bound the error. First, the basin cannot contain the other equilibrium q, which lies at distance
∣q∣=22/3+24/3=2.0266334527…
from the origin, so no estimate of this kind can exceed 2.0267. Second, the estimate 2 came from a crude bound on the cubic term; putting x=rcosθ, y=rsinθ gives exactly
and L˙<0 on the punctured disc of radius ρ exactly when ρ≤min{−a(θ)/b(θ):b(θ)>0}. Minimising that ratio numerically over θ gives
ρmax=2.0237302380…,attained at θ=0.9234330…rad=52.909∘.
So this Lyapunov function, used to the full, reaches to within 0.15% of the absolute ceiling set by q, while the one-line bound of Step 3 reaches 70% of the way. Both numbers were computed here, the first in closed form and the second by minimising the ratio over a fine grid in θ and refining.
Example(LaSalle where the linearization is silent)Original
Consider the nonlinearly restored, linearly damped oscillator
x˙=y,y˙=−x3−2y.
Step 1 (linearization fails). The only equilibrium is the origin, since y=0 and x3=0. Its Jacobian is (001−2), with T=−2 and D=0: an eigenvalue is zero, so the equilibrium is not hyperbolic and Linearization and the Hartman-Grobman theorem gives no verdict.
Step 2 (energy). The mechanical energy for the restoring force −x3 is
L(x,y)=4x4+2y2,
positive definite on the whole plane, with
L˙=x3y+y(−x3−2y)=−2y2≤0.
So L is a Lyapunov function and the origin is stable — but L˙ vanishes on the entire line y=0, so it is not strict and Lyapunov's theorem stops here.
Step 3 (LaSalle). Fix c>0 and let Pc={L≤c}. It is closed and bounded (because L(x,y)→∞ as ∣(x,y)∣→∞), it contains the origin in its interior, and it is positively invariant since L˙≤0. Suppose an entire solution in Pc kept L constant. Then L˙≡0 along it, so y≡0; differentiating, y˙≡0, and the second equation gives −x3−0≡0, so x≡0. The only such solution is the equilibrium. LaSalle therefore gives asymptotic stability with Pc inside the basin — for every c, so the origin is globally asymptotically stable.
Step 4 (a number). Start at (x,y)=(2,0), where L=24/4=4. Since L never increases, for all t≥0
Numerical integration of this orbit to t=60 (Runge-Kutta, tolerance 10−10) confirms the bounds and gives an actual maximum speed maxt∣y(t)∣=1.663, with L down to 3.8×10−5 at t=60 and x(60)=0.111. The decay is slow — algebraic, not exponential — because the equilibrium is not hyperbolic: near the origin the damping dominates and x˙≈−x3/2, which gives x(t)∼t−1/2. Asymptotic stability says nothing about the rate, and here there is no exponential one.
ApplicationOriginal
The oldest and still the largest use of Lyapunov functions is transient stability in power engineering. A synchronous generator tied to a grid obeys the swing equation Mδ¨+Bδ˙=Pm−Pmaxsinδ for the rotor angle δ, which is the damped pendulum that Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), §9.2, pp. 194-195, 199) analyse with E=21v2+1−cosθ, E˙=−bv2, and then with LaSalle to bound the basin. The engineering question after a fault — will the machine return to synchronism, or slip a pole? — is exactly the question of whether the post-fault state lies in the basin of attraction of the stable equilibrium, and the classical equal-area criterion is the statement that the energy at the moment the fault clears must not exceed the energy of the saddle: a sublevel set of a Lyapunov function, computed as in Example 1 above. Whole-network versions of the same construction, under the name transient energy function methods, let an operator certify stability without integrating the equations of a system with thousands of states.
The same pattern recurs wherever the dynamics cannot be solved but a quantity is known to be dissipated. In chemical kinetics the free energy, or the relative entropy against the equilibrium composition, decreases along reaction trajectories and so certifies that a complex network relaxes to its equilibrium. In control engineering, Lyapunov's criterion is run backwards: one designs a feedback law to make a chosen positive definite function decrease, which is how most nonlinear controllers and adaptive schemes are built and proved. And in dissipative mechanics the criterion explains why linearization can be dispensed with at all: the energy is there before any coordinates are chosen, it does not require the equilibrium to be hyperbolic, and, as Example 2 shows, it answers cases where the Jacobian has nothing to say.
Gives two ways to certify that a centre predicted by linearization survives the nonlinear terms - a conserved quantity with a local minimum, or a time-reversal symmetry - together with the quadrature that turns the conserved quantity into periods, and shows how the level set through a saddle produces homoclinic and heteroclinic connections.
Definition(Conserved quantities, conservative systems, reversibility)AdaptedAdapted from Strogatz (1994), §6.5, pp. 160-161, and §6.6, pp. 164, 166-167
Throughout, f=(f1,f2) has continuous first partial derivatives on an open set U⊆R2, and the system under study is
x˙=f1(x,y),y˙=f2(x,y).
By the existence and uniqueness theorem each point of U lies on exactly one maximal solution, so the phase portrait is a partition of U into orbits.
A conserved quantity for the system is a continuous real-valued function E on U that is constant along every solution and is not constant on any nonempty open subset of U. A system that has one is conservative. The second clause is what gives the notion content: without it the function E≡0 would be conserved by every system whatsoever.
When E is differentiable, "constant along every solution" is the single equation
∇E⋅f=∂x∂Ef1+∂y∂Ef2=0on U,
by the chain rule, so a conserved quantity is a function whose gradient is everywhere orthogonal to the field. Each orbit therefore lies inside one level set of E.
The prototype is a particle of mass m on a line under a force F depending on position alone, mx¨=F(x). Writing F=−V′ for a potential V and putting y=x˙ gives the planar system x˙=y, my˙=−V′(x), and
E(x,y)=21my2+V(x)
satisfies ∇E⋅f=my⋅(−V′(x)/m)+V′(x)y=0. So E is conserved, and it is nonconstant on every open set as long as V is nonconstant on every interval.
A reversal is a map R:U→U with R∘R=id (an involution). The system is reversible with respect to R if t↦R(φ(−t)) is a solution whenever t↦φ(t) is; equivalently, by differentiating,
f(R(p))=−DR(p)f(p)(p∈U).
The case used most often is the reflection R(x,y)=(x,−y) in the x-axis, for which the condition reads: f1 is odd in y and f2 is even in y. Every system x˙=y, y˙=F(x)/m coming from mx¨=F(x) is reversible in this sense, which is the statement that a film of the motion run backwards still shows a possible motion.
An orbit γ is homoclinic to a fixed point p if the solution along it tends to p both as t→+∞ and as t→−∞, and heteroclinic from p to q if it tends to p as t→−∞ and to q as t→+∞ with p=q. A heteroclinic orbit joining two saddles is also called a saddle connection. A homoclinic orbit is not a periodic orbit: it never arrives at p, since p is a fixed point and solutions through a point are unique.
Method(A conservative system has no attracting fixed point)AdaptedAdapted from Strogatz (1994), §6.5, Example 6.5.1, p. 160
Suppose the system has a conserved quantity E. Then no fixed point p is attracting; that is, there is no neighbourhood W of p such that every solution starting in W tends to p as t→∞.
Method(Nonlinear centres in a conservative system)AdaptedAdapted from Strogatz (1994), §6.5, Thm 6.5.1, p. 163
Let f be continuously differentiable on U⊆R2, let E be a conserved quantity for x˙=f(x), and let p∈U be an isolated fixed point, meaning that some neighbourhood of p contains no other fixed point. If p is a local minimum of E, then every trajectory sufficiently close to p is closed, so p is a centre.
The same conclusion holds at a local maximum of E, since −E is also a conserved quantity and turns maxima into minima.
What is proved below. Strogatz states the theorem for a conserved quantity E that is merely continuous, and his argument rests on the assertion that the level sets of E near a local extremum are closed curves, which he states and does not prove. The proof given below supplies that step, and therefore proves the theorem in full, under the extra hypothesis that E is twice continuously differentiable near p with positive definite Hessian D2E(p) (which already forces p to be an isolated local minimum). Everything in this unit's examples meets that hypothesis.
Method(Nonlinear centres in a reversible system)AdaptedAdapted from Strogatz (1994), §6.6, Thm 6.6.1, p. 164
Let f be continuously differentiable near the origin, with f(0)=0, and suppose that the system x˙=f(x) is reversible with respect to the reflection R(x,y)=(x,−y). If the linearization Df(0) has a pair of purely imaginary nonzero eigenvalues, so that the origin is a linear centre in the sense of the trace-determinant classification, then every trajectory sufficiently close to the origin is closed, and the origin is a centre of the nonlinear system as well.
Idea of the proofOriginal
All three results turn on the same observation: a conserved quantity or a symmetry forbids the kind of one-way drift that destroys a centre.
If E is conserved, every trajectory is trapped inside a level set of E. Near a local minimum those level sets are small loops around p, and a loop carrying no fixed point has nowhere for a trajectory to stop, so the trajectory must run all the way round it. The same trapping rules out attraction: if every nearby trajectory converged to p, then E would take the value E(p) on a whole open set, which a conserved quantity is forbidden to do.
Reversibility works differently and uses no conserved quantity at all. It says that reflecting a trajectory in the x-axis and reversing its arrow produces another trajectory. A linear centre makes the angle around the origin increase steadily, so a trajectory leaving the positive x-axis arrives on the negative x-axis. Its mirror image is a second trajectory with the same two endpoints and the opposite arrow; uniqueness of solutions forces the two to be one and the same orbit, which is therefore closed.
Proof(No attracting fixed point)AdaptedAdapted from Strogatz (1994), §6.5, Example 6.5.1, p. 160
Suppose p were attracting, with W a neighbourhood of p all of whose solutions tend to p. Fix q∈W and let φ be the solution with φ(0)=q.
E is constant along φ, by the definition of a conserved quantity; call that constant E(q).
The limit identifies the constant. Since φ(t)→p and E is continuous, E(q)=limt→∞E(φ(t))=E(p). Continuity of E is used here and nowhere else.
Contradiction. The point q was an arbitrary point of W, so E≡E(p) on W. A neighbourhood is a nonempty open set, and a conserved quantity is by definition nonconstant on every nonempty open set. □
Proof(Nonlinear centres in a conservative system)AdaptedAdapted from Strogatz (1994), §6.5, proof idea of Thm 6.5.1, p. 163
Write c=E(p) and choose r>0 so small that the closed disc Br(p) lies in U, contains no fixed point except p (possible because p is isolated) and satisfies E≥c on it (possible because p is a local minimum).
Trajectories are confined to level sets.E is constant along solutions, so the orbit through a point q lies inside Lq={x:E(x)=E(q)}.
Near p the relevant level sets are loops. We use the following fact, which Strogatz states and does not prove: for all small ε>0 the connected component of {E=c+ε} that meets every neighbourhood of p is a simple closed curve encircling p. The remark below proves exactly this when E is twice continuously differentiable with a positive definite Hessian at p, together with the extra information that ∇E=0 on the loop; that is the case in every example in this unit, and we use that version from here on.
The loop carries no fixed point. By the choice of r the only fixed point in the disc is p, and p is not on the loop because E(p)=c=c+ε. This is the step that needs p to be isolated: on a level set that passed through a second fixed point, a trajectory could approach that point and never return.
A trajectory in the loop is the whole loop. Let γ be the loop, let q∈γ and let φ be the solution with φ(0)=q; its orbit O lies in γ by step 1. Since ∇E=0 along γ, the implicit function theorem makes γ a one-dimensional C1 submanifold of the plane. In a chart of γ near φ(t) the curve φ has derivative φ′(t)=f(φ(t))=0, tangent to γ because γ is a level set of E and ∇E⋅f=0; so by the inverse function theorem φ carries a neighbourhood of t homeomorphically onto a relatively open subarc of γ. Hence O is relatively open in γ, and so is every other orbit lying in γ, which makes the union γ∖O relatively open too. As γ is connected and O=∅, we get O=γ. (The same local picture is the flow box used in the next unit.)
Hence the orbit is periodic. The solution t↦φ(t) through q is a surjection R→γ that is a local homeomorphism by step 4. If it were injective it would be a continuous bijection of R onto the compact set γ, and its inverse would be continuous because the map is a local homeomorphism; that is impossible, since R is not compact and compactness is preserved by homeomorphisms. So φ(t1)=φ(t2) for some t1<t2, and uniqueness of solutions makes φ periodic with period t2−t1.
Every point of a small enough punctured neighbourhood of p lies on such a loop, so every trajectory sufficiently close to p is closed. □
Remark(The closed-contour step, proved)Original
Suppose in addition that E is twice continuously differentiable near p and that its Hessian H=D2E(p) is positive definite; then p is automatically an isolated local minimum of E. Write x=p+ρu(θ) with u(θ)=(cosθ,sinθ) and g(ρ,θ)=E(p+ρu(θ))−c.
Because ∇E(p)=0 and E∈C2, Taylor's theorem in several variables gives ∇E(x)=H(x−p)+o(∣x−p∣), hence
∂ρ∂g(ρ,θ)=∇E(p+ρu)⋅u=ρ(Hu⋅u)+o(ρ)≥ρ(λmin+o(1)),
with λmin>0 the smallest eigenvalue of H. So there is ρ0>0 with ∂g/∂ρ>0 for 0<ρ≤ρ0 and all θ; in particular ∇E=0 on the punctured disc. For each θ the function ρ↦g(ρ,θ) is therefore strictly increasing from 0 to g(ρ0,θ)≥21λminρ02=:2ε0>0, so for every ε∈(0,ε0] there is exactly one radius ρ(θ) in (0,ρ0] with g=ε. The implicit function theorem (Implicit and inverse function theorems), applicable because ∂g/∂ρ=0, makes θ↦ρ(θ) continuously differentiable and 2π-periodic. Its graph is a simple closed curve encircling p, and it is the whole of {E=c+ε} inside the disc of radius ρ0. □
Proof(Nonlinear centres in a reversible system)AdaptedAdapted from Strogatz (1994), §6.6, proof idea of Thm 6.6.1, pp. 164-165
Strogatz's route is to reflect a trajectory in the x-axis and glue it to the original; he notes that the step in which the trajectory is said to swirl round the origin is not made rigorous there. We follow the route and supply that step from the hypothesis of a linear centre.
The angle increases. Write A=Df(0)=(acbd). In polar coordinates r2θ˙=xy˙−yx˙, and expanding f to first order,
xy˙−yx˙=Q(x,y)+O(r3),Q(x,y)=cx2+(d−a)xy−by2.
A linear centre has trA=a+d=0 and detA=ad−bc>0, so d=−a and bc<−a2. The discriminant of the quadratic form Q is then (d−a)2+4bc=4a2+4bc<0, so Q is definite: there are μ>0 and a sign ε=±1 with εQ(x,y)≥μr2 everywhere. Hence there is r0>0 with
εθ˙=r2εQ+O(r)≥21μ(0<r≤r0),
so the angle advances monotonically near the origin, counterclockwise if ε=+1 and clockwise if ε=−1.
This is the only place the linear-centre hypothesis is used. With trA=0 the discriminant works out to −4detA, so definiteness of Q is exactly detA>0: the argument fails at a linear saddle, where detA<0, and the calculation of the discriminant fails outright at a linear spiral, where trA=0.
Staying small. Since f is C1 and f(0)=0, there is K with ∣f(x)∣≤K∣x∣ for ∣x∣≤r0, so ∣r˙∣≤Kr and r(t)≤r(0)eK∣t∣ as long as r≤r0. Put Tmax=2π/μ. Since ∣θ˙∣≥μ/2 while r≤r0, a time Tmax is enough for θ to move by π, and a time 2Tmax is enough for it to move by 2π. Choose
r1=r0e−KTmax,r2=r1e−2KTmax.
A trajectory starting with r≤r1 then has r≤r0 for all ∣t∣≤Tmax, and one starting with r≤r2 has r≤r1≤r0 for all ∣t∣≤2Tmax.
From the positive to the negative x-axis. Let p=(x0,0) with 0<x0≤r1, and let φ be the solution with φ(0)=p, so θ(0)=0. By steps 1 and 2 the angle increases at rate at least μ/2 while r stays below r0, so there is a first time T∈(0,Tmax] with θ(T)=επ: the solution has reached a point q=(x1,0) with x1<0.
Reflect. Let S(x,y)=(x,−y). Reversibility says that u(t):=Sφ(T−t) is again a solution. It satisfies u(0)=Sφ(T)=Sq=q (as q is on the x-axis) and u(T)=Sφ(0)=Sp=p.
Glue by uniqueness. The two solutions u(t) and φ(T+t) agree at t=0, both being q, so by uniqueness of solutions they agree for all t. At t=T this gives φ(2T)=u(T)=p=φ(0), so φ is periodic with period 2T and its orbit is closed.
Every nearby trajectory. Let z satisfy 0<∣z∣≤r2=r1e−2KTmax. Run the flow backwards from z. By step 1 the angle moves at rate at least μ/2 in modulus while r≤r0, and in backward time it moves in the direction of decreasing εθ; starting from an arbitrary angle, εθ must be decreased by at most 2π to reach the positivex-axis, which takes time at most 2Tmax. This is why the budget is 2Tmax and not Tmax: a half-turn reaches the x-axis, but possibly its negative half, and the gluing of steps 3-5 starts on the positive half. By step 2 the trajectory keeps r≤r1 throughout those 2Tmax units of backward time, so it arrives at a point (x0,0) with 0<x0≤r1. By steps 3-5 the orbit through that point is closed, and z lies on it. □
Remark(Where each hypothesis is used, and what fails without it)Original
p isolated (conservative theorem). Without it the level set through p may carry other fixed points and the nearby orbits need not close up. Strogatz records this as the second of his two remarks on the theorem.
Local extremum of E (conservative theorem). At a saddle point of E the nearby level sets are hyperbola-like and unbounded, not loops; the system x˙=y, y˙=x has the conserved quantity E=21(y2−x2), whose only critical point is a saddle of E at the origin, and the origin is a saddle of the flow with no closed orbits near it.
Linear centre (reversible theorem). Step 1 of its proof needs the quadratic form Q to be definite, which is exactly trDf(0)=0 together with detDf(0)>0. For a system reversible under a reflection, at a fixed point on the symmetry axis the trace vanishes by itself: for R(x,y)=(x,−y) the function f1 is odd in y, so f1≡0 and hence ∂f1/∂x=0 on the axis, while f2 is even in y, so ∂f2/∂y=0 there. The real content of the hypothesis is therefore detDf(0)>0, and dropping it destroys the conclusion: the second example below is reversible under R(x,y)=(−x,y) and its fixed points (0,±1) lie on the symmetry set x=0, but there detDf=−2<0 and they are saddles, with no closed orbit anywhere near them. Reversibility by itself does not even forbid attraction: that same example has an attracting node at (−1,0).
E nonconstant on open sets. Drop it and the first theorem collapses: E≡0 is then "conserved" by x˙=−x, y˙=−y, which has an attracting fixed point.
Two meanings of "conservative". We follow Strogatz, for whom a conservative system is one with a conserved quantity. Hirsch, Smale and Devaney reserve the word for the Newtonian system X′=V, V′=−gradU(X)/m coming from a force field F=−gradU, and treat the general case under the heading of Hamiltonian systems, x˙=∂H/∂y, y˙=−∂H/∂x, for which H is a first integral. Every planar Hamiltonian system is conservative in Strogatz's sense as long as H is nonconstant on open sets; the converse fails, because a conserved quantity need not generate the field. The two books agree on the mathematics and differ only in which class the word names.
Reversible is weaker than conservative. A conservative system cannot have an attracting fixed point; a reversible one can, as the second example shows. What reversibility does force is a pairing: if p is a fixed point then so is R(p), and the two have mirror-image portraits with the arrows reversed, so an attracting fixed point off the symmetry set is always matched by a repelling one.
Method(Analysing a planar system by energy or by symmetry)Original
Look for a conserved quantity. If the system is x˙=y, my˙=F(x), take V(x)=−∫xF and E=21my2+V(x). Otherwise test a candidate E against ∇E⋅f=0, or check whether the field is Hamiltonian, f=(∂H/∂y,−∂H/∂x), which happens exactly when ∂f1/∂x+∂f2/∂y=0 on a simply connected domain.
Look for a reversal. Test the reflections R(x,y)=(x,−y), (−x,y) and (−x,−y) against f(R(p))=−DR(p)f(p). For R(x,y)=(x,−y) this is the parity test: f1 odd in y, f2 even in y.
Promote the predicted centres. Two certificates are available, and either suffices.
Energy. If a conserved quantity E exists and the fixed point is an isolated local minimum or maximum of E, it is a centre. This needs no information about the linearization at all.
Symmetry. If the system is reversible, the fixed point lies on the symmetry set and the linearization there is a centre, it is a centre. The last clause cannot be dropped: the saddles (0,±1) of the second example below lie on the symmetry set of a reversible system, and nothing near them is closed.
Draw the level curves of E and put arrows on them by evaluating f at one convenient point of each; continuity fixes the rest.
Read off the connections. The level set through a saddle is the separatrix structure: a loop in it that returns to the same saddle is a homoclinic orbit, an arc joining two different saddles is a heteroclinic orbit.
Get periods by quadrature. For x˙=y, my˙=−V′(x) a closed orbit at energy E runs between consecutive roots x−<x+ of V(x)=E, and y=±2(E−V(x))/m gives
T(E)=2∫x−x+2(E−V(x))/mdx.
Near a nondegenerate minimum of V this tends to the harmonic value 2πm/V′′(x∗).
Example(A cubic potential: centre, saddle, homoclinic loop, and two periods)Original
Take m=1 and
x˙=y,y˙=x−x2.
Conserved quantity. The force is F(x)=x−x2=−V′(x) with V(x)=−21x2+31x3, so
E(x,y)=21y2−21x2+31x3
is conserved, and it is nonconstant on every open set because V is a nonconstant polynomial. The system is also reversible under (x,y)↦(x,−y): f1=y is odd in y and f2=x−x2 is even in y.
Fixed points.y=0 and x−x2=0 give (0,0) and (1,0). The Jacobian is (01−2x10).
At (0,0): eigenvalues ±1, a saddle; E(0,0)=0.
At (1,0): eigenvalues ±i, a linear centre. Now V′′(x)=−1+2x and V′′(1)=1>0, so V has a strict local minimum at x=1 and E has a strict local minimum at (1,0), with Hessian diag(1,1), which is positive definite. The fixed point is isolated. By the conservative-centre theorem, and again by the reversible-centre theorem, (1,0) is a nonlinear centre, with E(1,0)=−21+31=−61.
The homoclinic loop. The saddle sits at level E=0, whose equation is
y2=x2−32x3=x2(1−32x).
For 0<x≤23 this defines a loop leaving the origin, reaching its rightmost point (23,0) and returning; for x<0 it defines an unbounded branch. The loop carries no fixed point except the origin at its ends. On its lower half y<0, so x˙=y<0 and x decreases monotonically to a limit L≥0; if L>0 the point would converge to (L,yL) with yL=0, contradicting x˙→0, so L=0. Hence the forward orbit tends to the saddle, and by the reflection symmetry so does the backward orbit: the loop is a homoclinic orbit. It is not periodic, because reaching the origin would take infinite time.
substituting u=2x/3. (Numerical quadrature of the first integral returns 1.200000.)
Two periods. Small oscillations about the centre have period 2π/V′′(1)=2π≈6.283185. For the closed orbit at energy E=−0.1, the turning points are the two roots of 31x3−21x2+0.1=0 that bracket x=1, namely
x−=0.5670689…,x+=1.3305409…,
and the quadrature of step 7 gives
T(−0.1)=2∫x−x+2(−0.1−V(x))dx=6.7364788…
(computed to 25 digits with tanh-sinh quadrature, and confirmed to ten digits by integrating the system numerically from (x−,0) and timing the return). The orbit is 7.2% slower than the harmonic estimate, which is the anharmonicity of the cubic well.
Example(Reversible but not conservative)Original
Consider
x˙=x2+y2−1,y˙=xy.
Reversibility. Take the reflection in the y-axis, R(x,y)=(−x,y), an involution. Here f1 is even in x and f2 is odd in x, which is exactly f(R(p))=−DR(p)f(p) with DR=diag(−1,1). Concretely, if (x(t),y(t)) is a solution then so is (−x(−t),y(−t)).
Fixed points.xy=0 forces x=0 or y=0; with x2+y2=1 this gives the four points (±1,0) and (0,±1). The Jacobian is (2xy2yx):
fixed point
Jacobian
eigenvalues
type
(1,0)
diag(2,1)
2,1
unstable node
(−1,0)
diag(−2,−1)
−2,−1
stable node
(0,1)
(0120)
±2
saddle
(0,−1)
(0−1−20)
±2
saddle
Conclusion. The point (−1,0) has both eigenvalues negative, hence is attracting by the linearization theorem. A conservative system has no attracting fixed point, so this system has no conserved quantity at all: reversibility is strictly weaker than conservation. The reversal is visible in the table: R exchanges (1,0) with (−1,0) and reverses the arrows, turning the unstable node into the stable one, and fixes each of the two saddles on the symmetry axis x=0.
Application(Radial motion in a central force, and why orbits are symmetric about their apsides)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §13.2, pp. 280-281, §13.3, pp. 282-285, §13.4, p. 286, and §9.4, Proposition, p. 208; Strogatz (1994), §6.6, pp. 163-164
A particle of mass m in the plane under a central potential V(r) has, by the canonical equations, the Hamiltonian
H=2mpr2+2mr2pθ2+V(r),
in which θ does not appear, so pθ=mr2θ˙=ℓ is constant: this is the angular momentum whose conservation Hirsch, Smale and Devaney derive directly from centrality, together with the fact that the motion stays in a plane in the first place. Fixing ℓ reduces the radial motion to exactly the systems of this unit,
r˙=mpr,p˙r=−Veff′(r),Veff(r)=V(r)+2mr2ℓ2,
a one-degree-of-freedom conservative system with conserved quantity E=pr2/2m+Veff(r). Everything in this unit then reads off at once:
A circular orbit is a fixed point of the radial system, that is, a critical point r∗ of Veff. The Jacobian there is (0−Veff′′(r∗)1/m0), with eigenvalues λ satisfying λ2=−Veff′′(r∗)/m — the two cases ±iλ and ±λ that Hirsch, Smale and Devaney's proposition allows at an equilibrium of a planar Hamiltonian system, and no others. The sign of Veff′′(r∗) decides which:
Local minimum,Veff′′(r∗)>0. Then E=pr2/2m+Veff(r) has an isolated local minimum at (r∗,0) and the centre theorem applies: nearby orbits oscillate periodically in r between a pericentre and an apocentre, so the motion is bounded and radially periodic. This is the stability of nearly circular orbits.
Local maximum,Veff′′(r∗)<0. Then (r∗,0) is not an extremum of E but a saddle of it — E increases in the pr direction and decreases in the r direction — so neither the theorem nor its −E form applies. The eigenvalues ±−Veff′′(r∗)/m are real and of opposite sign, and the circular orbit is an unstable one, a saddle of the radial flow.
The radial periodTr is the quadrature of step 7 with V replaced by Veff. Over one radial period the angle advances by Δθ=∫0Trℓ/(mr2)dt, which is in general not 2π; the orbit in the plane closes exactly when Δθ is a rational multiple of 2π, and otherwise it precesses, filling an annulus densely.
The radial system is reversible under pr→−pr, t→−t, since Veff depends on r only. The gluing argument of the reversible-centre proof is the statement that the trajectory is symmetric about every apsis: the outbound arc from pericentre to apocentre is the mirror image of the inbound arc. This is the standard apsidal symmetry of central-force orbits, and it is a consequence of the symmetry alone, not of any special form of V.
Finally, a local maximum of Veff is a saddle of the radial system, and the level set through it is a separatrix in the radial phase plane, exactly like the homoclinic loop of the worked example: inside it the particle oscillates in the well, outside it escapes or falls in. This happens for attractive potentials steeper than the inverse square — for V=−k/r3, say, Veff=−k/r3+ℓ2/(2mr2) has its single critical point at r∗=3km/ℓ2, where Veff′′(r∗)=−3k/r∗5<0: a maximum, hence an unstable circular orbit with a separatrix through it.
The Kepler problem is not of that kind, and it is worth saying why. For V=−k/r and ℓ=0,
Veff(r)=−rk+2mr2ℓ2,Veff′(r)=r2k−mr3ℓ2,
which vanishes at the single radius r∗=ℓ2/(mk), where Veff′′(r∗)=k/r∗3>0. So Veff has exactly one critical point, a minimum, and the radial system has no saddle at all: every circular Kepler orbit is stable. The familiar threshold E=0 between bound ellipses and unbounded hyperbolas is not a critical level of Veff but the non-critical limiting value Veff(r)→0 as r→∞; the parabolic orbit is the borderline case E=0, on which the particle reaches infinity with zero speed.
The same picture is what makes a pendulum's phase portrait intelligible, with librations inside the separatrix and whirling motions outside it, and what underlies the magnetic-mirror confinement of charged particles, where the conserved magnetic moment turns the motion along a field line into one-dimensional motion in an effective potential well.
Attaches to every loop in the phase plane the integer number of turns the vector field makes along it, shows that this integer is the sum of the indices of the fixed points inside and equals +1 on any closed orbit, and so proves from the shape of the field alone that a periodic solution must surround fixed points of total index one.
Definition(The index of a curve and of a fixed point)AdaptedAdapted from Strogatz (1994), §6.8, pp. 174-175, 178; do Carmo (2016), §4-5, p. 283; do Carmo (2016), §5-7, Thm 1 p. 400 and Remark 1 p. 402
Let f=(u,v) have continuous first partial derivatives on an open set U⊆R2. A point p with f(p)=0 is a fixed point of x˙=f(x) (do Carmo calls it a singular point of f); it is isolated if some neighbourhood of it contains no other.
Let γ:[0,1]→U be any continuous loop, γ(0)=γ(1), with fnowhere zero on it. Then
g=∣f∘γ∣f∘γ
is a loop in the unit circle S1, and the index of γ is
Iγ=2π1[ϕ]γ,[ϕ]γ=ϕ(1)−ϕ(0),
where ϕ:[0,1]→R is a continuous determination of the angle of f along γ: a continuous function with g(t)=(cosϕ(t),sinϕ(t)) for every t. Three things have to be said, and all three come from lifting paths through the covering R→S1, s↦(coss,sins), for which see Connectedness and simple connectivity: such a ϕ exists for every continuous g; any two of them differ by a constant in 2πZ, so [ϕ]γ does not depend on which is used; and g(1)=g(0) forces [ϕ]γ∈2πZ, so that Iγ is an integer. Only continuity of γ is needed, which matters below, because the deformations in (I1) run through loops that are neither smooth nor simple.
When C is a simple closed curve in U — it does not cross itself — avoiding the fixed points, the index of C, written IC, means Iγ for a parametrisation γ traversing C once counterclockwise. In words: IC counts how many net counterclockwise turns the arrows of the field make as we walk once counterclockwise round C. It is read off the field on C alone; nothing inside C is needed.
If in addition C is piecewise C1 and f is C1, a continuous determination of ϕ is C1 on each smooth piece of C and the net change is the line integral of dϕ, which gives the computational form
IC=2π1∮Cu2+v2udv−vdu=2π1∮Cd(arctanuv).
In that case IC is the winding number n(f∘γ,0) about the origin of The winding number, which is defined there by this integral, for piecewise differentiable cycles only; the lifting definition above is what carries the notion over to merely continuous loops. Do Carmo takes the integral as the definition, on a surface, with the angle measured from a coordinate field; in the plane with Cartesian coordinates his angle is ours. Piecewise smoothness is needed for this formula and for nothing else.
If p is an isolated fixed point, the index of p is
I(p)=ICfor any simple closed curve C around p enclosing no other fixed point,
a definition justified by property (I1) below.
A closed orbit (also called a periodic orbit, or a cycle) of x˙=f(x) is the image C=φ(R) of a nonconstant periodic solutionφ: a solution defined for all t, not constant, with φ(t+S)=φ(t) for all t and some S>0. Such a φ has a least periodT>0. Indeed the set P of its periods is a subgroup of R, and it is closed, because φ is continuous and P=⋂t{S:φ(t+S)=φ(t)}; a closed subgroup of R is {0}, all of R, or TZ with T its smallest positive element. The first is excluded because S∈P, the second because φ would then be constant. Throughout, T denotes the least period; the proof of (I4) shows that φ is then injective on [0,T), so that C is a simple closed curve and IC is defined for it.
Two classical facts are taken as external input. The arguments below use the Jordan curve theorem — a simple closed plane curve C separates the plane into a bounded inside and an unbounded outside, with C as their common boundary — and the Schoenflies theorem, that the closed region D bounded by C is carried onto a closed round disc by a homeomorphism of the plane. Neither is proved in this handbook and neither is in any unit, so neither appears in the prerequisites; they are used here in the same way as Poincare's theorem in the application below, as quoted classical results. For regular curves, which is all that (I4) and the examples need, do Carmo proves a differentiable Jordan curve theorem (§5-7, Theorem 1, p. 400) and records in Remark 1, p. 402, that the inside is then homeomorphic to an open disc, attributing that proof to Stoker; the closed-disc form used in step 2 of the additivity proof is the full Schoenflies statement, which he does not give.
Theorem(Properties of the index)AdaptedAdapted from Strogatz (1994), §6.8, pp. 177-178; do Carmo (2016), §4-5, pp. 270, 284
Let f be a C1 vector field on U and let all curves below be simple closed curves in U avoiding the fixed points of f.
(I1) Deformation. If C can be deformed continuously into C′ inside U without any intermediate loop meeting a fixed point, then IC=IC′. The intermediate loops are only required to be continuous, which is why the index was defined for continuous loops.
(I2) Nothing inside. If the region enclosed by C contains no fixed point, then IC=0.
(I3) Reversing time. Replacing f by −f, which reverses every arrow, leaves every index unchanged.
(I4) A closed orbit. If C is itself a closed orbit of x˙=f(x), then IC=+1.
(I5) Additivity. If C encloses exactly the isolated fixed points p1,…,pn and no others, then
IC=k=1∑nI(pk).Theorem((I6) The index of a nondegenerate fixed point)Original
If f(p)=0 and detDf(p)=0, then p is isolated and
I(p)=sign(detDf(p)).
In particular saddles have index −1, while nodes, stars, degenerate nodes, spirals and centres all have index +1: the index sees the sign of the determinant and is blind to stability.
Neither source states the determinant criterion in this form. Strogatz computes the indices of the elementary linear types one by one (Example 6.8.4, p. 178), and do Carmo pictures three of them in Figure 4-35, p. 285 — v=(−x,−y) and v=(−y,x) with I=1, and v=(−x,y) with I=−1 — alongside a fourth field of index −3 that is not linear and is precisely a case the criterion does not reach, since such an index is possible only when detDf(p)=0. The statement above packages those computations, and it is proved below from (I1) alone.
Theorem(Closed orbits enclose index one)AdaptedAdapted from Strogatz (1994), §6.8, Thm 6.8.2, p. 180
Let C be a closed orbit of a C1 planar system. Then IC=+1. If every fixed point enclosed by C is isolated, then there are only finitely many of them and their indices sum to +1:
k∑I(pk)=+1.
Consequently:
a closed orbit always encloses at least one fixed point;
a closed orbit cannot enclose exactly one fixed point if that point is a saddle;
a system with no fixed points at all in a simply connected region has no closed orbit lying in that region.
Consequence 1 needs no isolation hypothesis. Isolation is needed for the sum, and cannot be dropped: the field f(x,y)=max(r−1,0)2(−y,x), with r=x2+y2, is C1 on the plane and has the circle r=2 as a closed orbit, yet its fixed-point set is the whole closed disc r≤1, in which no point is isolated and no index is defined.
Idea of the proofOriginal
The index is a counting device of the same family as Gauss's law and the argument principle: an integral over a boundary that reports what is inside. The angle of the field is a function defined only up to multiples of 2π, so walking round a loop can return it to a different branch; the number of branches crossed is an integer, and an integer that varies continuously with the loop cannot vary at all. That single remark gives (I1), and (I2) follows because a loop enclosing no fixed point can be shrunk to a point.
Additivity is the familiar "cut the region into simply connected pieces" argument. On a simply connected set where the field never vanishes, the angle can be defined globally and single-valuedly, so every loop there contributes nothing. Slit the inside of C along disjoint bridges running to small circles about the fixed points; what remains is simply connected, so the total circuit of its boundary contributes zero, and since each bridge is walked once each way its contributions cancel. What is left is C against the small circles, which is the identity we want.
Property (I4) is the one genuinely topological input. On a closed orbit the field is tangent to the curve, so the index of the field counts the turns of the tangent to a simple closed plane curve, and that count is ±1 — the theorem of turning tangents, due to Hopf and quoted by do Carmo. Combined with additivity it says that a periodic solution must sit over fixed points whose indices balance out to exactly one, which is a strong restriction obtained from no computation at all.
Proof((I1), (I2) and (I3))AdaptedAdapted from Strogatz (1994), §6.8, pp. 177-178
Write Z={x∈U:f(x)=0} for the fixed-point set, so f:U∖Z→R2∖{0} is continuous.
(I1). Let Γ:[0,1]×[0,1]→U∖Z be the deformation, each Γ(s,⋅) a loop, with Γ(0,⋅)=γ and Γ(1,⋅)=γ′; the hypothesis that no intermediate curve meets a fixed point is exactly what makes Γ land in U∖Z. Put
gs=∣f∘Γ(s,⋅)∣f∘Γ(s,⋅),
a loop in S1 for each s, with IΓ(s,⋅) equal to the net change of a lift of gs, divided by 2π. It is enough to show that this integer is the same for s=0 and s=1.
The map (s,t)↦gs(t) is continuous on the compact square, hence uniformly continuous, so there is a subdivision 0=s0<s1<⋯<sN=1 with ∣gsj+1(t)−gsj(t)∣<2 for every t and every j: consecutive loops are never at antipodal points of the circle. Fix j and write gsj+1(t)=Rδ(t)gsj(t), where Rδ is the rotation by δ and δ(t) is the unique such angle in (−π,π) — unique and continuous in t exactly because the two points are never antipodal, which is where the subdivision is used. If ϕ is a continuous determination of the angle for gsj, then ϕ+δ is one for gsj+1; and δ(1)=δ(0), since both loops close up and δ is determined pointwise. So the two net changes are equal, and the index is unchanged across each step of the subdivision, hence IC=IC′. (For piecewise differentiable loops this is the homotopy invariance of the winding number, The winding number; the argument just given needs nothing but continuity, and will be used again in (I5).)
(I2). Let D be the closed region bounded by C. By the Jordan curve theorem together with the Schoenflies theorem — the classical facts that a simple closed plane curve bounds, and that D is homeomorphic to a closed disc — γ is null-homotopic inside D. The hypothesis is that D∩Z=∅, so that null-homotopy takes place in U∖Z, and (I1) gives IC=Ipoint=0: a constant loop has winding number zero. See Connectedness and simple connectivity for the notion of null-homotopy used here.
(I3). If ϕ is a continuous determination of the angle of f along C, then ϕ+π is one for −f. The two differ by a constant, so they have the same net change over the circuit, and IC is the same for f and for −f. □
Proof((I5) additivity)AdaptedAdapted from Strogatz (1994), §6.8, ideas behind Thm 6.8.1, p. 179
Step 1: a single-valued angle on a simply connected set. Let W⊆U∖Z be open and simply connected. Then f:W→R2∖{0} lifts: there is a continuous Φ:W→R with f=∣f∣(cosΦ,sinΦ), because a continuous map from a simply connected space into the circle (after normalising f/∣f∣) lifts through the covering R→S1. Consequently the net change of the angle round any loop in W is zero, since the angle is given there by the single-valued function Φ.
Step 2: cut the inside of C. Let D be the closed Jordan domain bounded by C. The fixed points inside are isolated, and D is compact, so there are finitely many of them, p1,…,pn; if there were infinitely many they would accumulate, and an accumulation point of Z in D would be a non-isolated fixed point. Choose pairwise disjoint small circles σk⊂intD about the pk, each enclosing only pk, and pairwise disjoint simple arcs β1,…,βn in D, meeting no fixed point, with βk running from a point of C to a point of σk and otherwise meeting neither C nor any σj. Such arcs exist: by the Schoenflies theorem, already used in (I2), there is a homeomorphism of the plane carrying D onto a closed round disc, and in that model one simply draws n disjoint radial-looking arcs from the boundary circle to the n small discs. This is the topological input that the balloon picture in Strogatz makes plausible.
Step 3: the slit region is simply connected. Let D∗ be D with the open discs bounded by the σk removed. It is compact and carries no fixed point, but it has n holes; the slits are what remove them. Put
S=int(D∗)∖(β1∪⋯∪βn),
an open subset of U∖Z. It is connected and simply connected: in the Schoenflies model of step 2 it is an open round disc with n smaller closed discs and n arcs joining them to the boundary deleted, and such a slit disc is homeomorphic to an open disc. (It is S, not D∗, that is simply connected: deleting the closed arcs from the compact D∗ leaves a set that is neither compact nor open, which is why the simply connected set is taken open from the start.) Step 1 therefore applies with W=S.
Step 4: the boundary circuit. Let Γ be the closed loop that traverses the boundary of the slit region: C once counterclockwise, broken at the n points where the arcs meet it; each arc βk traversed inwards and later outwards; and each circle σk traversed once clockwise. The loop Γ does not lie in S, since it runs along the slits and along the deleted circles themselves. But for small η>0 the parallel circuit Γη that keeps the distance η on the inside — just inside C, down one side of each slit, just outside each σk, back up the other side — does lie in S, and Γη is homotopic to Γ through loops in D∗, hence in U∖Z, by letting η shrink to 0. By step 1 the net change of the angle of f along Γη is 0, and by the homotopy invariance proved in (I1) — which was proved for arbitrary continuous loops, and is used here for a loop that is far from simple — the net change along Γ is 0 as well.
Step 5: cancel the bridges. Each βk is traversed once in each direction, and the net change of a continuous angle along a path and along its reverse are negatives of one another, so these contributions cancel in pairs. What survives is the change along C, which is 2πIC, and the changes along the clockwise circles, which are −2πI(pk) by the definition of the index of a point. Hence
0=2πIC−2πk=1∑nI(pk),
which is (I5). Taking n=1 shows in passing that Iσ is the same for every small circle around a single isolated fixed point, so the index of a point is well defined, as claimed in the definition. □
Proof((I4) a closed orbit has index one)AdaptedAdapted from do Carmo (2016), §4-5, Theorem of Turning Tangents, p. 270; do Carmo (2016), §5-7, Thm 2 and its proof, pp. 402-403; Strogatz (1994), §6.8, property 4, p. 178
Let C be a closed orbit, the image of a nonconstant periodic solution φ, and let T>0 be its least period, which exists by the definition block. The field f is C1, hence locally Lipschitz, so solutions of x˙=f(x) are unique: two solutions that agree at one time agree at every time of their common interval (Existence and uniqueness of solutions). It is used twice in the two paragraphs that follow.
C carries no fixed point, so it is regular. If φ(t0) were a fixed point, the constant function t↦φ(t0) would be a solution agreeing with φ at t0, so by uniqueness φ would be constant, which it is not. Hence φ′(t)=f(φ(t))=0 for every t.
C is a simple closed curve. It is enough that φ be injective on [0,T), for then t↦φ(Tt) maps [0,1] onto C, is injective on [0,1) and closes up, which is what a simple closed curve is. Suppose instead φ(t1)=φ(t2) with 0≤t1<t2<T. The two functions t↦φ(t+t1) and t↦φ(t+t2) both solve x˙=f(x) and agree at t=0, so by uniqueness they agree for every t; that says t2−t1 is a period of φ, and 0<t2−t1<T contradicts the minimality of T. Minimality is what the argument needs and is used here: for a non-minimal period kT (k≥2) the solution restricted to an interval of that length is not injective, and the loop it traces runs k times round C, so it is neither simple nor the once-round parametrisation that the definition of IC asks for.
For the parametrisation in the definition of IC take γ(t)=φ(εTt), t∈[0,1], with T the least period, so that γ runs round C exactly once, and with ε=+1 if the solution runs counterclockwise round C and ε=−1 if it runs clockwise. Either way γ traverses C once counterclockwise, as the definition of IC asks; and unlike a general continuous parametrisation it is C1 and regular, so its tangent may be differentiated and discussed.
The field is a multiple of the tangent, of the sign ε. By the chain rule γ′(t)=εTφ′(εTt)=εTf(γ(t)), so
f∘γ=Tεγ′,
a positive multiple of γ′ when ε=+1 and a negative one when ε=−1.
2. Scaling and sign do not change the winding. Multiplying a nonvanishing loop pointwise by a positive continuous function does not change its angle, hence not its winding number. So in the first case IC equals the winding number of γ′. In the second case IC is the winding number of −γ′, and that too equals the winding number of γ′: if ψ is a continuous determination of the angle of γ′, then ψ+π is one of the angle of −γ′, and the two have the same net change over the circuit because they differ by a constant. (This is the same one-line argument as in (I3), applied to the tangent loop rather than to the field of the system.)
3. The turning of the tangent. What is left is the winding number of γ′, the rotation index of C. The theorem of turning tangents says it is ±1: do Carmo states it for piecewise regular plane curves (§4-5, p. 270) and attributes it to Hopf, and he proves the regular case, which is ours, as Theorem 2 of §5-7, p. 402.
Which sign. Both of do Carmo's statements say only ±1, "depending on the orientation", without saying which orientation gives which sign, so that half-line has to be supplied. His proof of Theorem 2 supplies it. He takes a line L touching the curve at a point p=γ(t0) with the whole curve on one side of L — a support line, which exists because C is compact — and reaches the value +1 in the case where turning γ′(t0) counterclockwise by π, rather than clockwise, sweeps it to −γ′(t0) across the half-plane holding the curve: that is, when the curve lies to the left of its tangent at p. For the counterclockwise traversal the inside of C lies to the left of the tangent — that is what counterclockwise means for a simple closed curve — and near the support point the inside is on the same side of L as the curve. So the sign is +1 for the counterclockwise traversal. The round circle is the check: for γ(t)=(cos2πt,sin2πt) the tangent γ′(t)=2π(−sin2πt,cos2πt) has angle 2πt+π/2, whose net change over the circuit is +2π, rotation index +1.
Hence IC=+1. □
Proof((I6) the index of a nondegenerate fixed point)Original
Let f(p)=0 and A=Df(p) with detA=0. Translate so that p=0.
Isolation. Since f is differentiable at 0, f(x)=Ax+ε(x) with ∣ε(x)∣/∣x∣→0. With m=min∣w∣=1∣Aw∣>0 (positive because A is invertible and the unit circle is compact), ∣Ax∣≥m∣x∣, so ∣f(x)∣≥(m−∣ε(x)∣/∣x∣)∣x∣>0 for 0<∣x∣≤ρ and ρ small. So 0 is isolated and I(0)=ICρ for the circle Cρ of radius ρ.
The linear part already decides. On Cρ the homotopy fs(x)=Ax+sε(x), s∈[0,1], never vanishes, by the same estimate. By homotopy invariance of the winding number, I(0) is the index of the linear field x↦Ax on Cρ.
Only the sign of the determinant survives. The set of invertible 2×2 real matrices with positive determinant is path-connected: Gram-Schmidt applied to the columns of A gives A=QR with Q a rotation and R upper triangular with positive diagonal, and both factors are joined to the identity by paths staying invertible (rotate the angle to zero; move the diagonal entries to 1 and the off-diagonal entry to 0). Along such a path As the fields x↦Asx never vanish on Cρ, so homotopy invariance of the winding number keeps the index constant, and every A with detA>0 has the index of the identity. For the identity field x↦x the angle of the field equals the polar angle, which increases by 2π, so that index is +1. If detA<0, the same argument joins A to the reflection diag(1,−1), whose field (x,−y) has angle −θ at the point of polar angle θ, giving index −1.
The elementary types. For a planar matrix detA=λ1λ2. A saddle has eigenvalues of opposite sign, so detA<0 and index −1. Nodes, stars and degenerate nodes have eigenvalues of one sign, and spirals and centres have conjugate eigenvalues α±iβ with detA=α2+β2>0; in all these cases the index is +1. □
Proof(Closed orbits enclose index one)AdaptedAdapted from Strogatz (1994), §6.8, Thm 6.8.2, p. 180
Let C be a closed orbit. By (I4), IC=+1. If the enclosed fixed points are isolated, there are finitely many of them, by step 2 of the additivity proof — this is the step that uses isolation, through the fact that an infinite subset of the compact Jordan domain would accumulate — and (I5) gives ∑kI(pk)=IC=+1.
If no fixed point were enclosed, (I2) would give IC=0, contradicting IC=+1. So at least one fixed point is enclosed, and this argument uses no isolation hypothesis.
A saddle is nondegenerate, hence isolated by (I6); if it were the only enclosed fixed point the sum would be −1, not +1.
If a simply connected region contains no fixed point and a closed orbit C lay in it, then the Jordan domain of C would lie in that region as well, so IC would be 0 by (I2), contradicting (I4). □
Remark(Where the hypotheses are used)Original
f=0 on C. Without it the angle ϕ is undefined at some point of C and IC does not exist. This is why an index is attached to a curve avoiding the fixed points, never to a curve through one.
C simple. Additivity, (I2) and (I4) all use the Jordan domain of C. A figure-eight has no inside in that sense, and the tangent of a figure-eight has winding number 0, not ±1, so (I4) fails for it outright.
Counterclockwise orientation. The index changes sign with the orientation of C; the convention is fixed once and for all in the definition. In (I4) the orientation of C as a curve is counterclockwise whatever direction the orbit runs, which is why step 2 of that proof needs (I3).
Isolated fixed points. Additivity is a finite sum, and the finiteness came from isolation plus compactness. A field such as x˙=0, y˙=−y, whose fixed points fill the whole x-axis, has no index at any of them: every small circle around a point of the axis meets the fixed-point set.
Nondegeneracy in (I6) only. Degenerate fixed points have indices too, and they need not be 0 or ±1. The field x˙=x2−y2, y˙=2xy, which is z↦z2 under z=x+iy, has a single fixed point at the origin of index +2; its mirror x˙=x2−y2, y˙=−2xy has index −2. These are exactly the cases detDf(0)=0 that (I6) does not reach.
The index is not stability. By (I3) a stable node and an unstable node have the same index +1, and so does a centre. Index theory constrains which fixed points a closed orbit can surround, never whether the orbit attracts.
Method(Computing an index)Original
By the integral. Parametrise C and evaluate
IC=2π1∮Cu2+v2udv−vdu.
This is the method to use numerically; the answer must come out an integer, which is a built-in check on the quadrature.
2. By counting turns. Mark eight or so points round C, compute f at each, redraw the arrows from a common origin in the order met, and count the net turns. Cheap, and enough when the field is simple.
3. By the Jacobian. At a fixed point with detDf(p)=0, use (I6): the index is +1 if detDf(p)>0 and −1 if detDf(p)<0. No integration at all.
4. By the argument principle, when the field is analytic. If u+iv=F(z) with z=x+iy and F analytic, then IC is the number of zeros of F inside C counted with multiplicity, by the argument principle. Anti-analytic fields, u+iv=F(z), give minus that count.
5. To exclude closed orbits. Compute the indices of all fixed points, then ask which subsets of them can have index sum +1 and be encircled by a single simple loop. If no subset can, there is no closed orbit. If some can, the test is inconclusive and a further argument is needed, for instance Dulac's criterion.
Example(An analytic field: three indices and one exclusion)Original
Take
x˙=x2−y2−1,y˙=2xy.
Under z=x+iy this is z˙=z2−1, so the field is analytic.
Fixed points.z2=1 gives (1,0) and (−1,0). The Jacobian is (2x2y−2y2x), which is 2I at (1,0) and −2I at (−1,0): an unstable star and a stable star. Both have detDf=4x2+4y2=4>0, so by (I6)
I((1,0))=I((−1,0))=+1.
The index of a big circle, three ways. Let C be the circle ∣z∣=2, counterclockwise.
By additivity. Both fixed points lie inside, so IC=1+1=2 by (I5).
By the argument principle.F(z)=z2−1 has two zeros inside ∣z∣=2, counted with multiplicity, so IC=2. (Rouche settles the count at once: on ∣z∣=2 we have ∣−1∣=1<4=∣z2∣, so z2−1 has as many interior zeros as z2 does, that is two.)
By the integral. With u=x2−y2−1, v=2xy and x=2cost, y=2sint, numerical quadrature of ∮C(udv−vdu)/(u2+v2) returns
12.566370614…=4π,soIC=2π4π=2.
The same quadrature on the circle of radius 21 about the origin, which encloses no fixed point, returns a value below 10−15 in absolute value — 6×10−16 at the default tolerance and 2×10−16 with the tolerance tightened to 10−14. The digits of such a number are pure cancellation and are not reproducible between runs or libraries; what is reproducible is that the result is zero to the accuracy of the double-precision adaptive quadrature used, whose own error estimate is 2×10−14, as (I2) requires. On the circle of radius 21 about (1,0) it returns 6.283185307179586 and on the one about (−1,0) it returns 6.283185307179588: both are 2π=6.283185307179586… to within one unit in the last place, that is an index of +1 at each fixed point, as (I6) requires.
What this excludes. Suppose the system had a closed orbit C′. By the index theorem its enclosed fixed points have index sum +1. The available subsets are: none (sum 0), {(1,0)} (sum 1), {(−1,0)} (sum 1), both (sum 2). So no closed orbit can encircle both fixed points, and any closed orbit would have to surround exactly one of them. Index theory alone stops here; it has reduced an infinite search to two cases.
Example(A field with no fixed points)Original
The system
x˙=1+x2+y2,y˙=x−y
has x˙≥1>0 everywhere, so it has no fixed point in the plane. By consequence 3 of the index theorem it has no closed orbit whatsoever — no computation of any index was needed, only the observation that the sum over an empty set of fixed points is 0 and not 1.
Application(Topological defects, and why a sphere cannot be combed)AdaptedAdapted from do Carmo (2016), §4-5, Poincare's theorem, p. 286
The index is the mathematics of topological defects, and the reason physicists meet it far from dynamical systems.
A superfluid, a thin superconducting film or a two-dimensional magnet in the XY universality class is described by a complex order parameter ψ=∣ψ∣eiχ, or equivalently by the planar vector field (Reψ,Imψ). A vortex is a point where ψ=0 — a fixed point of the field — and its charge is precisely the index defined above: the number of turns χ makes round a small loop enclosing it. Three statements of this unit become three physical facts:
Charge is quantised and measurable from far away.IC is an integer determined by the field on C alone, which is why the circulation of a superfluid round any loop is an integer multiple of h/m regardless of how the loop is drawn.
Charge is conserved under deformation. By (I1) and (I5) the total charge inside a loop cannot change while the field evolves continuously and no defect crosses the loop. Vortices therefore cannot appear or disappear singly: they are created and annihilated in pairs of opposite index, which is exactly what is seen in vortex-antivortex annihilation and in the Berezinskii-Kosterlitz-Thouless picture of unbinding.
Curvature forces defects. Do Carmo's theorem of Poincare adds up the indices over a compact surface: for a field that is differentiable and whose zeros are all isolated, the total comes to the Euler-Poincare characteristic of the surface. He draws the corollary that on a surface homeomorphic to a sphere, where that number is 2, no differentiable tangent field is nowhere zero. So a spherical shell of ordered material must carry defects of total index 2, and a differentiable wind field on the globe must have at least one calm point. This is the "hairy ball" statement. (It is also true for merely continuous tangent fields, but that does not follow from the quoted theorem without first approximating a continuous field by a differentiable one, which we do not do here.) It is the reason ordered phases on curved closed surfaces are never defect-free.
Two cautions. First, the theorem quoted is about vector fields; the director of a nematic liquid crystal is a line field, defined only up to sign, and its defects carry half-integer indices, so the bookkeeping is similar but the theorem above does not cover it as stated. Second, for the dynamical-systems use the point of (I4) is the opposite one: there the index is used not to count defects but to forbid oscillations, as in the worked examples.
Settles the long-run behaviour of a planar flow: a trajectory trapped in a compact region free of fixed points must be, or spiral onto, a closed orbit, which gives a constructive test for self-sustained oscillation, while gradient structure and Dulac's criterion give two tests that rule oscillation out.
Definition(Closed orbits, limit sets, limit cycles, trapping regions)AdaptedAdapted from Strogatz (1994), §7.0, pp. 196-197; Hirsch, Smale & Devaney (2013), §10.1, p. 215, and §10.6, p. 225
Let f be a vector field with continuous first partial derivatives on an open set U⊆R2, and write ϕt for the flow of x˙=f(x), so that t↦ϕt(X) is the solution through X (Existence and uniqueness of solutions).
A closed orbit (or periodic orbit, or cycle) is the image of a nonconstant periodic solution. A limit cycle is a closed orbit that is isolated: some annular neighbourhood of it contains no other closed orbit. It is stable or attracting if all nearby trajectories approach it, unstable if all recede, and half-stable if those on one side approach and those on the other recede. Linear systems never have limit cycles: if x(t) solves x˙=Ax and is periodic, so is cx(t) for every c=0, so each closed orbit sits in a one-parameter family and none is isolated. A limit cycle is a genuinely nonlinear object.
The ω-limit set of X is
ω(X)={Y:ϕtn(X)→Y for some tn→+∞}=T≥0⋂{ϕt(X):t≥T},
defined when the forward solution exists for all t≥0; the α-limit setα(X) is the same with tn→−∞. A limit set means either of these.
A set D⊆U is positively invariant if ϕt(D)⊆D for all t≥0. A trapping region is a compact, connected, positively invariant set; in practice one is produced by exhibiting a compact region with piecewise smooth boundary on which f points strictly inward, that is, f⋅n<0 at every boundary point where the outward normal n is defined.
Finally, let X0 be a point with f(X0)=0. Let L be the straight line through X0perpendicular to f(X0), so that f(X0) is by construction not tangent to L: writing w=f(X0), the line is L={Y:Y⋅w=X0⋅w} and f(X0)⋅w=∣w∣2>0. Since f is continuous, f(Y)⋅w>0 for all Y in some neighbourhood of X0. A local section at X0 is an open segment S⊆L containing X0 and short enough to lie in that neighbourhood. Trajectories therefore cross S transversally, and all in the same direction, namely the direction of increasing Y⋅w.
Theorem(Gradient systems have no closed orbits)AdaptedAdapted from Strogatz (1994), §7.2, Thm 7.2.1, p. 199
Suppose x˙=−∇V on U, with V:U→R single-valued and continuously differentiable. Then the system has no closed orbit.
Theorem(Dulac's criterion)AdaptedAdapted from Strogatz (1994), §7.2, Dulac's criterion, p. 202
Suppose that R⊆R2 is simply connected, and that on R
the field f of the system x˙=f(x) has continuous first partial derivatives, and
some real-valued weightg, also with continuous first partial derivatives, makes the divergence of the weighted field take one strict sign: either ∇⋅(gf)>0 at every point of R, or ∇⋅(gf)<0 at every point of R, where
∇⋅(gf)=∂x∂(gf1)+∂y∂(gf2).
Then no periodic solution of the system has its whole trajectory inside R.
Strogatz phrases condition 2 as the requirement that the divergence of the weighted field have one sign on R. The weight is entirely at the user's disposal, and the criterion is worth only as much as the choice of g.
Limit-set form (the version proved here, following Hirsch, Smale and Devaney). Let f have continuous first partial derivatives on an open subset of the plane, and let Λ be a nonempty, closed and bounded limit set of the system that contains no fixed point. Then Λ is a closed orbit.
Trapping-region form (the version used in practice, following Strogatz). Suppose that
R⊆R2 is compact;
f has continuous first partial derivatives on some neighbourhood of R;
no point of R is a fixed point of f;
for some X∈R the forward trajectory never leaves R: ϕt(X)∈R for every t≥0. Call that trajectory C.
Then C is itself a closed orbit, or it spirals towards one as t→∞; in either case R contains a closed orbit.
The second form follows from the first, as shown below. Both are statements about the plane and both are false in R3.
Idea of the proofOriginal
Two exclusion results and one existence result, with three quite different mechanisms.
For a gradient system the quantity V decreases strictly along every nonstationary trajectory, while going once round a closed orbit must return it to its starting value. Dulac's criterion replaces "V decreases" by a flux computation: a closed orbit is tangent to the field, so the flux of gf through it vanishes, and Green's theorem turns that into the vanishing of the integral of ∇⋅(gf) over the inside — impossible if the integrand never changes sign.
The Poincare-Bendixson theorem is where planarity enters. Its engine is the fact that a trajectory and a transversal segment cannot be interleaved: if a trajectory crosses a short transversal at three points in time order, the second must lie between the first and the third along the segment, because the loop formed by the trajectory arc and the piece of segment between the first two crossings is a Jordan curve the trajectory can never re-enter. That one-dimensional monotonicity is what forbids a planar trajectory from wandering forever, and it is exactly what fails in three dimensions. Given it, if Y lies in a limit set then the trajectory of Y can meet a transversal only once; but if some Z lies in the limit set of Y then the trajectory of Y returns to a transversal at Z again and again, so two of those returns must be the same point and Y is periodic.
Proof(Gradient systems)AdaptedAdapted from Strogatz (1994), §7.2, proof of Thm 7.2.1, pp. 199-200
Suppose x(t) traversed a closed orbit, with period T>0.
Round trip.x(T)=x(0), so V(x(T))−V(x(0))=0. This is where V is required to be single-valued: a multivalued potential, such as the polar angle, may return to a different branch.
The same difference as an integral. By the chain rule and the equation of motion,
Strict sign. The integrand is continuous and nonnegative, so the integral vanishes only if x˙≡0, that is only if the orbit is a single fixed point. A closed orbit is by definition the image of a nonconstant periodic solution, so the integral is strictly positive and the difference strictly negative, contradicting step 1. □
Proof(Dulac's criterion)AdaptedAdapted from Strogatz (1994), §7.2, proof of Dulac's criterion, p. 202
Suppose a closed orbit C lay entirely in R.
The inside is available.C is a simple closed curve, so by the Jordan curve theorem it bounds a region A. Because R is simply connected and C⊂R, the whole of A lies in R as well; this is the only use of simple connectivity, and the remark below shows it cannot be dropped.
with n the outward unit normal.
3. The flux vanishes. On a closed orbit the field is tangent to the curve: f(x) is the velocity of the solution tracing C. This is the step that uses the assumption that C is an orbit and not merely a curve. Hence f⋅n=0 at every point of C and the line integral is 0.
4. Contradiction. So ∬A∇⋅(gf)dA=0. But A has positive area and ∇⋅(gf) is continuous with one strict sign on R⊇A, so the integral is strictly positive or strictly negative. □
Proof(Three lemmas on local sections)AdaptedAdapted from Hirsch, Smale & Devaney (2013), §10.2, pp. 216-218, and §10.4, pp. 220-222
Throughout, S is a local section at a point X0 with f(X0)=0.
Lemma 1 (flow box and continuous arrival time). There are σ>0 and an open neighbourhood V of X0 such that every X∈V satisfies ϕt(X)∈S for exactly one t∈(−σ,σ). More generally, if ϕt0(Z0)=X0 then Z0 has a neighbourhood N and there is a continuous τ:N→R with τ(Z0)=t0 and ϕτ(Z)(Z)∈S for every Z∈N.
Proof. Write w=f(X0)=0 and G(Z,t)=ϕt(Z)⋅w. The line containing S is the level set {Y:Y⋅w=X0⋅w}. Now G is continuously differentiable in (Z,t) by smooth dependence of solutions on their data, G(Z0,t0)=X0⋅w, and
∂t∂G(Z0,t0)=f(ϕt0(Z0))⋅w=f(X0)⋅w=∣w∣2=0.
Apply the implicit function theorem: on some neighbourhood N of Z0 it produces a continuously differentiable τ with τ(Z0)=t0 and G(Z,τ(Z))≡X0⋅w, which says that ϕτ(Z)(Z) lies on the line of S; shrinking N puts it in the segment S. Taking Z0=X0, t0=0 gives the first statement, with uniqueness in (−σ,σ) because t↦G(Z,t) has derivative near ∣w∣2>0, hence is strictly increasing, for Z near X0 and ∣t∣ small. □
Lemma 2 (monotone crossings). The times at which a given orbit meets S form a discrete set, so they can be listed in increasing order. Let Y0,Y1,Y2 be the points of three successive such crossings — the orbit does not meet S at any time strictly between them — and suppose they are distinct. Then Y1 lies strictly between Y0 and Y2 in the order along the segment S. Consequently the whole sequence of successive crossings is monotone along S, and so is every subsequence of it.
Proof. Discreteness is Lemma 1: around any crossing point there is a flow box in which the orbit meets S at exactly one time of (−σ,σ).
Σ is a Jordan curve. Let A be the trajectory arc from Y0 to Y1 and let T⊆S be the closed subsegment between Y0 and Y1. By uniqueness of solutions A meets itself nowhere, and because the two crossings are successiveA∩S={Y0,Y1}, so A meets T only at the two shared endpoints. Hence Σ=A∪T is a simple closed curve; let D be the bounded open region it encloses (Jordan curve theorem). Without successiveness A could re-cross the interior of T and Σ would not be simple; the general monotonicity statement is then recovered by induction over consecutive triples.
The two leftover pieces of S. Since Σ∩S=T, the set S∖T is the union of two connected pieces: I0, beyond Y0, and I1, beyond Y1. Neither meets Σ, so each lies entirely inside D or entirely outside D.
One-way crossing. Every trajectory crosses S in the same direction, so either all of them leave D across the relative interior of T, or all of them enter. Take the first case; the second is the mirror image, obtained by exchanging the words "inside" and "outside" throughout.
The complement of D is positively invariant. A forward trajectory outside D cannot enter D through T, by the choice of case in step 3, and cannot cross A, by uniqueness of solutions.
Four rays at Y1. Take the flow-box coordinates (s,u) of Lemma 1 near Y1, in which the flow is ∂/∂s, S is {s=0}, Y1 is the origin and T is locally {s=0,u≤0}. Near Y1 the curve Σ is exactly two rays out of the origin: the incoming orbit arc {u=0,s≤0} and the segment {s=0,u≤0}. They cut a small disc about Y1 into two pieces, the quadrant {s<0,u<0} and the three-quadrant sector containing both the outgoing orbit {u=0,s>0} and the piece I1={s=0,u>0}. So the forward orbit just after Y1 and the piece I1 lie in the same component of the disc minus Σ. In the case of step 3 the forward orbit leaves D at Y1, so that component is outside D; hence I1 lies outside D too.
Four rays at Y0. The same picture at Y0, where T is locally {s=0,u≥0} and A is the outgoing ray {u=0,s>0}: the two rays cut the disc into the quadrant {s>0,u>0} and a three-quadrant sector containing I0={s=0,u<0}. In the case of step 3 trajectories cross T out of D, that is, in the direction of increasing s, so D lies locally on the side s<0 of T and the quadrant {s>0,u>0} is not in D. Therefore D is locally the other sector, and I0⊂D by step 2.
Conclusion. By step 4 the forward orbit from Y1 never meets D again, and T∪D⊆D, so Y2∈/T and Y2∈/D. Thus Y2∈S∖T=I0∪I1 with Y2∈/D⊇I0, which leaves Y2∈I1: beyond Y1 as seen from Y0. In the mirror case of step 3 steps 5 and 6 give I1⊂D and I0 outside D, while D is now positively invariant, so again Y2∈I1.
Finally, a sequence in which every consecutive triple is monotone is monotone. □
Lemma 3 (one crossing). If Y∈ω(X), then the orbit of Y meets any local section in at most one point. The same holds for Y∈α(X).
Proof. Suppose Y1=Y2 were two points of S on the orbit of Y. Both lie in ω(X), since ω(X) is invariant: if ϕtn(X)→Y then ϕtn+s(X)→ϕs(Y) by continuity of the flow. Take disjoint subsegments J1∋Y1 and J2∋Y2 of S and flow boxes around them from Lemma 1. Because Y1,Y2∈ω(X), the trajectory of X enters each flow box at arbitrarily large times, hence crosses J1 and J2 infinitely often. Choose crossings alternately, a1∈J1, b1∈J2, a2∈J1, b2∈J2,…, in increasing time. These points form a subsequence of the sequence of successive crossings of S by the orbit of X, which Lemma 2 says is monotone along S; a subsequence of a monotone sequence is monotone. But this one alternates between two disjoint subsegments, which no monotone sequence can do. □
Let Λ=ω(X) be nonempty, closed, bounded and free of fixed points (the case Λ=α(X) follows by reversing time).
Part 1: Λ contains a closed orbit. Pick Y∈Λ. The forward orbit of Y stays in the compact set Λ, because Λ is closed and invariant, so ω(Y) is nonempty: it is an intersection of a nested family of nonempty compact sets (Compactness). Moreover ω(Y)⊆ω(X)=Λ. Pick Z∈ω(Y). Then f(Z)=0, since Λcontains no fixed point, so there is a local section S at Z, with a flow box V as in Lemma 1.
Since Z∈ω(Y), there are times tn→∞ with ϕtn(Y)→Z, so ϕtn(Y)∈V for all large n; by Lemma 1 each such point reaches S after a time of modulus less than σ. So the orbit of Y meets S at two different times s<r. But Y∈ω(X), so Lemma 3 says the orbit of Y meets S at most once: therefore ϕr(Y)=ϕs(Y), and by uniqueness ϕr−s(Y)=Y with r−s>0. As Y is not a fixed point, Y lies on a closed orbit γ⊆Λ, of some least period τ>0.
Part 2: Λ=γ. It suffices to show d(ϕt(X),γ)→0 as t→∞, since then every point of ω(X) is a limit of points at distance →0 from the compact set γ, hence lies in γ; the reverse inclusion is Part 1.
Fix Y∈γ, a local section S at Y and a flow box as in Lemma 1, with parameter ϵ>0 chosen so that every point of the box reaches S within time ϵ.
A monotone sequence of crossings. Since Y∈ω(X), the trajectory of X enters the flow box at arbitrarily large times, so there are t0<t1<⋯→∞ listing the successive times at which ϕt(X)∈S near Y. Put Xn=ϕtn(X). By Lemma 2 the Xn are monotone along S, and since Y is a limit of them, Xn→Y.
The gaps are bounded. Let τ be the period of γ, so ϕτ(Y)=Y. By continuity of the flow, for Xn close enough to Y the point ϕτ(Xn) lies in the flow box, hence meets S after a further time of modulus at most ϵ. Therefore tn+1−tn≤τ+ϵ for all large n.
Uniform closeness. Let β>0. By continuous dependence on initial conditions (Dependence on initial data and Gronwall's inequality) on the compact time interval [−(τ+ϵ),τ+ϵ], there is δ>0 such that ∣W−Y∣<δ implies ∣ϕs(W)−ϕs(Y)∣<β for all ∣s∣≤τ+ϵ. Choose n0 with ∣Xn−Y∣<δ for n≥n0.
Conclusion. Let t≥tn0 and pick n≥n0 with tn≤t≤tn+1, so 0≤t−tn≤τ+ϵ by step 2. Then
d(ϕt(X),γ)≤ϕt−tn(Xn)−ϕt−tn(Y)<β,
because ϕt−tn(Y)∈γ. As β was arbitrary, d(ϕt(X),γ)→0. □
Assume (1)-(4) and let X be a point of the confined trajectory C.
The forward orbit exists for all time and is bounded, since it stays in the compact set R.
ω(X) is nonempty and compact: it is the intersection of the nested nonempty compact sets {ϕt(X):t≥T}⊆R (Compactness).
ω(X)⊆R, because R is closed and the forward orbit lies in it.
ω(X) contains no fixed point, by hypothesis (3).
By the limit-set form, ω(X) is a closed orbit γ⊆R. If X∈γ then C=γ is itself closed; otherwise Part 2 of the previous proof says d(ϕt(X),γ)→0, which is the spiralling alternative. □
Remark(Where the hypotheses are used, and what fails without them)Original
Planarity. Lemma 2 is the only place dimension two is used, through the Jordan curve theorem, and it is indispensable. In R3 a trajectory may be confined to a compact fixed-point-free region forever without approaching any closed orbit. The simplest witness is the linear flow on an invariant 2-torus: embed the torus in R3 with angle coordinates (α,β) and let α˙=1, β˙=ϱ with ϱ irrational, extending the field smoothly off the torus. The torus is compact, invariant and carries no fixed point, every orbit on it is dense in it, and none is closed. Everything else in the proof works in any dimension.
Aperiodic bounded motion in three dimensions is of course also what the Lorenz system displays, and that is the example Strogatz gives on p. 210; but it is not a witness for the sharper statement above, because the Lorenz flow has an equilibrium at the origin which lies on the attractor, so no bounded absorbing region of that flow is free of fixed points.
No fixed point in Λ. Without it the conclusion is false: the limit set may be a fixed point together with one or more homoclinic orbits returning to it, which Hirsch, Smale and Devaney illustrate on pp. 224-225, or a cycle of heteroclinic connections joining several saddles, as on their p. 215. This is why trapping regions in practice are drawn as annuli, or as a region with the fixed point punctured out.
R compact in the trapping-region form. Compactness is what makes ω(X) nonempty in step 2, and it is boundedness that is doing the work, not closedness. On the closed but unbounded region R={y≥1} the system x˙=1, y˙=0 satisfies every other hypothesis — R is closed, the field is smooth and nowhere zero, and the trajectory through (0,1) never leaves R — yet its ω-limit set is empty and the system has no closed orbit anywhere.
Simple connectivity in Dulac's criterion. It was used to put the inside of the orbit into R, and it cannot be dropped. In polar coordinates let r˙=r−1, θ˙=1, that is
x˙=x−y−rx,y˙=x+y−ry,r=x2+y2.
Then ∇⋅f=2−1/r, which is strictly positive on the annulus R={0.6<r<2} with the weight g≡1; yet the unit circle r=1 is a closed orbit lying in R. The annulus is not simply connected, and the disc inside the orbit leaves R.
Single-valued V for gradient systems. Step 1 fails otherwise. Taking V=−θ, the polar angle, the "gradient system" x˙=−∇V is the rotation field (−y,x)/r2, every orbit of which is closed.
Two definitions of "limit cycle". We follow Strogatz: a limit cycle is an isolated closed orbit. Hirsch, Smale and Devaney instead call γ a limit cycle when γ⊆ω(X) or γ⊆α(X) for some X∈/γ, that is, when some other trajectory spirals onto it. These are genuinely different conditions. Take the field given in polar coordinates by θ˙=1 together with r˙=0 for r≤1 and r˙=−(r−1)2 for r>1, which is continuously differentiable (in Cartesian form it is (−y,x)+a(r)(x,y)/r with a(r)=−max(r−1,0)2, smooth near the origin because a vanishes there, and C1 across r=1 because a(1)=a′(1)=0). Outside the unit circle r−1 decays like 1/t, so every such trajectory spirals onto r=1 and that circle is a limit cycle in the Hirsch-Smale-Devaney sense; but every circle of radius r≤1 is also a closed orbit, so r=1 is not isolated and is not a limit cycle in Strogatz's sense. Both books agree on the statement that matters here, that a nonconstant trajectory confined in a compact fixed-point-free planar region approaches a closed orbit, and the theorem proved above does not depend on the choice. Where this unit says "limit cycle" without qualification, Strogatz's meaning is intended.
The test is one-way. Dulac's criterion and the gradient test can only rule orbits out; failing to find g proves nothing. Likewise Poincare-Bendixson produces a closed orbit but says nothing about how many there are or whether any is attracting.
Method(Deciding whether a planar system oscillates)Original
To rule closed orbits out, in increasing order of effort:
No fixed points? By index theory a closed orbit must enclose fixed points of index sum +1; a system with no fixed points in a simply connected region has no closed orbit there, and one whose only enclosed fixed point is a saddle has none either.
Gradient? Check whether ∂f1/∂y=∂f2/∂x on a simply connected domain. If it holds, f=−∇V for a single-valued V, and there are no closed orbits. On a star-shaped domain a potential is given outright by
V(p)=−∫01f(sp)⋅pds,
or in practice: integrate Vx=−f1 in x, then fix the remaining function of y from Vy=−f2.
3. A Liapunov function? See Lyapunov functions; a strict one forbids closed orbits in its domain.
4. Dulac. Try g=1 first; then g=xayb on a quadrant, eax or eay. Compute ∇⋅(gf) and look for one strict sign on a simply connected region.
To prove a closed orbit exists:
Build a compact region R with piecewise smooth boundary and check f⋅n<0 on every edge. Useful boundary pieces: the coordinate axes or lines x=δ, y=δ near them; horizontal or vertical lines far out, where one component of f has a definite sign; and a line of constant x+y, which works whenever x˙+y˙ has a sign there, because dtd(x+y)=f1+f2. Sketching the nullclines f1=0 and f2=0 first shows which sign each component has in which region.
Locate the fixed points in R. If there are none, Poincare-Bendixson applies at once.
If the only fixed point p in R is a repeller, puncture R with a Lyapunov ellipse. What is needed is a small open set B∋p on whose boundary the field points strictly out of B, for then no forward trajectory can enter B and R∖B is again compact and positively invariant, now with no fixed point in it. It is not enough that every trajectory eventually leaves B, and a round disc will not do in general: a repeller can have the field pointing into a circle around it over part of that circle. Build B as a level set of a quadratic form instead. Let J=Df(p) have both eigenvalues with positive real part, so that −J has both in the left half-plane and
P=∫0∞e−JTse−Jsds
converges; P is symmetric positive definite, and integrating dsd(e−JTse−Js) from 0 to ∞ gives the Lyapunov equationJTP+PJ=I (see Lyapunov functions). Put W(X)=(X−p)TP(X−p). Along the flow, with U=X−p and f(X)=JU+O(∣U∣2),
W˙=UT(JTP+PJ)U+O(∣U∣3)=∣U∣2+O(∣U∣3),
so W˙>0 on some punctured neighbourhood of p. Choose c>0 small enough that {W≤c} lies in that neighbourhood and inside R, and take B={W<c}. On ∂B the outward normal is parallel to ∇W and f⋅∇W=W˙>0, so the field points strictly out of B. Then R∖B contains a closed orbit.
8. Report what has been proved. The theorem gives existence, not uniqueness or stability; those need a Poincare map or a Hopf-bifurcation analysis.
Example(A trapping region for an autocatalytic reaction)Original
Consider the two-species kinetic system (the Brusselator, with rate constants a=1, b=3)
x˙=1−4x+x2y,y˙=3x−x2y,x,y>0.
Step 1: the fixed point. Adding the two equations gives the identity used throughout,
dtd(x+y)=1−x.
A fixed point needs x˙+y˙=0, hence x=1, and then y˙=3−y=0 gives y=3. So (1,3) is the only fixed point. Its Jacobian is
Step 2: the trapping region. Let Q be the quadrilateral with vertices
(81,0),(81,25),(2,25),(27,0),
and check the four edges in turn.
edge
field there
inward?
x=81, 0≤y≤25
x˙=1−21+64y≥21
yes, x˙>0 points right
y=0, 81≤x≤27
y˙=3x≥83
yes, y˙>0 points up
y=25, 81≤x≤2
y˙=3x−25x2=x(3−25x)≤−641
yes, y˙<0 points down
x+y=27, 2≤x≤27
dtd(x+y)=1−x≤−1
yes, x+y decreases
The worst case on the third edge is x=81, where x(3−25x)=81(3−825)=−641<0, and on the fourth edge x≥2 throughout, so 1−x≤−1. Every edge is entered strictly, so Q is positively invariant. (A numerical sweep of 20001 points per edge returns the minimum inward components 21, 83, 641 and 1, matching the table.)
Step 3: puncture with a Lyapunov ellipse. The fixed point (1,3) lies inside Q, so hypothesis (3) fails as it stands, and it must be cut out by a set the field cannot be entering. A round disc will not serve here: at the point (1,3+ρ) of the circle of radius ρ about (1,3) the field is (ρ,−ρ), whose component along the outward radius (0,1) is −ρ<0, so on every such circle there are points where the field points inward. Step 7 of the method supplies the right set. Writing u=x−1, v=y−3, the Lyapunov equation JTP+PJ=I for J=(2−31−1) has the unique solution P=21(11776), so with
W=2(u,v)P(u,v)T=11u2+14uv+6v2
(positive definite: 11>0 and 11⋅6−72=17>0, with eigenvalues λ±=21(17±221), so λ−=1.0669…) one computes exactly, substituting u˙=2u+v+3u2+2uv+u2v and v˙=−3u−v−3u2−2uv−u2v,
W˙=2(u2+v2)+(24u3+22u2v+4uv2)+(8u3v+2u2v2).
Each monomial of degree d is at most ∣U∣d in modulus, where ∣U∣2=u2+v2, so
W˙≥∣U∣2(2−50∣U∣−10∣U∣2),
the coefficient sums being 24+22+4=50 and 8+2=10. This is positive whenever ∣U∣≤301, since 2−3050−90010=0.3222…>0. Now W≤c forces ∣U∣≤c/λ−, so c=10−3 gives ∣U∣≤0.030615<301. Take
B={11u2+14uv+6v2<10−3}.
Then W˙>0 on ∂B — a direct numerical minimisation over ∂B gives minW˙=1.139×10−4, comfortably above zero and consistent with the crude bound — so the field points strictly out of B at every boundary point. The ellipse is tiny, ∣x−1∣≤0.0188 and ∣y−3∣≤0.0255, hence well inside Q.
The set Q∗=Q∖B is therefore compact, positively invariant (nothing leaves Q across ∂Q by step 2, and nothing enters B across ∂B) and contains no fixed point. By the Poincare-Bendixson theorem, Q∗ contains a closed orbit.
Step 4: what it looks like. Integrating the system numerically from (1.1,3) with absolute and relative tolerances 10−14 and 10−12, the solution settles onto a cycle on which successive upward crossings of the line x=1 are separated by
T=7.156920(six decimal places; the five successive gaps agree to 10−12),
and the cycle occupies 0.3706≤x≤3.7518, 0.8523≤y≤4.7221 — comfortably inside Q, as it must be. The existence of the orbit was proved before any of these numbers were computed; the numbers only say where it is.
Example(Dulac rules out oscillation in a competition model)Original
Two competing populations obey
x˙=x(4−x−3y),y˙=y(3−2x−y),x,y>0.
Take R to be the open first quadrant, which is simply connected, and the weight g(x,y)=1/(xy), continuously differentiable and nonzero there. Then
gx˙=y4−x−3y,gy˙=x3−2x−y,
so
∇⋅(gf)=∂x∂(y4−x−3y)+∂y∂(x3−2x−y)=−y1−x1<0
throughout the quadrant. By Dulac's criterion there is no closed orbit in the open first quadrant. (The interior fixed point is (1,1), where the Jacobian is (−1−2−3−1) with det=−5<0: a saddle. Index theory gives the same conclusion for any orbit enclosing only that point, but Dulac's computation excludes every closed orbit in the quadrant at once.)
Example(A gradient system)Original
For
x˙=3x2−3y,y˙=−3x+3y2,
we have ∂f1/∂y=−3=∂f2/∂x, so on the plane, which is simply connected, the field is a gradient:
f=−∇V,V(x,y)=−x3+3xy−y3,
as one checks by differentiating. Hence the system has no closed orbits, although it has two fixed points, (0,0) and (1,1) (solve y=x2 and x=y2). No trapping-region argument is needed and no index is computed.
Application(Chemical clocks, and why chaos needs a third dimension)AdaptedAdapted from Strogatz (1994), §7.3, No Chaos in the Phase Plane, p. 210; Hirsch, Smale & Devaney (2013), §10.7, pp. 228-229
Self-sustained oscillation. A system that returns to the same cycle after being disturbed, with an amplitude and period set by the equations rather than by the initial data, is a limit cycle; a closed orbit that is not isolated will not do, because a perturbation moves the system permanently to a neighbouring orbit. Poincare-Bendixson is the standard way to prove that a two-variable model has such a cycle without solving it: build a trapping region from the nullclines, check that the single interior fixed point is a repeller, puncture, conclude. This is how oscillation is established in models of glycolysis, in the two-variable reductions of the Belousov-Zhabotinsky and chlorine dioxide-iodine-malonic acid reactions discussed by Hirsch, Smale and Devaney, and in two-variable caricatures of pacemaker cells. The bonus in applications is that the construction gives where the oscillation lives: in the worked example the orbit was confined to a specific quadrilateral before any integration was performed, which is the kind of bound a numerical study cannot supply by itself.
A sharp limit on two-variable models. The theorem also says what a planar model can never do. Suppose the system has only finitely many equilibria, and consider a trajectory that stays in a bounded set. Its ω-limit set is then nonempty, closed and bounded, and it is either a closed orbit, by the theorem proved above, or it contains an equilibrium; in the latter case Hirsch, Smale and Devaney record (p. 215) that such a limit set is made up of equilibria together with the orbits joining them, so the trajectory converges to a fixed point, to a homoclinic loop, or to a cycle of heteroclinic connections. The hypothesis of finitely many equilibria is needed: a planar system whose fixed-point set is a whole curve can have a bounded trajectory whose limit set is an arc of that curve, which is none of the three. With it, there is no room for aperiodic bounded motion that settles to nothing. So no two-variable deterministic model, however nonlinear, can produce the irregular, sensitively dependent time series that turbulence, chemical chaos and simple nonlinear circuits display. Three dimensions are the minimum, and the Lorenz equations are the standard demonstration that three suffice. When a planar model of a real system fails to reproduce observed irregularity, this theorem says the fault is structural: another degree of freedom, a delay, forcing or noise must be added, and no amount of tuning the two-variable model will do instead.
Describes the one way a planar fixed point can lose stability without colliding with another fixed point - a conjugate pair of eigenvalues crossing the imaginary axis - and separates the benign supercritical case, in which a small stable cycle of amplitude proportional to the square root of the distance past the threshold grows out of the fixed point, from the dangerous subcritical case, in which an unstable cycle shrinks onto the fixed point and the state jumps to a distant attractor.
Let f(⋅,μ) be a family of vector fields on an open set of the plane, depending on a real parameter μ, with f and its first partial derivatives in x continuous in (x,μ). Suppose the family has a branch of fixed points x∗(μ), continuous in μ, and write
J(μ)=Dxf(x∗(μ),μ),λ±(μ)=α(μ)±iω(μ)
for the Jacobian there and its eigenvalues. The family has a Hopf bifurcation at μ=μc when
the eigenvalues at μc are a nonzero purely imaginary pair: α(μc)=0 and ω(μc)=0, and
they cross the imaginary axis as μ passes μc: α changes sign there, the usual sufficient condition being the transversality conditionα′(μc)=0.
Condition 1 is the only alternative to a zero eigenvalue. If a planar fixed point is stable, both eigenvalues lie in Reλ<0; since they are the roots of a real quadratic they are either both real or a conjugate pair; and to destabilise the fixed point one or both must reach the imaginary axis. A real eigenvalue reaching 0 gives the saddle-node, transcritical and pitchfork bifurcations of Bifurcations in one dimension, all of which involve a collision of fixed points. Condition 1 is the other case, and in it no new fixed point appears: what is created or destroyed is a closed orbit.
The bifurcation is called
supercritical when a stable limit cycle exists on the side of μc where x∗ is unstable, and shrinks to x∗ as μ→μc from that side;
subcritical when an unstable limit cycle exists on the side where x∗ is stable, and shrinks onto x∗ as μ→μc from that side, so that past μc the state must leave the neighbourhood altogether;
degenerate when conditions 1 and 2 hold but there is no limit cycle on either side.
"Limit cycle" means an isolated closed orbit, as in Limit cycles and the Poincare-Bendixson theorem. Nothing in the definition is special to the plane: Hopf bifurcations occur in any dimension n≥2, with the remaining n−2 eigenvalues staying off the imaginary axis. Everything below is planar.
Fix ω=0 and b∈R, and consider, in polar coordinates on R2∖{0} together with the origin,
r˙=μr−r3,θ˙=ω+br2.
In Cartesian coordinates this is the polynomial system
x˙=μx−ωy−(x2+y2)(x+by),y˙=ωx+μy+(x2+y2)(bx−y),
which has the origin as a fixed point, with Jacobian there (μω−ωμ) and eigenvalues μ±iω. So α(μ)=μ, ω(μ)=ω, and conditions 1 and 2 of the definition hold at μc=0 with α′(0)=1.
Then:
For μ≤0 the origin is the only fixed point; it is asymptotically stable and attracts every trajectory, and there is no closed orbit. This includes μ=0, where the linearization is a centre.
For μ>0, provided ω+bμ=0, the origin is again the only fixed point and is an unstable spiral, the circle r=μ is a closed orbit, and every trajectory other than the origin spirals onto it. It is therefore a stable limit cycle, and the only closed orbit of the system; its amplitude is exactly μ and its period is exactly 2π/∣ω+bμ∣.
The hypothesis ω+bμ=0 cannot be dropped, and its failure is not a bookkeeping nuisance. Suppose b=0 and μ=−ω/b>0. On the circle r=μ both r˙=r(μ−r2) and θ˙=ω+br2=ω+bμ vanish, so the Cartesian field is identically zero there: the circle is a curve of fixed points, not a closed orbit, the origin is not the only fixed point, and 2π/∣ω+bμ∣ is undefined. (With ω=1, b=−1, μ=1 this is the unit circle.) For fixed b=0 this happens at the single parameter value μ=−ω/b, and only if that value is positive.
The bifurcation at μ=0 is supercritical. There is at most one exceptional μ, and none at all when b=0 or when −ω/b≤0; when there is one it sits a fixed positive distance from 0. So for all sufficiently small μ>0 the second item applies and the small stable cycle is there.
Theorem(The subcritical model)AdaptedAdapted from Strogatz (1994), §8.2, pp. 251-252
With ω=0 and b∈R as before, consider
r˙=μr+r3−r5,θ˙=ω+br2,
whose Cartesian form again has the origin as a fixed point with Jacobian eigenvalues μ±iω. Write ρ=r2. As in the previous statement, assume
ω+bρ=0for every root ρ≥0 of ρ2−ρ−μ=0;
equivalently, if b=0, that −ω/b is not such a root. A fixed point away from the origin needs r˙=0 and θ˙=0 at once, which is exactly a root ρ of that quadratic with ω+bρ=0; so under this hypothesis the origin is the only fixed point, and the circles named below are genuine closed orbits with period 2π/∣ω+bρ∣. Without it the offending circle consists entirely of fixed points, just as in the supercritical model. Then:
For −41<μ<0 there are exactly two closed orbits, the circles of radii
r±=(21±1+4μ)1/2,
the inner one (r−) unstable and the outer one (r+) stable, and the origin is asymptotically stable. Three attracting or repelling objects coexist, and which one a trajectory meets depends on where it starts: inside r− it falls to the origin, outside it winds onto r+.
At μ=0 the inner cycle has shrunk to the origin: r−=0, r+=1, and the origin is unstable (weakly, with r˙=r3−r5>0 for small r>0).
For μ>0 the origin is an unstable spiral and the only closed orbit is r+=(21+1+4μ)1/2>1, which is stable.
At μ=−41 the two cycles coalesce at r=1/2=0.7071067… and for μ<−41 there is no closed orbit at all.
The bifurcation at μ=0 is subcritical. The event at μ=−41 is not a Hopf bifurcation — the origin takes no part in it — but a saddle-node bifurcation of cycles, treated in Global bifurcations of cycles.
Theorem(A cycle is born whenever the cubic coefficient has the right sign)Original
The two models are not special. Suppose that on {0<r≤r0}, for ∣μ∣≤μ0, a planar family has the polar form
r˙=μr+ar3+R(r,θ,μ),θ˙=ω+Θ(r,θ,μ),
with the origin a fixed point for every μ, where ω=0, a<0, the functions R,Θ are continuous and 2π-periodic in θ, and for some constant C
∣R(r,θ,μ)∣≤Cr4,∣Θ(r,θ,μ)∣≤Cr(0≤r≤r0,∣μ∣≤μ0).
Put
μ1=min{μ0,2∣a∣r02,32C2∣a∣3,2C2∣a∣ω2}.
Then:
For −μ0≤μ≤0 the origin is asymptotically stable, and no closed orbit meets the disc {r≤r1}, where r1=min{r0,∣a∣/(2C)}.
For 0<μ<μ1 the compact annulus
Aμ={2∣a∣μ≤r≤∣a∣2μ}
is positively invariant and contains no fixed point, so by the Poincare-Bendixson theorem it contains a closed orbit; every such orbit has all its radii between μ/(2∣a∣) and 2μ/∣a∣.
So the amplitude of the born cycle is pinned between two multiples of μ. The argument gives existence, not uniqueness: it does not by itself show the orbit is isolated, hence a limit cycle. For a>0 the mirror-image argument (run time backwards) produces an unstable cycle for −μ1<μ<0, which is the subcritical picture.
Theorem(Two rules of thumb)AdaptedAdapted from Strogatz (1994), §8.2, Rules of Thumb, p. 251
Strogatz (Strogatz (1994), p. 251) draws two generic consequences from the model above. He is careful to call them rules of thumb and not theorems, and they are reported here with that status.
Period at onset. Read off the imaginary part of the eigenvalue pair at μ=μc. At the instant the cycle is born this is its angular frequency exactly; a little past μc it is wrong by O(μ−μc), so
T=Imλ(μc)2π+O(μ−μc).
Amplitude just past onset. No cycle of finite size appears from nowhere: the amplitude rises continuously out of zero, and near μc it is a constant times μ−μc.
He also warns against reading two features of the model as general. Its cycle is an exact circle, whereas a cycle met in practice is elliptical near onset and distorts further as μ leaves μc — so the model is a faithful picture topologically, not geometrically. And here Imλ is the constant ω, so the eigenvalue pair slides along horizontal lines; a typical pair follows a curved path and meets the imaginary axis at an angle.
Idea of the proofOriginal
Everything turns on one observation: when the eigenvalues are ±iω with ω=0, the angular velocity near the fixed point is close to ω and in particular never zero, so trajectories must circulate. The radial motion is then a one-dimensional flow, r˙=μr+ar3+⋯, carried round and round.
That reduces the whole question to a bifurcation on the half-line. The linear term μr changes the stability of r=0 as μ crosses zero, and the first nonlinear term ar3 decides what the trajectory does once μr has pushed it away. If a<0 the cubic term fights the growth and the two balance at r≈μ/∣a∣: a stable zero of the radial field, which in the plane is a stable cycle. If a>0 both terms push outwards for μ>0 and nothing stops the trajectory nearby, while for μ<0 they balance at an unstable zero: an unstable cycle, which collapses onto the origin as μ↑0.
For the two models the radial equation is exact and this is a complete analysis. For a general family the higher-order terms are not zero, but they are O(r4), and at the radii of interest r∼μ that is O(μ2) against the O(μ3/2) of the terms that matter — small enough that the inward and outward inequalities on two circles survive. That is the content of the third statement: not a change of coordinates, just two circles and a trapping region.
using rcosθ=x, rsinθ=y and r2=x2+y2; the computation for y˙ is the same with y˙=r˙sinθ+rθ˙cosθ. Both right-hand sides are polynomials whose nonlinear parts are cubic, so they are differentiable at the origin and the Jacobian there is the matrix of the linear parts, (μω−ωμ). Its characteristic polynomial is (λ−μ)2+ω2, with roots μ±iω. Since θ˙=ω+br2=0 for r small, the origin is the only fixed point near it. Globally, a fixed point with r>0 would need r˙=0 and θ˙=0 at the same place, that is r2=μ and ω+br2=0, hence ω+bμ=0 — which the hypothesis of the statement excludes. So for μ>0 with ω+bμ=0, and for every μ≤0, the origin is the only fixed point in the plane.
2. The radial flow decides. Because r obeys the autonomous one-dimensional equation r˙=h(r):=r(μ−r2), the analysis of flows on the line applies to it verbatim on r≥0.
3. μ≤0. For r>0, h(r)=r(μ−r2)<0, so r(t) is strictly decreasing and bounded below by 0; its limit L≥0 satisfies h(L)=0 by continuity of h (if h(L)<0 then r would keep decreasing past L), and the only nonnegative zero is L=0. Hence every trajectory tends to the origin and no closed orbit exists, since on a closed orbit r would have to return to its initial value. At μ=0 this uses the cubic term: the linearization there is a centre and says nothing.
4. μ>0, the cycle.h(r)=r(μ−r2) vanishes on r>0 only at r=μ. The set r=μ is therefore invariant, and on it θ˙=ω+bμ, a constant that is nonzero by the hypothesis ω+bμ=0; so θ advances by 2π in time 2π/∣ω+bμ∣ and the circle is a closed orbit of exactly that period.
5. μ>0, stability and uniqueness.h>0 on (0,μ) and h<0 on (μ,∞), so r(t)→μ for every initial r>0, by the monotone argument of step 3 applied on each side. Hence every nonzero trajectory spirals onto the circle, which is consequently isolated among closed orbits — a limit cycle — and is the only closed orbit, because a closed orbit would need r to be periodic, while r is strictly monotone off r=μ. At the origin h(r)=r(μ−r2)>0 for small r>0 and the eigenvalues are μ±iω with μ>0: an unstable spiral. □
Proof(The subcritical model)AdaptedAdapted from Strogatz (1994), §8.2, pp. 251-252
Again the radial equation r˙=h(r)=r(μ+r2−r4) decides everything, and again θ˙=ω+br2 is nonzero for small r and, by the hypothesis of the statement, at each of the cycle radii as well, so every invariant circle below really is an orbit.
Where the cycles are. For r>0, h(r)=0 exactly when ρ=r2 solves ρ2−ρ−μ=0, that is ρ=21(1±1+4μ). Both roots are real iff μ≥−41, and both are positive iff in addition μ<0 (their product is −μ and their sum is 1). For μ≥0 the root with the minus sign is ≤0 and only r+ survives. This gives the four cases listed.
Their stability. At a positive zero r we have μ=r4−r2, so
h′(r)=μ+3r2−5r4=(r4−r2)+3r2−5r4=2r2(1−2r2).
Thus h′(r)>0 when r2<21 and h′(r)<0 when r2>21. Since r−2=21(1−1+4μ)<21<21(1+1+4μ)=r+2 whenever μ<0 and 1+4μ>0, the inner circle is a repelling zero of the radial field and the outer one an attracting zero: an unstable and a stable limit cycle respectively.
3. The origin. For −41<μ<0, h(r)=r(μ+r2−r4)<0 for 0<r<r− because μ+r2−r4 is negative between its roots... more directly, μ+ρ−ρ2 is a downward parabola in ρ vanishing at ρ=r±2, hence negative for 0≤ρ<r−2. So the origin attracts the disc r<r−. At μ=0, h(r)=r3(1−r2)>0 for 0<r<1, so the origin repels, although only algebraically: the eigenvalues ±iω are on the axis and the linearization again says nothing. This is the step where the destabilising sign of the cubic term is used.
4. The collision. At μ=−41 the discriminant 1+4μ vanishes, r−=r+=1/2, and h(r)=r(−41+r2−r4)=−r(r2−21)2≤0 with equality only at r=1/2: a single half-stable cycle. For μ<−41 the parabola μ+ρ−ρ2 has no real root, so h(r)<0 for all r>0 and there is no closed orbit; the origin attracts everything. □
Proof(A cycle is born whenever the cubic coefficient has the right sign)Original
where the middle inequality uses Cr3≤21∣a∣r2, which is exactly r≤∣a∣/(2C), and the last uses μ≤0. So r decreases strictly on (0,r1]; the disc {r≤r1} is positively invariant and, by the argument of step 3 of the first proof, every trajectory in it tends to the origin. A closed orbit meeting the disc would have to stay in it (positive invariance) and have r periodic, which strict decrease forbids.
2. Part 2, the outer circle. Let 0<μ<μ1 and ρ+=2μ/∣a∣. From μ≤∣a∣r02/2 we get ρ+≤r0, so the hypotheses apply there. Since ∣a∣ρ+2=2μ,
Now Cρ+3<μ/2: squaring, this is C2(2μ/∣a∣)3<μ2/4, that is μ<∣a∣3/(32C2), which is one of the four bounds defining μ1. Hence r˙≤−21μρ+<0 on the outer circle: the field points strictly inward.
3. The inner circle. Let ρ−=μ/(2∣a∣), so ∣a∣ρ−2=μ/2. Then
using ρ−<ρ+ and the bound of step 2 again. The field points strictly outward across the inner circle.
4. No fixed point in the annulus. On Aμ we have r≤ρ+, so ∣Θ∣≤Cρ+<∣ω∣ precisely when C2⋅2μ/∣a∣<ω2, that is μ<∣a∣ω2/(2C2) — the last of the four bounds. Hence θ˙=ω+Θ has the sign of ω and is nonzero, so the velocity field does not vanish anywhere on Aμ. This is the only place the hypothesis ω=0 is used, and it is indispensable: with ω=0 the annulus could contain fixed points and the next step would fail.
5. Conclusion.Aμ is compact and connected; by steps 2 and 3 the field points strictly into it on both boundary circles, so it is positively invariant; by step 4 it contains no fixed point. The Poincare-Bendixson theorem, in the trapping-region form of Limit cycles and the Poincare-Bendixson theorem, gives a closed orbit inside it. Any closed orbit lying in Aμ has all its radii in [ρ−,ρ+] by definition of Aμ; and no closed orbit can cross either boundary circle, because the field is strictly transverse there. □
Remark(Where the hypotheses are used, and what goes wrong without them)Original
ω(μc)=0. Used in step 4 of the last proof to keep the annulus free of fixed points, and implicitly everywhere else to make "the radial equation plus rotation" an honest description. With ω(μc)=0 the eigenvalue pair is 0,0 and the bifurcation is of a different family altogether.
Transversality α′(μc)=0. Not used in any proof above, because the models have α(μ)=μ by construction. It is what guarantees that the eigenvalues really do cross rather than touch the axis; without it, α(μ)=μ2 say, the fixed point is unstable on both sides and no transition occurs.
The sign of the cubic coefficient, and why linearization cannot supply it. The two models have the same linearization at the origin for every μ — eigenvalues μ±iω — and opposite behaviour. So no amount of eigenvalue computation decides whether a Hopf bifurcation is sub- or supercritical; the information lives in the cubic terms. Strogatz (Strogatz (1994), p. 253) reports that a closed-form criterion exists but is awkward enough that he sends the reader to the computer instead. The experiment has two halves, and both are needed. Push the parameter a little past μc and see what the trajectory settles onto: a cycle of small amplitude there is the supercritical signature, whereas a jump to something far away is not. Then bring the parameter back down: if the cycle shrinks away again and the fixed point is recovered at the values it had before, the first reading stands; if the two sweeps disagree over an interval, the bifurcation is subcritical. He is explicit that none of this proves anything — a numerical experiment is evidence and the numerics themselves need checking — and the worked example below is run in that spirit.
Degenerate case. Conditions 1 and 2 can hold with no cycle on either side. Strogatz's example is the pendulum with damping of either sign, x¨+μx˙+sinx=0, equivalently x˙=y, y˙=−sinx−μy. At the origin the Jacobian is (0−11−μ) with trace −μ and determinant 1, hence eigenvalues 21(−μ±μ2−4): a stable spiral for 0<μ<2, an unstable spiral for −2<μ<0, and ±i at μ=0. Both conditions hold. But at μ=0 the system is conservative, with E=21y2+1−cosx, and by Conservative and reversible systems the origin is a nonlinear centre: a continuous band of closed orbits, none of them isolated, hence no limit cycle, and none for μ=0 either. The usual signature of degeneracy is exactly this, a system that becomes conservative at the bifurcation value.
Isolation is not proved. The third statement produces a closed orbit in the annulus, not a limit cycle. Proving uniqueness needs more — a Poincare map, as in Poincare maps, or a Dulac function, or the normal-form reduction that turns the O(r4) remainder into an O(r5) one. In the two models uniqueness came free because the radial field was exactly autonomous.
The bound ∣R∣≤Cr4 is a hypothesis, not a consequence. A general smooth planar family with linear part (μω−ωμ) has r˙=μr+O(r2), with a quadratic term that depends on θ; it is only after a polynomial change of coordinates (the Poincare normal form) that the quadratic term and the θ-dependence of the cubic term are removed, leaving r˙=μr+ar3+O(r4). That reduction is not carried out here; the statement takes its output as a hypothesis and shows what follows from it.
Method(Detecting and classifying a Hopf bifurcation)AdaptedAdapted from Strogatz (1994), §8.2, pp. 248-253, and §8.3, Example 8.3.3, p. 260
Given a planar family x˙=f(x,μ):
Find the branch of fixed pointsx∗(μ). Solve f=0 symbolically if possible; by the implicit function theorem the branch is as smooth as f wherever Dxf is invertible, which it is throughout a Hopf bifurcation (the eigenvalues are ±iω=0, so the determinant is ω2>0).
Compute the trace and determinant of J(μ)=Dxf(x∗(μ),μ), not the eigenvalues. For a 2×2 real matrix the eigenvalues are 21(τ±τ2−4Δ) with τ=trJ, Δ=detJ, so
α=21τ,ω=214Δ−τ2when τ2<4Δ.
Solve τ(μ)=0 for candidate values μc, and at each check Δ(μc)>0. That is condition 1: then ω(μc)=Δ(μc)=0. (If Δ(μc)≤0 the eigenvalues are real at μc and this is not a Hopf point.) See Classification of linear planar systems for the trace-determinant picture: a Hopf bifurcation is the crossing of the positive Δ-axis.
Check transversality by differentiating: α′(μc)=21τ′(μc)=0.
Quote the period as T≈2π/ω(μc)=2π/Δ(μc), by rule of thumb 2.
Decide sub- or supercritical numerically, as in the remark: integrate at a few parameter values on the unstable side, measure the amplitude of whatever the trajectory settles onto, and test whether it goes to zero like μ−μc. Equivalently, check whether (amplitude)/μ−μc tends to a finite nonzero constant. Then sweep back and look for hysteresis.
Report the scaling, not just the threshold. It is the pair (amplitude ∝μ−μc, period →2π/Δ) that identifies the bifurcation as a Hopf event rather than one of the global bifurcations of Global bifurcations of cycles, which have quite different scalings.
Example(Enrichment destabilises a predator-prey equilibrium)Original
A predator y feeds on a prey x whose growth saturates at carrying capacity K, with a feeding rate that saturates at high prey density:
x˙=x(1−Kx)−1+xxy,y˙=y(1+xx−21),x,y>0.
The parameter is K. Everything is dimensionless; the predator's death rate has been fixed at 21.
Step 1: the coexistence equilibrium.y˙=0 with y>0 needs x/(1+x)=21, so x∗=1 for every K. Then x˙=0 gives 1−1/K=y/2, so
x∗=1,y∗=2(1−K1),
which is in the positive quadrant exactly when K>1.
Step 2: trace and determinant. Differentiating and substituting x=1, y=y∗,
Step 3: the Hopf point.τ=0 only at Kc=3, and there Δ=2/12=1/6>0. The eigenvalues at Kc are ±i1/6=±0.4082483i, purely imaginary and nonzero: condition 1 holds. For K<3 the trace is negative and the equilibrium is stable; for K>3 it is positive and the equilibrium is unstable. It is a spiral only when τ2<4Δ, which here reads (K−3)2<4K(K−1), that is 3K2+2K−9>0, that is K>31(27−1)=1.4305. So on 1<K<1.4305 the equilibrium is a stable node, and only past K=1.4305 does it spiral — comfortably below the bifurcation value, so the picture near K=3 is spiral on both sides.
Step 4: transversality.τ′(K)=2K23, so α′(3)=21τ′(3)=121>0=0. A genuine crossing.
Step 5: the predicted period. Rule of thumb 2 gives
T≈Δ(3)2π=2π6=15.390598.
Step 6: supercritical or subcritical? Simply integrating and waiting is not enough this close to the threshold: the cycle attracts at a rate of order ∣τ∣/2≈1.7×10−4 at K=3.002, so a transient of length 104 still leaves errors in the sixth decimal. Instead, locate the cycle exactly by shooting. Take the section x=1, x˙>0, which the field crosses at y<y∗, let Π(y) be the value of y at the next such crossing, and solve Π(y)=y to machine precision (integration at 3×10−14 relative, 10−16 absolute tolerance; the extrema of x are located as roots of x˙=0 rather than by sampling). This gives the cycle's half-range A=21(xmax−xmin) and its period T directly:
K
A
A/K−3
period T
3.002
0.0632572
1.414474
15.396477
3.005
0.1000461
1.414865
15.405297
3.010
0.1415514
1.415514
15.420003
3.020
0.2003665
1.416805
15.449436
3.050
0.3176604
1.420620
15.537896
3.100
0.4511915
1.426793
15.685860
3.200
0.6433024
1.438468
15.983709
(Recomputing the whole table with the tolerance relaxed to 10−12 reproduces every entry to nine decimals.) The third column decreases monotonically towards 1.41421…=2 as K↓3, rather than growing or collapsing, so the amplitude behaves like 2K−3: rule of thumb 1. And the cycle shrinks to nothing as K→3+. Starting at (x,y)=(1.05,y∗(K)), a distance 0.05 from the equilibrium, the trajectory returns to it for K=2.98 and K=2.99: the distance at t=5000 is 9.3×10−6 and 6.7×10−4 respectively, and still falling, in both cases close to the linear prediction 0.05eαt with α=(K−3)/(4K). So there is no cycle on the stable side and no hysteresis. The bifurcation is supercritical.
Step 7: check the period against the prediction. The last column tends to 15.3906, matching 2π6=15.390598 from step 5: at K=3.002 the measured period is 15.396477, which is 0.0382% high. The discrepancy falls very nearly in proportion to K−3, as rule of thumb 2 says it should: it is 5.8789×10−3 at K=3.002 and 1.472978×10−1 at K=3.05, a ratio of 25.06 for a parameter ratio of 25.
Reading. The equilibrium is destroyed by improving the environment. Below K=3 predator and prey settle to a steady coexistence; above it they oscillate, with an amplitude that grows like K−3, so that at K=3.5 the prey already swings between 0.29 and 2.37. Enrich further and the minimum prey density keeps falling, which in a finite population means extinction. This is the paradox of enrichment, and the calculation above locates it exactly: at K=3, with a period of 2π6.
Example(A subcritical bifurcation, and the jump it causes)Original
Take the subcritical model with ω=1, b=0:
r˙=μr+r3−r5,θ˙=1.
Sweep μ upwards from −0.5, starting each time from the state reached at the previous value.
At μ=−0.5<−41 the only invariant set is the origin and the system sits there.
At μ=−0.2 the cycles are at r−=0.5257311 and r+=0.8506508 (from r±2=21(1±0.2), with 0.2=0.4472136). The state is still at the origin, inside r−, so nothing happens — the system gives no warning that two cycles now exist a short distance away.
At μ=0+ the inner cycle reaches the origin and the state is thrown out to r+=1: an O(1) jump, from amplitude 0 to amplitude 1, triggered by an arbitrarily small change in μ.
Now sweep back down. At μ=−0.2 the state stays on the outer cycle at r=0.8506508, not at the origin: the two sweeps disagree over the whole interval −41<μ<0. The oscillation only stops at μ=−41, where the outer cycle collides with the inner one at r=1/2=0.7071068 and both vanish; the state then falls to the origin.
So the response curve is a hysteresis loop, 41 wide in μ: sweeping up, the amplitude jumps from 0 to 1 at μ=0; sweeping down, it drops from 1/2=0.7071068 to 0 at μ=−41. The two branches disagree over the whole interval −41<μ<0, and the amplitude at which the oscillation finally switches off, 0.7071068, is larger than the amplitude at which it switched on, which was 0. Compare the supercritical case, where the amplitude is a continuous single-valued function of μ and the two sweeps agree everywhere. In an engineering setting this is the whole difference between a vibration that announces itself by growing slowly and one that appears at full size without notice.
Oscillating reactions. The two-variable reduction of the chlorine dioxide-iodine-malonic acid reaction,
x˙=a−x−1+x24xy,y˙=bx(1−1+x2y),
is Strogatz's worked case: the single fixed point loses stability through a Hopf bifurcation as the parameter b is lowered past bc=3a/5−25/a, and the period of the resulting chemical clock is read off from Imλ exactly as in step 5 of the method, giving T=2π[(a2+25)/(15a2−625)]1/2, which tends to 2π/15 for large a; Strogatz quotes this limit as ≈1.63, and the value is 1.62231… The useful point for an experimenter is that the theory predicts the period at onset from the linearization alone, before anything is known about the amplitude.
Circuits. Hirsch, Smale and Devaney put a tunable resistor in a Lienard circuit, x˙=y−(x3−μx), y˙=−x, where μ might be the resistor's temperature. The origin is the only equilibrium; its Jacobian (μ−110) has eigenvalues 21(μ±μ2−4), a spiral sink for −1≤μ<0 and a spiral source for 0<μ≤1. For μ≤0 the resistor is passive and the circuit is dead: every solution decays to zero, including at μ=0 where the linearization is a centre. As μ becomes positive a unique attracting periodic solution appears and shrinks back to the origin as μ→0+. The circuit starts to oscillate of its own accord, and this is the mechanism behind every electronic oscillator; at μ=1 it is the van der Pol oscillator.
Why the sub- or supercritical distinction is the practical question. Strogatz lists three places where the subcritical case appears: the firing of nerve cells, aeroelastic flutter and other vibrations of aircraft wings, and instabilities of fluid flows. In all three, what matters is not that an oscillation exists but that it arrives at finite amplitude and without warning, and that reversing the parameter does not switch it off — the hysteresis of the second worked example. A wing that flutters subcritically cannot be recovered by slowing down to the speed at which flutter began. Supercritical onset, by contrast, is announced by small oscillations whose amplitude grows like a square root, which is why an experimenter who measures that square-root law has evidence about the structure of the underlying model and not merely about where it becomes unstable.
Ecology. The worked example is the standard predator-prey model with a saturating functional response, and the Hopf bifurcation there is the paradox of enrichment: raising the prey's carrying capacity, which looks like an unambiguous improvement, converts a stable coexistence into a limit cycle whose troughs take the prey arbitrarily close to zero.
Gives the two approximations that turn a known limit cycle into numbers: in the strongly damped limit the orbit splits into slow crawls along a nullcline and instantaneous jumps between its branches, so the period is a single quadrature; in the weakly nonlinear limit two-timing collapses the oscillator to a pair of averaged equations for a slowly drifting amplitude and phase.
Definition(Lienard systems, relaxation and weakly nonlinear limits)AdaptedAdapted from Strogatz (1994), §7.4, p. 210, and §7.6, p. 215; Hirsch, Smale & Devaney (2013), §12.2, p. 261
A Lienard equation is a second-order scalar equation
x¨+f(x)x˙+g(x)=0,
read as a unit mass carrying a position-dependent damping coefficient f and a restoring force −g. Written as a planar system in the usual way it becomes
x˙=y,y˙=−g(x)−f(x)y.(L)
Throughout, f and g are continuously differentiable on all of R and
F(x)=∫0xf(u)du
denotes the antiderivative of the damping coefficient that vanishes at 0.
The van der Pol oscillator is the case g(x)=x, f(x)=μ(x2−1):
x¨+μ(x2−1)x˙+x=0,μ>0.
Its damping is negative for ∣x∣<1 and positive for ∣x∣>1, so small motions are driven and large ones are braked. Two opposite limits of this one equation are the subject of this unit.
The relaxation limit is μ≫1. The cycle then consists of long intervals of slow drift separated by short, almost instantaneous transitions, and an oscillation with that waveform is called a relaxation oscillation: tension built up slowly is released suddenly, and the build-up begins again.
The weakly nonlinear limit is the opposite extreme. A weakly nonlinear oscillator is an equation
x¨+x+ϵh(x,x˙)=0,0<ϵ≪1,
with h smooth: a small perturbation of the harmonic oscillator of The damped and driven harmonic oscillator. Van der Pol with μ=ϵ small is the case h(x,x˙)=(x2−1)x˙, and the undamped Duffing oscillator is the case h(x,x˙)=x3.
One further piece of notation is needed. Hirsch, Smale and Devaney use the name Lienard system for the first-order pair
x˙=y−φ(x),y˙=−x,(L′)
which is what the circuit of their Chapter 12 produces directly, with φ the current-voltage characteristic of the resistor. The two forms are related by φ=F and g(x)=x: we verify this in the proof below, where (L′) appears as the natural coordinate system for the relaxation limit. Their van der Pol equation is (L′) with φ(x)=x3−x, which is our van der Pol equation with μ=1 after the rescaling x↦x/3. When the letter f appears below it always means the damping coefficient in (L), never the characteristic in (L′).
Method(Lienard's theorem)AdaptedAdapted from Strogatz (1994), §7.4, Lienard's Theorem, pp. 210-211
Impose on the coefficients of (L) the following five conditions.
Both f and g are continuously differentiable on the whole real line.
g is an odd function: g(−x)=−g(x) at every x.
On the positive axis g is strictly positive: g(x)>0 whenever x>0.
f is an even function: f(−x)=f(x) at every x.
The antiderivative F(x)=∫0xf(u)du, which conditions 1 and 4 make odd and continuously differentiable, vanishes at exactly one point a of the positive axis; it takes negative values throughout 0<x<a; on x>a it is positive and never decreasing; and F(x) is unbounded above as x grows.
Under conditions 1-5 the planar system (L) possesses exactly one limit cycle. That cycle encircles the origin, and it is stable.
This is Strogatz's conclusion, and it is weaker than it may look: limit cycle means isolated closed orbit, so what is asserted is that exactly one closed orbit is isolated, not that there is exactly one closed orbit of any kind. We do not strengthen it. No proof is given here; Strogatz states the theorem and points to the monograph literature for the argument, and we follow him in that, using the conclusion only as a licence to write the cycle rather than a cycle. In the van der Pol case conditions 1-4 need no checking, and F(x)=3μx(x2−3) satisfies condition 5 with a=3; so every μ>0 gives a system with a unique stable limit cycle.
A sharper statement is available in the single case μ=1. Hirsch, Smale and Devaney prove, for (L′) with φ(x)=x3−x, that there is exactly one nontrivial periodic solution and that every solution other than the equilibrium tends to it.
Method(The period of a relaxation oscillation)AdaptedAdapted from Strogatz (1994), §7.5, Examples 7.5.1 and 7.5.2, pp. 212-214
Take g(x)=x and a damping coefficient μf(x) with μ>0, so that the equation is
x¨+μf(x)x˙+x=0,
and assume
f is even and continuously differentiable, with f(0)<0;
f has a positive zero xB, with f<0 on (0,xB) and f>0 on (xB,∞), so that xB is its only positive zero;
F(x)=∫0xf→∞ as x→∞.
Then F is odd, has its only positive minimum at xB with F(xB)<0, and there is a unique xA>xB with F(xA)=−F(xB). As μ→∞ the period of the limit cycle satisfies
T=2μ∫xBxAxf(x)dx+o(μ)
The leading term is O(μ): the cycle is slow, and almost all of the time is spent creeping along the two slow branches.
For van der Pol, f(x)=x2−1 gives xB=1 and xA=2, and the integral is 23−ln2, so
T=μ(3−2ln2)+o(μ)≈1.6137μ.
Strogatz reports a refinement of this van der Pol formula, due to Cartwright, in which the next term is +2αμ−1/3 with α=2.338 the smallest root of Ai(−α)=0; the correction comes from the time taken to round the corner between a jump and a crawl. We neither prove nor generalise that term, but the example below checks numerically that the error in the boxed formula does decay like μ−1/3 for a damping coefficient other than the cubic one.
Method(The averaged equations)AdaptedAdapted from Strogatz (1994), §7.6, pp. 218-224, Eq. (53), p. 224
Let h be smooth and consider the weakly nonlinear oscillator
x¨+x+ϵh(x,x˙)=0,0<ϵ≪1.
Introduce the fast time τ=t and the slow time T=ϵt, treat them as independent variables, and seek
x(t,ϵ)=x0(τ,T)+ϵx1(τ,T)+O(ϵ2)
with x0 and x1 bounded in τ. Then necessarily
x0(τ,T)=r(T)cos(τ+ϕ(T)),
and, writing θ=τ+ϕ and
⟨⋅⟩=2π1∫02π(⋅)dθ,h=h(rcosθ,−rsinθ),
the slowly varying amplitude r and phase ϕ obey the averaged equations
dTdr=⟨hsinθ⟩,rdTdϕ=⟨hcosθ⟩,
in which r is held fixed while the average over θ is taken. The angular frequency of the resulting oscillation is
ω=dtd(τ+ϕ(T))=1+ϵdTdϕ+O(ϵ2).
Strogatz gives the following rule of thumb for the accuracy of the leading term, and refers the precise statements to the monograph literature: x0 stays within O(ϵ) of the true solution started from the same data for all times up to and including t∼O(1/ϵ), and, if the true solution is periodic, for all time.
Idea of the proofOriginal
Both halves of the unit are the same trick: find the one coordinate in which the motion is slow, and integrate along it.
In the relaxation limit the fast coordinate is x itself. Changing variables so that the damping is absorbed into the vertical coordinate turns the equation into a system in which the horizontal speed is O(μ) everywhere except on one curve, and O(1/μ) on it. A trajectory therefore shoots horizontally onto that curve, crawls along it until the curve turns over, and shoots across again. Only the crawls take any time, and on a crawl the motion is one-dimensional, so the period is an ordinary integral.
In the weakly nonlinear limit the fast coordinate is the phase of an almost harmonic oscillation, and the slow ones are its amplitude and phase offset. Declaring the two time scales independent and refusing to let the first correction grow without bound forces the slow variables to drift at exactly the rate at which the perturbation pumps energy and advances phase per cycle; that rate is a Fourier coefficient, and extracting it is an average.
Proof(The period of a relaxation oscillation)AdaptedAdapted from Strogatz (1994), §7.5, pp. 212-214
Step 1: a coordinate in which the damping is integrated. The two terms x¨+μf(x)x˙ are the time derivative of x˙+μF(x), because dtdF(x)=f(x)x˙ by the chain rule. So with
w=x˙+μF(x)
the equation x¨+μf(x)x˙+x=0 says exactly w˙=−x. Rescaling the new coordinate as y=w/μ gives the planar system
x˙=μ[y−F(x)],y˙=−μx.(∗)
This is the form (L′) of the definition up to the factors of μ, which are what carry the limit.
Step 2: everywhere off the curve y=F(x) the motion is horizontal. Suppose ∣y−F(x)∣ and ∣x∣ are both of order 1. Then (∗) gives ∣x˙∣∼μ and ∣y˙∣∼1/μ, so the velocity vector is longer horizontally than vertically by a factor μ2. The trajectory is therefore horizontal to within O(μ−2), and its direction is set by the sign of y−F(x): rightward above the curve, leftward below it.
Whether that horizontal run closes the gap depends on the branch it is heading for, and here the sign of f decides. Differentiating the gap along the flow and using (∗),
dtd(y−F(x))=y˙−f(x)x˙=−μx−μf(x)(y−F(x)),
whose dominant term shrinks the gap at rate μf(x) where f>0 and grows it at rate μ∣f(x)∣ where f<0. So the jump ends on one of the two outer branches of the curve, ∣x∣>xB, where f>0; the middle branch ∣x∣<xB repels horizontally and is never followed. (Strogatz argues the half of this that he needs, the approach to the outer branch from above; the sign test is what the statement here needs in the extra generality of a general even f.) Crossing a horizontal distance of order 1 at speed μ takes a time of order μ−1: this is the jump.
Step 3: on the curve the motion is slow and one-dimensional. The trajectory cannot stay off the curve, and it cannot cross it transversally at speed μ, so it must settle into the region where y−F(x) is small enough for the two components of (∗) to be comparable, namely y−F(x)=O(μ−2). To leading order the trajectory then follows y=F(x). Differentiating that relation along the trajectory and using the second equation of (∗),
f(x)x˙=dtdF(x)=y˙=−μx,sox˙=−μf(x)x.
The speed is O(μ−1): this is the crawl. Separating variables, the time to traverse a piece of the curve is
dt=−xμf(x)dx.
Step 4: the endpoints of a crawl. The branch of y=F(x) to the right of xB is where F increases (f>0 there), and it is the branch the flow descends: on it x>0, f(x)>0, so x˙<0. The crawl ends where the branch ends, at the local minimum x=xB of F, since F′(xB)=f(xB)=0. There the trajectory jumps horizontally, at the level y=F(xB), to the only other point of the curve at that height. Because F is odd, that point is x=−xA, where xA>xB is determined by
F(xA)=−F(xB).
Such an xA exists and is unique: on (xB,∞) the function F increases strictly from F(xB)<0 to +∞, so it takes the value −F(xB)>0 exactly once. By the oddness of F and x↦−x, y↦−y symmetry of (∗), the left branch is traversed from −xA to −xB in the same time, after which the trajectory jumps to +xA and the cycle repeats.
Step 5: the period. The two crawls contribute equally, and each takes
∫dt=∫xAxB(−xμf(x))dx=μ∫xBxAxf(x)dx,
which is positive because f>0 and x>0 on (xB,xA). The two jumps contribute O(μ−1) each. Hence
T=2μ∫xBxAxf(x)dx+o(μ),
as claimed. □
Proof(The averaged equations)AdaptedAdapted from Strogatz (1994), §7.6, pp. 219-224
Step 1: the chain rule in two times. With τ=t and T=ϵt regarded as independent, any function u(τ,T) evaluated along the solution has
The O(1) equation is the harmonic oscillator in the fast variable only, so its general solution is x0=rcosτ+ssinτ with r and s constant in τ — but free to depend on T, which is the whole point of the method. Writing the same thing in amplitude-phase form,
x0(τ,T)=r(T)cosθ,θ=τ+ϕ(T).
Step 3: what the right-hand side looks like. From x0=rcosθ we get ∂τx0=−rsinθ and ∂Tx0=r′cosθ−rϕ′sinθ, where ′ is d/dT; differentiating the latter in τ,
Step 4: kill the resonance. In the fast variable τ, with T frozen, the left-hand side is the harmonic oscillator of unit frequency, and θ differs from τ by the constant ϕ. A forcing term proportional to sinθ or cosθ is therefore resonant and produces a response growing like τ (Regular perturbation expansions); such an x1 is not bounded in τ, so the expansion assumed at the outset would be invalid once τ∼1/ϵ. We must make the first harmonic of the right-hand side vanish.
Since h(θ) is 2π-periodic in θ, expand it in a Fourier series (Fourier series),
These are the averaged equations. No other harmonic is resonant, so with this choice x1 can be solved for as a bounded 2π-periodic function of θ plus a harmonic term, and the expansion is consistent.
Step 5: the frequency. The argument of the cosine in x0 is θ=τ+ϕ(T)=t+ϕ(ϵt), so
ω=dtdθ=1+ϵϕ′(T),
correct through first order in ϵ; the O(ϵ2) frequency shift is invisible at this order and would require a third time scale t2=ϵ2t. □
Remark(Where the hypotheses are used, and what fails without them)Original
Relaxation. Evenness of f and oddness of g enter only in Step 4, to say that the two crawls are mirror images and so contribute equally; without the symmetry the same argument gives the period as a sum of two different integrals, one per branch. The requirement that f has a single positive zero is what makes y=F(x) have exactly one knee on each side, hence exactly one jump per half-cycle; a damping coefficient with two positive zeros produces a nullcline with two knees and a correspondingly more elaborate itinerary. The condition f(0)<0 makes the origin a source, so that the system is not simply damped to rest. Finally the estimate is asymptotic in μ only: the neglected corner regions contribute at order μ−1/3, which is larger than the O(μ−1) jumps and is the true next term.
Averaging. The derivation needs h smooth enough for the Fourier series of h(θ) to have well-defined first harmonics, which continuity is already enough for, and needs the unperturbed problem to be the harmonic oscillator, so that the kernel of ∂ττ+1 is spanned by cosθ and sinθ and "resonant" means "first harmonic". If the unperturbed restoring force is nonlinear, every harmonic of the unperturbed orbit can resonate, and the method must be rebuilt in action-angle variables.
The boundedness of x1 in τ is a requirement imposed by us, not a consequence of the equation; it is what selects the averaged equations out of the many ways of splitting the solution between x0 and x1. Drop it and the expansion still solves the equation order by order, but it reproduces the useless secular series of Regular perturbation expansions.
Nothing in the derivation justifies treating τ and T as independent, and Strogatz says plainly that this is hard to make rigorous. What is being asserted is only an asymptotic approximation on a stated time interval, and the rule of thumb quoted in the statement is reported as a rule of thumb. The honest check is the one carried out in the examples below: solve the full equation numerically and compare.
A warning that both halves share: the two limits are limits of the same equation, and the formulas do not meet in the middle. The relaxation formula predicts T≈1.6137μ for van der Pol and the averaging calculation predicts T≈2π; at μ=1 these read 1.61 and 6.28, and neither is right (the true period is about 6.66).
Method(Getting numbers out of a limit cycle)AdaptedAdapted from Strogatz (1994), §7.4, Example 7.4.1, p. 211, §7.5, pp. 212-214, and §7.6, pp. 223-226
Given x¨+f(x)x˙+g(x)=0, or any planar system suspected of oscillating:
Establish that a cycle exists before approximating it. Either check the five hypotheses of Lienard's theorem, or build a trapping region and apply Limit cycles and the Poincare-Bendixson theorem.
Decide which limit you are in. Measure the damping against the restoring force. If the damping coefficient is large compared with the frequency, you are in the relaxation limit; if it is small, in the weakly nonlinear limit. Introduce an explicit parameter (μ or ϵ) that makes this quantitative, rescaling x and t if necessary.
Relaxation limit (μ≫1, equation x¨+μf(x)x˙+x=0):
Form F(x)=∫0xf and pass to the Lienard plane x˙=μ[y−F(x)], y˙=−x/μ.
Locate the knee xB, the positive zero of f (there is only one, by the hypothesis on the sign of f), and the landing point xA>xB, the unique solution of F(xA)=−F(xB).
Evaluate T≈2μ∫xBxAf(x)dx/x. Expect a relative error of order μ−4/3.
Substitute x=rcosθ, x˙=−rsinθ into h and expand the result in powers of cosθ and sinθ.
Average term by term, using ⟨cos2⟩=⟨sin2⟩=21, ⟨cos4⟩=⟨sin4⟩=83, ⟨cos2sin2⟩=81, ⟨cos6⟩=⟨sin6⟩=165, and ⟨⋅⟩=0 for any monomial cosmθsinnθ with m+n odd.
Read off r′=⟨hsinθ⟩ and rϕ′=⟨hcosθ⟩. The limit cycles are the positive zeros of r′, their stability is that of the corresponding fixed points of the one-dimensional flow r′=⟨hsinθ⟩ (Flows on the line), their amplitude is the root r∗, and their frequency is 1+ϵ⟨hcosθ⟩/r∗.
For the transient, solve the scalar equation for r(T) and restore T=ϵt.
Example(A quartic relaxation oscillator, and the size of the error)Original
Take the damping coefficient f(x)=x4−1, so the equation is
x¨+μ(x4−1)x˙+x=0,μ≫1.
Lienard's theorem applies.f and g(x)=x are polynomials, g is odd and positive on x>0, and f is even. Here
F(x)=5x5−x=5x(x4−5),
whose only positive zero is a=51/4≈1.49535; F<0 on (0,a); F′=f>0 for x>1, so F is increasing, hence positive and nondecreasing, on (a,∞); and F→∞. So there is a unique stable limit cycle.
The knee and the landing point. The positive zero of f is xB=1, and F(1)=51−1=−54. We need xA>1 with F(xA)=54, that is
5x5−x=54⟺x5−5x−4=0⟺(x+1)(x4−x3+x2−x−4)=0.
The quartic factor is negative at x=1 (value −4) and positive at x=2 (value 16−8+4−2−4=+6), and its derivative 4x3−3x2+2x−1 is positive for x≥1, so it has exactly one root there. Newton's method from x=1.7 gives
Numerically xA4=7.4233183, so xA4/4=1.855830 and lnxA=0.501157, giving the integral 1.104673 and
T=2.209346μ+o(μ)
Check against the full equation. Integrating x˙=μ[y−F(x)], y˙=−x/μ numerically (DOP853/LSODA, relative tolerance 10−11) from (2,0) and measuring the time between successive upward crossings of x=0 on the settled orbit:
μ
T (numerical)
T/μ
T−2.209346μ
10
24.82190
2.482190
2.72844
100
222.22477
2.222248
1.29016
1000
2209.94354
2.209944
0.59745
T/μ converges to the predicted 2.209346, and the residual shrinks by factors 2.115 and 2.159 when μ is multiplied by 10 — against 101/3=2.15443. So the first correction to the boxed formula is of order μ−1/3 here too, exactly as in the van der Pol case for which Cartwright computed its coefficient.
Example(A Rayleigh-Duffing oscillator by averaging)Original
Consider
x¨+x+ϵ(x˙3−2x˙+x3)=0,0<ϵ≪1,
a cubically damped oscillator — negative damping at small speeds, positive at large ones — on a hardening spring. Here h(x,x˙)=x˙3−2x˙+x3.
Substitute. With x=rcosθ and x˙=−rsinθ,
h(θ)=−r3sin3θ+2rsinθ+r3cos3θ.
Average. Multiplying by sinθ and averaging, the term cos3θsinθ has odd total degree and averages to zero, so
Multiplying by cosθ instead, both sin3θcosθ and sinθcosθ average to zero and only the Duffing term survives:
rϕ′=⟨hcosθ⟩=r3⟨cos4θ⟩=83r3,soϕ′=83r2.
Read off the cycle. The scalar flow r′=83r(38−r2) on r>0 has the single fixed point
r∗=8/3=326=1.6329932…,
and r′ is positive below it and negative above it, so it is stable: a stable limit cycle of amplitude r∗+O(ϵ). On it ϕ′=83⋅38=1 exactly, so
ω=1+ϵ+O(ϵ2),T=1+ϵ2π+O(ϵ2).
The transient. Put u=r2. Then u′=2rr′=2u−43u2, a logistic equation with growth rate 2 and capacity 38, so
u(T)=1+Ce−2T8/3,C=3r(0)28−1,
and with r(0)=1 (so C=5/3) the envelope is r(t)=[(8/3)(1+35e−2ϵt)−1]1/2.
Check against the full equation. Integrating the second-order equation numerically (DOP853, relative tolerance 10−11) from x(0)=1, x˙(0)=0 and reading the amplitude and the period off the settled orbit:
ϵ
amplitude
error vs 1.632993
ω
error vs 1+ϵ
0.10
1.538603
−0.094390
1.082264
−0.017736
0.05
1.579771
−0.053222
1.044980
−0.005020
0.02
1.609904
−0.023089
1.019124
−0.000876
The amplitude error falls by roughly the factor by which ϵ falls — it is O(ϵ), as advertised. The frequency error falls by 3.5 and then by 5.7 as ϵ is halved and then divided by 2.5, against 4 and 6.25 for a quadratic: it is O(ϵ2), which is why the frequency prediction is the more accurate of the two.
Application(Nerve cells, bowed strings and clocks)AdaptedAdapted from Strogatz (1994), §7.5, p. 212, and §7.6, pp. 226-227; Hirsch, Smale & Devaney (2013), §12.1, pp. 257-260, and §12.3, p. 263
Van der Pol wrote his equation for a vacuum-tube oscillator, and the circuit that produces it is still the cleanest derivation: in a series circuit with one resistor whose current-voltage characteristic is not monotone, Kirchhoff's laws and the capacitor and inductor relations reduce the state to the inductor current and the capacitor voltage, and give exactly x˙=y−φ(x), y˙=−x. The non-monotone piece of φ is the part of the characteristic where the device acts as a source rather than a sink of energy, and it is what makes the circuit oscillate on its own instead of ringing down. A linear circuit can never behave this way, which is what the uniqueness theorem records: an isolated periodic orbit attracting all other solutions is impossible for a linear system.
The relaxation limit is the normal regime for excitable biological tissue. A neuron at rest integrates input slowly, fires in a millisecond, and resets: the slow-fast decomposition of this unit is the skeleton of the FitzHugh-Nagumo reduction of the Hodgkin-Huxley equations, with the cubic nullcline playing the part of the membrane's current-voltage curve and the knee playing the part of the firing threshold. The same picture, with mechanical names, describes stick-slip motion: a bowed violin string, a squealing brake, an earthquake fault. The practical payoff of the formula in this unit is that the period of such an oscillator is set by the slow variable alone, so it scales linearly with the large damping parameter and is almost independent of the details of the fast jump.
Averaging answers a different question, and the one experimenters usually ask about a nearly sinusoidal oscillator: how do amplitude and frequency depend on each other? For a linear oscillator the frequency does not depend on the amplitude at all; the term rϕ′=⟨hcosθ⟩ measures exactly the failure of that isochronism. In the example above it is a hardening spring that shifts the frequency up; a pendulum, whose restoring force softens, shifts it down, which is why a pendulum clock must be driven at a fixed amplitude to keep time. The same two averaged equations are the standard tool for the response of a driven weakly nonlinear oscillator, where fixed points of the averaged system are phase-locked oscillations and saddle-node bifurcations of those fixed points are the jumps of the familiar bent resonance curve (see Global bifurcations of cycles).
Names and analyses the three ways a planar limit cycle can be created or destroyed without shrinking into a fixed point - two cycles colliding, a fixed point appearing on the cycle, and the cycle swelling onto a saddle - and gives the amplitude and period scaling laws by which each can be told from the others in measured data.
Definition(Local and global bifurcations of cycles)AdaptedAdapted from Strogatz (1994), §8.4, pp. 260-263
Let x˙=f(x,μ) be a smooth family of vector fields on an open subset of the plane, depending on a real parameter μ. A bifurcation of cycles at μ=μc is a change, as μ crosses μc, in the number or stability of the closed orbits of the family.
The bifurcation is local if everything that happens is visible in an arbitrarily small neighbourhood of a single fixed point: the Hopf bifurcation of The Hopf bifurcation is local, because the cycle it creates shrinks onto the fixed point as μ→μc, and the eigenvalues of the Jacobian at that one point announce it. The bifurcation is global if it is not: the cycle involved stays a finite distance across at the bifurcation, so no amount of magnification around any one point reveals it, and no linearization detects it. Strogatz lists four common ways a planar limit cycle is created or destroyed; one is the Hopf bifurcation, and the other three are global. They are:
Saddle-node (fold) bifurcation of cycles. Two limit cycles, one stable and one unstable, approach each other, merge into a single half-stable cycle at μc, and annihilate. Read on a cross-section this is a saddle-node bifurcation of fixed points of the return map, which is where the name comes from. No fixed point of the flow takes part.
Infinite-period bifurcation. A fixed point of the flow appears on the cycle. The cycle develops a bottleneck - a short arc along which the speed drops towards zero - and the period grows without bound; at μc the orbit is no longer periodic, having become a fixed point together with the one trajectory that leaves it and returns to it, namely the rest of the former cycle. Past μc the fixed point splits into two, and the former cycle is then two fixed points joined by a pair of arcs, one running each way.
Homoclinic (saddle-loop) bifurcation. The cycle swells until it touches a saddle point lying outside it. At μc the cycle has become a homoclinic orbit of that saddle, a trajectory that approaches the same fixed point as t→+∞ and as t→−∞ (Conservative and reversible systems). Past μc the saddle connection breaks: the branch of the unstable manifold that used to close up now misses the stable manifold on one side, and no closed orbit remains nearby. The period again diverges, so this too is an infinite-period bifurcation in the literal sense; the names homoclinic and saddle-loop are used to keep the two cases apart.
The letter μ is used below in two ways, and they should be kept apart. In the three prototype systems it is the system's own parameter and μc is its bifurcation value, so that μc=−41 in the first prototype and μc=1 in the second, and the distance from the bifurcation is ∣μ−μc∣. In the scaling table, and there only, μ denotes that distance itself: a dimensionless quantity with 0<μ≪1.
Definition(Three prototypes)AdaptedAdapted from Strogatz (1994), §8.4, pp. 261-264
(a) Saddle-node of cycles. Let ω and b be real numbers such that ω+bρ=0 for every ρ∈[0,1] (it is enough that ω=0 and that b has the same sign as ω, or that ∣b∣<∣ω∣), and let
r˙=μr+r3−r5,θ˙=ω+br2.
This is the subcritical Hopf family of The Hopf bifurcation, now examined away from μ=0. For −41<μ<0 it has two limit cycles, of radii r±=(21(1±1+4μ))1/2, the outer one stable and the inner one unstable, and the origin is also stable. As μ↓−41 the two radii merge at
rc=1/2=0.7071067…,
leaving a single half-stable cycle, and for μ<−41 there is no cycle at all. The origin is a stable spiral throughout and takes no part. At the bifurcation the cycle has amplitude rc and period 2π/∣ω+b/2∣, both O(1).
The condition on ω+bρ is ours, not Strogatz's: he reads only the orders of magnitude off this example and needs nothing beyond ω=0. It cannot be dropped once a period is quoted. Taking b=−2ω, for instance, makes θ˙ vanish on the circle r2=21, which is then a circle of fixed points rather than a closed orbit, and 2π/∣ω+b/2∣ is not a period but a division by zero.
(b) Infinite-period. For μ>0 let
r˙=r(1−r2),θ˙=μ−sinθ.
The radial and angular equations are uncoupled. Every trajectory except the origin approaches the unit circle, which is invariant. On it the angular equation is the nonuniform oscillator of Flows on the circle. For μ>1 the circle is a limit cycle, traversed with period
T(μ)=∫02πμ−sinθdθ=μ2−12π∼μ−1π2(μ↓1).
At μ=1 a half-stable fixed point is born on the circle at θ=π/2, and for μ<1 it splits into a saddle and a stable node. The amplitude is exactly 1 throughout; only the period diverges.
(c) Homoclinic. Let f be smooth in x and in μ, let p(μ) be a hyperbolic saddle of f(⋅,μ) with eigenvalues λu>0>−λs depending smoothly on μ, and suppose that at μ=μc one branch of its unstable manifold coincides with one branch of its stable manifold, forming a homoclinic loop Γ. Fix a cross-section of Γ near the saddle and let u0(μ) be the signed separation of the two manifolds measured on it, counted positive when the returning branch of Wu(p) lies on the interior side of the loop, so that u0(μc)=0; assume the manifolds separate at a nonzero rate, u0′(μc)=0. Write λu and λs for the two eigenvalues at μ=μc and assume λu=λs.
Then there is at least one closed orbit near Γ for μ on one side of μc and none on the other, and which side that is turns on the comparison of the two eigenvalues as well as on the sign of the separation: the orbits live where u0(μ)>0 when λu<λs, and where u0(μ)<0 when λu>λs. Every closed orbit near Γ crosses the entry section the proof below works with at a distance from the invariant curve that that section crosses comparable to the separation of the two manifolds read on that same section; writing σ(μ)>0 for that separation, on the side of μc where the orbits exist, the crossing lies between σ and 2σ. The two cases differ in which section, which curve and which separation these are. When λu<λs the section is the edge of the proof's small square that crosses the stable manifold, the distance is measured from the stable manifold, and σ=u0. When λu>λs the section is, after the proof's reversal of time, the edge that crosses the unstable manifold, the distance is measured from the unstable manifold, and σ is the separation read on that edge, which Step 6 of the proof computes as u~0=−u0/a+O(u02) with a>0: comparable in size to ∣u0∣ and positive where the orbits exist, but of the opposite sign to u0, which is why the two cases carry the orbits on opposite sides of μc. On another cross-section of the same branch the bracket holds up to a positive factor bounded away from 0 and from ∞, but not on a section of the other branch: when λu>λs the cycle crosses the incoming branch at a distance from Ws of order ∣u0∣λu/λs, which is far smaller than ∣u0∣. In either case all the closed orbits lie in an annulus that shrinks onto Γ as μ→μc. The set of them attracts a one-sided neighbourhood of Γ if λu<λs, and repels one if λu>λs. Attracts means here that the forward orbit of every point of that neighbourhood converges to the set of closed orbits. The side in question is the interior of the loop, the quadrant u>0, s>0 in the coordinates of the proof below, in both cases; the proof controls the return map only there. It is not in general the side into which the returning branch of Wu(p) enters: when λu>λs that branch returns to the exterior, while the repelled neighbourhood is still the interior. Nothing is claimed about the other side, and nothing here distinguishes one closed orbit from another: a point beyond the stable manifold leaves along the second branch of Wu(p), and whether it ever comes back depends on global structure that the hypotheses above do not fix - in a figure-eight configuration, where the other two branches also close up, such a point can be recaptured. As μ→μc from the side where they exist, the period of each of them satisfies
Their amplitudes tend to the diameter of Γ, which is O(1). The crossing point on the section converges to Γ by the above, and continuous dependence on initial conditions over the compact outside arc of Step 2, which is traversed in a bounded time, carries that convergence to the part of the orbit lying outside the square N of the proof. Inside N continuous dependence is of no use, since the transit time there diverges; what holds the orbit near Γ inside N is the size of the square itself, a side δ that may be fixed as small as one pleases before μ is brought to μc.
Two words are deliberately absent. One is unique: the statement claims at least one closed orbit, not exactly one. The other is limit cycle, which means an isolated closed orbit. One short further argument, sketched at Step 4 of the proof, supplies both - uniqueness, isolation, and with them the stability or instability of the single cycle - by showing that the return map is a contraction; but it needs a C1 linearization at the saddle, which this handbook does not prove and which the proof below therefore avoids. Everything asserted above, the period law included, is proved without it.
Three things about this coefficient. It is the reciprocal of the smaller of the two rates, so it is 1/λu exactly in the stable case; the minimum has to be there because reversing time exchanges λu with λs and interchanges the two cases while leaving every period unchanged. It is ours, not Strogatz's: his table records only the order ln(1/μ) and refers the estimate to the literature. And the error term is o(1) relative, rather than the O(1) absolute error one often sees quoted, because the proof below uses only the invariant curves at the saddle; with a C1 linearization there - which in the plane is a theorem of Hartman for C2 fields, not something this handbook proves - the remainder is indeed O(1).
The scaling table. Writing μ for the dimensionless distance from the bifurcation, the generic laws for bifurcations of cycles in the plane are, as Strogatz tabulates them:
bifurcation
amplitude of the stable limit cycle
period
supercritical Hopf
O(μ1/2)
O(1)
saddle-node of cycles
O(1)
O(1)
infinite-period
O(1)
O(μ−1/2)
homoclinic
O(1)
O(lnμ)
The heading of the amplitude column is Strogatz's and the restriction in it is not decorative: the laws are for the stable cycle, and part (c) shows that an unstable cycle obeys the logarithm with the other eigenvalue in it. Strogatz states the homoclinic row without derivation, referring to the literature; the proof below supplies the estimate. The O(lnμ) entry means a divergence of that order as μ→0, that is, growth like ln(1/μ).
Idea of the proofOriginal
Each prototype is a two-dimensional problem with a one-dimensional bottleneck, and in each case the bifurcation happens in the one dimension that matters.
For (a) the radial equation is a flow on a half-line, uncoupled from the angle, so a collision of limit cycles is a collision of fixed points on that line. For (b) the angular equation is a flow on a circle, and the cycle's period is the time to go once around; a saddle-node on that circle makes the integrand blow up non-integrably at one point, and the familiar square-root law of the saddle-node bottleneck gives the rate. For (c) the slow region is a neighbourhood of the saddle, where the flow is linear to leading order; the time to pass a hyperbolic saddle entering at distance d from its stable manifold is λu−1ln(1/d), and the splitting between the two invariant manifolds, which is what sets d, shrinks linearly as the parameter approaches its bifurcation value. Whether the cycle really enters at distance comparable to the splitting is a question about the return map, and the answer depends on which rate at the saddle is the larger: if the saddle contracts faster than it expands, the return map lands every point of the section beyond the splitting and overshoots it by at most half that point's own distance from the stable manifold, which traps every closed orbit between the splitting and twice the splitting, and the answer is yes, giving the coefficient 1/λu; if not, one reads the same estimate backwards in time, which exchanges the two eigenvalues. Either way the coefficient belongs to the smaller eigenvalue. In all three the amplitude stays O(1) simply because the cycle at the bifurcation is a curve of finite size, which is exactly what makes the bifurcation global.
Proof((a) Saddle-node of cycles)AdaptedAdapted from Strogatz (1994), §8.4, pp. 261-262
Because r˙ does not involve θ, every circle r=const on which r˙ vanishes is invariant, and it is a single closed orbit rather than a curve of fixed points because the hypothesis ω+bρ=0 keeps θ˙ from vanishing anywhere on it. Conversely a closed orbit must have r constant, since r(t) is monotone whenever r˙=0. So the closed orbits are exactly the positive zeros of
r˙=r(μ+r2−r4),
and their stability as cycles is the stability of those zeros for the one-dimensional flow (Flows on the line). Substituting ρ=r2>0, the positive zeros solve ρ2−ρ−μ=0, that is
ρ±=21±1+4μ.
For −41<μ<0 both roots are real and positive, giving two cycles; at μ=−41 the discriminant vanishes and they coincide at ρ=21, so rc=1/2; for μ<−41 there is no real root and no cycle. This is precisely the saddle-node normal form of Bifurcations in one dimension in the variable ρ, since μ+ρ−ρ2=(μ+41)−(ρ−21)2 identically, which is the normal form with bifurcation parameter μ+41.
At μ=−41 the sign of r˙ is negative on both sides of rc, so the cycle is half-stable: attracting from outside, repelling from inside.
The origin plays no part: r˙≈μr near r=0, so it is a stable spiral for every μ<0, including at μc. Finally the period of either cycle is 2π/∣θ˙∣=2π/∣ω+bρ∣, which at ρ=21 is 2π/∣ω+b/2∣: finite and bounded away from 0 and ∞ as μ→μc. Hence both amplitude and period are O(1). □
Proof((b) Infinite-period)AdaptedAdapted from Strogatz (1994), §8.4, p. 262, and §4.3, pp. 98-100
The radial equation r˙=r(1−r2) is independent of θ, and on r>0 it has the single stable fixed point r=1; so every trajectory other than the origin tends to the unit circle, which is invariant.
On that circle the dynamics is θ˙=μ−sinθ, a flow on the circle. By the dichotomy of Flows on the circle, if μ>1 then μ−sinθ>0 everywhere, the motion is periodic, and its period is the integral of dθ/θ˙ once around, evaluated there as
T(μ)=∫02πμ−sinθdθ=μ2−12π.
Every trajectory off the circle spirals onto it, so the circle is a limit cycle. If μ<1 the equation sinθ=μ has two solutions on the circle; both are fixed points of the full planar system, lying on the unit circle, one a saddle and one a stable node, and no closed orbit survives. At μ=1 the two merge at θ=π/2.
For the rate, factor the square root: μ2−1=μ−1μ+1, and μ+1→2 as μ↓1, so
T(μ)=2μ−12π(1+O(μ−1))=μ−1π2+O((μ−1)1/2).
The exponent −21 is not an accident of this example. Near θ=π/2 write θ=π/2+ψ; then θ˙=μ−cosψ=(μ−1)+21ψ2+O(ψ4), the saddle-node normal form, and the time to cross the bottleneck is ∫dψ/((μ−1)+21ψ2)=π2/μ−1 when taken over the whole line - the same square-root law that governs a saddle-node ghost in Flows on the circle, and here it reproduces the leading term exactly. □
Proof((c) The logarithm at a homoclinic bifurcation)Original
Put the saddle at the origin. The period of the cycle is the time for one circuit, and the circuit has a slow part, a passage close to the saddle, and a fast part, the rest of the loop. Everything is in getting the slow part right, and for that we need coordinates near the saddle.
We do not take them from Linearization and the Hartman-Grobman theorem. The Hartman-Grobman theorem stated there supplies a homeomorphism, and although it preserves the time parametrisation, and so would give the passage time, the argument below also measures the separation of two invariant manifolds in the original coordinates and differentiates it in μ; a topological conjugacy, constructed afresh at each parameter value, carries no information of that kind. We therefore work with the invariant curves themselves.
Step 0: coordinates in which the manifolds are the axes. A hyperbolic planar saddle has a local stable curve and a local unstable curve through it, tangent to the two eigendirections, invariant under the flow, and for a smooth vector field they are smooth (Hirsch, Smale and Devaney (Hirsch, Smale & Devaney (2013), §8.3, the Stable Curve Theorem and the remarks following it, p. 169)). That they also depend smoothly on a parameter is standard and is used here without proof. Take the two curves as coordinate axes: for each μ near μc there is a smooth change of coordinates near the saddle, depending smoothly on μ, in which the unstable curve is the axis {s=0} and the stable curve is the axis {u=0}, oriented so that the loop leaves along u>0 and returns along s>0. Each axis is invariant, so the u-component of the field vanishes where u=0 and the s-component where s=0, and a smooth function vanishing on a coordinate axis is that coordinate times a continuous factor. Hence
u˙=uA(u,s),s˙=−sB(u,s),A(0,0)=λu,B(0,0)=λs,
with A and B continuous. Fix a tolerance ϵ with 0<2ϵ<∣λu−λs∣ and ϵ<min(λu,λs); the second constraint is what keeps λu−ϵ positive, so that the bound of Step 1 and the exponent ν+=(λs+ϵ)/(λu−ϵ) mean what they are meant to. Since A and B are continuous in (u,s) and in μ, one square
N={∣u∣≤δ,∣s∣≤δ}
serves for all μ in some interval around μc, with ∣A−λu∣≤ϵ and ∣B−λs∣≤ϵ on it. Both δ and the interval depend on ϵ, but not on μ; at the end we let ϵ→0.
Step 1: crossing the square. Let a trajectory enter N through the top edge at (uin,δ) with 0<uin≤δ. Then u grows and s decays, and s stays positive because a decaying exponential never reaches 0; so the trajectory leaves through the side edge u=δ, at a time tN that the bounds on A sandwich:
λu+ϵ1lnuinδ≤tN≤λu−ϵ1lnuinδ.
This is the only slow part of the circuit, and it is slow precisely because the velocity vanishes at the saddle. The exit height is sout=δexp(−∫0tNBdt), and putting the same bounds into the exponent gives
both exponents tending to ν=λs/λu as ϵ→0. The map uin↦sout is increasing, since trajectories of a planar flow cannot cross.
Step 2: the rest of the circuit takes a bounded time. Outside N the loop Γ is a compact arc on which f does not vanish, so the speed has a positive minimum there. The flow therefore carries a neighbourhood of the exit edge u=δ to the entry edge s=δ in a time tΓ bounded above and below by constants depending on δ and f but neither on μ nor on the exit point, and the resulting transition map is a diffeomorphism between two transversals:
uin=u0(μ)+asout+O(sout2),
with a positive and bounded away from 0 and from ∞. Its value at sout=0 is where the unstable curve comes back to the entry edge, which is the separation u0(μ) of the two manifolds measured there.
Step 3: the separation is linear in the parameter. At μ=μc the manifolds coincide, so u0(μc)=0, and the transversality assumed in the statement says u0′(μc)=0. Hence
u0(μ)=c(μ−μc)+O((μ−μc)2),c=0.
Only ln∣u0∣ enters below, so the constant c, and equally the bounded distortion that the change of coordinates of Step 0 applies to any such separation, contribute O(1) and may be ignored.
Which cross-section the separation is read on does not matter either, and this is worth saying once because the statement fixes a section in advance while the proof has measured u0 on the entry edge s=δ of a square chosen only in Step 0. Between any two cross-sections of Γ near the saddle the flow furnishes a transition map that is a diffeomorphism, depends smoothly on μ, preserves orientation, and carries the point where the stable curve meets the first section to the point where it meets the second. A separation measured on the first section is therefore the separation measured on the second times a positive factor, smooth in μ and bounded away from 0 and ∞. Hence the vanishing at μc, the sign of u0 - which, together with the comparison of the two eigenvalues (Step 6), is what decides the side of μc on which the cycle lives - and the non-degeneracy u0′(μc)=0 all transfer unchanged from one section to any other.
Step 4: where the cycle enters, from the return map. A closed orbit near Γ meets the top edge in a fixed point of the return map R=(Step 2)∘(Step 1), and by those two steps R is continuous and increasing with
u0(μ)+a−δ(δu)ν+≤R(u)≤u0(μ)+a+δ(δu)ν−
for constants 0<a−≤a+ independent of μ.
Suppose first that λu<λs. The choice 2ϵ<λs−λu makes ν−>1, so the nonlinear term is small compared with its argument: there is u1>0 with a+δ(u/δ)ν−≤u/2 for 0<u≤u1. Take the side of μc on which b:=u0(μ)>0, and let μ be close enough to μc that 2b≤u1. Then R(b)≥b and R(2b)≤b+b=2b, so by the intermediate value theorem R has a fixed point u∗ with
u0(μ)≤u∗≤2u0(μ):
a cycle exists. The same inequality shows that every fixed point in (0,u1] lies in [u0,2u0], since R(u)≥u0 always and R(u)≤u0+u/2<u once u>2u0: whatever closed orbits sit near Γ, they all enter the square at a distance comparable to the separation of the manifolds. That they enter it at all, in that interval, is worth a line: a closed orbit confined to a tube of radius r≪δ about Γ must cross the entry edge, since Γ crosses it transversally at (0,δ) and meets it nowhere else, and must cross it at some u with ∣u∣<r; it cannot cross at u=0, a point of the stable curve, whose orbit tends to the saddle and is not closed, nor at u<0, for there ∣u∣ grows and the trajectory leaves N through the edge u=−δ, following the other branch of the unstable curve. Shrink δ once more, if need be, so that Γ meets N only in the two positive semi-axes; that edge is then a compact set disjoint from Γ, hence at a positive distance from it, and a tube thinner than that distance cannot contain the trajectory. Taking r smaller than both that distance and u1 puts the crossing in (0,u1], where the two bounds apply. That is the step that cannot be taken on trust, and it is where the comparison of the two eigenvalues enters.
The same two inequalities make R map (0,u1] into itself, and R is increasing, so the iterates of any point of that interval are monotone and converge to a fixed point: the closed orbits near Γ, confined to the thin annulus just described, attract as a set every trajectory that meets the entry edge in (0,u1], whether it meets it inside that annulus or outside it. Those points include a whole one-sided neighbourhood of Γ, namely the interior side - which in this case, u0>0, is also the side into which the unstable curve returns. For a point lying close enough to the arc of Γ outside N, on its interior side, is carried to the entry edge at a small positive u by the transition map of Step 2; and a point of N on the interior side has u>0 and s>0, so it leaves through the exit edge u=δ and is then carried to the entry edge in the same way. Making the neighbourhood thin enough puts every one of those landing points in (0,u1]. That is the one-sided attraction asserted in (c), and it is exactly as much as these bounds give. They do not give isolation. An increasing continuous self-map of an interval may have a repelling fixed point, with R(u)<u just below it and R(u)>u just above, and the two bounds u0≤R(u)≤u0+u/2 do not forbid such a configuration inside [u0,2u0]; nor do they forbid a whole interval of fixed points there. Uniqueness is nevertheless true, and with it isolation and the stability of the single cycle, but it needs an ingredient we have avoided: with a C1 linearization at the saddle one computes R′(u)=O(uν−1)→0, so that R is a contraction on (0,u1] and has one fixed point, attracting. The period estimate of Step 5 uses none of this - only u∗∈[u0,2u0], which is established above for every fixed point. On the other side of μc, where u0<0, the upper bound gives R(u)≤u0+u/2<u for every 0<u≤u1, so there is no fixed point and no cycle.
Step 5: add the two times. With T(μ)=tN+tΓ, where tΓ depends on the exit point as well as on μ and only its uniform two-sided bound from Step 2 is used, Steps 3 and 4 give ln(δ/u∗)=ln(1/∣μ−μc∣)+O(1), and Steps 1 and 2 then give
with O(1) terms that depend on ϵ but not on μ. Dividing by ln(1/∣μ−μc∣) and letting μ→μc confines the limit of T(μ)/ln(1/∣μ−μc∣) to [1/(λu+ϵ),1/(λu−ϵ)]; since ϵ was any sufficiently small positive number, the limit is 1/λu. This is the case λu<λs, where min(λu,λs)=λu.
Step 6: the other case, by reversing time. Let λu>λs and run Steps 0-5 for the reversed field −f. Reversing time exchanges the stable and the unstable curve of the saddle, hence exchanges the two eigenvalues, so the reversed family has unstable rate λs and stable rate λu and falls under the case just treated. It has the same loop at the same μc, and a separation of the same order - though not the same separation, since reversing time changes the edge the gap is read on and the curve it is read from. That edge is now u=δ, which has become the entry edge of the square, and it is joined to s=δ only by the global map Ψ(sout)=u0+asout+O(sout2) of Step 2, with a>0; this is not a pair of sections to which the section-independence of Step 3 applies, because Ψ does not carry the crossing of one curve to the crossing of the same curve. The reversed unstable curve is the old stable curve, which meets s=δ at u=0, so the reversed separation is u~0=Ψ−1(0)=−u0/a+O(u02): the magnitude is comparable and the vanishing at μc and the non-degeneracy survive, but the sign reverses. That is why the cycle lives on the opposite side of μc in the two cases. What Step 3 supplies is then applied to u~0: it vanishes at μc and does so linearly, so ∣u~0∣ is comparable to ∣μ−μc∣, which is all Step 5 reads from it. Closed orbits and their periods are untouched by reversing time, while attraction becomes repulsion. Hence the set of closed orbits repels the one-sided neighbourhood of Γ instead of attracting it, and
ln(1/∣μ−μc∣)T(μ)⟶λs1=min(λu,λs)1.
The two cases together are the statement. What is determined is the coefficient of the logarithm, and only in the limit: the argument leaves open an amount that is o(ln(1/∣μ−μc∣)), because the rates A and B were controlled only to within ϵ. With an exactly linear flow on N the same steps give u∗=u0(μ)(1+O(u0ν−1)) and leave open a bounded amount instead, depending on the whole loop through tΓ, δ and c. □
Remark(Where the hypotheses are used, and when the table fails)AdaptedAdapted from Strogatz (1994), §8.4, Example 8.4.1 and the closing remarks, pp. 264-265; Hirsch, Smale & Devaney (2013), §16.4, pp. 377-380
Hyperbolicity and transversality in (c). The coefficient comes from the linear part at the saddle, so it needs the eigenvalues to be nonzero: at a non-hyperbolic fixed point the passage time is a power of the distance rather than a logarithm, which is exactly case (b). Transversality, u0′(μc)=0, is what makes the separation linear in μ−μc; without it the logarithm acquires a different multiple, since ln∣μ−μc∣k=kln∣μ−μc∣.
Both eigenvalues matter, and three different things depend on them. It is tempting to say that only the unstable eigenvalue appears, since it alone sets the passage time, but that is wrong: the passage time is measured from an entry point, and the entry point is set by the return map, whose exponent is the ratio ν=λs/λu. If ν>1 the return map obeys u0≤R(u)≤u0+u/2 near the saddle, which is enough to confine its fixed points to [u0,2u0] without any claim that it contracts: the cycle enters at a distance comparable to the separation of the manifolds, and the coefficient is 1/λu; if ν<1 it does not, and reading the circuit backwards in time gives 1/λs instead. So the sign of the saddle quantityλu−λs decides both which eigenvalue appears in the period and whether the cycle is stable or unstable, while the question of which side of μc carries the cycle is settled by the sign of u0′(μc), the direction in which the manifolds separate (Step 3 of the proof), together with the sign of λu−λs: for a fixed signu0′(μc) the cycle sits on one side of μc when λu<λs and on the other when λu>λs, because reversing time reverses the sign of the separation (Step 6). So the saddle quantity decides three things, this among them. The same comparison of λu with λs matters a great deal more in dimension three, where it separates a tame homoclinic bifurcation from one that creates infinitely many periodic orbits.
Genericity. The table reports the generic laws, and a symmetry or a coincidence can defeat them. Strogatz's instance is the van der Pol oscillator x¨+ϵ(x2−1)x˙+x=0: its eigenvalues cross the imaginary axis at ϵ=0, yet the cycle appears with amplitude 2 rather than O(ϵ1/2), because the nonlinear term is multiplied by the same ϵ that drives the crossing and so vanishes with it. Setting u=ϵ1/2x removes the coincidence and the amplitude of the cycle in u is 2ϵ1/2, obeying the Hopf law after all. The lesson is that a bifurcation diagnosed from a scaling law should first be checked for a degeneracy of this kind.
Dimension. All four laws persist in higher dimensions, but the list stops being exhaustive, and the homoclinic case changes character completely. In the plane, breaking a saddle loop destroys one cycle. In three dimensions a homoclinic orbit can have a return map with infinitely many fixed points, and breaking the connection unfolds infinitely many horseshoes, so that the family passes through infinitely many bifurcations as the parameter crosses μc; Hirsch, Smale and Devaney exhibit this for a perturbation of the double-scroll system. The planar analysis above is protected from any of that by the Poincare-Bendixson theorem (Limit cycles and the Poincare-Bendixson theorem).
Why local tests miss these. In (a) and (c) the Jacobian at every fixed point of the family passes through μc without doing anything: in (a) the only fixed point is the origin, a stable spiral throughout, and in (c) the saddle stays a hyperbolic saddle with eigenvalues depending smoothly on μ. An eigenvalue test run at the fixed points therefore reports nothing at all. Case (b) is the exception, and it is instructive: there a fixed point with a zero eigenvalue is born at μc, so a purely local scan does flag that parameter value - as a saddle-node bifurcation of fixed points. What the scan cannot report is that the new fixed point was born on a closed orbit and has destroyed it, nor that the period was diverging as μ↓μc; those are facts about a curve of O(1) size, invisible to any computation carried out at the fixed points alone. What changes is the global arrangement of invariant manifolds, and detecting it requires following trajectories over O(1) distances: numerically, by continuing the cycle or by tracking the unstable manifold of the saddle, which is the practical criterion Strogatz uses.
Method(Telling the four apart from measurements)AdaptedAdapted from Strogatz (1994), §8.4, pp. 260-265
Suppose an experiment or a simulation shows a stable oscillation that disappears as a control parameter p is varied, and you want to know which bifurcation did it.
Record amplitude and period as functions of p as close to the transition pc as the data allow, and plot both against μ=∣p−pc∣ on logarithmic axes. Treat pc as a fitted parameter, not a measured one: it is usually the largest source of error.
Does the amplitude vanish? If it falls like μ1/2, the bifurcation is a supercritical Hopf and is local; confirm by checking that the oscillation becomes sinusoidal and that a fixed point survives at its centre.
If the amplitude stays finite, look at the period. A finite period means a saddle-node of cycles; an O(1) jump in the oscillation's basin, and hysteresis on reversing p, confirm it. A diverging period means one of the two infinite-period cases.
Separate the two divergences. Fit T against μ−1/2 and against ln(1/μ). A power-law divergence is the infinite-period case; a logarithmic one is homoclinic. The waveform distinguishes them too: at an infinite-period bifurcation the trajectory stalls at one point of an otherwise unchanged orbit, while at a homoclinic bifurcation it stalls at a saddle that sits at a corner of the orbit.
Extract a number from the logarithm. If the divergence is logarithmic, the slope of T against ln(1/μ) is 1/min(λu,λs), the reciprocal of the smaller of the two rates at the saddle - a quantity of the model that the fit measures directly. For an oscillation seen in data the cycle is the stable one, so the slope is 1/λu and measures the unstable eigenvalue; an unstable cycle, which has to be found by continuation or by integrating backwards in time, gives 1/λs instead. Checking the measured slope against both eigenvalues of a candidate model is a sharper test of the model than the logarithm alone.
Check for degeneracy before trusting a fit, by asking whether a symmetry or a shared parameter makes a nonlinear coefficient vanish at pc, as in the van der Pol case of the remark.
Example(A homoclinic bifurcation located exactly by energy balance)Original
Consider the family
x˙=y,y˙=x−x2+ϵy(α−x),
with 0<ϵ≪1 and α the bifurcation parameter. We locate its homoclinic bifurcation in closed form.
The unperturbed system. At ϵ=0 the system is x¨=x−x2, conservative with
H(x,y)=21y2−21x2+31x3.
Its fixed points are a saddle at (0,0) (eigenvalues ±1) and a centre at (1,0), with H=0 and H=−61 respectively. The level set H=0 contains a homoclinic loop Γ to the saddle, closing at x=23, where indeed H(23,0)=−89+89=0. On Γ,
The loop survives the perturbation only if H returns to its starting value after one circuit. To leading order in ϵ we may evaluate the integral along the unperturbed loop, giving the condition
using ∫−11σ2(1−σ2)dσ=2(31−51)=154 and ∫−11σ2(1−σ2)2dσ=2(31−52+71)=10516. Hence
αc=6/536/35=76=0.857142857…
Which side carries the cycle. For α slightly larger than 6/7 the integral is positive, so to leading order in ϵ the unstable manifold returns at a higher level of H than it left and falls outside the stable manifold. Meanwhile the fixed point (1,0) has Jacobian (0−11ϵ(α−1)), with determinant 1 and trace ϵ(α−1)<0 for α<1: a stable spiral. An outward-going unstable manifold and an inward-attracting fixed point must be separated by a closed orbit, so an unstable limit cycle exists for αc(ϵ)<α<1, where αc(ϵ) is the exact bifurcation value of the family, equal to 6/7 to leading order in ϵ and located numerically below: at ϵ=0.1, for instance, αc=0.8572178, so there is no cycle in 6/7<α<0.8572178. The cycle is destroyed at α=αc(ϵ) by the homoclinic bifurcation and at α=1 by a subcritical Hopf bifurcation, where the trace changes sign. Being unstable, it is the case λu>λs of part (c), as the eigenvalues below confirm.
The predicted period scaling. At ϵ=0.1 and α=6/7 the Jacobian at the origin is (011ϵα), whose eigenvalues solve λ2−353λ−1=0; numerically λu=1.0437751 and −λs=−0.9580608. The saddle quantity is the trace, λu−λs=ϵα=3/35=0.0857143>0, so this is the case λu>λs of the statement - which is only to be expected, since the cycle here has already been identified as the unstable one. The coefficient of the logarithm is therefore
min(λu,λs)1=λs1=1.0437751,
and not1/λu=0.9580608; the two differ by 9%. They are easy to confuse in this family, because detJ=−1 at the saddle forces λuλs=1, so that 1/λs and 1/λu are the two eigenvalues over again.
Check against the full system. We located αc(ϵ) numerically by integrating the unstable manifold forward and the stable manifold backward from 10−10 along the respective eigendirections (DOP853, relative tolerance 10−13), recording the x-coordinate at which each first meets the line y=0, and solving for the α that makes the two coincide:
ϵ
αc (numerical)
αc−6/7
(αc−6/7)/ϵ2
0.10
0.8572178176
7.4960×10−5
0.0074960
0.05
0.8571615988
1.8742×10−5
0.0074967
0.02
0.8571458559
2.9987×10−6
0.0074968
0.01
0.8571436068
7.4969×10−7
0.0074969
The last two columns are formed from αc at full working precision, not from the ten decimals printed in the second: rounding αc at the tenth decimal would already disturb the fifth figure of the ratio, since at ϵ=0.01 that ratio magnifies an error in αc by ϵ−2=104. The last column is constant to three significant figures, and to four from ϵ=0.05 down, still creeping upwards in the fifth, consistent with a limit near 0.0074969, so the most that should be read off it is
αc(ϵ)=76+0.007497ϵ2+O(ϵ3).
The energy balance therefore gets the bifurcation value right with an error of order ϵ2, one order better than the O(ϵ) accuracy the derivation guarantees. We record the extra order as an observation about this family and do not account for it; the derivation above claims only the leading term.
Measuring the coefficient of the logarithm. The cycle is unstable, so we integrate the system backwards in time, where it attracts, start from (1.1,0) at α=αc(ϵ)+d, discard the transient and time the return to the section y=0, x>1 (DOP853, relative tolerance 10−13). Two members of the family are worth doing. At ϵ=0.1 the saddle sits at αc=0.8572178176, where λu=1.0437790 and λs=0.9580572, so the predicted slope of T against ln(1/d) is 1/λs=1.0437790. The second member is ϵ=1, far outside the perturbative regime that produced 6/7, but within reach of the same shooting, which gives αc=0.8645452474. That number is worth pausing on: (x,y,t)↦(x,−y,−t) turns this member of the family into x˙=y, y˙=μy+x−x2+xy with μ=−α, whose homoclinic bifurcation Strogatz (Strogatz (1994), §8.4, p. 263) locates numerically at μc≈−0.8645, so the shooting reproduces his value and we may trust it here. At this αc the saddle has λu=1.5217035 and λs=0.6571582, so the prediction is 1/λs=1.5217035 against the alternative 1/λu=0.6571582: a factor 2.32 apart. The slope column is the increment of T divided by ln10.
d
T at ϵ=0.1
slope
T at ϵ=1
slope
10−2
9.966954
-
10.250481
-
10−3
12.773439
1.218841
13.728902
1.510659
10−4
15.412528
1.146142
17.227064
1.519232
10−5
17.969326
1.110403
20.730098
1.521348
10−6
20.47992
1.09034
24.23385
1.52166
10−7
22.9616
1.07779
27.7377
1.52170
10−8
25.4239
1.06934
31.2416
1.5217
The last three rows are quoted to fewer figures than the others, because at small d the measured time inherits the uncertainty of αc: an error δ in αc is an error δ in the true distance from the bifurcation, and so moves T by about δ/(λsd). The shooting pins αc down to some 10−12, which at d=10−8 is 1.5×10−4 in T and 6×10−5 in the slope - the last digit printed in that row may be off by a unit or two - while at d≥10−5 it is below 10−6 and all six decimals stand. The ϵ=1 column settles on 1.5217, which is 1/λs=1.5217035 to the five figures the shooting supports, and is nowhere near 1/λu=0.6571582. The ϵ=0.1 column converges much more slowly - the two rates there are nearly equal, and the relative error in the cycle's entry point decays only like dλu/λs−1=d0.0895 - but it falls monotonically towards 1.0437790, its distance from that value shrinking by a factor of roughly 0.7 per decade of d across the range shown, and every entry in it stays well above the wrong candidate 0.9580572. The reciprocal of the smaller eigenvalue is what the data measure.
Example(Reading a period off an infinite-period bifurcation)Original
Take prototype (b) at μ=1.01, a distance 0.01 past the bifurcation. The exact period is
against 2π/μ=6.2209756, the time the same circle would take at the uniform angular speed μ: the bottleneck slows the orbit down by a factor 7.124, although its shape is exactly the same unit circle. The square-root law predicts
μ−1π2=0.14.4428829=44.428829,
which overestimates the exact value by 0.25% - the relative error being O(μ−1)=O(10−2), as the expansion in the proof says. Halving the distance, to μ=1.005, gives T=62.753460, a factor 1.41598 against 2=1.414214; that factor is the signature to look for in data: a period that doubles when the distance from the transition is divided by four, with the waveform developing a flat stretch but keeping its amplitude.
Application(Hysteresis in a driven pendulum, and two kinds of neuron)AdaptedAdapted from Strogatz (1994), §8.5, pp. 265-272
A pendulum driven by a constant torque, and equally a Josephson junction driven by a constant current, obeys
ϕ′′+αϕ′+sinϕ=I
on the cylinder, with α>0 the dimensionless damping and I the dimensionless drive. For I<1 there are two fixed points, a saddle and a sink - the pendulum can hang at rest with gravity balancing the torque - and they merge in a saddle-node at I=1. For I>1 no rest state exists and there is instead a unique stable cycle on which the pendulum whirls over the top. The interesting regime is weak damping, where the whirling cycle and the rest state coexist over a range of I, so that the state reached depends on the history of the drive: raising I past one threshold makes the junction switch to the voltage-carrying whirling state, and lowering it again does not switch it back until a different, lower threshold. Only the lower threshold belongs to this unit. The upper one is the saddle-node bifurcation of fixed points at I=1, which destroys the rest state and which Strogatz lists separately from the two bifurcations of cycles. Which bifurcation of cycles supplies the lower threshold is decided by the damping. At strong damping the cycle is destroyed as soon as a fixed point appears on it, which happens at I=1 itself: an infinite-period bifurcation, and no hysteresis at all, because the two thresholds then coincide. At weak damping the cycle instead survives below I=1 and is destroyed later, at a lower drive Ic<1, by running into the saddle: a homoclinic bifurcation, with the junction bistable throughout Ic<I<1. Strogatz (Strogatz (1994), §8.5, p. 272) reports from Guckenheimer and Holmes, who obtain it by Melnikov's method, that the curve Ic(α) is tangent to the line I=4α/π as α→0. So the bistable range is at its widest at the weakest damping, where it is nearly the whole of 0<I<1, and it narrows as α grows until the homoclinic curve reaches I=1 and the infinite-period scenario takes over. That width is the usable switching margin of the device.
The same distinction organises a textbook fact about nerve cells. Strogatz's section is about the pendulum and the junction and says nothing about neurons; what follows is our own reading of the scaling table, and we give it because it is the clearest everyday instance of the table being used as a diagnostic. A neuron that begins to fire at an arbitrarily low rate as the injected current is raised is one whose resting state loses stability through a fixed point disappearing from a cycle: the period diverges at threshold, so arbitrarily long interspike intervals are available. A neuron that begins to fire abruptly at a nonzero rate is one whose rest state loses stability through a Hopf bifurcation, where the period stays O(1) and only the amplitude grows. Measuring the firing rate just above threshold therefore identifies the bifurcation, and with it which reduced model can reproduce the cell - exactly the inference the scaling table is for.