Moving interface, PINNs, Extended variable technique, Error analysis

1 Introduction↩︎

Interface problems arise in many scientific and engineering models that involve multiple materials with different chemical or physical properties (see [1], [2] and references therein). In these models, the computational domain is composed of subdomains that are separated by smooth curves (or surfaces) known as interfaces. The interface geometry itself may be dynamic; that is, the interface between different subdomains varies in time. The low global regularity of the solutions and the evolving interface pose challenges for numerical simulations, particularly when the interface undergoes large deformations.

Let \(\Omega \subset \mathbb{R}^d\), \(d\in \mathbb{N}\), be a fixed bounded domain and let \(\Omega^+(t)\) and \(\Omega^-(t)\) be two subdomains of \(\Omega\) separated by the evolving interface \(\Gamma(t)\) on a time interval \([0,T_{\text{end}}]\). In this paper, we consider the following parabolic equation: \[\label{eq1461} \left\{ \begin{align} \partial_t u - \nabla \cdot \left( \beta \nabla u\right) &= f \quad &&\mathrm{in}\; \Omega = \Omega^+(t) \cup \Omega^-(t), \quad t \in [0, T_{\text{end}}], \\ u\left( \cdot, t\right) &= g \quad &&\mathrm{on}\; \partial \Omega, \quad t \in [0, T_{\text{end}}], \\ u\left( \cdot, 0\right) &= u_0 \quad &&\mathrm{in}\; \Omega = \Omega^+(0) \cup \Omega^-(0), \end{align}\right.\tag{1}\] with the jump conditions on the moving interface \(\Gamma(t)\): \[\begin{align} [u]_{\Gamma(t)} &:= u^+|_{\Gamma(t)} - u^-|_{\Gamma(t)} = h_D, \quad t\in [0, T_{\text{end}}] \tag{2},\\ [\beta \nabla u \cdot \mathbf{n}]_{\Gamma(t)}&:= \beta^+ \nabla u^+ \cdot \mathbf{n}|_{\Gamma(t)} - \beta^- \nabla u^- \cdot \mathbf{n}|_{\Gamma(t)} = h_N, \quad t\in [0, T_{\text{end}}] \tag{3}, \end{align}\] where \(f \in L^2(0,T_{\text{end}}; L^2(\Omega))\), and \(\mathbf{n}\) is the unit normal vector to \(\Gamma(t)\) pointing from \(\Omega^-(t)\) to \(\Omega^+(t)\). The restrictions of \(u\) on \(\Omega^+(t)\) and \(\Omega^-(t)\) are denoted by \(u^+ = u|_{\Omega^+(t)}\) and \(u^- = u|_{\Omega^-(t)}\). The coefficient function \(\beta\) is a piecewise positive constant defined as follows: \[\beta(\boldsymbol{x}, t) = \begin{cases} \beta^+, & \text{for } \boldsymbol{x} \in \Omega^+(t), \\ \beta^-, & \text{for } \boldsymbol{x} \in \Omega^-(t). \end{cases} \label{eq1466}\tag{4}\]

Furthermore, we assume that there is a prescribed velocity field \(\mathcal{V}(\boldsymbol{x}, t)\) that governs the movement of the interface, i.e., \[\dfrac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t} = \mathcal{V}(\boldsymbol{x}, t), \quad \boldsymbol{x} \in \Gamma(t). \label{eq1467}\tag{5}\]

The parabolic moving interface problem (1 )–(5 ) appears in many applications, such as the Stefan problem [3] and the Burton–Cabrera–Frank-type model [4]. For instance, \(u\) represents the temperature and \(\mathcal{V}\) is determined by the heat flux across the interface in Stefan problems. Under appropriate assumptions on the initial, boundary, and jump conditions (see Section 3 for detailed discussions), it can be proved that Eqs. (1 )–(5 ) admit a unique solution \(u\) belonging to the space \(L^2(0, T_{\text{end}}; H^1(\Omega))\).

It is well known that moving interface problems pose significant challenges for numerical simulations because the computational domain itself evolves. When employing traditional finite element methods (FEMs), one must ensure the generation of body-fitted meshes that conform to the interface; otherwise, the accuracy of the numerical solution will be severely degraded (cf. [5]). Chen and Zou [3] studied the linear FEM on nearly fitted quasi-uniform meshes and proved optimal-order error estimates up to some logarithmic factors. In general, mesh regeneration should be minimized or even avoided whenever possible, as mesh generation is notoriously laborious and time-consuming, particularly for complex interface geometries and higher-dimensional PDEs. Consequently, arbitrary Lagrangian–Eulerian methods [6] and space-time finite element methods [7] were proposed to overcome this difficulty.

To completely eliminate the need for mesh generation in solving PDEs, unfitted methods have attracted significant attention over recent decades. Unfitted methods allow the interface to cut through mesh elements, while special techniques are required to incorporate the jump conditions across the interface. To handle interface-cutting elements, one class of methods enforces the jump conditions through computational schemes such as CutFEM [8], [9] and IPFEM [10]. Another class of methods constructs specialized basis functions within interface elements that are designed to satisfy the interface jump conditions, such as GFEM [11], [12], XFEM [13], and IFEM [14], [15]. Guo [16] presented the first discrete analysis and optimal error estimates for the backward Euler IFEM applied to parabolic moving interface problems. All the aforementioned methods are effective for solving interface problems and eliminate the need for body-fitted mesh regeneration at each time step. However, they still require either careful handling of small cut elements [17] near the interface or the construction of specialized basis functions [15] within interface elements that satisfy the jump conditions.

Solving interface problems via deep learning methods has attracted extensive attention in recent years. As mesh-free methods, neural networks are particularly well suited for addressing problems in computational domains with complex geometries and higher-dimensional settings. Deep learning methods effectively mitigate the curse of dimensionality and eliminate the need for handling solution meshes. The well-known PINNs [18] and DeepONet [19], along with various improved methods [20]–[25], have been successfully used to solve elliptic interface problems, yielding promising results. To address more challenging moving interface problems, Wang et al. [26] proposed a multi-network framework for Stefan problems. In their approach, two distinct networks approximate the temperature fields in the two phases, and an additional network captures the interface evolution, serving essentially to classify the domain points. Lin et al. [27] extended the discontinuous ELM method by advancing the solution in time with an implicit–explicit discretization scheme. PF-PINNs [28] were developed for two-dimensional immiscible incompressible two-phase flows.

In this work, we develop an improved PINN method, termed the Extended Interface PINN (XI-PINN), for solving moving interface problems. The method is built upon the extended variable technique (EVT) [20], [21]. Compared with existing numerical methods for moving interface problems, this approach offers the following main advantages:

• As a meshfree method, it avoids the difficulties that traditional grid-based methods face when dealing with high-dimensional problems or large deformations of the interface geometry.

• Unlike existing multi-network approaches [26], [29], XI-PINN treats both subdomains with a single neural network. This means the network can adaptively distribute its parameters according to the complexity of the solution in each region, leading to a more uniform error distribution across subdomains. In contrast, multi-network methods that assign independent networks to different subdomains may suffer from unsuitable a priori parameter allocation, which can degrade the accuracy of the coupled solution for interface problems.

• By designing the extension variable, XI-PINN can accurately capture various types of interface conditions, circumventing the limitations inherent in diffuse interface models.

• When the solution exhibits no jump across the interface, XI-PINN can exactly enforce the zero-jump condition [30]. This is particularly important for two-phase flow problems, as it helps preserve the correct physical behavior.

One limitation of the EVT is that it presupposes an explicitly given level set function to describe the interface. In practice, obtaining this function typically requires solving a level set evolution equation, which is a nontrivial task when an analytical expression is unavailable. To overcome this difficulty, we represent the level set function using an additional neural network. By employing a tailored temporal training strategy, this network can accurately track the moving interface, even under large deformations.

We further present a comprehensive error analysis for the XI-PINN framework. Following recent developments [31], [32], we decompose the total error into three distinct components: the approximation error associated with the neural network space, the statistical error, and the optimization error. The bounds for the approximation error are guaranteed by the universal approximation properties of deep neural networks [33], while the optimization error vanishes in the over-parameterized regime, as established in the literature [34]. Consequently, our primary theoretical contribution lies in providing a detailed proof of the statistical error specifically tailored to parabolic moving interface problems, thereby completing the overarching analytical framework.

We evaluate the performance of XI-PINN through a series of numerical experiments, including high-dimensional problems, the Oseen equations, and test cases characterized by severe interface deformations. Leveraging second-order optimization algorithms, XI-PINN delivers highly accurate solutions. In cases where the exact solution is continuous across the interface, or when employing a small network architecture, XI-PINN exhibits better accuracy compared to multi-network approaches. Furthermore, for a benchmark problem with the interface condition \(h_D=0\), we show the spectrum of the neural tangent kernel (NTK). The results indicate that XI-PINN yields a more favorable spectral distribution than the Vanilla-PINN, demonstrating that within a single-network framework, XI-PINN effectively mitigates the spectral bias of the Vanilla-PINN and successfully captures sharp interface features.

The rest of this paper is organized as follows. In Section 2, we first present the fundamental framework of Deep Neural Networks (DNNs), then introduce the XI-PINN method and propose our novel approach for solving the level set function. In Section 3, we conduct a detailed error analysis of the XI-PINN method. Numerical results are presented in Section 4 to demonstrate the accuracy and efficacy of XI-PINN. Finally, we give some concluding remarks in Section 5.

2 Extended Interface PINN Method↩︎

In this section, we first review preliminaries on DNNs and the EVT. Subsequently, we introduce the XI-PINN method for solving parabolic interface problems. This method efficiently handles problems with geometrically complex interfaces and can be seamlessly extended to higher-dimensional scenarios.

2.1 DNNs↩︎

We consider a standard fully connected feedforward neural network consisting of multiple linear transformations and nonlinear activation functions. Let the network have depth \(L \ge 2\), \(L \in \mathbb{Z}\), and let \(\{ n_l \}_{l=0}^L\) denote the layer dimensions, where \(n_0\) is the input dimension and \(n_L\) is the output dimension. For each layer \(l = 1, \dots, L\), define an affine transformation \(\mathbf{T}^l: \mathbb{R}^{n_{l-1}} \to \mathbb{R}^{n_l}\) as follows: \[\mathbf{T}^l(\mathbf{X}^{l-1}):= \mathbf{A}^l \mathbf{X}^{l-1} + \mathbf{b}^l,\] where \(\mathbf{A}^l \in \mathbb{R}^{n_l \times n_{l-1}}\) is the weight matrix and \(\mathbf{b}^l \in \mathbb{R}^{n_l}\) is the bias vector. The \(l\)-th hidden layer (\(l = 1, \dots, L-1\)) applies an element-wise nonlinear activation function \(\sigma: \mathbb{R} \to \mathbb{R}\): \[\mathbf{X}^{l} = f^{l}(\mathbf{X}^{l-1}) := \sigma\bigl( \mathbf{T}^l( \mathbf{X}^{l-1} ) \bigr).\] The final layer (\(l = L\)) produces the output via a purely affine transformation without activation, i.e., \(\mathbf{X}^L = \mathbf{T}^L ( \mathbf{X}^{L-1})\). The complete \(L\)-layer DNN is then the composition \[\mathcal{NN}(\boldsymbol{x}; \theta) := \bigl( \mathbf{T}^L \circ f^{L-1} \circ \cdots \circ f^1 \bigr) (\mathbf{X}^0), \quad \mathbf{X}^0 = \boldsymbol{x},\] where \(\boldsymbol{x} \in \mathbb{R}^{n_0}\) is the input of the neural network, and the set of trainable parameters \(\{ \mathbf{A}^l, \mathbf{b}^l\}_{l=1}^L\) is collected into a parameter vector \(\theta\). The integers \(L\) and \(W:=\max \{ n_l,\;l=1,\dots,L\}\) are called the depth and width of the DNN, respectively.

Following the definition in [35], we denote by \(\mathcal{H}_\sigma(L,N_\theta,B_\theta)\) the collection of DNN functions of depth \(L\), with at most \(N_\theta\) nonzero trainable parameters, where each parameter is bounded in absolute value by \(B_\theta\), i.e., \[\begin{align} \mathcal{H}_\sigma( L,N_\theta,B_\theta) := \bigl\{ \text{L-layer } \mathcal{NN}(\boldsymbol{x}; \theta): |\theta|_{\ell^0} \le N_\theta,\;|\theta|_{\ell^\infty} \le B_\theta \bigr\}. \end{align}\] Here \(|\cdot|_{\ell^0}\) and \(|\cdot|_{\ell^\infty}\) denote the number of nonzero entries and the maximum norm of a vector, respectively. For brevity, we often write \(\mathcal{H}\) to mean \(\mathcal{H}_\sigma(L,N_\theta,B_\theta)\).

2.2 EVT↩︎

The EVT has demonstrated remarkable success in solving elliptic interface problems, enabling highly accurate approximations of solutions with low regularity across interfaces using a single neural network (see [20], [21], [23], [25]).

We first introduce two important functions: the level set function \(\phi: \Omega \times [0,T_{\text{end}}] \to \mathbb{R}\) and the indicator function \(\chi: \Omega \times [0,T_{\text{end}}] \to \mathbb{R}\). The level set function \(\phi\) characterizing the interface satisfies [36] \[\phi(\boldsymbol{x}, t) \begin{cases} <0, & \text{if } \boldsymbol{x} \in \Omega^-(t), \\ =0, & \text{if } \boldsymbol{x} \in \Gamma(t), \\ >0, & \text{if } \boldsymbol{x} \in \Omega^+(t), \end{cases}\] and the indicator function is defined as \[\chi(\boldsymbol{x}, t)= \begin{cases} -1, & \text{if } \boldsymbol{x} \in \Omega^-(t), \\ 1, & \text{if } \boldsymbol{x} \in \Omega^+(t). \end{cases}\] It is clear that the indicator function can be written in terms of the level set function. Consequently, we define an extension function \(\widetilde{u}(\boldsymbol{x}, t, z)\) on \(\widetilde{\Omega}:= \Omega \times [0, T_{\text{end}}] \times \mathbb{R} \subset \mathbb{R}^{d+2}\) such that \[\label{eq2465} \widetilde{u}(\boldsymbol{x}, t, z(\boldsymbol{x}, t)) = u(\boldsymbol{x}, t), \quad \boldsymbol{x} \in \Omega, \quad t\in[0,T_{\text{end}}],\tag{6}\] where \(z\) represents the extension variable. Depending on the interface jump condition \(h_D\), the extension variable is defined as \[z( \boldsymbol{x}, t) = \begin{cases} \chi( \boldsymbol{x}, t), & \text{if } h_D \neq 0,\\[4pt] \Phi( \boldsymbol{x}, t):= | \phi( \boldsymbol{x}, t) |, & \text{if } h_D = 0, \end{cases}\] As demonstrated in [30], setting \(z=\Phi(\boldsymbol{x}, t)\) for the case \(h_D=0\) improves approximation accuracy, since the neural network component \(u_{\mathcal{NN}} \in \mathcal{H}\) naturally enforces continuity across the entire domain \(\Omega\). It can be observed that, by the definitions of \(\widetilde{u}\) and the extended variable \(z(\boldsymbol{x}, t)\), \(\widetilde{u}\) naturally exhibits discontinuities or jumps in its derivatives at the interface \(\Gamma(t)\).

2.3 XI-PINN↩︎

As discussed in Section 2.2, although \(u\) (or its derivative) is discontinuous across the interface, the EVT lifts it to a continuous function \(\widetilde{u}\) defined on \(\widetilde{\Omega} \subset \mathbb{R}^{d+2}\). We therefore approximate the solution of Eqs. 1 –5 by a DNN \(\widetilde{u}_\theta\in\mathcal{H}\) with input dimension \(d+2\), trained through the following discrete optimization problem: \[\label{eq2467} \operatorname*{arg\,min}_{\widetilde{u}_{\theta} \in \mathcal{H}} \ell = \sum_{i=1}^{5} \ell_i(\widetilde{u}_{\theta}),\tag{7}\] with \[\begin{align} \ell_1(\widetilde{u}_{\theta}) &= \dfrac{1}{N_\Omega} \sum_{\boldsymbol{s}_j \in \mathcal{S}_{\Omega}} \bigl| \partial_t \widetilde{u}_{\theta}(\boldsymbol{s}_j, z(\boldsymbol{s}_j)) -\nabla \cdot \bigl( \beta \nabla \widetilde{u}_{\theta}(\boldsymbol{s}_j, z(\boldsymbol{s}_j)) \bigr) - f(\boldsymbol{s}_j) \bigr|^2, \tag{8} \\ \ell_2(\widetilde{u}_{\theta}) &= \dfrac{1}{N_{\partial \Omega}} \sum_{\boldsymbol{s}_j \in \mathcal{S}_{\partial \Omega}} \bigl| \widetilde{u}_{\theta}(\boldsymbol{s}_j, z(\boldsymbol{s}_j)) - g(\boldsymbol{s}_j) \bigr|^2, \tag{9} \\ \ell_3(\widetilde{u}_{\theta}) &= \dfrac{1}{N_{\Omega_0}} \sum_{\boldsymbol{s}_j \in \mathcal{S}_{\Omega_0}} \bigl| \widetilde{u}_{\theta}(\boldsymbol{s}_j, z(\boldsymbol{s}_j)) - u_0(\boldsymbol{s}_j) \bigr|^2, \tag{10}\\ \ell_4(\widetilde{u}_{\theta}) &= \dfrac{1}{N_{\Gamma}} \sum_{\boldsymbol{s}_j \in \mathcal{S}_{\Gamma}} \bigl| \bigl[ \beta \nabla \widetilde{u}_{\theta}(\boldsymbol{s}_j, z(\boldsymbol{s}_j)) \cdot \mathbf{n} \bigr]_{\Gamma(t)} - h_N(\boldsymbol{s}_j) \bigr|^2, \tag{11}\\ \ell_5(\widetilde{u}_{\theta}) &= \dfrac{1}{N_{\Gamma}} \sum_{\boldsymbol{s}_j \in \mathcal{S}_{\Gamma}} \bigl| \bigl[ \widetilde{u}_{\theta}(\boldsymbol{s}_j, z(\boldsymbol{s}_j)) \bigr]_{\Gamma(t)} - h_D(\boldsymbol{s}_j) \bigr|^2, \tag{12} \end{align}\] where \(\mathcal{S}_\Omega\), \(\mathcal{S}_{\partial\Omega}\), \(\mathcal{S}_{\Omega_0}\), and \(\mathcal{S}_\Gamma\) denote the sets of uniformly sampled collocation points within the respective domains, with \(N_\Omega, N_{\partial \Omega}, N_{\Omega_0}\), and \(N_{\Gamma}\) denoting the total number of points in each corresponding set.

2.4 Construction of the Level Set Function↩︎

The level set function plays a pivotal role in the EVT, as it determines the location of points relative to the subdomains, thereby defining the indicator function. Furthermore, when \(h_D=0\), it is used to directly construct the extended variable. Consequently, obtaining the time-evolving level set function is essential for the proposed framework. Assuming the velocity field \(\mathcal{V}\) is continuously extended to the entire domain \(\Omega\) such that the trajectory of every point \(\boldsymbol{x} \in \Omega\) is governed by the ODE (5 ), the level set function \(\phi(\boldsymbol{x}, t)\) satisfies the following advection equation: \[\begin{align} \label{advection32equation} \frac{\partial \phi(\boldsymbol{x}, t)}{\partial t} + \mathcal{V} \cdot \nabla \phi(\boldsymbol{x}, t) &= 0, \quad \boldsymbol{x} \in \Omega,\; t \in [0, T_{\text{end}}],\\ \phi(\boldsymbol{x}, 0) &= \phi_0(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega, \end{align}\tag{13}\] where \(\phi_0(\boldsymbol{x})\) denotes the prescribed initial level set function at \(t=0\). The governing equation may be solved using the method of characteristics. For any point \(\boldsymbol{x}_0 \in \Omega\) at \(t=0\), provided that its position at \(t=T\) admits an explicit representation \(\boldsymbol{x}_T = \boldsymbol{x}_0 + \mathcal{F}(\boldsymbol{x}_0, T)\), where \(\mathcal{F}\) is the displacement function (i.e., \(\mathcal{F}(\boldsymbol{x}_0, T) = \boldsymbol{x}_T - \boldsymbol{x}_0\)), the level set function can be constructed as \[\phi(\boldsymbol{x}, T) = \phi_0\bigl( \boldsymbol{x} - \mathcal{F}(\boldsymbol{x}, T) \bigr).\] Alternatively, to bypass solving the evolution equation entirely, we can employ a neural network to directly construct an appropriate level set function. For a given reference time \(t_0 \in [0,T_{\text{end}}]\), we define the flow map as follows: \[\label{eq24616} \mathbf{X}: \Gamma(t_0) \times [t_0, T_{\text{end}}] \to \bigcup_{t\in[t_0,T_{\text{end}}]}\Gamma(t), \quad \boldsymbol{x}_t = \mathbf{X}(\boldsymbol{x}, t; t_0), \quad t\in[t_0,T_{\text{end}}],\tag{14}\] where \(\boldsymbol{x} \in \Gamma(t_0)\) and \(\boldsymbol{x}_t \in \Gamma(t)\). The flow map can be obtained using the solution of the initial-value problem (5 ): \[\frac{\mathrm{d}}{\mathrm{d}t} \mathbf{X}(\boldsymbol{x},t; t_0) = \mathcal{V}\bigl( \mathbf{X}(\boldsymbol{x},t; t_0), t \bigr), \quad \mathbf{X}(\boldsymbol{x},t_0; t_0) = \boldsymbol{x}, \quad \boldsymbol{x} \in \Gamma(t_0).\] The inverse of \(\mathbf{X}(\boldsymbol{x},t; t_0)\) is denoted by \(\widetilde{\mathbf{X}}(\boldsymbol{x}, t; t_0): \Gamma(t) \times [t_0 ,T_{\text{end}}] \to \Gamma(t_0)\). We extend the mapping \(\widetilde{\mathbf{X}}\) to the whole domain and define the extension mapping \[\widehat{\mathbf{X}} : \Omega \times [t_0,T_{\text{end}}] \to \mathbb{R}^d,\] which is required to satisfy the following properties: \[\label{flow32map32properties} \small \widehat{\mathbf{X}}(\mathbf{x},t;t_0)|_{\mathbf{x} \in \Gamma(t)} = \widetilde{\mathbf{X}}(\mathbf{x},t;t_0) \text{ and } \widehat{\mathbf{X}}(\mathbf{x},t;t_0) \text{ is orientation-preserving and injective.}\tag{15}\] Therefore, the level set function at time \(t\) can be defined as \[\label{eq24620} \phi( \boldsymbol{x}, t) = \phi_0\bigl( \widehat{\mathbf{X}}(\boldsymbol{x},t;0) \bigr).\tag{16}\] The gradient of this function can be computed via the chain rule: \[\nabla_{\boldsymbol{x}} \phi(\boldsymbol{x},t) = \nabla_{\widehat{\mathbf{X}}} \phi_0\bigl( \widehat{\mathbf{X}} \bigr) \, \nabla_{\boldsymbol{x}} \widehat{\mathbf{X}}(\boldsymbol{x},t;0),\] where \(\nabla_{\widehat{\mathbf{X}}} \phi_0\) is the gradient vector and \(\nabla_{\boldsymbol{x}} \widehat{\mathbf{X}}\) is the Jacobian matrix.

Next, we employ neural networks to approximate the mapping \(\widehat{\mathbf{X}}\). Denote by \(\mathcal{F}_\theta\) the neural network function, and define the approximating function for \(\widehat{\mathbf{X}}\) as \[\label{eq24621} \widehat{\mathbf{X}}_\theta( \boldsymbol{x}, t; 0) = \mathcal{F}_\theta( \boldsymbol{x}, t) + \boldsymbol{x}.\tag{17}\] Let \(\mathbf{X}_{\mathrm{RK}}\) be the approximation of \(\mathbf{X}\) computed using the 4th-order explicit Runge–Kutta method [37]. We then determine the approximate displacement function \(\mathcal{F}_\theta\) and the resulting flow map \(\widehat{\mathbf{X}}_\theta\) by minimizing the following loss function: \[\ell_{\mathrm{LSF}} = \frac{1}{N_{\Omega}} \sum_{i=1}^{N_{\Omega}} \bigl| \mathcal{F}_\theta( \boldsymbol{x}_{i, \Omega}, 0) \bigr|^2 + \frac{1}{N_t N_\Gamma} \sum_{j=1}^{N_t} \sum_{i=1}^{N_\Gamma} \bigl| \widehat{\mathbf{X}}_\theta\bigl( \mathbf{X}_{\mathrm{RK}}(\boldsymbol{x}_{i,\Gamma}, t_j; 0), t_j; 0 \bigr) - \boldsymbol{x}_{i,\Gamma} \bigr|^2,\] where \(\{ \boldsymbol{x}_{i,\Omega} \}_{i=1}^{N_{\Omega}}\) and \(\{ \boldsymbol{x}_{i,\Gamma} \}_{i=1}^{N_{\Gamma}}\) are training points sampled from the domain \(\Omega\) and the initial interface \(\Gamma(0)\), respectively, and \(\{ t_j \}_{j=1}^{N_{t}}\) forms a sufficiently fine temporal discretization of \([0,T_{\text{end}}]\).

Notably, we require \(\widehat{\mathbf{X}}_\theta\) to satisfy the properties 15 . If the approximate flow map \(\widehat{\mathbf{X}}_\theta\) fails to be injective or orientation-preserving, the level set function \[\phi_\theta(\boldsymbol{x}, t) := \phi_0\bigl( \widehat{\mathbf{X}}_\theta(\boldsymbol{x}, t; 0) \bigr)\] loses its validity as an implicit interface representation. Lack of injectivity produces global folding: two distinct points \(\boldsymbol{x}_1 \neq \boldsymbol{x}_2\) in the current domain can be pulled back to the same reference point \(\boldsymbol{y}\), forcing \(\phi_\theta(\boldsymbol{x}_1, t) = \phi_\theta(\boldsymbol{x}_2, t) = \phi_0(\boldsymbol{y})\). If \(\boldsymbol{y}\) lies on the initial interface, the zero level set spuriously bifurcates, creating non-physical internal boundaries and corrupting the subdomain indicator function. Violation of orientation preservation, on the other hand, leads to local Jacobian degeneracy or sign reversal; the map can then erroneously assign points that belong to \(\Omega^+(t)\) to \(\Omega^-(0)\), destroying the inside/outside sign convention. These defects are the direct analogues of mesh entanglement and element inversion in finite element computations.

To mitigate this issue, we adopt a temporal partitioning strategy, dividing \([0, T_{\text{end}}]\) into smaller sub-intervals: \[[0, T_{\text{end}}] = \bigcup_{k=1}^{K}[T_{k-1}, T_k], \quad T_0 = 0,\; T_K = T_{\text{end}},\] where the set \(\{ T_k \}_{k=1}^{K} \subseteq \{ t_j \}_{j=1}^{N_t}\) can be adaptively chosen by monitoring the Jacobian determinant of \(\widehat{\mathbf{X}}_\theta\). For any \(t\) in the \(k\)-th sub-interval \([T_{k-1}, T_k]\), we can decompose \(\widehat{\mathbf{X}}(\boldsymbol{x}, t; 0)\) into a composition of sub-mappings: \[\label{eq24624} \widehat{\mathbf{X}}(\boldsymbol{x}, t;0) = \bigl[ \widehat{\mathbf{X}}(\cdot, t;T_{k-1}) \circ \widehat{\mathbf{X}}(\cdot, T_{k-1}; T_{k-2}) \circ \cdots \circ \widehat{\mathbf{X}}(\cdot, T_{1}; 0) \bigr](\boldsymbol{x}).\tag{18}\] For each short time sub-interval \([T_{k-1}, T_k]\), we compute a high-precision approximate sub-mapping \(\widehat{\mathbf{X}}_{\theta, k}\) using a small neural network. Note that employing an excessively large network may lead to overfitting, which in turn could introduce undesirable mapping entanglement.

Moreover, the displacement of interface points, \(\Delta \boldsymbol{x}\), remains moderate over short time intervals, effectively preventing mesh entanglement, particularly for problems involving large interface deformations. Correspondingly, we employ a series of sub-networks \(\widehat{\mathbf{X}}_{\theta, k}\) to approximate \(\widehat{\mathbf{X}}(\cdot, t; T_{k-1})\) for \(t \in [T_{k-1}, T_k]\), with the loss function adjusted as follows: \[\label{eq24625} \begin{align} \ell_{\mathrm{LSF}}^k = &\frac{1}{N_{\Omega}} \sum_{i=1}^{N_{\Omega}} \bigl| \mathcal{F}_\theta(\boldsymbol{x}_{i, \Omega}, T_{k-1}) \bigr|^2 + \frac{1}{N_{t,k} N_{\Gamma}} \sum_{j=1}^{N_{t,k}} \sum_{i=1}^{N_{\Gamma}} \\ & \qquad \bigl| \widehat{\mathbf{X}}_\theta\bigl( \mathbf{X}_{\mathrm{RK}}(\boldsymbol{x}_{i,\Gamma}, t_j^k; 0) , t_j^k; T_{k-1} \bigr) - \mathbf{X}_{\mathrm{RK}}(\boldsymbol{x}_{i,\Gamma}, T_{k-1}; 0) \bigr|^2, \end{align}\tag{19}\] where \(\{ t_j^k \}_{j=1}^{N_{t,k}}\) denotes a discrete time series within \([T_{k-1}, T_k]\), sampled from the global set \(\{ t_j \}_{j=1}^{N_t}\).

To detect potential mesh entanglement (indicated by negative Jacobian determinants), we locally linearize the mapping via a first-order Taylor expansion. This approach employs uniform sampling over \(\Omega\) with a sufficiently small grid spacing \(h\), which ensures the validity of the local linear approximation: \[\label{eq24626} \widehat{\mathbf{X}}_{\theta, k}(\boldsymbol{x},t;T_{k-1}) \approx \widehat{\mathbf{X}}_{\theta, k}(\mathbf{x}_0,t;T_{k-1}) + \mathbf{J}(\boldsymbol{x}_0)(\boldsymbol{x} - \boldsymbol{x}_0), \quad \forall \boldsymbol{x} \in B_\varepsilon(\boldsymbol{x}_0),\tag{20}\] where \(\mathbf{J}(\boldsymbol{x}_0):=\nabla_{\boldsymbol{x}} \widehat{\mathbf{X}}_{\theta,k}\) denotes the Jacobian matrix of the mapping and \(B_\varepsilon \subset \Omega\) is a small neighborhood of the sample point \(\boldsymbol{x}_0 \in \Omega\). A negative Jacobian determinant (\(\det(\mathbf{J}) < 0\)) indicates that the mapping fails to preserve orientation. It is worth noting that, for a smooth map between two compact simply connected domains, orientation preservation (positive Jacobian determinant everywhere) guarantees global injectivity. To enforce this condition numerically, we impose a positivity constraint by maintaining \(\det(\mathbf{J}) > \delta\), where \(\delta>0\) is a prescribed threshold. A Jacobian determinant falling below \(\delta\) (i.e., \(\det(\nabla_{\boldsymbol{x}} \widehat{\mathbf{X}}_{\theta,k}) < \delta\)) serves as a topological early warning signal. At this critical point, we initialize a new neural network \(\widehat{\mathbf{X}}_{\theta,k+1}\) to learn the deformation mapping for the subsequent interval \([T_k, T_{k+1}]\).

We therefore develop an adaptive time-stepping algorithm, which consists of (1) a sequence of sub-mappings \(\{ \widehat{\mathbf{X}}_{\theta, k} \}_{k=1}^{K}\) and (2) a corresponding partition of the time domain into sub-intervals \(\{ [T_{k-1}, T_k] \}_{k=1}^{K}\) based on the discrete temporal grid \(\{ t_j \}_{j=1}^{N_t}\) of \([0,T_{\text{end}}]\). This adaptive network growth strategy ensures that \(\phi_\theta\) remains valid over the entire time domain. We describe the complete time-adaptive training procedure in Algorithm 1 as follows:

3 Error analysis↩︎

In this section, we present the error analysis of the XI-PINN framework introduced in Section 2. We assume that the domain \(\Omega \subseteq \mathbb{R}^d\) is bounded with a Lipschitz boundary \(\partial \Omega\), and \(\Gamma(t)\) is a closed smooth curve (or surface) within \(\Omega\). To simplify the analysis, we assume \(h_D = 0\). This implies that the neural network approximation \(\widetilde{u}_{\mathcal{NN}} \in \mathcal{H}\) remains continuous across the interface.

The space-time domain is denoted by \(Q:=\Omega \times [0,T_{\text{end}}]\) and \(\Sigma:= \partial \Omega \times [0,T_{\text{end}}]\). The extended space-time domain is defined as \[\widetilde{Q}:= Q \times \bigl[ \min_{(\boldsymbol{x}, t)\in Q} z(\boldsymbol{x}, t), \max_{(\boldsymbol{x}, t)\in Q} z(\boldsymbol{x}, t) \bigr].\] For each subdomain \(\omega \subseteq \Omega\), we denote by \(H^k(\omega)\) the standard Sobolev space equipped with the norm \(\|\cdot\|_{H^k(\omega)}\), and define the time-dependent Bochner space \(H^l(0, T; H^k(\omega))\) as the space of all measurable functions \(u: [0,T] \to H^k(\omega)\) with the norm \(\|\cdot\|_{H^l(0, T; H^k(\omega))}\). Here, we let \(H^k(\Omega^+(t) \cup \Omega^-(t))\) denote the space of functions that belong to \(H^k\) on each subdomain separately, equipped with the piecewise norm \[\|v\|_{H^k(\Omega^+(t) \cup \Omega^-(t))}:= \|v\|_{H^k(\Omega^+(t))} + \|v\|_{H^k(\Omega^-(t))}.\] We also employ the Besov space \(B_{2,1}^{1/2}(\omega):=\bigl( L^2(\omega), H^1(\omega) \bigr)_{1/2, 1}\) obtained by real interpolation (see [38]). The piecewise Besov space \(B_{2,1}^{1/2}(\Omega^+(t) \cup \Omega^-(t))\) and its norm \(\|v\|_{B_{2,1}^{1/2}(\Omega^+(t) \cup \Omega^-(t))}\) are defined analogously. Moreover, we have the embedding \(H^{1}(\omega) \hookrightarrow B_{2,1}^{1/2}(\omega)\). Let \(r\) and \(s\) be two non-negative real numbers, and define the space \(H^{r,s}(Q) := L^2(0,T; H^{r}(\Omega)) \cap H^s(0, T; L^2(\Omega))\) with the norm \(\|\cdot\|_{H^{r,s}}\); this is a Hilbert space (see [39]). For a function \(v(t)\), we write \(v' := \frac{\mathrm{d}}{\mathrm{d} t} v\) for its time derivative (or \(\partial_t v\) if \(v\) depends on space as well) and let \(v_0:=v(0)\) denote its initial value. For simplicity, we use the shorthand notation \(A \lesssim B\) for the inequality \(A \le C B\), where \(C>0\) is a generic constant independent of \(N_\Omega\), \(N_{\partial \Omega}\), \(N_{\Omega_0}\), and \(N_\Gamma\), the numbers of training points specified in the loss functions (8 )–(11 ).

3.1 Error analysis of the XI-PINN method↩︎

We first establish an energy estimate for the solution \(u\) to Eqs. (1 )-(5 ).

For subsequent error analysis, we define the continuous loss function as: \[\label{eq34612} {\mathcal{L}}(\widetilde{u}_{\theta}) = \mathop{\sum}_{i=1}^4 {\mathcal{L}_i}(\widetilde{u}_{\theta}),\tag{24}\] \[\begin{align} &{\mathcal{L}_1}(\widetilde{u}_{\theta}) = \int_{0}^{T_{\text{end}}}\int_{\Omega} \left| \partial_t \widetilde{u}_{\theta}({\boldsymbol{x}}, t, z({\boldsymbol{x}}, t)) - \nabla \cdot (\beta \nabla\widetilde{u}_{\theta}({\boldsymbol{x}}, t, z({\boldsymbol{x}}, t))) - f({\boldsymbol{x}}, t) \right|^2 \mathrm{d}{\boldsymbol{x}} \mathrm{d}t, \\ &\begin{aligned} {\mathcal{L}_2}(\widetilde{u}_{\theta}) = &\int_{0}^{T_{\text{end}}}\int_{\partial \Omega} { \mathop{\sum}_{|\alpha|=0}^2 }\left| \partial^\alpha_{\boldsymbol{x}} \left( \widetilde{u}_{\theta}({\boldsymbol{x}}, t, z({\boldsymbol{x}}, t)) - g({\boldsymbol{x}}, t)\right) \right|^2 \mathrm{d}{\boldsymbol{x}} \mathrm{d}t\\ &+\int_{0}^{T_{\text{end}}}\int_{\partial \Omega} \left| \partial_t \widetilde{u}_{\theta}({\boldsymbol{x}}, t, z({\boldsymbol{x}}, t)) - \partial_t g({\boldsymbol{x}}, t) \right|^2 \mathrm{d}{\boldsymbol{x}} \mathrm{d}t \end{aligned} \\ &{\mathcal{L}_3}(\widetilde{u}_{\theta}) = \int_{\Omega} {\mathop{\sum}_{|\alpha| = 0}^1 }\left| \partial^\alpha_{\boldsymbol{x}} \left( \widetilde{u}_{\theta}({\boldsymbol{x}}, 0, z({\boldsymbol{x}}, 0)) - u_0({\boldsymbol{x}}) \right) \right|^2 \mathrm{d}{\boldsymbol{x}}, \\ &{\mathcal{L}_4}(\widetilde{u}_{\theta}) = \int_{0}^{T_{\text{end}}}\int_{\Gamma(t)} \left| [\beta \nabla \widetilde{u}_{\theta}({\boldsymbol{x}}, t, z({\boldsymbol{x}}, t)) \cdot \mathbf{n}]_{\Gamma(t)} - h_N({\boldsymbol{x}}, t)\right|^2 \mathrm{d}{\boldsymbol{x}} \mathrm{d}t, \end{align}\] where \(\alpha = (\alpha_1, \alpha_2, \dots, \alpha_d) \in \left( \mathbb{N}\cup \left\lbrace 0\right\rbrace \right) ^d\) is a multi-index for spatial coordinates with \(|\alpha| = \sum_{k=1}^d \alpha_k\), and \(\partial_{{\boldsymbol{x}}}^\alpha = \dfrac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1} \cdots \partial x_d^{\alpha_d}}\). Further, the corresponding discrete form of loss function 24 is defined by \[\label{eq34617} {\mathcal{L}}^d(\widetilde{u}_{\theta}) = \mathop{\sum}_{i=1}^4 \widehat{\ell}_i(\widetilde{u}_{\theta}),\tag{25}\] with \(\widehat{\ell}_1(\widetilde{u}_{\theta})=\left|\Omega(t)\right| {\ell}_1(\widetilde{u}_{\theta})\), \(\widehat{\ell}_4(\widetilde{u}_{\theta}) =\left|\Gamma(t)\right| {\ell}_4(\widetilde{u}_{\theta})\), and \[\begin{align} &\begin{aligned} \widehat{\ell}_2(\widetilde{u}_{\theta}) &= \dfrac{\left| \partial \Omega(t)\right| }{N_{\partial \Omega}} \sum_{\boldsymbol{s}_j \in \mathcal{S}_{\partial \Omega}} \left[ { \mathop{\sum}_{|\alpha|=0}^2 }\left| \partial^\alpha_{\boldsymbol{x}} \left( \widetilde{u}_{\theta}(({\boldsymbol{s}}_j), z({\boldsymbol{s}}_j)) - g(\boldsymbol{s}_j)\right) \right|^2 \right. \\ &\left. \qquad \qquad \qquad \qquad + \left| \left( \partial_t \widetilde{u}_{\theta}((\boldsymbol{s}_j), z(\boldsymbol{s}_j)) - \partial_t g(\boldsymbol{s}_j)\right) \right|^2\right] , \end{aligned} \\ &\widehat{\ell}_3(\widetilde{u}_{\theta}) = \dfrac{\left| \Omega(0) \right|}{N_{\Omega_0}} \mathop{\sum}_{\boldsymbol{s}_j \in \mathcal{S}_{\Omega_0}} { \mathop{\sum}_{|\alpha| = 0}^1 } \left| \partial^\alpha_{\boldsymbol{x}}\left( \widetilde{u}_{\theta}\left(\boldsymbol{s}_j,z\left(\boldsymbol{s}_j\right) \right) - u_0\left(\boldsymbol{s}_j\right)\right) \right| ^2, \end{align}\]

Compared to (7 )-(11 ), the term \(\ell_5\) has been eliminated by exploiting the level-set function’s continuity, while derivative information has been incorporated into both the initial condition (10 ) and boundary term (9 ) inspired by the energy estimates. We remark that although the error analysis motivates including derivative terms in the loss function, our neural network implementation retains the original loss function defined in (7 )-(11 ).

Let \(\widetilde{u}^*\) and \(\widetilde{u}^*_d\) denote the best approximations in the function class \(\mathcal{H}\) according to formulas (24 ) and (25 ), respectively. Namely, \[\widetilde{u}^* := \mathop{\arg \min}_{v_\theta \in \mathcal{H}} {\mathcal{L}}(v_\theta), \quad \widetilde{u}_d^* := \mathop{\arg \min}_{v_\theta \in \mathcal{H}} {\mathcal{L}}^d(v_\theta).\] Further denote by \(\widetilde{u}^*_{d\mathscr{A}}\) the numerical solution obtained via optimization procedures based on the loss function (25 ). We then have the following result:

Following Theorem 1, we decompose the error into three distinct components:

• Approximation error \(\mathcal{E}_\text{approx}\): Best approximation error of the target function \(\widetilde{u} \in H^4(\widetilde{Q})\) in the neural network function class \(\mathcal{H}\) under the Sobolev norm on bounded domain \(\widetilde{Q}\);

• Statistical error \(\mathcal{E}_\text{stat}\): Error induced by finite sampling;

• Optimization error \(\mathcal{E}_\text{opt}\): Error from numerical optimization procedures.

Under over-parameterized assumptions, the optimization error can be assumed arbitrarily small, as established in the literature [34]. Therefore, we assume that the optimization error \(\mathcal{E}_{\text{opt}}\) can be arbitrarily small. For approximation error \(\mathcal{E}_\text{approx}\) in Sobolev space, the following approximation theorem holds [33]:

Theorem 2 demonstrates that the approximation error \(\mathcal{E}_{\text{approx}}\) of the neural network can be effectively controlled. Therefore, for \(\varepsilon>0\), there exists a DNN \(v_\theta \in \mathcal{H}\) such that \(\|v_\theta - \widetilde{u}\|_{H^4(\widetilde{Q})} \le \varepsilon\).

Next, we bound the statistical error via the Rademacher complexity \(\mathfrak{R}_n\) which measures the complexity of a collection of function by the correlation between function values with Rademacher random variables, i.e., with probability \(P\left( \omega =1\right) = P\left( \omega = -1\right) = \frac{1}{2}\) (cf. [41]).

To derive an upper bound estimate for the statistical error \(\mathcal{E}_\text{stat}\), we invoke the following several lemmas:

As established in Section 2.4, we have presented a method for approximating the level set function using neural networks, which is denoted by \(\phi_\theta({\boldsymbol{x}}, t)\). The zero level set of \(\phi_\theta({\boldsymbol{x}}, t)\) represents the approximate interface \[\Gamma_\theta(t) :=\left\{ {\boldsymbol{x}} | \phi_\theta({\boldsymbol{x}}, t) = 0, {\boldsymbol{x}} \in \Omega\right\}, \quad t\in[0, T_{\text{end}}].\] Although \(\Gamma_\theta(t)\) deviates from the true interface \(\Gamma(t)\), we establish that the solution remains controllable relative to the original problem’s solution under suitable assumptions with sufficiently small interface discrepancy. We quantify this interface deviation via the Hausdorff distance: \(\delta = \mathop{\max}_{t \in [0, T_{\text{end}}]} d_H\left( \Gamma(t), \Gamma_\theta(t)\right),\) where \(d_H\) denotes Hausdorff distance: \[d_H(A, B) := \mathop{\max}\left\lbrace \mathop{\sup}_{a\in A} \mathop{\inf}_{b\in B} \mathrm{dist}\left( a, b\right) , \mathop{\sup}_{b\in B} \mathop{\inf}_{a\in A} \mathrm{dist}\left( a, b\right) \right\rbrace.\]

We define a tubular neighborhood of width \(\delta > 0\) about the interface \(\Gamma(t)\) as follows: \[S_{\delta}\left( t\right) :=\left\lbrace {\boldsymbol{x}}| \mathrm{dist}({\boldsymbol{x}}, \Gamma(t)) < \delta \right\rbrace, \quad S_\delta:= \mathop{\bigcup}_{t \in [0, T_{\text{end}}]} S_\delta(t).\]

By use of the Lemma 2.1 in [45], we have:

Finally, we conclude the following result:

4 Numerical experiments↩︎

In this section, we present several numerical examples to demonstrate the effectiveness and accuracy of the XI-PINN method for moving interface problems, particularly emphasizing its ability to capture solution discontinuities across the interface. In all test cases, we use the hyperbolic tangent activation function \(\sigma(x) = \tanh(x)\). To minimize the loss function, we adopt the improved Levenberg–Marquardt (LM) method [46]. Compared to commonly used first-order optimizers (e.g., gradient descent), second-order methods have been shown to effectively mitigate the ill-conditioned loss landscapes of PINNs [47] and alleviate the gradient conflicts arising in multi-objective optimization [48]. Furthermore, while training PINNs via gradient descent often requires adaptively adjusting the penalty parameters of different residual terms to ensure convergence [49], [50]. Since we adopt a second-order optimizer throughout this work, we simply fix all penalty parameters to \(1\) for all neural network methods, exactly as in the loss function defined in 7 .

In all the following experiments, the optimization is terminated after \(2000\) epochs. We compare the proposed XI-PINN with the Vanilla-PINN [18] and with X-PINNs [29]. Vanilla-PINN, like XI-PINN, employs a single neural network to represent the solution in both subdomains, while X-PINNs use two separate networks to approximate the solution in \(\Omega^+(t)\) and \(\Omega^-(t)\), respectively. For Vanilla-PINN, the interface jump condition is approximated by a finite difference: \[\bigl[ \beta \nabla u \cdot \mathbf{n} \bigr]_{\Gamma(t)} \approx \beta^+ \nabla u(\boldsymbol{x} + \varepsilon \mathbf{n}, t) \cdot \mathbf{n} - \beta^- \nabla u(\boldsymbol{x}, t) \cdot \mathbf{n},\] where we take \(\varepsilon = 10^{-3}\). For X-PINNs, one only needs to enforce the interface conditions between the two networks (see, e.g., the treatment of elliptic interface problems in [22]). When a neural network is used to construct the level set function \(\phi_\theta\) within the XI-PINN framework, we denote the resulting scheme as \(\text{XI-PINN}^*\) and set the threshold \(\delta = 0.2\) in Algorithm 1. In the present work, we consider the case where the interface velocity field \(\mathcal{V}\) is known; therefore, the level set function \(\phi_\theta\) can be constructed independently prior to the training of XI-PINN. For a neural network approximation \(u_\theta\), we define the relative errors in the \(L^2(0,T;L^2(\Omega))\) and \(L^2(0,T;H^1(\Omega))\) norms as \[E_{L^2}^{\text{Rel}}(u_\theta) = \frac{\|u - u_\theta\|_{L^2(0,T;L^2(\Omega))}}{\|u\|_{L^2(0,T;L^2(\Omega))}}, \qquad E_{H^1}^{\text{Rel}}(u_\theta) = \frac{\|u - u_\theta\|_{L^2(0,T;H^1(\Omega))}}{\|u\|_{L^2(0,T;H^1(\Omega))}}.\] All experiments were run on an NVIDIA V100-SXM2-32GB GPU.

4.1 Parabolic equations↩︎

For this example we conduct a thorough comparison of Vanilla-PINN, X-PINNs, XI-PINN, and \(\text{XI-PINN}^*\) using networks of different sizes. In all the comparative experiments, the numbers of training points are set to \((N_\Omega, N_{\partial\Omega}, N_{\Omega_0}, N_\Gamma)=(10\text{K}, 3\text{K}, 0.3\text{K}, 2\text{K})\), and we keep the total numbers of parameters of networks of different architectures as close as possible. For \(\text{XI-PINN}^*\), the level set function \(\phi_\theta\) is constructed by a single fully connected network of architecture \((L,W)=(2,64)\), whose training takes \(59.2\) seconds. The decay of the loss \(\ell\) and the relative \(L^2(0,T;L^2(\Omega))\) error during training are shown in Fig. 2, and the final numerical results are reported in Table 1.

The results show that Vanilla-PINN fails to yield an accurate solution because it cannot capture discontinuities at the interface. Although the computational costs (measured by training time and including the overhead for training \(\widehat{\mathbf{X}}_\theta\) in \(\text{XI-PINN}^*\)) are similar across the three methods, XI-PINN achieves slightly higher accuracy than X-PINNs. As shown by the error distributions in Fig. 3 and Fig. 4, XI-PINN produces a smooth error field without interfacial jumps. This is because its input variables, \((\boldsymbol{x}, t, \Phi(\boldsymbol{x}, t))\), are constructed to be globally continuous. Furthermore, its single-network architecture shares parameters across subdomains, allowing for adaptive parameter allocation and leading to a more uniform error distribution. In contrast, X-PINNs use independent networks for each subdomain. This multi-network architecture only approximates the interface conditions and does not exactly satisfy continuity, which introduces localized error jumps. Moreover, due to the strongly coupled interface conditions, an error in one network easily propagates into the interior of the other subdomains. This error propagation becomes much more obvious when smaller networks are used (see Fig. 4).

For \(\text{XI-PINN}^*\), the zero level set of the learned function \(\phi_\theta\) has a slight geometric deviation from the exact interface. This geometric error causes a small drop in numerical accuracy, which is consistent with the theoretical bounds established in Theorem 4. Finally, we note a trade-off in the optimization process. While parameter sharing in XI-PINN improves accuracy, it results in a dense Jacobian matrix in the Levenberg-Marquardt (LM) algorithm. Conversely, the independent networks in X-PINNs naturally yield a block-diagonal Jacobian matrix, which provides a more favorable structure for computation.

For the XI-PINN, we investigate the effect of statistical error on the numerical accuracy by varying the number of training points. The numerical results in Table 2 show that, as the number of training points decreases, the approximation error grows even though the loss \(\ell\) remains very small. This indicates that reducing training points increases the statistical error, which eventually degrades the accuracy—a behavior consistent with the analysis in Section 3.

Both Vanilla-PINN and XI-PINN are able to enforce continuity across the interface, but Vanilla-PINN fails to capture the normal flux discontinuity and the overall training fails. To understand this, we examine the NTK matrices of the two methods. For this comparison we use a single-hidden-layer network with \(512\) neurons for both models, and set the numbers of training points to \((N_\Omega, N_{\partial\Omega}, N_{\Omega_0}, N_\Gamma)=(1000, 500, 200, 300)\). The eigenvalue spectra of the NTK matrices are plotted in Fig. 5. We observe that the eigenvalues of the XI-PINN NTK decay more slowly than those of Vanilla-PINN. Since small eigenvalues of the NTK are typically associated with high-frequency components of the solution[49], [51], XI-PINN is better able to capture the interface discontinuity while still maintaining continuity across the interface. We do not compare with X-PINNs here, because X-PINNs deal with a piecewise smooth problem, which is fundamentally different from the setting of Vanilla-PINN and XI-PINN.

4pt

For problems where the solution is discontinuous across the interface, i.e., \(h_D \neq 0\), we choose the indicator function as the extension variable. In this example we compare XI-PINN with X-PINNs under three different network sizes; the numerical results are shown in Fig. 6 and Table 3. The numbers of training points are set to \((N_\Omega, N_{\partial\Omega}, N_{\Omega_0}, N_\Gamma) = (20\text{K}, 6\text{K}, 0.6\text{K}, 4\text{K})\). For small network architectures, XI-PINN achieves higher accuracy by sharing parameters across the two subdomains. Unlike Example 1, for large networks XI-PINN does not gain an accuracy advantage over X-PINNs. This is because, when each subnetwork of X-PINNs has sufficiently many parameters, it can represent the solution in its own subdomain with high accuracy. In contrast, a pre-assigned small network may suffer from insufficient capacity in one subdomain, which in turn degrades the approximation in the other. Furthermore, when \(h_D \neq 0\), XI-PINN no longer benefits from exactly enforcing the zero-jump condition. In such scenarios, both XI-PINN and X-PINNs must learn the interface jump conditions, although XI-PINN uniquely encodes this discontinuity via the extended variable. Consequently, for high-accuracy computations utilizing large networks in the presence of non-homogeneous interface jumps, X-PINNs emerge as a more favorable alternative. The numerical solutions and absolute error distributions on the hyperplane \(z=0\) at \(t=0.5\) for both methods are given in Fig. 7. These results also highlight the advantage of neural networks for high-dimensional PDEs; handling such moving interface problems via grid-based space-time finite element methods is a nontrivial task.

4.2 Oseen equations↩︎

The XI-PINN method can also be applied to solve more complex time-varying interface problems such as Oseen equations [52]. The model is described as follows: \[\label{eq4461} \left\{\begin{align} \partial_t \boldsymbol{u}^{\pm} + (\mathcal{V}\cdot\nabla)\boldsymbol{u}^{\pm} - \nu^{\pm} \Delta \boldsymbol{u}^{\pm} + \nabla p^{\pm} &= \boldsymbol{f}^{\pm}, \qquad \operatorname{div}\boldsymbol{u}^{\pm} = 0, && \text{in } \Omega^{\pm}(t), \\ \boldsymbol{u}^{\pm}\big|_{t=0} &= \boldsymbol{u}_{0}^{\pm}, && \text{in } \Omega^{\pm}(0), \\ [\nu \partial_{\mathbf{n}}\boldsymbol{u} - p \mathbf{n}]_{\Gamma(t)} &= \boldsymbol{h}_N, \quad [\boldsymbol{u}]_{\Gamma(t)} = \boldsymbol{h}_D, && \text{on } \Gamma(t),\\ \boldsymbol{u}^{+} &= \boldsymbol{g}, && \text{on } \partial\Omega, \end{align}\right.\tag{28}\] where the viscosities \(\nu\) is piece-wise positive constants, which is defined as \[\label{eq4465} \nu(\mathbf{x}, t) = \begin{cases} \nu^+, & \text{for } \boldsymbol{x} \in \Omega^+(t), \\ \nu^-, & \text{for } \boldsymbol{x} \in \Omega^-(t). \end{cases}\tag{29}\]

3.5pt

In this example, since the pressure \(p\) is discontinuous across the interface, XI-PINN requires two separate neural networks to approximate the velocity \(\boldsymbol{u}\) and the pressure \(p\), respectively. For the network approximating the velocity, the extension variable is \(z = \Phi(\boldsymbol{x},t)\), while for the network approximating the pressure we take \(z = \chi(\boldsymbol{x},t)\). To ensure a fair comparison, X-PINNs also use different networks for the velocity and for the pressure; hence both methods involve a total of four neural networks. In both XI-PINNs and X-PINNs, the neural network approximating \(\boldsymbol{u}\) yields two outputs, corresponding to the vector components \((u,v)\). In this experiment the numbers of training points are set to \((N_\Omega, N_{\partial\Omega}, N_{\Omega_0}, N_\Gamma) = (10\text{K}, 3\text{K}, 0.3\text{K}, 2\text{K})\). As before, \(\text{XI-PINN}^*\) refers to the variant where the level set function is constructed by a neural network \(\phi_\theta\); this network has architecture \((L,W)=(2,64)\) and its training takes \(59.7\) seconds.

The decay of the loss and the relative errors during training for the different methods are shown in Fig. 8, and the final numerical results together with the training times are reported in Table 4. The detailed network architectures are listed in Table 5. The numerical behavior is similar to that observed in the scalar case (Example 1): XI-PINN enforces continuity at the interface by construction and therefore achieves better accuracy than X-PINNs. From the error distributions at \(t=0.5\) (Fig. 9) we can see that the error of XI-PINN is more uniformly distributed, whereas X-PINNs exhibit noticeable discrepancies between the interior and exterior of the interface. The cross‑section plots further indicate that such differences may be caused by the error accumulated at the interface. XI-PINN accurately captures the discontinuity in the normal derivative of \(\boldsymbol{u}\) across the interface, and its advantage becomes more pronounced when the size of network is small.

We use this example to illustrate the advantages of neural networks for moving interface problems with large deformations. Since both the velocity \(\boldsymbol{u}\) and the pressure \(p\) are discontinuous across the interface, for XI-PINN a single network suffices, with outputs \((u,v,p)\). For X-PINNs, two networks are needed, one for each subdomain, to approximate the velocity and the pressure. As an analytical level set function is not available in this example, we construct \(\phi_\theta\) using Algorithm 1. Four fully connected networks with architecture \((L,W)=(4,32)\) are employed to build the composite inverse mapping \(\widehat{\mathbf{X}}_\theta\); training these four networks takes \(366.4\) seconds in total. For the X-PINN framework, we similarly utilize the learned level set function \(\phi_\theta\) to classify the collocation points within the global domain \(\Omega\). The numbers of training points used for the PDE solution are the same as in Example 3.

The numerical results are presented in Fig. 10 and Table 6. Here XI-PINN is tested with two network sizes: \((L,W)=(3,32)\) and \((3,64)\); for X-PINNs we use \((3,23)\) and \((3,45)\). With networks of comparable size, both models achieve high accuracy and the numerical results show no significant difference. Because the exact solution \((\boldsymbol{u},p)\) is fairly smooth, both methods reach relative errors on the order of \(10^{-7}\) even with small networks containing just around \(1300\) parameters. The numerical solutions and error distributions of \(\left| \boldsymbol{u}\right|\) at different times are shown in Fig. 11. These results demonstrate that PINN-based methods can handle moving interface problems with large deformations in a simple and effective manner, fully exploiting the advantages of mesh-free approaches.

5 Conclusion↩︎

In this paper, we have presented the XI‑PINN framework for solving parabolic moving interface problems. The method employs a single neural network combined with an extended variable technique driven by a level set function, allowing it to capture the evolving interface and the associated jump conditions without requiring separate networks for each subdomain or body‑fitted mesh generation. By analyzing the error structure, we have decomposed the total error into approximation, statistical, and optimization components, and provided an a priori error estimate that confirms the method’s reliability under appropriate regularity assumptions. 

Extensive numerical experiments, including two‑ and three‑dimensional examples and the Oseen equations, demonstrate the accuracy and robustness of XI‑PINN. The method consistently outperforms both the standard Vanilla‑PINN, which fails to resolve interface discontinuities, and the multi‑network X‑PINN when the solution is continuous across the interface or when network capacities are limited. The single‑network design enables adaptive parameter sharing, leading to a more uniform error distribution and a favorable spectral behavior of the neural tangent kernel, which alleviates the spectral bias commonly observed in physics‑informed neural networks. Moreover, the time‑adaptive neural‑network construction of the level set function successfully handles large interface deformations while preserving topological validity. 

While the present work focuses on moving interface problems with prescribed velocity fields, promising directions for future research include the extension to fully coupled multiphysics problems, such as two‑phase incompressible Navier–Stokes flows [53], where the interface motion is intrinsically linked to the solution fields. In such settings, XI‑PINN’s ability to naturally enforce velocity continuity across the interface will be particularly advantageous, but it also calls for the development of stable and efficient iterative algorithms that alternately update the PDE solution and the evolving interface representation.

References↩︎

---

1. School of Mathematics, Nanjing University, Nanjing 210093, People’s Republic of China (, , ).↩︎

2. Submitted to the editors DATE.↩︎