#!/usr/bin/env python3 """Audit the confining modified-Mathieu/DCHE spectral realization. The dimensionless mechanical convention is H_kappa = -d^2/dx^2 + 2*kappa*cosh(2*x), kappa > 0, H_kappa Phi = A Phi, Phi in L^2(R). With s=exp(x), zeta=-2*sqrt(kappa)*s, and Phi = s^(1/2) exp[-sqrt(kappa)*(s+s^(-1))] Y(zeta), the equation becomes the DLMF doubly confluent Heun equation Y'' + (delta/zeta^2 + gamma/zeta + 1)Y' + (alpha*zeta-q_D)/zeta^2 Y = 0 on the negative zeta ray, with alpha=1, gamma=2, delta=-4*kappa, q_D=2*kappa-A-1/4. Two numerically independent representations are compared: 1. parity-reduced uniform-grid Jacobi matrices, refined by step halving and second-order Richardson extrapolation; 2. direct DCHE shooting from the algebraic formal solution at zeta=infinity, refined independently in endpoint distance, formal-series order, and Runge--Kutta step. The program also starts the algebraic solution at zeta=0 and checks Abel's matching-point invariant away from the spectrum. All reported shifts are empirical floating-point diagnostics, not certified enclosures. Requirements: CPython 3.9+ and NumPy 1.24+. """ from __future__ import annotations import argparse import csv from dataclasses import dataclass import math from pathlib import Path import platform from typing import Iterable, Sequence import numpy as np MIN_KAPPA = 1.0 MAX_KAPPA = 9.0 DEFAULT_KAPPA = 4.0 def require(condition: bool, message: str) -> None: """Raise an optimization-safe exception when an audit fails.""" if not condition: raise RuntimeError(message) def potential(x: np.ndarray | float, kappa: float) -> np.ndarray | float: """Return 2*kappa*cosh(2*x).""" return 2.0 * kappa * np.cosh(2.0 * np.asarray(x)) @dataclass(frozen=True) class Settings: """All cutoffs for one reproducibility profile.""" grid_coarse_cells: int domain: float series_coarse_order: int series_fine_order: int s_coarse_min: float s_coarse_max: float s_fine_min: float s_fine_max: float zeta_coarse_step_factor: float zeta_fine_step_factor: float @dataclass(frozen=True) class GridSpectrum: """Richardson-refined parity spectra and mesh diagnostics.""" even: np.ndarray odd: np.ndarray fine_even: np.ndarray fine_odd: np.ndarray maximum_shift: float sturm_width: float @dataclass(frozen=True) class SpectralRow: """One ordered eigenvalue in the two numerical representations.""" index: int parity: str grid: float dche: float dche_refinement_shift: float @property def cross_residual(self) -> float: return abs(self.grid - self.dche) def settings_for_mode(mode: str) -> Settings: """Return frozen default or high-refinement settings.""" if mode == "default": return Settings( grid_coarse_cells=2400, domain=5.0, series_coarse_order=10, series_fine_order=20, s_coarse_min=0.050, s_coarse_max=20.0, s_fine_min=0.035, s_fine_max=40.0, zeta_coarse_step_factor=0.005, zeta_fine_step_factor=0.0025, ) if mode == "high": return Settings( grid_coarse_cells=4800, domain=5.0, series_coarse_order=20, series_fine_order=28, s_coarse_min=0.035, s_coarse_max=40.0, s_fine_min=0.025, s_fine_max=60.0, zeta_coarse_step_factor=0.0025, zeta_fine_step_factor=0.00125, ) raise ValueError(f"unknown mode {mode!r}") def dche_parameters(kappa: float, energy: float) -> tuple[float, float, float, float]: """Return (alpha, gamma, delta, q_D) in the DLMF convention.""" return 1.0, 2.0, -4.0 * kappa, 2.0 * kappa - energy - 0.25 def grid_parity_blocks( kappa: float, domain: float, cells: int, ) -> tuple[ tuple[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray], ]: """Return even and odd Jacobi blocks on the half-line grid.""" require(domain >= 4.0, "grid domain is too short for the calibrated range") require(cells >= 400, "coordinate grid is too coarse") step = domain / cells kinetic_diagonal = 2.0 / step**2 kinetic_off_diagonal = -1.0 / step**2 even_x = step * np.arange(cells, dtype=float) even_diagonal = kinetic_diagonal + potential(even_x, kappa) even_off = np.full(cells - 1, kinetic_off_diagonal) # Orthonormal reduction of the full-line couplings adjacent to x=0. even_off[0] *= math.sqrt(2.0) odd_x = step * np.arange(1, cells, dtype=float) odd_diagonal = kinetic_diagonal + potential(odd_x, kappa) odd_off = np.full(cells - 2, kinetic_off_diagonal) return ( (np.asarray(even_diagonal), even_off), (np.asarray(odd_diagonal), odd_off), ) def sturm_count( diagonal: np.ndarray, off_diagonal: np.ndarray, value: float, ) -> int: """Count Jacobi eigenvalues strictly below value by an LDL sequence.""" scale = max(1.0, float(np.max(np.abs(diagonal)))) pivot_floor = 8.0 * np.finfo(float).eps * scale pivot = float(diagonal[0] - value) if abs(pivot) < pivot_floor: pivot = -pivot_floor count = int(pivot < 0.0) for index in range(1, diagonal.size): pivot = float( diagonal[index] - value - off_diagonal[index - 1] ** 2 / pivot ) if abs(pivot) < pivot_floor: pivot = -pivot_floor count += int(pivot < 0.0) return count def tridiagonal_bounds( diagonal: np.ndarray, off_diagonal: np.ndarray, ) -> tuple[float, float]: """Return Gershgorin bounds for a real symmetric tridiagonal matrix.""" radii = np.zeros_like(diagonal) radii[:-1] += np.abs(off_diagonal) radii[1:] += np.abs(off_diagonal) scale = max(1.0, float(np.max(np.abs(diagonal)))) pad = 8.0 * np.finfo(float).eps * scale return ( float(np.min(diagonal - radii) - pad), float(np.max(diagonal + radii) + pad), ) def lowest_tridiagonal_eigenvalues( diagonal: np.ndarray, off_diagonal: np.ndarray, count: int, ) -> tuple[np.ndarray, float]: """Return the first ``count`` Jacobi eigenvalues by Sturm bisection.""" require( diagonal.size == off_diagonal.size + 1, "invalid tridiagonal dimensions", ) require(1 <= count < diagonal.size, "invalid requested eigenvalue count") global_lower, global_upper = tridiagonal_bounds(diagonal, off_diagonal) eigenvalues: list[float] = [] largest_width = 0.0 for target_index in range(count): lower = global_lower upper = global_upper for _ in range(96): midpoint = 0.5 * (lower + upper) if midpoint == lower or midpoint == upper: break if sturm_count(diagonal, off_diagonal, midpoint) <= target_index: lower = midpoint else: upper = midpoint eigenvalues.append(0.5 * (lower + upper)) largest_width = max(largest_width, upper - lower) return np.asarray(eigenvalues), largest_width def grid_spectrum( kappa: float, domain: float, coarse_cells: int, count_each: int, ) -> GridSpectrum: """Refine both parity spectra by grid step halving and Richardson.""" coarse_blocks = grid_parity_blocks(kappa, domain, coarse_cells) fine_blocks = grid_parity_blocks(kappa, domain, 2 * coarse_cells) coarse_even, width_ce = lowest_tridiagonal_eigenvalues( *coarse_blocks[0], count_each ) coarse_odd, width_co = lowest_tridiagonal_eigenvalues( *coarse_blocks[1], count_each ) fine_even, width_fe = lowest_tridiagonal_eigenvalues( *fine_blocks[0], count_each ) fine_odd, width_fo = lowest_tridiagonal_eigenvalues( *fine_blocks[1], count_each ) even = (4.0 * fine_even - coarse_even) / 3.0 odd = (4.0 * fine_odd - coarse_odd) / 3.0 maximum_shift = max( float(np.max(np.abs(even - fine_even))), float(np.max(np.abs(odd - fine_odd))), ) return GridSpectrum( even=even, odd=odd, fine_even=fine_even, fine_odd=fine_odd, maximum_shift=maximum_shift, sturm_width=max(width_ce, width_co, width_fe, width_fo), ) def zero_algebraic_series( zeta: float, kappa: float, energy: float, order: int, ) -> tuple[float, float]: """Evaluate the truncated algebraic DCHE series at zeta=0.""" _, _, delta, q_d = dche_parameters(kappa, energy) require(zeta < 0.0, "the physical DCHE ray has zeta < 0") require(order >= 2, "endpoint series order must be at least two") coefficients = [1.0] previous = 0.0 current = 1.0 for n in range(order - 1): following = -( (n * (n + 1.0) - q_d) * current + n * previous ) / (delta * (n + 1.0)) coefficients.append(following) previous, current = current, following value = 0.0 derivative = 0.0 power = 1.0 for n, coefficient in enumerate(coefficients): value += coefficient * power if n: derivative += n * coefficient * power / zeta power *= zeta return value, derivative def infinity_algebraic_series( zeta: float, kappa: float, energy: float, order: int, ) -> tuple[float, float]: """Evaluate the truncated algebraic DCHE series at zeta=infinity.""" _, _, delta, q_d = dche_parameters(kappa, energy) require(zeta < 0.0, "the physical DCHE ray has zeta < 0") require(order >= 2, "endpoint series order must be at least two") coefficients = [1.0] previous = 0.0 current = 1.0 for n in range(order - 1): following = ( (n * (n + 1.0) - q_d) * current - delta * n * previous ) / (n + 1.0) coefficients.append(following) previous, current = current, following value = 0.0 derivative = 0.0 for n, coefficient in enumerate(coefficients): value += coefficient * zeta ** (-n - 1) derivative -= (n + 1.0) * coefficient * zeta ** (-n - 2) return value, derivative def dche_rhs( zeta: float, value: float, derivative: float, delta: float, q_d: float, ) -> tuple[float, float]: """Return the first-order form of the specialized DLMF DCHE.""" second = -( 1.0 + 2.0 / zeta + delta / zeta**2 ) * derivative - ( 1.0 / zeta - q_d / zeta**2 ) * value return derivative, second def rk4_segment( start: float, state: tuple[float, float], stop: float, max_step: float, delta: float, q_d: float, ) -> tuple[float, float]: """Propagate one DCHE solution over a real-zeta segment.""" require(start < 0.0 and stop < 0.0, "DCHE propagation cannot cross zero") steps = max(1, int(math.ceil(abs(stop - start) / max_step))) step = (stop - start) / steps zeta = start value, derivative = state for _ in range(steps): k1v, k1d = dche_rhs(zeta, value, derivative, delta, q_d) k2v, k2d = dche_rhs( zeta + 0.5 * step, value + 0.5 * step * k1v, derivative + 0.5 * step * k1d, delta, q_d, ) k3v, k3d = dche_rhs( zeta + 0.5 * step, value + 0.5 * step * k2v, derivative + 0.5 * step * k2d, delta, q_d, ) k4v, k4d = dche_rhs( zeta + step, value + step * k3v, derivative + step * k3d, delta, q_d, ) value += step * (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0 derivative += step * (k1d + 2.0 * k2d + 2.0 * k3d + k4d) / 6.0 zeta += step require( math.isfinite(value) and math.isfinite(derivative), "DCHE propagation overflowed; retune the calibrated range", ) return value, derivative def propagate_to_targets( start: float, state: tuple[float, float], targets: Sequence[float], max_step: float, delta: float, q_d: float, ) -> dict[float, tuple[float, float]]: """Propagate successively to monotone target points.""" if not targets: return {} direction = math.copysign(1.0, targets[0] - start) ordered = sorted(targets, reverse=direction < 0.0) require( all(math.copysign(1.0, target - start) == direction for target in ordered), "targets must lie on one side of the initial point", ) current = start current_state = state result: dict[float, tuple[float, float]] = {} for target in ordered: current_state = rk4_segment( current, current_state, target, max_step, delta, q_d, ) current = target result[target] = current_state return result def dche_parity_residual( energy: float, kappa: float, parity: str, order: int, s_max: float, step_factor: float, ) -> float: """Return a normalized midpoint residual from the infinity endpoint.""" _, _, delta, q_d = dche_parameters(kappa, energy) root = math.sqrt(kappa) zeta_start = -2.0 * root * s_max zeta_midpoint = -2.0 * root state = infinity_algebraic_series( zeta_start, kappa, energy, order, ) value, derivative = rk4_segment( zeta_start, state, zeta_midpoint, max_step=2.0 * root * step_factor, delta=delta, q_d=q_d, ) scaled_derivative = zeta_midpoint * derivative norm = math.hypot(value, scaled_derivative) require(norm > 0.0, "vanishing DCHE midpoint state") if parity == "odd": return value / norm if parity == "even": return (scaled_derivative + 0.5 * value) / norm raise ValueError(f"unknown parity {parity!r}") def bracketed_root( function, guess: float, initial_width: float, ) -> float: """Bracket a nearby simple root and refine it by bisection.""" width = initial_width lower = guess - width upper = guess + width f_lower = function(lower) f_upper = function(upper) for _ in range(8): if f_lower == 0.0: return lower if f_upper == 0.0: return upper if f_lower * f_upper < 0.0: break width *= 1.6 lower = guess - width upper = guess + width f_lower = function(lower) f_upper = function(upper) else: raise RuntimeError(f"could not bracket DCHE root near A={guess:.8g}") for _ in range(72): midpoint = 0.5 * (lower + upper) if midpoint == lower or midpoint == upper: break f_midpoint = function(midpoint) if f_midpoint == 0.0: return midpoint if f_lower * f_midpoint < 0.0: upper = midpoint f_upper = f_midpoint else: lower = midpoint f_lower = f_midpoint return 0.5 * (lower + upper) def dche_root( guess: float, kappa: float, parity: str, order: int, s_max: float, step_factor: float, ) -> float: """Find one parity root of the DCHE endpoint connection problem.""" def residual(energy: float) -> float: return dche_parity_residual( energy, kappa, parity, order, s_max, step_factor, ) return bracketed_root( residual, guess, initial_width=max(0.15, 0.20 * math.sqrt(kappa)), ) def ordered_grid_levels( grid: GridSpectrum, levels: int, ) -> list[tuple[float, str]]: """Merge the even and odd grid spectra into spectral order.""" merged = [(float(value), "even") for value in grid.even] merged += [(float(value), "odd") for value in grid.odd] merged.sort(key=lambda item: item[0]) return merged[:levels] def build_rows( kappa: float, grid: GridSpectrum, settings: Settings, levels: int, ) -> list[SpectralRow]: """Compute refined DCHE roots from the independent grid guesses.""" rows: list[SpectralRow] = [] for index, (guess, parity) in enumerate(ordered_grid_levels(grid, levels)): coarse = dche_root( guess, kappa, parity, settings.series_coarse_order, settings.s_coarse_max, settings.zeta_coarse_step_factor, ) fine = dche_root( guess, kappa, parity, settings.series_fine_order, settings.s_fine_max, settings.zeta_fine_step_factor, ) rows.append( SpectralRow( index=index, parity=parity, grid=guess, dche=fine, dche_refinement_shift=abs(fine - coarse), ) ) return rows def abel_invariant_spread( kappa: float, energy: float, order: int, s_min: float, s_max: float, step_factor: float, ) -> tuple[float, list[tuple[float, float]]]: """Check the Abel-normalized Wronskian at three matching points.""" _, gamma, delta, q_d = dche_parameters(kappa, energy) root = math.sqrt(kappa) zeta_left = -2.0 * root * s_min zeta_right = -2.0 * root * s_max targets = [-1.6 * root, -2.0 * root, -2.5 * root] max_step = 2.0 * root * step_factor left_state = zero_algebraic_series( zeta_left, kappa, energy, order, ) right_state = infinity_algebraic_series( zeta_right, kappa, energy, order, ) left_values = propagate_to_targets( zeta_left, left_state, targets, max_step, delta, q_d, ) right_values = propagate_to_targets( zeta_right, right_state, targets, max_step, delta, q_d, ) invariants: list[tuple[float, float]] = [] for zeta in targets: left_value, left_derivative = left_values[zeta] right_value, right_derivative = right_values[zeta] wronskian = ( left_value * right_derivative - left_derivative * right_value ) invariant = ( math.exp(zeta - delta / zeta) * zeta**gamma * wronskian ) invariants.append((zeta, invariant)) values = [value for _, value in invariants] scale = max(max(abs(value) for value in values), np.finfo(float).tiny) spread = (max(values) - min(values)) / scale return abs(spread), invariants def deep_well_energy(index: int, kappa: float) -> float: """Return the fixed-index large-kappa expansion through O(1).""" return ( 2.0 * kappa + (4.0 * index + 2.0) * math.sqrt(kappa) + (2.0 * index * index + 2.0 * index + 1.0) / 4.0 ) def turning_points(energy: float, kappa: float) -> tuple[float, float]: """Return the two real turning points for energy above the minimum.""" require(energy > 2.0 * kappa, "turning points require A > 2*kappa") positive = 0.5 * math.acosh(energy / (2.0 * kappa)) return -positive, positive def write_main_csv(path: Path, rows: Iterable[SpectralRow]) -> None: """Write the two-method spectral table.""" path.parent.mkdir(parents=True, exist_ok=True) with path.open("w", newline="", encoding="utf-8") as stream: writer = csv.writer(stream) writer.writerow( [ "index", "parity", "grid_richardson", "dche_shooting", "absolute_difference", "dche_refinement_shift", ] ) for row in rows: writer.writerow( [ row.index, row.parity, f"{row.grid:.16e}", f"{row.dche:.16e}", f"{row.cross_residual:.16e}", f"{row.dche_refinement_shift:.16e}", ] ) def write_figure_data( prefix: Path, kappa: float, rows: Sequence[SpectralRow], ) -> None: """Export the exact data embedded in the companion TikZ figure.""" prefix.parent.mkdir(parents=True, exist_ok=True) potential_path = Path(f"{prefix}-potential.csv") levels_path = Path(f"{prefix}-levels.csv") ray_path = Path(f"{prefix}-ray.csv") with potential_path.open("w", newline="", encoding="utf-8") as stream: writer = csv.writer(stream) writer.writerow(["x", "potential"]) for x in np.linspace(-1.15, 1.15, 93): writer.writerow([f"{x:.8f}", f"{float(potential(x, kappa)):.12f}"]) with levels_path.open("w", newline="", encoding="utf-8") as stream: writer = csv.writer(stream) writer.writerow(["index", "parity", "energy"]) for row in rows[:4]: writer.writerow([row.index, row.parity, f"{row.dche:.12f}"]) ground = rows[0].dche x_left, x_right = turning_points(ground, kappa) root = math.sqrt(kappa) with ray_path.open("w", newline="", encoding="utf-8") as stream: writer = csv.writer(stream) writer.writerow(["point", "x", "s", "zeta"]) for label, x in ( ("left_turning", x_left), ("midpoint", 0.0), ("right_turning", x_right), ): s_value = math.exp(x) writer.writerow( [label, f"{x:.12f}", f"{s_value:.12f}", f"{-2.0*root*s_value:.12f}"] ) def parse_arguments() -> argparse.Namespace: """Parse and validate the command line.""" parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--mode", choices=("default", "high"), default="default") parser.add_argument("--kappa", type=float, default=DEFAULT_KAPPA) parser.add_argument("--levels", type=int, default=6) parser.add_argument("--csv", type=Path) parser.add_argument( "--figure-data", type=Path, help="prefix for potential, level, and DCHE-ray CSV exports", ) arguments = parser.parse_args() if not MIN_KAPPA <= arguments.kappa <= MAX_KAPPA: parser.error( f"--kappa must lie in the calibrated interval " f"[{MIN_KAPPA:g}, {MAX_KAPPA:g}]" ) if not 2 <= arguments.levels <= 8: parser.error("--levels must lie between 2 and 8") return arguments def main() -> None: """Run the complete two-representation calibration.""" arguments = parse_arguments() settings = settings_for_mode(arguments.mode) count_each = (arguments.levels + 1) // 2 grid = grid_spectrum( arguments.kappa, settings.domain, settings.grid_coarse_cells, count_each, ) rows = build_rows( arguments.kappa, grid, settings, arguments.levels, ) audit_energy = 0.5 * (rows[0].dche + rows[1].dche) abel_spread, invariants = abel_invariant_spread( arguments.kappa, audit_energy, settings.series_fine_order, settings.s_fine_min, settings.s_fine_max, settings.zeta_fine_step_factor, ) alpha, gamma, delta, _ = dche_parameters(arguments.kappa, rows[0].dche) print("Modified-Mathieu / doubly confluent Heun calibration") print("H = -d^2/dx^2 + 2*kappa*cosh(2x), Phi in L^2(R)") print( f"Python {platform.python_version()}; NumPy {np.__version__}; " f"calibrated {MIN_KAPPA:g} <= kappa <= {MAX_KAPPA:g}" ) print( f"profile={arguments.mode}; kappa={arguments.kappa:g}; " f"grid={settings.grid_coarse_cells}->{2*settings.grid_coarse_cells} " f"cells on [0,{settings.domain:g}]" ) print() print("Exact DCHE passport") print(f" alpha_D={alpha:g}, gamma_D={gamma:g}, delta_D={delta:g}") print(" q_D=2*kappa-A-1/4") print(" physical ray: zeta in (-infinity,0)") print() print("Ordered spectrum") print( f"{'n':>3s} {'parity':>7s} {'grid/Richardson':>19s} " f"{'DCHE shooting':>19s} {'|difference|':>14s} {'DCHE shift':>12s}" ) for row in rows: print( f"{row.index:3d} {row.parity:>7s} {row.grid:19.12f} " f"{row.dche:19.12f} {row.cross_residual:14.3e} " f"{row.dche_refinement_shift:12.3e}" ) print() print("Independent diagnostics") print(f" grid Richardson maximum shift: {grid.maximum_shift:.3e}") print( " Sturm bisection stagnation width (not an error bound): " f"{grid.sturm_width:.3e}" ) print(f" maximum grid/DCHE energy residual: {max(row.cross_residual for row in rows):.3e}") print(f" maximum DCHE refinement shift: {max(row.dche_refinement_shift for row in rows):.3e}") print(f" Abel-invariant relative spread: {abel_spread:.3e}") for zeta, invariant in invariants: print(f" zeta={zeta:8.3f}: invariant={invariant:+.12e}") print() print("Deep-well calibration through O(1)") for index, row in enumerate(rows[:4]): approximation = deep_well_energy(index, arguments.kappa) print( f" n={index}: exact={row.dche:.10f}, " f"asymptotic={approximation:.10f}, " f"difference={row.dche-approximation:+.3e}" ) print() print("Caveat: refinement shifts and method gaps are empirical double-precision diagnostics.") print("The program does not evaluate a Painleve tau function or an NS period.") if arguments.csv is not None: write_main_csv(arguments.csv, rows) print(f"Wrote {arguments.csv}") if arguments.figure_data is not None: write_figure_data(arguments.figure_data, arguments.kappa, rows) print(f"Wrote {arguments.figure_data}-potential.csv") print(f"Wrote {arguments.figure_data}-levels.csv") print(f"Wrote {arguments.figure_data}-ray.csv") if __name__ == "__main__": main()