Theory¶
This document proves every mathematical fact the Twindex checker relies on. Each VALID verdict is a consequence of the results below plus exact rational arithmetic; no floating-point computation is involved in checking.
Notation: I = [0, 1], ‖f‖ = max_{x∈I} |f(x)|.
1. The problem¶
A model fixes a coefficient degree d ≥ 0, a polynomial source g,
boundary values α, β (rational; default 0), and a nonlinearity
f(u) = Σ_{m ∈ M} a_m u^m, every m ∈ M odd, every a_m ≥ 0.
For θ ∈ ℝ^{d+1} let
q(x; θ) = Σ_{j=0}^{d} θ_j B_j(x), B_j(x) = C(d, j) x^j (1 − x)^{d−j},
and consider
(P_θ) −u''(x) + q(x; θ) f(u(x)) = g(x) on (0, 1), u(0) = α, u(1) = β.
The original (v1) model is d = 3, f(u) = u³, g = 1, α = β = 0.
Two facts are used throughout.
- The
B_jare nonnegative onIand sum to 1. Hence θ ≥ 0 (componentwise) implies q ≥ 0 on I. Every prior box Twindex accepts haslo_j ≥ 0. - For odd
m, defineh_m(a, b) = Σ_{i=0}^{m−1} a^i b^{m−1−i}, so thata^m − b^m = (a − b) h_m(a, b). Becauset ↦ t^mis nondecreasing for oddm,h_m(a, b) ≥ 0whenevera ≠ b; by continuity it is also≥ 0whena = b. SettingH(a, b) = Σ_m a_m h_m(a, b) ≥ 0gives f(a) − f(b) = (a − b) H(a, b) with H ≥ 0 continuous.
2. A one-dimensional comparison lemma¶
Lemma 1. Let c ∈ C(I) with c ≥ 0, and let w ∈ C²(0,1) ∩ C(I) satisfy
−w'' + c w ≥ 0 on (0,1) and w(0) ≥ 0, w(1) ≥ 0. Then w ≥ 0 on I.
Proof. Suppose w(x₀) < 0 for some x₀. Let (a, b) be the connected
component of the open set {x ∈ I : w(x) < 0} containing x₀. At an endpoint
inside (0,1) we have w = 0 by continuity. At an endpoint equal to 0 or 1 we
have w ≥ 0 by hypothesis and w ≤ 0 as a limit of negative values, so again
w = 0. On (a, b), c w ≤ 0, so −w'' ≥ −c w ≥ 0: w is concave on [a, b].
A concave function vanishing at both ends of an interval is ≥ 0 on it,
contradicting w(x₀) < 0. ∎
3. Existence and uniqueness¶
Theorem 1. If q ∈ C(I) and q ≥ 0, then (P_θ) has exactly one solution
u ∈ C²(I).
Proof.
Existence. Let F(t) = Σ_m a_m t^{m+1}/(m+1). Each m + 1 is even, so F is
convex and F ≥ 0. Let ℓ(x) = α(1 − x) + βx and write u = ℓ + v with
v ∈ H¹₀(0,1). Consider
E(v) = ∫₀¹ ( ½ (v' + ℓ')² + q F(ℓ + v) − g (ℓ + v) ) dx. In one dimension
H¹₀ ⊂ C(I) with ‖v‖_{L²} ≤ ‖v‖ ≤ ‖v'‖_{L²}, so E is well defined and
E(v) ≥ ½‖v'‖²_{L²} − (|ℓ'| + ‖g‖_{L²})‖v'‖_{L²} − C for a constant C, which
is coercive. E is strictly convex
(strictly convex ½‖v'‖², plus the convex ∫ q F(v) since q ≥ 0, plus a
linear term). It is weakly lower semicontinuous: the first term is a norm
squared, and the other two are continuous under the compact embedding
H¹₀ ⊂ C(I). So E has a unique minimiser u, whose Euler–Lagrange equation is
∫ u' φ' + q f(u) φ − g φ = 0 for all φ ∈ H¹₀, with u = ℓ + v. Since q f(u) − g ∈ C(I), the
weak second derivative u'' = q f(u) − g is continuous, and u ∈ C²(I) solves
(P_θ) classically, with u(0) = α, u(1) = β.
Uniqueness. If u₁, u₂ are classical solutions, e = u₁ − u₂ satisfies
−e'' + c e = 0 with c = q H(u₁, u₂) ≥ 0 (§1), and e(0) = e(1) = 0.
Lemma 1 applied to e and to −e gives e = 0. ∎
Special case (v1 model, f(u) = u³, g = 1). Here also
0 ≤ u ≤ x(1 − x)/2. Lower bound: at a negative interior minimum x₀,
−u''(x₀) + q u(x₀)³ ≤ 0 < 1, a contradiction. Upper bound: v = x(1−x)/2 − u
satisfies −v'' = q u³ ≥ 0 with zero boundary values, so v ≥ 0 by Lemma 1.
The checker does not use these bounds.
4. The barrier bound (from a candidate to the true solution)¶
Theorem 2. Let q ≥ 0 on I and let u solve (P_θ). Let p be a
polynomial with p(0) = α and p(1) = β, set r = −p'' + q f(p) − g, and suppose
‖r‖ ≤ ρ. Then for every x ∈ I
|u(x) − p(x)| ≤ ρ · x(1 − x) / 2.
Proof. Put e = u − p and c = q H(u, p). By §1, c ≥ 0, and it is
continuous. Subtracting the equations, and using f(u) − f(p) = (u − p) H(u, p),
−e'' + c e = (−u'' + q f(u)) − (−p'' + q f(p)) = g − (r + g) = −r.
Let b = ρ x(1−x)/2, so −b'' = ρ and b ≥ 0. For w± = b ∓ e,
−w±'' + c w± = ρ + c b ± r ≥ ρ − |r| ≥ 0,
and w±(0) = w±(1) = 0, because e = u − p vanishes at both ends. Lemma 1
gives w± ≥ 0, i.e. |e| ≤ b. ∎
Notes.
* The bound needs no assumption about how p was obtained: no convergence
of a numerical method, no spectral tail estimate and no inverse-operator norm.
It also holds pointwise; the checker uses the pointwise value x_i(1−x_i)/2 at
each sensor, which is never larger than 1/8.
* The hypotheses the checker must verify are exactly: p(0) = 0, p(1) = 0
(exact rational equality), θ ≥ 0 (via the prior box), and ‖r‖ ≤ ρ (§6).
With inhomogeneous boundary values the first two become p(0) = α, p(1) = β.
5. Uniform bounds over a parameter box¶
Theorem 3. Fix p. For θ, c ∈ ℝ^{d+1},
r(x; θ) = r(x; c) + Σ_j (θ_j − c_j) B_j(x) f(p(x)). Consequently, if
|θ_j − c_j| ≤ h_j for all j,
‖r(·; θ)‖ ≤ ‖r(·; c)‖ + Σ_j h_j ‖B_j f(p)‖.
Proof. q is linear in θ and r depends on θ only through q f(p). The
inequality is the triangle inequality. ∎
So one candidate p and one exact computation give a residual bound ρ_B valid
for every θ in a box B = c ± h. If B lies in a prior box with lo ≥ 0,
Theorem 2 then applies to every θ ∈ B with the same ρ_B.
Theorem 3′ (affine candidate families). Fix polynomials p₀, s₀, …, s_d
with p₀(0) = α, p₀(1) = β and every s_j(0) = s_j(1) = 0, a centre c and
half-widths h ≥ 0. For Δ ∈ D = Π_j [−h_j, h_j] put
p_Δ = p₀ + Σ_j Δ_j s_j, θ = c + Δ.
The residual r(x; Δ) = −p_Δ'' + q(x; c + Δ) f(p_Δ) − g is a polynomial in
(x, Δ). Write it as Σ_α Δ^α R_α(x), a finite sum over multi-indices α.
Then:
sup_{x∈I, Δ∈D} |r(x; Δ)| ≤ ρ_D := Σ_α h^α ‖R_α‖, whereh^α = Π_j h_j^{α_j};- if
c + D ⊆ [0, ∞)^{d+1}, then for everyθ = c + ΔwithΔ ∈ D,|u_θ(x) − p_Δ(x)| ≤ ρ_D x(1−x)/2; - at each
x,{p_Δ(x) : Δ ∈ D} = [p₀(x) − Σ_j h_j |s_j(x)|, p₀(x) + Σ_j h_j |s_j(x)|].
Proof. (1) |Δ^α| ≤ h^α on D, then the triangle inequality. (2) For a fixed
Δ, p_Δ is a polynomial with p_Δ(0) = α, p_Δ(1) = β, and its residual is bounded by
ρ_D by (1). Since θ ≥ 0, q ≥ 0, so Theorem 2 applies to the candidate p_Δ.
(3) p_Δ(x) is affine in Δ, so over the box it attains exactly the interval
p₀(x) ± Σ_j h_j |s_j(x)|. ∎
With all s_j = 0 this is exactly Theorem 3. The gain comes from choosing s_j
close to the sensitivities ∂u/∂θ_j: the first-order terms R_{e_j} then
nearly cancel, leaving a residual that is O(|h|²). The checker expands r
exactly (rational arithmetic, at most 20 000 monomials in Δ) and computes every
‖R_α‖ with the certified bound of §6.
Sensor tests over the box (§7) use the exact interval from (3). Every θ ∈ c + D
is compatible if max(|lo_i − y_i|, |hi_i − y_i|) + s_i ≤ ε at every sensor.
No θ is compatible if lo_i − y_i > ε + s_i or y_i − hi_i > ε + s_i at some
sensor. Here s_i = ρ_D x_i(1−x_i)/2.
5.1 Monotone two-corner boxes¶
Lemma 2 (positivity). If q ≥ 0, g ≥ 0 on I and α, β ≥ 0, then the
solution of (P_θ) satisfies u ≥ 0.
Proof. Write f(u) = u φ(u) with φ(u) = Σ_m a_m u^{m−1}. Every m − 1 is
even and every a_m ≥ 0, so φ ≥ 0, and φ(u(x)) is continuous. Then
−u'' + c u = g ≥ 0 with c = q φ(u) ≥ 0, and u(0) = α ≥ 0, u(1) = β ≥ 0.
Lemma 1 gives u ≥ 0. ∎
Theorem 6 (monotone dependence on θ). Under the hypotheses of Lemma 2, if
0 ≤ θ ≤ θ' componentwise, then u_θ ≥ u_θ' on I.
Proof. Let w = u_θ − u_θ' and write q = q(·; θ), q' = q(·; θ').
Subtracting the two equations,
−w'' + q (f(u_θ) − f(u_θ')) = (q' − q) f(u_θ').
By §1, f(u_θ) − f(u_θ') = H(u_θ, u_θ') w with H ≥ 0, so the left side is
−w'' + c w with c = q H ≥ 0. On the right, q' − q = Σ_j (θ'_j − θ_j) B_j ≥ 0,
and f(u_θ') ≥ 0 because u_θ' ≥ 0 (Lemma 2) and f is nondecreasing with
f(0) = 0. Since w(0) = w(1) = 0, Lemma 1 gives w ≥ 0. ∎
Corollary 7 (two-corner box). Let B = Π_j [lo_j, hi_j] ⊆ [0, ∞)^{d+1} and
assume the hypotheses of Lemma 2. Let p_lo and p_hi be candidates for
θ = lo and θ = hi with residual bounds ρ_lo, ρ_hi (Theorem 2). Then for
every θ ∈ B and every x ∈ I,
p_hi(x) − ρ_hi x(1−x)/2 ≤ u_hi(x) ≤ u_θ(x) ≤ u_lo(x) ≤ p_lo(x) + ρ_lo x(1−x)/2.
So B ⊆ C_y if every sensor interval lies in [y_i − ε, y_i + ε], and
B ∩ C_y = ∅ if some sensor interval misses it entirely.
Proof. lo ≤ θ ≤ hi, so Theorem 6 gives the two middle inequalities. Theorem 2
gives the outer ones. ∎
Corollary 7 needs only two ordinary residual bounds per box, independent of the
box width. Its sensor intervals enclose the true range [u_hi(x), u_lo(x)] up to
the two corner errors, so they are nearly sharp. The checker certifies g ≥ 0
exactly: every Bernstein coefficient of every piece of a dyadic subdivision is
≥ 0, which implies g ≥ 0 (§6, Fact A). It reports INVALID if some attained
value is negative, and INCONCLUSIVE if neither is decided.
6. Certified sup norms of polynomials¶
Let p have degree at most n, and write it in the Bernstein basis of degree
n on I:
p(x) = Σ_{i=0}^{n} b_i C(n, i) x^i (1 − x)^{n−i},
b_i = Σ_{k=0}^{i} [C(i, k) / C(n, k)] a_k (a_k = monomial coefficients).
Fact A (enclosure). The basis functions are nonnegative on I and sum to 1,
so p(x) is a convex combination of the b_i and |p(x)| ≤ max_i |b_i|.
Fact B (endpoints). b_0 = p(0) and b_n = p(1).
Fact C (subdivision). The de Casteljau algorithm at t = 1/2 maps the
coefficients b of p on [α, β] to the exact Bernstein coefficients of p
on [α, m] and [m, β], where m = (α + β)/2. Each step averages adjacent
entries; the implementation keeps integer numerators and multiplies the common
denominator by 2ⁿ, so the computation is exact.
Algorithm (sup_abs_bound). Start from the piece I. Keep a best attained
value L: the maximum of |b_0| and |b_n| over every piece generated so far,
each of which is |p| at a point by Fact B. Repeatedly take the piece with the
largest enclosure U = max |b_i|. If U ≤ L, stop: by Fact A every remaining
piece satisfies |p| ≤ U ≤ L on it. Otherwise, if the piece has reached the
depth limit, record U; if not, split it by Fact C. Return the maximum of L
and all recorded U.
Theorem 4. The returned value R satisfies R ≥ ‖p‖.
Proof. The pieces form a partition of I into dyadic intervals; every point
lies in some piece that was either recorded (bound U ≤ R), discarded when its
enclosure fell to U ≤ L ≤ R, or never split because the loop stopped (its
enclosure is at most the popped U, which is ≤ L ≤ R). In each case Fact A
bounds |p| on that piece by a quantity ≤ R. ∎
L is also a lower bound on ‖p‖, so R − L measures how tight the result is.
7. From sensor inequalities to the claims¶
Fix a problem instance (sensors x_i, rational data y_i, tolerance ε,
prior box Θ ⊆ [0, ∞)^{d+1}, target index k), and let
C_y = {θ ∈ Θ : max_i |u_θ(x_i) − y_i| ≤ ε}.
Suppose Theorem 2 (or Theorem 3 with Theorem 2) gives |u_θ(x_i) − p(x_i)| ≤ s_i.
By the triangle inequality:
- if
|p(x_i) − y_i| + s_i ≤ εfor everyi, thenθ ∈ C_y(certified compatible); - if
|p(x_i) − y_i| − s_i > εfor somei, thenθ ∉ C_y(certified incompatible); - otherwise these sufficient conditions decide nothing (inconclusive).
For a box B ⊆ Θ with a uniform bound, the same three cases certify
B ⊆ C_y, B ∩ C_y = ∅, or nothing.
Corollary (pair certificate). If θ, θ' ∈ C_y are certified compatible,
then D_Q(y) := sup_{φ,φ'∈C_y} |φ_k − φ'_k| ≥ |θ_k − θ'_k|. Moreover every
estimator T (any function of the data) satisfies
max( |T(y) − θ_k|, |T(y) − θ'_k| ) ≥ |θ_k − θ'_k| / 2,
because both θ and θ' are consistent with the same observed y at noise
level ε. This is the triangle inequality. So no method, PINN or otherwise, can
recover θ_k from these data with worst-case error below half the certified gap.
Theorem 5 (paving certificate). Let a binary tree be given whose root box is
Θ and in which every internal node splits its box B on an axis j at a point
m with lo_j(B) < m < hi_j(B) into the two closed boxes
{θ ∈ B : θ_j ≤ m} and {θ ∈ B : θ_j ≥ m}. Suppose every leaf labelled
inside is certified ⊆ C_y, and every leaf labelled outside is certified
∩ C_y = ∅, by §5 and the case analysis above. Let In be the inside leaves and
Out⁺ the inside and unresolved leaves. Then
max_{B∈In} hi_k(B) − min_{B∈In} lo_k(B) ≤ D_Q(y) ≤ max_{B∈Out⁺} hi_k(B) − min_{B∈Out⁺} lo_k(B),
where the left side requires In ≠ ∅. If Out⁺ = ∅ then C_y = ∅.
Proof. Coverage: by induction on the tree, the leaves cover Θ, because each
split writes its box as the union of its two children. Upper bound: let
θ ∈ C_y ⊆ Θ. Then θ lies in some leaf, and that leaf is not an outside leaf
(those are disjoint from C_y), so it belongs to Out⁺ and
θ_k ∈ [min_{Out⁺} lo_k, max_{Out⁺} hi_k]. Hence any two points of C_y differ
in coordinate k by at most the stated width. If Out⁺ = ∅ no leaf can contain
θ, so C_y = ∅. Lower bound: pick an inside leaf A attaining the maximum
of hi_k and an inside leaf A' attaining the minimum of lo_k. Then choose
θ ∈ A with θ_k = hi_k(A) and θ' ∈ A' with θ'_k = lo_k(A'). Both lie in
C_y because inside leaves are subsets of C_y. ∎
The checker verifies the split condition at every node, recomputes every leaf box exactly from the root and the splits, re-certifies every inside and outside leaf, and compares the claimed bounds with the two expressions above.
Theorem 8 (oriented paving). Fix a rational origin o, an invertible
rational matrix A and a root box R in coordinates φ, and set
θ = o + A φ. Suppose R contains the interval hull of A⁻¹(Θ − o), and take a
bisection tree over R whose leaves are φ-boxes P. Write θ(P) for the
parallelotope o + A P and H(P) for its exact interval hull. Suppose:
- every leaf labelled outside either satisfies
H(P) ∩ Θ = ∅, or carries a candidate familyp_Δ = p₀ + Σ_k Δ_k s_k(withθ = o + A(c + Δ),cthe centre ofP,Δ ∈ P − c) that the sensor tests of Theorem 3′ prove incompatible; - every leaf labelled inside satisfies
H(P) ⊆ Θand is proved compatible in the same way.
Let In be the inside leaves and Out⁺ the inside and unresolved leaves with
H(P) ∩ Θ ≠ ∅. Then
max_{In} hi_k(H(P)) − min_{In} lo_k(H(P)) ≤ D_Q(y) ≤
max_{Out⁺} min(hi_k(H(P)), hi_k(Θ)) − min_{Out⁺} max(lo_k(H(P)), lo_k(Θ)).
Proof. Every θ ∈ Θ has φ = A⁻¹(θ − o) ∈ R, hence lies in θ(P) for some
leaf P. A leaf with H(P) ∩ Θ = ∅ contains no point of Θ. For θ ∈ θ(P) ∩ Θ
we have θ ≥ 0, so q ≥ 0. Since q(c_θ + AΔ) = q(c_θ) + Σ_k Δ_k E_k with
E_k = Σ_j A_{jk} B_j, the expansion of Theorem 3′ with directions E_k bounds
the residual of the family member at Δ = φ − c, and Theorem 2 applies. Hence
outside leaves contain no point of C_y, and inside leaves (which lie in Θ)
contain only points of C_y. The upper bound follows as in Theorem 5, using
H(P) ∩ Θ as an enclosure of each remaining leaf. For the lower bound, θ_k over
the parallelotope θ(P) attains both ends of H(P)_k at vertices. ∎
Oriented pavings let the leaves follow a thin, slanted compatible set, for
example by choosing the columns of A along the singular directions of the
sensitivity matrix. Monotone corners (Corollary 7) are not used in oriented
pavings, because the componentwise order is not preserved by A.
7.1 Two space dimensions (unit square)¶
Let Ω = (0, 1)², q(x, y; θ) = Σ_{j,k} θ_{jk} B_j(x) B_k(y) with θ ≥ 0 (so
q ≥ 0), f as in §1, g a polynomial, and consider
(P²_θ) −Δu + q f(u) = g in Ω, u = 0 on ∂Ω.
Well-posedness. E(v) = ∫_Ω ½|∇v|² + q F(v) − g v is coercive, strictly convex
and weakly lower semicontinuous on H¹₀(Ω). Since H¹₀ ⊂ L^r for every r < ∞ in
two dimensions, ∫ q F(v) is finite. The unique minimiser is the unique weak
solution u. Monotone truncation (Stampacchia) gives u ∈ L^∞, and interior
elliptic regularity makes u continuous, in fact smooth, inside Ω.
Theorem 9 (2-D barrier). Let p be a polynomial vanishing on ∂Ω with
|−Δp + q f(p) − g| ≤ ρ on Ω̄. Let w ∈ C²(Ω̄) satisfy −Δw ≥ 1 in Ω and
w ≥ 0 on ∂Ω. Then |u − p| ≤ ρ w in Ω.
Proof. e = u − p ∈ H¹₀(Ω) is a weak solution of −Δe + c e = −r with
c = q H(u, p) ∈ L^∞, c ≥ 0 (§1). The functions v± = ρw ∓ e satisfy
−Δv± + c v± ≥ ρ + cρw ± r ≥ 0 weakly, because ρw ≥ 0 on Ω̄ (by the weak
maximum principle, −Δw ≥ 0 and w ≥ 0 on ∂Ω give w ≥ 0), and v± ≥ 0 on
∂Ω. The weak maximum principle for −Δ + c with c ≥ 0 (Gilbarg–Trudinger,
Theorem 8.1) gives v± ≥ 0 a.e. Since u and p are continuous in Ω, the
bound holds pointwise. ∎
The checker uses, pointwise at each sensor, the minimum of the three barriers
x(1−x)/2, y(1−y)/2 and (x(1−x) + y(1−y))/4. Each satisfies −Δw = 1 and
vanishes or is nonnegative on ∂Ω, so each bound holds separately. The residual
bound is certified by the tensor-product Bernstein enclosure (twindex/poly2.py),
whose Facts A–C hold coordinatewise. Corner compatibility: a
polynomial vanishing on ∂Ω has Δp = 0 at the corners, so its residual there is
−g(corner). If g is nonzero at a corner, the solution is not C² up to that
corner and ρ ≥ |g(corner)|. This floor can block certification at small ε,
though a large enough ε can still close. For small residuals, use sources that
vanish at the corners. "p vanishes on ∂Ω" is checked exactly on
the coefficient grid: the first row and column vanish, and the row and column
sums vanish.
7.2 Time-dependent problems (heat equation)¶
Let Q = (0, 1) × (0, 1] (time scaled to [0, 1]), q(x; θ) = Σ_j θ_j B_j(x) with
θ ≥ 0, f as in §1, g(x, t) and u₀(x) polynomials with u₀(0) = u₀(1) = 0, and
consider
(H_θ) u_t − u_xx + q f(u) = g in Q, u(0, t) = u(1, t) = 0, u(x, 0) = u₀(x).
Well-posedness. f is nondecreasing with f(0) = 0 and q ≥ 0, so the
comparison principle gives uniqueness and an a-priori L^∞ bound. Standard theory
for semilinear parabolic equations with monotone lower-order terms gives a solution
that is C^{2,1} in Q and continuous on Q̄. With the compatibility below it is
smooth up to the corners to the order used. See Lieberman, Second Order Parabolic
Differential Equations, Ch. 2 and 5.
Theorem 10 (parabolic barrier). Let p(x, t) be a polynomial that matches the
data exactly: p(0, t) = p(1, t) = 0 and p(x, 0) = u₀(x). Suppose
|p_t − p_xx + q f(p) − g| ≤ ρ on Q̄. Then
|u − p| ≤ ρ · min(t, x(1 − x)/2) in Q.
Proof. e = u − p vanishes on the parabolic boundary. It satisfies
e_t − e_xx + c e = −r, where r is the residual of p and
c = q H(u, p) ≥ 0 is bounded (§1). Let w be either t or x(1 − x)/2. Both
satisfy w_t − w_xx = 1 and w ≥ 0 on Q̄. The functions v± = ρw ∓ e then
satisfy v±_t − v±_xx + c v± = ρ(1 + c w) ± r ≥ 0 and v± ≥ 0 on the parabolic
boundary. The weak maximum principle for ∂_t − ∂_xx + c with bounded c ≥ 0
gives v± ≥ 0, so |e| ≤ ρw. Taking the smaller barrier pointwise finishes the
proof. ∎
The checker verifies the data conditions exactly on the coefficient grid
c[i][j] of x^i t^j:
- row 0 vanishes (
p(0, t) = 0); - every column sum vanishes (
p(1, t) = 0); - column 0 equals the coefficients of
u₀(p(x, 0) = u₀).
The residual bound uses the tensor Bernstein enclosure of §7.1. Sensors must have
0 < x < 1 and 0 < t ≤ 1.
Corner compatibility. At a corner (c, 0) with c ∈ {0, 1}, p(c, t) = 0
forces p_t(c, 0) = 0, so the residual there is −u₀″(c) − g(c, 0). Any data with
u₀″(c) + g(c, 0) ≠ 0 therefore force ρ to be at least that value. Polynomial
candidates converge quickly only when the next order also holds. For g = 0,
differentiating the equation in time at the lateral boundary gives u_xxt = 0.
Consequently the equation differentiated twice in space requires
u₀⁗(c) = (q f(u₀))″(c). When f′(0) = 0, this reduces to u₀⁗(c) = 0
(using u₀(c) = 0 and the odd-power form of f). A linear reaction term can
make the right-hand side nonzero, so vanishing fourth derivatives alone need
not suffice. The default cubic nonlinearity and
u₀ = x − (5/3)x³ + x⁵ − (1/3)x⁶ satisfy both orders. In this setting,
space-time Chebyshev candidates with n = 16 reach ρ ≈ 6.5·10⁻³, and
n = 24 reaches ρ ≈ 4·10⁻⁴. The certified twins in tests/test_parabolic.py
are θ₀ = 5 and θ₀ = 10 at noise level ε = 10⁻³.
7.3 Polynomial ODE systems¶
Let F(t, u; θ) be a polynomial vector field on ℝ × ℝⁿ and u₀(θ) a polynomial
initial value. Signs are unrestricted. Consider
(O_θ) u′ = F(t, u; θ), u(0) = u₀(θ), t ∈ [0, T],
with observables y = c·u + c₀ (affine in the state) measured at times tᵢ. F is
smooth, so (O_θ) has a unique maximal solution, which may blow up. Write |·| for
the max norm on ℝⁿ. Write μ(A) = maxᵢ (Aᵢᵢ + Σ_{j≠i} |Aᵢⱼ|) for its logarithmic
norm, which is subadditive and positively homogeneous, hence convex.
A candidate consists of:
- breaks
0 = b₀ < … < b_K ≥ T; - on each segment, polynomials
p_k(τ)inτ = (t − b_k)/h_k; - a tube radius
δ > 0.
On segment k, let ρ_k ≥ sup |p_k′ − F(t, p_k)|, and let
μ_k ≥ sup { μ(∂F/∂u(t, v)) : t ∈ segment, |v − p_k(t)| ≤ δ }.
Theorem 11 (tube validation). Let E₀ = |p₀(0) − u₀|. For k ≥ 1, let
E_k = Ē_{k−1} + |p_k(0) − p_{k−1}(1)|, where Ē_{k−1} is the end-of-segment
bound below. Suppose that on every segment μ_k⁺ h_k < 1, where x⁺ = max(x, 0),
and
(E_k + ρ_k h_k) / (1 − μ_k⁺ h_k) < δ.
Then the solution exists on [0, b_K], and for t = b_k + s in segment k
|u(t) − p_k(t)| ≤ E_k / (1 − μ_k s) + ρ_k s / (1 − μ_k⁺ s) =: B_k(s),
with Ē_k = B_k(h_k).
Proof. Fix segment k and suppose |e(b_k)| ≤ E_k with e = u − p_k. As long as
|e| ≤ δ, the segment from p_k to u lies in the tube. Hence
e′ = A(t) e − r with A = ∫₀¹ ∂F/∂u(t, p_k + σe) dσ, and by convexity
μ(A) ≤ μ_k. The upper Dini derivative then satisfies
D⁺|e| ≤ μ(A)|e| + |r| ≤ μ_k|e| + ρ_k (Dahlquist 1958; see Söderlind, BIT 46 (2006) 631–652). Comparison
with the scalar equation gives
|e(b_k + s)| ≤ E_k e^{μ_k s} + ρ_k (e^{μ_k s} − 1)/μ_k.
Two inequalities bound the right-hand side:
e^x ≤ 1/(1 − x)for everyx < 1, including negativex.(e^{μs} − 1)/μ = ∫₀ˢ e^{μσ} dσ ≤ s · e^{μ⁺s} ≤ s/(1 − μ⁺s).
Together they give B_k.
The supremum of B_k over [0, h_k] is at most (E_k + ρ_k h_k)/(1 − μ_k⁺ h_k) < δ.
So the first time |e| reached δ would contradict the bound. The solution therefore
stays in the tube, which also rules out blow-up, and the estimate holds on the whole
segment. Induction over the segments, together with the triangle inequality at each
jump, finishes the proof. ∎
Weighted norms. Given exact scales w_i > 0, apply the theorem to
z_i = u_i / w_i. That system has right-hand side F_i / w_i and Jacobian
J_ij w_j / w_i. In practice:
- the residual and the initial or jump errors are divided by
w_i; - the tube becomes
|v_i − p_i| ≤ δ w_i; μ_w = max_i (J_ii + Σ_{j≠i} |J_ij| w_j / w_i);- an observable
c·upicks upΣ_j |c_j| w_j B.
Nothing else changes. Scaling matters for models whose states carry very
different units (tests/test_ode.py::test_state_scales_rescue_badly_scaled_units).
Corollary. At a measurement (tᵢ, y = c·u + c₀, yᵢ) the checker compares
|c·p(tᵢ) + c₀ − yᵢ| ± ‖c‖₁ B(tᵢ) with ε. The verdict rules are those of §7: a
witness is VALID when every sum is ≤ ε, and INVALID when some difference
exceeds ε. INVALID is sound because (O_θ) has a unique solution, which the
theorem places within B of p.
Certified enclosures. Both ρ_k and the zero-order diagonal term of μ_k are
univariate polynomials in τ, bounded by Bernstein enclosures with de Casteljau
refinement (§6). The diagonal term uses a signed upper bound, so dissipative
systems keep μ_k < 0.
Every other contribution is bounded through the tube expansion. Write
v = p_k(τ) + δs with |s| ≤ 1 and expand
∂F/∂u(t, v) = Σ_α c_α(τ) s^α. Then |s^α| ≤ 1 gives
Jᵢᵢ + Σ_{j≠i} |Jᵢⱼ| ≤ sup c⁽ⁱⁱ⁾₀ + Σ_{α≠0} sup |c⁽ⁱⁱ⁾_α| + Σ_{j≠i} Σ_α sup |c⁽ⁱʲ⁾_α|.
Model strings are parsed by a closed grammar (twindex/mpoly.py): numbers, names,
+ − *, division by nonzero constants, and integer powers. Decimal literals are
read exactly, so "0.1" means 1/10.
Example (tests/test_ode.py). Consider SIR (S′ = −βSI, I′ = βSI − γI,
S₀ = 0.99, I₀ = 0.01), with I observed at t = 0.5, 1, …, 4 and ε = 10⁻³.
Twindex certifies both (β, γ) = (0.5, 0.2) and (0.72, 0.4178) as data-compatible.
These share the early growth rate βS₀ − γ, yet their reproduction numbers
β/γ are 2.5 and about 1.72. Early-epidemic data of this kind provably cannot
determine R₀.
7.4 Elementary right-hand sides (certified polynomialization)¶
Right-hand sides may use + − *, division by any expression, integer powers
(including negative ones), exp, sin and cos. This covers Michaelis–Menten
and Hill kinetics, Monod growth, Holling responses, frequency-dependent
incidence βSI/N, Arrhenius rates exp(−E/T) and periodic forcing.
Construction. Atoms are processed innermost first. Each non-polynomial atom
a whose argument depends on time or on a state becomes an auxiliary state.
Write ḋ = ∂_t d + Σ_x ∂_x d · F̃_x for the derivative of an argument d along
the lifted vector field F̃; the sum runs over every lifted state x. Then:
| atom | auxiliary state | equation | initial value |
|---|---|---|---|
1/D |
z |
z′ = −z² Ḋ |
1/D(0) (exact rational) |
exp(g) |
z |
z′ = z ġ |
exp(g(0)) (exact rational interval) |
sin(g), cos(g) |
s, c |
s′ = c ġ, c′ = −s ġ |
sin g(0), cos g(0) (intervals) |
A reciprocal 1/D whose argument involves only parameters is a derived
constant 1/D(θ). It is computed exactly per θ, and its value is INVALID if
D(θ) = 0. The lifted right-hand side F̃ is a polynomial.
Theorem 12. Let v = (u, z) solve the lifted system on [0, T] with the
exact (possibly irrational) initial values above. Then u solves the original
system on [0, T], and every denominator D(u(t)) is nonzero there.
Proof. Take the atoms in creation order. Each atom's argument involves only states and earlier atoms, which are already known to equal their intended expressions.
- Reciprocal. Let
w = D z. Thenw′ = Ḋ z + D(−z² Ḋ) = Ḋ z (1 − w), andw(0) = 1. This is a linear equation inwwith continuous coefficient, sow ≡ 1. Hencez = 1/DandD ≠ 0. - Exponential.
w = z e^{−g}hasw′ = z ġ e^{−g} − z e^{−g} ġ = 0, sow ≡ 1. - Sine and cosine. Let
φ = g(u(t)), so thatφ′ = ġ. Both(s, c)and(sin φ, cos φ)solve the linear systems′ = c ġ,c′ = −s ġwith the same initial value. By uniqueness they coincide, and the phase is pinned, not only the circles² + c² = 1.
So each atom equals its intended expression, F̃(v) = F(u), and u′ = F(u). ∎
Certification. Theorem 11 applies to the lifted polynomial system with one
change: an initial value given only as an interval [lo, hi] contributes
max(|p(0) − lo|, |p(0) − hi|) to E₀. This bounds |p(0) − v(0)| for the true
initial value whatever it is. Nothing is chosen inside the interval, so a
candidate that assumes a wrong value of, say, exp(1) pays for the difference.
The enclosures of exp, sin and cos at a rational point x with |x| ≤ 64
are exact Taylor sums with Lagrange remainder |x|^k/k!, multiplied by
3^{⌈x⌉} for exp with x > 0. They are rounded outward to multiples of 2⁻¹⁶⁰.
Initial arguments that involve earlier atoms are evaluated in exact interval
arithmetic.
The checker derives the lifting itself, deterministically, from the original
model strings in the certificate. The witness supplies candidates for the
auxiliary states as well. Observables refer to the user's states only. A model
whose true solution reaches a singularity, for example u′ = 1/(1 − u) at
t = 1/2, can never be certified past it, because Theorem 11 would prove that
the lifted solution exists there.
Prior work. Polynomialization of elementary ODEs is classical; see Hemery, Fages and Soliman (2021) and the quadratization literature. What is new here is that the lifting runs inside an exact, replayable certificate checker.
8. What a certificate does not say¶
- It concerns the declared data
yand the bounded-noise modelmax_i |noise_i| ≤ ε. It is not a statistical statement and does not assert non-identifiability for exact (noise-free) data. - A pair certificate gives a lower bound on
D_Q(y). Upper bounds need a covering ofΘby a paving certificate (Theorem 5), and a paving is only as tight as its unresolved leaves allow. - INCONCLUSIVE means the sufficient inequalities did not close. It is not evidence that the claim is false.
9. Related work¶
The ingredients are classical; Twindex's contribution is putting them together in an exact, independently replayable checker.
- Optimal recovery and the radius of information: Micchelli & Rivlin (1977); Foucart & Liao, arXiv 2111.02601 and arXiv 2401.11112.
- Optimal recovery for parametric PDEs: Binev et al., SIAM/ASA JUQ 5 (2017); Cohen, Dahmen, Mula & Nichols, SIAM/ASA JUQ 10 (2022), arXiv 2009.02687.
- Guaranteed (set-membership) parameter estimation and set inversion: Milanese & Vicino (1991); Jaulin & Walter, Automatica 29 (1993); Raïssi, Ramdani & Candau, Automatica (2004).
- Certified reduced-basis possibility regions: Sen et al., J. Comput. Phys. (2006); Nguyen, Rozza, Huynh & Patera (2011).
- Maximum-principle a posteriori analysis for semilinear problems: Nochetto, Schmidt, Siebert & Veeser.
- Computer-assisted proofs for PDEs: Nakao, Plum & Watanabe (2019).
- Error certification for physics-informed networks: Eiras et al., ICML 2024.