Saddle-Node Homoclinic Continuation: the `mtn` Demo

This gallery reproduces auto-07p’s demos/mtn homoclinic-bifurcation tutorial — M. Scheffer’s predator-prey model with a saddle-node homoclinic orbit — using PyCoBi’s ODESystem.continue_homoclinic() entry point. It’s the small, well-conditioned test case for PyCoBi’s HomCont support (auto-07p Ch. 20): a 2-D system with a single saddle-node, a pre-computed near-homoclinic orbit shipped as mtn.dat next to this script, and sign-flip detection of the PSI(15) / PSI(16) test functions that identify non-central homoclinic to saddle-node bifurcations [1].

The model is:

\[\begin{split}\dot u &= 0.5\, u\, (1 - u/K) - w \cdot \psi_Z(u) + d K, \\ \dot w &= 0.6\, \psi_Z(u)\, w - 0.15\, w - G_f\, \psi_F(w),\end{split}\]

with the Holling-II functional responses

\[\psi_Z(u) = \frac{0.4\, u}{0.6 + u}, \qquad \psi_F(w) = \frac{w^2}{w^2 + H_z^2}.\]

The seed orbit lives at \((K, G_f, d, H_z) = (6.0, 0.0673, 0.01, 0.5)\) with the saddle equilibrium at \((u^\star, w^\star) = (5.7386, 0.5108)\) and a truncation interval \(\text{PAR}(11) = 1046.18\). Continuing in \((K, G_f)\) while monitoring the PSI(15)/(16) test functions, the first run walks the curve of homoclinics to a saddle-node and flags a UZ point at \(K \approx 6.61,\ G_f \approx 0.069\) where \(\psi_{15} = 0\) (matching auto-07p’s mtn.1 step 26 to all significant figures). The UZ marks the codim-2 non-central SNIC in the sense of Nechyporenko, Ashwin & Tsaneva-Atanasova (2026, arXiv:2412.12298) — where the homoclinic orbit’s return trajectory comes in along the saddle-node’s stable manifold rather than along the central (zero-eigenvalue) direction — not the codim-1 SNIC of a periodic orbit.

References

Step 1: build the model + ship the seed alongside

The model is small enough to define inline with a PyRates OperatorTemplate — no YAML required. We rename a couple of parameters from the auto-07p Fortran (GF → Gf, HZ → Hz) because GF collides with sympy’s Galois-field constructor and the equation parser then fails to sympify the RHS.

The starting orbit lives in mtn.dat next to this script — auto-07p pre-computed it for the demo; we copy the file verbatim into the working directory before the HomCont run so continue_homoclinic() can read it (ISTART=1 requires the seed on disk).

import shutil
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
from pyrates import CircuitTemplate, NodeTemplate, OperatorTemplate

from pycobi import ODESystem

here = Path(__file__).resolve().parent
work = here  # write the c.* / .f90 / output files alongside the script
shutil.copy(here / 'mtn.dat', work / 'mtn_seed.dat')

op = OperatorTemplate(
    name='mtn_op',
    equations=[
        "u' = 0.5*u*(1.0 - u/K) - w*0.4*u/(0.6 + u) + d*K",
        "w' = 0.6*0.4*u/(0.6 + u)*w - 0.15*w - Gf*w*w/(w*w + Hz*Hz)",
    ],
    variables={
        'u': 'output(5.738626)', 'w': 'variable(0.5108401)',
        'K': 6.0, 'Gf': 0.06729762, 'd': 0.01, 'Hz': 0.5,
    },
)
node = NodeTemplate(name='mtn_node', operators=[op])
circuit = CircuitTemplate(name='mtn', nodes={'p': node})

Step 2: instantiate the ODESystem and request the ‘hom’ scenario

auto_constants=('hom',) generates a c.hom file pre-configured for auto-07p’s IPS=9 mode with sensible HomCont defaults. PyRates emits the analytical Jacobian (analytical_jacobian=True); HomCont consumes it through the same DFDU path as the equilibrium and limit-cycle scenarios.

ode = ODESystem.from_template(
    template=circuit,
    working_dir=str(work),
    auto_dir="~/PycharmProjects/auto-07p",
    init_cont=False,
    analytical_jacobian=True,
    auto_constants=('hom',),
)

Step 3: continue the saddle-node homoclinic curve (auto-07p’s mtn.1)

This is the heart of the demo. Mode-B continue_homoclinic() (omit origin, pass dat_basename) skips the LC seed-extraction step and reads mtn_seed.dat straight off disk. The PAR={...} override supplies the seed orbit’s parameter values, including PAR(11) = 1046.18 (the truncation interval) and PAR(12), PAR(13) (the saddle’s equilibrium values that IEQUIB=2 projects the unstable manifold from).

NUNSTAB = NSTAB = -1 tells HomCont to derive the eigenvalue split from NDIM; for this 2-D model it picks the saddle-node-compatible split internally. IPSI=(15, 16) selects the non-central saddle-node homoclinic test functions; auto-07p stores them at PAR(35), PAR(36). UZR={35: 0.0, 36: 0.0} flags zero crossings of either as a UZ label — that’s where the saddle-node homoclinic codim-2 phenomenon happens.

sols_fwd, _ = ode.continue_homoclinic(
    dat_basename='mtn_seed',
    ICP=['K', 'Gf'],
    NUNSTAB=-1, NSTAB=-1,
    IEQUIB=2, ITWIST=0, IPSI=(15, 16),
    name='snhom_fwd',
    PAR={11: 1046.178, 12: 5.738626, 13: 0.5108401},
    NTST=100, NCOL=4, IAD=1, ISP=0, ILP=0,
    NMX=30, NPR=10, IID=2, ITMX=10, ITNW=7, NWTN=3,
    DS=0.01, DSMIN=0.001, DSMAX=0.05,
    UZR={35: 0.0, 36: 0.0},
)

print(f"\nForward HomCont branch ({len(sols_fwd)} points)")
print(f"  bifurcations: {dict(sols_fwd['bifurcation'].value_counts())}")
sel = sols_fwd[[c for c in sols_fwd.columns
                if (isinstance(c, tuple) and c[0] in
                    ('K', 'Gf', 'PAR(35)', 'PAR(36)', 'bifurcation'))]]
print(sel.to_string())
print(
    "Reference (auto-07p mtn.1, step 26):  K = 6.61046, Gf = 0.069325, "
    "PSI(15) = 5.129e-9 → UZ"
)

Step 4: continue in the opposite direction (auto-07p’s mtn.2)

The reverse-direction half of the saddle-node-homoclinic curve uses a negative DS so auto-07p walks toward smaller K instead. The resulting branch picks up the other SNIC point — where the second PSI test function (PSI(16) = non-central homoclinic to saddle-node in the unstable manifold) flips sign. We restart from the original seed rather than chaining off snhom_fwd because the .dat file is the cleanest entry point in PyCoBi’s mode-B workflow.

sols_rev, _ = ode.continue_homoclinic(
    dat_basename='mtn_seed',
    ICP=['K', 'Gf'],
    NUNSTAB=-1, NSTAB=-1,
    IEQUIB=2, ITWIST=0, IPSI=(15, 16),
    name='snhom_rev',
    PAR={11: 1046.178, 12: 5.738626, 13: 0.5108401},
    NTST=50, NCOL=4, IAD=1, ISP=0, ILP=0,
    NMX=40, NPR=10, IID=2, ITMX=10, ITNW=7, NWTN=3,
    DS=-0.01, DSMIN=0.001, DSMAX=0.05,
    UZR={35: 0.0, 36: 0.0},
)

print(f"\nReverse HomCont branch ({len(sols_rev)} points)")
print(f"  bifurcations: {dict(sols_rev['bifurcation'].value_counts())}")

Step 5: visualise the saddle-node-homoclinic curve in \((K, G_f)\)

def _col(df, name):
    """Pull a column out of a possibly-MultiIndexed PyCoBi summary."""
    for c in df.columns:
        if (isinstance(c, tuple) and c[0] == name) or c == name:
            return df[c].to_numpy(dtype=float)
    raise KeyError(name)

fig, ax = plt.subplots(figsize=(7, 5))
for sols, lbl, col in [(sols_fwd, 'forward', '#1F77B4'),
                       (sols_rev, 'reverse', '#FF7F0E')]:
    K = _col(sols, 'K')
    Gf = _col(sols, 'Gf')
    ax.plot(K, Gf, '-', color=col, lw=1.5, label=f'SN-homoclinic ({lbl})')

# Overlay each UZ — auto-07p flags these where the PSI test functions
# pass through zero, i.e. where the homoclinic touches the saddle-node
# line (the codim-2 SNIC point).
for sols, col in [(sols_fwd, '#1F77B4'), (sols_rev, '#FF7F0E')]:
    bif = sols[('bifurcation', '')] if isinstance(sols.columns, type(sols.columns)) \
        and (('bifurcation', '') in sols.columns) else sols['bifurcation']
    K = _col(sols, 'K')
    Gf = _col(sols, 'Gf')
    for i, b in enumerate(bif):
        if str(b).strip() == 'UZ':
            ax.scatter(K[i], Gf[i], marker='*', s=200, c=col,
                       edgecolor='k', linewidth=0.8, zorder=10,
                       label='SNIC (PSI=0)' if i == 0 else None)

ax.set_xlabel('K (carrying capacity)')
ax.set_ylabel(r'$G_f$ (predator gain)')
ax.set_title('Scheffer predator-prey: saddle-node homoclinic curve')
ax.legend(loc='best')
plt.tight_layout()
plt.show()

Gallery generated by Sphinx-Gallery

Saddle-Node Homoclinic Continuation: the mtn Demo