diff --git a/test/segmenter-recurrence-oracle.test.ts b/test/segmenter-recurrence-oracle.test.ts new file mode 100644 index 00000000..097034a9 --- /dev/null +++ b/test/segmenter-recurrence-oracle.test.ts @@ -0,0 +1,83 @@ +import { describe, expect, it } from 'vitest'; +import { makeSpringValueSampler } from '../src/internal/solver.js'; +import type { SpringParams } from '../src/spring.js'; + +/** + * RED-оракул #227 (шаг 1-2 из issue): root-matched recurrence обязана + * воспроизводить аналитическую сетку с дрейфом на порядки ниже бюджета + * реконструкции (tolerance/2). Референс НЕ использует production-коэффициенты: + * S и P выводятся здесь независимо из полюсов (q₁+q₂, q₁q₂), а сравнение идёт + * против makeSpringValueSampler — двух независимых форм одной ОДУ. + */ +const SWEEP: readonly [label: string, p: SpringParams, v0: number][] = [ + ['underdamped wobbly', { mass: 1, stiffness: 180, damping: 12 }, 0], + ['underdamped v0=3', { mass: 1, stiffness: 180, damping: 12 }, 3], + ['near-critical снизу', { mass: 1, stiffness: 100, damping: 19.99 }, 0], + ['critical', { mass: 1, stiffness: 100, damping: 20 }, 1.5], + ['near-critical сверху', { mass: 1, stiffness: 100, damping: 20.01 }, 0], + ['overdamped ζ=2', { mass: 1, stiffness: 100, damping: 40 }, -2], + ['stiff', { mass: 1, stiffness: 210, damping: 20 }, 0], +]; + +function recurrenceGrid(p: SpringParams, v0: number, T: number, n: number): number[] { + const omega0 = Math.sqrt(p.stiffness / p.mass); + const zeta = p.damping / (2 * p.mass * omega0); + const h = T / n; + let S: number; + let P: number; + if (zeta < 1) { + const zw = zeta * omega0; + const omegaD = omega0 * Math.sqrt(1 - zeta * zeta); + S = 2 * Math.exp(-zw * h) * Math.cos(omegaD * h); + P = Math.exp(-2 * zw * h); + } else if (zeta === 1) { + const q = Math.exp(-omega0 * h); + S = 2 * q; + P = q * q; + } else { + const sqrtTerm = Math.sqrt(zeta * zeta - 1); + const r2 = -omega0 * (zeta + sqrtTerm); + const r1 = (omega0 * omega0) / r2; + const q1 = Math.exp(r1 * h); + const q2 = Math.exp(r2 * h); + S = q1 + q2; + P = q1 * q2; + } + const sample = makeSpringValueSampler(p, v0); + const ys = new Array(n + 1); + ys[0] = 0; + ys[1] = sample(h); + // y = x − 1 удовлетворяет однородной ОДУ ⇒ рекуррентность точна для y. + let prev = ys[0]! - 1; + let curr = ys[1]! - 1; + for (let i = 2; i <= n; i++) { + const next = S * curr - P * prev; + prev = curr; + curr = next; + ys[i] = 1 + next; + } + return ys; +} + +describe('root-matched recurrence — независимый оракул против аналитики (#227)', () => { + for (const [label, p, v0] of SWEEP) { + it(`${label}: дрейф на сетках 256/1024/4096 << tolerance/2`, () => { + const omega0 = Math.sqrt(p.stiffness / p.mass); + const zeta = p.damping / (2 * p.mass * omega0); + const slow = zeta >= 1 ? zeta - Math.sqrt(zeta * zeta - 1) : zeta; + const T = Math.log(100) / (omega0 * Math.max(slow, 1e-6)); + const sample = makeSpringValueSampler(p, v0); + for (const n of [256, 1024, 4096]) { + const ys = recurrenceGrid(p, v0, T, n); + let worst = 0; + for (let i = 0; i <= n; i++) { + const exact = sample((i / n) * T); + const err = Math.abs(ys[i]! - exact); + if (err > worst) worst = err; + } + // Бюджет реконструкции ≥ 1/400/2 = 1.25e-3; требуем запас ≥ 4 порядка. + expect(worst).toBeLessThan(1e-7); + } + }); + } +});