Physics-informed neural network for the heat equation under imperfect contact conditions and its error analysis

1.   Introduction

Composite media often contain imperfect contacts at the material interfaces, resulting in discontinuous temperatures across the interface. Since the heat transfer and energy dissipation through imperfectly matched interfaces can significantly influence the temperature profile in internal domains ^[1,2,3,4], neglecting the effects of thermal resistance occurring from imperfect contact in modeling can deteriorate the predictions of thermal behavior in the domain of interest. Kapitza's thermal resistance is one of the widely adopted heat transfer models in material science ^[5,6,7], where the amount of jump in the temperature variable along the interface is proportional to the normal heat flux.

One of the difficulties in predicting the thermal profile in imperfectly contacted materials is that the temperature jump across interfaces is implicitly given. In ^[8,9], the jump is treated as independent variables; the jump profiles are updated in an iterative way, where in each iteration the elliptic interface problem is solved with the approximate jump. On the other hand, some people applied finite element method (FEM)-based algorithms, where interface conditions are implicitly contained in bilinear forms ^[10,11]. In addition, there are extended approaches based on the finite element method (XFEM) to solve the imperfect contact problem where additional degrees of freedom appear near the interface to fit implicit jump conditions ^[7,12,13,14]. However, the numerical methods discussed above are inherently complex to implement, as they require additional degrees of freedom in the case of XFEM ^[13,14] or involve intricate modifications of the bilinear form ^[10,11] or demand iterative updates for the jump profile ^[8,9]. Therefore, there arises a need for convenient and efficient methods for the imperfect contact heat transfer equation.

Recently, neural network-based methods have emerged as successful alternatives for solving partial differential equations (PDEs) in the scientific computing community. In the physics-informed neural network (PINN) method, neural network functions serve as a surrogate model for solving PDEs, with its parameters optimized by minimizing objective functions that include the strong form of the governing equations. For an excellent review of PINNs, we refer to ^[15]. With the advancement of highly efficient automatic differentiation techniques such as those provided by PyTorch ^[16], PINNs have been successfully applied to a wide range of problems, including Peridynamics ^[17], thermal-fluid dynamics ^[18,19], electromagnetism ^[20], and Poisson-Boltzmann ^[22].

In this work, we develop a new PINN-based method for the heat equation under the imperfect contact condition. One of the restrictions of conventional PINN in solving PDE lies in the fact that neural networks belong to a continuous family of functions. Therefore, when PDE contains some interface with possibly discontinuous solutions along the interface, it is hard to expect an accurate reconstruction of solutions by PINN. One remedy is to extend the piecewise continuous function to a continuous function of higher dimension. In ^[21,22], the axis augmentation techniques are proposed, where an additional variable labels the locations of subdomains. Following the way of ^[21,22], we first employ continuous Sobolev extensions for the piecewise $H^1$-temperature variable. Then, together with the auxiliary variable that labels the locations of subregions, we develop a continuous representation of the heat variable in an augmented variable. Finally, we define the physics-informed neural network in an augmented variable by the minimizer of the loss functional, which includes the implicit jump condition across the interface. In this way, our method offers a user-friendly way to solve heat transfer equations with imperfect contact conditions, without the need for mesh generation or iterative methods to determine jump profiles across the interface.

Aside from the ease of implementation, another advantage of using a continuous representation of solutions in augmented variables is that it allows error analysis in the space of smooth functions. We provide an error analysis of the proposed method. Our main theorem (see Theorem 3.3) states that the difference between the exact solution and the neural network approximation in an energy-like norm is bounded by the physics-informed loss functional. Here, the loss functional can be made arbitrarily small by increasing the parameters in the neural network, such as the number of nodes and number of training samples (see Theorem 3.6). The numerical tests in the result section support the error analysis. To the best of the authors' knowledge, the proposed method is the first successful application of PINN for solving interface problems with implicit jump conditions, accompanied by an error analysis. Let us summarize the advantages of the proposed work: 1) by establishing the continuous representation of solutions in the augmented variable, our method benefits from the universal approximation property of the neural networks; 2) since our methods are formulated on smooth function space, we were able to derive the error estimates for the proposed PINN; 3) our method does not require mesh generation or linearization procedures, in contrast to the conventional FEM-based method.

The rest of the paper is organized as follows. In the next section, we describe the governing equation and develop our version of PINN. The error estimates are carried out in Section 3. In Section 4, we document the performance of the proposed algorithm. The conclusion follows in Section 5.

2.   Method

In this section, we develop a PINN method for a heat equation involving imperfect contact. We employ continuous Sobolev extensions for the temperature variable to represent the solution as a continuous function in an augmented variable. The model equation is described in Subsection 2.1, and our version of the PINN method is proposed in Subsection 2.2.

We introduce some definitions. Suppose $G$ is any domain in $\mathbb{R}^d$. For integer $d \geq 1$, $s\geq 0$, and $1 \leq p \leq \infty$, we denote by $W^{s, p}(G)$ the usual Sobolev space associated with norm $\|\cdot\|_{W^{s, p}(G)}$ and seminorm $|\cdot|_{W^{s, p}(G)}$. In particular, if $p = 2$, we then write $W^{s, 2}(G) = H^s(G)$, with $\| {\cdot} \|_{W^{s, 2}(G)} = \| {\cdot} \|_{s, G}$ and $|\cdot|_{W^{s, 2}(G)} = |\cdot|_{s, G}$. We also write $W^{0, p}(G) = L^p(G)$ and $\| {\cdot} \|_{W^{0, p}(G)} = \| {\cdot} \|_{L^p(G)}$. For the inner product on $G$, we use the notation $(\cdot, \cdot)_G$. Given an integer $k \geq 0$, let $C^k(G)$ denote the space of all $C^k$ functions in the closure of $G$. Given a real-valued function $\psi$, its support is denoted by $\text{supp}(\psi)$. We denote by $\chi_G$ the indicator function of $G$.

2.1.   Heat equation under imperfect contact condition

We consider a (bounded) composite material domain $\Omega \subset \mathbb{R}^d$ ($d = 2, 3$) divided by material interface $\Gamma$, i.e., $\Omega = \Omega^+\cup\Omega^-\cup\Gamma$ (see Figure 1). The heat flux variable is defined as $\mathbf{q} = -k\nabla T$, where $T$ and $k$ are the temperature variable and thermal conductivity constant, respectively. In an imperfect contact situation, the temperature jump across the material interface is proportional to the heat flux in the normal direction of the interface:

where $\alpha$ the is Kapitza constant ^[5,6,7]. Now, let us summarize the steady-state heat equation on $\Omega$ under the imperfect condition as follows:

where $f$ is the source term and $\mathbf{n}_\Gamma = \mathbf{n}_{\Omega^-}$. Also, $[\cdot]_\Gamma$ denotes a jump operator across the interface, i.e.,

Here, we assume that $k$ is piecewise constant, that is, $k = k^+\chi_{\Omega^+} + k^-\chi_{\Omega^-}$ for some positive constants $k^+$ and $k^-$. The problem (2.1) has a unique solution ^[23,24].

2.2.   Physics-informed neural network method

In this subsection, we describe the PINN method for the governing equation. The overall process is illustrated in Figure 2.

One challenge is that the solution and its derivative are discontinuous across $\Gamma$, whereas conventional neural network functions are inherently smooth on the entire domain. As a result, traditional neural network functions may be inefficient in capturing the discontinuities in the governing equation. To remedy this, we adopt the axis augmentation technique introduced in ^[21]. The key idea is to introduce an augmented variable $z$ to extend the discontinuous function $T$ to an augmented function $T^{aug}$ that is continuous over $\Omega\times[-1, 1]$. Let us assume that $T|_{\Omega^-}$ and $T|_{\Omega^+}$ are sufficiently smooth. We apply the Sobolev extension theorem (see, e.g., ^[25,26]), which guarantees the existence of smooth extensions $T^s$, ($s = \pm$), satisfying

With the additional variable $z$, we define axis-augmented function

Then, it follows that

In this way, we have derived a continuous function $T^{aug}$ in an augmented variable, where the variable $z$ labels the locations of subregions, i.e., $z = 1$ for $\Omega^+$ and $z = -1$ for $\Omega^-$.

The advantage of introducing such a continuous representation is that it enables the design of neural network surrogate functions in the augmented variable space, which retains the universal approximation property ^[27]. The neural network functions in the augmented variable are defined as follows: Given integer $n\geq 1$, we introduce the set $\mathcal{N}_n^{aug}$ of all augmented neural network functions $\phi_\theta^{aug}$ having $n$ hidden layers:

where $L_i$'s are linear affine functions and $\sigma$ is a tanh activation function. Here, $\theta$ denotes the collection of all training parameters. The advantage of defining such an augmented neural network function is that, by restricting $z$ to subdomains, we can obtain the piecewise continuous function, i.e.,

which can efficiently approximate the desired piecewise continuous solution $T$. It remains to impose loss functionals for the functions in

Since the residual of the governing equation can be evaluated on discrete points, we need to define the collection of training points in the variable $(\mathbf{x}, z)$. Let $\{(\mathbf{x}_i^+, 1)\}_{i = 1}^{M^+}$ and $\{(\mathbf{x}_i^-, -1)\}_{i = 1}^{M^-}$ be the collections of training points in $\Omega^{+}\times\{+1\}$ and $\Omega^{-}\times\{-1\}$, respectively. Let $\{ \mathbf{x}_i^{\Gamma}\}_{i = 1}^{M^{\Gamma}}$ be the collection of training points on $\Gamma$. We also let $\{ \mathbf{x}_i^{\partial}\}_{i = 1}^{M^{\partial}}$ be the collection of points on $\partial\Omega$. Some choices of the training points can be found in ^[21,22]. Motivated by the governing Eqs. (2.1a)–(2.1d), the loss functionals are defined as follows: For any piecewise smooth function $\phi_\theta \in \mathcal{N}_n^\pm$,

The detailed calculations of loss functionals are provided in Appendix A. Finally, our version of the PINN method for (2.1) is formulated as follows: Find $T_{\theta}\in\mathcal{N}_n^{\pm}$ such that

Once $T_\theta$ is obtained, the function inherently satisfies the implicit jump conditions (2.1b)–(2.1c). Consequently, we introduce an efficient and user-friendly method to solve the imperfect contact problem without requiring mesh generation or iterative procedures to determine the jump profile.

3.   Error estimates

In this section, we derive the error estimates of the proposed PINN method. We define the residuals $\mathcal{R}^\Omega, \mathcal{R}^\Gamma$, and $\mathcal{R}^\partial$ as follows: Given a real-valued function $\phi$ on $\Omega$,

Also, let $\mathcal{R} : = \mathcal{R}^\Omega + \mathcal{R}^{\Gamma} + \mathcal{R}^\partial$. Then the main result in this section is summarized as follows.

● (Theorem 3.3) If $T$ is the solution of (2.1) which is sufficiently smooth in $\Omega^+$ and $\Omega^-$, then for any $T_{\theta}\in \mathcal{N}_n^{\pm}$, the piecewise $H^1$-error between $T$ and $T_{\theta}$ can be bounded above by the residual $\mathcal{R}(T_{\theta})$.

● (Theorem 3.6) If $T$ is sufficiently smooth in $\Omega^+$ and $\Omega^-$, then the residual can be made arbitrarily small. For any $\varepsilon > 0$ there exists $T_{\theta}\in \mathcal{N}_n^{\pm}$ such that $\mathcal{R}(T_{\theta}) < \varepsilon$.

Note that, for a sufficiently large number of training points in the definition of the loss functional $\mathcal{L}$, the errors between the loss functional $\mathcal{L}$ and the residual $\mathcal{R}$ can be made arbitrarily small (see, e.g., ^{[28,29,30,31]}). Therefore, we assume that the following holds.

Assumption 3.1. It holds that $\mathcal{L}^s(\phi) \approx \mathcal{R}^s(\phi)$ for $s = \Omega, \Gamma, \partial$.

This assumption, together with the results in Theorem 3.3 and Theorem 3.6, implies that the error between the exact solution and the neural network solution can be made arbitrarily small by controlling the parameters in the neural network.

Finally, we introduce some broken function spaces

and associated norms (respectively)

3.1.   Error bound by residuals

To prove Theorem 3.3, we need the Poincaŕe inequality for the piecewise $H^1$-functions which can be found in [22, Lemma 1].

Lemma 3.2. It holds that

where $C$ is a positive constant depending only on $\Omega$ and $\Gamma$.

Theorem 3.3. Suppose that the solution $T$ of (2.1) satisfies $T\in H^2(\Omega^\pm)\cap C^1(\Omega^\pm)$. Then we have

where $C$ is a positive constant depending only on $\Omega$, $\Gamma$, $k$, $\alpha$, $\| { T} \|_{C^1(\Omega^{\pm})}$ and $\| { T_{\theta}} \|_{C^1(\Omega^{\pm})}$.

Proof. Let $v = T - T_\theta$. By Lemma 3.2 and by the fact that $k^-, k^+$, and $\alpha$ are positive constants, we have

where $\widetilde{C} = C_0 {\rm max} \{(k^-)^{-1}, (k^+)^{-1}, \alpha^{-1}\}$. By the boundary condition (2.1d), we have

Using integration by parts and the fact that $- k\Delta T = f$ on each $\Omega^-\cup\Omega^+$, we have

We remark that $\mathbf{n}_\Gamma = \mathbf{n}_{\Omega^-} = - \mathbf{n}_{\Omega^+}$. Using the identity that

we can rewrite the last three terms of (3.2) as

Here, by the interface condition (2.1b)–(2.1c), we have

Therefore, (3.3) becomes

By (3.2)–(3.5), we obtain

We estimate each of the terms. By the Cauchy-Schwarz and Young's inequalities, we have for $I_1$, $I_2$,

Here, $\delta_1$ and $\delta_2$ are positive constants related to the Young's inequality, which will be determined later.

Next, let us estimate $I_3$. By the fact that $T, T_\theta \in C^1(\Omega^\pm)$ and $T = g$ on $\partial\Omega$, we have

By the Cauchy-Schwarz and Young's inequalities (with $\delta_3, \, \delta_4 > 0$), we have, for $I_4$ and $I_5$,

Since $T, T_\theta \in C^1(\Omega^\pm)$ and $[k \nabla \cdot \mathbf{n}_\Gamma] = 0$, we have, for $I_6$,

Then, plugging the estimates of $I_1$, $\cdots$, $I_6$ into (3.6), together with (3.1), we have

Now, by the choice

we obtain the desired inequality. □

Remark 3.4. As stated in Theorem 3.3, the constant $C$ may depend on the piecewise $C^1$-norm of the neural network $T_{\theta}$. In practice, to prevent this norm from becoming excessively large, we can include an $L^2$-regularization term in the loss functional $\mathcal{L}$, following the approach in ^[29,31].

3.2.   Estimates on the residual and the training error

As proved in the previous subsection, the piecewise $H^1$-error can be controlled by the residuals $\mathcal{R}^{\Omega}$, $\mathcal{R}^{\partial}$ and $\mathcal{R}^{\Gamma}$. Given some regularity assumption, we will prove that if $n > 1$, then the residuals can be made arbitrarily small, i.e., for any $\varepsilon > 0$, there exists $T_{\theta}\in\mathcal{N}_n^{\pm}$ such that $\mathcal{R}(T_{\theta}) < \varepsilon$.

Lemma 3.5. Suppose that $n > 1$ and the solution $T$ of (2.1) satisfies $T|_{\Omega^s}\in W^{3, \infty}(\Omega^s)$ for $s = \pm$. For any $\varepsilon > 0$, there exists a neural network $T_{\theta}^{\pm}\in \mathcal{N}_n^{\pm}$ such that

Proof. By the argument given in the beginning of Subsection 2.2, for sufficiently large $R > 0$, one can construct an extension $T^{aug}\in W^{3, \infty}([-R, R]^{d+1})$ of $T$ such that $\Omega\times[-1, 1]\subset [-R, R]^{d+1}$ and

Then, according to the approximation theorem of neural network functions with the tanh activation function ^[32], there exists $T_{\theta}^{aug} \in \mathcal{N}_n^{aug}$ such that

Now the conclusion follows from the Sobolev embedding theorem (cf. ^[26]). □

Using the lemma above, we obtain the desired theorem.

Theorem 3.6. Suppose that $n > 1$ and the solution $T$ of (2.1) satisfies $T|_{\Omega^s}\in W^{3, \infty}(\Omega^s)$ for $s = \pm$. For any $\varepsilon > 0$, there exists $T_{\theta}\in \mathcal{N}_n^{\pm}$ such that

where $C$ is a positive constant depending only on $\Omega$, $\Gamma$, $k$, $\alpha$ and $T$.

Proof. Let $\varepsilon > 0$ be chosen arbitrarily. By Lemma 3.5, there exists $T_{\theta}\in\mathcal{N}_n^{\pm}$ such that

Therefore, we have

Since $T$ is the solution of (2.1), we obtain from (3.7) and the trace inequality that

This completes the proof. □

Remark 3.7. In the proof of Lemma 3.5, the approximation property of the tanh neural network functions ^[32] was employed. This validates our choice of tanh activation functions in (2.3).

4.   Numerical results

In this section, we present some examples. The domain $\Omega = [-1, 1]^2$ is separated by zeros of a level set function $L(x, y)$, i.e., $\Omega^- = \{(x, y)\in \Omega : L(x, y) < 0\}$ and $\Omega^+ = \{(x, y)\in \Omega : L(x, y) > 0\}$. For all examples, we assume that $n = 2$ (that is, the number of hidden layers is two) and the tangent hyperbolic function is adopted for the activation function. In addition, we assign $N$-number of nodes for each hidden layer. The weights in the neural network are initialized using the Kaiming uniform method ^[33], with biases set to zero.

We present three examples, each featuring a different interface shape: circle, line, and perturbed circle. In Example 4.1–4.2, the $L^2$- and $L^\infty$-errors against the exact solutions are computed, where we observe reasonable convergence as the parameters increase. In Example 4.3, a unidirectional heat flux generated by some boundary condition is considered. Since the analytic solution is unknown in this case, we validate the PINN-predicted solution by comparing it with the FEM solution. All examples are performed on a single core of an Intel(R) Core(TM) i7-13700 CPU.

4.1.   Choices of optimizers and training sample sets

In this subsection, we justify our choice of optimizer and the training sampling strategy. We compare the results obtained by two optimizers in Example 4.1. The first is the Levenberg–Marquardt (LM) optimizer ^[34] with 3,000 epochs. The LM optimizer is configured with an initial damping parameter of $10^5$, which is updated at every iteration. The damping parameter is divided by 1.3 when the process is successful and is multiplied by 3 otherwise. The second is the Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimizer ^[35] with 10,000 epochs, a learning rate of 0.1, and a strong Wolfe line search option.

To describe the training points, let us introduce the order-related parameter $m$, which determines the total number of training points. The inner points ($m^2$-number) are determined either by Chebyshev–Gauss points, Gauss–Legendre points, or uniformly spaced grid points. For example, Chebyshev–Gauss points of order $m$ is defined by

The interface ($4m$-number) and boundary points ($4m$-number) are uniformly sampled at the interface and boundary, respectively. In this way, the total number of training points corresponding to the parameter $m$ is $m^2+8m$. We compare the results obtained by different types of training points in Example 4.1.

Example 4.1 (Circle shaped interface). In this example, the level set function is $L(x, y) = x^2+y^2-r_0^2$ and the exact solution is

where $r_0 = 0.5$ and $\alpha = 1$. We consider two different thermal conductivity contrasts: $(k^-, k^+) = (1, 10)$, or $(100, 1)$.

We compare loss curves with increasing training epochs for the LM and L-BFGS optimizer in Figure 3. The loss function decreases to the lower level by the LM optimizer. Consequently, we adopted the LM optimizer for all experiments in this study. Next, we compare the $L^2$ and $L^\infty$ errors computed using Chebyshev–Gauss, Gauss–Legendre, and uniform grids type training points. We report errors in Figure 4. Overall, the smallest errors are achieved with the choice of Chebyshev–Gauss points and with $N = 20$. Therefore, we use Chebysev–Gauss type training points for the rest of the Section.

Now, let us present numerical results by the proposed PINN method. Figure 5 presents the numerical solutions for the two conductivity cases. Notably, no spurious oscillation appears near the interface in any case. We report the errors $L^2$ and $L^\infty$ that are computed at 10,000 uniformly sampled test points. The errors and CPU time for the two conductivity contrasts, that is, $(k^+, k^-) = (1, 10)$ or $(100, 1)$, are documented in Tables 1 and 2, respectively.

Let us discuss the computational complexities with respect to $N$ and $m$. One of the bottlenecks of the LM optimizer is the inversion of the Hessian matrix associated with the loss functional. In our case, the number of entries in the Hessian matrix is $\mathcal{O}(N^2m^2)$. As a results, the CPU time may dramatically increases when $N$ and $m$ increase. Fortunately, the smallest errors were achieved at relatively small parameters, specifically $N = 20$ and $m = 32$. In particular, with this parameter choice, the $L^2$ errors remain below $4\cdot 10^{-7}$ in Tables 1 and 2 with corresponding CPU times of 261 (s) and 274 (s), respectively.

Before closing this example, we compare the proposed method with the FEM-type method. One of the authors has the experience of implementing FEM for elliptic interface problems with implicit jump conditions, as described in ^[10], where the $L^2$ error is of $\mathcal{O}(h^{-2})$. To achieve comparable accuracy, the FEM approach requires solving an algebraic system with over $10^5$ unknowns, resulting in substantial computational cost (see e.g., Table 2 in ^[10]). Although it is difficult to compare directly, our method achieves a comparable accuracy with less CPU time.

4.2.   Comparison with other methods

In this subsection, we compare the proposed PINN with other methods. In Example 4.2, we compare our method with the piecewise PINN method ^[36], while in Example 4.3, we compare it with a FEM-type method ^[11].

Example 4.2 (Line interface). The level set function is $L(x, y) = x+2y+r_0$ and the exact solution is

where $\alpha = 0.5$ and $r_0 = 0.1$. The thermal conductivity is $(k^-, k^+) = (1, 10)$. The graphs of $T_\theta$ and $|T-T_\theta|$ are presented in Figure 6. The $L^2$ and $L^\infty$ errors are reported on the left side of Table 3. Notably, the $L^2$ errors are below $2\times 10^{-7}$ when $N = 10$ and $m = 16$. Therefore, the proposed PINN successfully approximates the exact solution.

We compare our version of PINN with piecewise PINN introduced in ^[36]. In ^[36], the authors propose two separate neural network functions for each subdomain. As a result, the total number of parameters is approximately doubled when the same number of nodes is used in each hidden layer. We report the $L^2$ and $L^\infty$ errors and CPU time on the right side of Table 3. We observe that $L^2$ and $L^\infty$ errors obtained by the proposed PINN and piecewise PINN are similar. However, in terms of CPU time, our PINN is more efficient, which is due to the increased number of parameters in piecewise PINN.

Example 4.3 (Perturbed-circle shaped interface). In this example, we consider the unidirectional heat flow given by the homogeneous outer source ($f = 0$) and the boundary condition that

where $(k^-, k^+) = (1, 2)$ and $\alpha = 0.5$. We consider the perturbed circle shape interface:

Since the exact solution is unknown, we compare the solution obtained by PINN with that by FEM ^[11], denoted as $T_{FEM}$. In Appendix B, we provide a brief description of FEM for the completeness of the presentation.

Figure 7 presents the comparison between $T_{FEM}$ and $T_\theta$, showing the overall agreement between the two solutions. We also calculate $L^2$ and $L^\infty$ differences between $T_{\theta}$ and $T_{FEM}$:

which are reasonably small. Thus, we conclude that the solutions obtained by PINN are comparable to those from the well-established FEM. However, in terms of implementation, the proposed method is simpler than FEM type method, without the needs for the mesh generation or iteration techniques to determine jump profiles.

4.3.   Discussion

In this subsection, we provide some discussion regarding the numerical experiments. We begin by examining the relation between the choice of the parameter $N$ (number of nodes in a hidden layer) and its influences on the accuracy. In Figure 8, we plot the $L^2$ errors of the reconstructed solutions obtained by our methods with varying number of $N$, with $m = 16$ fixed both for the Example 4.1 and 4.2. Here, we do not observe clear decaying trends in errors as $N$ increases. This could result from difficulties in finding the global minimizer using optimization algorithms, as discussed in ^[37,38]. The convergence of the optimizer, closly related to the optimization landscape, depends on various factors including the interface shapes, PDE-related parameters (e.g., $k^\pm$), etc. Notably, in Example 4.2, where both the interface shape and the solutions are relatively simple, the smallest $L^2$ error ($2.3\text{E}-7$) was observed with just $N = 10$. A similar phenomenon is observed for the $m$ parameter, which determines the number of training points. In Figure 4, we do not observe the clear decay of errors as $m$ increases. Based on our parameter study, we suggest selecting moderate values such as $N = 10, 20$ and $m = 16$ rather than significantly large values.

We now turn to the discussion of the stability of the optimization process. Through our experiments, we observed that the optimization process becomes unstable when there is a large contrast of $k^+$ and $k^-$ or the parameters $N$ and $m$ are too large. In particular, when the damping parameter in LM becomes close to zero, the ill-conditioned Hessian matrix may cause the loss function to diverge. To remedy this, we modified the LM optimizer slightly by introducing an additional rule: if current loss exceeds twice the previous loss, increase the damping parameter by a factor of $10^5$. This adjustment improved the robustness of the optimization process.

Finally, let us discuss the regularity assumptions used in our error estimates. In the statements of Lemma 3.5 and Theorem 3.6, we assumed that $T|_{\Omega^s} \in W^{3, \infty}(\Omega^s)$ for $s = \pm$, which is a strong regularity condition. However, this assumption may not hold in practical applications, especially when material interfaces induce singularities. In such cases, the errors may not be controlled by increasing the number of parameters since the universal approximation fails under low regularity conditions.

5.   Conclusions

In this work, we developed a new PINN method to solve heat transfer equations involving imperfect contact conditions. There are mainly two difficulties in solving imperfect contact problems: 1) the solutions are discontinuous across the interfaces, and 2) the jumps are implicitly determined. To remedy these, we adopted the axis-augmentation techniques introduced in ^[21] to establish a continuous representation of the solution in the augmented variable. Then, we defined the neural network function in an augmented variable with the physics-informed loss functional, which includes the implicit jump conditions. After the loss function is minimized, the resulting solutions obtained by PINN naturally satisfy the desired jump conditions. One of the advantages of such a continuous representation formulation was the availability of error analysis. We carried out the error analysis; our main theorem states that given a sufficient number of parameters, the loss between exact and PINN solutions can be made arbitrarily small. The numerical experiments that support our analysis are reported in the result section. Let us discuss the limitations of the proposed method. First, the number of parameters in the neural network increases with the introduction of a new axis. Second, while our analysis assumes strong regularity conditions, these conditions may not always hold in the real-world applications.

Appendix

A.   Calculation of loss functional

We provide details regarding the calculation of loss functionals presented in (2.5)–(2.9). Recall first that $\phi_{\theta}(\mathbf{x}) = \phi_{\theta}^{aug}(\mathbf{x}, 1)\chi_{\Omega^+}(\mathbf{x}) + \phi_{\theta}^{aug}(\mathbf{x}, -1)\chi_{\Omega^-}(\mathbf{x})$. While the loss functionals are written in the function $\phi_\theta$, the augmented neural network function $\phi_\theta^{aug}(\mathbf{x}, z)$ is employed during the implementation stage. Specifically, the parameter $z\in \{1, -1 \}$, indicates whether the sampling points belong to $\Omega^-$ or $\Omega^+$, respectively. Below, we summarize how operators such as $\nabla$, $\nabla$, and $[\cdot]$ acting on $\phi_\theta$ are calculated using the augmented neural network functions:

where

B.   FEM for heat equation involving imperfect contact condition

We briefly describe the FEM introduced in ^[11], where the so-called immersed finite element (IFE) method was employed. In IFE methodology, uniform grids are employed for the interface problem, allowing the interface to cut through the element, which contrasts with the usage of interface-fitted grids. Instead, the basis function is modified so that the they satisfy the interface conditions.

Let $\mathcal{T}_h$ be a triangulation of the domain $\Omega$. When $E\in\mathcal{T}_h$ is not cut through by the interface, the local space is a standard linear space, denoted by $S_h(E)$. Suppose $E$ is cut through by the interface; the basis function is defined as a piecewise polynomial of the form:

Here, the coefficients in (B.1) are determined by the interface conditions (2.1b)-(2.1c), i.e.,

In this way, the local space $\widehat{S}_h(E)$ is defined with nodes and degrees of freedom. The global space $\widehat{S}_h(\Omega)$ is defined by patching local spaces, imposing the continuity at the nodes. Let us remark that the dimension of $\widehat{S}_h(\Omega)$ corresponds to the number of nodes in $\mathcal{T}_h$. Finally, FEM solution is obtained by solving the weak problem : find $T_h \in \widehat{S}_h$ such that satisfies

where

The detailed implementation can be found in ^[11].

Author contributions

Hansaem Oh: Formal analysis, Methodology, Writing-original draft; Gwanghyun Jo: Conceptualization, Methodology, Writing-review & editing.

Acknowledgements

The first author (H. Oh) is supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NO. RS-2024-00463063). We thank Dr. Hyeokjoo Park for his valuable advice on the error analysis.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

All authors declare no conflicts of interest in this paper.