Ordinary Differential Equations (7): Stability Theory
Chen Kai BOSS
  2019-05-07 14:00:00 2019-05-07 14:00      5.8k Words  36 Mins  

What happens to a system after a small perturbation? Does it return to its original state, or drift further and further away? The answer to this question determines whether a bridge will collapse in the wind, whether an ecosystem can maintain balance, and whether economic cycles will spiral out of control. Stability theory gives us the mathematical tools to answer these questions. In this chapter, we'll start from intuition and gradually build Lyapunov stability theory, seeing its wide applications in engineering, biology, and physics.

Starting with an Inverted Pendulum

Imagine a vertical stick that you're balancing on your fingertip.

Two equilibrium states: 1. Stick pointing down (center of mass below support point): Give it a slight push, it will swing a few times and return to vertical — this is stable 2. Stick pointing up (center of mass above support point): The slightest deviation and it will fall — this is unstable

Mathematically, both equilibrium points satisfy the same condition (derivative equals zero), but their nature is completely different. The task of stability theory is to distinguish these cases.

Precise Definitions of Stability

Lyapunov Stability

Consider an autonomous system: Letbe an equilibrium point, i.e.,.

Definition (Lyapunov stable): The equilibrium pointis stable if for any, there existssuch that when, we havefor all.

Intuition: Starting near the equilibrium, the trajectory never goes too far.

Definition (Asymptotically stable): Ifis stable, and there existssuch that when, then.

Intuition: Not only doesn't go far, but eventually returns to the equilibrium.

Definition (Unstable): Ifis not stable, it is called unstable.

Basin of Attraction

Definition: The basin of attraction of equilibrium pointis the set of all initial valuessatisfying.

The basin of attraction can be the entire space (globally asymptotically stable), or just a region around the equilibrium (locally asymptotically stable).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Visualizing different types of stability
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# System 1: Stable but not asymptotically stable (harmonic oscillator)
def harmonic(x, t):
    return [x[1], -x[0]]

ax1 = axes[0]
t = np.linspace(0, 20, 1000)
for r in [0.5, 1, 1.5, 2]:
    for theta0 in np.linspace(0, 2*np.pi, 4, endpoint=False):
        x0 = [r*np.cos(theta0), r*np.sin(theta0)]
        sol = odeint(harmonic, x0, t)
        ax1.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.5, alpha=0.7)

ax1.plot(0, 0, 'ko', markersize=8)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('Stable (not asymptotically)\nCenter', fontsize=12, fontweight='bold')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# System 2: Asymptotically stable (damped oscillator)
def damped(x, t):
    return [x[1], -x[0] - 0.3*x[1]]

ax2 = axes[1]
for r in [0.5, 1, 1.5, 2]:
    for theta0 in np.linspace(0, 2*np.pi, 4, endpoint=False):
        x0 = [r*np.cos(theta0), r*np.sin(theta0)]
        sol = odeint(damped, x0, t)
        ax2.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.8, alpha=0.7)

ax2.plot(0, 0, 'ko', markersize=8)
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('y', fontsize=12)
ax2.set_title('Asymptotically Stable\nStable Spiral', fontsize=12, fontweight='bold')
ax2.set_aspect('equal')
ax2.grid(True, alpha=0.3)

# System 3: Unstable (saddle point)
def saddle(x, t):
    return [x[0], -x[1]]

ax3 = axes[2]
t_short = np.linspace(0, 2, 500)
t_back = np.linspace(0, -2, 500)

initial_conditions = [
    [0.1, 0.5], [0.1, -0.5], [-0.1, 0.5], [-0.1, -0.5],
    [0.5, 0.1], [0.5, -0.1], [-0.5, 0.1], [-0.5, -0.1]
]

for x0 in initial_conditions:
    sol = odeint(saddle, x0, t_short)
    sol_back = odeint(saddle, x0, t_back)
    ax3.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.8, alpha=0.7)
    ax3.plot(sol_back[:, 0], sol_back[:, 1], 'b-', linewidth=0.8, alpha=0.7)

# Stable and unstable manifolds
ax3.arrow(0, 0, 0, 1.5, head_width=0.1, head_length=0.1, fc='green', ec='green', linewidth=2)
ax3.arrow(0, 0, 0, -1.5, head_width=0.1, head_length=0.1, fc='green', ec='green', linewidth=2)
ax3.arrow(0, 0, 1.5, 0, head_width=0.1, head_length=0.1, fc='red', ec='red', linewidth=2)
ax3.arrow(0, 0, -1.5, 0, head_width=0.1, head_length=0.1, fc='red', ec='red', linewidth=2)

ax3.plot(0, 0, 'ko', markersize=8)
ax3.set_xlabel('x', fontsize=12)
ax3.set_ylabel('y', fontsize=12)
ax3.set_title('Unstable\nSaddle Point', fontsize=12, fontweight='bold')
ax3.text(1.2, 0.2, 'Unstable', fontsize=10, color='red')
ax3.text(0.1, 1.2, 'Stable', fontsize=10, color='green')
ax3.set_aspect('equal')
ax3.grid(True, alpha=0.3)
ax3.set_xlim(-2, 2)
ax3.set_ylim(-2, 2)

plt.tight_layout()
plt.savefig('stability_types.png', dpi=150, bbox_inches='tight')
plt.show()

Linearization Method

Basic Idea

For a nonlinear system, perform Taylor expansion near equilibrium point:Since, let:whereis the Jacobian matrix.

Jacobian Matrix

Hartman-Grobman Theorem

Theorem: If all eigenvalues ofhave nonzero real parts (hyperbolic equilibrium point), then the qualitative behavior of the nonlinear system nearis the same as the linearized system.

Precise meaning of "same qualitative behavior": There exists a homeomorphism (continuous bijection with continuous inverse) that maps trajectories of the nonlinear system to trajectories of the linear system, preserving time direction.

Important corollaries: - If all eigenvalues ofhave negative real parts → origin is asymptotically stable - Ifhas any eigenvalue with positive real part → origin is unstable

Warning: Ifhas purely imaginary eigenvalues, the linearization method fails! Higher-order analysis is needed.

Example: Damped Pendulum

Pendulum equation (without small angle approximation):Writing as a first-order system (,):$ $

Equilibrium points:and, i.e.,, Jacobian matrix:At:Eigenvalues:If, both eigenvalues have negative real parts → asymptotically stable (hanging position)

At:Eigenvalues:One positive, one negative → saddle point → unstable (inverted position)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Damped pendulum
gamma = 0.3  # Damping coefficient
omega0 = 1.0  # Natural frequency

def pendulum(X, t):
    x, y = X
    dxdt = y
    dydt = -gamma * y - omega0**2 * np.sin(x)
    return [dxdt, dydt]

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Phase space
ax1 = axes[0]
t = np.linspace(0, 50, 2000)

# Trajectories from different initial values
for x0 in np.linspace(-3*np.pi, 3*np.pi, 15):
    for y0 in np.linspace(-3, 3, 5):
        X0 = [x0, y0]
        sol = odeint(pendulum, X0, t)
        ax1.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.3, alpha=0.5)

# Mark equilibrium points
for n in range(-2, 3):
    if n % 2 == 0:  # Stable equilibrium
        ax1.plot(n*np.pi, 0, 'go', markersize=10, label='Stable' if n == 0 else '')
    else:  # Unstable equilibrium (saddle)
        ax1.plot(n*np.pi, 0, 'rx', markersize=10, mew=2, label='Unstable' if n == 1 else '')

ax1.set_xlabel('θ (angle)', fontsize=12)
ax1.set_ylabel('ω (angular velocity)', fontsize=12)
ax1.set_title('Phase Portrait of Damped Pendulum', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-3*np.pi, 3*np.pi)
ax1.set_ylim(-4, 4)

# Time evolution
ax2 = axes[1]
initial_conditions = [
    [0.5, 0, 'Near stable eq'],
    [np.pi - 0.3, 0, 'Near unstable eq'],
    [0, 2, 'Large velocity']
]

for x0, y0, label in initial_conditions:
    sol = odeint(pendulum, [x0, y0], t)
    ax2.plot(t, sol[:, 0], linewidth=2, label=label)

ax2.axhline(y=0, color='g', linestyle='--', alpha=0.5)
ax2.axhline(y=np.pi, color='r', linestyle='--', alpha=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('θ (angle)', fontsize=12)
ax2.set_title('Time Evolution', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('damped_pendulum.png', dpi=150, bbox_inches='tight')
plt.show()

# Calculate eigenvalues
print("Equilibrium point analysis:")
for eq_name, x_eq in [("Hanging position (0, 0)", 0), ("Inverted position (π, 0)", np.pi)]:
    A = np.array([[0, 1], [-omega0**2 * np.cos(x_eq), -gamma]])
    eigenvalues = np.linalg.eigvals(A)
    print(f"\n{eq_name}:")
    print(f"  Jacobian matrix:\n{A}")
    print(f"  Eigenvalues: {eigenvalues}")
    if all(ev.real < 0 for ev in eigenvalues):
        print(f"  Conclusion: Asymptotically stable")
    elif any(ev.real > 0 for ev in eigenvalues):
        print(f"  Conclusion: Unstable")

Lyapunov's Direct Method

Inspiration from Energy

The stability of physical systems is usually related to energy: - States with minimum energy are stable - If energy keeps decreasing, the system will tend toward a stable state

Lyapunov's direct method mathematizes this idea.

Lyapunov Functions

Definition: For system, if there exists a functionsatisfying: 1.$V(^) = 0V() > 0 ^ = V V$is called a Lyapunov function.

Stability Theorems

Theorem (Lyapunov stability theorem): 1. If a Lyapunov function exists with, thenis stable 2. Iffor all, thenis asymptotically stable 3. If there exists a functionwith, thenis unstable

Intuition:is like "generalized energy." If it always decreases, the system will eventually reach the lowest energy point (equilibrium).

How to Construct Lyapunov Functions?

This is an art; there's no universal formula. But there are some commonly used techniques:

Technique 1: Physical energy

For mechanical systems, total energy(kinetic + potential) is often a good candidate.

Technique 2: Quadratic forms

For linear systems, try, whereis a positive definite matrix.If we can findsuch that(positive definite), then→ asymptotically stable.

This is called the Lyapunov equation.

Technique 3: Trial and error

Start from simple forms (like), check the sign of.

Example: Stability of Van der Pol Oscillator

Van der Pol equation:This describes an electronic tube oscillator with nonlinear damping.

Writing as a first-order system:$ $

Equilibrium point: Jacobian matrix:Eigenvalues:When,→ unstable!

But there's an interesting phenomenon: although the origin is unstable, the system has a stable limit cycle.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Van der Pol oscillator
def van_der_pol(X, t, mu):
    x, y = X
    dxdt = y
    dydt = mu * (1 - x**2) * y - x
    return [dxdt, dydt]

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

t = np.linspace(0, 50, 2000)

# Phase portraits for different μ values
mu_values = [0.5, 1.0, 3.0]

for ax, mu in zip(axes, mu_values):
    # Starting from inside
    for r in [0.1, 0.3, 0.5]:
        X0 = [r, 0]
        sol = odeint(van_der_pol, X0, t, args=(mu,))
        ax.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.8, alpha=0.7)
    
    # Starting from outside
    for r in [3, 4, 5]:
        X0 = [r, 0]
        sol = odeint(van_der_pol, X0, t, args=(mu,))
        ax.plot(sol[:, 0], sol[:, 1], 'r-', linewidth=0.8, alpha=0.7)
    
    ax.plot(0, 0, 'ko', markersize=8)
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title(f'Van der Pol Oscillator (μ = {mu})', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')

plt.tight_layout()
plt.savefig('van_der_pol.png', dpi=150, bbox_inches='tight')
plt.show()

# Time evolution showing limit cycle
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

mu = 1.0
t_long = np.linspace(0, 100, 5000)

# Starting from inside
ax1 = axes[0]
X0_in = [0.1, 0]
sol_in = odeint(van_der_pol, X0_in, t_long, args=(mu,))
ax1.plot(t_long, sol_in[:, 0], 'b-', linewidth=1)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('x(t)', fontsize=12)
ax1.set_title('Starting Inside Limit Cycle', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# Starting from outside
ax2 = axes[1]
X0_out = [4, 0]
sol_out = odeint(van_der_pol, X0_out, t_long, args=(mu,))
ax2.plot(t_long, sol_out[:, 0], 'r-', linewidth=1)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('x(t)', fontsize=12)
ax2.set_title('Starting Outside Limit Cycle', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('van_der_pol_time.png', dpi=150, bbox_inches='tight')
plt.show()

LaSalle's Invariance Principle

The Problem

Whenbut not strictly less than zero, the standard Lyapunov theorem can only guarantee stability, not asymptotic stability.

LaSalle's Invariance Principle

Theorem: Letbe a Lyapunov function with. Define the setExtra close brace or missing open braceE = \{\mathbf{x} : \dot{V}(\mathbf{x}) = 0}.

Letbe the largest invariant set contained in(i.e., if a trajectory starts in, it stays inforever).

Then all bounded trajectories approach.

Application: Ifcontains only the equilibrium point, thenis asymptotically stable.

Example

Consider a pendulum with friction:Take the energy function:Compute: Extra close brace or missing open braceE = \{\dot{\theta} = 0}.

On ,is constant, so. Substituting into the equation:→.

For, this is indeed an invariant set in. Sois asymptotically stable.

Global Stability

Global Asymptotic Stability

Definition: If the basin of attraction of equilibrium pointis the entire space, it is called globally asymptotically stable.

Theorem: If: 1.for all$ V() || < 0 ^*$is globally asymptotically stable.

Example: Logistic EquationEquilibrium points:and.

For, take.Whenand,. So on,is asymptotically stable.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Logistic equation
r = 1.0  # Growth rate
K = 100  # Carrying capacity

def logistic(N, t):
    return r * N * (1 - N/K)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

t = np.linspace(0, 15, 500)

# Solutions from different initial values
ax1 = axes[0]
N0_values = [10, 30, 50, 80, 120, 150, 200]

for N0 in N0_values:
    sol = odeint(logistic, N0, t)
    ax1.plot(t, sol, linewidth=2, label=f'N ₀ = {N0}')

ax1.axhline(y=K, color='r', linestyle='--', linewidth=2, label=f'K = {K}')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Population N', fontsize=12)
ax1.set_title('Logistic Growth: Global Stability of K', fontsize=14, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# Lyapunov function contours
ax2 = axes[1]
N_range = np.linspace(1, 200, 200)
V = 0.5 * (N_range - K)**2

ax2.plot(N_range, V, 'b-', linewidth=2)
ax2.axvline(x=K, color='r', linestyle='--', linewidth=2, label=f'N = K = {K}')
ax2.set_xlabel('Population N', fontsize=12)
ax2.set_ylabel('V(N) = ½(N - K)²', fontsize=12)
ax2.set_title('Lyapunov Function', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# Plot dV/dt on the same figure
ax2_twin = ax2.twinx()
dV = -(r/K) * N_range * (N_range - K)**2
ax2_twin.plot(N_range, dV, 'g-', linewidth=2, alpha=0.7)
ax2_twin.set_ylabel('$\\dot{V}$', fontsize=12, color='g')
ax2_twin.tick_params(axis='y', labelcolor='g')

plt.tight_layout()
plt.savefig('logistic_stability.png', dpi=150, bbox_inches='tight')
plt.show()

Lyapunov Analysis of Linear Systems

Lyapunov Equation

For a linear system, let.

Theorem: If all eigenvalues ofhave negative real parts, then for any positive definite matrix, there exists a unique positive definite matrixsatisfying:This is called the continuous Lyapunov equation.

Corollary: If we can findsuch thatis negative definite, the system is asymptotically stable.

Numerical Solution

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import numpy as np
from scipy.linalg import solve_continuous_lyapunov
import matplotlib.pyplot as plt

# Lyapunov analysis of linear system
A = np.array([[-2, 1], 
              [-1, -3]])

# Check stability (eigenvalues)
eigenvalues = np.linalg.eigvals(A)
print("Eigenvalues:", eigenvalues)
print("All eigenvalues have Re < 0:", all(ev.real < 0 for ev in eigenvalues))

# Solve Lyapunov equation A'P + PA = -Q
Q = np.eye(2)  # Take identity matrix
P = solve_continuous_lyapunov(A.T, -Q)

print("\nQ =")
print(Q)
print("\nP =")
print(P)

# Verify P is positive definite
eigenvalues_P = np.linalg.eigvals(P)
print("\nEigenvalues of P:", eigenvalues_P)
print("P is positive definite:", all(ev > 0 for ev in eigenvalues_P))

# Visualize Lyapunov function
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Contour plot
ax1 = axes[0]
x1 = np.linspace(-2, 2, 100)
x2 = np.linspace(-2, 2, 100)
X1, X2 = np.meshgrid(x1, x2)

# V(x) = x'Px
V = P[0,0]*X1**2 + (P[0,1]+P[1,0])*X1*X2 + P[1,1]*X2**2

contour = ax1.contour(X1, X2, V, levels=20, cmap='viridis')
ax1.clabel(contour, inline=True, fontsize=8)

# Plot trajectories
from scipy.integrate import odeint

def linear_system(x, t):
    return A @ x

t = np.linspace(0, 5, 500)
for angle in np.linspace(0, 2*np.pi, 8, endpoint=False):
    x0 = [1.5*np.cos(angle), 1.5*np.sin(angle)]
    sol = odeint(linear_system, x0, t)
    ax1.plot(sol[:, 0], sol[:, 1], 'r-', linewidth=1, alpha=0.7)

ax1.plot(0, 0, 'ko', markersize=8)
ax1.set_xlabel('$x_1$', fontsize=12)
ax1.set_ylabel('$x_2$', fontsize=12)
ax1.set_title('Lyapunov Function V(x) = x\'Px\nwith Trajectories', fontsize=14, fontweight='bold')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# V along trajectory
ax2 = axes[1]
x0 = [1.5, 1.0]
sol = odeint(linear_system, x0, t)

V_traj = []
for state in sol:
    V_traj.append(state @ P @ state)

ax2.plot(t, V_traj, 'b-', linewidth=2)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('V(x(t))', fontsize=12)
ax2.set_title('Lyapunov Function Decreases Along Trajectory', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lyapunov_linear.png', dpi=150, bbox_inches='tight')
plt.show()

Classification of Two-Dimensional Systems

Complete Classification

For two-dimensional linear systems, classification based on eigenvalues:

Type Eigenvalue Condition Stability Phase Portrait
Stable node (real) Asymptotically stable Trajectories approach origin
Unstable node (real) Unstable Trajectories leave origin
Saddle point (real) Unstable Saddle-type trajectories
Stable spiral ,(complex) Asymptotically stable Spiral toward origin
Unstable spiral ,(complex) Unstable Spiral away from origin
Center ,(purely imaginary) Stable (not asymptotic) Closed orbits
Degenerate node , lacks eigenvectors Depends on sign Star/degenerate

Classification Using Trace and Determinant

Let,.

Discriminant: -: Saddle point -,,: Stable node -,,: Unstable node -,,: Stable spiral -,,: Unstable spiral -,: Center -: Non-isolated equilibria (a line)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
import numpy as np
import matplotlib.pyplot as plt

# Trace-determinant plane classification diagram
fig, ax = plt.subplots(1, 1, figsize=(10, 8))

tau = np.linspace(-4, 4, 400)
delta = np.linspace(-4, 4, 400)
Tau, Delta = np.meshgrid(tau, delta)

# Boundary lines
# 1. Delta = 0 (x-axis)
ax.axhline(y=0, color='black', linewidth=2)

# 2. tau = 0 (y-axis, in region Delta > 0)
ax.axvline(x=0, color='black', linewidth=2)

# 3. Discriminant D = tau^2 - 4*Delta = 0, i.e., Delta = tau^2/4
tau_curve = np.linspace(-4, 4, 200)
delta_curve = tau_curve**2 / 4
ax.plot(tau_curve, delta_curve, 'k-', linewidth=2)

# Fill regions
# Stable node: Delta > 0, tau < 0, tau^2 > 4*Delta
ax.fill_between(tau[tau < 0], 0, tau[tau < 0]**2/4, 
                alpha=0.3, color='blue', label='Stable Node')

# Stable spiral: Delta > tau^2/4, tau < 0
tau_neg = tau[tau < 0]
ax.fill_between(tau_neg, tau_neg**2/4, 4, 
                alpha=0.3, color='cyan', label='Stable Spiral')

# Unstable node: Delta > 0, tau > 0, tau^2 > 4*Delta
ax.fill_between(tau[tau > 0], 0, tau[tau > 0]**2/4, 
                alpha=0.3, color='red', label='Unstable Node')

# Unstable spiral: Delta > tau^2/4, tau > 0
tau_pos = tau[tau > 0]
ax.fill_between(tau_pos, tau_pos**2/4, 4, 
                alpha=0.3, color='orange', label='Unstable Spiral')

# Saddle: Delta < 0
ax.fill_between(tau, -4, 0, alpha=0.3, color='yellow', label='Saddle')

# Center (tau = 0, Delta > 0) marked with line
ax.plot([0, 0], [0, 4], 'g-', linewidth=3, label='Center (τ=0, Δ>0)')

# Labels
ax.text(-2.5, 2, 'Stable\nSpiral', fontsize=12, ha='center', fontweight='bold')
ax.text(-2.5, 0.3, 'Stable\nNode', fontsize=12, ha='center', fontweight='bold')
ax.text(2.5, 2, 'Unstable\nSpiral', fontsize=12, ha='center', fontweight='bold')
ax.text(2.5, 0.3, 'Unstable\nNode', fontsize=12, ha='center', fontweight='bold')
ax.text(0, -2, 'Saddle', fontsize=12, ha='center', fontweight='bold')

ax.set_xlabel('τ = tr(A)', fontsize=14)
ax.set_ylabel('Δ = det(A)', fontsize=14)
ax.set_title('Classification of 2D Linear Systems', fontsize=16, fontweight='bold')
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.legend(loc='upper left', fontsize=10)
ax.grid(True, alpha=0.3)

# Add mathematical relation
ax.text(2.8, 3.5, '$D = \\tau^2 - 4\\Delta = 0$', fontsize=10)
ax.annotate('', xy=(2.5, 1.56), xytext=(2.8, 3.3),
            arrowprops=dict(arrowstyle='->', color='black'))

plt.tight_layout()
plt.savefig('trace_det_classification.png', dpi=150, bbox_inches='tight')
plt.show()

Introduction to Bifurcation Theory

What is Bifurcation?

When system parameters change, the number or stability of equilibrium points may undergo sudden changes. Such changes are called bifurcations.

Saddle-Node Bifurcation

Consider a one-dimensional system: - When: Two equilibrium points - When: One equilibrium point(non-hyperbolic) - When: No equilibrium points

At, two equilibrium points "collide and disappear."

Transcritical BifurcationEquilibrium points:and

  • When:stable,unstable
  • When:unstable,stable

At, two equilibrium points "exchange stability."

Pitchfork Bifurcation

Supercritical pitchfork bifurcation: - When: One stable equilibrium point - When:becomes unstable, two new stable equilibrium pointsappear

Subcritical pitchfork bifurcation:Bifurcation direction is reversed, may lead to "jump" behavior.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Three types of bifurcation
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Saddle-node bifurcation
ax1 = axes[0]
r_range = np.linspace(-2, 2, 200)

# Plot equilibrium points
r_neg = r_range[r_range < 0]
x_eq_pos = np.sqrt(-r_neg)
x_eq_neg = -np.sqrt(-r_neg)

ax1.plot(r_neg, x_eq_pos, 'b-', linewidth=2, label='Stable')
ax1.plot(r_neg, x_eq_neg, 'r--', linewidth=2, label='Unstable')
ax1.plot(0, 0, 'ko', markersize=8)

ax1.set_xlabel('r', fontsize=12)
ax1.set_ylabel('x*', fontsize=12)
ax1.set_title('Saddle-Node Bifurcation\n$\\dot{x} = r + x^2$', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-2, 2)
ax1.set_ylim(-2, 2)
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.axvline(x=0, color='k', linewidth=0.5)

# Transcritical bifurcation
ax2 = axes[1]
ax2.plot(r_range, np.zeros_like(r_range), 'b-', linewidth=2, label='x* = 0')
ax2.plot(r_range, r_range, 'r-', linewidth=2, label='x* = r')

# Stability markers
ax2.plot(r_range[r_range < 0], np.zeros_like(r_range[r_range < 0]), 'b-', linewidth=3)
ax2.plot(r_range[r_range > 0], np.zeros_like(r_range[r_range > 0]), 'b--', linewidth=3)
ax2.plot(r_range[r_range < 0], r_range[r_range < 0], 'r--', linewidth=3)
ax2.plot(r_range[r_range > 0], r_range[r_range > 0], 'r-', linewidth=3)

ax2.plot(0, 0, 'ko', markersize=8)

ax2.set_xlabel('r', fontsize=12)
ax2.set_ylabel('x*', fontsize=12)
ax2.set_title('Transcritical Bifurcation\n$\\dot{x} = rx - x^2$', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-2, 2)
ax2.set_ylim(-2, 2)
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=0, color='k', linewidth=0.5)
ax2.text(-1.5, -1.5, 'Unstable', fontsize=10, color='red')
ax2.text(1.5, 1.5, 'Stable', fontsize=10, color='red')
ax2.text(-1.5, 0.2, 'Stable', fontsize=10, color='blue')
ax2.text(1.5, 0.2, 'Unstable', fontsize=10, color='blue')

# Pitchfork bifurcation (supercritical)
ax3 = axes[2]
r_pos = r_range[r_range >= 0]

ax3.plot(r_range, np.zeros_like(r_range), 'k-', linewidth=2)
ax3.plot(r_range[r_range <= 0], np.zeros_like(r_range[r_range <= 0]), 'b-', linewidth=3)
ax3.plot(r_range[r_range >= 0], np.zeros_like(r_range[r_range >= 0]), 'r--', linewidth=3)

# New stable branches
x_branch = np.sqrt(r_pos)
ax3.plot(r_pos, x_branch, 'b-', linewidth=3)
ax3.plot(r_pos, -x_branch, 'b-', linewidth=3)

ax3.plot(0, 0, 'ko', markersize=8)

ax3.set_xlabel('r', fontsize=12)
ax3.set_ylabel('x*', fontsize=12)
ax3.set_title('Supercritical Pitchfork\n$\\dot{x} = rx - x^3$', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.set_xlim(-2, 2)
ax3.set_ylim(-2, 2)
ax3.axhline(y=0, color='k', linewidth=0.5)
ax3.axvline(x=0, color='k', linewidth=0.5)
ax3.text(-1, 0.2, 'Stable', fontsize=10, color='blue')
ax3.text(1, 0.2, 'Unstable', fontsize=10, color='red')
ax3.text(1.5, 1.2, 'Stable', fontsize=10, color='blue')

plt.tight_layout()
plt.savefig('bifurcations.png', dpi=150, bbox_inches='tight')
plt.show()

Hopf Bifurcation

Birth of Limit Cycles

When the real part of complex eigenvalues changes from negative to positive, a Hopf bifurcation occurs: the equilibrium changes from a stable spiral to an unstable spiral, and a limit cycle appears.

Supercritical Hopf bifurcation: The limit cycle is stable Subcritical Hopf bifurcation: The limit cycle is unstable

Example$

In polar coordinates:$

- When: The origin is a stable spiral - When: The origin is an unstable spiral, a stable limit cycle exists at 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Hopf bifurcation example
def hopf_system(X, t, mu):
    x, y = X
    r2 = x**2 + y**2
    dxdt = mu*x - y - x*r2
    dydt = x + mu*y - y*r2
    return [dxdt, dydt]

fig, axes = plt.subplots(2, 3, figsize=(15, 10))

mu_values = [-0.5, 0, 0.5]
t = np.linspace(0, 50, 2000)

for i, mu in enumerate(mu_values):
    ax1 = axes[0, i]
    ax2 = axes[1, i]
    
    # Phase portrait
    # Starting from different initial values
    for r0 in [0.3, 0.6, 1.0, 1.5]:
        for theta0 in [0, np.pi/2, np.pi, 3*np.pi/2]:
            x0 = [r0*np.cos(theta0), r0*np.sin(theta0)]
            sol = odeint(hopf_system, x0, t, args=(mu,))
            ax1.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.5, alpha=0.7)
    
    # Draw limit cycle (if exists)
    if mu > 0:
        theta = np.linspace(0, 2*np.pi, 100)
        r_lim = np.sqrt(mu)
        ax1.plot(r_lim*np.cos(theta), r_lim*np.sin(theta), 'r-', linewidth=2, label='Limit cycle')
    
    ax1.plot(0, 0, 'ko', markersize=8)
    ax1.set_xlabel('x', fontsize=12)
    ax1.set_ylabel('y', fontsize=12)
    ax1.set_title(f'μ = {mu}', fontsize=12, fontweight='bold')
    ax1.set_aspect('equal')
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(-2, 2)
    ax1.set_ylim(-2, 2)
    if mu > 0:
        ax1.legend(fontsize=10)
    
    # Time evolution
    x0 = [1.5, 0]
    sol = odeint(hopf_system, x0, t, args=(mu,))
    r_traj = np.sqrt(sol[:, 0]**2 + sol[:, 1]**2)
    
    ax2.plot(t, r_traj, 'b-', linewidth=1.5)
    if mu > 0:
        ax2.axhline(y=np.sqrt(mu), color='r', linestyle='--', linewidth=2, label=f'$r = \\sqrt{{ \\mu }} = {np.sqrt(mu):.2f}$')
        ax2.legend(fontsize=10)
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('r(t)', fontsize=12)
    ax2.set_title(f'Radius vs Time (μ = {mu})', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)

plt.suptitle('Supercritical Hopf Bifurcation', fontsize=16, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('hopf_bifurcation.png', dpi=150, bbox_inches='tight')
plt.show()

# Bifurcation diagram
fig, ax = plt.subplots(figsize=(10, 6))

mu_range = np.linspace(-1, 1, 200)

# Origin stability
ax.plot(mu_range[mu_range <= 0], np.zeros(np.sum(mu_range <= 0)), 'b-', linewidth=3, label='Stable eq.')
ax.plot(mu_range[mu_range >= 0], np.zeros(np.sum(mu_range >= 0)), 'r--', linewidth=3, label='Unstable eq.')

# Limit cycle amplitude
mu_pos = mu_range[mu_range > 0]
ax.plot(mu_pos, np.sqrt(mu_pos), 'g-', linewidth=3, label='Stable limit cycle')
ax.plot(mu_pos, -np.sqrt(mu_pos), 'g-', linewidth=3)

ax.axvline(x=0, color='k', linestyle=':', linewidth=1)
ax.set_xlabel('μ (bifurcation parameter)', fontsize=14)
ax.set_ylabel('Amplitude', fontsize=14)
ax.set_title('Hopf Bifurcation Diagram', fontsize=16, fontweight='bold')
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('hopf_diagram.png', dpi=150, bbox_inches='tight')
plt.show()

Application: Stability of Ecosystems

Lotka-Volterra Predator-Prey Model$

Equilibrium points: 1.: Both species extinct 2.: Coexistence

Jacobian matrix at the coexistence equilibrium:Eigenvalues:(purely imaginary)

Conclusion: The coexistence equilibrium is a center (non-hyperbolic), linearization fails. In fact, the system has infinitely many closed orbits, with populations oscillating periodically.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Lotka-Volterra model
a = 1.0  # Prey growth rate
b = 0.5  # Predation rate
c = 0.5  # Predator death rate
d = 0.2  # Conversion efficiency

def lotka_volterra(X, t):
    x, y = X
    dxdt = a*x - b*x*y
    dydt = -c*y + d*x*y
    return [dxdt, dydt]

# Equilibrium point
x_eq = c/d
y_eq = a/b

print(f"Coexistence equilibrium: ({x_eq:.2f}, {y_eq:.2f})")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Phase space
ax1 = axes[0]
t = np.linspace(0, 50, 2000)

# Trajectories from different initial values
initial_conditions = [
    [1, 1], [2, 1], [3, 1], [4, 1], [5, 1],
    [1, 2], [1, 3], [2, 3], [3, 2]
]

for x0 in initial_conditions:
    sol = odeint(lotka_volterra, x0, t)
    ax1.plot(sol[:, 0], sol[:, 1], linewidth=1, alpha=0.7)

ax1.plot(x_eq, y_eq, 'ro', markersize=10, label=f'Equilibrium ({x_eq:.1f}, {y_eq:.1f})')
ax1.plot(0, 0, 'ko', markersize=8)
ax1.set_xlabel('Prey x', fontsize=12)
ax1.set_ylabel('Predator y', fontsize=12)
ax1.set_title('Lotka-Volterra: Closed Orbits (Center)', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Time evolution
ax2 = axes[1]
x0 = [4, 1]
sol = odeint(lotka_volterra, x0, t)

ax2.plot(t, sol[:, 0], 'b-', linewidth=2, label='Prey x(t)')
ax2.plot(t, sol[:, 1], 'r-', linewidth=2, label='Predator y(t)')
ax2.axhline(y=x_eq, color='b', linestyle='--', alpha=0.5)
ax2.axhline(y=y_eq, color='r', linestyle='--', alpha=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Population', fontsize=12)
ax2.set_title('Population Dynamics: Periodic Oscillations', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lotka_volterra.png', dpi=150, bbox_inches='tight')
plt.show()

# Energy function (conserved quantity)
# V = d*x - c*ln(x) + b*y - a*ln(y) = const
fig, ax = plt.subplots(figsize=(10, 8))

x_range = np.linspace(0.1, 8, 200)
y_range = np.linspace(0.1, 6, 200)
X, Y = np.meshgrid(x_range, y_range)

V = d*X - c*np.log(X) + b*Y - a*np.log(Y)

contour = ax.contour(X, Y, V, levels=30, cmap='viridis')
ax.clabel(contour, inline=True, fontsize=8)

ax.plot(x_eq, y_eq, 'ro', markersize=10)
ax.set_xlabel('Prey x', fontsize=12)
ax.set_ylabel('Predator y', fontsize=12)
ax.set_title('Lotka-Volterra: Conserved Quantity (Level Sets)', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lotka_volterra_energy.png', dpi=150, bbox_inches='tight')
plt.show()

Application: Control System Design

Feedback Stabilization

Consider an unstable system, how can we stabilize it through feedback?

State feedback:Closed-loop system: Goal: Choosesuch that all eigenvalues ofhave negative real parts.

Pole placement theorem: Ifis controllable, then the eigenvalues ofcan be arbitrarily assigned.

Example: Inverted Pendulum Control

Linearized inverted pendulum equations:whereis the force applied to the cart.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
import numpy as np
from scipy.linalg import solve_continuous_are
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Inverted pendulum parameters
M = 1.0  # Cart mass
m = 0.1  # Pendulum mass
l = 0.5  # Pendulum length
g = 9.81  # Gravitational acceleration

# State space matrices (linearized around θ=0)
A = np.array([
    [0, 1, 0, 0],
    [0, 0, -m*g/M, 0],
    [0, 0, 0, 1],
    [0, 0, (M+m)*g/(M*l), 0]
])

B = np.array([[0], [1/M], [0], [-1/(M*l)]])

print("Open-loop eigenvalues:", np.linalg.eigvals(A))
print("System is unstable (has positive eigenvalue)")

# LQR design
Q = np.diag([10, 1, 100, 1])  # State weights
R = np.array([[0.1]])  # Control weight

# Solve Riccati equation
P = solve_continuous_are(A, B, Q, R)
K = np.linalg.inv(R) @ B.T @ P

print("\nFeedback gain K:", K.flatten())
print("Closed-loop eigenvalues:", np.linalg.eigvals(A - B @ K))
print("System is stable (all eigenvalues have Re<0)")

# Simulation
def inverted_pendulum(X, t, K, A, B):
    u = -K @ X
    return (A @ X + B @ u).flatten()

def inverted_pendulum_open(X, t, A):
    return (A @ X).flatten()

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

t = np.linspace(0, 5, 500)
X0 = [0, 0, 0.2, 0]  # Initial angle deviation of 0.2 rad

# Open-loop response
sol_open = odeint(inverted_pendulum_open, X0, t, args=(A,))

# Closed-loop response
sol_closed = odeint(inverted_pendulum, X0, t, args=(K, A, B))

# Angle
ax1 = axes[0, 0]
ax1.plot(t, sol_open[:, 2] * 180/np.pi, 'r--', linewidth=2, label='Open-loop')
ax1.plot(t, sol_closed[:, 2] * 180/np.pi, 'b-', linewidth=2, label='Closed-loop')
ax1.set_xlabel('Time (s)', fontsize=12)
ax1.set_ylabel('Angle θ (degrees)', fontsize=12)
ax1.set_title('Pendulum Angle', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# Cart position
ax2 = axes[0, 1]
ax2.plot(t, sol_open[:, 0], 'r--', linewidth=2, label='Open-loop')
ax2.plot(t, sol_closed[:, 0], 'b-', linewidth=2, label='Closed-loop')
ax2.set_xlabel('Time (s)', fontsize=12)
ax2.set_ylabel('Position x (m)', fontsize=12)
ax2.set_title('Cart Position', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# Angular velocity
ax3 = axes[1, 0]
ax3.plot(t, sol_open[:, 3] * 180/np.pi, 'r--', linewidth=2, label='Open-loop')
ax3.plot(t, sol_closed[:, 3] * 180/np.pi, 'b-', linewidth=2, label='Closed-loop')
ax3.set_xlabel('Time (s)', fontsize=12)
ax3.set_ylabel('Angular velocity (deg/s)', fontsize=12)
ax3.set_title('Angular Velocity', fontsize=14, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# Control input
ax4 = axes[1, 1]
u_closed = np.array([-K @ sol_closed[i] for i in range(len(t))]).flatten()
ax4.plot(t, u_closed, 'b-', linewidth=2)
ax4.set_xlabel('Time (s)', fontsize=12)
ax4.set_ylabel('Control Force u (N)', fontsize=12)
ax4.set_title('Control Input', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.suptitle('Inverted Pendulum Stabilization via LQR', fontsize=16, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('inverted_pendulum.png', dpi=150, bbox_inches='tight')
plt.show()

Summary

Key Points of This Chapter

  1. Lyapunov stability:
    • Stable: doesn't go far
    • Asymptotically stable: eventually returns to equilibrium
    • Unstable: will go far
  2. Linearization method:
    • Jacobian matrix
    • Eigenvalues determine stability
    • Hartman-Grobman theorem
  3. Lyapunov's direct method:
    • Lyapunov functions (generalized energy)
    • LaSalle's invariance principle
  4. Bifurcation theory:
    • Saddle-node, transcritical, pitchfork bifurcations
    • Hopf bifurcation and limit cycles
  5. Applications:
    • Ecosystems, circuits, control

Preview of Next Chapter

In "Numerical Methods," we will: - Learn various ODE numerical solving algorithms - Understand step size selection and error control - Handle stiff equations - Implement adaptive methods

Exercises

Basic Problems

  1. Determine the stability of the origin for the following systems: -, -,(near origin)

  2. For, use Lyapunov functionto analyze the stability of the origin.

  3. Find all bifurcation points of the equationand draw the bifurcation diagram.

  4. For the damped pendulum, prove that total energyis a Lyapunov function.

  5. Calculate the trace and determinant of the following matrix and determine the equilibrium type:

Advanced Problems

  1. Bistable system: Analyze all equilibrium points and their stability for. Draw the-graph and trajectories.

  2. Hysteresis phenomenon: For,

    • Draw the-graph for different values of - Find bifurcation points
    • Explain why hysteresis occurs

  3. Prove that ifis a Lyapunov function with, then all bounded trajectories approach the set of equilibria.

  4. Competitive exclusion: Lotka-Volterra competition model: Analyze stability whenand when.

Programming Problems

  1. Write a program to automatically determine the equilibrium type of 2D linear systems (given matrix).

  2. Implement a Lyapunov equation solver and use it to verify stability of linear systems.

  3. Simulate the Van der Pol oscillator, creating an animation showing the formation of the limit cycle.

  4. Design a PID controller for the inverted pendulum and compare performance with the LQR controller.

Thought Questions

  1. Why does the Hartman-Grobman theorem require eigenvalues to have nonzero real parts? Give a counterexample showing how linearization can fail in the center case.

  2. Is the energy of a physical system always a good Lyapunov function candidate? When isn't it?

  3. What is "critical slowing down" near bifurcation points? What is its practical significance?

References

  1. Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos. CRC Press.
  2. Khalil, H. K. (2002). Nonlinear Systems. Prentice Hall.
  3. Perko, L. (2001). Differential Equations and Dynamical Systems. Springer.
  4. Guckenheimer, J., & Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer.
  5. Slotine, J.-J. E., & Li, W. (1991). Applied Nonlinear Control. Prentice Hall.

Next Chapter: Nonlinear Systems and Phase Portraits →

This is Chapter 7 of the "World of Ordinary Differential Equations" series.

  • Post title：Ordinary Differential Equations (7): Stability Theory
  • Post author：Chen Kai
  • Create time：2019-05-07 14:00:00
  • Post link：https://www.chenk.top/ode-chapter-07-systems-and-phase-plane/
  • Copyright Notice：All articles in this blog are licensed under BY-NC-SA unless stating additionally.
 Comments