Übersetzung möglicherweise veraltet Die englische Quellseite ist neuer; diese deutsche Seite ist möglicherweise noch nicht vollständig synchronisiert.
Englische Quelle öffnen

Insel-Jacobian-Analyse: periodische Orbits und Monodromiematrizen#

Dieses Notebook demonstriert die vollständige Jacobian-/Monodromieanalyse für magnetische Inselketten in einem analytischen Stellarator-Gleichgewicht, portiert aus einem Julia-basierten Workflow.

Wissenschaftlicher Hintergrund#

Ein periodischer Feldlinienorbit oder Zyklus ist eine Feldlinie, die nach genau \(m\) toroidalen Umläufen zu ihrem Startpunkt zurückkehrt: \(X(\phi + 2\pi m) = X(\phi)\).

Das sind Fixpunkte der \(m\)-ten Iteration der Poincaré-Karte \(P^m\). Für eine Resonanz \(q = m/n\) (mit \(q = B_\phi r / (B_{pol} R)\)) beträgt die Orbitperiode \(m\) Umläufe.

Die Jacobian-Matrix \(DX(\phi)\) entwickelt sich gemäß:

\[\frac{dDX}{d\phi} = A(r(\phi), z(\phi), \phi) \cdot DX, \quad DX(\phi_0) = I\]

wobei \(A_{ij} = \partial(R B_{pol,i}/B_\phi) / \partial x_j\).

Die Monodromiematrix \(M = DP^m(\phi_0)\) liefert die linearisierte \(m\)-Umlauf-Karte. Eigenwerte von \(M\):

  • \(|\lambda| = 1\): elliptisch (O-Punkt, stabil)

  • \(|\lambda| > 1\): hyperbolisch (X-Punkt, instabil)

[1]:

import numpy as np
import matplotlib.pyplot as plt

from pyna.toroidal.equilibrium.stellarator import StellaratorSimple
from pyna.topo.toroidal_cycle import (
    poincare_map_n, poincare_map_n_trajectory,
    jacobian_of_poincare_map, find_cycle, find_all_cycles_near_resonance,
    ToroidalPeriodicOrbitTrace as PeriodicOrbit,
)
from pyna.topo.monodromy import (
    evolve_DPm_along_cycle, build_A_matrix_func, build_delta_b_pol_func,
    orbit_shift_under_perturbation, monodromy_change_under_perturbation,
    CycleVariationalData,
)

plt.rcParams['figure.figsize'] = (10, 7)
plt.rcParams['font.size'] = 12
print('Importe OK')

Importe OK

1. Einrichtung: StellaratorSimple mit q=2/1-Inselkette aufbauen#

Wir verwenden m_h=2, n_h=1, epsilon_h=0.02, um eine q=2-Resonanz zu erzeugen. In diesem Modell gilt q = m_h/n_h = 2, und die Inselkette hat Periode 2 (2 toroidale Umläufe).

Das q-Profil: q(r) = q0 + (q1-q0)*psi = q0 + (q1-q0)*(r/r0)²

  • Bei r=0: q = q0 = 1.5

  • Bei r=r0: q = q1 = 3.5

  • q=2-Fläche bei: psi = (q-q0)/(q1-q0) = 0.5/2.0 = 0.25, also r/r0 = 0.5

[2]:

# Build the stellarator
epsilon_h = 0.02  # helical ripple amplitude
stellarator = StellaratorSimple(
    R0=3.0, r0=0.35, B0=1.0,
    q0=1.5, q1=3.5,
    m_h=2, n_h=1, epsilon_h=epsilon_h,
)

R0, r0 = stellarator.R0, stellarator.r0
field_func = stellarator.field_func

# Domain limits
RZlimit = (R0 - r0*1.5, R0 + r0*1.5, -r0*1.5, r0*1.5)

# Find q=2/1 resonante Fläche
psi_list = stellarator.resonant_psi(2, 1)
psi_res = psi_list[0]
r_res = np.sqrt(psi_res) * r0
print(f"q=2/1 resonance at psi={psi_res:.3f}, r_res={r_res:.4f} m")

# q profile plot
psi_arr = np.linspace(0, 1, 100)
q_arr = [stellarator.q_of_psi(p) for p in psi_arr]

fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.plot(psi_arr, q_arr, 'b-', linewidth=2)
ax.axhline(2.0, color='r', linestyle='--', label='q=2-Resonanz')
ax.axvline(psi_res, color='r', linestyle=':', alpha=0.5)
ax.set_xlabel('Normalisierter Fluss ψ')
ax.set_ylabel('Sicherheitsfaktor q')
ax.set_title('q-Profil von StellaratorSimple')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('q_profile.png', dpi=100, bbox_inches='tight')
plt.show()
print('q-Profil-Plot gespeichert')

q=2/1 resonance at psi=0.250, r_res=0.1750 m
q-Profil-Plot gespeichert

2. Poincaré-Schnitt und Visualisierung der Inselkette#

[3]:

# Compute Poincaré section
print('Berechne Poincaré-Schnitt (das kann ~30 s dauern)...')

# Trace Feldlinien from various starting Punkte
n_field_lines = 8
n_poincare_turns = 100

poincare_Punkte = []  # (R, Z) lists per field line

for i, r_frac in enumerate(np.linspace(0.15, 0.95, n_field_lines)):
    r_start = r_frac * r0
    R_start = R0 + r_start
    Z_start = 0.0

    Rs, Zs = [R_start], [Z_start]
    R_curr, Z_curr = R_start, Z_start

    for _ in range(n_poincare_turns):
        R_next, Z_next = poincare_map_n(
            field_func, [R_curr, Z_curr, 0.0], n_turns=1, dt=0.15, RZlimit=RZlimit
        )
        if np.isnan(R_next):
            break
        Rs.append(R_next)
        Zs.append(Z_next)
        R_curr, Z_curr = R_next, Z_next

    poincare_Punkte.append((np.array(Rs), np.array(Zs)))
    print(f'  Feldlinie {i+1}/{n_field_lines}: {len(Rs)} Punkte')

print('Poincaré-Schnitt berechnet')

Berechne Poincaré-Schnitt (das kann ~30 s dauern)...
  Feldlinie 1/8: 101 Punkte
  Feldlinie 2/8: 101 Punkte
  Feldlinie 3/8: 101 Punkte
  Feldlinie 4/8: 101 Punkte
  Feldlinie 5/8: 101 Punkte
  Feldlinie 6/8: 101 Punkte
  Feldlinie 7/8: 101 Punkte
  Feldlinie 8/8: 101 Punkte
Poincaré-Schnitt berechnet

[4]:

# Find O-Punkte und X-Punkte using Newton-Raphson
print('Suche X/O-Punkte...')

# Known approximate locations from scan (for epsilon_h=0.02, q=2/1)
# X-Punkte: near theta=±0.79 and ±2.36 at r=r_res
# O-Punkts: near theta=±2.32 at r=r_res (for small epsilon)
cycle_seeds = [
    # X-Punkt candidates
    np.array([R0 + r_res*np.cos(0.79), r_res*np.sin(0.79), 0.0]),
    np.array([R0 + r_res*np.cos(-0.79), r_res*np.sin(-0.79), 0.0]),
    np.array([R0 + r_res*np.cos(2.36), r_res*np.sin(2.36), 0.0]),
    np.array([R0 + r_res*np.cos(-2.36), r_res*np.sin(-2.36), 0.0]),
    # O-Punkt candidates (scan found distance minimum near theta=-2.02, r=0.1625)
    np.array([R0 + r_res*np.cos(3.02), r_res*np.sin(3.02), 0.0]),
    np.array([R0 + r_res*np.cos(-2.32), r_res*np.sin(-2.32), 0.0]),
    np.array([R0 + r_res*np.cos(2.32), r_res*np.sin(2.32), 0.0]),
    np.array([R0 + r_res*np.cos(0.25), r_res*np.sin(0.25), 0.0]),
]

found_orbits = []
for seed in cycle_seeds:
    orbit = find_cycle(
        field_func, seed, n_turns=2, dt=0.15, RZlimit=RZlimit,
        max_iter=30, tol=1e-8, n_fallback_seeds=4, fallback_radius=0.02,
    )
    if orbit is None:
        continue
    dist_axis = np.sqrt((orbit.rzphi0[0]-R0)**2 + orbit.rzphi0[1]**2)
    if dist_axis < 0.05:
        continue  # skip axis fixed point
    # Deduplicate
    dup = any(np.linalg.norm(orbit.rzphi0[:2]-fo.rzphi0[:2]) < 1e-4 for fo in found_orbits)
    if not dup:
        found_orbits.append(orbit)

o_Punkte = [o for o in found_orbits if o.is_stable]
x_Punkte = [o for o in found_orbits if not o.is_stable]
print(f'Gefunden: {len(o_Punkte)} O-Punkte und {len(x_Punkte)} X-Punkte')
for o in o_Punkte:
    print(f'  O-Punkt: ({o.rzphi0[0]:.4f}, {o.rzphi0[1]:.4f}), k={o.stability_index:.4f}')
for x in x_Punkte:
    print(f'  X-Punkt: ({x.rzphi0[0]:.4f}, {x.rzphi0[1]:.4f}), k={x.stability_index:.4f}')

Suche X/O-Punkte...
Gefunden: 2 O-Punkte und 4 X-Punkte
  O-Punkt: (3.1237, -0.1237), k=-0.1045
  O-Punkt: (2.8763, -0.1237), k=0.3026
  X-Punkt: (3.1237, 0.1237), k=0.3024
  X-Punkt: (2.8763, 0.1237), k=-0.1042
  X-Punkt: (2.8383, 0.0169), k=2.0324
  X-Punkt: (3.1587, 0.0546), k=2.0320

[5]:

# Plot Poincaré section with X/O Punkte marked
fig, ax = plt.subplots(figsize=(10, 8))

colors = plt.cm.tab10(np.linspace(0, 0.9, n_field_lines))
for i, (Rs, Zs) in enumerate(poincare_Punkte):
    ax.scatter(Rs, Zs, s=0.5, color=colors[i], alpha=0.6)

# Mark O-Punkts (blue circles)
for o in o_Punkte:
    ax.plot(o.rzphi0[0], o.rzphi0[1], 'bo', markersize=10,
            label='O-Punkt' if o is o_Punkte[0] else '', zorder=5)

# Mark X-Punkte (red x's)
for x in x_Punkte:
    ax.plot(x.rzphi0[0], x.rzphi0[1], 'rx', markersize=12, markeredgewidth=2,
            label='X-Punkt' if x is x_Punkte[0] else '', zorder=5)

# Mark resonante Fläche
theta_arr = np.linspace(0, 2*np.pi, 200)
ax.plot(R0 + r_res*np.cos(theta_arr), r_res*np.sin(theta_arr), 'k--',
        linewidth=1, alpha=0.5, label=f'q=2 surface (r={r_res:.3f})')

ax.set_xlabel('R (m)')
ax.set_ylabel('Z (m)')
ax.set_title(f'Poincaré-Schnitt - q=2/1-Inselkette\n(epsilon_h={epsilon_h})')
ax.set_aspect('equal')
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('poincare_section.png', dpi=100, bbox_inches='tight')
plt.show()
print('Poincaré-Schnitt gespeichert')

# Note: For interactive Plotly version, use:
# import plotly.graph_objects as go
# fig = go.Figure()
# for Rs, Zs in poincare_Punkte:
#     fig.add_trace(go.Scatter(x=Rs, y=Zs, Modus='markers', marker=dict(size=1)))
# for o in o_Punkte:
#     fig.add_trace(go.Scatter(x=[o.rzphi0[0]], y=[o.rzphi0[1]], Modus='markers',
#                              marker=dict(size=10, symbol='circle', color='blue')))
# fig.show()

Poincaré-Schnitt gespeichert

3. Monodromieanalyse für den O-Punkt#

Die Monodromiematrix \(M = DP(\phi_0 + 2\pi n)\) kodiert die Stabilität des periodischen Orbits. Für einen O-Punkt (elliptisch) liegen die Eigenwerte auf dem Einheitskreis: \(|\lambda| = 1\).

Greene-Residuum: \(R_G = (2 - \text{Tr}(M))/4\)

  • \(0 < R_G < 1\): elliptisch (stabil)

  • \(R_G < 0\): hyperbolisch (instabil, Standardfall)

  • \(R_G > 1\): hyperbolisch mit Reflexion

[6]:

# Pick the first O-Punkt for analysis
if not o_Punkte:
    print('No O-Punkts found; using first found orbit for demo')
    opoint = found_orbits[0] if found_orbits else None
else:
    opoint = o_Punkte[0]

if opoint is None:
    print('No orbits found to analyze')
else:
    print(f'Analysiere O-Punkt bei ({opoint.rzphi0[0]:.4f}, {opoint.rzphi0[1]:.4f})')
    print(f'  stability_index = {opoint.stability_index:.6f}')
    print(f'  eigenvalues = {opoint.eigenvalues}')

    # Compute full monodromy analysis
    print('Berechne Monodromie (Variationsgleichungen)...')
    monodromy_O = evolve_DPm_along_cycle(
        field_func, opoint, dt_output=0.1, rtol=1e-8, atol=1e-9
    )

    print(f'  Monodromiematrix M:')
    print(f'    {monodromy_O.DPm}')
    print(f'  det(M) = {np.linalg.det(monodromy_O.DPm):.8f} (sollte sein ~1)')
    print(f'  Tr(M)/2 = {monodromy_O.stability_index:.6f}')
    print(f"  Greene's residue = {monodromy_O.Greene_residue:.6f}")
    print(f'  Eigenwerte = {monodromy_O.eigenvalues}')

Analysiere O-Punkt bei (3.1237, -0.1237)
  stability_index = -0.104472
  eigenvalues = [-0.10447198+0.99401547j -0.10447198-0.99401547j]
Berechne Monodromie (Variationsgleichungen)...
  Monodromiematrix M:
    [[-1.95431405  1.21064084]
 [-3.64518107  1.74639533]]
  det(M) = 1.00000011 (sollte sein ~1)
  Tr(M)/2 = -0.103959
  Greene's residue = 0.551980
  Eigenwerte = [-0.10395936+0.9945816j -0.10395936-0.9945816j]

[7]:

if opoint is not None and 'monodromy_O' in dir():
    phi_arr = monodromy_O.phi_arr

    # Compute eigenvalues and det along orbit
    DP_eigvals = np.array([np.linalg.eigvals(DP) for DP in monodromy_O.DPm_arr])
    DP_dets = np.array([np.linalg.det(DP) for DP in monodromy_O.DPm_arr])

    fig, axes = plt.subplots(3, 1, figsize=(10, 10))

    # Plot |eigenvalues| of DPm along orbit
    ax = axes[0]
    ax.plot(phi_arr/(2*np.pi), np.abs(DP_eigvals[:, 0]), 'b-', label='|lambda(phi)|')
    ax.plot(phi_arr/(2*np.pi), np.abs(DP_eigvals[:, 1]), 'r--', label='|lambda(phi)|')
    ax.axhline(1.0, color='k', linestyle=':', alpha=0.5)
    ax.set_xlabel('φ/(2π)')
    ax.set_ylabel('|eigenvalue of DP|')
    ax.set_title('O-Punkt: Eigenwertentwicklung entlang des Orbits')
    ax.legend()
    ax.grid(True, alpha=0.3)

    # Plot det(DPm) - should stay near 1
    ax = axes[1]
    ax.plot(phi_arr/(2*np.pi), DP_dets, 'g-', linewidth=2)
    ax.axhline(1.0, color='k', linestyle=':', alpha=0.5, label='det=1 (area-preserving)')
    ax.set_xlabel('φ/(2π)')
    ax.set_ylabel('det(DPm(φ))')
    ax.set_title('Flächenerhaltungsprüfung: det(DPm) ~ 1')
    ax.legend()
    ax.grid(True, alpha=0.3)

    # Plot Greene residue evolution (Tr(J)/2)
    ax = axes[2]
    Tr_half = np.array([np.trace(J)/2 for J in monodromy_O.DPm_arr])
    ax.plot(phi_arr/(2*np.pi), Tr_half, 'm-', linewidth=2)
    ax.axhline(0.0, color='k', linestyle=':', alpha=0.3)
    ax.axhline(monodromy_O.stability_index, color='r', linestyle='--',
               label=f'Final Tr(M)/2 = {monodromy_O.stability_index:.4f}')
    ax.set_xlabel('φ/(2π)')
    ax.set_ylabel('Tr(DPm(φ))/2')
    ax.set_title('Entwicklung des Stabilitätsindex entlang des O-Punkt-Orbits')
    ax.legend()
    ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('monodromy_opoint.png', dpi=100, bbox_inches='tight')
    plt.show()
    print('Monodromie-O-Punkt-Plot gespeichert')

Monodromie-O-Punkt-Plot gespeichert

4. Monodromieanalyse für den X-Punkt#

Für einen X-Punkt (hyperbolisch) sind die Eigenwerte reell und \(|\lambda| > 1\). Nahe Orbits laufen exponentiell auseinander; der zugehörige Ljapunow-Exponent wird aus \(\lambda_{max}\) bestimmt.

[8]:

# Pick the first X-Punkt
if not x_Punkte:
    print('No X-Punkte found')
    xpoint = None
else:
    xpoint = x_Punkte[0]
    print(f'Analysiere X-Punkt bei ({xpoint.rzphi0[0]:.4f}, {xpoint.rzphi0[1]:.4f})')
    print(f'  stability_index = {xpoint.stability_index:.6f}')
    print(f'  eigenvalues = {xpoint.eigenvalues}')

    monodromy_X = evolve_DPm_along_cycle(
        field_func, xpoint, dt_output=0.1, rtol=1e-8, atol=1e-9
    )

    print(f'  det(M) = {np.linalg.det(monodromy_X.DPm):.8f}')
    print(f"  Greene's residue = {monodromy_X.Greene_residue:.6f} (< 0 -> hyperbolisch)")

Analysiere X-Punkt bei (3.1237, 0.1237)
  stability_index = 0.302381
  eigenvalues = [0.30238099+0.95326451j 0.30238099-0.95326451j]
  det(M) = 0.99999969
  Greene's residue = 0.348880 (< 0 -> hyperbolisch)

[9]:

if xpoint is not None and 'monodromy_X' in dir():
    phi_arr_X = monodromy_X.phi_arr
    DP_eigvals_X = np.array([np.linalg.eigvals(DP) for DP in monodromy_X.DPm_arr])
    DP_dets_X = np.array([np.linalg.det(DP) for DP in monodromy_X.DPm_arr])

    fig, axes = plt.subplots(2, 1, figsize=(10, 7))

    ax = axes[0]
    ax.semilogy(phi_arr_X/(2*np.pi), np.abs(DP_eigvals_X[:, 0]), 'b-', label='|lambda(phi)|')
    ax.semilogy(phi_arr_X/(2*np.pi), np.abs(DP_eigvals_X[:, 1]), 'r--', label='|lambda(phi)|')
    ax.axhline(1.0, color='k', linestyle=':', alpha=0.5)
    ax.set_xlabel('φ/(2π)')
    ax.set_ylabel('|eigenvalue of DP| (log scale)')
    ax.set_title('X-Punkt: hyperbolisches Eigenwertwachstum entlang des Orbits')
    ax.legend()
    ax.grid(True, alpha=0.3)

    ax = axes[1]
    ax.plot(phi_arr_X/(2*np.pi), DP_dets_X, 'g-', linewidth=2)
    ax.axhline(1.0, color='k', linestyle=':', alpha=0.5, label='det=1')
    ax.set_xlabel('φ/(2π)')
    ax.set_ylabel('det(DPm(φ))')
    ax.set_title('Flächenerhaltungsprüfung für X-Punkt-Orbit')
    ax.legend()
    ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('monodromy_xpoint.png', dpi=100, bbox_inches='tight')
    plt.show()

5b. Stabile und instabile Mannigfaltigkeiten vom X-Punkt#

Die stabile Mannigfaltigkeit \(W^s\) (blau, kühl) enthält Trajektorien, die bei Vorwärtsiteration zum X-Punkt konvergieren. Die instabile Mannigfaltigkeit \(W^u\) (orange, warm) enthält Trajektorien, die vom X-Punkt divergieren. Zusammen bilden sie das Separatrix-Gerüst der Inselkette.

[10]:

# === Stable/Unstable Manifold Visualization ===
from pyna.topo.variational import PoincaréMapVariationalEquations
from pyna.topo.manifold_improve import StableManifold, UnstableManifold
from pyna.toroidal.visual.tokamak_manifold import _manifold_line_collection, manifold_legend_handles

# field_func_2d wrapper: stellarator.field_func([R,Z,phi]) returns unit tangent
# We need (R,Z,phi) -> (dR/dphi, dZ/dphi)
def field_func_2d(R, Z, phi):
    tang = stellarator.field_func(np.array([R, Z, phi]))  # [dR/ds, dZ/ds, dphi/ds]
    dphi_ds = tang[2]
    if abs(dphi_ds) < 1e-15:
        return np.array([0.0, 0.0])
    return np.array([tang[0] / dphi_ds, tang[1] / dphi_ds])

if xpoint is not None and 'monodromy_X' in dir():
    R_xpt, Z_xpt = xpoint.rzphi0[0], xpoint.rzphi0[1]
    print(f'Lasse Mannigfaltigkeiten vom X-Punkt wachsen ({R_xpt:.4f}, {Z_xpt:.4f})')

    # Use the Jacobian from monodromy_X (2x2 monodromy matrix)
    M = monodromy_X.DPm
    eigvals = np.linalg.eigvals(M)
    print(f'Monodromie-Eigenwerte: {eigvals}')
    print(f'|lambda_1|={abs(eigvals[0]):.4f}, |lambda_2|={abs(eigvals[1]):.4f}  (product~1: {abs(np.prod(eigvals)):.4f})')

    # Grow manifolds
    sm = StableManifold([R_xpt, Z_xpt], M, field_func_2d)
    um = UnstableManifold([R_xpt, Z_xpt], M, field_func_2d)
    sm.grow(n_turns=12, init_length=5e-5, n_init_pts=6, both_sides=True)
    um.grow(n_turns=12, init_length=5e-5, n_init_pts=6, both_sides=True)
    print(f'Stabile Segmente: {len(sm.segments)}, Instabile Segmente: {len(um.segments)}')

    # Combined figure: Poincaré + manifolds
    fig_mf, ax_mf = plt.subplots(figsize=(9, 8))
    ax_mf.set_facecolor('white')

    # Replot Poincaré scatter (from existing results_n variable)
    try:
        R_pts = results_n[:, 0]; Z_pts = results_n[:, 1]
        psi_n = np.clip(((R_pts - R0)**2 + Z_pts**2) / r_res**2, 0, 1.4)
        from matplotlib.colors import Normalize
        import matplotlib.cm as cm
        colors_sc = cm.plasma(np.clip(psi_n / 1.4, 0.05, 0.95))
        ax_mf.scatter(R_pts, Z_pts, s=0.5, c=colors_sc, rasterized=True, alpha=0.5, zorder=2)
    except NameError:
        pass  # Poincaré results not available, skip

    # Manifolds with arc-length coloring
    for seg in sm.segments:
        if len(seg) > 2:
            lc, _ = _manifold_line_collection(seg, cmap='GnBu')
            lc.set_linewidth(1.3); lc.set_alpha(0.92); lc.set_zorder(6)
            ax_mf.add_collection(lc)

    for seg in um.segments:
        if len(seg) > 2:
            lc, _ = _manifold_line_collection(seg, cmap='Oranges')
            lc.set_linewidth(1.3); lc.set_alpha(0.92); lc.set_zorder(6)
            ax_mf.add_collection(lc)

    # X-Punkt marker
    ax_mf.plot(R_xpt, Z_xpt, 'r+', ms=14, mew=2.5, zorder=8, label='X-Punkt')

    # Also plot all found O/X Punkte
    for op in o_Punkte:
        ax_mf.plot(op.rzphi0[0], op.rzphi0[1], 'go', ms=7, mew=1.5, zorder=7)
    for xp in x_Punkte:
        ax_mf.plot(xp.rzphi0[0], xp.rzphi0[1], 'r+', ms=12, mew=2.5, zorder=8)

    # Resonant surface circle
    theta_c = np.linspace(0, 2*np.pi, 300)
    ax_mf.plot(R0 + r_res*np.cos(theta_c), r_res*np.sin(theta_c),
               '--', color='tomato', lw=0.8, alpha=0.6, label='$q=2/1$ surface')
    ax_mf.plot(R0 + stellarator.r0*np.cos(theta_c), stellarator.r0*np.sin(theta_c),
               'k-', lw=1.2, label='LCFS')

    # Legend + labels
    mfld_handles = manifold_legend_handles('Oranges', 'GnBu')
    ax_mf.legend(handles=mfld_handles + [
        plt.Line2D([0],[0], marker='+', color='r', ms=10, mew=2, lw=0, label='X-Punkt'),
        plt.Line2D([0],[0], marker='o', color='g', ms=7, lw=0, label='O-Punkt'),
        plt.Line2D([0],[0], color='k', lw=1.2, label='LCFS'),
    ], loc='upper right', fontsize=9, framealpha=0.9)

    ax_mf.set_xlim(R0 - 1.3*stellarator.r0, R0 + 1.3*stellarator.r0)
    ax_mf.set_ylim(-1.3*stellarator.r0, 1.3*stellarator.r0)
    ax_mf.set_xlabel('$R$ (m)', fontsize=12)
    ax_mf.set_ylabel('$Z$ (m)', fontsize=12)
    ax_mf.set_title('Stabile ($W^s$, blau) & instabile ($W^u$, orange) Mannigfaltigkeiten - $q=2/1$-X-Punkt', fontsize=12)
    ax_mf.set_aspect('equal')
    plt.tight_layout()
    plt.savefig('island_manifolds.png', dpi=150, bbox_inches='tight')
    plt.show()
    print('Mannigfaltigkeitsabbildung fertig.')
else:
    print('Kein X-Punkt oder keine Monodromie verfügbar; Mannigfaltigkeitsvisualisierung wird übersprungen.')

Lasse Mannigfaltigkeiten vom X-Punkt wachsen (3.1237, 0.1237)
Monodromie-Eigenwerte: [0.3022398+0.95323176j 0.3022398-0.95323176j]
|lambda_1|=1.0000, |lambda_2|=1.0000  (product~1: 1.0000)
Stabile Segmente: 2, Instabile Segmente: 2
Mannigfaltigkeitsabbildung fertig.

5. Wirkung einer Störung: Orbitverschiebung unter \(\delta B\)#

Wir fügen eine kleine zusätzliche helikale Störung \(\delta B\) mit Amplitude epsilon_pert hinzu. Die Orbitverschiebung erfüllt:

\[\frac{dX_{cyc}}{d\phi} = A \cdot X_{cyc} + \delta b_{pol}\]

wobei \(\delta b_{pol} = [R\delta B_R/B_\phi - R B_R \delta B_\phi/B_\phi^2; \text{entsprechend für Z}]\).

[11]:

epsilon_pert = 0.005  # perturbation amplitude

# Build perturbed stellarator (slightly different helical phase to create delta B)
stellarator_pert = StellaratorSimple(
    R0=3.0, r0=0.35, B0=1.0,
    q0=1.5, q1=3.5,
    m_h=2, n_h=1, epsilon_h=epsilon_h + epsilon_pert,
)

def delta_field_func(rzphi):
    """Perturbation field: difference between perturbed and unperturbed."""
    f0 = np.asarray(field_func(rzphi))
    f1 = np.asarray(stellarator_pert.field_func(rzphi))
    # Note: field_func returns unit tangent vectors, not raw B fields.
    # For the orbit displacement, we use the difference in the phi-parameterized field.
    # The perturbation in terms of dR/dphi, dZ/dphi:
    g0 = np.array([f0[0]/f0[2], f0[1]/f0[2], 1.0])
    g1 = np.array([f1[0]/f1[2], f1[1]/f1[2], 1.0])
    return g1 - g0

if opoint is not None and 'monodromy_O' in dir():
    print(f'Berechne Orbitverschiebung für O-Punkt unter epsilon_pert={epsilon_pert}...')
    orbit_shift_O = orbit_shift_under_perturbation(
        field_func, delta_field_func, opoint, monodromy_O
    )

    print(f'  Anfängliche Orbitverschiebung: dR={orbit_shift_O[0,0]:.6f}, dZ={orbit_shift_O[0,1]:.6f}')
    print(f'  Max. |Verschiebung| entlang des Orbits: {np.max(np.linalg.norm(orbit_shift_O, axis=1)):.6f} m')

if xpoint is not None and 'monodromy_X' in dir():
    print(f'Berechne Orbitverschiebung für X-Punkt unter epsilon_pert={epsilon_pert}...')
    orbit_shift_X = orbit_shift_under_perturbation(
        field_func, delta_field_func, xpoint, monodromy_X
    )
    print(f'  Anfängliche Orbitverschiebung: dR={orbit_shift_X[0,0]:.6f}, dZ={orbit_shift_X[0,1]:.6f}')

Berechne Orbitverschiebung für O-Punkt unter epsilon_pert=0.005...
  Anfängliche Orbitverschiebung: dR=-0.000000, dZ=-0.000000
  Max. |Verschiebung| entlang des Orbits: 0.000000 m
Berechne Orbitverschiebung für X-Punkt unter epsilon_pert=0.005...
  Anfängliche Orbitverschiebung: dR=0.000000, dZ=-0.000000

[12]:

if opoint is not None and 'orbit_shift_O' in dir():
    phi_arr_O = monodromy_O.phi_arr

    fig, axes = plt.subplots(2, 1, figsize=(10, 7))

    ax = axes[0]
    ax.plot(phi_arr_O/(2*np.pi), orbit_shift_O[:, 0], 'b-', label='δR(φ)')
    ax.plot(phi_arr_O/(2*np.pi), orbit_shift_O[:, 1], 'r-', label='δZ(φ)')
    ax.set_xlabel('φ/(2π)')
    ax.set_ylabel('Orbitverschiebung (m)')
    ax.set_title(f'O-Punkt-Orbitverschiebung unter δε={epsilon_pert}')
    ax.legend()
    ax.grid(True, alpha=0.3)

    if xpoint is not None and 'orbit_shift_X' in dir():
        phi_arr_X = monodromy_X.phi_arr
        ax = axes[1]
        ax.plot(phi_arr_X/(2*np.pi), orbit_shift_X[:, 0], 'b-', label='δR(φ)')
        ax.plot(phi_arr_X/(2*np.pi), orbit_shift_X[:, 1], 'r-', label='δZ(φ)')
        ax.set_xlabel('φ/(2π)')
        ax.set_ylabel('Orbitverschiebung (m)')
        ax.set_title(f'X-Punkt-Orbitverschiebung unter δε={epsilon_pert}')
        ax.legend()
        ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('orbit_shift.png', dpi=100, bbox_inches='tight')
    plt.show()

6. Stabilitätsdiagramm: Greene-Residuum gegen Störungsamplitude#

Wir scannen epsilon_pert von 0 bis 0.02. Für jede Amplitude:

  • O- und X-Zyklen finden

  • Greene-Residuum \(R_G = (2 - \text{Tr}(M))/4\) berechnen

Der Übergang \(R_G: 0 \to 1\) markiert die Inselzerstörung (KAM-ähnlicher Übergang).

[13]:

print('Berechne Stabilitätsdiagramm (Scan über epsilon_h)...')
epsilon_arr = np.linspace(0.005, 0.04, 12)

greene_O = []  # O-Punkt Greene residues
greene_X = []  # X-Punkt Greene residues
valid_eps = []

for eps in epsilon_arr:
    st_eps = StellaratorSimple(
        R0=3.0, r0=0.35, B0=1.0,
        q0=1.5, q1=3.5,
        m_h=2, n_h=1, epsilon_h=eps,
    )
    ff = st_eps.field_func

    # Find X-Punkt near known seed (2.8763, -0.1237)
    orbit_x = find_cycle(
        ff, np.array([2.8763, -0.1237, 0.0]),
        n_turns=2, dt=0.15, RZlimit=RZlimit,
        max_iter=30, tol=1e-8, n_fallback_seeds=6, fallback_radius=0.02,
    )

    # Find O-Punkt near known seed
    psi_r = st_eps.resonant_psi(2, 1)[0]
    r_r = np.sqrt(psi_r) * st_eps.r0
    opoint_seed = np.array([st_eps.R0 + r_r*np.cos(-2.32), r_r*np.sin(-2.32), 0.0])
    orbit_o = find_cycle(
        ff, opoint_seed,
        n_turns=2, dt=0.15, RZlimit=RZlimit,
        max_iter=30, tol=1e-8, n_fallback_seeds=8, fallback_radius=0.02,
    )

    rg_x = None
    rg_o = None

    if orbit_x is not None:
        M_x = orbit_x.DPm
        rg_x = (2.0 - np.trace(M_x)) / 4.0

    if orbit_o is not None:
        dist_axis = np.sqrt((orbit_o.rzphi0[0]-st_eps.R0)**2 + orbit_o.rzphi0[1]**2)
        if dist_axis > 0.05:
            M_o = orbit_o.DPm
            rg_o = (2.0 - np.trace(M_o)) / 4.0

    if rg_x is not None or rg_o is not None:
        valid_eps.append(eps)
        greene_X.append(rg_x)
        greene_O.append(rg_o)
        rg_x_str = f'{rg_x:.4f}' if rg_x is not None else 'N/A'
        rg_o_str = f'{rg_o:.4f}' if rg_o is not None else 'N/A'
        print(f'  eps={eps:.3f}: R_G(X)={rg_x_str}, R_G(O)={rg_o_str}')

print(f'Berechnet {len(valid_eps)} Datenpunkte')

Berechne Stabilitätsdiagramm (Scan über epsilon_h)...
  eps=0.005: R_G(X)=0.0024, R_G(O)=0.0024
  eps=0.008: R_G(X)=0.0325, R_G(O)=0.0325
  eps=0.011: R_G(X)=0.0852, R_G(O)=0.0852
  eps=0.015: R_G(X)=0.1612, R_G(O)=0.1612
  eps=0.018: R_G(X)=0.2615, R_G(O)=0.2615
  eps=0.021: R_G(X)=0.3873, R_G(O)=0.3873
  eps=0.024: R_G(X)=0.5402, R_G(O)=0.5402
  eps=0.027: R_G(X)=0.7220, R_G(O)=0.7220
  eps=0.030: R_G(X)=0.9348, R_G(O)=0.9348
  eps=0.034: R_G(X)=1.1812, R_G(O)=1.1812
  eps=0.037: R_G(X)=1.4640, R_G(O)=1.4640
  eps=0.040: R_G(X)=1.7865, R_G(O)=1.7865
Berechnet 12 Datenpunkte

[14]:

fig, ax = plt.subplots(figsize=(9, 5))

# Filter valid entries
eps_X = [e for e, g in zip(valid_eps, greene_X) if g is not None]
rg_X = [g for g in greene_X if g is not None]
eps_O = [e for e, g in zip(valid_eps, greene_O) if g is not None]
rg_O = [g for g in greene_O if g is not None]

if eps_X:
    ax.plot(eps_X, rg_X, 'rx-', linewidth=2, markersize=8, label="Greene-Residuum des X-Punkts")
if eps_O:
    ax.plot(eps_O, rg_O, 'bo-', linewidth=2, markersize=8, label="Greene-Residuum des O-Punkts")

ax.axhline(0.0, color='k', linestyle='--', alpha=0.5, label='R_G=0 (elliptische Grenze)')
ax.axhline(1.0, color='k', linestyle=':', alpha=0.5, label='R_G=1 (hyperbolische Grenze)')
ax.fill_between([0, 0.05], 0, 1, alpha=0.1, color='green', label='Elliptischer Bereich')

ax.set_xlabel('Amplitude des helikalen Ripples ε_h')
ax.set_ylabel("Greene-Residuum R_G")
ax.set_title("Stabilitätsdiagramm: Greene-Residuum gegen Störungsamplitude\n"
             f"(q=2/1 island chain, m_h=2, n_h=1)")
ax.legend(loc='best')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 0.05)

plt.tight_layout()
plt.savefig('stability_diagram.png', dpi=100, bbox_inches='tight')
plt.show()
print('Stabilitätsdiagramm gespeichert')

Stabilitätsdiagramm gespeichert

Zusammenfassung#

Dieses Notebook hat gezeigt:

  1. Suche periodischer Orbits mit Newton-Raphson auf der Poincaré-Karte

  2. Jacobian-Entwicklung entlang des Orbits über Variationsgleichungen

  3. Monodromiematrix DPm der \(m\)-Umlauf-Poincaré-Karte mit det(M) ~ 1 (Flächenerhaltung)

  4. Greene-Residuum \(R_G\) als Stabilitätsindikator

  5. Orbitverschiebung unter Störung mit periodischer Randbedingung

  6. Stabilitätsdiagramm: wie sich die Stabilität einer Inselkette mit der Störungsamplitude ändert

Wichtige Ergebnisse:

  • O-Punkte haben \(0 < R_G < 1\) (elliptisch, stabil)

  • X-Punkte haben \(R_G < 0\) oder \(R_G > 1\) (hyperbolisch, instabil)

  • det(DPm(phi)) ~ 1 durchgehend; das bestätigt Flächenerhaltung (Hamiltonsche Struktur)

Diese Portierung einer Jacobian-/Monodromieanalyse aus einem Julia-Workflow nach Python (analytischer Stellarator) demonstriert dieselben physikalischen Algorithmen ohne proprietäre Gleichgewichtsdaten.