Continuation and bifurcations¶

Branches of solutions followed through folds by pseudo-arclength, with branch points, Hopf points, cycles and codimension-two curves located on the way.

Pseudo-arclength continuation follows a curve of solutions \(G(x, p) = 0\) by stepping along its tangent and correcting with Newton on the bordered system, so that folds are passed rather than hit. Events (folds, branch points, Hopf points) are located on the curve by bisection on test functions.

A nonlinear boundary-value problem¶

The Bratu problem \(u'' + \lambda e^u = 0\), \(u(\pm 1) = 0\), written natively and compiled at a declared resolution:

import numpy as np

import conatus as cn
from conatus.expert import ContinuationSettings, NewtonSolver, NLBVPContinuation

m = cn.Model(cn.Box(x=cn.walls(-1, 1)), name="Bratu", lam=0.3)
m.scalar("u")
m.equations("dx(dx(u)) + lam*exp(u) = 0")
m.at("walls", u=0)
compiled = m.compile(48, dealias=1)        # a nonlinear boundary-value problem, N = 48
problem, u = compiled.problem, compiled.fields["u"]

dealias=1 is declared here: exp(u) is not a polynomial nonlinearity, so no padding makes it exact. Newton works on the exact Jacobian of the discretised problem:

newton = NewtonSolver(problem, tol=1e-13, damping=1.0, max_iterations=20,
                      params={"lam": 0.3})
sol = newton.solve({u: lambda x: 0 * x})                # initial guess u = 0
print(sol.explain())                                    # residual per iteration, orders

Following the branch through the fold¶

settings = ContinuationSettings(ds=0.2, ds_min=1e-8, ds_max=0.6, adapt=1.5,
                                fast_iterations=3, tolerance=1e-12, max_newton=15,
                                max_points=12, dp_fd=1e-6)
branch = NLBVPContinuation(newton, "lam", settings).run(sol.coeffs, 0.3, +1,
                                                        p_bounds=(0.3, 1.0))
(fold,) = branch.folds()
print(branch.explain())
print(f"fold at lam = {fold.p:.12f}")

Every step-size rule is in ContinuationSettings: initial, minimum and maximum steps, adaptation factor, Newton budget and tolerance, the finite-difference step of \(\partial_p G\) (or None when exact derivatives exist).

Bifurcations of a dynamical system¶

For a system of ODEs, or a PDE reduced to one by ir_system, Conatus computes the exact Jacobian, \(\partial_p G\) and the multilinear forms used by the normal forms.

from conatus.expert import continue_equilibria, ir_system, ode_problem

lorenz = ode_problem(["x", "y", "z"], {"sigma": 10.0, "r": 0.5, "b": 8 / 3},
                     lambda X, P: [P["sigma"] * (X["y"] - X["x"]),
                                   P["r"] * X["x"] - X["y"] - X["x"] * X["z"],
                                   X["x"] * X["y"] - P["b"] * X["z"]], name="Lorenz")
L = ir_system(lorenz, ["r"])            # G(x; r), exact Jacobian, dG/dr, multilinear forms
cfg = ContinuationSettings(ds=0.1, ds_min=1e-9, ds_max=2.0, adapt=1.3, fast_iterations=4,
                           tolerance=1e-11, max_newton=12, max_points=60, dp_fd=None)
origin = continue_equilibria(L, np.zeros(3), "r", cfg, marginal_tol=1e-9, bounds=(0, 3))
(bp,) = origin.events_of("branch point")

Branch switching, Hopf points and cycles¶

from conatus.expert import CollocationMesh, branch_switch, normal_form

S = np.diag([-1.0, -1.0, 1.0])          # (x, y, z) -> (-x, -y, z), a symmetry of Lorenz
cplus = branch_switch(origin, bp, cfg, step=0.05, symmetry=S, bounds=(0, 30))
(hopf,) = cplus.events_of("hopf")
hf = normal_form(L, hopf.point.z[:-1], kind="hopf", par={"r": hopf.point.z[-1]},
                 parameter="r")
print(hf.explain())                     # first Lyapunov coefficient: sub- or supercritical

cfo = ContinuationSettings(ds=0.2, ds_min=1e-9, ds_max=1.0, adapt=1.3, fast_iterations=4,
                           tolerance=1e-9, max_newton=12, max_points=12, dp_fd=None)
cycles = branch_switch(cplus, hopf, cfo, amplitude=0.2,
                       mesh=CollocationMesh.uniform(40, 4), tests=("lpc",))
print(cycles.explain())                 # periodic orbits born at the Hopf point

At a symmetry-breaking pitchfork the symmetry is declared (symmetry=S) and the branch point is classified. At a Hopf point the first Lyapunov coefficient decides between sub- and supercritical; the cycles born there are continued by orthogonal collocation on a declared mesh. A branch of cycles can also start from any orbit given as a collocation state z0 (the mesh values, the period and the parameter): continue_orbits(PeriodicOrbitProblem(system, "r", mesh), z0, settings), with the same tests ("lpc", "pd", "ns", "bpc") and explain().

Codimension two¶

from conatus.expert import continue_codim2

bt = ir_system(ode_problem(["x", "y"], {"b1": -0.5, "b2": -1.0},
                           lambda X, P: [X["y"], P["b1"] + P["b2"] * X["x"]
                                         + X["x"] * X["x"] - X["x"] * X["y"]],
                           name="Bogdanov-Takens"), ["b1", "b2"])
cfg2 = ContinuationSettings(ds=0.05, ds_min=1e-9, ds_max=0.3, adapt=1.3, fast_iterations=3,
                            tolerance=1e-11, max_newton=12, max_points=40, dp_fd=None)
eq = continue_equilibria(bt, np.array([(1 - np.sqrt(3)) / 2, 0.0]), "b1", cfg2,
                         marginal_tol=1e-9, bounds=(-2, 2))
fold_pt = eq.events_of("fold")[0]
folds = continue_codim2(bt, fold_pt.point.z[:-1], "fold", ("b1", "b2"), cfg2,
                        par={"b1": fold_pt.point.z[-1]}, bounds={"b2": (-1.5, 1.5)})
print(folds.explain())                  # the fold curve b1 = b2^2 / 4 in the (b1, b2) plane

Fold and Hopf curves are continued in two parameters; codimension-two points on them (Bogdanov–Takens, cusp, generalized Hopf) are reported as events.

Amplitude equations (amplitude_equation, ginzburg_landau) reduce a PDE near onset to its Landau or Ginzburg–Landau form; tutorial 10. Bifurcations: branch switching, cycles, codimension 2, normal forms checks them on Swift–Hohenberg.