pycont-lite 0.4.0

Lightweight matrix-free arclength continuation with bifurcation detection and branch switching.

PyCont-Lite

PyCont-Lite is a lightweight, matrix-free arclength continuation library for solving nonlinear parametric systems

$$ G(u, p) = 0 $$

with automatic bifurcation detection, branch switching, and stability analysis.

Key Features

• Matrix-free: only need to implement $G(u, p)$, no Jacobians required.
• Adaptive continuation: robust predictor–corrector with adjustable step sizes.
• Bifurcation detection: automatically localizes and classifies folds and branch points.
• Branch switching: continues all discovered branches recursively.
• Stability analysis: computes leading eigenvalues to distinguish stable vs. unstable segments.
• Lightweight plotting: quickly visualize bifurcation diagrams with plotBifurcationDiagram. Plots automatically distinguish stable vs. unstable branches.
• Structured output: branches (with stability) and events (folds, bifurcations) are returned for further analysis.”

Installation & Requirements

PyCont-Lite requires only:

• numpy
• scipy
• matplotlib

Install via PyPI (pip):

pip install pycont-lite

Quick Examples

Pitchfork Bifurcation

import numpy as np
import pycont

# Define the pitchfork function
def G(u, p):
    return u**3 - p*u

# Initial guess
u0 = np.array([-3.0])
p0 = 0.0

# Run continuation
ds_max = 0.01
ds_min = 1.e-6
ds = 0.001
n_steps = 1000
solver_parameters = {"tolerance": 1e-10}
continuation_result = pycont.arclengthContinuation(G, u0, p0, ds_min, ds_max, ds, n_steps, solver_parameters=solver_parameters)

# Plot the solution curve
pycont.plotBifurcationDiagram(continuation_result)

Fold Bifurcation

import numpy as np
import pycont

# Define the fold function
G = lambda x, r: r + x**2

# Initial guess
u0 = np.array([-5.0])
p0 = -u0[0]**2

# Run continuation
ds_max = 0.01
ds_min = 1.e-6
ds = 0.001
n_steps = 5000
solver_parameters = {"tolerance": 1e-10}
continuation_result = pycont.arclengthContinuation(G, u0, p0, ds_min, ds_max, ds, n_steps, solver_parameters=solver_parameters)

# Plot the curves
pycont.plotBifurcationDiagram(continuation_result)

Bratu: PDE Example

The nonlinear boundary value problem

$$ u_{xx} + \lambda e^u = 0, \quad u(0)=u(1)=0 $$

discretized with finite differences:

import numpy as np
import pycont

N = 101 # total number of points
x = np.linspace(0.0, 1.0, N)
dx = x[1] - x[0]

def G(u: np.ndarray, lam: float) -> np.ndarray:
    u_full = np.zeros(N, dtype=float)
    u_full[1:-1] = u
    
    u_xx = (u_full[:-2] - 2.0 * u_full[1:-1] + u_full[2:]) / (dx * dx)
    r = u_xx + lam * np.exp(u_full[1:-1])
    return r

# We know that u = 0 for lambda = 0 - otherwise we must solve G(u, lambda0) = 0.
lam0 = 0.0
u0 = np.zeros(N-2)

# Do continuation
ds_max = 0.01
ds_min = 1e-6
ds0 = 1e-4
n_steps = 2000
solver_parameters = {"tolerance": 1e-10}
continuation_result = pycont.arclengthContinuation(G, u0, lam0, ds_min, ds_max, ds0, n_steps, solver_parameters=solver_parameters)

# Plot the bifurcation diagram (lambda, max(u))
u_transform = lambda u : np.sign(u[50]) * np.max(np.abs(u))
pycont.plotBifurcationDiagram(continuation_result, u_transform=u_transform, p_label=r'$\lambda$', u_label=r'$\text{sign}(u) ||u||_{\infty}$')

This produces the classical S-shaped bifurcation curve with a fold near $\lambda \approx 3.51$.

Allen-Cahn

The nonlinear boundary value problem

$$ \varepsilon \phi_{xx} - \varepsilon^{-1} \phi \left(\phi^2 - 1\right) = 0, \quad \phi_x(-1) = \phi_x(1)=0 $$

discretized with finite differences:

import numpy as np
import pycont

N = 100
x = np.linspace(-1.0, 1.0, N)
dx = (x[-1] - x[0]) / (N-1)

def laplace_neumann(phi, dx):
    phi_ext = np.hstack([phi[1], phi, phi[-2]])
    phi_l = np.roll(phi_ext, -1)[1:-1]
    phi_r = np.roll(phi_ext,  1)[1:-1]
    return (phi_l - 2.0*phi + phi_r) / dx**2
def G(phi : np.ndarray, eps : float):
    phi_xx = laplace_neumann(phi, dx)
    rhs = eps * phi_xx - phi * (phi**2 - 1.0) / eps
    return rhs

# Initial Point
eps0 = 0.6
phi0 = np.zeros(N)

# Do continuation
tolerance = 1e-9
ds_max = 1e-2
ds_min = 1e-6
ds0 = 1e-4
n_steps = 1000
solver_parameters = {"tolerance" : tolerance, "param_min" : 0.22, "param_max" : 0.7}
continuation_result = pycont.arclengthContinuation(G, phi0, eps0, ds_min, ds_max, ds0, n_steps, solver_parameters=solver_parameters, verbosity='verbose')

# Plot the bifurcation diagram eps versus phi(x=-1)
u_transform = lambda phi: phi[0]
pycont.plotBifurcationDiagram(continuation_result, u_transform=u_transform, p_label=r'$\varepsilon$', u_label=r'$\phi(x=-1)$')

This reproduces the many bifurcation points as $\varepsilon \to 0$.

Hopf Bifurcation

Since version 0.4.0, PyCont-Lite is also able to detect Hopf bifurcations. A Hopf bifurcation occurs when two complex-conjugated eigenvalues of the Jacobian $G_u(u,p)$ cross the imaginary axis. The prototypical example is the Hopf normal form

$$ \begin{align} \mu x - y - (x^2 + y^2) x &= 0 \ x - \mu y - (x^2 + y^2) y &= 0 \ \end{align} $$

which has a Hopf point at $(x,y) = (0,0)$ with $\mu = 0$. PyCont-Lite can find this Hopf bifurcation point (see examples/NormalHopf.py):

import numpy as np
import pycont

def G(u : np.ndarray, mu : float) -> np.ndarray:
    x = u[0]; y = u[1]
    Gx = mu*x - y - (x**2 + y**2) * x
    Gy = x + mu*y - (x**2 + y**2) * y
    return np.array([Gx, Gy])
mu0 = -1.0
u0 = np.array([0.0, 0.0])

ds_max = 0.01
ds_min = 1.e-6
ds = 0.01
n_steps = 200
solver_parameters = {"tolerance": 1e-10, 'hopf_detection' : True}
continuation_result = pycont.arclengthContinuation(G, u0, mu0, ds_min, ds_max, ds, n_steps, solver_parameters)

pycont.plotBifurcationDiagram(continuation_result, p_label=r'$\mu$')

which produces the (trivial) bifurcation diagram

For now, Hopf bifurcation detection is disabled by default, so the user must supply it via 'hopf_detection' : True in the solver parameters.

Advanced Hopf Example: The Fitzhugh-Nagumo PDEs

For a more interesting system that exhibits a Hopf bifurcation, consider the Fitzhugh-Nagumo PDEs $$ \begin{align} u_{xx} + u - u^3 - v = 0 \ \delta v_{xx} + \varepsilon \left(u - a_1 v - a_0\right) = 0 \end{align} $$ with $\delta = 4$, $a_1 = 2$ and $a_0 = -0.03$. This system exhibits a Hopf bifurcation point near $\varepsilon \approx 0.018$, and a fold point near $\varepsilon \approx 0.94$. See examples/FitzhughNagumoPDEs.py for the code:

N = 100
L = 20.0
x = np.linspace(0.0, L, N)
dx = L / (N-1)

# Build the FHN objective function through finite differences
a0 = -0.03
a1 = 2.0
delta = 4.0
def G(z : np.ndarray, eps : float):
    u, v = z[:N], z[N:]
    u_xx = laplace_neumann(u, dx)
    v_xx = laplace_neumann(v, dx)
    u_rhs = u_xx + u - u**3 - v
    v_rhs = delta * v_xx + eps * (u - a1*v - a0)
    return np.concatenate((u_rhs, v_rhs))

# Do continuation.
tolerance = 1e-9
ds_max = 0.01
ds_min = 1e-6
ds0 = 1e-3
n_steps = 1000
solver_parameters = {"tolerance" : tolerance, "param_min" : 0.01, "hopf_detection" : True}
continuation_result = pycont.arclengthContinuation(G, z0, eps0, ds_min, ds_max, ds0, n_steps, solver_parameters)

# Plot the bifurcation diagram `eps` versus <u>
u_transform = lambda z: np.average(z[:N])
pycont.plotBifurcationDiagram(continuation_result, u_transform=u_transform, p_label=r'$\varepsilon$', u_label=r'$<u>$')

Note that we limit exploration to a minimal parameter value of $\varepsilon = 0.01$ for plotting purposes. PyCont-Lite detects the fold and Hopf points nicely

Plotting Bifurcation Diagrams

PyCont-Lite includes a helper function plotBifurcationDiagram for quick visualization. Stable segments are shown as solid lines and unstable segments as dashed lines, just like in AUTO/MATCONT. By default, it plots the parameter value on the horizontal axis and the transformed state variable on the vertical axis.

import pycont

# After running continuation...
pycont.plotBifurcationDiagram(continuation_result)

For multi-dimensional systems, you can specify how to reduce $u$ to a single scalar for plotting via the u_transform argument.

• Default behavior:
  • If $u$ has dimension 1 → plot $u$ directly.
  • If $u$ has dimension >1 → plot the Euclidean norm $||u||$.
• Custom transform (example: plot the maximum component of $u$):

pycont.plotBifurcationDiagram(result, u_transform=lambda u: u.max())

Solver Parameters

You can fine-tune the solver by passing a dictionary:

solver_parameters = {
    "rdiff": 1e-8,              # finite-difference step
    "nk_maxiter": 20,           # Newton-Krylov iterations
    "tolerance": 1e-10,         # nonlinear tolerance
    "bifurcation_detection": True,
    "analyze_stability": True,  # compute leading eigenvalue
    "initial_directions": "both"  # 'both', 'increase_p', 'decrease_p'
}

Verbosity Options

Control how much progress info PyCont-Lite prints during continuation. Three levels are supported:

Level	What you see
off	No progress messages (errors only).
info	One-line progress per step + event summaries (recommended default).
verbose	Solver details: Newton–Krylov iterations, step rejections, preconditioner notes.

Pass the level to arclengthContinuation via the verbosity argument. You can use a string (case-insensitive) or the enum:

from pycont import arclengthContinuation, Verbosity

# String (case-insensitive)
arclengthContinuation(G, u0, p0, ..., verbosity="info")

# Enum
arclengthContinuation(G, u0, p0, ..., verbosity=Verbosity.VERBOSE)

Default is info. See the Allen–Cahn example above for a typical verbose run.

Output Format

arclengthContinuation returns a ContinuationResult object with:

• branches: list of Branch objects (u_path, p_path, stability flag, etc.)
• events: list of Event objects (start points, folds, bifurcations)

This makes it easy to explore and plot bifurcation diagrams programmatically.

Planned Features

The following features are under active consideration for future releases:

• Hopf bifurcation detection
• Limit cycle continuation
• More stable bifurcation detection by employing two independent (but fast) tests

License

This project is licensed under the MIT License — see the LICENSE file for details.

Acknowledgement

I started this project because most continuation software is either legacy Fortran (AUTO), Matlab-only (MATCONT, COCO), or heavyweight. PyCont-Lite is meant to be modern, lightweight, and useful both for industry and academia!

For feature requests or contributions, feel free to open an issue or reach out: hannesvdc[at]gmail[dot]com.