back to table · edit · history · where entries came from · files · download
6915 bytes, as of the version from 2026-10-27 01:01 (current). Recorded here, not run.
"""Independent Mathieu endpoint ODE shooting and normalization checks.
python3 -m pip install numberdb mpmath
python3 mathieu_checks.py --output /path/to/elliptic/run
Taylor solutions and their derivatives with respect to the eigenvalue are
continued from 0 to pi/2. The boundary Wronskian gives the squared norm.
Two step/order/precision profiles must agree. These are numerical agreement
checks, not certified ODE remainder estimates. No website writes occur.
"""
import argparse
import json
from decimal import Decimal
from pathlib import Path
import mpmath as mp
from mathieu_reference import (DIAGNOSTICS, ORDERS, PARAMETERS, PROFILES,
MathieuReference, finite_section)
def shoot(kind, eigenvalue, coupling, digits, steps, order):
with mp.workdps(digits):
step = mp.pi / (2 * steps)
position = mp.mpf(0)
state = [mp.mpf(1 if kind == "ce" else 0),
mp.mpf(0 if kind == "ce" else 1), mp.mpf(0), mp.mpf(0)]
for interval in range(steps):
cosine = [mp.cos(2 * position), -2 * mp.sin(2 * position)]
for index in range(2, order):
cosine.append(-4 * cosine[index - 2] / (index * (index - 1)))
potential = [2 * coupling * value for value in cosine]
potential[0] -= eigenvalue
solution = state[:2]
variation = state[2:]
for index in range(order - 1):
denominator = (index + 1) * (index + 2)
solution.append(mp.fsum(potential[column] * solution[index - column]
for column in range(index + 1)) / denominator)
variation.append((mp.fsum(potential[column] * variation[index - column]
for column in range(index + 1))
- solution[index]) / denominator)
state = [mp.polyval(list(reversed(solution)), step),
mp.polyval(list(reversed([index * solution[index]
for index in range(1, len(solution))])), step),
mp.polyval(list(reversed(variation)), step),
mp.polyval(list(reversed([index * variation[index]
for index in range(1, len(variation))])), step)]
position += step
return state
def normalized(kind, order, initial, profile):
digits, steps, degree = profile
with mp.workdps(digits):
neumann = order % 2 == (0 if kind == "ce" else 1)
eigenvalue = +initial
for iteration in range(2):
value, slope, variation, variation_slope = shoot(
kind, eigenvalue, mp.mpf(1), digits, steps, degree)
residual = slope if neumann else value
jacobian = variation_slope if neumann else variation
eigenvalue -= residual / jacobian
value, slope, variation, variation_slope = shoot(
kind, eigenvalue, mp.mpf(1), digits, steps, degree)
norm = slope * variation - value * variation_slope
if not mp.isfinite(norm) or norm <= 0:
raise ArithmeticError("Invalid independent norm")
scale = mp.sqrt(mp.pi / (4 * norm))
return {"eigenvalue": +eigenvalue, "0": +scale if kind == "ce" else mp.mpf(0),
"1/2": +(scale * value) if neumann else mp.mpf(0),
"boundary_residual": +(slope if neumann else value)}
def main(output):
output.mkdir(parents=True, exist_ok=True)
ode_profiles = ((160, 12, 140), (190, 16, 170))
with mp.workdps(170):
for kind in ("ce", "se"):
for order in (1, 2, 3):
value, slope, variation, variation_slope = shoot(
kind, mp.mpf(order * order), mp.mpf(0), 160, 12, 140)
expected_value = mp.cos(order * mp.pi / 2) if kind == "ce" else mp.sin(order * mp.pi / 2) / order
expected_slope = -order * mp.sin(order * mp.pi / 2) if kind == "ce" else mp.cos(order * mp.pi / 2)
expected_norm = mp.pi / 4 if kind == "ce" else mp.pi / (4 * order * order)
if max(abs(value - expected_value), abs(slope - expected_slope),
abs(slope * variation - value * variation_slope - expected_norm)) > mp.mpf("1e-130"):
raise ArithmeticError("Zero-coupling ODE control failed")
independent = {}
for kind in ("ce", "se"):
orders = range(11) if kind == "ce" else range(1, 7)
for order in orders:
initial = finite_section(kind, order % 2, *PROFILES[-1])[order]["eigenvalue"]
profiles = [normalized(kind, order, initial, profile) for profile in ode_profiles]
with mp.workdps(210):
for key in ("eigenvalue", "0", "1/2"):
if abs(profiles[0][key] - profiles[1][key]) > mp.mpf("1e-125"):
raise ArithmeticError("Independent ODE profiles disagree")
independent[kind, order] = profiles[-1]
print("ODE", kind, order, "done", flush=True)
for tid in PARAMETERS:
generator = MathieuReference(tid)
rows = []
for params in generator.enumerate():
written = str(generator.value(params, 100))
kind = "se" if tid == "T438" else "ce"
order = int(params["n"])
key = "eigenvalue" if tid == "T436" else params[generator.parameters[-1]]
with mp.workdps(210):
value = independent[kind, order][key]
unit = mp.power(10, Decimal(written).as_tuple().exponent)
if written == "0":
if value != 0:
raise ArithmeticError("Exact-zero mismatch")
discrepancy = mp.mpf(0)
else:
discrepancy = abs(value - mp.mpf(written)) / unit
if discrepancy >= 1:
raise ArithmeticError("Independent ODE outside output interval")
rows.append({"parameters": params, "written": written,
"ode_error_ulps": mp.nstr(discrepancy, 20)})
report = {"tid": tid, "rows": rows, "checked": len(rows),
"fourier_profiles": PROFILES, "ode_profiles": ode_profiles,
"rigour": "heuristic (agreement-checked)", "controls_passed": True}
directory = output / tid
directory.mkdir(parents=True, exist_ok=True)
(directory / "mathieu-checks.json").write_text(json.dumps(report, indent=2) + "\n")
(output / "mathieu-diagnostics.json").write_text(json.dumps(DIAGNOSTICS, indent=2) + "\n")
if __name__ == "__main__":
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--output", type=Path, required=True)
arguments = parser.parse_args()
main(arguments.output)