Finite-Time Dynamic Bounds

• Finite-Time Dynamic Bounds are a numerical approach that enforces physical invariants in time-dependent PDEs using variational inequality formulations and Bernstein bases.
• The methodology extends standard finite element and collocation methods to uniformly preserve solution bounds in both space and time.
• It employs a transformation linking traditional stage values to Bernstein coefficients, ensuring properties like nonnegative concentrations and bounded order parameters throughout simulations.

The enforcement of bounds constraints in finite element (FE) approximations of time-dependent partial differential equations (PDEs) is critical in applications where physical or mathematical invariants dictate that solutions must remain within prescribed intervals. Standard FE and collocation-based Runge–Kutta time discretizations generally do not preserve such invariants automatically. The referenced work develops a comprehensive methodology to enforce finite-time dynamic bounds by extending variational inequality formulations to high-order collocation-type time integrators and by employing a Bernstein basis reformulation for both the spatial FE and the time collocation polynomial. This permits bounds enforcement either at finite-dimensional degrees of freedom or, more stringently, uniformly throughout the spatial–temporal domain.

1. Motivation and Problem Setting

Many models—especially those in reaction–diffusion, phase separation, or ecological and biological systems—feature natural bounds on their state variables. For example, concentrations must be nonnegative, and order parameters in the Cahn–Hilliard equation lie within (–1, 1). Standard FE spatial discretizations, when combined with classical time-stepping schemes (explicit or implicit Runge–Kutta, collocation), cannot guarantee that the approximate solution lies within these bounds at all times, which may lead to nonphysical behaviors such as negative concentrations or divergence due to singular free energy terms.

To address this, the methodology in (Kirby et al., 20 Jun 2025) replaces the standard variational equality with a variational inequality (VI) at the fully discrete level and leverages specific FE bases (notably Bernstein polynomials) for both space and time, enforcing bounds at the level of basis coefficients. The approach is designed to yield formally high-order schemes that honor the bounds constraint uniformly in time.

2. Space–Time Discretization and Collocation Schemes

After semi-discretization in space (e.g., by Galerkin projection into a FE space $\mathcal{V}_h \subset \mathcal{V}$), the evolution equation $(y', v) = \langle F_t(y), v\rangle$ is discretized in time using an $s$-stage implicit Runge–Kutta method. The corresponding stage equations are: $$(Y_i, v_i) = (y^n, v_i) + k \sum_{j=1}^s A_{ij} \langle F_{t^n + c_jk}(Y_j), v_i \rangle, \qquad \forall v_i\in \mathcal{V}$$ where $Y = (Y_1, ..., Y_s)$ are the stage values, $A_{ij}$ and $c_j$ denote the Butcher tableau, and $k$ is the time step.

Such schemes are equivalent to collocation methods: a collocating polynomial $u(t)$ is sought on $[t^n, t^{n+1}]$ so that

$$u(t^n) = y^n, \qquad (u'(t^n + c_ik), v) = \langle F_{t^n + c_ik}(u(t^n + c_ik)), v\rangle.$$

A standard representation uses Lagrange interpolation at collocation points

$$u(t^n + \tau k) = \sum_{i=0}^s \bar{Y}_i \ell_i(\tau)$$

where $\ell_i(\tau)$ are Lagrange basis polynomials, and $\bar{Y}_i \in \mathcal{V}_h$.

3. Bounds Enforcement via Variational Inequalities

To ensure boundedness, the method replaces the variational equality with a VI: $$\langle \mathcal{F}(Y_h), V_h - Y_h \rangle \geq 0\quad \forall V_h \in \mathcal{K}_h^s,$$ where $\mathcal{K}_h$ is a convex set encoding the pointwise or coefficient-wise bounds (e.g., $[m,M]$-intervals in Bernstein bases). For spatial FE using Lagrange basis, constraints are enforced at the nodal degrees of freedom. For Bernstein basis, because of the convex hull property, if Bernstein coefficients of $p(x) = \sum_i c_i b_i(x)$ all lie in $I$, then $p(x)\in I$ for all $x$ on the element.

However, in time-collocation, enforcing constraints at the stage nodes (in Lagrange representation) does not guarantee that intermediate values respect the bounds due to possible overshoots. This motivates the Bernstein basis reformulation in time.

4. Bernstein Reformulation and Uniformly-In-Time Bounds

By representing the collocation polynomial in the time interval $[t^n, t^{n+1}]$ as

$u(t^n + k\tau) = \sum_{j=0}^s \bar{Z}_j b_j(\tau)$

where $b_j(\tau)$ are Bernstein polynomials and $\bar{Z}_0 = y^n$, one can enforce bounds directly on the $\bar{Z}_j$ coefficients. The mapping between traditional stage values $\bar{Y}_i$ and Bernstein coefficients is captured by a transformation (Bernstein–Vandermonde matrix).

Imposing bounds on the $\bar{Z}_j$ ensures that $u(t)$ stays within $[m,M]$ for all $\tau\in[0,1]$ (i.e., for the entire $[t^n, t^{n+1}]$ interval), achieving uniform-in-time bounds-preservation. This property is independent of the collocation nodes' placement and holds regardless of the method's stiff accuracy.

The complete bounds-constrained method is, at each step, to solve the variational inequality for the Bernstein coefficients $\bar{Z}$, using the transformation to map to stage values as needed for evaluating the collocation scheme.

5. Numerical Examples and Performance

(a) Phytoplankton Growth Model

A stiff nonlinear ODE model for phytoplankton growth, with known invariants (e.g., conservation of total carbon/nitrogen), is used. Without bounds enforcement, integrators can yield negative concentrations, causing subsequent blow-up. Both Lagrange and Bernstein-based VIs maintain nonnegativity, but only the Bernstein form ensures the solution is positive at all intermediate times. This is illustrated by plotting the collocation polynomial and its coefficients; in the Bernstein formulation, the minimum over time never breaches zero.

(b) Heat Equation

For parabolic problems (e.g., with Dirichlet or Neumann conditions), negative overshoots when simulating temperatures are avoided by enforcing bounds. The method demonstrates optimal convergence rates across polynomial degree and mesh refinement. Comparisons show that Bernstein-based methods, with bounds on coefficients, guarantee nonnegativity pointwise, not just at nodal points, and perform robustly with only a modest increase in computational cost (e.g., in Newton iterations).

(c) Cahn–Hilliard Equation

For the Flory–Huggins Cahn–Hilliard model, the order parameter must satisfy $c \in (-1, 1)$. Standard methods can overshoot, leading to failure due to nonphysical logarithms. Using quadratic Bernstein FE in space and Bernstein collocation in time, the full solution remains strictly within $(–1, 1)$, as shown by both solution snapshots and time series of the extremal coefficients. The scheme also preserves the qualitative behavior, such as monotonic free energy decay.

6. Practical Implications and Extensions

• The combination of FE discretization, collocation-type implicit Runge–Kutta schemes, and variational inequality enforcement using Bernstein bases in both space and time creates a framework for arbitrarily high-order, bounds-preserving time integration of PDEs with invariant regions.
• The approach is applicable not only to classical parabolic or elliptic PDEs but also to nonlinear systems with singularities or conservation laws.
• A critical distinction exists between bounds enforcement at discrete nodes (which does not guarantee uniform boundedness in between) and bounds enforcement via the Bernstein basis (ensuring global-in-time–and–space invariance).
• Newton-based solution of the resulting nonlinear VI systems converges efficiently in practice, and the additional cost of bounds-enforcement is moderate.

7. Summary Table: Basis Choice and Bounds Enforcement

Basis	Bounds Enforcement	Uniform-in-Time/Space?
Lagrange FE (space), Lagrange (time)	At nodal points only	No
Bernstein FE (space or time)	Coefficient constraints	Yes

8. Conclusion

The bounds-constrained FE approach via variational inequalities and Bernstein representation yields formally high-order, stable, and physically meaningful numerical solutions for time-dependent PDEs. By enforcing bounds uniformly through the basis coefficients, this method is especially suited for problems with strict invariants, offering a general, computationally tractable, and easily extensible framework for preserving finite-time dynamic bounds in a variety of applications (Kirby et al., 20 Jun 2025).