Analyse du Jacobien des îlots : orbites périodiques et matrices de monodromie#
Ce notebook démontre l’analyse complète du Jacobien et de la monodromie des chaînes d’îlots magnétiques dans un équilibre de stellarator analytique, portée depuis un workflow Julia de type MCF_scripts.
Contexte scientifique#
Une orbite périodique de ligne de champ ou cycle est une ligne de champ qui revient à son point de départ après exactement \(n\) tours toroïdaux : \(X(\phi + 2\pi n) = X(\phi)\).
Ce sont des points fixes de la \(n\)-ième itération de la carte de Poincaré \(P^n\). Pour une résonance \(q = m/n\) (avec \(q = B_\phi r / (B_{pol} R)\)), la période orbitale est de \(m\) tours.
La matrice jacobienne \(DX(\phi)\) évolue selon :
\[\frac{dDX}{d\phi} = A(r(\phi), z(\phi), \phi) \cdot DX, \quad DX(\phi_0) = I\]
où \(A_{ij} = \partial(R B_{pol,i}/B_\phi) / \partial x_j\).
La matrice de monodromie \(M = DP(\phi_0 + 2\pi n)\) donne la carte linéarisée à \(n\) tours. Valeurs propres de \(M\) :
\(|\lambda| = 1\) : elliptique (point O, stable)
\(|\lambda| > 1\) : hyperbolique (point X, instable)
[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('Imports OK')
Imports OK
1. Mise en place : construire StellaratorSimple avec une chaîne d’îlots q=2/1#
Nous utilisons m_h=2, n_h=1, epsilon_h=0.02 pour créer une résonance q=2. Dans ce modèle, q = m_h/n_h = 2, et la chaîne d’îlots a une période 2 (2 tours toroïdaux).
Le profil q : q(r) = q0 + (q1-q0)*psi = q0 + (q1-q0)*(r/r0)²
À r=0 : q = q0 = 1.5
À r=r0 : q = q1 = 3.5
Surface q=2 à : psi = (q-q0)/(q1-q0) = 0.5/2.0 = 0.25, donc 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 surface résonante
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='résonance q=2')
ax.axvline(psi_res, color='r', linestyle=':', alpha=0.5)
ax.set_xlabel('Flux normalisé ψ')
ax.set_ylabel('Facteur de sécurité q')
ax.set_title('profil q de 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('tracé du profil q enregistré')
q=2/1 resonance at psi=0.250, r_res=0.1750 m
tracé du profil q enregistré
2. Section de Poincaré et visualisation de la chaîne d’îlots#
[3]:
# Compute Poincaré section
print('Calcul de la section de Poincaré (peut prendre ~30 s)...')
# Trace lignes de champ from various starting points
n_field_lines = 8
n_poincare_turns = 100
poincare_points = [] # (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_points.append((np.array(Rs), np.array(Zs)))
print(f' Ligne de champ {i+1}/{n_field_lines}: {len(Rs)} points')
print('Section de Poincaré calculée')
Calcul de la section de Poincaré (peut prendre ~30 s)...
Ligne de champ 1/8: 101 points
Ligne de champ 2/8: 101 points
Ligne de champ 3/8: 101 points
Ligne de champ 4/8: 101 points
Ligne de champ 5/8: 101 points
Ligne de champ 6/8: 101 points
Ligne de champ 7/8: 101 points
Ligne de champ 8/8: 101 points
Section de Poincaré calculée
[4]:
# Find points O et points X using Newton-Raphson
print('Recherche des points X/O...')
# Known approximate locations from scan (for epsilon_h=0.02, q=2/1)
# points X: near theta=±0.79 and ±2.36 at r=r_res
# point Os: near theta=±2.32 at r=r_res (for small epsilon)
cycle_seeds = [
# point X 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]),
# point O 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_points = [o for o in found_orbits if o.is_stable]
x_points = [o for o in found_orbits if not o.is_stable]
print(f'Trouvé : {len(o_points)} points O et {len(x_points)} points X')
for o in o_points:
print(f' point O: ({o.rzphi0[0]:.4f}, {o.rzphi0[1]:.4f}), k={o.stability_index:.4f}')
for x in x_points:
print(f' Point X : ({x.rzphi0[0]:.4f}, {x.rzphi0[1]:.4f}), k={x.stability_index:.4f}')
Recherche des points X/O...
Trouvé : 2 points O et 4 points X
point O: (3.1237, -0.1237), k=-0.1045
point O: (2.8763, -0.1237), k=0.3026
Point X : (3.1237, 0.1237), k=0.3024
Point X : (2.8763, 0.1237), k=-0.1042
Point X : (2.8383, 0.0169), k=2.0324
Point X : (3.1587, 0.0546), k=2.0320
[5]:
# Plot Poincaré section with X/O points 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_points):
ax.scatter(Rs, Zs, s=0.5, color=colors[i], alpha=0.6)
# Mark point Os (blue circles)
for o in o_points:
ax.plot(o.rzphi0[0], o.rzphi0[1], 'bo', markersize=10,
label='point O' if o is o_points[0] else '', zorder=5)
# Mark points X (red x’s)
for x in x_points:
ax.plot(x.rzphi0[0], x.rzphi0[1], 'rx', markersize=12, markeredgewidth=2,
label='point X' if x is x_points[0] else '', zorder=5)
# Mark surface résonante
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'Section de Poincaré - chaîne d’îlots q=2/1\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('Section de Poincaré enregistrée')
# Note: For interactive Plotly version, use:
# import plotly.graph_objects as go
# fig = go.Figure()
# for Rs, Zs in poincare_points:
# fig.add_trace(go.Scatter(x=Rs, y=Zs, mode='markers', marker=dict(size=1)))
# for o in o_points:
# fig.add_trace(go.Scatter(x=[o.rzphi0[0]], y=[o.rzphi0[1]], mode='markers',
# marker=dict(size=10, symbol='circle', color='blue')))
# fig.show()
Section de Poincaré enregistrée
3. Analyse de monodromie pour le point O#
La matrice de monodromie \(M = DP(\phi_0 + 2\pi n)\) encode la stabilité de l’orbite périodique. Pour un point O (elliptique), les valeurs propres sont sur le cercle unité : \(|\lambda| = 1\).
Résidu de Greene : \(R_G = (2 - \text{Tr}(M))/4\)
\(0 < R_G < 1\) : elliptique (stable)
\(R_G < 0\) : hyperbolique (instable, standard)
\(R_G > 1\) : hyperbolique avec réflexion
[6]:
# Pick the first point O for analysis
if not o_points:
print('No point Os found; using first found orbit for demo')
opoint = found_orbits[0] if found_orbits else None
else:
opoint = o_points[0]
if opoint is None:
print('No orbits found to analyze')
else:
print(f'Analyse du point O en ({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('Calcul de la monodromie (équations variationnelles)...')
monodromy_O = evolve_DPm_along_cycle(
field_func, opoint, dt_output=0.1, rtol=1e-8, atol=1e-9
)
print(f' Matrice de monodromie M:')
print(f' {monodromy_O.DPm}')
print(f' det(M) = {np.linalg.det(monodromy_O.DPm):.8f} (devrait être ~1)')
print(f' Tr(M)/2 = {monodromy_O.stability_index:.6f}')
print(f" Greene’s residue = {monodromy_O.Greene_residue:.6f}")
print(f' Valeurs propres = {monodromy_O.eigenvalues}')
Analyse du point O en (3.1237, -0.1237)
stability_index = -0.104472
eigenvalues = [-0.10447198+0.99401547j -0.10447198-0.99401547j]
Calcul de la monodromie (équations variationnelles)...
Matrice de monodromie M:
[[-1.95431405 1.21064084]
[-3.64518107 1.74639533]]
det(M) = 1.00000011 (devrait être ~1)
Tr(M)/2 = -0.103959
Greene’s residue = 0.551980
Valeurs propres = [-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('Point O : évolution des valeurs propres le long de l’orbite')
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('Contrôle de préservation d’aire: 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('Évolution de l’indice de stabilité le long de l’orbite du point O')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('monodromy_opoint.png', dpi=100, bbox_inches='tight')
plt.show()
print('tracé de monodromie du point O enregistré')
tracé de monodromie du point O enregistré
4. Analyse de monodromie pour le point X#
Pour un point X (hyperbolique), les valeurs propres sont réelles et \(|\lambda| > 1\). Les orbites voisines divergent exponentiellement ; le multiplicateur instable est \(\lambda_{max}\) et l’exposant associé est \(\log|\lambda_{max}|\) par période.
[8]:
# Pick the first point X
if not x_points:
print('No points X found')
xpoint = None
else:
xpoint = x_points[0]
print(f'Analyse du point X en ({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 -> hyperbolique)")
Analyse du point X en (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 -> hyperbolique)
[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('Point X : croissance hyperbolique des valeurs propres le long de l’orbite')
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('Contrôle de préservation d’aire pour l’orbite du point X')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('monodromy_xpoint.png', dpi=100, bbox_inches='tight')
plt.show()
5b. Variétés stables et instables depuis le point X#
La variété stable \(W^s\) (bleu, froid) contient les trajectoires qui convergent vers le point X sous l’itération directe. La variété instable \(W^u\) (orange, chaud) contient les trajectoires qui divergent du point X. Ensemble, elles forment le squelette de séparatrice de la chaîne d’îlots.
[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'Croissance des variétés depuis le point X ({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'Valeurs propres de monodromie : {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'Segments stables: {len(sm.segments)}, Segments instables: {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)
# point X marker
ax_mf.plot(R_xpt, Z_xpt, 'r+', ms=14, mew=2.5, zorder=8, label='point X')
# Also plot all found O/X points
for op in o_points:
ax_mf.plot(op.rzphi0[0], op.rzphi0[1], 'go', ms=7, mew=1.5, zorder=7)
for xp in x_points:
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='point X'),
plt.Line2D([0],[0], marker='o', color='g', ms=7, lw=0, label='point O'),
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('Variétés stable ($W^s$, bleu) et instable ($W^u$, orange) - point X $q=2/1$', fontsize=12)
ax_mf.set_aspect('equal')
plt.tight_layout()
plt.savefig('island_manifolds.png', dpi=150, bbox_inches='tight')
plt.show()
print('Figure des variétés terminée.')
else:
print('Aucun point X ou monodromie disponible ; visualisation des variétés ignorée.')
Croissance des variétés depuis le point X (3.1237, 0.1237)
Valeurs propres de monodromie : [0.3022398+0.95323176j 0.3022398-0.95323176j]
|lambda_1|=1.0000, |lambda_2|=1.0000 (product~1: 1.0000)
Segments stables: 2, Segments instables: 2
Figure des variétés terminée.
5. Effet d’une perturbation : déplacement orbital sous \(\delta B\)#
Nous ajoutons une petite perturbation hélicoïdale supplémentaire \(\delta B\) d’amplitude epsilon_pert. Le déplacement de l’orbite vérifie :
\[\frac{dX_{cyc}}{d\phi} = A \cdot X_{cyc} + \delta b_{pol}\]
où \(\delta b_{pol} = [R\delta B_R/B_\phi - R B_R \delta B_\phi/B_\phi^2; \text{idem pour 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'Calcul du déplacement orbital du point O sous epsilon_pert={epsilon_pert}...')
orbit_shift_O = orbit_shift_under_perturbation(
field_func, delta_field_func, opoint, monodromy_O
)
print(f' Déplacement orbital initial: dR={orbit_shift_O[0,0]:.6f}, dZ={orbit_shift_O[0,1]:.6f}')
print(f' Max |déplacement| le long de l’orbite: {np.max(np.linalg.norm(orbit_shift_O, axis=1)):.6f} m')
if xpoint is not None and 'monodromy_X' in dir():
print(f'Calcul du déplacement orbital du point X sous epsilon_pert={epsilon_pert}...')
orbit_shift_X = orbit_shift_under_perturbation(
field_func, delta_field_func, xpoint, monodromy_X
)
print(f' Déplacement orbital initial: dR={orbit_shift_X[0,0]:.6f}, dZ={orbit_shift_X[0,1]:.6f}')
Calcul du déplacement orbital du point O sous epsilon_pert=0.005...
Déplacement orbital initial: dR=-0.000000, dZ=-0.000000
Max |déplacement| le long de l’orbite: 0.000000 m
Calcul du déplacement orbital du point X sous epsilon_pert=0.005...
Déplacement orbital initial: 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('Déplacement orbital (m)')
ax.set_title(f'Déplacement orbital du point O sous δε={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('Déplacement orbital (m)')
ax.set_title(f'Déplacement orbital du point X sous δε={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. Diagramme de stabilité : résidu de Greene en fonction de l’amplitude de perturbation#
Nous balayons epsilon_pert de 0 à 0.02. Pour chaque amplitude :
trouver les cycles O et X
calculer le résidu de Greene \(R_G = (2 - \text{Tr}(M))/4\)
La transition \(R_G: 0 \to 1\) marque la destruction de l’îlot (transition de type KAM).
[13]:
print('Calcul du diagramme de stabilité (balayage de epsilon_h)...')
epsilon_arr = np.linspace(0.005, 0.04, 12)
greene_O = [] # point O Greene residues
greene_X = [] # point X 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 point X 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 point O 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'Calculé {len(valid_eps)} points de données')
Calcul du diagramme de stabilité (balayage de 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
Calculé 12 points de données
[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="résidu de Greene du point X")
if eps_O:
ax.plot(eps_O, rg_O, 'bo-', linewidth=2, markersize=8, label="résidu de Greene du point O")
ax.axhline(0.0, color='k', linestyle='--', alpha=0.5, label='R_G=0 (frontière elliptique)')
ax.axhline(1.0, color='k', linestyle=':', alpha=0.5, label='R_G=1 (frontière hyperbolique)')
ax.fill_between([0, 0.05], 0, 1, alpha=0.1, color='green', label='région elliptique')
ax.set_xlabel('Amplitude d’ondulation hélicoïdale ε_h')
ax.set_ylabel("Résidu de Greene R_G")
ax.set_title("Diagramme de stabilité : résidu de Greene vs amplitude de perturbation\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('Diagramme de stabilité enregistré')
Diagramme de stabilité enregistré
Résumé#
Ce notebook a démontré :
Recherche d’orbites périodiques avec Newton-Raphson sur la carte de Poincaré
Évolution jacobienne le long de l’orbite via les équations variationnelles
Matrice de monodromie \(M = J(\phi_0 + 2\pi n)\) avec \(\det(M) \simeq 1\) (préservation d’aire)
Résidu de Greene \(R_G\) comme indicateur de stabilité
Déplacement orbital sous perturbation avec la condition aux limites périodique
Diagramme de stabilité : variation de la stabilité de la chaîne d’îlots avec l’amplitude de perturbation
Résultats clés :
Les points O ont \(0 < R_G < 1\) (elliptiques, stables)
Les points X ont \(R_G < 0\) ou \(R_G > 1\) (hyperboliques, instables)
\(\det(DPm(\phi)) \simeq 1\) partout, ce qui confirme la préservation d’aire (structure hamiltonienne)
Ce port de l’analyse du Jacobien et de la monodromie depuis un workflow Julia vers Python (stellarator analytique) démontre les mêmes algorithmes physiques sans données d’équilibre privées.