Ordinary Differential Equations (10): Bifurcation Theory
Chen Kai BOSS
  2019-05-24 09:15:00 2019-05-24 09:15      5k Words  30 Mins  

Why do lakes sometimes shift abruptly from clear to turbid? Why do stock markets crash without warning? Why do climate systems seem to undergo dramatic changes at certain "tipping points"? Behind these seemingly unrelated phenomena lies a common mathematical principle —bifurcation. When system parameters change slowly, the qualitative behavior of a system can suddenly transform: stable equilibria disappear, periodic solutions emerge, and chaos erupts. Bifurcation theory is the mathematical tool for studying such "qualitative changes," enabling us to predict and understand various abrupt phenomena in nature.

What Is Bifurcation?

Starting with Everyday Examples

Imagine boiling water. Below 100° C, water remains liquid; at 100° C, it starts boiling and becomes vapor. This is a type of "bifurcation"— the qualitative state of the system fundamentally changes at a critical parameter value.

Now imagine a bridge. Within design limits, the bridge is stable; but when the weight exceeds a critical value, the bridge suddenly collapses. This transition from stability to instability is a manifestation of bifurcation.

Mathematical Definition

Bifurcation refers to the phenomenon where the topological structure of a system undergoes a qualitative change when a parameter crosses a critical value.

Consider a parameterized dynamical system: whereis the parameter. Asvaries: - The number of equilibrium points may change - The stability of equilibrium points may change - Periodic solutions may appear or disappear - Chaos may emerge

The parameter valueat which such changes occur is called the bifurcation value or critical point.

Bifurcation Diagrams: A Panoramic View

The bifurcation diagram is the most intuitive tool for understanding bifurcations. It plots the parameteron the horizontal axis and the system's steady states (equilibrium points, amplitudes of periodic solutions, etc.) on the vertical axis, showing how system behavior changes with the parameter.

Bifurcations in One-Dimensional Systems

One-dimensional systemsare the ideal starting point for studying bifurcations. Though simple, they already exhibit the core features of bifurcation.

Saddle-Node Bifurcation

The most basic type of bifurcation— two equilibrium points "collide" and annihilate.

Normal form:Analysis: - When: No real equilibrium points (has no solution) - When: One equilibrium point at (semi-stable) - When: Two equilibrium points atStability analysis: -: (stable) -: (unstable)

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
import numpy as np
import matplotlib.pyplot as plt

def saddle_node_bifurcation():
    """Visualization of saddle-node bifurcation"""
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Left: Phase lines for different mu values
    ax1 = axes[0]
    x = np.linspace(-2, 2, 500)
    
    mu_values = [-0.5, 0, 0.5, 1.0]
    colors = ['blue', 'green', 'orange', 'red']
    
    for mu, color in zip(mu_values, colors):
        f = mu - x**2
        ax1.plot(x, f, color=color, linewidth=2, label=f'μ = {mu}')
        
        # Mark equilibrium points
        if mu >= 0:
            x_eq = np.sqrt(mu)
            ax1.plot(x_eq, 0, 'o', color=color, markersize=10)
            ax1.plot(-x_eq, 0, 's', color=color, markersize=10)
    
    ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
    ax1.set_xlabel('x', fontsize=12)
    ax1.set_ylabel('dx/dt = μ - x ²', fontsize=12)
    ax1.set_title('Phase Lines for Different μ', fontsize=13, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(-2, 2)
    ax1.set_ylim(-2, 2)
    
    # Right: Bifurcation diagram
    ax2 = axes[1]
    mu_range = np.linspace(-1, 2, 300)
    
    # Stable branch (x = sqrt(mu))
    mu_pos = mu_range[mu_range >= 0]
    x_stable = np.sqrt(mu_pos)
    ax2.plot(mu_pos, x_stable, 'b-', linewidth=2.5, label='Stable')
    ax2.plot(mu_pos, -x_stable, 'r--', linewidth=2.5, label='Unstable')
    
    # Bifurcation point
    ax2.plot(0, 0, 'ko', markersize=12, zorder=5)
    ax2.annotate('Bifurcation point', xy=(0, 0), xytext=(0.5, 0.8),
                fontsize=11, arrowprops=dict(arrowstyle='->', color='black'))
    
    ax2.axvline(x=0, color='gray', linestyle=':', alpha=0.7)
    ax2.set_xlabel('μ', fontsize=12)
    ax2.set_ylabel('x* (equilibria)', fontsize=12)
    ax2.set_title('Saddle-Node Bifurcation Diagram', fontsize=13, fontweight='bold')
    ax2.legend(fontsize=11)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(-1, 2)
    ax2.set_ylim(-2, 2)
    
    plt.tight_layout()
    plt.savefig('saddle_node_bifurcation.png', dpi=150, bbox_inches='tight')
    plt.close()

saddle_node_bifurcation()

Physical meaning: Saddle-node bifurcation describes "threshold phenomena"— below the critical value there is no equilibrium state, but one suddenly appears when the threshold is crossed.

Real-world examples: - Lasers: No laser output below threshold current, sudden emission above threshold - Neurons: No firing below threshold voltage, action potential generated above threshold

Transcritical Bifurcation

Two equilibrium points exchange stability, but neither disappears.

Normal form:Analysis: - Equilibrium points are alwaysand - Stability: - When:stable,unstable - When:unstable,stable

At, the two equilibrium points "exchange" stability.

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
import numpy as np
import matplotlib.pyplot as plt

def transcritical_bifurcation():
    """Visualization of transcritical bifurcation"""
    
    fig, ax = plt.subplots(figsize=(10, 7))
    
    mu_range = np.linspace(-2, 2, 400)
    
    # x* = 0 branch
    ax.plot(mu_range[mu_range <= 0], np.zeros(sum(mu_range <= 0)), 
            'b-', linewidth=2.5, label='x*=0 (stable when μ<0)')
    ax.plot(mu_range[mu_range >= 0], np.zeros(sum(mu_range >= 0)), 
            'b--', linewidth=2.5, label='x*=0 (unstable when μ>0)')
    
    # x* = μ branch
    ax.plot(mu_range[mu_range <= 0], mu_range[mu_range <= 0], 
            'r--', linewidth=2.5, label='x*=μ (unstable when μ<0)')
    ax.plot(mu_range[mu_range >= 0], mu_range[mu_range >= 0], 
            'r-', linewidth=2.5, label='x*=μ (stable when μ>0)')
    
    # Bifurcation point
    ax.plot(0, 0, 'ko', markersize=15, zorder=5)
    ax.annotate('Transcritical\nbifurcation', xy=(0, 0), xytext=(0.8, -1.2),
                fontsize=11, arrowprops=dict(arrowstyle='->', color='black'))
    
    ax.set_xlabel('μ', fontsize=12)
    ax.set_ylabel('x*', fontsize=12)
    ax.set_title('Transcritical Bifurcation: Exchange of Stability', 
                fontsize=13, fontweight='bold')
    ax.legend(fontsize=10, loc='upper left')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.axhline(y=0, color='gray', linestyle='-', linewidth=0.5)
    ax.axvline(x=0, color='gray', linestyle='-', linewidth=0.5)
    
    plt.tight_layout()
    plt.savefig('transcritical_bifurcation.png', dpi=150, bbox_inches='tight')
    plt.close()

transcritical_bifurcation()

Physical meaning: Transcritical bifurcation is common in systems with a "zero state"— the zero state is always an equilibrium, but its stability changes with the parameter.

Real-world examples: - Epidemic models:is the critical point -: Disease-free equilibrium is stable (disease dies out) -: Disease-free equilibrium is unstable (epidemic outbreak) - Population growth: Species goes extinct below a certain environmental quality, survives above it

Pitchfork Bifurcation

One equilibrium point splits into three— the typical bifurcation for symmetric systems.

Supercritical pitchfork bifurcation:Analysis: - When: Only(stable) - When: Three equilibrium points:(unstable),(stable)

Subcritical pitchfork bifurcation: - When: Only(unstable) - When: Three equilibrium points:(stable),(unstable)

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
import numpy as np
import matplotlib.pyplot as plt

def pitchfork_bifurcations():
    """Supercritical and subcritical pitchfork bifurcations"""
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    mu_range = np.linspace(-2, 2, 400)
    mu_pos = mu_range[mu_range >= 0]
    mu_neg = mu_range[mu_range <= 0]
    
    # Left: Supercritical
    ax1 = axes[0]
    
    # x = 0 branch
    ax1.plot(mu_range[mu_range <= 0], np.zeros(sum(mu_range <= 0)), 
             'b-', linewidth=2.5)
    ax1.plot(mu_range[mu_range >= 0], np.zeros(sum(mu_range >= 0)), 
             'r--', linewidth=2.5)
    
    # x = ± sqrt(mu) branch
    ax1.plot(mu_pos, np.sqrt(mu_pos), 'b-', linewidth=2.5)
    ax1.plot(mu_pos, -np.sqrt(mu_pos), 'b-', linewidth=2.5)
    
    ax1.plot(0, 0, 'ko', markersize=12)
    ax1.set_xlabel('μ', fontsize=12)
    ax1.set_ylabel('x*', fontsize=12)
    ax1.set_title('Supercritical Pitchfork Bifurcation\n(dx/dt = μ x - x ³)', 
                  fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(-2, 2)
    ax1.set_ylim(-2, 2)
    ax1.legend(['Stable', 'Unstable'], fontsize=10, loc='upper left')
    
    # Right: Subcritical
    ax2 = axes[1]
    
    # x = 0 branch
    ax2.plot(mu_range[mu_range >= 0], np.zeros(sum(mu_range >= 0)), 
             'r--', linewidth=2.5)
    ax2.plot(mu_range[mu_range <= 0], np.zeros(sum(mu_range <= 0)), 
             'b-', linewidth=2.5)
    
    # x = ± sqrt(-mu) branch
    ax2.plot(mu_neg, np.sqrt(-mu_neg), 'r--', linewidth=2.5)
    ax2.plot(mu_neg, -np.sqrt(-mu_neg), 'r--', linewidth=2.5)
    
    ax2.plot(0, 0, 'ko', markersize=12)
    ax2.set_xlabel('μ', fontsize=12)
    ax2.set_ylabel('x*', fontsize=12)
    ax2.set_title('Subcritical Pitchfork Bifurcation\n(dx/dt = μ x + x ³)', 
                  fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(-2, 2)
    ax2.set_ylim(-2, 2)
    ax2.legend(['Unstable', 'Stable'], fontsize=10, loc='upper left')
    
    plt.tight_layout()
    plt.savefig('pitchfork_bifurcations.png', dpi=150, bbox_inches='tight')
    plt.close()

pitchfork_bifurcations()

Physical meaning: - Supercritical: The system smoothly transitions to a new state after bifurcation ("soft" bifurcation) - Subcritical: The system may suddenly jump to a distant state before bifurcation ("hard" bifurcation — dangerous!)

Real-world examples: - Euler buckling: A vertical slender rod under compression bends after exceeding the critical load (left or right, breaking symmetry) - Lasers: Above threshold, the optical field phase breaks symmetry

Summary of Bifurcation Types

Bifurcation Type Normal Form Characteristics Danger Level
Saddle-node Equilibria created/destroyed Medium
Transcritical Stability exchange Low
Supercritical pitchfork Symmetric splitting, soft transition Low
Subcritical pitchfork Symmetric splitting, hard jump High

Bifurcations in Two-Dimensional Systems

Two-dimensional systems can exhibit richer bifurcations, particularly those related to periodic solutions.

Hopf Bifurcation: Birth of Periodic Solutions

Hopf bifurcation is the most important two-dimensional bifurcation — a stable focus becomes an unstable focus while simultaneously giving birth to a limit cycle.

Normal form (polar coordinates):Analyzing theequation (same as supercritical pitchfork): - When:is stable (focus) - When:is unstable,is stable (limit cycle)

Limit cycle radius: Cartesian form: 

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
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def hopf_system(state, t, mu, omega):
    """Hopf bifurcation normal form"""
    x, y = state
    r_sq = x**2 + y**2
    dxdt = mu * x - omega * y - x * r_sq
    dydt = omega * x + mu * y - y * r_sq
    return [dxdt, dydt]

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
omega = 1.0
t = np.linspace(0, 50, 5000)

mu_values = [-0.5, 0, 0.5]
titles = ['μ = -0.5 (Stable Focus)', 'μ = 0 (Bifurcation Point)', 
          'μ = 0.5 (Limit Cycle)']

for ax, mu, title in zip(axes, mu_values, titles):
    # Multiple initial conditions
    initial_conditions = [(0.1, 0), (0.5, 0.5), (1.5, 0), (0, 1.5), (-1, -1)]
    
    for x0, y0 in initial_conditions:
        sol = odeint(hopf_system, [x0, y0], t, args=(mu, omega))
        ax.plot(sol[:, 0], sol[:, 1], linewidth=0.7, alpha=0.7)
        ax.plot(x0, y0, 'go', markersize=5)
    
    # Draw limit cycle if present
    if mu > 0:
        theta = np.linspace(0, 2*np.pi, 100)
        r = np.sqrt(mu)
        ax.plot(r * np.cos(theta), r * np.sin(theta), 'r-', 
                linewidth=2.5, label=f'Limit cycle (r = {r:.2f})')
        ax.legend(fontsize=9)
    
    ax.plot(0, 0, 'ko', markersize=8)
    ax.set_xlabel('x', fontsize=11)
    ax.set_ylabel('y', fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_aspect('equal')

plt.tight_layout()
plt.savefig('hopf_bifurcation.png', dpi=150, bbox_inches='tight')
plt.close()

Supercritical vs. Subcritical Hopf Bifurcation

Supercritical Hopf (as above): - Limit cycle radius grows continuously from zero after bifurcation - System transitions smoothly, relatively safe

Subcritical Hopf: - Unstable limit cycle exists before bifurcation - System may jump to a distant state after bifurcation - Hysteresis phenomenon exists

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
import numpy as np
import matplotlib.pyplot as plt

def hopf_bifurcation_diagram():
    """Comparison of Hopf bifurcation diagrams"""
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    mu_range = np.linspace(-1, 1, 200)
    mu_pos = mu_range[mu_range >= 0]
    mu_neg = mu_range[mu_range <= 0]
    
    # Left: Supercritical
    ax1 = axes[0]
    
    # Origin stability
    ax1.plot(mu_neg, np.zeros(len(mu_neg)), 'b-', linewidth=2.5, label='Stable equilibrium')
    ax1.plot(mu_pos, np.zeros(len(mu_pos)), 'r--', linewidth=2.5, label='Unstable equilibrium')
    
    # Limit cycle
    ax1.plot(mu_pos, np.sqrt(mu_pos), 'b-', linewidth=2.5, label='Stable limit cycle')
    ax1.plot(mu_pos, -np.sqrt(mu_pos), 'b-', linewidth=2.5)
    
    ax1.fill_between(mu_pos, -np.sqrt(mu_pos), np.sqrt(mu_pos), alpha=0.2, color='blue')
    
    ax1.plot(0, 0, 'ko', markersize=12)
    ax1.set_xlabel('μ', fontsize=12)
    ax1.set_ylabel('Amplitude', fontsize=12)
    ax1.set_title('Supercritical Hopf Bifurcation', fontsize=13, fontweight='bold')
    ax1.legend(fontsize=9, loc='upper left')
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(-1, 1)
    ax1.set_ylim(-1.5, 1.5)
    
    # Right: Subcritical (with hysteresis)
    ax2 = axes[1]
    
    # Origin stability
    ax2.plot(mu_neg, np.zeros(len(mu_neg)), 'b-', linewidth=2.5)
    ax2.plot(mu_pos, np.zeros(len(mu_pos)), 'r--', linewidth=2.5)
    
    # Unstable limit cycle (exists in mu<0 region)
    mu_unstable = mu_range[(mu_range >= -0.5) & (mu_range <= 0)]
    r_unstable = np.sqrt(-mu_unstable + 0.5)
    ax2.plot(mu_unstable, r_unstable, 'r--', linewidth=2.5, label='Unstable limit cycle')
    ax2.plot(mu_unstable, -r_unstable, 'r--', linewidth=2.5)
    
    # Stable limit cycle (large amplitude)
    mu_stable = mu_range[mu_range >= -0.5]
    r_stable = np.sqrt(mu_stable + 0.5)
    ax2.plot(mu_stable, r_stable, 'b-', linewidth=2.5, label='Stable limit cycle')
    ax2.plot(mu_stable, -r_stable, 'b-', linewidth=2.5)
    
    # Hysteresis region
    ax2.fill_between(mu_neg[mu_neg >= -0.5], 
                     -np.sqrt(-mu_neg[mu_neg >= -0.5] + 0.5),
                     np.sqrt(-mu_neg[mu_neg >= -0.5] + 0.5),
                     alpha=0.2, color='red', label='Bistable region')
    
    ax2.plot(0, 0, 'ko', markersize=12)
    ax2.plot(-0.5, 0, 'ko', markersize=12)
    ax2.set_xlabel('μ', fontsize=12)
    ax2.set_ylabel('Amplitude', fontsize=12)
    ax2.set_title('Subcritical Hopf Bifurcation\n(with hysteresis)', 
                  fontsize=13, fontweight='bold')
    ax2.legend(fontsize=9, loc='upper left')
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(-1, 1)
    ax2.set_ylim(-1.5, 1.5)
    
    plt.tight_layout()
    plt.savefig('hopf_comparison.png', dpi=150, bbox_inches='tight')
    plt.close()

hopf_bifurcation_diagram()

The Hopf Bifurcation Theorem

Theorem (Hopf, 1942): Letbe an equilibrium point of the systemat, with the Jacobian matrix having a pair of complex conjugate eigenvalues. If:

1.(real part crosses zero) 2.(imaginary part is nonzero) 3.(transversality condition)

then a Hopf bifurcation occurs at, with a branch of periodic solutions emerging with period approximately.

Global Bifurcations

The bifurcations discussed above are local bifurcations— occurring near equilibrium points or limit cycles. Global bifurcations involve structural changes over larger regions of phase space.

Homoclinic Bifurcation

A homoclinic orbit is a trajectory that starts from a saddle point and returns to the same saddle point.

When parameter changes cause a homoclinic orbit to appear or disappear, a homoclinic bifurcation occurs.

Special features of this bifurcation: - The orbit's period approaches infinity (since it spends infinite time near the saddle) - May lead to chaos (Shilnikov theorem)

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
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def homoclinic_example(state, t, mu):
    """Example system demonstrating homoclinic orbits"""
    x, y = state
    dxdt = y
    dydt = x - x**3 + mu * y  # Adjustable damping
    return [dxdt, dydt]

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
t = np.linspace(0, 50, 5000)

# Different mu values
mu_values = [-0.3, -0.05, 0.3]
titles = ['μ = -0.3: Stable limit cycles', 
          'μ ≈ 0: Near homoclinic bifurcation',
          'μ = 0.3: Trajectories escape to infinity']

for ax, mu, title in zip(axes, mu_values, titles):
    # Equilibrium points
    ax.plot(0, 0, 'rs', markersize=10, label='Saddle at origin')
    ax.plot(1, 0, 'bo', markersize=8, label='Stable/unstable node')
    ax.plot(-1, 0, 'bo', markersize=8)
    
    # Multiple trajectories
    initial_conditions = [
        (0.5, 0.2), (-0.5, 0.2), (0.5, -0.2), (-0.5, -0.2),
        (0.1, 0.5), (-0.1, 0.5), (0.1, -0.5), (-0.1, -0.5),
        (1.5, 0), (-1.5, 0)
    ]
    
    for x0, y0 in initial_conditions:
        try:
            sol = odeint(homoclinic_example, [x0, y0], t, args=(mu,))
            # Only plot bounded portion
            mask = (np.abs(sol[:, 0]) < 3) & (np.abs(sol[:, 1]) < 3)
            if np.any(mask):
                ax.plot(sol[mask, 0], sol[mask, 1], 'g-', linewidth=0.7, alpha=0.7)
        except:
            pass
        ax.plot(x0, y0, 'ko', markersize=4)
    
    ax.set_xlabel('x', fontsize=11)
    ax.set_ylabel('y', fontsize=11)
    ax.set_title(title, fontsize=11, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-2, 2)
    ax.set_ylim(-1.5, 1.5)
    ax.legend(fontsize=8, loc='upper right')

plt.tight_layout()
plt.savefig('homoclinic_bifurcation.png', dpi=150, bbox_inches='tight')
plt.close()

Heteroclinic Bifurcation

A heteroclinic orbit connects two different saddle points.

When multiple saddle points exist, the connections between them may change with parameters, leading to heteroclinic bifurcation.

SNIC Bifurcation

SNIC (Saddle-Node on Invariant Circle)— a saddle-node bifurcation occurring on a limit cycle.

Result: The limit cycle "opens" into an invariant curve through a saddle-node.

Characteristic: Period diverges as.

Codimension and Universality

The Concept of Codimension

Codimension is the minimum number of independent parameters that must be tuned to cause a bifurcation.

  • Codimension-1 bifurcations: Occur by tuning one parameter (saddle-node, transcritical, Hopf, etc.)
  • Codimension-2 bifurcations: Require simultaneously tuning two parameters (Bogdanov-Takens, cusp, etc.)

Higher codimension means the bifurcation is more "unstable"— harder to encounter in parameter space.

Universal Unfolding and Normal Forms

The normal form is the "simplest representative" of a bifurcation.

Through coordinate transformations, any system satisfying certain bifurcation conditions can be reduced to its normal form (near the bifurcation point).

This is the power of bifurcation theory: although specific systems vary greatly, the types of bifurcations are finite!

Applications in Natural Sciences

Ecology: Species Extinction and Recovery

Consider a population model with the Allee effect:whereis the Allee threshold — below this value, the population heads toward extinction.

When the environment deteriorates (decreases), a saddle-node bifurcation may occur: the positive equilibrium disappears, and the population suddenly collapses!

This explains why fishery resources can "suddenly" deplete — not a linear decline, but collapse after crossing a critical point.

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
import numpy as np
import matplotlib.pyplot as plt

def allee_effect_bifurcation():
    """Bifurcation caused by the Allee effect"""
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Parameters
    r = 1.0
    A = 10.0  # Allee threshold
    
    K_values = [50, 30, 20, 15]  # Different carrying capacities
    colors = ['green', 'blue', 'orange', 'red']
    
    N = np.linspace(0, 60, 500)
    
    # Left: Phase diagram
    ax1 = axes[0]
    
    for K, color in zip(K_values, colors):
        dNdt = r * N * (1 - N/K) * (N/A - 1)
        ax1.plot(N, dNdt, color=color, linewidth=2, label=f'K = {K}')
        
        # Mark equilibrium points
        from scipy.optimize import brentq
        try:
            # Find positive equilibria
            if K > A:
                N_stable = brentq(lambda n: r*n*(1-n/K)*(n/A-1), A+1, K-1)
                N_unstable = brentq(lambda n: r*n*(1-n/K)*(n/A-1), 1, A-1)
                ax1.plot(N_stable, 0, 'o', color=color, markersize=8)
                ax1.plot(N_unstable, 0, 's', color=color, markersize=8)
        except:
            pass
    
    ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
    ax1.set_xlabel('Population N', fontsize=12)
    ax1.set_ylabel('dN/dt', fontsize=12)
    ax1.set_title('Allee Effect: Phase Lines for Different K', 
                  fontsize=13, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(0, 60)
    ax1.set_ylim(-5, 10)
    
    # Right: Bifurcation diagram
    ax2 = axes[1]
    
    K_range = np.linspace(10, 60, 200)
    
    # Stable high-density branch
    N_stable = []
    N_unstable = []
    K_valid_s = []
    K_valid_u = []
    
    for K in K_range:
        # Stable: root between A and K
        if K > A:
            try:
                ns = brentq(lambda n: r*n*(1-n/K)*(n/A-1), A+0.1, K-0.1)
                N_stable.append(ns)
                K_valid_s.append(K)
            except:
                pass
            try:
                nu = brentq(lambda n: r*n*(1-n/K)*(n/A-1), 0.1, A-0.1)
                N_unstable.append(nu)
                K_valid_u.append(K)
            except:
                pass
    
    ax2.plot(K_valid_s, N_stable, 'b-', linewidth=2.5, label='Stable (high density)')
    ax2.plot(K_valid_u, N_unstable, 'r--', linewidth=2.5, label='Unstable (Allee threshold)')
    ax2.plot(K_range, np.zeros(len(K_range)), 'b-', linewidth=2.5, label='Stable (extinction)')
    
    # Bifurcation point
    ax2.axvline(x=A, color='gray', linestyle=':', linewidth=2)
    ax2.text(A+1, 50, 'Saddle-node\nbifurcation', fontsize=10)
    
    ax2.set_xlabel('Carrying capacity K', fontsize=12)
    ax2.set_ylabel('Equilibrium population N*', fontsize=12)
    ax2.set_title('Bifurcation Diagram: Population Collapse', 
                  fontsize=13, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(10, 60)
    ax2.set_ylim(-5, 60)
    
    plt.tight_layout()
    plt.savefig('allee_bifurcation.png', dpi=150, bbox_inches='tight')
    plt.close()

allee_effect_bifurcation()

Climate Science: Tipping Points and Irreversible Change

Climate systems may have multiple stable states (e.g., "greenhouse Earth" and "icehouse Earth").

When greenhouse gas concentrations cross a critical point, an irreversible state transition may occur — these are the "tipping points" that climate scientists warn about.

Examples: - Arctic sea ice: Albedo-temperature positive feedback may cause sudden disappearance - Amazon rainforest: Drought-fire-forest degradation may lead to savannification - Atlantic circulation: Freshwater injection may cause circulation collapse

These are all manifestations of saddle-node bifurcations or subcritical Hopf bifurcations.

Neuroscience: Neuronal Firing Patterns

Neurons can switch between "resting" and "firing" states.

The Hodgkin-Huxley model (simplified as FitzHugh-Nagumo) exhibits multiple bifurcations: - Saddle-node bifurcation: Class I neurons (frequency continuously adjustable from zero) - Hopf bifurcation: Class II neurons (frequency has a lower bound)

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
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def fitzhugh_nagumo(state, t, I, a, b, tau):
    """FitzHugh-Nagumo model"""
    v, w = state
    dvdt = v - v**3/3 - w + I
    dwdt = (v + a - b*w) / tau
    return [dvdt, dwdt]

# Behavior under different external currents
I_values = [0, 0.3, 0.5, 0.8]
a, b, tau = 0.7, 0.8, 12.5

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
t = np.linspace(0, 200, 4000)

for ax, I in zip(axes.flatten(), I_values):
    # Phase diagram
    v_range = np.linspace(-2.5, 2.5, 20)
    w_range = np.linspace(-1.5, 2, 20)
    V, W = np.meshgrid(v_range, w_range)
    dV = V - V**3/3 - W + I
    dW = (V + a - b*W) / tau
    
    # Vector field
    M = np.sqrt(dV**2 + dW**2)
    M[M == 0] = 1
    ax.quiver(V, W, dV/M, dW/M, M, cmap='Greys', alpha=0.3)
    
    # Nullclines
    v_null = np.linspace(-2.5, 2.5, 100)
    w_vnull = v_null - v_null**3/3 + I
    w_wnull = (v_null + a) / b
    
    ax.plot(v_null, w_vnull, 'b-', linewidth=2, label='v-nullcline')
    ax.plot(v_null, w_wnull, 'r-', linewidth=2, label='w-nullcline')
    
    # Trajectory
    sol = odeint(fitzhugh_nagumo, [-1, 0], t, args=(I, a, b, tau))
    ax.plot(sol[:, 0], sol[:, 1], 'g-', linewidth=1.5, alpha=0.7)
    ax.plot(sol[0, 0], sol[0, 1], 'go', markersize=8)
    
    ax.set_xlabel('v (membrane potential)', fontsize=11)
    ax.set_ylabel('w (recovery variable)', fontsize=11)
    ax.set_title(f'I = {I}', fontsize=12, fontweight='bold')
    ax.legend(fontsize=9, loc='upper right')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-2.5, 2.5)
    ax.set_ylim(-1.5, 2)

plt.suptitle('FitzHugh-Nagumo Model: Hopf Bifurcation with Increasing Current', 
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('fhn_bifurcation.png', dpi=150, bbox_inches='tight')
plt.close()

Engineering: Structural Stability

Euler buckling is a classic pitchfork bifurcation:

A slender column under axial load: -: Upright state is stable -: Bifurcation point (critical load) -: Upright state unstable, bent state is stable

Critical load(: elastic modulus,: moment of inertia,: length)

This is a supercritical pitchfork bifurcation— the column smoothly bends to one side.

Economics: Market Crashes

Financial markets may have multiple "equilibria": - High-confidence, high-investment prosperity state - Low-confidence, low-investment recession state

When certain indicators cross a critical value, a saddle-node bifurcation may occur — markets suddenly jump from prosperity to recession (crash).

This explains why market crashes tend to be sudden rather than gradual.

Period Doubling and the Route to Chaos

Period-Doubling Cascade

As parameters change, periodic solutions may undergo period-doubling bifurcations: - Period-1 → Period-2 → Period-4 → Period-8 → ... → Chaos

Feigenbaum constant:This constant is universal— regardless of the specific system, as long as it undergoes period doubling, the ratio converges to this constant!

The Logistic Map: A Simple Model of Chaos

Logistic map:This is the simplest model for studying chaos. Asincreases:

range Behavior
All orbits tend to 0
Tends to nonzero fixed point
Period-2
Period-4
... Period-doubling cascade
Chaos (with periodic windows) 
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
import numpy as np
import matplotlib.pyplot as plt

def logistic_bifurcation_diagram():
    """Bifurcation diagram of the logistic map"""
    
    r_values = np.linspace(2.5, 4.0, 2000)
    iterations = 1000
    last_points = 100
    
    fig, ax = plt.subplots(figsize=(14, 8))
    
    for r in r_values:
        x = 0.5  # Initial value
        
        # Iterate to steady state
        for _ in range(iterations - last_points):
            x = r * x * (1 - x)
        
        # Collect last few points
        x_values = []
        for _ in range(last_points):
            x = r * x * (1 - x)
            x_values.append(x)
        
        ax.scatter([r] * len(x_values), x_values, 
                   c='blue', s=0.1, alpha=0.5)
    
    ax.set_xlabel('r', fontsize=12)
    ax.set_ylabel('x*', fontsize=12)
    ax.set_title('Bifurcation Diagram of Logistic Map:$x_{n+1} = rx_n(1-x_n)$', 
                 fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # Mark key points
    ax.axvline(x=3, color='red', linestyle='--', alpha=0.5)
    ax.text(3.02, 0.1, 'Period-2\nbifurcation', fontsize=10, color='red')
    
    ax.axvline(x=3.5699, color='green', linestyle='--', alpha=0.5)
    ax.text(3.59, 0.1, 'Onset of\nchaos', fontsize=10, color='green')
    
    plt.tight_layout()
    plt.savefig('logistic_bifurcation.png', dpi=150, bbox_inches='tight')
    plt.close()

logistic_bifurcation_diagram()

Periodic Windows and Crises

Even in the chaotic regime, there are periodic windows— parameter intervals where the system returns to periodic motion.

The most famous is the period-3 window near.

Li-Yorke Theorem: "Period 3 implies chaos"— if a one-dimensional map has a period-3 orbit, then it has orbits of all periods, and there exist uncountably many chaotic orbits.

Numerical Methods: Continuation and Detection

Continuation Methods

How do we numerically "track" equilibrium points or periodic solutions as parameters vary?

Pseudo-arclength continuation: 1. Define an arclength parameter along the solution curve 2. Advance in the arclength direction while adjusting both parameter and state 3. Can track "fold-back" branches (at saddle-nodes)

Common software: - AUTO (classic tool) - MATCONT (MATLAB toolbox) - PyDSTool (Python) - XPPAut

Bifurcation Detection

Detecting saddle-node bifurcations: Monitor eigenvalues of the Jacobian matrix; alert when an eigenvalue crosses zero.

Detecting Hopf bifurcations: Monitor real parts of eigenvalues; Hopf bifurcation occurs when a pair of complex conjugate eigenvalues have their real parts cross zero.

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
import numpy as np
from scipy.linalg import eigvals

def detect_bifurcations(jacobian_func, param_range, x_eq_func):
    """
    Detect bifurcations along a parameter path
    
    jacobian_func: Function to compute Jacobian matrix J(x, mu)
    param_range: Sequence of parameter values
    x_eq_func: Function to compute equilibrium x*(mu)
    """
    
    bifurcations = []
    prev_eigenvalues = None
    
    for mu in param_range:
        x_eq = x_eq_func(mu)
        J = jacobian_func(x_eq, mu)
        eigenvalues = eigvals(J)
        
        if prev_eigenvalues is not None:
            # Detect real eigenvalue crossing zero (saddle-node, transcritical, pitchfork)
            prev_real = np.sort(np.real(prev_eigenvalues))
            curr_real = np.sort(np.real(eigenvalues))
            
            for pr, cr in zip(prev_real, curr_real):
                if pr * cr < 0:  # Sign change
                    bifurcations.append({
                        'type': 'Real eigenvalue crossing',
                        'mu': mu,
                        'eigenvalues': eigenvalues
                    })
                    break
            
            # Detect Hopf bifurcation (complex eigenvalue real part crosses zero)
            complex_eigs = eigenvalues[np.abs(np.imag(eigenvalues)) > 1e-10]
            if len(complex_eigs) > 0:
                prev_complex = prev_eigenvalues[np.abs(np.imag(prev_eigenvalues)) > 1e-10]
                if len(prev_complex) > 0:
                    pr = np.real(prev_complex[0])
                    cr = np.real(complex_eigs[0])
                    if pr * cr < 0:
                        bifurcations.append({
                            'type': 'Hopf bifurcation',
                            'mu': mu,
                            'eigenvalues': eigenvalues
                        })
        
        prev_eigenvalues = eigenvalues
    
    return bifurcations

Summary

In this chapter, we systematically studied the core content of bifurcation theory:

  1. One-dimensional bifurcations:
    • Saddle-node bifurcation (equilibrium creation/annihilation)
    • Transcritical bifurcation (stability exchange)
    • Pitchfork bifurcation (symmetry breaking)
  2. Two-dimensional bifurcations:
    • Hopf bifurcation (birth of periodic solutions)
    • Global bifurcations (homoclinic, heteroclinic)
  3. Codimension and universality:
    • Bifurcation classification depends on codimension
    • Normal forms provide unified descriptions
  4. Applications:
    • Ecology (species extinction)
    • Climate science (tipping points)
    • Neuroscience (firing patterns)
    • Engineering (structural instability)
    • Economics (market crashes)
  5. Route to chaos:
    • Period-doubling cascade
    • Feigenbaum constant

Bifurcation theory provides a mathematical framework for understanding "abrupt changes" in nature. It tells us: gradual causes can produce abrupt effects— understanding bifurcation is the key to understanding critical phenomena.

Exercises

Conceptual Problems

  1. What is bifurcation? What is the special significance of the bifurcation value?

  2. Explain the difference between supercritical and subcritical bifurcations. Why are subcritical bifurcations more "dangerous"?

  3. What is the fundamental difference between Hopf bifurcation and saddle-node bifurcation?

  4. What is codimension? Why are codimension-2 bifurcations harder to observe than codimension-1 bifurcations?

Computational Problems

  1. Analyze the bifurcation behavior of the system. Find all bifurcation points and classify them.

  2. For the system

  1. Find the equilibrium points
  2. Analyze the type of bifurcation at (c) Calculate the value ofat which Hopf bifurcation occurs (if it exists)
  1. For the logistic map:
    1. Find the fixed points and analyze their stability
    2. Prove that period-doubling bifurcation occurs at (c) Find an explicit expression for the period-2 orbit
  2. Consider the logistic equation with time delay:Qualitatively analyze how the delayaffects system stability.

Analytical Problems

  1. Lake ecosystems may exist in two stable states: "clear water" and "turbid water". Use bifurcation theory to explain:
    1. Why lakes may "suddenly" shift from clear to turbid
    2. Why reducing pollutants may not restore lake clarity
    3. What is the "hysteresis effect"
  2. Explain why period-doubling cascade ultimately leads to chaos. What does the "universality" of the Feigenbaum constant mean?

Programming Problems

  1. Write a program to plot the complete bifurcation diagram for the following system:Mark all bifurcation points.

  2. Implement an AUTO-type continuation algorithm to track the equilibrium branches of, including fold-back portions.

  3. For the FitzHugh-Nagumo model, numerically calculate the location of the Hopf bifurcation point and plot the bifurcation diagram (limit cycle amplitude vs. external current).

  4. Verify the Feigenbaum constant: For the logistic map, numerically compute the period-doubling bifurcation values, calculate the ratio sequence, and observe its convergence to approximately 4.669.

References

  1. Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos. CRC Press.

  2. Kuznetsov, Y. A. (2004). Elements of Applied Bifurcation Theory. Springer.

  3. Guckenheimer, J., & Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer.

  4. Wiggins, S. (2003). Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer.

  5. Seydel, R. (2010). Practical Bifurcation and Stability Analysis. Springer.

  6. Feigenbaum, M. J. (1978). "Quantitative universality for a class of nonlinear transformations." Journal of Statistical Physics, 19(1), 25-52.

  7. May, R. M. (1976). "Simple mathematical models with very complicated dynamics." Nature, 261, 459-467.

  8. Scheffer, M., et al. (2009). "Early-warning signals for critical transitions." Nature, 461, 53-59.

Previous Chapter: ← Ordinary Differential Equations (IX): Chaos Theory and the Lorenz System

Next Chapter: → Ordinary Differential Equations (XI): Numerical Methods for Differential Equations

This article is Chapter 10 of "The World of Ordinary Differential Equations" series.

  • Post title：Ordinary Differential Equations (10): Bifurcation Theory
  • Post author：Chen Kai
  • Create time：2019-05-24 09:15:00
  • Post link：https://www.chenk.top/ode-chapter-10-bifurcation-theory/
  • Copyright Notice：All articles in this blog are licensed under BY-NC-SA unless stating additionally.
 Comments