Convergent Operator-Splitting Scheme for Viscosity Solutions: A Foundation for Learning Domain-to-Solution Maps

This work introduces and rigorously analyzes a novel operator-splitting finite element scheme for approximating viscosity solutions of a broad class of constrained second-order partial differential equations. By decoupling the primary PDE evolution from the enforcement of constraints, the proposed method combines a stabilized finite element method for spatial discretization with an efficient semi-implicit time-stepping strategy. The cornerstone of our analysis is a proof that the scheme satisfies a discrete comparison principle. We demonstrate that under a mild time-step restriction and with appropriate stabilization, the discrete operator yields an M-matrix, which is sufficient to guarantee the scheme’s monotonicity and consequent $L^{\mathrm{\infty }}$-stability. These properties – consistency, stability, and monotonicity – are shown to be sufficient to prove convergence of the numerical approximation to the unique viscosity solution within the celebrated Barles–Souganidis framework. For solutions with enhanced regularity, we further establish an optimal-order error estimate of $O(\mathrm{\Delta }t+h^2)$. The rigorously established stability of the scheme provides a blueprint for a novel Physics-Constrained Neural Operator (PCNO) architecture. We prove that by emulating the scheme’s structure, the PCNO can provably break the curse of dimensionality for the challenging class of domain-to-solution mapping problems with complex topological variations, a problem for which standard learning approaches often fail. Numerical experiments for both a Hamilton-Jacobi equation with state constraints and a controlled reaction-diffusion system are presented to validate the theoretical findings and demonstrate the scheme’s effectiveness.