1 Introduction↩︎

The control theory of partial differential equations has been the subject of sustained and intensive investigation over the last several decades. For wave, heat, and dispersive equations, a large body of work has established exact, null, and approximate controllability via internal and boundary controls under a variety of geometric, analytic, and structural hypotheses; see, for instance, [1]–[4] and the references therein. In the case of wave equations, the Geometric Control Condition (GCC) plays a central role: in multidimensional domains, exact controllability is closely tied to the requirement that every generalized ray of geometric optics enters the control region within the prescribed time interval [5].

By contrast, the numerical reconstruction of controls remains much less developed. Even when controllability is established at the analytical level, it is often difficult to compute a reliable control that steers the system between prescribed states. Nonlinearities may also introduce loss of compactness, nonconvexity, stiffness, and instability in the numerical reconstruction procedure. Thus, while the existence of controls is well understood in many settings, the design of stable, accurate, and computationally tractable methods for approximating such controls remains a challenging and comparatively underdeveloped area.

A number of recent works have begun to address this issue from a constructive and numerical perspective. Least-squares, space-time, and fixed-point approaches have been developed for partial differential equations, including distributed and boundary control problems; see, for example, [6]–[10]. These works show that numerical control reconstruction is a subtle problem even in low dimensions. In higher dimensions, the difficulty is compounded by the cost of mesh-based discretizations and by the need to resolve interactions among the control region, high-frequency modes, and the conservative structure of the equations.

Machine learning provides a natural complementary set of tools for this problem, although most existing scientific machine learning methods were not designed specifically for control reconstruction. Physics-informed neural networks (PINNs), introduced in [11]–[13], approximate solutions of differential equations by training neural networks against PDE residuals, boundary conditions, and available data. PINNs have become a central tool in scientific machine learning for forward and inverse PDE problems, especially in data-sparse settings. However, the direct reconstruction of controls is a fundamentally different task: the control is an unknown input that must be chosen so that the corresponding state satisfies terminal steering constraints. Thus, the control problem imposes constraints not only on the state equation but also on the mechanism by which the unknown input enters the dynamics.

Several PINN variants have been proposed to address optimization difficulties, gradient flow stiffness, spectral bias, and gradient imbalance among different loss components; see, for example, [14]–[21]. These developments substantially improve the performance of PINNs for many forward and inverse problems. Nevertheless, they do not by themselves resolve the particular difficulty of control reconstruction, where one must recover an admissible control that steers the PDE dynamics to a prescribed terminal state and whose effect may be localized, bilinear, or strongly coupled to the state. For connections between Control Theory and Machine Learning, we refer to the discussions in [22].

In the work [23], the authors introduced WeightedPINN, an operator-decomposed PINN framework designed specifically for the numerical reconstruction of controls in nonlinear PDEs. The main idea of WeightedPINN is to assign separate adaptive weights to the different components of the controlled PDE residual, such as the time derivative, spatial differential operators, nonlinear response, and control term, as well as to the boundary, initial, and terminal constraints. This produces a control-aware residual metric that is better adapted to gradient imbalance at the operator level and to the mechanism through which the unknown control enters the equation. In this sense, WeightedPINN represents a viable attempt to adapt the PINN methodology to the control-reconstruction setting.

The present paper develops a different approach. Rather than modifying the space-time PINN loss, we exploit the spectral structure of the wave equation. The wave operator is skew-adjoint in the energy space, and its spectrum is purely imaginary. On rectangular domains, the corresponding Dirichlet eigenfunctions provide an explicit orthonormal modal basis. We use this structure to reduce the control reconstruction problem to the learning of time-dependent spectral coefficients. Both the state and the control are represented as finite sums of spatial eigenfunctions, while the modal coefficients are parameterized by one-dimensional neural networks in time. We call this a Neural-Spectral method.

This spectral viewpoint is especially natural for wave equations. The conservative dynamics are encoded mode by mode, and the projected controlled equation becomes a system of forced oscillators. The unknown control enters through the modal forcing coefficients, or, in the case of a localized actuator, through the corresponding actuator overlap matrix. Thus, the method does not ask a neural network to rediscover the spatial oscillatory structure of the wave equation. Instead, the spatial structure is built into the approximation, and the learning task is reduced to reconstructing the temporal coefficients that steer the retained modes.

The purpose of this paper is to develop and analyze a structure-preserving neural approximation framework for wave equation control reconstruction. The method is designed to exploit the modal structure of the wave equation and to provide a computationally useful control ansatz whose training loss is directly aligned with the projected controlled dynamics.

The main contributions of the paper are as follows.

• We formulate a Neural-Spectral framework for approximating controls of wave equations. The method represents the state and the control in the Dirichlet spectral basis and parameterizes the resulting temporal modal coefficients by shallow neural networks.

• We prove approximation results for the Neural-Spectral method. Under the standing assumption that an exact control exists in the relevant energy framework, we show that the control can be approximated in \(L^2((0,T);\mathcal{E})\). We then prove residual and solution approximation estimates showing that exact controlled trajectories can be approximated uniformly in time, in the energy norm, by Neural-Spectral control-state pairs.

• We state a conditional computable error estimate that separates the principal numerical error sources: spectral truncation, neural-network approximation, quadrature or collocation error, and optimization error. This estimate is not a global convergence theorem for the nonconvex optimizer, but it clarifies how the computable training loss relates to the forward-validation error.

• We identify a structural obstruction faced by standard time-stepping schemes for wave equations. Explicit Euler produces spurious amplification, implicit Euler introduces artificial dissipation, and Crank–Nicolson, although energy-preserving, compresses high-frequency phases and fails to converge to the continuous unitary wave group in the operator norm. This motivates a formulation that works directly at the level of spectral coefficients rather than relying on a direct time-stepping discretization of the full PDE.

• We provide numerical experiments in one, two, and three space dimensions. In one dimension, we test a strongly nonlinear semilinear wave equation and compare with a Galerkin–Schauder HUM fixed-point implementation. In the reported strongly nonlinear regime, this HUM-S iteration does not converge, while the Neural-Spectral method produces a forward-validated steering control. In two dimensions, we compare against a standard soft-constrained PINN on a GCC-compatible boundary-collar actuator. We also add a linear Gramian reference benchmark, where the finite-dimensional minimum-energy Galerkin control is computed explicitly and used as a classical reference solution. This benchmark verifies that the Neural-Spectral formulation can recover a control close to the Galerkin HUM control in a setting where an exact finite-dimensional reference is available. In three dimensions, we compare against a DeepONet forward surrogate combined with gradient-based control inversion. In all cases, the recovered controls are evaluated by independent forward simulation of the projected dynamics, rather than only by the training loss.

The numerical results illustrate the distinction between the present method and several natural computational baselines. In the one-dimensional strongly nonlinear test, a Galerkin–Schauder HUM fixed-point implementation loses its effective contraction property: the iterates fail to stabilize, and the resulting control does not steer the nonlinear semi-discrete dynamics accurately. By contrast, the Neural-Spectral method directly minimizes the projected residual and endpoint mismatch over a joint state-control parameterization; in the reported test, the training converges, and the recovered control gives a small forward-validation terminal error. This comparison should not be interpreted as a failure of HUM theory itself, but rather as evidence that the fixed-point implementation can become ineffective in strongly nonlinear regimes where a residual-based Neural-Spectral optimization remains useful. In the two-dimensional comparison, the PINN is trained with PDE, boundary, initial, and terminal penalties in the physical variables \((t,x,y)\). For the architecture, weighting, and optimization schedule used in this benchmark, the recovered PINN control has small amplitude and the forward-simulation terminal error remains close to the free-evolution baseline. The Neural-Spectral method, by contrast, optimizes the projected modal dynamics used in the validation step and gives a smaller terminal error on this test. This comparison highlights the potential advantage of aligning the training loss with the finite-dimensional controlled dynamics, while leaving open the possibility that more specialized PINN variants or different penalty strategies could improve the PINN baseline. In the three-dimensional comparison, the DeepONet experiment illustrates a cost–accuracy trade-off between per-instance optimization and amortized operator learning. The Neural-Spectral method is trained separately for each control problem and gives smaller forward-validated terminal errors on the two reported tests. DeepONet, on the other hand, shifts part of the computational cost to an offline training phase and therefore provides faster online inversion once the surrogate has been trained. The comparison should therefore be interpreted as evidence that the proposed method can provide higher per-instance accuracy in the tested setting, whereas operator-learning approaches may be preferable in many-query regimes where rapid online evaluation is more important than solving a single instance as accurately as possible.

The additional two-dimensional linear Gramian benchmark complements these comparisons by providing a classical non-machine-learning reference. In that test, the exact minimum-energy control of the projected Galerkin system is computed from the controllability Gramian, and the Neural-Spectral reconstruction is compared directly against this reference control. The results show close agreement with the Gramian/HUM control and small forward-validation terminal error. This benchmark should be viewed as a consistency check for the spectral control formulation in a setting where the finite-dimensional reference solution is available, rather than as a general superiority claim over all possible control solvers.

The paper is organized as follows. Section 2.1 introduces the controlled wave equation, the energy formulation, the spectral decomposition, and the Neural-Spectral ansatz. We then prove the control, residual, and solution approximation results in Section 2.2. Section 2.3 states the computable conditional error estimate. Section 3 discusses the obstruction created by standard time-stepping schemes for conservative wave dynamics. Finally, Section 4 presents the numerical experiments: a strongly nonlinear one-dimensional test, a two-dimensional comparison with a standard PINN, a two-dimensional linear Gramian reference benchmark, and a three-dimensional comparison with a DeepONet forward surrogate.

2 Controlling Wave Equations↩︎

2.1 The Neural-Spectral framework↩︎

Let \(\Omega=(0,1)^N\), and let \(\omega\subset \Omega\) be a nonempty open subset. For a given final time \(T>0\), we consider the following controlled wave equation with homogeneous Dirichlet boundary conditions: \[\label{Wave1} \begin{cases} \partial_{tt}{u} - \Delta u + g(u) + \alpha u + f\,\mathbf{1}_{\omega} = 0 & \text{in } (0,T)\times \Omega, \\[1mm] u(t,x) = 0 & \text{on } (0,T)\times \partial\Omega, \\[1mm] u(0,\cdot) = u^0,\quad \partial_t{u}(0,\cdot) = u^1 & \text{in } \Omega . \end{cases}\tag{1}\]

Here, \(u=u(t,x)\) denotes the state variable, and \(f=f(t,x)\) is an interior control acting on the subdomain \(\omega\), with \(\mathbf{1}_{\omega}\) denoting the characteristic function of \(\omega\). The constant \(\alpha \ge 0\) is real, and \(g:\mathbb{R}\to\mathbb{R}\) represents a nonlinear term. Subscripts in time denote differentiation with respect to \(t\). The objective is to construct an interior control \(f\), supported in \(\omega\), that steers the solution from the prescribed initial state \((u^0,u^1)\) to the prescribed terminal state \((u_T^0,u_T^1)\).

Control problem. Given initial data \(\left(u^0,u^1\right)\) and target data \(\left(u_T^0,u_T^1\right)\), find a control \(f\) such that \[\label{Wave2} u(T,\cdot) = u_T^0, \qquad \partial_t{u}(T,\cdot) = u_T^1 \quad \text{in } \Omega .\tag{2}\]

According to [24], under the assumption that \(\omega\) satisfies appropriate geometric conditions (e.g., being a neighborhood of a suitable part of the boundary) and \(T\) is sufficiently large, if the nonlinearity \(g\) is globally Lipschitz or satisfies an appropriate growth condition at infinity (e.g., \(\lim\limits_{|s|\to\infty} \dfrac{|g(s)|}{|s|\log^2|s|} = 0\) [25]), the system 1 is exactly controllable. That is, for any initial data \[\left(u^0,u^1\right)\in H_0^1(\Omega)\times L^2(\Omega)\] and any target data \[\left(u_T^0,u_T^1\right)\in H_0^1(\Omega)\times L^2(\Omega),\] there exists a control \(f\) such that 2 holds.

Introducing the auxiliary variable \(z = \partial_t{u}\), system 1 can be rewritten as the first-order system \[\label{Wave3} \left\{ \begin{align} &u' = z, \\ &z' = \Delta u - \alpha u -g(u)- f\,\mathbf{1}_{\omega}, \\ &u(0) = u^0,\quad z(0) = u^1 . \end{align} \right.\tag{3}\]

Next, denoting \(\Phi = (u,z)\), \(G(\Phi) = (0,g(u))\), \(\mathbf{F} = (0,f)\), and \(\Phi_0=(u^0,u^1),\) and \(\Phi_T=(u^0_T,u^1_T),\) system 3 can be written in the abstract Cauchy form \[\label{Wave4} \left\{ \begin{align} &\Phi' + A\Phi + \mathbf{F}\,\mathbf{1}_{\omega} +G(\Phi) = 0, \\ &\Phi(0,\cdot) = \Phi_0,\;\;\; \Phi(T,\cdot) = \Phi_T . \end{align} \right.\tag{4}\]

We work in the energy space \[\mathcal{E}:=H_0^1(\Omega)\times L^2(\Omega),\] equipped with the inner product \[\left\langle (u_1,z_1),(u_2,z_2)\right\rangle_{\mathcal{E}} = \int_\Omega \left(\nabla u_1\cdot \overline{\nabla u_2} +\alpha u_1\overline{u_2}\right)\,dx + \int_\Omega z_1\overline{z_2}\,dx .\] Throughout this section, complex Hilbert inner products are taken to be linear in the first argument and conjugate-linear in the second. We define the operator \(A:D(A)\subset \mathcal{E}\to \mathcal{E}\) by \[D(A):=\left(H^2(\Omega)\cap H_0^1(\Omega)\right)\times H_0^1(\Omega),\] and \[A(u,z):= \begin{pmatrix} -z\\ -\Delta u+\alpha u \end{pmatrix}.\] With this realization, \(A\) is skew-adjoint on \(\mathcal{E}\) (see [26]).

To make the spectral representation compatible with neural operator learning frameworks, it is natural and convenient to represent the states using a one-dimensional sequence rather than a complex multi-index. We establish this rigorous transition in the following lemma. The proof of this lemma is presented in Appendix 5.

With the one-dimensional re-indexing \((\lambda_n)_{n\in\mathbb{N}^*}\), \((\phi_n)_{n\in\mathbb{N}^*}\) established above, we now formulate the learning problem. The key idea is to represent both the state \(\Phi\) and the control \(\mathbf{F}\,\mathbf{1}_{\omega}\) as finite linear combinations of eigenfunctions of \(A\), with time-dependent coefficients parameterized by neural networks. Specifically, for a truncation order \(P\), we project the problem onto the finite-dimensional subspace \[\label{Wave7} \mathscr{V}_P := \operatorname{span}\left\{\phi^+_{q(j)},\, \phi^-_{q(j)} : j = 1, \ldots, P\right\} = \operatorname{span}\left\{\phi_{2j-1},\, \phi_{2j} : j = 1, \ldots, P\right\} = \{\phi_i\}_{i=1}^{2P} .\tag{7}\] We seek an approximate solution of the form \(\widetilde{\Phi}_{\theta_a} = \left(\widetilde{u}_{\theta_a}, \partial_t\widetilde{u}_{\theta_a}\right)\), where \[\widetilde{u}_{\theta_a}(t,x) = \sum_{j=1}^{P} a_j(t;\theta_a)e_j(x).\] where \[a_j(t;\theta_a)=A_j^a\cdot \sigma\left(W_j^a t+B_j^a\right),\qquad A_j^a,W_j^a,B_j^a\in\mathbb{R}^{m_a}, \qquad j=1,\ldots,P.\] Similarly, we approximate the control \(\mathbf{F}\,\mathbf{1}_{\omega}=\left(0,f\,\mathbf{1}_\omega\right)\) by \(\widetilde{\mathbf{F}}_{\theta_c}\,\mathbf{1}_{\omega}=\left(0,\widetilde{f}_{\theta_c}\mathbf{1}_\omega\right)\), where \[\widetilde{f}_{\theta_c}(t,x) =\sum_{j=1}^{P} c_j(t;\theta_c)e_j(x),\] and \[c_j(t;\theta_c) = A_j^c\cdot \sigma\left(W_j^c t+B_j^c\right), \qquad A_j^c,W_j^c,B_j^c\in\mathbb{R}^{m_c}, \qquad j=1,\ldots,P.\] Here \(\sigma\) acts component-wise and “\(\cdot\)” denotes the Euclidean inner product. Thus \[\theta_a = \left(A_i^{a},W_i^{a},B_i^{a}\right)_{i=1}^{P} \in (\mathbb{R}^{m_a})^{3P}, \qquad \theta_c = \left(A_j^c,W_j^c,B_j^c\right)_{j=1}^{P} \in (\mathbb{R}^{m_c})^{3P}.\] The full parameter vector is then \[\theta=(\theta_a,\theta_c).\] Accordingly, the state and control approximations are written as \[\widetilde{\Phi}_{\theta_a}(t,x) = \sum_{j=1}^{P}{\left(a_j(t;\theta_a), a_j'(t;\theta_a)\right)e_j(x)},\] and \[\widetilde{\mathbf{F}}_{\theta_c}(t,x)\,\mathbf{1}_\omega(x) = \left(0,\,\mathbf{1}_\omega(x)\sum_{j=1}^{P} c_j(t;\theta_c)e_j(x)\right).\] Suppose that \(\widetilde{\Phi}_{\theta}(t,x) =\displaystyle\sum_{j=1}^{P}{\left(a_j(t;\theta_a), a_j'(t;\theta_a)\right)e_j(x)}\) is a solution of 4 projected onto \(\mathscr{V}_P\) and controlled by \(\widetilde{\mathbf{F}}_{\theta_c}\). Restricting the problem 4 to \(\mathscr{V}_P\), taking the \(\mathcal{E}\)-inner product of 4 with \(\phi_i\) and noting that \[\left\langle A\widetilde{\Phi}_{\theta}(t,\cdot), \phi_i \right\rangle = \left\langle \widetilde{\Phi}_{\theta}(t,\cdot), A^*\phi_i \right\rangle = \left\langle \widetilde{\Phi}_{\theta}(t,\cdot), -A\phi_i\right\rangle = \left\langle \widetilde{\Phi}_{\theta}(t,\cdot), -\lambda_i\phi_i\right\rangle = -\overline{\lambda_i}\left\langle \widetilde{\Phi}_{\theta}(t,\cdot), \phi_i \right\rangle,\] we arrive at \[\begin{align} \left\langle \widetilde{\Phi}'_{\theta}(t,\cdot), \phi_i\right\rangle +\lambda_i\left\langle \widetilde{\Phi}_{\theta}(t,\cdot), \phi_i \right\rangle + \left\langle G\left(\widetilde{\Phi}_{\theta}(t,\cdot)\right),\phi_i\right\rangle_{\mathcal{E}} + \left\langle\widetilde{\mathbf{F}}_{\theta_c}(t,\cdot),\phi_i\right\rangle_{\mathcal{E}}=0,\quad t \in (0,T),\quad i=1,\ldots, 2P. \end{align}\] Note that we have used \(\overline{\lambda_i}=-\lambda_i\), which follows from \(\lambda_i \in \mathbf{i}\mathbb{R}\). Moreover, observe that \[\begin{align} \left\langle \widetilde{\Phi}_{\theta}(t,\cdot), \phi_{2j-1} \right\rangle_{\mathcal{E}} = \dfrac{1}{\sqrt{2}}\left(\dfrac{a_j(t;\theta)}{\overline{\lambda_{2j-1}}}\|e_j\|^2_{H_0^1(\Omega)}-a_j'(t;\theta)\|e_j\|^2_{L^2(\Omega)}\right) = \dfrac{-1}{\sqrt{2}}\left(a_j'(t;\theta)-\mathbf{i}\sqrt{\mu_j+\alpha}\,a_j(t;\theta)\right) \tag{8}\\ \left\langle \widetilde{\Phi}_{\theta}(t,\cdot), \phi_{2j} \right\rangle_{\mathcal{E}} = \dfrac{1}{\sqrt{2}}\left(\dfrac{a_j(t;\theta)}{\overline{\lambda_{2j}}}\|e_j\|^2_{H_0^1(\Omega)}-a_j'(t;\theta)\|e_j\|^2_{L^2(\Omega)}\right) = \dfrac{-1}{\sqrt{2}}\left(a_j'(t;\theta)+\mathbf{i}\sqrt{\mu_j+\alpha}\,a_j(t;\theta)\right). \tag{9} \end{align}\] Furthermore, since \(G\left(\widetilde{\Phi}_{\theta}\right)=\left(0,g\left(\widetilde{u}_{\theta}\right)\right)\), \(\widetilde{\mathbf{F}}_{\theta_c}=\left(0,\widetilde{f}_{\theta_c}\right)\), we have

\[\begin{align} \left\langle G\left(\widetilde{\Phi}_{\theta}(t,\cdot)\right),\phi_{2j-1}\right\rangle_{\mathcal{E}} &= \left\langle G\left(\widetilde{\Phi}_{\theta}(t,\cdot)\right),\phi_{2j}\right\rangle_{\mathcal{E}}=\dfrac{-1}{\sqrt{2}}\left\langle g\left(\widetilde{u}_{\theta}(t,\cdot)\right),e_j\right\rangle_{L^2(\Omega)} \\ & = \dfrac{-1}{\sqrt{2}}\left\langle g\left(\sum_{\ell=1}^{P} a_\ell(t;\theta_a)e_\ell\right),e_j\right\rangle_{L^2(\Omega)} \end{align}\] and \[\begin{align} \left\langle\widetilde{\mathbf{F}}_{\theta_c}(t,\cdot)\,\mathbf{1}_\omega,\phi_{2j-1}\right\rangle_{\mathcal{E}}&=\left\langle\widetilde{\mathbf{F}}_{\theta_c}(t,\cdot)\,\mathbf{1}_\omega,\phi_{2j}\right\rangle_{\mathcal{E}}=\dfrac{-1}{\sqrt{2}}\left\langle\widetilde{f}_{\theta_c}(t,\cdot)\,\mathbf{1}_\omega,e_j\right\rangle_{L^2(\Omega)}\\ &=\dfrac{-1}{\sqrt{2}}\left\langle\mathbf{1}_\omega\,\sum_{\ell=1}^{P} c_\ell(t;\theta_c)e_\ell,e_j\right\rangle_{L^2(\Omega)}=\dfrac{-1}{\sqrt{2}}\sum_{\ell=1}^P M_{j\ell}\,c_\ell(t;\theta_c), \end{align}\] where \[M_{j\ell} = \int_\omega e_j(x)e_\ell(x)\,dx.\] Therefore, considering the two cases \(i=2j-1\) and \(i=2j\), we obtain the same system as follows \[\begin{align} \label{Wave8} a_j''(t;\theta)+(\mu_j+\alpha)a_j(t;\theta) +\left\langle g\left(\sum_{\ell=1}^{P} a_\ell(t;\theta)e_\ell\right),e_j\right\rangle_{L^2(\Omega)} +\sum_{\ell=1}^P M_{j\ell}\,c_\ell(t;\theta_c)=0,\notag\\ t \in (0,T),\,j=1,\ldots, P. \end{align}\tag{10}\]

However, since the neural network \(a_j(t;\theta)\) is determined by the parameters \(\theta\) rather than by solving 10 , the left-hand side of 10 does not vanish in general. For each \(t \in (0,T)\) and \(j=1,\ldots, P\), we define the residual \(R_j(t;\theta)\) as this remaining quantity \[\begin{align} \label{Wave9} &R_j(t;\theta)=a_j''(t;\theta)+(\mu_j+\alpha)a_j(t;\theta) +\left\langle g\left(\sum_{\ell=1}^{P} a_\ell(t;\theta)e_\ell\right),e_j\right\rangle_{L^2(\Omega)} +\sum_{\ell=1}^P M_{j\ell}\,c_\ell(t;\theta_c). \end{align}\tag{11}\] The goal is to find \(\theta\) such that \(R_j(t;\theta) \approx 0\) for all \(j=1,\ldots, P\) and the terminal condition \(\widetilde{\Phi}_\theta(T,\cdot) \approx \Phi_T\) is satisfied. This motivates us to minimize \[\label{Wave10} \min_\theta \left[ \sum_{j=1}^P \|R_j(\cdot;\theta)\|_{L^2(0,T)}^2 + \left\| \sum_{j=1}^P \left(a_j(0;\theta), a_j'(0;\theta)\right)e_j-\Phi_0 \right\|_{\mathcal{E}}^2 + \left\| \sum_{j=1}^P \left(a_j(T;\theta),a_j'(T;\theta)\right)e_j-\Phi_T \right\|_{\mathcal{E}}^2 \right].\tag{12}\]

A minimizer \(\theta^*\) of 12 simultaneously yields an approximate control \(\widetilde{\mathbf{F}}_{\theta^*}\) and an approximate solution \(\widetilde{\Phi}_{\theta^*}\).

2.2 Approximation results↩︎

Having established the Neural-Spectral formulation, we now turn to its theoretical guarantees. In this section, we prove that the proposed architecture possesses sufficient representational capacity to approximate exact control-state pairs. We begin by specifying the regularity requirements for the neural network activation function, which are essential for establishing the density of the approximation space.

This structural identity, together with the strong convergence of spectral Galerkin projections [27], allow us to split the total approximation error into a spectral truncation contribution and a neural-network approximation contribution.

Theorem 5 guarantees that \(\widetilde{\mathbf{F}}_{\theta_c}\mathbf{1}_\omega\) approximates the exact control \(\mathbf{F}\,\mathbf{1}_\omega\) in \(L^2((0,T);\mathcal{E})\). A natural next question is whether replacing \(\mathbf{F}\,\mathbf{1}_\omega\) by \(\widetilde{\mathbf{F}}_{\theta_c}\mathbf{1}_\omega\) in system 4 produces a solution that remains close to the exact solution \(\Phi\). Thus, we aim to estimate the difference between the exact solution \(\Phi\) and the approximate solution \(\widetilde{\Phi}_{\theta}\) generated by the neural operator. Let the initial data \(\Phi_0 = (u_0, u_1) \in \mathcal{E} = H^1_0(\Omega) \times L^2(\Omega)\) and the internal control \(\mathbf{F}\,\mathbf{1}_\omega = (0, f\,\mathbf{1}_\omega)\) with \(f\in L^2((0,T) \times \omega)\) be given. Under the appropriate asymptotic growth conditions of the nonlinearity \(g\), it is a well-known result from the well-posedness theory of semi-linear wave equations (see, e.g., [25], [32]) that the system 1 admits a unique solution \[u \in C\!\left([0,T]; H^1_0(\Omega)\right) \cap C^1\!\left([0,T]; L^2(\Omega)\right).\] Consequently, the system 4 has a unique solution \[\Phi = \left(u, \partial_t u\right) \in C([0,T]; \mathcal{E}).\]

Let \(\widetilde{\Psi}_{\theta_c}=\left(\widetilde{\psi}_{\theta_c},\partial_t\widetilde{\psi}_{\theta_c}\right) \in C([0,T]; \mathcal{E})\) be the solution of the approximate problem obtained by replacing the control \(\mathbf{F}\,\mathbf{1}_\omega\) of 4 with the neural operator \(\widetilde{\mathbf{F}}_{\theta_c}\mathbf{1}_\omega\) defined as in Theorem 5, i.e., \[\label{Wave495approx} \partial_t{\widetilde{\Psi}_{\theta_c}}+A\widetilde{\Psi}_{\theta_c}+\widetilde{\mathbf{F}}_{\theta_c}\mathbf{1}_\omega+G\left(\widetilde{\Psi}_{\theta_c}\right)=0, \quad \text{with } \widetilde{\Psi}_{\theta_c}(0,\cdot)=\Phi(0,\cdot).\tag{17}\] For each \(j=1,\dots,P\), denoting \(a_j^*(t;\theta_c)=\left\langle \widetilde{\psi}_{\theta_c}(t,\cdot),e_j \right\rangle_{L^2(\Omega)}\), we prove the regularity of \(a_j^*\) in the following lemma.

Having established the \(W^{2,2}\)-regularity of the Galerkin coefficients \(a_j^*(t;\theta_c)\) under the neural control \(\widetilde{f}_{\theta_c}\,\mathbf{1}_\omega\) in Lemma 5, we are now ready to bound the residual generated by the neural approximation.

With the approximation properties for both the control input (Theorem 5) and the PDE residual (Theorem 7) independently established, we now synthesize these results. The forthcoming main theorem demonstrates that by employing a structural unification technique, a single, global neural operator architecture can be rigorously constructed to bound the overall solution error.

2.3 A conditional computable error estimate↩︎

We now record a conditional error estimate that separates the four sources of error in the fully computable Neural-Spectral method: spectral truncation, neural approximation, quadrature/collocation, and optimization error. This result is intended as a numerical analysis interpretation of the training procedure, rather than as a global convergence theorem for the nonconvex optimizer.

Let \(\Pi_P:\mathcal{E}\to\mathscr V_P\) denote the orthogonal projection onto the first \(P\) wave modes, and let \(\Phi=(u,\partial_t{u})\) be the exact controlled trajectory corresponding to an exact control \(\mathbf{F}=(0,f\mathbf{1}_\omega)\). For a trained Neural-Spectral state-control pair \[\widetilde{\Phi}_{\theta}^{P,m,Q}(t,x) = \sum_{j=1}^{P} \left(a_j(t;\theta_a),a_j'(t;\theta_a)\right)e_j(x), \qquad \widetilde{\mathbf{F}}_{\theta}^{P,m,Q}(t,x)\,\mathbf{1}_\omega =\left(0\,,\mathbf{1}_\omega\sum_{j=1}^{P}c_j(t;\theta_c)e_j(x)\right),\] we define the continuous residual \[\label{def-res} R_j(t;\theta) = a_j''(t;\theta_a)+\omega_j^2a_j(t;\theta_a) + Q_j(t;\theta_a) + \sum_{\ell=1}^P M_{j\ell}\,c_\ell(t;\theta_c), \qquad j=1,\ldots,P,\tag{30}\] where \(\omega_j=\sqrt{\mu_j+\alpha}\), \(M_{j\ell} = \displaystyle\int_\omega e_j(x)e_\ell(x)\,dx.\) and \(Q_j(\cdot;\theta)\) denotes the projected nonlinear term as follows \[Q_j(t;\theta):=\left\langle g\left(\sum_{\ell=1}^P{a_\ell(t;\theta)e_\ell}\right),e_j\right\rangle_{L^2(\Omega)}.\] In the linear case, \(Q_j\equiv0\).

The computable discrete loss is \[\label{eq:discrete95loss95PQ} \begin{align} \mathcal{L}_{P,Q}(\theta) :={}& \sum_{q=1}^{Q} w_q \sum_{j=1}^{P} \left|R_j(t_q;\theta)\right|^2 \\ &+ \lambda_{\mathrm{IC}} \left\| \widetilde{\Phi}_{\theta}^{P,m,Q}(0,\cdot)-\Pi_P\Phi_0 \right\|_{\mathcal{E}}^2 + \lambda_{\mathrm{TC}} \left\| \widetilde{\Phi}_{\theta}^{P,m,Q}(T,\cdot)-\Pi_P\Phi_T \right\|_{\mathcal{E}}^2 . \end{align}\tag{31}\] Here \(\{t_q,w_q\}_{q=1}^{Q}\) is the quadrature rule used to approximate the \(L^2(0,T)\)-norm of the residual.

We introduce the following error quantities: \[\begin{align} E_{\mathrm{spec}}(P) &:=\sup_{t\in[0,T]}\left\| \Phi(t,\cdot)-\Pi_P\Phi(t,\cdot) \right\|_{\mathcal{E}}, \tag{32} \\ E_{\mathrm{net}}(m) &:= \inf_{\theta\in\Theta_{P,m}} \left[ \sum_{j=1}^{P} \|R_j(\cdot;\theta)\|_{L^2(0,T)}^2 + \left\| \widetilde{\Phi}_\theta^{P,m,Q}(0,\cdot)-\Pi_P\Phi_0 \right\|_{\mathcal{E}}^2 + \left\| \widetilde{\Phi}_\theta^{P,m,Q}(T,\cdot)-\Pi_P\Phi_T \right\|_{\mathcal{E}}^2 \right]^{1/2}, \tag{33} \\ E_{\mathrm{quad}}(Q) &:= \sup_{\theta\in\Theta_{P,m}} \left| \sum_{j=1}^{P}\|R_j(\cdot;\theta)\|_{L^2(0,T)}^2 - \sum_{q=1}^{Q}w_q\sum_{j=1}^{P}|R_j(t_q;\theta)|^2 \right|^{1/2}, \tag{34} \\ E_{\mathrm{opt}}(\theta) &:= \left[ \mathcal{L}_{P,Q}(\theta) - \inf_{\vartheta\in\Theta_{P,m}}\mathcal{L}_{P,Q}(\vartheta) \right]_+^{1/2}. \tag{35} \end{align}\] Here \(\Theta_{P,m}\) denotes the class of Neural-Spectral ansatz functions with \(P\) modes and width \(m\).

3 Spectral Obstructions for Standard Time Discretizations↩︎

The preceding sections formulate the control reconstruction problem directly in spectral variables. We now explain why this choice is natural for wave equations. The point of this section is not to claim that classical numerical methods for wave equations fail in general. Many structure-preserving, filtered, symplectic, spectral, and discrete-HUM methods are available and can be highly effective when designed appropriately. Rather, the purpose is to isolate a simple obstruction faced by standard unfiltered one-step time-discretization schemes when they are asked to approximate the full conservative wave group uniformly over its unbounded spectrum.

For the homogeneous linear wave equation, the first-order operator \(A\) is skew-adjoint on the energy space \(\mathcal{E}\). Consequently, the exact evolution \(e^{-tA}\) is a unitary group. In spectral variables, every mode has purely imaginary frequency and evolves by a phase factor of modulus one. Thus, the conservative structure of the wave equation is encoded mode by mode.

A standard one-step time discretization replaces the exact multiplier \(e^{-\Delta t\lambda}\) by a rational amplification factor \(r(\Delta t \lambda)\). For low and moderate frequencies, this may give an accurate approximation. However, because the wave spectrum is unbounded, the behavior of \(r(i\Delta t\omega)\) as \(|\omega|\to\infty\) becomes relevant for uniform-in-frequency accuracy and, in particular, for high-frequency controllability mechanisms. The following elementary computations show three possible distortions: spurious amplification, artificial dissipation, and high-frequency phase compression.

This motivates the neural-spectral formulation used in this paper: rather than training through a direct time-stepping discretization of the full wave equation, we work with the modal dynamics themselves and learn the temporal coefficients of the state and control.

3.1 The Fundamental Obstruction: Purely Imaginary Spectrum↩︎

We begin by identifying a fundamental obstruction that affects standard time-stepping schemes when they are applied to the controlled wave system 4 . The difficulty stems from the spectral properties of the operator \(A\) established in Lemma 1: all eigenvalues \(\lambda_k=i\omega_k\) are purely imaginary, with \(|\omega_k|\to\infty\) as \(k\to\infty\). Thus, the spectrum lies on the imaginary axis and has no dissipative component. The purely imaginary nature of the spectrum of \(A\) is a direct consequence of its skew-adjointness. By Stone’s theorem (see, e.g., [37]), the skew-adjointness property of \(A\) is equivalent to \(e^{-tA}\) being a strongly continuous unitary group on \(\mathcal{E}\), which satisfies \[\label{eq:unitary95group} \left\|e^{-tA}\right\| = 1, \qquad \text{for all } t \geq 0.\tag{41}\] The above property reflects the conservative nature of the wave equation: the energy \(\left\|\Phi(t,\cdot)\right\|^2_H = \left\|\Phi(0,\cdot)\right\|^2_H\) is an exact constant of motion for the unforced linear system, for all initial data. This identity holds mode by mode: for each eigenfunction \(\phi_k\), \[|\langle e^{-tA}\Phi(0,\cdot), \phi_k\rangle_{\mathcal{E}}| = |\langle\Phi(0,\cdot),\phi_k\rangle_{\mathcal{E}}|.\] Indeed, applying 23 to \(\Psi=\Phi(0,\cdot)\), we obtain \[e^{-tA}\Phi(0,\cdot) = \sum_{k=1}^{\infty}{e^{-\lambda_kt}\left\langle \Phi(0,\cdot),\phi_k \right\rangle}\phi_k,\] which implies that for each \(k \in \mathbb{N}^*\), \[|\langle e^{-tA}\Phi(0,\cdot), \phi_k\rangle_{\mathcal{E}}| = \left|e^{-\lambda_kt}\left\langle \Phi(0,\cdot),\phi_k \right\rangle_{\mathcal{E}}\right| = \left|e^{-\lambda_kt}\right|\left|\left\langle \Phi(0,\cdot),\phi_k \right\rangle_{\mathcal{E}}\right|=\left|\left\langle \Phi(0,\cdot),\phi_k \right\rangle_{\mathcal{E}}\right|.\] The last equality follows from the fact that \(\left|e^{-\lambda_kt}\right| = \left|e^{-i\omega_kt}\right|=1\). Classical finite-difference time-stepping schemes approximate the action of the semi-group \(e^{-tA}\) over one time step \(\Delta{t} > 0\) by a rational function \(r(z)\) of \(z = \Delta{t} A\) (see, e.g., [33]). Specifically, the numerical solution after \(n\) steps takes the form \(\widetilde{\Phi}^n = r(\Delta{t} A)^n \Phi(0,\cdot)\), and the amplification factor of mode \(\phi_k\) under one step is \[r(\Delta{t}\lambda_k) = r(i\Delta{t}\omega_k) \in \mathbb{C}.\] For the scheme to faithfully reproduce 41 across the entire spectrum, one would need \(|r(i\Delta{t}\omega_k)| = 1\) for all \(k \in \mathbb{N}^*\). If instead \(|r(i\Delta{t}\omega_k)| > 1\), the amplitude of mode \(k\) grows as \(|r(i\Delta{t}\omega_k)|^n \to \infty\) as \(n \to \infty\), causing spurious exponential growth (instability). If \(|r(i\Delta{t}\omega_k)| < 1\), the amplitude decays as \(|r(i\Delta{t}\omega_k)|^n \to 0\) as \(n \to \infty\), introducing artificial dissipation that removes energy from that mode irreversibly. In either case, the scheme fails to replicate the conservative dynamics of the exact semi-group. A direct computation shows that this condition is satisfied only by the Crank–Nicolson scheme among the three classical choices:

• Explicit Euler: \(r(z) = 1 - z\), so \(|r(i\Delta{t}\omega_k)| = \sqrt{1 + \Delta{t}^2\omega_k^2} > 1\);

• Implicit Euler: \(r(z) = (1+z)^{-1}\), so \(|r(i\Delta{t}\omega_k)| = (1 + \Delta{t}^2\omega_k^2)^{-1/2} < 1\);

• Crank–Nicolson: \(r(z) = (1-z/2)(1+z/2)^{-1}\), so \(|r(i\Delta{t}\omega_k)| = 1\).

The explicit and implicit Euler schemes fail this conservative test in opposite ways: the former introduces spurious amplification, whereas the latter introduces artificial dissipation. Crank–Nicolson preserves the modulus of every spectral coefficient, but it replaces the exact phase \(\Delta t\,\omega_k\) by \[2\arctan\left(\frac{\Delta t\,\omega_k}{2}\right).\] Thus, high frequencies experience a large phase distortion, since the numerical phase saturates while the exact phase grows linearly in \(|\omega_k|\). Therefore, it still introduces a frequency-dependent phase error, which may affect high-frequency controllability.

3.2 Spectral Obstructions in Standard Time Discretization Schemes↩︎

We consider the unforced linear homogeneous case \(F=0\), \(G=0\) in system 4 . Any viable numerical method must first correctly capture the continuous evolution of the linear system \(\Phi' + A\Phi = 0\) over its unbounded spectrum before it can accurately resolve the controlled or nonlinear dynamics. Recall that \[\Phi(t,\cdot)=e^{-tA}\Phi(0,\cdot),\] where the semi-group \(e^{-tA}\) is defined by \[e^{-tA}\Psi=\sum_{j=1}^{\infty}{e^{-\lambda_jt}\left\langle \Psi, \phi_j \right\rangle}\phi_j.\] Since all the eigenvalues \(\lambda_j\) of \(A\) are purely imaginary, it follows that \(\left|e^{-t\lambda_j}\right|=1\) for all \(j \in \mathbb{N}^*\), and thus \[\left\|e^{-tA}\right\|=\sup\limits_{j\in\mathbb{N}^*}\left|e^{-t\lambda_j}\right|=1.\] Let \(0:=t_0<t_1<t_2<\ldots<t_N:=T\) be a uniform partition of the interval \([0,T]\) with time step \(\Delta{t}=\dfrac{T}{N}\). We approximate \(\Phi(t_n)\) by a numerical solution \(\widetilde{\Phi}^n\) obtained by a classical finite difference time-stepping scheme.

3.2.1 Unconditional instability of explicit Euler↩︎

Applying an explicit Euler scheme for time discretization to \(\Phi'+A\Phi=0\), we obtain \[\dfrac{\widetilde{\Phi}^n-\widetilde{\Phi}^{n-1}}{\Delta{t}}+A\widetilde{\Phi}^{n-1}=0,\quad\text{with }\widetilde{\Phi}^0=\Phi(0,\cdot),\] or equivalently, \[\widetilde{\Phi}^n=(I-\Delta{t}A)\widetilde{\Phi}^{n-1}.\] Consequently, \[\widetilde{\Phi}^n=(I-\Delta{t}A)^n\widetilde{\Phi}^0=(I-\Delta{t}A)^n\Phi(0,\cdot).\]

3.2.2 Artificial dissipation of implicit Euler↩︎

Applying an implicit Euler scheme for time discretization to \(\Phi'+A\Phi=0\), we obtain \[\dfrac{\widetilde{\Phi}^n-\widetilde{\Phi}^{n-1}}{\Delta{t}}+A\widetilde{\Phi}^{n}=0,\quad\text{with }\widetilde{\Phi}^0=\Phi(0,\cdot),\] or equivalently, \[\widetilde{\Phi}^n=(I+\Delta{t}A)^{-1}\widetilde{\Phi}^{n-1}.\] Consequently, \[\widetilde{\Phi}^n=(I+\Delta{t}A)^{-n}\widetilde{\Phi}^0=(I+\Delta{t}A)^{-n}\Phi(0,\cdot).\]

3.2.3 Lack of uniform convergence for Crank–Nicolson↩︎

Applying a Crank-Nicolson scheme for time discretization to \(\Phi'+A\Phi=0\), we obtain \[\dfrac{\widetilde{\Phi}^n-\widetilde{\Phi}^{n-1}}{\Delta{t}}+A\dfrac{\widetilde{\Phi}^{n}+\widetilde{\Phi}^{n-1}}{2}=0, \quad\text{with }\widetilde{\Phi}^0=\Phi(0,\cdot),\] or equivalently, \[\widetilde{\Phi}^n=\left(I+\dfrac{\Delta{t}}{2}A\right)^{-1}\left(I-\dfrac{\Delta{t}}{2}A\right)\widetilde{\Phi}^{n-1}.\] Consequently, \[\label{CN-discrete} \widetilde{\Phi}^n=\left(I+\dfrac{\Delta{t}}{2}A\right)^{-n}\left(I-\dfrac{\Delta{t}}{2}A\right)^n\widetilde{\Phi}^0=\left(I+\dfrac{\Delta{t}}{2}A\right)^{-n}\left(I-\dfrac{\Delta{t}}{2}A\right)^n\Phi(0,\cdot).\tag{45}\] For each \(k \in \mathbb{N}^*\) and \(\Delta{t}>0\), applying 21 to \(h:=\left(I + \dfrac{\Delta{t}}{2}A\right)^{-1} \left(I - \dfrac{\Delta{t}}{2}A\right)^1\) and \(\Psi:=\phi_k\), we have \[\left(I + \dfrac{\Delta{t}}{2}A\right)^{-1} \left(I - \dfrac{\Delta{t}}{2}A\right)\phi_k = \sum\limits_{j=1}^{\infty}{\dfrac{1 - \dfrac{\lambda_j\Delta{t}}{2}}{1 + \dfrac{\lambda_j\Delta{t}}{2}}\left\langle \phi_k, \phi_j\right \rangle_{\mathcal{E}} \phi_j}= \dfrac{1 - \dfrac{\lambda_k\Delta{t}}{2}}{1 + \dfrac{\lambda_k\Delta{t}}{2}}\phi_k = \dfrac{1 - \dfrac{i\omega_k\Delta{t}}{2}}{1 + \dfrac{i\omega_k\Delta{t}}{2}}\phi_k\] Thus, the amplification factor of mode \(\phi_k\) under one step of the scheme is the following Cayley transform \[\mu_k := r_{\mathrm{CN}}(i\omega_k\Delta{t}) := \dfrac{1 - \dfrac{i\omega_k\Delta{t}}{2}}{1 + \dfrac{i\omega_k\Delta{t}}{2}}.\] Since \(\omega_k \in \mathbb{R}\), we have \(1-\dfrac{i\omega_k\Delta{t}}{2}\) and \(1+\dfrac{i\omega_k\Delta{t}}{2}\) are complex conjugates, which leads to \(|\mu_k| = 1\). By 22 , it follows that \[\left\|r_{\mathrm{CN}}(\Delta{t} A)^n\right\| = \sup_{k \in \mathbb{N}^*}|\mu_k|^n = 1, \qquad \text{for all } n \geq 1,\] so the Crank–Nicolson operator is an isometry on \(\mathcal{E}\), which is compatible with 41 . Furthermore, \[\arg(\mu_k)=\arg\left(1-\dfrac{i\omega_k\Delta{t}}{2}\right)-\arg\left(1+\dfrac{i\omega_k\Delta{t}}{2}\right)=-2\arctan\left(\dfrac{\omega_k\Delta{t}}{2}\right).\] We write \(\mu_k = e^{-i\omega_k^h\Delta{t}}\), where \[\label{eq:num95freq} \omega_k^h := -\frac{1}{\Delta{t}}\arg(\mu_k) = \frac{2}{\Delta{t}}\arctan\!\left(\dfrac{\omega_k\Delta{t}}{2}\right).\tag{46}\] Thus, Crank–Nicolson preserves amplitudes exactly, but replaces the exact frequency \(\omega_k\) by the modified frequency \(\omega_k^h\). Recall that \(A\) is \(m\)-dissipative with dense domains. By [26], we obtain the formula for the exact solution \(\Phi\) of the linear system \(\Phi' + A\Phi = 0\) as follows \[\Phi(t,\cdot)=e^{-tA}\Phi(0,\cdot),\] where the semi-group \(e^{-tA}\) is defined by \[e^{-tA}\Psi=\sum_{j=1}^{\infty}{e^{-\lambda_jt}\left\langle \Psi, \phi_j \right\rangle}\phi_j.\] Consequently, it follows that \[\Phi(T,\cdot) = \Phi(N\Delta{t}) = e^{-TA}\Phi(0,\cdot),\] and from 45 , \[\widetilde{\Phi}^N = \left(I+\dfrac{\Delta{t}}{2}A\right)^{-N}\left(I-\dfrac{\Delta{t}}{2}A\right)^N\Phi(0,\cdot).\] Despite being an isometry, Crank–Nicolson fails for a more subtle reason:

4 Numerical Tests↩︎

4.1 A numerical experiment with strong semilinear nonlinearity↩︎

The purpose of this section is to test the Neural-Spectral framework in a strongly nonlinear regime in which the Galerkin–Schauder HUM fixed-point implementation considered here no longer converges. The experiment is not intended as a comparison with an exact semilinear HUM control, nor as a statement about all possible HUM-based numerical realizations. Rather, it compares the proposed residual-based Neural-Spectral method with one natural fixed-point HUM implementation in a regime where this implementation loses its effective contraction property.

Problem setup↩︎

We consider 1 with \[N=1, \qquad \Omega=(0,1), \qquad \alpha=0,\] and the globally Lipschitz nonlinearity \[g(u)=\beta\sin(u), \qquad \beta=20.\] Thus, \(g\) satisfies the Lipschitz hypothesis of Theorem 7 with constant \[L_g=\beta=20.\]

The control acts on the boundary collar \[\omega_\delta = \left\{ x\in\Omega: \min\{x,1-x\}<\delta \right\}, \qquad \delta=0.3,\] and we take \(T=2.5\). Equivalently, \[\omega_\delta = \Omega\setminus[\delta,1-\delta]=(0,\delta)\cup(1-\delta,\,1).\] For \(\delta=0.3\), this gives \[|\omega_\delta| = 1-(1-2\delta) = 1-0.4 = 0.6,\] so the controlled domain occupies \(60\%\) of \(\Omega\). In addition, note that the collar \(\omega_\delta\) contains a neighbourhood of both endpoints \(0\) and \(1\). Hence, every unit-speed ray of geometric optics reflecting on \(\partial\Omega\) must meet \(\omega_\delta\) within a maximum time of \[1 - 2\delta = 0.4 < T,\] implying that \((\omega_\delta, T)\) satisfies the one-dimensional Geometric Control Condition (GCC).

The initial data are \[u^0(x) = \sin(\pi x) + \frac{1}{2}\sin(2\pi x) - \frac{3}{10}\sin(3\pi x), \qquad u^1\equiv0,\] and the target state is \[\left(u_T^0,u_T^1\right)=(0,0).\] The energy norm of the initial state, which serves as the baseline for evaluating the performance of the control, is \(\left\|\left(u^0,u^1\right)\right\|_{\mathcal{E}}=3.72\).

We use the spectral truncation order \(P=8\). The \(L^2(\Omega)\)-orthonormal Dirichlet basis is \[e_j(x)=\sqrt{2}\sin(j\pi x), \qquad j=1,\ldots,P,\] with frequencies \[\omega_j=j\pi .\]

Galerkin–Schauder HUM iteration↩︎

For the linear problem \(g\equiv0\), the Hilbert Uniqueness Method (HUM) introduced by J.-L. Lions in [2] gives, via duality, a constructive minimum-energy control in the corresponding Galerkin space. For the semilinear problem, there is no analogous closed-form formula. The controllability proof in [25] is based on a Schauder/Leray–Schauder fixed-point argument, and a direct numerical implementation leads naturally to a Picard-type iteration on a sequence of linearized HUM subproblems. We use this implementation as a reference and refer to it as the Galerkin–Schauder HUM iteration, abbreviated HUM-S.

By HUM duality for the linearized problem, the adjoint state \(\phi\) satisfies the homogeneous wave equation \(\partial_{tt}\phi - \Delta\phi = 0\). Using the same truncation space \(V_P\) as in (9), we project \(\phi\) onto \(E_P=\operatorname{span}\{e_j\}_{j=1}^P\), yielding the harmonic oscillator ODE \(p_k''(t) + \omega_k^2 p_k(t) = 0\) for each temporal mode. Thus, we can parameterize the adjoint state by \[\phi(t,x) = \sum_{k=1}^{P}p_k(t)e_k(x), \qquad p_k(t) = A_k\cos(\omega_k t)+B_k\sin(\omega_k t),\] where \((A_k,B_k)_{k=1}^P\in\mathbb{R}^{2P}\) are the unknown coefficients. By HUM duality, the candidate control is \[f(t,x)=\phi(t,x)\mathbf{1}_\omega(x).\] Its projection onto \(E_P=\operatorname{span}\{e_j\}_{j=1}^P\) has modal coefficients \[c_j(t) = \sum_{k=1}^{P}M_{jk}p_k(t), \qquad M_{jk} := \int_\omega e_j(x)e_k(x)\,dx .\] The matrix \(M\in\mathbb{R}^{P\times P}\) is computed explicitly from elementary trigonometric identities.

Let \[a_j(t):=\langle u(t,\cdot),e_j\rangle_{L^2(\Omega)} .\] The projected state satisfies the nonlinear ODE system \[\label{eq:hum95s95ode} a_j''(t)+\omega_j^2a_j(t)+S_j(t)+c_j(t)=0, \qquad j=1,\ldots,P,\tag{47}\] where \[S_j(t) := \int_\Omega g\left(\sum_{\ell=1}^{P}a_\ell(t)e_\ell(x)\right)e_j(x)\,dx .\] The HUM-S iteration proceeds by freezing the nonlinear source \(S\) from the current iterate. Specifically, at iteration \(n+1\), we seek a control \(c^{(n+1)}\) and state coefficients \(a^{(n+1)}\) satisfying the linearized ODE system \[\label{eq:HUM95iteration} \frac{d^2}{dt^2} a_j^{(n+1)}(t) + \omega_j^2 a_j^{(n+1)}(t) + S_j^{(n)}(t) + c_j^{(n+1)}(t) = 0, \quad j = 1, \dots, P,\tag{48}\] where \(S_j^{(n)}(t) := \displaystyle\int_{\Omega} g\left(\sum_{\ell=1}^P a_\ell^{(n)}(t)e_\ell(x)\right)e_j(x) \, dx\). The system 48 is a linear control problem, which we solve for \(\left(A^{(n+1)}, B^{(n+1)}\right)\) by least-squares inversion of the control-to-final-state map. The state is then updated by forward integration of 47 . All ODE integrations are performed with the eighth-order Dormand–Prince method DOP853 with relative tolerance \(10^{-11}\).

For small Lipschitz constants \(L_g\), this fixed-point procedure is observed to converge geometrically. At \(\beta=20\), however, the iteration is no longer contractive in practice. Figure 1 shows the successive change \[\left\|u^{(n+1)}-u^{(n)}\right\|_{L^2((0,T);L^2(\Omega))}\] as a function of the outer iteration index \(n\). After a short transient, the difference grows rapidly and then settles into a persistent oscillatory regime around \(3.66\). Thus, the HUM-S iteration does not converge for this test. Moreover, when the final iterate’s control is inserted into a high-accuracy DOP853 integration of the full nonlinear semi-discrete system, the terminal error in the energy norm is \(2.90\), which is comparable to the size of the initial energy. Hence, this nonconverged HUM-S iterate does not provide a useful steering control in this regime.

Neural-Spectral controller↩︎

The Neural-Spectral method of Section 2.1 applies without modification. We use the ansatz \[\widetilde{u}_{\theta_a}(t,x) = \sum_{j=1}^{P}a_j(t;\theta_a)e_j(x), \qquad \widetilde{f}_{\theta_c}(t,x)\,\mathbf{1}_\omega = \sum_{j=1}^{P}c_j(t;\theta_c)e_j(x)\,\mathbf{1}_\omega,\] with shallow per-mode networks of widths \[m_a=m_c=64.\] The residual 11 becomes \[R_j(t;\theta) = a_j''(t;\theta_a) + \omega_j^2a_j(t;\theta_a) + Q_j(t;\theta_a) +\sum_{\ell=1}^P M_{j\ell}\,c_\ell(t;\theta_c), \qquad j=1,\ldots,P,\] where \[Q_j(t;\theta_a) := \int_\Omega g\left( \sum_{\ell=1}^{P}a_\ell(t;\theta_a)e_\ell(x) \right)e_j(x)\,dx .\] The nonlinear Galerkin terms \(Q_j\) are evaluated by composite trapezoidal quadrature on \(257\) spatial nodes. Training uses \(250\) collocation points, endpoint weights \[\lambda_{\mathrm{IC}}=\lambda_{\mathrm{TC}}=500,\] \(3000\) Adam steps, followed by \(12\) blocks of \(400\) L-BFGS iterations. The total wall time is approximately \(90\) seconds on a single CPU.

Result at \(\beta=20\)↩︎

A crucial point of interpretation is that exact controllability is fundamentally a feasibility statement. Even in a finite-dimensional Galerkin space, there may be many valid state-control pairs satisfying both the projected dynamics and the required endpoint conditions. Indeed, note that the overlap matrix \(M_P = \left(M_{j\ell}\right)_{j,\ell=1}^P\) is positive definite. For any nonzero vector \(\xi\in\mathbb{R}^P\), we have \[\xi^T M_P \xi = \int_{\omega_\delta} \left|\sum_{j=1}^P \xi_j e_j(x,y)\right|^2 \,dx\,dy > 0.\] The strict inequality holds because a nontrivial finite linear combination of Dirichlet eigenfunctions cannot vanish identically on the open set \(\omega_\delta\). Consequently, \(M_P\) is invertible. Thus, given a sufficiently smooth modal path \[\widetilde{u}(t,x) = \sum_{j=1}^{P}a_j(t)e_j(x)\] satisfying the prescribed initial and terminal conditions, the equation \[\sum_{\ell=1}^P M_{j\ell}\,c_\ell(t;\theta_c) = -\left(a_j''(t)+\omega_j^2a_j(t)+Q_j(t)\right), \qquad j=1,\ldots,P,\] uniquely defines a control vector that forces the projected residual to vanish. Different choices of the path \(\left\{a_j(t)\right\}_{j=1}^P\) therefore generate different admissible controls. The loss 12 does not contain a minimum-energy selection principle; it seeks a dynamically consistent control-state pair satisfying the endpoint constraints. Thus, the learned control should be judged by whether it actually steers the nonlinear dynamics to the target, rather than by its proximity to the nonconverged HUM-S iterate.

In the reported run, the Neural-Spectral optimization converges to a small residual and yields a forward-validated steering control. The final diagnostics are reported in Table 1. The self-reported terminal error from the network is exceptionally small \(\left(3.13 \times 10^{-3}\right)\). More importantly, the forward-validation error, obtained by independently integrating the full nonlinear semi-discrete system using the recovered control with a high-accuracy DOP853 solver, is \(4.33 \times 10^{-1}\). Combined with a low accumulated residual sum \(\left(5.89 \times 10^{-1}\right)\), this demonstrates that the recovered control genuinely steers the physical system. This forward-validation error is nearly one order of magnitude smaller than the terminal error \(2.90\) obtained from the nonconverged HUM-S iterate.

The learned control is shown in Figure 2. Panel (a) displays the spectral coefficients \(c_j(t)\), \(j=1,\ldots,P\). Panel (b) shows the corresponding spatial profile \[\widetilde{f}_{\theta_c}(t,x) = \sum_{j=1}^{P}c_j(t)e_j(x)\] at three intermediate times. The shaded band marks the actuator region \(\omega_{\delta}=(0,0.3) \cup (0.3,0.7)\).

4.2 2D wave equation: Neural-Spectral method vs.standard PINN↩︎

In the previous sections, the Neural-Spectral method was compared against the classical Galerkin–Schauder HUM iteration, which is the natural reference within the spectral framework. Here we step outside the spectral framework and compare against a different, widely used family of PDE-constrained learning methods: physics-informed neural networks (PINNs) in the joint \((t,x,y)\) variables. The question is operational: at comparable parameter count and wall time, how does a standard soft-constrained PINN perform on the controlled wave equation in a genuinely two-dimensional spatial domain?

Problem↩︎

We consider 1 on \[\Omega=(0,1)^2\] with the boundary-collar actuator domain \[\omega_\delta = \left\{ (x,y)\in\Omega: \min\{x,1-x,y,1-y\}<\delta \right\}, \qquad \delta=0.3.\] Equivalently, \[\omega_\delta = \Omega\setminus[\delta,1-\delta]^2 = \Omega\setminus[0.3,0.7]^2.\] Thus, \(\omega_\delta\) is the shell of width \(0.3\) adjacent to \(\partial\Omega\), and \[|\omega_\delta| = 1-(1-2\delta)^2 = 1-0.4^2 = 0.84.\] Hence, the actuator domain occupies \(84\%\) of \(\Omega\). We take \[T=7,\qquad \alpha=0,\] and no nonlinearity, \[g\equiv0.\] The controlled equation is therefore \[u_{tt}-\Delta u+f\,\mathbf{1}_{\omega_\delta}=0 \quad\text{on }(0,T)\times\Omega, \qquad u|_{\partial\Omega}=0.\]

The actuator \(\omega_\delta\) is compatible with the Geometric Control Condition (GCC). Indeed, a unit-speed ray that has not entered \(\omega_\delta\) must remain inside the central square \([\delta,1-\delta]^2\). Let \(v = (v_1, v_2)\) be the direction vector of this unit-speed ray, which implies its Euclidean norm is \(\|v\|_{\mathbb{R}^2} = \sqrt{v_1^2 + v_2^2} = 1\). Since \(v_1^2 + v_2^2 \le 2\max\left\{|v_1|^2, |v_2|^2\right\}\), it immediately follows that \[\max\{|v_1|,|v_2|\}\ge \frac{1}{\sqrt2}.\] Moreover, the maximum distance the ray can travel along any coordinate axis while remaining inside the central square is exactly its side length, \(1-2\delta\). Therefore, the total time \(T_\delta\) the ray can spend inside the central square is bounded by the time it takes to exit along the fastest coordinate direction, that is, \[T_\delta \le \min\left\{ \frac{1-2\delta}{|v_1|}, \frac{1-2\delta}{|v_2|} \right\} = \frac{1-2\delta}{\max\{|v_1|,|v_2|\}} \le \sqrt2(1-2\delta).\] For \(\delta=0.3\), this gives \[T_\delta\le0.4\sqrt2\approx0.566.\] Thus, the choice \(T=7\) is safely above this sufficient GCC time.

The initial data are \[u^0(x,y) = \sin(\pi x)\sin(\pi y) + \frac{1}{2}\sin(\pi x)\sin(2\pi y), \qquad u^1\equiv0,\] and the target state is \[\left(u_T^0,u_T^1\right) \equiv (0,0).\] The \(L^2(\Omega)\)-orthonormal Dirichlet sine basis is \[e_{j,k}(x,y) = 2\sin(j\pi x)\sin(k\pi y), \qquad j,k\in\mathbb{N}^*,\] with eigenvalues \[\mu_{j,k} = \left(j^2+k^2\right)\pi^2.\] We use the spectral truncation \[j,k\in\{1,2,3,4\},\] that is, \(P_{\max}^2=16\) retained modes for the spectral method. The PINN itself is not spectrally truncated during training.

Let \(m=(j,k)\) and \(n=(j',k')\) denote multi-indices in \(\{1,2,3,4\}^2\). If the lifted control field is written as \[\widetilde{f}(t,x,y) = \sum_n c_n(t)e_n(x,y),\] then the actual applied control is \[\mathbf{1}_{\omega_\delta}\widetilde{f}.\] The projected dynamics are \[\label{eq:test6-spec-ode} \ddot a_m(t)+\mu_m a_m(t)+\sum_n M_{mn}c_n(t)=0, \qquad m\in\{1,2,3,4\}^2,\tag{49}\] where \[M_{mn} = \int_{\omega_\delta} e_m(x,y)e_n(x,y)\,dx\,dy.\] For \(\delta=0.3\), the overlap matrix \(M\) has eigenvalues in \[[0.085,1.000]\] and condition number approximately \(11.7\), so the retained finite-dimensional control problem is well-conditioned.

The energy norm of the initial data is \[\left\|\left(u^0,u^1\right)\right\|_{\mathcal{E}}^2 = 2\pi^2\left(\frac{1}{2}\right)^2 + 5\pi^2\left(\frac{1}{4}\right)^2 \approx 8.02,\] and therefore \[\left\|\left(u^0,u^1\right)\right\|_{\mathcal{E}}\approx2.83.\] Since the target is zero and the uncontrolled linear wave equation conserves energy, any control that has a negligible effect on the dynamics will produce a forward-simulation terminal error close to this value.

Method 1: Neural-Spectral method↩︎

For each retained mode \((j,k)\in\{1,\ldots,4\}^2\), we use two independent shallow \(\tanh\) networks of width \(64\), \[a_{j,k}(t;\theta_a), \qquad c_{j,k}(t;\theta_c),\] for the state coefficient and the lifted-control coefficient. The parameters are trained by minimizing the unsupervised projected residual loss \[\begin{align} &\mathcal{L}_{\mathrm{spec}}(\theta_a,\theta_c) = \int_0^T \sum_m \left| \ddot a_m(t;\theta_a) + \mu_m a_m(t;\theta_a) + \sum_n M_{mn}c_n(t;\theta_c) \right|^2 \,dt \\ &+\lambda\left(\left\| \sum_{m} \left(a_m(0;\theta_a), \dot{a}_m(0;\theta_a)\right)e_j-\left(u^0, u^1\right) \right\|_{\mathcal{E}}^2 + \left\| \sum_{m}\left(a_m(T;\theta_a),\dot{a}_m(T;\theta_a)\right)e_j-\left(u_T^0, u_T^1\right) \right\|_{\mathcal{E}}^2\right), \end{align}\] with \(\lambda=100\). The initial and terminal mismatch terms are measured in the energy norm. The residual is the projected PDE residual of 49 ; consequently, driving this residual to zero together with the endpoint constraints yields a true solution of the retained truncated system.

In the present experiment, \[u^1\equiv0, \qquad (u_T^0,u_T^1)=(0,0),\] so this loss reduces to \[\begin{align} \mathcal{L}_{\mathrm{spec}}(\theta_a,\theta_c) =&{} \int_0^T \sum_m \left| \ddot a_m(t;\theta_a) + \mu_m a_m(t;\theta_a) + \sum_n M_{mn}c_n(t;\theta_c) \right|^2 \,dt \\ &+\lambda\left(\left\| \sum_{m} \left(a_m(0;\theta_a), \dot{a}_m(0;\theta_a)\right)e_j-\left(u^0, 0\right) \right\|_{\mathcal{E}}^2 + \left\| \sum_{m}\left(a_m(T;\theta_a),\dot{a}_m(T;\theta_a)\right)e_j \right\|_{\mathcal{E}}^2\right). \end{align}\]

The training schedule uses \(250\) time-collocation points, \(2000\) Adam steps, followed by \(25\) blocks of \(300\) L-BFGS iterations. The total number of trainable parameters is \(6144\).

Method 2: PINN in \((t,x,y)\)↩︎

The PINN baseline uses two fully connected feedforward networks with \(\tanh\) activations, width \(40\), and depth \(4\): \[u_\theta:\mathbb{R}^3\to\mathbb{R}, \qquad f_\theta:\mathbb{R}^3\to\mathbb{R}.\] The first network represents the state, and the second represents the control field. The actuator is imposed through the indicator \[\chi_{\omega_\delta}(x,y) = \mathbf{1}_{\{\min\{x,\,1-x,\,y,\,1-y\}<\delta\}}.\]

Let \[\mathcal{G}_{\mathrm{int}} = \{(t_q,x_q,y_q)\}_{q=1}^{N_{\mathrm{int}}} \subset (0,T)\times\Omega\] be the interior collocation set for the PDE residual, \[\mathcal{G}_{\mathrm{bc}} = \{(t_q,x_q,y_q)\}_{q=1}^{N_{\mathrm{bc}}} \subset (0,T)\times\partial\Omega\] be the boundary collocation set for the homogeneous Dirichlet condition, and \[\mathcal{G}_{\mathrm{ep}} = \{(x_q,y_q)\}_{q=1}^{N_{\mathrm{ep}}} \subset \Omega\] be the endpoint grid used at \(t=0\) and \(t=T\). To solve this control problem using the standard PINN framework (as originally introduced by Raissi et al. [13]), we reformulate the constrained PDE system into an unconstrained multi-objective optimization problem. The physical laws and boundary constraints are enforced as soft penalty terms in the empirical loss function. The explicit PINN loss is \[\begin{align} \mathcal{L}_{\mathrm{PINN}}(\theta) ={}& \Lambda_{\mathrm{PDE}} \frac{1}{N_{\mathrm{int}}} \sum_{q=1}^{N_{\mathrm{int}}} \left| \partial_{tt}u_\theta(t_q,x_q,y_q) - \partial_{xx}u_\theta(t_q,x_q,y_q) - \partial_{yy}u_\theta(t_q,x_q,y_q) + f_\theta(t_q,x_q,y_q)\chi_{\omega_\delta}(x_q,y_q) \right|^2 \\ &+ \Lambda_{\mathrm{BC}} \frac{1}{N_{\mathrm{bc}}} \sum_{q=1}^{N_{\mathrm{bc}}} \left| u_\theta(t_q,x_q,y_q) \right|^2 \\ &+ \Lambda_{\mathrm{IC}} \frac{1}{N_{\mathrm{ep}}} \sum_{q=1}^{N_{\mathrm{ep}}} \left[ \left| u_\theta(0,x_q,y_q)-u^0(x_q,y_q) \right|^2 + \left| \partial_t u_\theta(0,x_q,y_q)-u^1(x_q,y_q) \right|^2 \right] \\ &+ \Lambda_{\mathrm{TC}} \frac{1}{N_{\mathrm{ep}}} \sum_{q=1}^{N_{\mathrm{ep}}} \left[ \left| u_\theta(T,x_q,y_q)-u_T^0(x_q,y_q) \right|^2 + \left| \partial_t u_\theta(T,x_q,y_q)-u_T^1(x_q,y_q) \right|^2 \right]. \end{align}\] In the present experiment, \[u^1\equiv0, \qquad (u_T^0,u_T^1)=(0,0),\] so this loss reduces to \[\begin{align} \mathcal{L}_{\mathrm{PINN}}(\theta) ={}& \Lambda_{\mathrm{PDE}} \frac{1}{N_{\mathrm{int}}} \sum_{q=1}^{N_{\mathrm{int}}} \left| \partial_{tt}u_\theta(t_q,x_q,y_q) - \partial_{xx}u_\theta(t_q,x_q,y_q) - \partial_{yy}u_\theta(t_q,x_q,y_q) + f_\theta(t_q,x_q,y_q)\chi_{\omega_\delta}(x_q,y_q) \right|^2 \\ &+ \Lambda_{\mathrm{BC}} \frac{1}{N_{\mathrm{bc}}} \sum_{q=1}^{N_{\mathrm{bc}}} \left| u_\theta(t_q,x_q,y_q) \right|^2 \\ &+ \Lambda_{\mathrm{IC}} \frac{1}{N_{\mathrm{ep}}} \sum_{q=1}^{N_{\mathrm{ep}}} \left[ \left| u_\theta(0,x_q,y_q)-u^0(x_q,y_q) \right|^2 + \left| \partial_t u_\theta(0,x_q,y_q) \right|^2 \right] \\ &+ \Lambda_{\mathrm{TC}} \frac{1}{N_{\mathrm{ep}}} \sum_{q=1}^{N_{\mathrm{ep}}} \left[ \left| u_\theta(T,x_q,y_q) \right|^2 + \left| \partial_t u_\theta(T,x_q,y_q) \right|^2 \right]. \end{align}\] We use \[\Lambda_{\mathrm{PDE}}=50, \qquad \Lambda_{\mathrm{BC}} = \Lambda_{\mathrm{IC}} = \Lambda_{\mathrm{TC}} = 100.\] The total number of trainable parameters is \(6962\).

After training, the recovered PINN control \(f_\theta(t,x,y)\) is evaluated by forward validation. For this purpose, we project the applied control \(f_\theta\mathbf{1}_{\omega_\delta}\) onto the retained sine basis: \[d_{j,k}(t) = \int_{\omega_\delta} f_\theta(t,x,y)e_{j,k}(x,y)\,dx\,dy, \qquad j,k\in\{1,2,3,4\}.\] The validation dynamics are then \[\ddot a_{j,k}(t)+\mu_{j,k}a_{j,k}(t)+d_{j,k}(t)=0, \qquad j,k\in\{1,2,3,4\}.\] Thus, both methods are evaluated by inserting their recovered controls into the corresponding projected wave dynamics and measuring the terminal energy error.

Results↩︎

The single-problem comparison is summarized in Table 2, and the error-versus-wall-time trajectories are shown in Figure 3.

What we found↩︎

(F1) The Neural-Spectral method is more accurate on this test. It reduces the forward-simulation terminal error from the free-evolution baseline \[\left\|u^0\right\|_{\mathcal{E}}\approx2.83\] to \[4.84\times10^{-1}.\] The PINN ends at \[2.78,\] which is only a modest reduction from the free-evolution baseline. Thus, on this benchmark, the Neural-Spectral method improves the terminal error by a factor of approximately \(5.7\) relative to the PINN.

(F2) The standard PINN baseline returns a weak control in this test. The recovered PINN control has applied \(L^2\)-norm \(0.21\), about \(57\) times smaller than that recovered by the Neural-Spectral method. The resulting learned pair behaves close to a free-evolution solution: it satisfies the wave equation, boundary condition, and initial condition reasonably well, but the recovered control is not strong enough to steer the state close to the terminal target.

(F3) The Neural-Spectral loss is aligned with the validated dynamics. The Neural-Spectral method optimizes the projected residual of 49 , which is the same finite-dimensional dynamics used in the forward-validation step. Thus, for the retained modal system, a small projected residual together with a small endpoint mismatch is directly relevant to the reported forward-validation error. By contrast, the PINN baseline uses a soft multi-objective loss in the physical variables, where the PDE residual, boundary condition, initial condition, and terminal condition enter as competing penalty terms. The observed weak-control solution is consistent with known optimization difficulties for PINNs, including gradient imbalance among loss components [21]. We do not claim that all PINN variants would behave similarly; rather, this experiment shows that the standard soft-constrained PINN used here is less effective than the Neural-Spectral formulation on this benchmark.

4.3 3D wave equation: Neural-Spectral method vs. DeepONet↩︎

The previous sections compared the Neural-Spectral method against either the classical Galerkin–Schauder HUM iteration in one-dimensional spatial domain or, in the two-dimensional case, a standard PINN. We now compare it against a different family of neural PDE solvers: operator-learning methods, with DeepONet [38] as a representative example. The test is performed on the linear three-dimensional wave equation. The purpose is operational: we compare a per-problem solver, namely the Neural-Spectral method, with an operator-learning surrogate that is trained once and then reused. Thus, the two methods occupy different positions in the cost-accuracy trade-off. We first compare them on a single three-dimensional control problem, and then examine what changes when the same trained DeepONet is reused on a second related problem.

Problem↩︎

We consider 1 on the three-dimensional unit cube \[\Omega=(0,1)^3,\] with no nonlinearity, \(g\equiv0\), and with \(\alpha=0\). The initial velocity is zero, \(u^1\equiv0\), and the target state is \[\left(u_T^0,u_T^1\right)=(0,0).\] The control acts on the boundary collar \[\omega_\delta = \left\{ (x,y,z)\in\Omega: \min\{x,1-x,y,1-y,z,1-z\}<\delta \right\}, \qquad \delta=0.3,\] and we take \(T=5\). Equivalently, \[\omega_\delta = \Omega\setminus[\delta,1-\delta]^3.\] For \(\delta=0.3\), this gives \[|\omega_\delta| = 1-(1-2\delta)^3 = 1-0.4^3 = 0.936,\] so the controlled domain occupies approximately \(93.6\%\) of the cube.

This choice is compatible with the Geometric Control Condition (GCC). Indeed, a unit-speed ray that has not yet entered \(\omega_\delta\) must remain inside the central cube \([\delta,1-\delta]^3\). If \(v=(v_1,v_2,v_3)\) is its direction, then \(\|v\|_{\mathbb{R}^3}=1\). As \(v_1^2+v_2^2+v_3^2 \le 3\max\limits_{1\le i\le3}|v_i|^2\), it holds that \[\max\limits_{1\le i\le3}|v_i|\ge \dfrac1{\sqrt3}.\] Moreover, the maximum distance the ray can travel along any coordinate axis while remaining inside the central square is exactly its side length, \(1-2\delta\). Therefore, the total time \(T_\delta\) the ray can spend inside the central square is strictly bounded by the time it takes to exit along the fastest coordinate direction, that is, \[T_\delta \le \min\limits_{1\le i\le3}\frac{1-2\delta}{|v_i|} = \frac{1-2\delta}{\max\limits_{1\le i\le3}|v_i|} \le \sqrt3(1-2\delta).\] For \(\delta=0.3\), this gives \[T_\delta\le 0.4\sqrt3\approx0.693.\] Thus, the choice \(T=5\) is safely above this sufficient GCC time.

The \(L^2(\Omega)\)-orthonormal Dirichlet sine basis is \[e_{j,k,\ell}(x,y,z) = \left(\sqrt2\right)^3 \sin(j\pi x)\sin(k\pi y)\sin(\ell\pi z), \qquad j,k,\ell\in\mathbb{N}^*,\] with eigenvalues \[\mu_{j,k,\ell} = \left(j^2+k^2+\ell^2\right)\pi^2.\] We use the spectral truncation \[1\le j,k,\ell\le P_{\max}=3,\] which gives \[n_{\mathrm{mode}}=P_{\max}^3=27\] retained modes.

Let \(m=(j,k,\ell)\) and \(n=(j',k',\ell')\) denote multi-indices in \(\{1,2,3\}^3\). If the lifted control field is written as \[\widetilde{f}(t,x,y,z) = \sum_n c_n(t)e_n(x,y,z),\] then the actual applied control is \[\mathbf{1}_{\omega_\delta}\widetilde{f}.\] The projected dynamics are \[\label{eq:test8-spec-ode} \ddot a_m(t)+\mu_m a_m(t)+\sum_n M_{mn}c_n(t)=0, \qquad m\in\{1,2,3\}^3,\tag{50}\] where \[M_{mn} = \int_{\omega_\delta} e_m(x,y,z)e_n(x,y,z)\,dx\,dy\,dz.\] For this actuator, the overlap matrix \(M\) has eigenvalues in \[[0.125,0.999]\] and condition number approximately \(8.0\), so the retained finite-dimensional control problem is well-conditioned. We use two test problems, each with nonzero spectral content in four low modes: \[\text{Test \#1:}\quad \left\|u^0\right\|_{\mathcal{E}}=1.15, \qquad \text{Test \#2:}\quad \left\|u^0\right\|_{\mathcal{E}}=1.06.\]

Method 1: Neural-Spectral method↩︎

For each retained spectral mode \(m\), we represent both the state coefficient \(a_m(t)\) and the lifted-control coefficient \(c_m(t)\) by shallow \(\tanh\) networks of width \(64\). The parameters \(\theta_a\) and \(\theta_c\) are trained by minimizing the projected residual loss \[\begin{align} \mathcal{L}_{\mathrm{spec}}(\theta_a,\theta_c) ={}& \int_0^T \sum_{j=1}^{3}\sum_{k=1}^{3}\sum_{\ell=1}^{3} \left| \ddot a_{j,k,\ell}(t;\theta_a) + \left(j^2+k^2+\ell^2\right)\pi^2 a_{j,k,\ell}(t;\theta_a) \right.\\ &\left. + \sum_{j'=1}^{3}\sum_{k'=1}^{3}\sum_{\ell'=1}^{3} M_{(j,k,\ell),(j',k',\ell')} c_{j',k',\ell'}(t;\theta_c) \right|^2 \,dt \\ &+ \Lambda \sum_{j=1}^{3}\sum_{k=1}^{3}\sum_{\ell=1}^{3} \left[ \left(j^2+k^2+\ell^2\right)\pi^2 \left|a_{j,k,\ell}(0;\theta_a)-u^0_{j,k,\ell}\right|^2 + \left|\dot{a}_{j,k,\ell}(0;\theta_a)\right|^2 \right] \\ &+ \Lambda \sum_{j=1}^{3}\sum_{k=1}^{3}\sum_{\ell=1}^{3} \left[ \left(j^2+k^2+\ell^2\right)\pi^2 \left|a_{j,k,\ell}(T;\theta_a)\right|^2 + \left|\dot{a}_{j,k,\ell}(T;\theta_a)\right|^2 \right], \end{align}\] where the initial and terminal mismatch terms are measured in the energy norm and \(\Lambda=100\). We use \(250\) time-collocation points, \(2000\) Adam steps, followed by \(18\) blocks of \(250\) L-BFGS iterations. The total number of trainable parameters is \(10\,368\). This method is trained from scratch for each test problem.

Method 2: DeepONet forward surrogate and gradient-descent control↩︎

The DeepONet is trained to approximate the finite-dimensional control-to-terminal-state map \[\mathcal{G}: c \longmapsto \left(u(T,\cdot),u_t(T,\cdot)\right),\] where \(c\) denotes the vector of piecewise-linear modal control coefficients and the initial data are zero. Since the wave equation is linear in this experiment, the terminal state for arbitrary initial data decomposes as \[(u(T,\cdot),u_t(T,\cdot)) = \left(u_{\mathrm{free}}(T,\cdot),v_{\mathrm{free}}(T,\cdot)\right) + \mathcal{G}(c).\] Thus, the same trained DeepONet can be reused for different initial and target states on the same domain. No HUM solver is used to generate the training labels; all labels are obtained from forward simulations of random controls.

4.3.0.1 Control parameterization.

The modal control is parameterized as a piecewise-linear interpolant on \[N_{t,\mathrm{sub}}=11\] time knots, with one curve per retained spectral mode. Thus, the branch input is \[b = \left(c_{j,k,\ell}(t_n)\right)_{ 1\le j,k,\ell\le3,\;0\le n\le10} \in\mathbb{R}^{297}.\]

4.3.0.2 Network architecture.

We use a two-channel DeepONet of the form \[\mathcal{G}_\theta(c)(x,y,z) = \left( \beta_\theta(c)\cdot \tau_\theta^u(x,y,z), \, \beta_\theta(c)\cdot \tau_\theta^v(x,y,z) \right).\] The branch network \[\beta_\theta:\mathbb{R}^{297}\to\mathbb{R}^q\] is a four-layer \(\tanh\) MLP of width \(128\), with \(q=160\). The trunk network \[\tau_\theta:\mathbb{R}^3\to\mathbb{R}^{2q}\] is a five-layer \(\tanh\) MLP of width \(128\). The total number of trainable parameters is \(183\,138\), about \(18\) times that of the Neural-Spectral method.

4.3.0.3 Offline training phase.

We generate \(1500\) random control samples by drawing \[c_{j,k,\ell}(t_n)\sim\mathcal{N}\left(0,0.8^2\right)\] independently for each retained mode and time knot. Each sampled control is forward-simulated from zero initial data under 50 to obtain \[\left(a_m(T),\dot{a}_m(T)\right).\] The terminal fields \[u(T,\cdot), \qquad u_t(T,\cdot)\] are then reconstructed by spectral expansion at \(100\) random query points in \(\Omega\). Since the forcing is piecewise linear in time, the forced linear oscillator response can be computed in closed form; no ODE solver is needed for this data-generation step. Generating all \(1500\) training samples takes approximately \(0.4\) seconds.

The DeepONet is trained by minimizing the mean-squared error of the two terminal channels at the query points. The training schedule consists of \(10\,000\) Adam steps with batch size \(128\), initial learning rate \(2\cdot10^{-3}\), and decay factor \(0.5\) every \(1500\) steps. The offline training time is approximately \(17\) seconds.

4.3.0.4 Online inversion phase.

Given a test problem \[\left(u^0,u^1,u_T^0,u_T^1\right),\] the free-evolution terminal state is computed in closed form and evaluated on a fixed \(9^3=729\)-point grid in \(\Omega\). The control parameter vector \[c\in\mathbb{R}^{297}\] is initialized at zero and optimized by \(1500\) Adam steps with initial learning rate \(0.05\), with learning-rate decay. The loss is the mean-squared mismatch between \(\mathcal{G}_\theta(c)\) and the required terminal offset on the evaluation grid. During this phase, the DeepONet weights are frozen, and gradients are taken only with respect to the control vector \(c\). The online inversion time per test problem is approximately \(6\)–\(8\) seconds.

4.3.0.5 Forward validation.

For both methods, the recovered control is inserted into the true truncated modal system 50 , starting from the actual initial data. The system is then integrated with DOP853 to a relative tolerance \(10^{-12}\). We report the model-independent terminal error \[\left\|\Phi(T,\cdot)-\Phi_T\right\|_{\mathcal{E}} = \left( \sum_m \mu_m \left(a_m(T)-u^0_{T,m}\right)^2 + \left(\dot{a}_m(T)-u^1_{T,m}\right)^2 \right)^{1/2}.\] This is the same forward-validation metric used in the previous numerical tests.

Main result: single-problem accuracy↩︎

The main comparison is summarized in Table 3 and in Figure 4.

On Test #1, the Neural-Spectral method gives a terminal energy-norm error \[4.39\times10^{-1},\] whereas DeepONet gives \[1.36\times10^{0}.\] Thus, on this single problem, the Neural-Spectral method gives the smaller forward-validation terminal error, whereas DeepONet is faster in online wall time because part of the computation has been shifted into an offline-trained forward surrogate.

One possible source of the DeepONet error in this test is surrogate accuracy in the region of control space reached by the online inversion. During the online phase, the surrogate mismatch \[\|\mathcal{G}_\theta(c)-\text{target offset}\|^2\] on the evaluation grid is reduced to approximately \(10^{-2}\). However, when the resulting control is inserted into the true projected wave dynamics 50 , the forward-validation terminal error is \(1.36\times10^{0}\). This indicates that the optimized control is effective for the learned surrogate but less effective for the true finite-dimensional dynamics in this experiment. Such behavior is a known concern when learned forward operators are used inside inverse or control loops: the inversion step may move toward regions of input space where the surrogate is less accurate than its average training error suggests.

The Neural-Spectral method avoids this particular surrogate-error mechanism because it optimizes the residual of the same projected dynamics used in the forward-validation step. Its remaining errors come from finite spectral truncation, neural parameterization, numerical quadrature, and nonconvex optimization.

Additional observation: amortization across multiple problems↩︎

The above comparison treats the two methods on the same footing: each is asked to solve a single control problem. The natural regime for operator-learning methods, however, is the many-query setting, where the same trained surrogate is reused across multiple problems on the same domain. We therefore ran a second test, Test #2, with different initial data satisfying \[\left\|u^0\right\|_{\mathcal{E}}=1.06.\] The initial data places mass on the cross modes \[(2,2,1),\qquad (1,2,2),\qquad (3,1,1),\] as well as on \((1,1,1)\).

With two test problems, the Neural-Spectral method gives the smaller forward-validation terminal error on each individual problem, while DeepONet has the smaller cumulative wall time because its offline training cost is amortized. In general, if \(W_{\mathrm{spec}}\) denotes the average wall time of a single Neural-Spectral solve, \(W_{\mathrm{off}}\) the one-time DeepONet offline training cost, and \(W_{\mathrm{inf}}\) the DeepONet online inference cost per new problem, then the two cumulative costs over \(N\) related problems are approximately \[N W_{\mathrm{spec}} \qquad\text{and}\qquad W_{\mathrm{off}}+N W_{\mathrm{inf}}.\] Thus, DeepONet becomes faster when \[N> \frac{W_{\mathrm{off}}}{W_{\mathrm{spec}}-W_{\mathrm{inf}}}.\] With the present timings, this threshold is below one problem, so DeepONet is already faster in wall time on the first test. The observed trade-off in these experiments is therefore between per-instance forward-validation accuracy and amortized online speed.

4.4 A two-dimensional linear Gramian reference benchmark↩︎

In this section, we include a linear reference benchmark for which the finite-dimensional Galerkin control can be computed explicitly through a controllability Gramian. The purpose of this test is different from the previous PINN and DeepONet comparisons. Here, we compare the Neural-Spectral reconstruction with a classical Gramian/HUM reference control for the same projected linear dynamics. This provides a non-machine-learning benchmark against which the accuracy of the recovered control can be assessed.

We consider the linear wave equation \[\label{eq:linear952d95wave95test} \begin{cases} u_{tt}-\Delta u+f\,\mathbf{1}_{\omega_\delta}=0, & (t,x,y)\in (0,T)\times \Omega,\\ u=0, & (t,x,y)\in (0,T)\times \partial\Omega,\\ u(0,x,y)=u^0(x,y),\qquad u_t(0,x,y)=u^1(x,y), & (x,y)\in \Omega, \end{cases}\tag{51}\] where \[\Omega=(0,1)^2,\qquad g\equiv 0,\qquad \alpha=0.\] The target is the null state \[u(T,\cdot)=0,\qquad u_t(T,\cdot)=0.\] We choose the boundary-collar actuator \[\omega_\delta = \Omega\setminus [\delta,1-\delta]^2 = \{(x,y)\in \Omega:\min(x,1-x,y,1-y)<\delta\},\] with \[\delta=0.2,\qquad T=4.\] The uncontrolled region is the central square \[[\delta,1-\delta]^2=[0.2,0.8]^2.\] For this collar geometry, any geometric ray traveling at unit speed that has not yet entered the control region must remain inside the central uncontrolled square. Thus, the distance such a ray can cover before hitting the collar is bounded by the diameter of this square, namely \((1-2\delta)\sqrt{2}\). Thus, the time before the unit-speed ray reaches the collar is also bounded by \[(1-2\delta)\sqrt 2=0.6\sqrt 2\approx 0.849.\] Since \(T=4\), the pair \((\omega_\delta,T)\) satisfies the Geometric Control Condition (GCC) with a large safety margin. This elementary bound is used only to justify the benchmark geometry; no claim of sharpness is needed here.

For the initial condition, we take a finite superposition of Dirichlet modes: \[u^0(x,y) = E_{1,2}(x,y) + 0.7E_{3,1}(x,y) + 0.4E_{2,3}(x,y), \qquad u^1(x,y)=0,\] where \[E_{m,n}(x,y)=2\sin(m\pi x)\sin(n\pi y)\] is the \(L^2(\Omega)\)-normalized Dirichlet eigenfunction associated with the eigenvalue \[\mu_{m,n}=\pi^2\left(m^2+n^2\right).\] All three active modes are already included at the smallest truncation \(P=8\). Indeed, the largest active eigenvalue is \(\mu_{2,3}=13\pi^2\), and the first eight Dirichlet modes on the square contain all modes with \(m^2+n^2\le 13\). Hence, the initial data is represented exactly at each truncation level \(P\in\{8,16,32\}\) used below, and the projected initial energy norm \(\|Y_0\|_{\mathcal{E}_P}\approx 10.87\) is independent of \(P\).

Using the spectral re-indexing bijection \(q: \mathbb{N}^* \to \mathbb{N}^* \times \mathbb{N}^*\) established in the proof of Lemma  2 (see Appendix 5), we map the two-dimensional eigenvalues \(\mu_{m,n}\) in ascending order. Let \(e_j(x,y) := E_{q(j)}(x,y)\) denote the corresponding sorted sequence of the first \(P\) Dirichlet eigenfunctions. We then write the Galerkin approximation as \[u_P(t,x,y)=\sum_{j=1}^P a_j(t)e_j(x,y), \qquad f_P(t,x,y)=\sum_{j=1}^P c_j(t)e_j(x,y).\] Projecting 51 onto \(\operatorname{span}\{e_1,\ldots,e_P\}\) gives \[a_j''(t)+\mu_j a_j(t) + \sum_{\ell=1}^P M_{j\ell}c_\ell(t)=0, \qquad j=1,\ldots,P,\] where \[M_{j\ell} = \int_{\omega_\delta} e_j(x,y)e_\ell(x,y)\,dx\,dy.\] In vector form, with \[a(t)=(a_1(t),\ldots,a_P(t))^T,\qquad c(t)=(c_1(t),\ldots,c_P(t))^T,\] and \[\Lambda_P=\operatorname{diag}(\mu_1,\ldots,\mu_P), \qquad M_P=\left(M_{j\ell}\right)_{j,\ell=1}^P,\] the projected system is \[a''(t)+\Lambda_P a(t)+M_Pc(t)=0.\] Equivalently, setting \[Y(t)= \begin{pmatrix} a(t)\\ a'(t) \end{pmatrix},\] we obtain the first-order controlled system \[\label{eq:linear952d95first95order95system} Y'(t)=A_PY(t)+B_Pc(t),\tag{52}\] where \[A_P= \begin{pmatrix} 0&I_P\\ -\Lambda_P&0 \end{pmatrix}, \qquad B_P= \begin{pmatrix} 0\\ -M_P \end{pmatrix}.\] For this boundary-collar actuator, \(M_P\) is positive definite. Indeed, for any nonzero vector \(\xi\in\mathbb{R}^P\), \[\xi^TM_P\xi = \int_{\omega_\delta} \left|\sum_{j=1}^P \xi_j e_j(x,y)\right|^2\,dx\,dy>0,\] because a nontrivial finite linear combination of Dirichlet eigenfunctions cannot vanish identically on the open set \(\omega_\delta\). Therefore, \(M_P\) is invertible. Consequently, \[\operatorname{rank}\,[B_P,A_PB_P]=2P,\] and by the Kalman’s controllability rank condition (see, e.g., [39]), the finite-dimensional Galerkin system is controllable.

We compute \(M_P\) in closed form by exploiting orthonormality on the full square. Since \(\displaystyle\int_\Omega e_je_\ell = \delta_{j\ell}\) and \(\omega_\delta\) is the complement of the tensor-product square \([\delta,1-\delta]^2\), \[M_P = I_P - O_x\odot O_y, \qquad (O_x)_{j\ell} = \int_\delta^{1-\delta} 2\sin(m_j\pi x)\sin(m_\ell\pi x)\,dx,\] and similarly for \(O_y\) in the \(n\)-indices, with \(\odot\) denoting entry-wise multiplication. Numerically, the smallest eigenvalue of \(M_P\) decreases as \(P\) grows, reflecting the weaker interaction of some high modes with the collar. For \(P=8,16,32\), we find \[\lambda_{\min}(M_P) \approx 2.5\times 10^{-2},\quad 3.5\times 10^{-3},\quad 6.2\times 10^{-5},\] respectively. Thus, \(M_P\) remains positive definite over the truncation levels used in this benchmark, although the increasing ill-conditioning is visible.

The finite-dimensional Gramian reference control is computed with respect to the modal coefficient norm \[\int_0^T |c(t)|_{\mathbb{R}^P}^2\,dt.\] Let \[Y_0= \begin{pmatrix} a(0)\\ a'(0) \end{pmatrix}, \qquad Y_T=0,\] and define \[\label{eq:linear952d95gramian} W_T = \int_0^T e^{A_P(T-s)}B_PB_P^T e^{A_P^T(T-s)}\,ds.\tag{53}\] Since the pair \((A_P,B_P)\) is controllable, \(W_T\) is symmetric positive definite. The corresponding minimum-energy control steering \(Y_0\) to \(0\) is \[\label{eq:linear952d95hum95control} c_P^{\textrm{HUM}}(t) = B_P^T e^{A_P^T(T-t)}W_T^{-1} \left( -Y_{\textrm{free}}(T) \right), \qquad Y_{\textrm{free}}(T)=e^{A_PT}Y_0.\tag{54}\] Equivalently, \[Y(T) = e^{A_PT}Y_0 + \int_0^T e^{A_P(T-s)}B_Pc_P^{\textrm{HUM}}(s)\,ds = 0.\] In our implementation, \(W_T\) is assembled by Simpson quadrature applied to 53 , using matrix exponentials for the homogeneous flow. The computed reference control 54 steers the Galerkin system to the origin to within \[\left\|Y^{\textrm{HUM}}(T)\right\|_{\mathcal{E}_P}\sim 10^{-8},\] so it provides a numerically accurate reference for the projected system.

4.4.0.1 Neural-Spectral control.

We compare this reference with the Neural-Spectral control \(c_P^{\textrm{NS}}\). For the finite-dimensional linear system, the formula 54 implies that each component of \(c_P^{\textrm{HUM}}\) is a finite linear combination of temporal harmonics \[\cos(\sqrt{\mu_k}\,t), \qquad \sin(\sqrt{\mu_k}\,t), \qquad k=1,\ldots,P.\] This motivates the following frequency-informed temporal representation: \[\label{eq:linear952d95ns95ansatz} c_{\theta,j}(t) = \Theta_{j,0} + \sum_{k=1}^P \Bigl[ \Theta_{j,2k-1}\cos(\sqrt{\mu_k}\,t) + \Theta_{j,2k}\sin(\sqrt{\mu_k}\,t) \Bigr],\tag{55}\] with learnable weights \(\Theta\in\mathbb{R}^{P\times(2P+1)}\). This dictionary should be viewed as a frequency-informed specialization of the Neural-Spectral temporal representation for this linear benchmark: the modal spatial representation and the reduced control objective are the same, but the generic trainable temporal network is replaced by a fixed feature map matched to the known wave frequencies.

Unlike the fully unconstrained Neural-Spectral formulation in previous sections (where the PDE residual is penalized as a soft constraint), here we employ a hard-constrained approach. Given \(c_\theta\), the state \(Y_\theta\) is propagated from \(Y_0\) through 52 using the same matrix-exponential discretization used for the reference computation. Thus, the projected dynamics and the initial condition are enforced through the finite-dimensional propagation. The control is therefore obtained by solely minimizing the reduced objective \[\label{eq:linear952d95ns95loss} \mathcal{L}_{\textrm{NS}}(\theta) = \bigl\|Y_\theta(T)\bigr\|_{\mathcal{E}_P}^2 + \lambda_{\textrm{reg}}\,\|c_\theta\|_{L^2(0,T;\mathbb{R}^P)}^2 , \qquad \lambda_{\textrm{reg}}=3\times 10^{-4}.\tag{56}\] The first term penalizes the terminal mismatch, while the second term biases the optimization toward lower-energy controls. For small terminal error and small regularization parameter, this penalized formulation approximates the constrained minimum-energy Galerkin control problem solved explicitly by the Gramian formula. Optimization uses Adam followed by L-BFGS. After training, the recovered control is validated by solving 52 forward with \(c_P^{\textrm{NS}}\) using a high-order adaptive integrator.

We report the forward-validation terminal error \[\left\|Y^{\textrm{NS}}(T)\right\|_{\mathcal{E}_P} = \left( \sum_{j=1}^P \mu_j \left|a_j^{\textrm{NS}}(T)\right|^2 + \sum_{j=1}^P \left|\left(a_j^{\textrm{NS}}\right)'(T)\right|^2 \right)^{1/2},\] the distance to the finite-dimensional HUM reference control, \[\left\|c_P^{\textrm{NS}}-c_P^{\textrm{HUM}}\right\|_{L^2(0,T;\mathbb{R}^P)},\] and the relative control norm \[\dfrac{\left\|c_P^{\textrm{NS}}\right\|_{L^2(0,T;\mathbb{R}^P)}} {\left\|c_P^{\textrm{HUM}}\right\|_{L^2(0,T;\mathbb{R}^P)}}.\]

4.4.0.2 Discussion.

The results in Table 5 and Figure 5 show that the frequency-informed Neural-Spectral reconstruction is close to the classical finite-dimensional Gramian/HUM control in this linear reference benchmark. The relative control norm matches that of \(c_P^{\textrm{HUM}}\) to within roughly one part in a thousand across the reported truncation levels, and the modal control trajectories are visually close to the Gramian reference. The forward-validated terminal error is reduced by approximately three orders of magnitude relative to the free evolution \[\|Y_0\|_{\mathcal{E}_P}\approx 10.87,\] which is independent of \(P\) because the initial data is captured exactly at every truncation level.

Two structural points should be noted. First, this benchmark is favorable to the modal Fourier ansatz 55 . Since the finite-dimensional reference control is a linear combination of the temporal harmonics \[\left\{\cos\left(\sqrt{\mu_k}t\right),\sin\left(\sqrt{\mu_k}t\right):k=1,\ldots,P\right\},\] a feature map built from these frequencies has a small representation error for this linear problem. A generic time-network, by contrast, must learn the relevant oscillatory time scales during training. Second, the regularization in 56 plays an important selection role. The projected system is controllable, so many admissible controls may steer \(Y_0\) to the origin. The regularized objective biases the optimization toward lower-energy controls, and in this benchmark, the resulting control is close to the Gramian/HUM reference.

Unlike the PINN and DeepONet comparisons, this benchmark does not depend on another learned solver. The comparison is made against the finite-dimensional Gramian/HUM control of the same projected linear dynamics. The close agreement with \(c_P^{\textrm{HUM}}\), together with a small forward-validation terminal error, provides a useful consistency check for the Neural-Spectral control formulation in a setting where a classical reference solution is available.

4.4.0.3 Effect of the temporal ansatz.

The convergence study above uses the modal Fourier dictionary 55 for the temporal representation, whereas the general Neural-Spectral formulation in this paper represents each modal coefficient with a shallow, trainable time-network. To quantify the effect of using frequency-informed temporal features, we repeat the benchmark at fixed \(P=16\) with both ansatzes. Both variants are trained with the same reduced objective 56 , the same optimizer schedule, and the same finite-dimensional state propagation. For the tanh ansatz, each modal control coefficient is represented as \[c_{\theta,j}(t) = \sum_{r=1}^{m} A_{j,r}\tanh\bigl(W_{j,r}\,t+B_{j,r}\bigr), \qquad m=64,\] with all of \(A,W,B\) trained. The results are reported in Table 6 and Figure 6, where the relative control error is \[\dfrac{\left\|c_P^{\textrm{NS}}-c_P^{\textrm{HUM}}\right\|_{L^2}} {\left\|c_P^{\textrm{HUM}}\right\|_{L^2}}.\]

Both variants drive the projected state close to the target, as reflected by their comparable terminal errors. They differ, however, in the control selected by the optimization. In this experiment, the tanh time-network produces an admissible steering control whose trajectory differs substantially from the minimum-energy Gramian/HUM reference. The modal Fourier ansatz incorporates the frequencies \(\sqrt{\mu_k}\) appearing in the propagator \(e^{A_P^T(T-t)}\), and therefore has a much smaller representation error for the reference control in this linear setting. With the same reduced objective, the Fourier-feature version recovers a control much closer to \(c_P^{\textrm{HUM}}\), with a relative control error below \(10^{-2}\).

The Fourier dictionary should therefore be interpreted as a frequency-informed specialization of the Neural-Spectral strategy for the linear wave equation, where the modal time scales are known analytically. It uses the same modal representation of the dynamics and the same reduced control objective, but replaces the generic trainable tanh time-network by a fixed feature map matched to the spectral time scales of the problem.

5 1D Re-indexing of Multi-dimensional Spectrum↩︎

As established in Lemma 1, the spectrum of the operator \(A\) is naturally indexed by pairs \[(k,\sigma)\in(\mathbb{N}^*)^N\times\{+,-\}.\] For the purposes of constructing neural operator approximations, it is convenient to replace this multi-index notation with a one-dimensional enumeration. This avoids notational and dimensional ambiguities when the spectral coefficients are represented as finite input and output vectors. The purpose of this appendix is to prove Lemma 2 by constructing such an explicit enumeration.

We first enumerate \((\mathbb{N}^*)^N\). For each \(k=(k_1,\dots,k_N)\in(\mathbb{N}^*)^N\), define \[s(k):=\sum_{j=1}^N k_j^2\in\mathbb{N}^*.\] For each \(m\in\mathbb{N}^*\), set \[S_m:=\left\{k\in(\mathbb{N}^*)^N:\;s(k)=m\right\}.\] Each set \(S_m\) is finite. Indeed, if \(k\in S_m\), then \[1\le k_j\le \sqrt m, \qquad j=1,\ldots,N.\] Hence, each coordinate \(k_j\) can take at most \(\lfloor\sqrt m\rfloor\) values, and therefore \[|S_m|\le \lfloor\sqrt m\rfloor^N<\infty.\] Moreover, since \(s(k)\in\mathbb{N}^*\) for every \(k\in(\mathbb{N}^*)^N\), we have \(k\in S_{s(k)}\). Thus \[(\mathbb{N}^*)^N=\bigcup_{m=1}^\infty S_m.\] The union is disjoint, because if \(k\in S_m\cap S_\ell\), then \(s(k)=m=\ell\).

For each nonempty set \(S_m\), we order its elements lexicographically and write \[S_m=\left\{k_1^m,k_2^m,\ldots,k_{|S_m|}^m\right\}.\] Define \[n_0:=0, \qquad n_m:=\sum_{\ell=1}^m |S_\ell|, \qquad m\ge1.\] Then \((n_m)_{m\ge0}\) is nondecreasing, and \(n_m\to\infty\) as \(m\to\infty\). To see this, for each \(r\in\mathbb{N}^*\), consider \[k^r:=(r,1,\ldots,1)\in(\mathbb{N}^*)^N.\] These vectors are pairwise distinct and satisfy \[s(k^r)=r^2+(N-1).\] Thus \(S_{r^2+(N-1)}\neq\emptyset\) for infinitely many values of \(r\). Consequently, the partial sums \[n_m=\sum_{\ell=1}^m |S_\ell|\] cannot remain bounded, and hence \(n_m\to\infty\).

We now define a map \[q:\mathbb{N}^*\to(\mathbb{N}^*)^N.\] Given \(j\in\mathbb{N}^*\), since \(n_m\to\infty\), there exists at least one \(m\ge1\) such that \(j\le n_m\). Let \(m\) be the smallest such integer. Then \[n_{m-1}<j\le n_m.\] In particular, \(S_m\neq\emptyset\), and we define \[\label{def:q} q(j):=k^m_{\,j-n_{m-1}}.\tag{57}\] This gives a well-defined map \[q:\mathbb{N}^*\longrightarrow(\mathbb{N}^*)^N, \qquad j\longmapsto q(j).\]

We now prove that \(q\) is a bijection. Suppose first that \(q(j_1)=q(j_2)\). Let \(m_1,m_2\) be the corresponding integers such that \[n_{m_1-1}<j_1\le n_{m_1}, \qquad n_{m_2-1}<j_2\le n_{m_2}.\] Then \[q(j_1)\in S_{m_1}, \qquad q(j_2)\in S_{m_2}.\] Since the sets \(S_m\) are pairwise disjoint and \(q(j_1)=q(j_2)\), we must have \(m_1=m_2\). Within the ordered list of \(S_{m_1}\), the equality \(q(j_1)=q(j_2)\) then implies that the two elements have the same position, namely \[j_1-n_{m_1-1}=j_2-n_{m_1-1}.\] Hence \(j_1=j_2\), and \(q\) is injective.

To prove surjectivity, let \(k\in(\mathbb{N}^*)^N\). Set \(m:=s(k)\). Then \(k\in S_m\), so there exists an index \(r\), with \(1\le r\le |S_m|\), such that \[k=k_r^m.\] By construction, \[q(n_{m-1}+r)=k_r^m=k.\] Therefore \(q\) is surjective. Hence \(q\) is a bijection.

Finally, to incorporate the two branches \(\{+,-\}\), we define \[\kappa:\mathbb{N}^*\to(\mathbb{N}^*)^N\times\{+,-\}\] by \[\kappa(2j-1):=(q(j),+), \qquad \kappa(2j):=(q(j),-), \qquad j\in\mathbb{N}^*.\] Since \(q\) is a bijection, it follows immediately that \(\kappa\) is also a bijection. Indeed, every pair \((k,+)\) is attained by \(2q^{-1}(k)-1\), and every pair \((k,-)\) is attained by \(2q^{-1}(k)\).

Using \(\kappa\), we may re-index the eigenvalues and eigenfunctions of \(A\) by a single integer. For each \(n \in \mathbb{N}^*\), we define \[\begin{align} \lambda_{2n-1}:=\lambda_{q(n)}^{+}, \qquad \lambda_{2n}:=\lambda_{q(n)}^{-} \\ \phi_{2n-1}:=\phi_{q(n)}^{+} \qquad \phi_{2n}:=\phi_{q(n)}^{-}. \end{align}\] Thus, the full spectral family \[\left\{\left(\lambda_k^\pm,\phi_k^\pm\right)\right\}_{k\in(\mathbb{N}^*)^N}\] is equivalently represented as the one-dimensional sequence \[\{(\lambda_n,\phi_n)\}_{n\in\mathbb{N}^*}.\]

References↩︎