From 2D dynamics to 3D dynamics - The Poincaré–Bendixson Theorem, and chaos.sagews

What sorts of long-term behavior can you
possibly get from a single (1D)
differential equation?

What sorts of long-term behavior can you
possibly get from a single (1D)
differential equation?

Answer:

• Constant solution (start on an equilibrium point, stay there forever)

• Approach a stable equilibrium point

• Grow without bound (increase to $\infty$ or decrease to $-\infty$)

What sorts of long-term behavior can you
possibly have in a 2D system of
differential equations?

state_vars = list(var("X, Y"))
system = (
    4 - Y,
    2 - X,
)
field(X, Y) = system
t_range = srange(0, 10, 0.1)
ics = {
    (2, 4): t_range,
}
p1 = plot_phase_portrait(field, (X, -2, 4), (Y, -1, 5), ics,
                         nullcline_colors=(None, None), aspect_ratio=1)
p = p1 + text("Constant solution (single point trajectory)", (1, 6), fontsize=12)
p.show(axes_labels=("$X$", "$Y$"))

state_vars = list(var("X, Y"))
system = (
    -X + 0.4*Y,
    -2*Y - X,
)
field(X, Y) = system
t_range = srange(0, 100, 0.1)
ics = {
    (2, 2): t_range,
}
p2 = plot_phase_portrait(field, (X, -3, 3), (Y, -3, 3), ics,
                         nullcline_colors=(None, None), aspect_ratio=1)
p = p2 + text("Trajectory approaches stable eq. point", (0, 4), fontsize=12)
p.show(axes_labels=("$X$", "$Y$"))

state_vars = list(var("X, Y"))
system = (
    X - 0.5*Y,
    2*Y + X,
)
field(X, Y) = system
t_range = srange(0, 15.55, 0.1)
ics = {
    (0.1, -0.1): t_range,
}
p3 = plot_phase_portrait(field, (X, -1e9, 1e9), (Y, -1e9, 1e9), ics,
                         nullcline_colors=(None, None), aspect_ratio=1)
p = p3 + text('Trajectory is unbounded ("goes to infinity")', (0, 1.35E9), fontsize=12)
p.show(axes_labels=("$X$", "$Y$"))

state_vars = list(var("X, Y"))
system = (
    Y,
    -X - (Y^3 - 0.1*Y),
)
field(X, Y) = system
t_range = srange(0, 10, 0.1)
ics = {
    (0.365, 0): t_range,
}
p4 = plot_phase_portrait(field, (X, -1, 1), (Y, -1, 1), ics,
                         nullcline_colors=(None, None), aspect_ratio=1)
p = p4 + text("Trajectory is periodic (closed loop orbit)", (0, 1.35), fontsize=12)
p.show(axes_labels=("$X$", "$Y$"))

state_vars = list(var("X, Y"))
system = (
    Y,
    -X - (Y^3 - 0.1*Y),
)
field(X, Y) = system
t_range = srange(0, 100, 0.1)
ics = {
    (1, -0.75): t_range,
}
p5 = plot_phase_portrait(field, (X, -1, 1), (Y, -1, 1), ics,
                         nullcline_colors=(None, None), aspect_ratio=1)
p = p5 + text("Trajectory approaches a periodic cycle (limit cycle attractor)", (0, 1.35), fontsize=12)
p.show(axes_labels=("$X$", "$Y$"))

state_vars = list(var("X, Y"))
system = (
    (Y - 0.25*X)*(X^2 - 1),
    X*(1 - Y^2),
)
field(X, Y) = system
t_range = srange(0, 100, 0.1)
ics = {
    (0.02, 0.02): t_range,
}
p6 = plot_phase_portrait(field, (X, -2, 2), (Y, -2, 2), ics,
                         aspect_ratio=1)
p = p6 + text('Multiple saddle points form a "polygon" of other', (0, 2.9), fontsize=12)
p += text('trajectories, and this trajectory is trapped inside!', (0, 2.7), fontsize=12)
p.show(axes_labels=("$X$", "$Y$"))

These are all of the long-term behaviors
that can happen to a trajectory in 2D!

(This fact is called the Poincaré–Bendixson Theorem.)

ga = graphics_array(((p1, p2, p3), (p4, p5, p6)))
ga.show()

In 3D (and higher dimensions), there's a
new kind of behavior that does not fit
any of these patterns...

state_vars = list(var("X, Y, Z"))
system = (
    X*(1 - X) - 5/3*X/(1/3 + X)*Y,
    5/3*X/(1/3 + X)*Y - 0.4*Y - 0.05*Y/(0.5 + Y)*Z,
    0.05*Y/(0.5 + Y)*Z - 0.01*Z,
)
field(X, Y, Z) = system
t_range = srange(0, 10000, 0.1)
solution = desolve_odeint(field, (1, 1, 1), t_range, state_vars)
start = 6000/0.1
#p = axes3d((X, 0, 1), (Y, 0, 0.7), (Z, 0, 12), grid=False)
p = Graphics()
p += line3d(((0,0,7), (1,  0, 7)), color="black") + text3d("X", (1.1,0,7), color="black")
p += line3d(((0,0,7), (0,0.5, 7)), color="black") + text3d("Y", (0,0.6,7), color="black")
p += line3d(((0,0,7), (0,  0,11)), color="black") + text3d("Z", (0,0, 11.1), color="black")
p += list_plot(solution[start::10,:], plotjoined=True, color="red")
p.show(frame=False, aspect_ratio=(3,6,1))

3D rendering not yet implemented

start = int(9000/0.1)
p  = list_plot(zip(t_range[start:], solution[start:,0]), plotjoined=True, color="red",   legend_label="$X$")
p += list_plot(zip(t_range[start:], solution[start:,1]), plotjoined=True, color="green", legend_label="$Y$")
p += list_plot(zip(t_range[start:], solution[start:,2]), plotjoined=True, color="blue",  legend_label="$Z$")
show(p, ymin=0, figsize=5, axes_labels=("$t$", "Populations"))

It never precisely repeats itself, so it's not oscillating. But it also does not go to equilibrium, nor does it go to infinity.

This is called chaotic behavior, and this will be the topic of Chapter 5.