#!/usr/bin/env python3 """Test actual relative equilibria with a frozen lower-depth nuclear support.""" import hashlib import importlib.util import json from datetime import datetime, timezone from pathlib import Path import numpy as np from scipy.integrate import solve_ivp import sympy as sp ROOT = Path(__file__).parent SPEC = importlib.util.spec_from_file_location('canonical_orbit', ROOT/'canonical-orbit-response-checks.py') BASE = importlib.util.module_from_spec(SPEC) SPEC.loader.exec_module(BASE) CHECKS = [] def zero(name, expression): values = list(expression) if isinstance(expression, sp.MatrixBase) else [expression] if any(sp.simplify(value) != 0 for value in values): raise AssertionError((name, expression)) CHECKS.append(name) def symbolic(): w, n, M, g = sp.symbols('w n M g', positive=True) A = sp.diag(sp.Matrix([[0, -w], [w, 0]]), sp.Matrix([[0, -n], [n, 0]])) q = sp.Matrix([1, 0, 1, 0]) v = A*q D = 1-M*(w*w+n*n) denominator = (1-M*(w+n)**2)*(1-M*(w-n)**2) c = g*n*(-1+M*(3*w*w+n*n))/denominator d = g*w*(1-M*(w*w+3*n*n))/denominator B = sp.Matrix([[0, -g*w, 0, c], [g*w, 0, d, 0], [0, -d, 0, -g*n], [-c, 0, g*n, 0]]) zero('complete rotating source response', B+M*(A*(A*B-B*A)-(A*B-B*A)*A)+g*(q*v.T-v*q.T)) F = (-(w*w+n*n)+M*(w**4+n**4+6*w*w*n*n))/denominator zero('mixed response contraction', n*c-w*d-g*F) force = (A*A+g*(A*B+B*A))*q zero('visible force balance', force[0]-(-w*w*(1+2*g*g)+g*g*F)) zero('unsupported force balance', force[2]-(-n*n*(1+2*g*g)+g*g*F)) zero('equal rotation cannot bind unsupported circle', F.subs(n, w)+2*w*w) L, y = sp.symbols('L y', positive=True) f = (-(L*L+1)+y*(L**4+6*L*L+1))/((1-y*(L+1)**2)*(1-y*(L-1)**2)) zero('response as sum and difference clocks', F.subs({w: L*n, M: y/n**2})-n**2*f) C = 2+1/g**2 polynomial = C*(L*L-1)**2*y*y-(2*C*(L*L+1)+L**4+6*L*L+1)*y+C+L*L+1 zero('equal-radius source equation', (C-f)*(1-y*(L+1)**2)*(1-y*(L-1)**2)-polynomial) q5 = sp.Matrix(sp.symbols('Q0:5')) support = sp.diag(1, 1, 1, 0, 0) radius = sp.sqrt(q5.dot(support*q5)) force = -support*q5/radius**3 potential = -1/radius zero('supported-space potential gradient', force+sp.Matrix([sp.diff(potential, z) for z in q5])) torque = q5*force.T-force*q5.T zero('supported rotational charge', torque[:3, :3]) zero('unsupported rotational charge', torque[3:, 3:]) return {'ratio_polynomial': str(polynomial), 'response': str(F), 'source_resonances': 'M*(w+n)**2=1 or M*(w-n)**2=1'} def rhs(time, state, depth, coupling, inertia, A): q, p, B, P = BASE.unpack(state, depth) support = np.zeros_like(q) support[:min(3, depth)] = q[:min(3, depth)] v = p-coupling*B@q qd = v-A@q pd = -coupling*B@v-support/(support@support)**1.5-A@p Bd = P/inertia-A@B+B@A Pd = -B-coupling*(np.outer(q, v)-np.outer(v, q))-A@P+P@A return np.concatenate([qd, pd, BASE.upper(Bd), BASE.upper(Pd)]) def equilibrium(depth, ratio, high, coupling=1.0): C = 2+1/coupling**2 qa = C*(ratio**2-1)**2 qb = -(2*C*(ratio**2+1)+ratio**4+6*ratio**2+1) qc = C+ratio**2+1 discriminant = qb*qb-4*qa*qc y = (-qb+np.sqrt(discriminant))/(2*qa) if high else 2*qc/(-qb+np.sqrt(discriminant)) nu = np.sqrt(1/((1+2*coupling**2)*(ratio**2-1))) omega = ratio*nu inertia = y/(nu*nu) A = np.zeros((depth, depth)) A[0, 1], A[1, 0], A[3, 4], A[4, 3] = -omega, omega, -nu, nu q = np.zeros(depth) q[0], q[3] = 1, 1 v = A@q basis = [] for i, j in BASE.pairs(depth): E = np.zeros((depth, depth));E[i, j], E[j, i] = 1, -1 basis.append(E) operator = np.column_stack([BASE.upper(E+inertia*(A@(A@E-E@A)-(A@E-E@A)@A)) for E in basis]) b = np.linalg.solve(operator, -coupling*BASE.upper(np.outer(q, v)-np.outer(v, q))) B = BASE.skew(b, depth) P = inertia*(A@B-B@A) p = (A+coupling*B)@q state = np.concatenate([q, p, BASE.upper(B), BASE.upper(P)]) return state, A, inertia, operator, omega, nu def energy(state, depth, coupling, inertia): q, p, B, P = BASE.unpack(state, depth) v = p-coupling*B@q radius = np.linalg.norm(q[:min(3, depth)]) return v@v/2+np.sum(BASE.upper(P)**2)/(2*inertia)+np.sum(BASE.upper(B)**2)/2-1/radius if __name__ == '__main__': symbols = symbolic() print(json.dumps(symbols, indent=2), flush=True) cases = [] structural = [{'depth': d, 'nuclear_support_rank': min(3, d), 'unsupported_rotation_planes': (d-min(3, d))//2, 'two-plane_branch_available': d >= 5} for d in (2, 3, 4, 5, 6, 10)] growth_runs = [] conservation = [] for depth in (2, 3, 4, 5, 6, 10): rng = np.random.default_rng(20261011+depth) size = depth*(depth-1)//2 state = rng.normal(size=2*depth+2*size)*0.1 state[0] += 1 derivative = rhs(0, state, depth, 0.7, 1.3, np.zeros((depth, depth))) step = 2e-6 high, low = state+step*derivative, state-step*derivative energy_defect = abs((energy(high, depth, 0.7, 1.3)-energy(low, depth, 0.7, 1.3))/(2*step)) def charge(y): q, p, B, P = BASE.unpack(y, depth) return np.outer(q, p)-np.outer(p, q)-(B@P-P@B) torque = (charge(high)-charge(low))/(2*step) rank = min(3, depth) allowed_defect = max(np.linalg.norm(torque[:rank, :rank]), np.linalg.norm(torque[rank:, rank:])) if max(energy_defect, allowed_defect) > 2e-8: raise AssertionError(('energy or applicable rotational charge', depth, energy_defect, allowed_defect)) conservation.append({'depth': depth, 'energy_derivative_defect': energy_defect, 'applicable_charge_derivative_defect': allowed_defect, 'full_charge_derivative_norm': float(np.linalg.norm(torque))}) for ratio in (1.5, 2.0, 2.5, 3.0, 4.0): for high in (False, True): for depth in (5, 6, 10): state, A, inertia, operator, omega, nu = equilibrium(depth, ratio, high) residual = np.max(np.abs(rhs(0, state, depth, 1.0, inertia, A))) if residual > 2e-10: raise AssertionError(('equilibrium', ratio, high, depth, residual)) jac = np.column_stack([np.imag(rhs(0, state.astype(complex)+1e-25j*np.eye(len(state))[i], depth, 1.0, inertia, A))/1e-25 for i in range(len(state))]) direction = np.random.default_rng(20261011+depth).normal(size=len(state)) direction /= np.linalg.norm(direction) step = 1e-6 finite = (rhs(0, state+step*direction, depth, 1.0, inertia, A)- rhs(0, state-step*direction, depth, 1.0, inertia, A))/(2*step) jacobian_defect = np.max(np.abs(finite-jac@direction)) if jacobian_defect > 2e-8: raise AssertionError(('independent Jacobian', ratio, high, depth, jacobian_defect)) values, vectors = np.linalg.eig(jac) growth = float(max(values.real)) cases.append({'depth': depth, 'ratio': ratio, 'branch': 'high-inertia' if high else 'low-inertia', 'source_inertia': inertia, 'orbit_clock': omega, 'unsupported_clock': nu, 'source_operator_condition': float(np.linalg.cond(operator)), 'equilibrium_residual': residual, 'max_real_eigenvalue': growth, 'independent_jacobian_defect': jacobian_defect, 'phase_space_dimension': len(state), 'eigenvalues': [[float(z.real), float(z.imag)] for z in values]}) if ratio == 2 and depth == 6 and growth > 1e-5: idx = int(np.argmax(values.real)) mode = vectors[:, idx] initial = mode.real/np.linalg.norm(mode.real) duration = 4/growth amplitude = 1e-8 solution = solve_ivp(lambda t, y: rhs(t, y, depth, 1.0, inertia, A), (0, duration), state+amplitude*initial, method='DOP853', rtol=1e-12, atol=1e-13) if not solution.success: raise AssertionError(solution.message) # Compare with the complete real linear flow, including any imaginary frequency. from scipy.linalg import expm expected = expm(jac*duration)@initial actual = (solution.y[:, -1]-state)/amplitude relative = np.linalg.norm(actual-expected)/np.linalg.norm(expected) if relative > 2e-3: raise AssertionError(('nonlinear unstable mode', high, relative)) growth_runs.append({'branch': 'high-inertia' if high else 'low-inertia', 'depth': depth, 'ratio': ratio, 'duration': duration, 'linear_amplification': float(np.linalg.norm(expected)), 'relative_growth_error': float(relative), 'energy_error': abs(energy(solution.y[:, -1], depth, 1, inertia)- energy(solution.y[:, 0], depth, 1, inertia)), 'rhs_evaluations': solution.nfev}) path = Path(__file__) stamp = datetime.now(timezone.utc).strftime('%Y%m%dT%H%M%S.%fZ') report = {'created_utc': stamp, 'sources': { 'experiments/cylindrical-orbit-response-checks.py': hashlib.sha256(path.read_bytes()).hexdigest(), 'experiments/canonical-orbit-response-checks.py': hashlib.sha256((ROOT/'canonical-orbit-response-checks.py').read_bytes()).hexdigest()}, 'exact_checks': CHECKS, 'symbolic': symbols, 'structural_controls': structural, 'equilibria': cases, 'nonlinear_growth': growth_runs, 'conservation': conservation} destination = path.with_name('cylindrical-orbit-response-results-'+stamp+'.json') destination.write_text(json.dumps(report, indent=2)+'\n') print(json.dumps({'report': str(destination), 'exact_checks': len(CHECKS), 'structural_controls': structural, 'cases': len(cases), 'resolved_unstable_cases': sum(c['max_real_eigenvalue'] > 1e-5 for c in cases), 'max_equilibrium_residual': max(c['equilibrium_residual'] for c in cases), 'max_independent_jacobian_defect': max(c['independent_jacobian_defect'] for c in cases), 'nonlinear_growth': growth_runs}, indent=2), flush=True)