1 Introduction↩︎

Partial differential equations (PDEs) arise in many scientific branches. In some cases problems, the PDE can be derived from a functional associated to a variational problem. In these cases, the solutions of the PDE problem can be recognized as critical points of the functional [1]. This forms the foundation of the calculus of variations [2]-[3], and its solutions involve the well-known Euler-Lagrange equations, which allow for a transition to the PDE approach. This approach is frequently found in physics [4]. In fact, there are many versions of this approach depending on the different fields of physics (classical mechanics [5],[6]; quantum mechanics [7],[8]; particle physics [9],[10]; general relativity [11],[12], etc.). These versions of the approach are known as Action Principles and constitute the core of the different fields of physics. An Action Principle starts with a scalar function called a Lagrangian, which describes the physical system. The accumulated value of this Lagrangian function between two states of the system is called the action. The Action Principles assume that, among all possibilities, the actual evolution of the physical system produces a critical point of the action.

When exact solutions are not possible, approximation techniques are usually employed, such as the Finite Element Method (FEM) [13], [14], Galerkin Method [15], and Ritz Method [16], [17]. The Ritz Method is an interesting approximation technique which consists of choosing a finite number of admissible trial functions such that the solution \(y(x)\) can be determined as a linear combination of these finite trial functions. Recently, novel versions of the Ritz Method integrated with deep learning have been explored. This is known as the Deep Ritz Method [18]-[19]. This method is based on the Universal Approximation Theorem, which essentially states that a neural network with a hidden layer is dense in subsets of a Sobolev space [20]-[21]. However, the boundary conditions of the variational problem are not satisfied by the neural network. To address this non-fulfillment of the boundary conditions, the Deep Ritz Method redefines the action by adding a penalty term.

On one hand, the added penalty term can cause the action to become a non-convex function. This is a drawback, as finding local minima and saddle points for the action is a non-trivial task. On the other hand, this penalty term is weighted by a hyperparameter that must be adjusted to assign more or less importance to the boundary conditions in the optimization process. This adjustment makes the convergence rate strongly dependent on this hyperparameter.

In recent works [22], [23], several ansatzes have been proposed to enforce the boundary conditions of admissible functions in PINNs (Physics-Informed Neural Networks) and VPINNs (Variational Physics-Informed Neural Networks), thereby avoiding the need for penalty terms. However, these works only provide numerical experiments and do not include theoretical results supporting the proposed ansatzes. In this work, we revisit the idea of using ansatzes that exactly satisfy the boundary conditions, but we introduce a simpler construction together with the corresponding theoretical results that justify the proposed approach in Sobolev spaces. This ansatz consists of one term that enforces the boundary behaviour and another term given by the product of a neural network and a function \(g\), chosen so that their ratio remains differentiable inside the domain. The function \(g\) can be adjusted to vanish at the boundary of the domain, ensuring that the boundary conditions are satisfied. For the rest of the domain, the behavior of the ansatz is dominated by the product of the neural network and the function \(g\). It is important to point out that, since the Lagrangian usually depends of the admisible functions and their derivatives, Sobolev spaces are the natural framework to demonstrate the theorems that support the proposed ansatzes.

Throughout this work, we demonstrate that the product of a polynomial and a single-hidden-layer neural network is dense in the relevant Sobolev spaces, provided that the polynomial is chosen so that their ratio remains differentiable inside the domain. This property allows us to consider the proposed ansatz as an admissible function to approximate the solution of the variational problem using the Ritz method approach. It is important to point out that the conditions ensuring the validity of using this type of ansatz had not been theoretically established before. Moreover, we present numerical experiments showing that the proposed ansatz significantly reduces the complexity of the neural network and prevents misleading outcomes in the optimization of the action functional (for instance, obtaining lower action values at the cost of not satisfying the boundary conditions exactly).

The remainder of the paper is organized as follows: The problem definition, theorems, and lemmas required are presented in Section 2. The ansatz and proposed method are described in Section 3. Several examples and comparative results are shown to compare the performance of the proposed method with the well-established Deep Ritz Method. Finally, concluding remarks are presented in Section 7.

2 Theoretical Foundations: Universal Approximation with Polynomial Factors↩︎

Formally, a Lagrangian function is defined as follows:

\[\begin{align} \label{eq1} \mathcal{L}:\left[\mathcal{C}^1(\Omega)\right]\times\mathbb{R} \rightarrow\mathbb{R} \end{align}\tag{1}\]

Where we assume that \(\Omega\subset\mathbb{R}^n\) is compact. The action is defined as the accumulated value of the Lagrangian over a domain \(\Omega\):

\[\begin{align} \label{eq2} S=\int_{\Omega}{\mathcal{L}\left(y(x),\frac{\partial y(x)}{\partial x},x\right)dx} \end{align}\tag{2}\]

where \(y(x)\in \mathcal{C}^1(\Omega)\) belongs to the set of admissible functions such that \(y(x)=b(x)\) \(\forall x\in\partial\Omega\subset\mathbb{R}^{n}\) and \(\left(\frac{\partial y(x)}{\partial x}\right)_{j}=\frac{\partial y}{\partial x_j}\) with \(j=1,\ldots,n\). The action principles assume that the real evolution \(y(x)\) produces a stationary critical point of the action \(S\) (i.e., \(\delta S=0\), where \(\delta S\) denotes the functional derivative). This constitutes the class of variational problems, and their solutions involve the well-known Euler-Lagrange equations, which allow for a transition to the PDE approach [1].

When exact solutions are not possible, approximation techniques are usually employed. One interesting approximation technique is the Ritz Method [13]. Recently, a novel version of the Ritz Method integrated with deep learning has been explored, known as the Deep Ritz Method [18]. This method is based on the Universal Approximation Theorem [20] and its extension [24]. In [20], it is established that a single-hidden-layer neural network is dense in the space \(\mathcal{C}(\Omega)\) of continuous functions. The work [24] provides a generalization, showing that single-hidden-layer neural networks with differentiable activation functions are dense in \(\mathcal{C}^{m}(\Omega)\) for compact domains \(\Omega\subset\mathbb{R}^{n}\). Consequently, this extended version of the theorem can also be formulated in terms of the Sobolev norms [25], [26]:
Definition 1. Let \(m,\mathrm{p}\in\mathbb{N}\) and \(1\le p<\infty\). Let \(f\in\mathcal{C}^m(\Omega)\). Then, Sobolev norms are defined as: \[\begin{align} \label{eq3} \|f\|_{m,\mathrm{p},\Omega} = \left( \sum_{|\alpha|\le m} \|D^\alpha f\|_{L^\mathrm{p}(\Omega)}^{\mathrm{p}} \right)^{1/\mathrm{p}} \end{align}\tag{3}\] \[\begin{align} &\label{eq3b} \|f\|_{m,\infty,\Omega} = \max_{|\alpha|\le m} \|D^\alpha f\|_{L^\infty(\Omega)} \end{align}\tag{4}\] where \(D^\alpha\) denotes multi-index notation for partial derivatives and \(||\cdot||_{L^\mathrm{p}(\Omega)}\) is \(L^\mathrm{p}\)-norm in \(\Omega\).
According to [24] the universal approximation theorem can be stated by using Sobolev norms as follows:
Universal Approximation Theorem for Sobolev Spaces [24]. Let \(m\in\mathbb{N}\) and \(1\le p<\infty\). Let \(\sigma\in\mathcal{C}^m(\mathbb{R})\) be a non-constant and bounded function, and let \(f\in \mathcal{C}^m(\Omega)\), where \(\Omega\subset\mathbb{R}^n\) is compact. Then, for any \(\varepsilon>0\), there exist \(N\in\mathbb{N}\), \(w_l\in\mathbb{R}^n\), \(b_l\in\mathbb{R}\), \(\alpha_l\in\mathbb{R}\), \(l=1,\ldots,N\), such that \[\begin{align} \label{eq4} \bigg\|f(x)-\sum_{l=1}^{N}\alpha_l \sigma(w_l^T x +b_l)\bigg\|_{m,p,\Omega}<\varepsilon, \end{align}\tag{5}\] where \(A^T\) denotes the transpose of \(A\).
According to [24], the Universal Approximation Theorem also holds in Sobolev spaces, where the approximation error is measured with the norm \(\|\cdot\|_{m,\mathrm{p},\Omega}\) instead of the supremum norm \(\|\cdot\|_{\infty}\) as it is stated in [20]. This point is important because, for variational problems with Dirichlet boundary conditions, the norm \(\|\cdot\|_{m,\mathrm{p},\Omega}\) is the natural norm to measure approximation errors. In particular, since the Lagrangian function depends of \(y(x)\) and \(\frac{\partial y(x)}{\partial x}\), the approximation \(y(x)\approx \sum_{l=1}^{N}{\alpha_l \sigma(w_l^T x +b_l)}\) made in the Deep Ritz Method should be measured by the norm \(\|y(x)-\sum_{l=1}^{N}{\alpha_l \sigma(w_l^T x +b_l)}\|_{1,\mathrm{p},\Omega}\) for a fixed value of \(N\), and searches for \(w_l\in\mathbb{R}^n\), \(b_l\in \mathbb{R}\), and \(\alpha_l\in\mathbb{R}\) such that \(\delta S=0\) should be done using a suitable optimization algorithm.
Note that the boundary condition \(y(x)=b(x)\) \(\forall x\in\partial\Omega\subset\mathbb{R}^{n}\) is not guaranteed, since the parameters of the neural network \(w_l\in\mathbb{R}^n\), \(b_l\in \mathbb{R}\), and \(\alpha_l\in\mathbb{R}\) are only optimized to ensure \(\delta S=0\). To address the absence of the boundary condition, the Deep Ritz Method redefines the action using a penalty method as \[\begin{align} \label{eq5} S_{\mathrm{DR}}=\int_{\Omega}{\mathcal{L}\left(y(x),\frac{\partial y(x)}{\partial x},x\right)dx}+\beta\int_{\partial\Omega}{(y(x)-b(x))^2dx} \end{align}\tag{6}\] where \(\beta\) is an adjustable parameter that weights the importance of the boundary condition in the optimization process to achieve \(\delta S_{\mathrm{DR}}=0\). The main drawback of the penalty method is that the boundary condition is not exactly satisfied but only approximated.

In order to ensure the exact satisfaction of the boundary conditions, we propose an ansatz that includes the product of a single-hidden-layer neural network and a suitably chosen function. With this construction, we present the following theorem, which extends the Universal Approximation Theorem to products of functions and single-hidden-layer neural networks:
Theorem 1. Let \(m\in\mathbb{N}\) and \(1\le p<\infty\). Let \(\sigma\in\mathcal{C}^m(\mathbb{R})\) be a non-constant and bounded function and \(f,g\in\mathcal{C}^m(\Omega)\), such that \(f/g\in\mathcal{C}^m(\Omega)\) where \(\Omega\subset\mathbb{R}^n\) is compact. Then, for any \(\varepsilon'>0\), there exist \(N\in\mathbb{N}\), \(w_l\in\mathbb{R}^n\), \(b_l\in \mathbb{R}\), and \(\alpha_l\in\mathbb{R}\) with \(l=1,\ldots,N\) such that \[\begin{align} \label{eq6} &\bigg\|f(x)-g(x)\sum_{l=1}^{N}{\alpha_l \sigma(w_l^T x +b_l)}\bigg\|_{m,\mathrm{p},\Omega}<\varepsilon' \end{align}\tag{7}\]

Remark 1. Note that 26 is a Hölder-type product estimate in Sobolev spaces
In other words, the product of a function \(g\in\mathcal{C}^m(\Omega)\) and a single-hidden-layer neural network is dense in the set \(\{\,f\in\mathcal{C}^m(\Omega) : f/g \in \mathcal{C}^{m}(\Omega)\,\}\), where \(\Omega \subset \mathbb{R}^{n}\) is a compact set.
Let us particularize Theorem 1 to the case of variational problems with boundary conditions, imposing that the function \(f\), and therefore \(g\), vanish at the boundary \(\partial\Omega\).
Corollary 1. Let \(m\in\mathbb{N}\) and \(1\le p<\infty\). Let \(\sigma\in\mathcal{C}^m(\mathbb{R})\) be a non-constant and bounded function. Let \(\Omega\subset\mathbb{R}^n\) be compact, and let \(f,g\in \mathcal{C}^m(\Omega)\) be such that \[\begin{align} \label{eq:boundary95zero} f(x)=0,\quad g(x)=0 \quad \forall x\in\partial\Omega, \end{align}\tag{29}\] and assume that \[\begin{align} \label{eq:ratio95condition} \frac{f}{g}\in \mathcal{C}^m(\Omega). \end{align}\tag{30}\] Then, for any \(\varepsilon'>0\), there exist \(N\in\mathbb{N}\), \(w_l\in\mathbb{R}^n\), \(b_l\in\mathbb{R}\), and \(\alpha_l\in\mathbb{R}\), \(l=1,\ldots,N\), such that \[\begin{align} \label{eq:corollary95result} \bigg\|f(x)-g(x)\sum_{l=1}^{N}{\alpha_l \sigma(w_l^T x +b_l)}\bigg\|_{m,\mathrm{p},\Omega}<\varepsilon'. \end{align}\tag{31}\]
This corollary requires that, if \(f(x)=0\) for all \(x\in\partial\Omega\), then it is also required that \(g(x)=0\) for all \(x\in\partial\Omega\), and furthermore, \(f\) must vanish at least at the same rate as \(g\) in \(\partial\Omega\) in order to satisfy Theorem 1.

To enforce the boundary condition \(y(x)=b(x)\) for all \(x\in\partial\Omega\), with \(b\in \mathcal{C}^{1}(\Omega)\), Corollary 1 suggests the following approximation: \[\begin{align} \label{eq14} &y(x)\approx B(x)+g(x)\sum_{l=1}^{N}\alpha_l\,\sigma(w_l^{T}x+b_l), \end{align}\tag{32}\] where \(g(x)=0\) and \(B(x)=b(x)\) for all \(x\in\partial\Omega\). For simplicity, we choose \(g(x)\) to be a polynomial \(p(x)\), since polynomial functions can be easily constructed so that \(p(x)=0\) for all \(x\in\partial\Omega\), although other choices are possible. Note that, since \(b\in\mathcal{C}^{1}(\Omega)\) and \(g(x)=p(x)\) is a polynomial, choosing \(B\in \mathcal{C}^{1}(\Omega)\) ensures that the proposed ansatz is an admissible function for the problem given in 2 .

Although all theorems and corollaries are stated for \(\Omega\) being a generic compact set, we assume that \(\Omega\) is (piecewise) diffeomorphic to \([0,1]^{n}\) or to the unit ball \(B_1^n\subset\mathbb{R}^n\) (i.e., there exists a piecewise smooth bijection between \(\Omega\subset\mathbb{R}^{n}\) and either \([0,1]^{n}\) or \(B_1^n\)). This assumption is motivated by the fact that common differentiable activation functions exhibit significant variation when their argument lies in \([0,1]\), while they tend to saturate outside this interval, as occurs with the sigmoid or the hyperbolic tangent. Moreover, several strategies are available to construct such a diffeomorphism (see Appendix C in [1]).

In the remainder of this work, and for simplicity, we assume that \(\Omega\) is (piecewise) diffeomorphic to \([0,1]^{n}\). This assumption facilitates the choice of a polynomial \(p(x)\) satisfying \(p(x)=0\) for all \(x\in\partial[0,1]^{n}\), such as \[\begin{align} \label{eq18} p(x) = \prod_{j=1}^{n} x_j (1 - x_j). \end{align}\tag{33}\] Note that this choice of polynomial ensures that \(p(x)\to 0\) linearly as \(x_j\to 0\) or \(x_j\to 1\) (i.e., \(p(x)\) vanishes linearly as \(x\) approaches the \(j\)-th face of the hypercube). Since \(p(x)\) vanishes linearly near \(\partial[0,1]^{n}\), the solution of the variational problem \(y(x)\) must verify \((y(x)-b(x))\rightarrow 0\) at least at the same rate as \(p(x)\) when approaching the boundary, as required to satisfy Theorem 1. This last condition is fulfilled, since \(y(x)\) and \(b(x)\) belong to \(\mathcal{C}^1(\Omega)\) and therefore satisfy Taylor theorem.

The case of \(B_1^n\) could be easily handled choosing \(g(x)=1-\|x\|^2_2\), but this case will not be considered in this work.

3 Ansatz Construction and Variational Minimization Procedure↩︎

Let us consider the differentiable transformation \(T:[0,1]^n\rightarrow\Omega\) with \(\Omega\subset\mathbb{R}^n\) such that if \(u\in\partial[0,1]^{n}\) then \(T(u)=x\in\partial\Omega\), where \(x=(x_1,\ldots,x_n)\) and \(u=(u_1,\ldots,u_n)\), such that \(x_j=T_j(u_1,\ldots,u_n)\) with \(T_j:[0,1]^n\rightarrow\mathbb{R}\). Applying this transformation to 2 leads to the following result: \[\begin{align} \label{eq15} \hat{S}=\int_{[0,1]^n}{\mathcal{L}\left(\hat{y}(u),(J_T(u)^{-1})^T\frac{\partial \hat{y}(u)}{\partial u},T(u)\right)|\mathrm{det}(J_{T}(u))|du} \end{align}\tag{34}\] where \(\left(J_{T}(u)\right)_{ij}=\frac{\partial T_i}{\partial u_j}\). In 34 , the expressions \(y(x)=\hat{y}(u)\), \(\frac{\partial y(x)}{\partial x}=(J_T(u)^{-1})^T\frac{\partial \hat{y}(u)}{\partial u}\) and \(x=T(u)\) have been taken into account. For \(y(x)\in \mathcal{C}^1\) to be an admissible function, it must satisfy \(y(x)=b(x)\) \(\forall x\in\partial\Omega\). Therefore, the \(T\) differentiable transformation has to ensure \(\hat{y}(u)=b(T(u))\) \(\forall u\in\partial[0,1]^{n}\). Let us denote \(\hat{b}(u)=b(T(u))\). According to the previous section, we propose the following ansatz: \[\begin{align} \label{eq16} &\hat{y}(u,\theta)= B(u)+p(u)N_{\mathrm{net}}(u,\theta) \end{align}\tag{35}\] where \(\theta=(\alpha_{1},\ldots,\alpha_{N},w_{1},\ldots,w_{N},b_{1},\ldots,b_{N})\), and \(N_{\mathrm{net}}(u,\theta)\) represents single-hidden-layer neural networks with \(N\) neurons: \[\begin{align} \label{eq17} &N_{\mathrm{net}}(u,\theta)= \sum_{l=1}^{N}{\alpha_{l} \tanh\left(w_{l}^T u +b_{l}\right)} \end{align}\tag{36}\] where \(N\in\mathbb{N}\), \(w_{l}\in\mathbb{R}^n\), \(b_{l}\in \mathbb{R}\), \(\alpha_{l}\in\mathbb{R}\), the polynomial \(p(u)\) is defined as 33 and where the term \(B(u)\) satisfies \(B(u)=\hat{b}(u)\) \(\forall u\in\partial[0,1]^n\). This can be achieved through transfinite interpolation [29],[30], following the inclusion-exclusion principle [31]. For instance, in the one-dimensional case, given \(\hat{b}(0)\) and \(\hat{b}(1)\), this leads to: \[\begin{align} \label{eq19} &B(u)=(1-u)\hat{b}(0)+u\hat{b}(1) \end{align}\tag{37}\] and for the bidimensional case (also known as the Coons patch [32]), given \(\hat{b}(0,u_2)\), \(\hat{b}(1,u_2)\), \(\hat{b}(u_1,0)\), \(\hat{b}(u_1,1)\), \(\hat{b}(0,0)\), \(\hat{b}(0,1)\), \(\hat{b}(1,0)\), and \(\hat{b}(1,1)\), this leads to: \[\begin{align} \nonumber &B(u)=\\ \nonumber &\left[(1-u_1)\hat{b}(0,u_2)+u_1\hat{b}(1,u_2)+(1-u_2)\hat{b}(u_1,0)+u_2\hat{b}(u_1,1)\right]-\\ \label{eq20} &\left[(1-u_1)(1-u_2)\hat{b}(0,0)+u_1(1-u_2)\hat{b}(1,0)+\right.\left.(1-u_1)u_2\hat{b}(0,1)+u_1u_2\hat{b}(1,1)\right] \end{align}\tag{38}\]

In what follows, for simplicity, we will consider a convex Lagrangian with a convex domain. Under this consideration, the condition of null variation of the action (i.e., \(\delta \hat{S}=0\)) is equivalent to minimizing the action \(\hat{S}\) [33],[3]. Therefore, we will aim to solve the following problem: \[\begin{align} \label{eq21} &\min_{\hat{y}(u)}\left\{\hat{S}(\hat{y}(u)) : \hat{y}(u)=\hat{b}(u),\;\forall u\in\partial[0,1]^{n}\right\} \end{align}\tag{39}\]

If we use the proposed ansatz given by 35 , the equivalent problem given by 39 becomes a minimization problem with respect to the parameters of the neural network \(\theta\). The advantage of this approach over the Deep Ritz method is that penalization is not necessary for the boundary condition, as it is inherently satisfied by the proposed ansatz (i.e., the proposed ansatz given by 35 constitutes a set of admissible functions).

The minimization of the action will be carried out using the gradient descent method, where the parameters \(\theta\) are updated iteratively as follows: \[\begin{align} \label{eq22} &\theta^{k+1}=\theta^{k}-\eta \frac{\partial\hat{S}(\theta^k)}{\partial \theta} \end{align}\tag{40}\] where \(\theta^k=(\theta_1^k,\ldots,\theta_{3N}^k)\) contains the neural network parameters \(\alpha_{1},\ldots,\alpha_{N},w_{1},\ldots,w_{N},b_{1},\ldots,b_{N}\) for the \(k\)th iteration, \(\eta\in\mathbb{R}\) is a learning rate (along results section, \(\eta=0.001\) will be used), and \(\frac{\partial\hat{S}}{\partial \theta}=(\frac{\partial\hat{S}}{\partial \theta_1},\ldots,\frac{\partial\hat{S}}{\partial \theta_{3N}})\) where \[\begin{align} \label{eq23} \frac{\partial\hat{S}}{\partial \theta}=\int_{[0,1]^{n}}{\left[\frac{\partial\mathcal{L}}{\partial \hat{y}}\cdot\frac{\partial\hat{y}}{\partial \theta}+\frac{\partial\mathcal{L}}{\partial \left(\frac{\partial\hat{y}}{\partial u}\right)}\cdot\frac{\partial^2\hat{y}}{\partial \theta \partial u}\right]|\mathrm{det}(J_{T}(u))|du} \end{align}\tag{41}\] where \(\left(\frac{\partial\hat{y}}{\partial \theta}\right)_{s}=\frac{\partial\hat{y}}{\partial \theta_s}\), \(\left(\frac{\partial\mathcal{L}}{\partial \left(\frac{\partial\hat{y}}{\partial u}\right)}\right)_{j}=\frac{\partial\mathcal{L}}{\partial \left(\frac{\partial\hat{y}}{\partial u_j}\right)}\), \(\left(\frac{\partial^2\hat{y}}{\partial \theta \partial u}\right)_{sj}=\frac{\partial^2\hat{y}}{\partial \theta_s \partial u_j}\). The derivatives of \(\hat{y}\) are calculated from the proposed ansatz given by 35 as follows: \[\begin{align} \tag{42} &\frac{\partial \hat{y}}{\partial \theta_s}=p(u)\frac{\partial N_{\mathrm{net}}}{\partial\theta_s}\\ \tag{43} &\frac{\partial^2\hat{y}}{\partial \theta_s \partial u_j}=\frac{\partial p}{\partial u_j}\frac{\partial N_{\mathrm{net}}}{\partial \theta_s}+p(u)\frac{\partial^2 N_{\mathrm{net}}}{\partial \theta_s\partial u_j} \end{align}\] where \[\begin{align} \label{eq26} &\frac{\partial p}{\partial u_j}=(1-2u_j)\prod_{i\neq j}^n u_i(1-u_i) \end{align}\tag{44}\]

The computation of the integral given by 34 and 41 is approximated using conventional quadrature rules, such as Gauss-Legendre, or specialized quadrature rules for cases involving singular integrands.

4 Results↩︎

In this section, we present some examples where the results obtained using the proposed method are compared with those achieved by the Deep Ritz method, as shown in [18]. All results were obtained using custom-developed code written in Python and PyTorch [34], leveraging the efficient neural network library to access the required evaluation functions and derivatives.

4.1 Example 1↩︎

In this section, we demonstrate the application of the proposed method for a simple action, in order to illustrate didactically how to proceed. In this example, we aim to solve the following problem: \[\begin{align} \label{eq27} \min_{y(x)}\left\{ S=\int_{a}^{b}{\left[y(x)+\left(\frac{dy}{dx}\right)^2\right]dx} : y(a)=A,\;y(b)=B\right\} \end{align}\tag{45}\] where \(x\in[a,b]\) and \(y(x)\in\mathbb{R}\). Let us consider the linear transformation \(T:[0,1]\rightarrow [a,b]\) defined as \(x=(b-a)u+a\), where \(u\in[0,1]\). Applying this linear transformation to the problem defined in 45 leads to: \[\begin{align} \label{eq28} \min_{\hat{y}(u)}\left\{ \hat{S}=\int_{0}^{1}{\left[(b-a)\hat{y}(u)+\frac{1}{b-a}\left(\frac{d\hat{y}}{du}\right)^2\right]du} : \hat{y}(0)=A,\;\hat{y}(1)=B\right\} \end{align}\tag{46}\] In order to minimize the action with a set of admissible functions, we propose the following ansatz, in accordance with 35 , and following 37 : \[\begin{align} \label{eq29} &\hat{y}(u,\theta)=A(1-u)+Bu+u(1-u)N_{\mathrm{net}}(u,\theta) \end{align}\tag{47}\] where \(N_{\mathrm{net}}(u,\theta)\) is defined as in 36 . Note that \(\hat{y}(0,\theta)=A\) and \(\hat{y}(1,\theta)=B\), using the proposed ansatz in 47 (i.e., the set of \(\hat{y}(u,\theta)\in \mathcal{C}^1\) provided by 47 constitutes a set of admissible functions). According to [18], the minimization problem via the Deep Ritz method is reduced to minimizing the action given by: \[\begin{align} \nonumber &\hat{S}_{\mathrm{DR}}=\int_{0}^{1}{\left[(b-a)N_{\mathrm{net}}(u,\theta)+\frac{1}{b-a}\left(\frac{\partial N_{\mathrm{net}}(u,\theta)}{\partial u}\right)^2\right]du}\\ \label{eq30} &+\beta\left[\left(N_{\mathrm{net}}(0,\theta)-A\right)^2+\left(N_{\mathrm{net}}(1,\theta)-B\right)^2\right] \end{align}\tag{48}\]

In the following, we will consider \(a=A=0\) and \(b=B=10\). Since the Lagrangian is convex with a convex domain, the exact solution can be obtained by solving the Euler-Lagrange equations. This exact solution, for these values of \(a\), \(b\), \(A\), and \(B\), is given by: \[\begin{align} \label{eq31} \hat{y}_{\mathrm{exact}}(u)=25u^2-15u \end{align}\tag{49}\]

Figure 1 shows the values of \(\hat{S}\) with respect to the iteration number of the gradient descent method. Three cases have been considered: case \(N=2\) when \(\hat{S}\) is minimized using the ansatz given by 47 ; case \(N=2\) when \(\hat{S}_{\mathrm{DR}}\) is minimized; and case \(N=30\) when \(\hat{S}_{\mathrm{DR}}\) is minimized. The value of \(\beta=50\) has been used for the minimization of \(\hat{S}_{\mathrm{DR}}\). With this value of \(\beta\), the boundary values of the \(N_{\mathrm{net}}\) reached in the last iteration are \(N_{\mathrm{net}}(0,\theta)=-0.04\) and \(N_{\mathrm{net}}(1,\theta)=9.94\) for \(N=2\), and \(N_{\mathrm{net}}(0,\theta)=-0.03\) and \(N_{\mathrm{net}}(1,\theta)=9.93\) for \(N=30\). The initial values of the \(\theta\) parameters were randomly sampled from a Gaussian distribution with a mean of 10 and a standard deviation of 0.1, and these initial values were identical for all cases considered.

The obtained results in Figure 1 show that faster convergence is achieved when \(\hat{S}\), using the ansatz given in 47 , is minimized. Specifically, four significant digits are reached at iteration 4200. In comparison, using the Deep Ritz method, four significant digits are reached at iterations 29200 and 29400 for \(N=2\) and \(N=30\), respectively.

It is worth noting that the values of \(\hat{S}\) obtained with the Deep Ritz method for \(N=30\) are the lowest at the final iteration. This might suggest that the solution closest to the exact solution is provided by the Deep Ritz method with \(N=30\). However, as shown in subsequent figures, the solution closest to the exact solution is actually provided by the proposed method.

This observation may appear counterintuitive, since the Deep Ritz method with \(N=30\) achieves the lowest value of \(\hat{S}\). Nevertheless, it is important to consider that the boundary values obtained using the Deep Ritz method with \(N=30\) are close to the exact values (\(N_{\mathrm{net}}(0,\theta)=-0.03\) and \(N_{\mathrm{net}}(1,\theta)=9.93\)), but not exact. In fact, with these boundary values, the lowest values of \(\hat{S}\) achieved by the Deep Ritz method with \(N=30\) are very close to the values of \(\hat{S}\) obtained by the proposed method at the final iteration. This indicates that a failure to exactly satisfy the boundary conditions can lead to this type of misleading result.

Figure [fig2](a) illustrates the \(\hat{y}(u)\) values with respect to \(u\in[0,1]\) for the \(\theta\) values obtained in the final iteration. Figure [fig2](b) displays the absolute error relative to the exact solution given by 49 . The three cases are once again considered: the case of the proposed method using the ansatz given in 47 , the case of the Deep Ritz method with \(N=2\), and the case of the Deep Ritz method with \(N=30\).

Although the boundary condition is approximately satisfied by the solution obtained using the Deep Ritz method for \(N=2\), achieving \(N=30\) is necessary to approximate the solution to the exact solution (i.e., the differences between the Deep Ritz solution and the exact solution are unacceptable when \(N=2\)). In fact, the absolute error achieved by the Deep Ritz method is approximately 0.01 when \(N=30\), whereas the error obtained by the proposed method is less than 0.0001 when \(N=2\). These results validate the use of the ansatz, as the absolute errors of the solution obtained are 100 times smaller than those of the Deep Ritz method, despite employing a very simple neural network with a single hidden layer and \(N=2\).

4.2 Example 2↩︎

In this section, we aim to solve the following problem: \[\begin{align} \label{eq32} &\min_{y(x_1,x_2)}\left\{\begin{matrix} S=\int_{\Omega}{\left[||\nabla y(x_1,x_2)||^2-f(x_1,x_2)y(x_1,x_2)\right]dx}\\ \\ y(x_1,x_2)=b(x_1,x_2), \forall (x_1,x_2)\in\partial\Omega \end{matrix} \right\} \end{align}\tag{50}\] where \(\Omega=[0,1]\times[0,1]\), the norm \(||\cdot||\) is the Euclidean norm, and \(f(x_1,x_2)\) is defined as follows: \[\begin{align} \label{eq33} &f(x_1,x_2)=e^{-x_1}(2-x_1-x_2^3-6x_2) \end{align}\tag{51}\] and where: \[\begin{align} \label{eq34} &b(x_1,x_2)=\left\{\begin{matrix}x_2^3\mathrm{if}x_1=0\\ e^{-1}(1+x_2^3)\mathrm{if}x_1=1\\x_1e^{-x_1}\mathrm{if}x_2=0\\ (x_1+1)e^{-x_1}\mathrm{if}x_2=1\\ 0\mathrm{otherwise} \end{matrix} \right.\; \end{align}\tag{52}\] The exact solution is provided as follows: \[\begin{align} \label{eq35} &y_{\mathrm{exact}}(x_1,x_2)=e^{-x_1}\left(x_1+x_2^3\right) \end{align}\tag{53}\] In order to minimize the action with a set of admissible functions, we propose the following ansatz, in accordance with 35 : \[\begin{align} \label{eq36} &y(x_1,x_2,\theta)=B(x_1,x_2)+x_1(1-x_1)x_2(1-x_2)N_{\mathrm{net}}(x_1,x_2,\theta) \end{align}\tag{54}\] where \(B(x_1,x_2)\) is directly derived from 38 by substituting \(x_1,x_2\) in place of \(u_1,u_2\). It is important to note that, in this case, we work directly with variables \(x_1\) and \(x_2\) to simplify the \(T\) transformation described in the previous section. According to [18], the minimization problem in the Deep Ritz method is reduced to minimizing the action given by: \[\begin{align} \nonumber &S_{\mathrm{DR}}=\int_{\Omega}{\left[||\nabla N_{\mathrm{net}}(x_1,x_2,\theta)||^2-f(x_1,x_2)N_{\mathrm{net}}(x_1,x_2,\theta)\right]dx}\\ \label{eq37} &+\beta\int_{\partial\Omega}{\left[N_{\mathrm{net}}(x_1,x_2,\theta)-b(x_1,x_2)\right]^2dx} \end{align}\tag{55}\]

Figure 3 illustrates the values of \(S\) with respect to the iteration number of the gradient descent method. Two cases have been considered: the case of \(N=5\), where \(S\) is minimized using the ansatz given by 54 , and the case of \(N=10\), with an additional hidden layer, where \(S_{\mathrm{DR}}\) is minimized. This additional hidden layer was necessary to achieve acceptable results in the minimization of \(S_{\mathrm{DR}}\). A value of \(\beta=30\) was used during the minimization of \(S_{\mathrm{DR}}\). The initial values of the \(\theta\) parameters were randomly sampled from a Gaussian distribution with a mean of 1 and a standard deviation of 0.1, with these initial values being identical across all considered cases.

The obtained results in Figure 3 show that faster convergence is achieved when \(S\), using the ansatz given in 54 , is minimized. Specifically, four significant digits are reached at iteration 2800, whereas with the Deep Ritz method, four significant digits are reached at iteration 9800.

It is worth noting that the lowest value of \(S\) is achieved by the Deep Ritz method at approximately 2000 iterations. However, this is accomplished at the expense of deviating the solution’s boundary values from those required by the boundary condition in the optimization process. This demonstrates another advantage of using an ansatz that exactly satisfies the boundary conditions, as it prevents such misleading outcomes.

We would like to point out that, as observed in the figure 3, the value of \(S\) obtained with the Deep Ritz method is not monotonically decreasing with respect to the iteration number, whereas the value of \(S\) produced by the proposed method is. This behavior is expected, since the Deep Ritz method does not minimize \(S\) directly but rather the functional \(S_{\mathrm{DR}}\) defined in 55 , which includes the penalization term enforcing the boundary conditions. In contrast, the functional \(S_{\mathrm{DR}}\) itself is monotonically decreasing along the iterations, as expected from the application of gradient descent.

Figure [fig4](a) illustrates the \(y(x_1,x_2)\) values with respect to \(x_1\) and \(x_2\) for the \(\theta\) values obtained in the final iteration of the proposed method, whereas Figure [fig4](b) shows the \(y(x_1,x_2)\) values for the \(\theta\) values obtained in the final iteration of the Deep Ritz method.

Although the results shown in Figures [fig4](a) and [fig4](b) are qualitatively similar, significant differences can be observed when compared to the exact solution given by 53 . These differences are highlighted in Figures [fig5](a) and [fig5](b), which display the absolute errors relative to the exact solution 53 for the proposed method and the Deep Ritz method, respectively.

The absolute errors obtained by the proposed method are less than 0.0005, whereas those obtained by the Deep Ritz method are bounded above by 0.09. Note that in the proposed method the error is zero at the boundaries, whereas in the Deep Ritz method the maximum error is attained in \((x_1,x_2)=(0, 1)\). It is also worth highlighting that the results achieved by the proposed method were obtained using a neural network with a single hidden layer and \(N=5\), whereas the Deep Ritz method required two hidden layers with \(N=10\) to achieve comparable results. These findings further underscore the value of employing an ansatz that precisely satisfies the boundary conditions, as it eliminates the need for additional complexity while ensuring accuracy.

4.3 Example 3↩︎

In this section, we address the task of finding the minimum eigenvalue for the quantum oscillator problem under Dirichlet boundary conditions. For this type of problem, the eigenvalue can be approximated as follows: \[\begin{align} \label{eq38} &\lambda[y(x)]=\frac{\int_{\Omega}{\left[||\nabla y(x)||^2+||x||^2y(x)^2\right]dx}}{\int_{\Omega}{y(x)^2dx}} \end{align}\tag{56}\] where \(y(x)^2\) is interpreted as the probability density of locating a particle at a given position. This interpretation leads to the normalization condition \(\int_{\Omega}{y(x)^2dx}=1\). Note that expression 56 is the Rayleigh quotient [26] associated with the energy operator \(\hat{E}=-\nabla^2 y(x)+||x||^2y(x)\). It is important to point out that, in quantum mechanics, the square norm is associated with the product of a function and its complex conjugate, since the wave function (in our case, \(y(x)\)) is generally a complex function. Therefore, \(y(x)\) is determined by solving the following problem: \[\begin{align} \label{eq39} &\min_{y(x_1,x_2)}\left\{\begin{matrix} S=\frac{\int_{\Omega}{\left[||\nabla y(x)||^2+||x||^2y(x)^2\right]dx}}{\int_{\Omega}{y(x)^2dx}}\\ \\ y(x)=0, \forall x\in\partial\Omega \end{matrix} \right\} \end{align}\tag{57}\] In the case of Deep Ritz method, the minimization problem is given by: \[\begin{align} \label{eq40} &S_{\mathrm{DR}}=\frac{\int_{\Omega}{\left[||\nabla y(x)||^2+||x||^2y(x)^2\right]dx}}{\int_{\Omega}{y(x)^2dx}}+\beta\int_{\partial\Omega}{y(x)^2dx} \end{align}\tag{58}\] Where \(y(x)=N_{\mathrm{net}}(x,\theta)\) in 58 . In [18], an example of finding the minimum eigenvalue for a quantum oscillator was presented for \(\Omega=[-3,3]^d\) by adding a penalization term \(\left(\int_{\Omega}{y(x)^2dx}-1\right)^2\). However, after implementing this minimization with the added penalization, we observed that although the penalization improves the stabilization of the optimization process, the convergence remains roughly similar to the case of minimizing \(S_{\mathrm{DR}}\) without the penalization term \(\left(\int_{\Omega}{y(x)^2dx}-1\right)^2\):

So, we compare the results obtained from the minimization of \(S_{\mathrm{DR}}\) with those from the minimization of \(S\), as defined by 55 , and the following proposed ansatz: \[\begin{align} \label{eq41} &y(x)=p(x)N_{\mathrm{net}}(x,\theta) \end{align}\tag{59}\] where: \[\begin{align} \label{eq42} &p(x)=\prod_{j=1}^d\left(1-\left(\frac{x_i}{3}\right)^2\right) \end{align}\tag{60}\] We consider the case where \(d=3\), for which the exact minimum eigenvalue is \(\lambda_{\mathrm{exact}}=3\).

Figure 6 illustrates the values of \(\lambda\) via 56 with respect to the iteration number of the gradient descent method. Two cases have been considered: one where \(S\) is minimized using the ansatz provided in 59 , and another where \(S_{\mathrm{DR}}\) is minimized. In both cases, \(N=50\) is used, with \(\beta=15\). The initial values of the \(\theta\) parameters were randomly sampled from a Gaussian distribution with a mean of 0 and a standard deviation of 0.1, and these initial values were identical for both considered cases.

We observe that convergence is achieved more rapidly with the proposed method than with the Deep Ritz method. Specifically, a value of \(\lambda=3.0096\) (i.e., a relative error of 0.3% with respect to the exact value) is obtained using the proposed method, whereas the Deep Ritz method reaches \(\lambda=3.5939\) (i.e., a relative error of 19.8% with respect to the exact value) in the final iteration. Moreover, spurious oscillations appear in the values of \(\lambda\) for the Deep Ritz method. To further investigate this, we performed numerical simulations using the Deep Ritz method with an added penalization term \(\left(\int_{\Omega}{y(x)^2dx}-1\right)^2\). While this adjustment eliminates the oscillations, the overall convergence behavior remains similar.

5 Conclusions↩︎

An ansatz satisfying the boundary conditions for the variational problem has been proposed. This ansatz eliminates the need for a penalty term to account for the boundary conditions, as required in the Deep Ritz method. The ansatz comprises a term representing boundary behavior and another term formed by the product of a neural network and a polynomial. The polynomial is designed to vanish at the boundaries of the domain, ensuring that the boundary conditions are satisfied. Within the rest of the domain, the behavior of the ansatz is primarily governed by the product of the neural network and the polynomial.

A central contribution of this work is that all supporting theorems and corollaries have been established in Sobolev spaces. This is particularly relevant because the natural notion of approximation in variational problems involves not only the function itself but also its derivatives. Accordingly, proving density and approximation properties in Sobolev spaces provides a mathematically rigorous justification for the admissibility and effectiveness of the proposed ansatz within the Ritz framework. In particular, we have shown that the product of a polynomial that vanishes linearly near the domain boundary and a single-hidden-layer neural network is dense in the corresponding Sobolev space. This result guarantees that the proposed ansatz is expressive enough to approximate the solution of the variational problem with respect to the natural Sobolev norm.

Several drawbacks have been addressed by the proposed ansatz in various examples, demonstrating its effectiveness in overcoming key challenges:

• The dependence on hyperparameters is unnecessary for enforcing the boundary conditions.

• The proposed ansatz effectively reduces complexity while maintaining accuracy, highlighting its practical advantage in addressing variational problems.

• While the Deep Ritz method achieves the lowest action values, the proposed method offers the solution closest to the exact one, emphasizing the critical role of accurately satisfying boundary conditions to prevent misleading outcomes in the optimization of the action.

Declarations↩︎

• Funding: This work was supported by the ‘Plan Estatal de Investigación’, under the project PID2022-138144NB-I00 entitled ‘Simetría y coherencia en sistemas críticos cuánticos’.

• Conflict of interest/Competing interests: Not applicable

• Ethics approval and consent to participate: Not applicable

• Consent for publication: Not applicable

• Data availability: Not applicable

• Materials availability: Not applicable

• Code availability: Available under a reasonable request to the authors.

• Author contribution: Not applicable

References↩︎