.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/homoclinic_mtn.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_homoclinic_mtn.py: 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 :meth:`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: .. math:: \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), with the Holling-II functional responses .. math:: \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 :math:`(K, G_f, d, H_z) = (6.0, 0.0673, 0.01, 0.5)` with the saddle equilibrium at :math:`(u^\star, w^\star) = (5.7386, 0.5108)` and a truncation interval :math:`\text{PAR}(11) = 1046.18`. Continuing in :math:`(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 :math:`K \approx 6.61,\ G_f \approx 0.069` where :math:`\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 ^^^^^^^^^^ .. [1] M. Scheffer (1991) *Should we expect strange attractors behind plankton dynamics — and if so, should we bother?*, J. Plankton Res. 13, 1291-1305; reproduced as auto-07p's ``demos/mtn`` (Doedel et al.). .. [2] AUTO manual, Chapter 20 (Homoclinic Bifurcation Analysis), §20.5 (the PSI test-function table) and §20.6 (the ``mtn`` walkthrough). .. GENERATED FROM PYTHON SOURCE LINES 52-65 Step 1: build the model + ship the seed alongside ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The model is small enough to define inline with a PyRates :class:`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 :meth:`continue_homoclinic` can read it (``ISTART=1`` requires the seed on disk). .. GENERATED FROM PYTHON SOURCE LINES 65-93 .. code-block:: Python 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}) .. GENERATED FROM PYTHON SOURCE LINES 94-102 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. .. GENERATED FROM PYTHON SOURCE LINES 102-112 .. code-block:: Python ode = ODESystem.from_template( template=circuit, working_dir=str(work), auto_dir="~/PycharmProjects/auto-07p", init_cont=False, analytical_jacobian=True, auto_constants=('hom',), ) .. GENERATED FROM PYTHON SOURCE LINES 113-131 Step 3: continue the saddle-node homoclinic curve (auto-07p's mtn.1) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This is the heart of the demo. Mode-B :meth:`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. .. GENERATED FROM PYTHON SOURCE LINES 131-156 .. code-block:: Python 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" ) .. GENERATED FROM PYTHON SOURCE LINES 157-167 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. .. GENERATED FROM PYTHON SOURCE LINES 167-184 .. code-block:: Python 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())}") .. GENERATED FROM PYTHON SOURCE LINES 185-187 Step 5: visualise the saddle-node-homoclinic curve in :math:`(K, G_f)` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 187-222 .. code-block:: Python 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() .. _sphx_glr_download_auto_examples_homoclinic_mtn.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: homoclinic_mtn.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: homoclinic_mtn.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: homoclinic_mtn.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_