Tutorial: Monodromy Matrix and Orbit Shift on an Analytic X-type Saddle Cycle#
This notebook demonstrates the core pyna monodromy workflow on an analytically constructed saddle cycle (X-point of a Poincaré map) embedded in a 3-D toroidal geometry.
What you will learn:
How to define an analytic field with a known m=3/n=1 X-cycle in cylindrical (R, Z, φ) coordinates
How
evolve_DPm_along_cycleintegrates the variational equationdDX_pol/dφ = A · DX_polto obtain the monodromy matrixDPmWhy the DPm evolution (commutator equation) needs
DPm(φ₀) = DX_pol(φ_end), not the identity ?and what that means geometricallyHow to compute the linear orbit shift under a perturbation field using
orbit_shift_under_perturbationHow to compare the analytic formula with the pyna API result
1. Analytic field construction#
We construct a field whose exact m=3 X-cycle is known analytically.
1.1 The cycle and its eigendirections#
The cycle lives on the ellipse
with \(\iota = 1/3\), so the orbit closes after \(m=3\) toroidal turns.
The Poincaré map at the X-point has eigenvalues \(\lambda_u = e^{+1/5}\) (unstable) and \(\lambda_s = e^{-1/5}\) (stable). The eigendirections rotate along the orbit following a prescribed angle \(\theta(\varphi)\).
[1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# ── Cycle parameters ────────────────────────────────────────────────────────
Rax, Zax = 1.0, 0.0
Rell, Zell = 0.3, 0.5
phi0_cycle = 0.0
BPhiax = 2.5
m = 3 # toroidal period (orbit closes after m turns)
n = 1
iota = n / m # rotational transform
# Eigenvalues of the m-turn Poincaré map
lam_u = np.e ** (1/5) # unstable multiplier
lam_s = np.e ** (-1/5) # stable multiplier
print(f"λ_u = {lam_u:.6f}, λ_s = {lam_s:.6f}, product = {lam_u * lam_s:.10f} (should be 1)")
# Rotating eigendirections along the orbit
def theta_u(phi):
"""Angle of unstable eigenvector at toroidal position phi."""
return phi/3 + np.pi/2 + np.pi/9
def theta_s(phi):
"""Angle of stable eigenvector at toroidal position phi."""
return phi/3 + np.pi/2 - np.pi/9
dtheta_dphi = 1.0/3 # same rotation rate as the orbit
# ── The cycle position ──────────────────────────────────────────────────────
def cycle_RZ(phi):
"""Return (R, Z) of the X-cycle at toroidal angle phi."""
return np.array([
Rell * np.cos(iota * phi + phi0_cycle) + Rax,
Zell * np.sin(iota * phi + phi0_cycle) + Zax,
])
# ── Toroidal field (axisymmetric model: R·B_φ = const) ─────────────────────
def BPhi(phi, RZ):
"""Toroidal field B_φ(R) = B_φ^ax · R_ax / R (free-force-balance model)."""
return BPhiax * Rax / RZ[0]
print("Cycle at phi=0:", cycle_RZ(0.0))
print("Cycle at phi=2π:", cycle_RZ(2*np.pi))
print("Cycle at phi=6π (closes):", cycle_RZ(6*np.pi))
λ_u = 1.221403, λ_s = 0.818731, product = 1.0000000000 (should be 1)
Cycle at phi=0: [1.3 0. ]
Cycle at phi=2π: [0.85 0.4330127]
Cycle at phi=6π (closes): [ 1.3000000e+00 -1.2246468e-16]
[2]:
# ── Poloidal field on the cycle ──────────────────────────────────────────────
#
# The poloidal field on the exact cycle is just B_pol = (dR_c/dl, dZ_c/dl)
# re-parameterised by phi. We compute the φ-derivatives:
# dR_c/dφ = -ι R_ell sin(ι φ + φ₀)
# dZ_c/dφ = +ι Z_ell cos(ι φ + φ₀)
def BR0_BZ0_on_cycle(phi):
"""Return (B_R, B_Z) of the unperturbed field on the X-cycle."""
Rc = cycle_RZ(phi)[0]
Bp = BPhi(phi, cycle_RZ(phi))
dRc = -iota * Rell * np.sin(iota * phi + phi0_cycle)
dZc = iota * Zell * np.cos(iota * phi + phi0_cycle)
return np.array([dRc / Rc * Bp, dZc / Rc * Bp])
# ── Full field B(R, Z, φ) with exact X-cycle ────────────────────────────────
#
# We construct B_pol(X, φ) so that:
# 1. B_pol(X_c(φ), φ) = B_pol^0(φ) (matches cycle tangent)
# 2. The Jacobian A(φ) of g = R·B_pol/B_φ along X_c has the prescribed
# eigenstructure (λ_u, λ_s with rotating eigenvectors θ_u, θ_s).
#
# Construction: let V(φ) = [e_u | e_s] be the eigenvector matrix and
# Λ = diag(log|λ_u|, log|λ_s|) / (2mπ) (per-φ growth rates).
# Then A(φ) = V Λ V^{-1} + dV/dφ · V^{-1} (accounts for rotating frame).
# The field is B_pol(X, φ) = B_pol^0(φ) + R·A(φ)·(X - X_c(φ)) / R_c
def _eigvec_matrix(phi):
"""2x2 matrix of eigenvectors [e_u | e_s] at phi."""
return np.array([
[np.cos(theta_u(phi)), np.cos(theta_s(phi))],
[np.sin(theta_u(phi)), np.sin(theta_s(phi))],
])
def _A_matrix(phi):
"""A(φ) = Jacobian of g = R·B_pol/B_φ w.r.t. (R,Z) along cycle."""
V = _eigvec_matrix(phi)
# per-radian growth rates
mu_u = np.log(np.abs(lam_u)) / (2 * m * np.pi)
mu_s = np.log(np.abs(lam_s)) / (2 * m * np.pi)
Lam = np.diag([mu_u, mu_s])
# rotating-frame correction: dV/dφ · V^{-1}
dV = np.array([
[-dtheta_dphi * np.sin(theta_u(phi)), -dtheta_dphi * np.sin(theta_s(phi))],
[ dtheta_dphi * np.cos(theta_u(phi)), dtheta_dphi * np.cos(theta_s(phi))],
])
return V @ Lam @ np.linalg.inv(V) + dV @ np.linalg.inv(V)
def field_g_pol(phi, X_pol):
"""φ-parameterised field direction g = R·B_pol/B_φ (units: m)."""
Xc = cycle_RZ(phi)
g0 = Xc # Xc gives the cycle tangent when dividing by R·B_phi/R_B_phi^ax
# Actually g = dX_c/dphi + A·(X - X_c)
dXc = np.array([
-iota * Rell * np.sin(iota * phi + phi0_cycle),
iota * Zell * np.cos(iota * phi + phi0_cycle),
])
A = _A_matrix(phi)
return dXc + A @ (X_pol - Xc)
# Verify: g on the cycle equals dX_c/dphi
phi_test = 0.7
Xc_test = cycle_RZ(phi_test)
g_on_cycle = field_g_pol(phi_test, Xc_test)
dXc_expected = np.array([
-iota * Rell * np.sin(iota * phi_test + phi0_cycle),
iota * Zell * np.cos(iota * phi_test + phi0_cycle),
])
print("g on cycle :", g_on_cycle)
print("dX_c/dφ :", dXc_expected)
print("Match:", np.allclose(g_on_cycle, dXc_expected))
g on cycle : [-0.02312218 0.16215018]
dX_c/dφ : [-0.02312218 0.16215018]
Match: True
2. Wrapping as a pyna field function#
pyna expects field functions with signature field_func(rzphi) ?(dR/dl, dZ/dl, dφ/dl), i.e. arc-length parameterised. We convert from the φ-parameterisation using \(d\varphi/dl = B_\varphi / (R |B|)\).
[3]:
def pyna_field(rzphi):
"""pyna-compatible field function: rzphi=(R,Z,φ) ?(dR/dl, dZ/dl, dφ/dl)."""
R, Z, phi = rzphi[0], rzphi[1], rzphi[2]
X_pol = np.array([R, Z])
g = field_g_pol(phi, X_pol) # (dR/dφ, dZ/dφ)
Bp = BPhi(phi, X_pol) # B_φ at (R, Z)
# dφ/dl = B_φ / (R · |B|), |B|² = B_R² + B_Z² + B_φ²
# B_R = g[0] · B_φ / R, B_Z = g[1] · B_φ / R
# So |B|² = B_φ² · (g²/R² + 1) ? dphi/dl = 1/sqrt(R² + |g|²)
dphi_dl = 1.0 / np.sqrt(R**2 + np.dot(g, g))
dR_dl = g[0] * dphi_dl
dZ_dl = g[1] * dphi_dl
return np.array([dR_dl, dZ_dl, dphi_dl])
# Quick sanity: integrate the cycle once and check it closes
def rhs_phi(phi, XZ):
f = pyna_field(np.array([XZ[0], XZ[1], phi]))
dphi_dl = f[2]
return np.array([f[0] / dphi_dl, f[1] / dphi_dl]) # dX/dphi
X0 = cycle_RZ(0.0)
sol = solve_ivp(rhs_phi, [0, 2 * m * np.pi], X0, rtol=1e-10, atol=1e-12, max_step=0.01)
X_end = sol.y[:, -1]
print(f"Cycle start: {X0}")
print(f"Cycle end : {X_end}")
print(f"Closure error: {np.linalg.norm(X_end - X0):.2e} (should be < 1e-6)")
Cycle start: [1.3 0. ]
Cycle end : [1.30000000e+00 3.20923843e-17]
Closure error: 4.45e-16 (should be < 1e-6)
3. Computing the monodromy matrix with pyna.topo.evolve_DPm_along_cycle#
The function evolve_DPm_along_cycle integrates two coupled equations along the orbit:
Phase 1 ?variational equation for DX_pol:
After one full period (\(\varphi_0 \to \varphi_0 + 2m\pi\)), we get \(\mathtt{DPm} = \mathtt{DX\_pol}(\varphi_{\rm end})\) ?the monodromy matrix.
Phase 2 ?commutator equation for the DPm evolution:
Note the initial condition is not the identity ?it is the monodromy matrix just computed in Phase 1. This ensures that DPm(φ) tracks how the full-period map Jacobian would look if we started the period at φ instead of φ₀.
[4]:
from pyna.topo.monodromy import evolve_DPm_along_cycle
from pyna.topo.toroidal_cycle import ToroidalPeriodicOrbitTrace as PeriodicOrbit
# Build a PeriodicOrbit object that pyna needs
rzphi0 = np.array([cycle_RZ(0.0)[0], cycle_RZ(0.0)[1], 0.0])
orbit = PeriodicOrbit(
rzphi0=rzphi0,
period_m=m,
trajectory=np.zeros((1, 3)), # placeholder; evolve_DPm_along_cycle re-integrates
DPm=np.eye(2), # placeholder
)
# Run monodromy computation
mono = evolve_DPm_along_cycle(
field_func=pyna_field,
orbit=orbit,
dt_output=2 * np.pi / 200, # ~200 output points per turn
rtol=1e-10, atol=1e-12,
)
print("DPm (monodromy matrix):")
print(mono.DPm)
print()
print("Eigenvalues :", mono.eigenvalues)
print("Stability index (Tr/2) :", mono.stability_index)
print("Greene residue :", mono.Greene_residue)
DPm (monodromy matrix):
[[ 1.02006674 -0.07328026]
[-0.55316615 1.02006674]]
Eigenvalues : [0.81873081 1.22140267]
Stability index (Tr/2) : 1.0200667408183737
Greene residue : -0.010033370409186837
[5]:
# ── Analytic verification ────────────────────────────────────────────────────
#
# The analytic DPm at phi=0 should be V(2mπ) · diag(λ_u, λ_s) · V^{-1}(0)
# where V(φ) is the eigenvector matrix.
V_end = _eigvec_matrix(2 * m * np.pi)
V_start = _eigvec_matrix(0.0)
DPm_analytic = V_end @ np.diag([lam_u, lam_s]) @ np.linalg.inv(V_start)
print("DPm analytic:")
print(DPm_analytic)
print()
print("DPm from pyna:")
print(mono.DPm)
print()
print("Max error:", np.max(np.abs(mono.DPm - DPm_analytic)))
DPm analytic:
[[ 1.02006676 -0.07328031]
[-0.55316612 1.02006676]]
DPm from pyna:
[[ 1.02006674 -0.07328026]
[-0.55316615 1.02006674]]
Max error: 5.4677180783002655e-08
4. DX_pol along the orbit and the DPm initial condition#
Here we visualise DX_pol(φ) ?the state transition matrix at each φ ?and explain why the Phase-2 DPm equation starts from DX_pol(φ_end) rather than I.
Geometric meaning: DPm(φ) is the monodromy of the “shifted” orbit that starts at φ instead of φ₀. Because the full-period map changes when you shift the Poincaré section, DPm(φ) ?DX_pol(φ_end) in general ?it evolves via the commutator equation, and its initial value must equal DX_pol(φ_end) to be consistent with the definition.
[6]:
phi_arr = mono.phi_arr
fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True)
fig.suptitle("DX_pol(φ) components along the m=3 X-cycle", fontsize=14)
labels = [["DX_pol[0,0]", "DX_pol[0,1]"], ["DX_pol[1,0]", "DX_pol[1,1]"]]
for i in range(2):
for j in range(2):
axes[i, j].plot(phi_arr / np.pi, mono.DX_pol_arr[:, i, j], color='steelblue')
axes[i, j].set_ylabel(labels[i][j])
axes[i, j].axvline(0, color='k', lw=0.5)
axes[i, j].axvline(2 * m, color='gray', lw=0.5, linestyle='--', label='φ_end')
# Mark the initial value of DPm (= DX_pol at phi_end)
axes[i, j].scatter([2 * m], [mono.DX_pol_arr[-1, i, j]],
color='red', zorder=5, label='DPm IC')
axes[0, 0].legend(fontsize=9)
for ax in axes[1]:
ax.set_xlabel("φ / π")
plt.tight_layout()
plt.savefig("DX_pol_components.png", dpi=120)
plt.show()
print("DX_pol(φ_end) = DPm_init (Phase-2 IC):")
print(mono.DX_pol_arr[-1])
print()
print("DPm (Phase-1 result = monodromy):")
print(mono.DPm)
print()
print("They are the same:", np.allclose(mono.DX_pol_arr[-1], mono.DPm))
DX_pol(φ_end) = DPm_init (Phase-2 IC):
[[ 1.02006674 -0.07328026]
[-0.55316615 1.02006674]]
DPm (Phase-1 result = monodromy):
[[ 1.02006674 -0.07328026]
[-0.55316615 1.02006674]]
They are the same: True
5. Linear orbit shift under a perturbation field#
Now we apply a small perturbation δB and compute the linear shift of the X-cycle using orbit_shift_under_perturbation.
The shift δX(φ) satisfies:
with the periodic boundary condition δX(φ₀ + 2mπ) = δX(φ₀), giving
We use a simple analytic perturbation: a uniform vertical shift δB_Z = ε.
[7]:
from pyna.topo.monodromy import orbit_shift_under_perturbation
epsilon = 0.01 # perturbation amplitude
def delta_field(rzphi):
"""Perturbation: uniform δB_Z = ε (arc-length parameterised)."""
R, Z, phi = rzphi[0], rzphi[1], rzphi[2]
X_pol = np.array([R, Z])
Bp = BPhi(phi, X_pol)
# dφ/dl from unperturbed field
g = field_g_pol(phi, X_pol)
dphi_dl = 1.0 / np.sqrt(R**2 + np.dot(g, g))
# δB_Z = ε ?δ(dZ/dl) = ε · |B| / |B| = ε / |B| · |B| = ε dphi_dl · R / Bp · Bp = ...
# Simple: δdZ/dl = ε · dphi_dl / dphi_dl = ε (unit B)
return np.array([0.0, epsilon * dphi_dl, 0.0])
delta_X = orbit_shift_under_perturbation(
field_func=pyna_field,
delta_field_func=delta_field,
orbit=orbit,
monodromy_analysis=mono,
)
print(f"Orbit shift at φ=0: δR={delta_X[0, 0]:.6f} m, δZ={delta_X[0, 1]:.6f} m")
print(f"Periodicity check: |δX(end) - δX(start)| = "
f"{np.linalg.norm(delta_X[-1] - delta_X[0]):.2e} (should be < 1e-6)")
Orbit shift at φ=0: δR=-0.029622 m, δZ=0.000000 m
Periodicity check: |δX(end) - δX(start)| = 1.36e-15 (should be < 1e-6)
[8]:
# ── Visualise the orbit shift in 3D ─────────────────────────────────────────
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
Rc = mono.trajectory[:, 0]
Zc = mono.trajectory[:, 1]
phi_plot = mono.phi_arr
# Unperturbed cycle (xyz)
ax.plot(Rc * np.cos(phi_plot), Rc * np.sin(phi_plot), Zc,
color='black', linewidth=2, label='Unperturbed X-cycle')
# Shifted cycle
R_shifted = Rc + delta_X[:, 0]
Z_shifted = Zc + delta_X[:, 1]
ax.plot(R_shifted * np.cos(phi_plot), R_shifted * np.sin(phi_plot), Z_shifted,
color='tomato', linewidth=2, linestyle='--', label=f'Shifted cycle (ε={epsilon})')
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_zlabel('Z [m]')
ax.legend()
ax.set_title('Linear orbit shift under uniform δB_Z perturbation')
plt.tight_layout()
plt.savefig("orbit_shift_3d.png", dpi=120)
plt.show()
6. Key takeaways#
Concept |
Symbol |
pyna API |
|---|---|---|
Variational matrix along orbit |
|
|
Full-period Poincaré map Jacobian |
|
|
Phase-2 initial condition |
|
automatic in |
Eigenvalues of DPm |
|
|
Stability index |
|
|
Greene residue |
|
|
Linear orbit shift |
|
|
_pol makes clear this matrix lives in cylindrical (polar) (R, Z) space, not Cartesian. The variable name carries physics meaning at a glance.DPm = Derivative of Poincaré map (m-turn period). A single letter M gives no hint of this rich structure.