HomPINNs: homotopy physics-informed neural networks for solving the
  inverse problems of nonlinear differential equations with multiple solutions
Due to the complex behavior arising from non-uniqueness, symmetry, and
bifurcations in the solution space, solving inverse problems of nonlinear
differential equations (DEs) with multiple solutions is a challenging task. To
address this, we propose homotopy physics-informed neural networks (HomPINNs),
a novel framework that leverages homotopy continuation and neural networks
(NNs) to solve inverse problems. The proposed framework begins with the use of
NNs to simultaneously approximate unlabeled observations across diverse
solutions while adhering to DE constraints. Through homotopy continuation, the
proposed method solves the inverse problem by tracing the observations and
identifying multiple solutions. The experiments involve testing the performance
of the proposed method on one-dimensional DEs and applying it to solve a
two-dimensional Gray-Scott simulation. Our findings demonstrate that the
proposed method is scalable and adaptable, providing an effective solution for
solving DEs with multiple solutions and unknown parameters. Moreover, it has
significant potential for various applications in scientific computing, such as
modeling complex systems and solving inverse problems in physics, chemistry,
biology, etc.
Haoyang Zheng, Yao Huang, Ziyang Huang, Wenrui Hao, Guang Lin
HomPINNs: homotopy physics-informed neural networks for solving the inverse problems of nonlinear differential equations with multiple solutions

###### Abstract

Due to the complex behavior arising from non-uniqueness, symmetry, and bifurcations in the solution space, solving inverse problems of nonlinear differential equations (DEs) with multiple solutions is a challenging task. To address this issue, we propose homotopy physics-informed neural networks (HomPINNs), a novel framework that leverages homotopy continuation and neural networks (NNs) to solve inverse problems. The proposed framework begins with the use of a NN to simultaneously approximate known observations and conform to the constraints of DEs. By utilizing the homotopy continuation method, the approximation traces the observations to identify multiple solutions and solve the inverse problem. The experiments involve testing the performance of the proposed method on one-dimensional DEs and applying it to solve a two-dimensional Gray-Scott simulation. Our findings demonstrate that the proposed method is scalable and adaptable, providing an effective solution for solving DEs with multiple solutions and unknown parameters. Moreover, it has significant potential for various applications in scientific computing, such as modeling complex systems and solving inverse problems in physics, chemistry, biology, etc.

keywords: Machine learning

Physics-informed neural networks

Nonlinear DEs

Multiple solutions

Homotopy continuation method

## 1 Introduction

Inverse problems are popular in various fields of computational physics, such as fluid dynamics [1], electromagnetics [2], and geophysics [3]. These problems often involve recovering a set of parameters given certain observations or measurements. Traditional techniques for solving inverse problems include optimization-based methods, such as adjoint methods [4], Bayesian inversion [5], and sampling-based methods [6], etc. While the above-mentioned approaches have proven successful in specific research areas, they come with several obstacles, such as high computational costs, a need for significant domain knowledge, or issues with ill-posed problems.

In the last few years, deep learning techniques, specifically physics-informed neural networks (PINNs) [7; 8], have shown potential for solving inverse problems more accurately and efficiently. By integrating prior physical knowledge into neural network (NN) architectures, PINNs enable to solve inverse problems with superior generalization ability and a reduced need for the number of observations. The concept of PINNs was first introduced by Raissi _et al._[8] in a series of papers that demonstrated their potential for solving forward and inverse problems in computational physics. The intuition behind PINNs is the incorporation of prior physical knowledge into the loss function. By minimizing this loss function, the NN learns to approximate the solution to the inverse problem while satisfying the
underlying physics. Following the idea of PINNs and with the help of automatic differentiation, a Python library [9] was developed by Lu _et al._ to solve inverse problems in the area of computational science and engineering.

Since then, PINNs and their variants have been successfully applied to solve inverse problems in many different areas. Building upon PINNs, Pang _et al._[10] developed fractional PINNs (fPINNs) for solving space-time fractional advection-diffusion equations from scattered and noisy data. Kharazmi _et al._[11] proposed hp-Variational PINNs (hp-VPINNs) with domain decomposition to construct both local and global neural network approximations. Liu _et al._[12] introduced Bayesian PINNs (B-PINNs) for solving partial differential equations (PDEs), providing accurate predictions and quantifying aleatoric uncertainty with noisy data. Jagtap _et al._[13] tried to enforce flux continuity and average solution at sub-domain interfaces and proposed conservative PINNs (cPINNs) on discrete domains for nonlinear conservation laws. Jagtap _et al._[14] extended a generalized space-time domain decomposition approach, extended PINNs (XPINNs), for solving nonlinear PDEs on complex-geometry domains. Wang _et al._[15] investigate the limitations of PINNs in approximating high-frequency or multi-scale functions, and propose novel architectures using spatio-temporal and multi-scale random Fourier features to improve the robustness and accuracy of PINNs. Wang _et al._[16] also investigated the training behavior and failures of PINNs using the Neural Tangent Kernel (NTK) framework, identifying a discrepancy in the convergence rate of different loss components, and proposed a novel gradient descent algorithm that adaptively calibrates the convergence rate. Additionally, PINNs have been implemented in various research domains, including fluid dynamics [17; 18], solid mechanics [19], molecular dynamics[20], quantum chemistry [21], etc. More papers related to solve inverse problems with the use of PINNs can be found in [9; 22].

Despite the advancements and vast potential of PINNs for solving inverse problems demonstrated in the aforementioned works, they may not be well-suited for addressing differential equations (DEs) with multiple solutions since they rely on gradient descent optimization algorithms that converge to a single local minimum. In problems with multiple solutions, PINNs may only converge to one of the solutions, but they fail to capture other solutions. Without any prior knowledge of multiple solutions of DEs, PINNs may get trapped in a local minimum instead of exploring other minima and resulting in failure to solve inverse problems.

In this work, we consider the combination of PINNs and homotopy continuation methods to address this issue by enabling the exploration of multiple solutions for DEs. Homotopy continuation [23] is a numerical technique that traces the solution path of a given problem as the parameter changes from an initial to a final value. Hao _et al._[24] utilized the WENO scheme in homotopy continuation to resolve steady-state problems associated with hyperbolic conservation laws. Hao _et al._[25] presented a homotopy continuation-based method using domain decomposition to solve large polynomial systems resulting from discretizing nonlinear algebraic DEs. Hao _et al._[26] proposed an adaptive step-size homotopy tracking method for computing bifurcation points in nonlinear systems. Hao _et al._[27] developed a novel approach to compute multiple solutions of nonlinear DEs by adaptively constructing a spectral approximation space using a greedy algorithm. These homotopy methods have been successfully applied to solving pattern formation problems arising from math biology [28].

In this context, we proposed homotopy physics-informed neural networks (HomPINNs), where the homotopy continuation method is incorporated into the PINNs to find multiple solutions by gradually transforming the problem from a simple problem with given conditions to the target deformed problem of interest. The transformation process is accomplished by adjusting the scale of constraints in the loss function during NN training: In the training process, we consider two constraints in the loss function: one term aims to minimize the distance between the NN approximation and the given observations, while the other term enforces the approximation's conformity with the governing equations. We first ensure that the NN approximation satisfies both two constraints. As training proceeds, we adjust the magnitude of a homotopy tracking parameter to gradually focus more on the first constraint. As the homotopy tracking parameter changes, both DEs with multiple solutions and the inverse problems can be solved. This strategy makes it possible to explore the solution space more effectively and identify various local minima corresponding to different solutions of the DEs.

In contrast to traditional methods used for solving inverse problems of DEs with multiple solutions, our proposed method is advantageous over previous methods in the following aspects:

1. The proposed method provides both adaptable and flexible solutions for inverse problems.
2. The proposed method is ideal for inverse problems that are high-dimensional or computationally intense.
3. The proposed method is able to solve inverse problems even when only a limited number of DE solutions are available in the observations.

The remainder of this paper is structured as follows: The problem description of the inverse problems and the proposed HomPINN method are presented in Section 2, along with a comprehensive explanation of the construction
details. The effectiveness of our proposed HomPINNs is demonstrated through several numerical examples in Section 3. The concluding remarks and discussion of our proposed method are presented in Section 4.

## 2 The proposed framework

In this section, we first define the inverse problems associated with DEs. Next, we detail the technique we proposed to identify the correct parameters within DEs whose solution is not unique, based on observations \(\{\mathbf{x}_{i},u_{i}\}_{i=1}^{N_{o}}\). As the observations contain different solutions of DEs and are unclassified, it is imperative to discuss how to classify those observations to different solutions of DEs. We finally delve further into inferring unknown parameters in DEs from the given observations.

### Problem description

We consider the following nonlinear DE:

\[\begin{cases}\mathcal{L}u\left(\mathbf{x}\right)=f\left(u,\mathbf{x};\mathbf{\lambda} \right),\ \ \mathbf{x}\in\Omega,\\ \mathcal{B}u\left(\mathbf{x}\right)=b\left(\mathbf{x}\right),\ \ \mathbf{x}\in\partial \Omega.\end{cases} \tag{1}\]

Here \(\mathcal{L}\) is a linear differential operator, \(f(u,\mathbf{x};\mathbf{\lambda})\) is a nonlinear forcing term with an unknown parameter vector \(\mathbf{\lambda}\), solution \(u\), and \(\Omega\in\mathbb{R}^{d}\) (\(d\) denotes the dimension of the input \(\mathbf{x}\)) is the domain of interest. \(\mathcal{B}\) is the operator for the boundary condition, and \(b\left(\mathbf{x}\right)\) is the force function at the boundary. Even under the homogeneous boundary condition, i.e., \(b(\mathbf{x})=0\), it is still possible that the solution of (1) is not unique and that the exact number of solutions is not known. Given the observation set, \(\{\mathbf{x}_{i},u_{i}\}_{i=1}^{N_{o}}\), where \(u_{i}\) is randomly observed from one of the solutions satisfying (1) and \(N_{o}\) is the total number of observations, the present study aims to determine the value of \(\mathbf{\lambda}\) in (1).

### Neural networks to solve differential equations with multiple solutions

By leveraging the observations \(\{\mathbf{x}_{i},u_{i}\}_{i=1}^{N_{o}}\), collocation points \(\{\mathbf{x}_{j}\}_{j=1}^{N_{o}}\) (\(N_{c}\) is the total number of collocation points), and governing equation (1), we are going a PINN to find the multiple solutions of (1). A general framework of the NN can be found in Fig. 1. Referring to [8], we are constructing a NN that can yield an approximation \(\hat{\mathbf{u}}(\mathbf{x}_{i};\mathbf{\theta},\mathbf{\lambda})\) from input \(\mathbf{x}_{i}\), network parameters \(\mathbf{\theta}\), and an unknown parameter vector \(\mathbf{\lambda}\). Note that \(\hat{\mathbf{u}}(\mathbf{x}_{i};\mathbf{\theta},\mathbf{\lambda})\) is actually a concatenation of \(\hat{u}_{m}(\mathbf{x}_{i};\mathbf{\theta},\mathbf{\lambda})\) (\(1\leq m\leq M\)), where \(\hat{u}_{m}(\mathbf{x}_{i};\mathbf{\theta},\mathbf{\lambda})\) attempts to approximate one of the solutions of (1), and \(M\) is an estimated total number of solutions for (1).

**Remark 1**.: _One concern is that the estimated \(M\) does not match the true quantity of solutions for (1). We will explore how to accurately identify its correct value in the experimental parts (see 3.1 and 3.2)._

Figure 1: Schematic of one PINN for the inverse problems: given the observations \(\{\mathbf{x}_{i},u_{i}\}_{i=1}^{N_{o}}\) and collocation points \(\{\mathbf{x}_{j}\}_{j=1}^{N_{o}}\). PINN makes prediction \(\hat{\mathbf{u}}\), where \(\hat{\mathbf{u}}\) is a concatenation of \(\hat{u}_{m}\) (\(m=1,2,\cdots,M\)). Here \(M\) is an estimated total number of solutions of (1), \(\hat{u}_{m}(\mathbf{x}_{i};\mathbf{\theta},\mathbf{\lambda})\) approximates of one of the solutions of (1). Two constraints (purple rectangular boxes) will be considered here to constrain PINNs during the optimization.
For DEs with multiple solutions within the same domain, a follow-up question is to design an appropriate loss function to classify observations into different individual solutions. Merely decreasing the distance (which is measured by traditional loss functions such as mean squared errors, mean absolute errors, etc.) between the NN predictions and observations is insufficient, as the observations from different solutions of (1) are unclassified. It could lead to a poor generalization ability of the NNs. A more intelligent strategy would be to classify the observations while reducing the errors between the NN predictions and observations, namely,

\[L(\mathbf{\theta},\mathbf{\lambda})=\sum_{i=1}^{N_{v}}\min\|\hat{u}_{m}(\mathbf{x}_{i};\bm {\theta},\mathbf{\lambda})-u_{i}\|^{2}+\alpha\sum_{j=1}^{N_{v}}\sum_{m=1}^{M}\big{(} \mathcal{L}\hat{u}_{m}(\mathbf{x}_{j};\mathbf{\theta},\mathbf{\lambda})-f(\hat{u}_{m}(\mathbf{ x}_{j};\mathbf{\theta},\mathbf{\lambda}),\mathbf{x}_{j};\mathbf{\lambda})\big{)}^{2}. \tag{2}\]

Here \(\alpha\) is a Lagrange multiplier to balance the weight between the two loss terms. By assuming that the observations are from \(M\) or fewer smooth and continuous solutions of (1), \(\hat{u}_{m}(\cdot)\) is an approximation to one of the solutions. The first term in (2) quantifies the distance between \(\hat{u}_{m}(\cdot)\) and the observations, and determines, at a specific observation point \(\mathbf{x}_{i}\), which \(\hat{u}_{m}(\mathbf{x}_{i})\) is closest to the observations at \(\mathbf{x}_{i}\). In the optimization process, we optimize \(\mathbf{\theta}\) (parameters in NNs) and \(\mathbf{\lambda}\) (parameters in DEs). Such processes enable the NN to improve its ability to classify observations into different \(\hat{u}_{m}(\cdot)\). One concern is that incorrectly classifying an observation to a specific \(\hat{u}_{m}(\cdot)\) at the beginning can cause nonsmooth and discontinuous function approximation. It may lead to a buildup of errors, which can be hard to correct later. A possible solution is to introduce the second loss term in (2) which evaluates the residual of (1) using \(\hat{u}_{m}(\cdot)\) at collocation points. By properly weighting the two terms in (2), the observations can be correctly classified.

### Physics-informed neural networks with homotopy process

While a well-trained PINN may work for some simple inverse problems, recovering unknown parameters \(\mathbf{\lambda}\) in (1) with multiple solutions is still challenging. A bad initial guess of \(\mathbf{\lambda}\) can fail the network training. By leveraging the idea of homotopy continuation [24] to start with a simpler and related problem with known solutions, track the original problem to the deformed one, and finally identify unknown parameters in (1), we developed HomPINNs to solve the inverse problem of DEs with multiple solutions. The magnitude of \(\alpha\) in the loss function (2) is utilized to trace from the constraint for satisfying both observations and DEs to the constraint only for observations. During the homotopy process, the unknown parameters \(\mathbf{\lambda}\) are adjusted to gradually match the observations until the correct \(\mathbf{\lambda}\) will make the NN approximation fit the observations perfectly. Fig. 2 illustrates a general workflow for executing HomPINNs. HomPINNs feature a shared network structure such that the parameters obtained in the previous homotopy step become the initial guess of the next step. In step \(k\) (\(1\leq k\leq K\)), the loss function becomes

\[L_{k}(\mathbf{\theta},\mathbf{\lambda})=\frac{1}{N_{o}}\sum_{i=1}^{N_{o}}\min\|\hat{u}_ {m}(\mathbf{x}_{i};\mathbf{\theta},\mathbf{\lambda})-u_{i}\|^{2}+\frac{\alpha_{k}}{MN_{c}} \sum_{j=1}^{N_{v}}\sum_{m=1}^{M}\big{(}\mathcal{L}\hat{u}_{m}(\mathbf{x}_{j};\mathbf{ \theta},\mathbf{\lambda})-f(\hat{u}_{m}(\mathbf{x}_{j};\mathbf{\theta},\mathbf{\lambda}),\mathbf{ x}_{j};\mathbf{\lambda})\big{)}^{2}, \tag{3}\]

Figure 2: Schematic of HomPINNs: a NN at step \(k\) (\(k=1,2,\cdots,K\)) is constructed as Fig. 1, where the network parameters are initialized by a previous well-trained NN at step \(k-1\). Then the observations and collocation points are used to train the NN at each homotopy step. The homotopy process starts from step \(1\) and ends in step \(K\). The unknown parameter \(\mathbf{\lambda}^{\prime}\), the optimal parameters \(\mathbf{\theta}^{\prime}\), and the solutions of (1) will be obtained finally.
where \(\alpha_{k}\) is a homotopy tracking parameter in step \(k\). The value of \(\alpha_{k}\) decreases in a monotonic fashion, i.e., \(\alpha_{0}>\alpha_{1}>\cdots>\alpha_{k}>\cdots>\alpha_{K-1}>\alpha_{K}\). As \(\alpha_{k}\) decreases, the NN will track more of the given observations, and \(\lambda\) will also be adjusted to better fit the observations. Besides, we further claim that the tracking \(\alpha_{k}\) in this study is structured to decrease exponentially:

\[\alpha_{k}=\alpha_{0}r^{k-1}. \tag{4}\]

Here \(\alpha_{0}\) serves as an initial value to start the process, and \(r\) is the decay rate during the homotopy process. The training process of HomPINNs is summarized in Fig. 2 and algorithm 1. With the given observations, collocation points, and parameter initializations, HomPINNs consist of \(K\) homotopy steps inside which the optimization is performed. The present study uses the Adam optimizer [29] with a learning rate of \(10^{-3}\) and beta values of \(0.9\) and \(0.99\) in each homotopy step. Upon finishing \(K\) homotopy steps, the optimized parameters are produced as the final output.

```
1:Input Observations \(\{\mathbf{x}_{i},\mathbf{u}_{i}\}_{i=1}^{N_{k}}\).
2:Input Collocations \(\{\mathbf{x}_{j}\}_{j=1}^{N_{k}}\).
3:Initialize DE parameter \(\mathbf{\lambda}_{0}\).
4:Initialize PINN parameters \(\mathbf{\theta}_{0}\).
5:Initialize Number of outputs \(M\).
6:Initialize Homotopy tracking parameters \(\alpha_{0}\), \(r\)
7:for\(k=1,2,\cdots,\)\(K\)do
8: Update homotopy tracking parameter \(\alpha_{k}\) from (4).
9: Initialize \(\mathbf{\theta}=\mathbf{\theta}_{k-1}\), \(\mathbf{\lambda}=\mathbf{\lambda}_{k-1}\),
10: Optimize parameters \([\mathbf{\theta}_{k},\mathbf{\lambda}_{k}]=\operatorname*{argmin}_{\mathbf{\theta},\mathbf{ \lambda}}L_{k}(\mathbf{\theta},\mathbf{\lambda})\), where \(L_{k}(\mathbf{\theta},\mathbf{\lambda})\) is computed from (3).
11:endfor
12:Output DE parameters \(\mathbf{\lambda}_{K}\).
13:Output PINN parameters \(\mathbf{\theta}_{K}\).
```

**Algorithm 1** The training process of the proposed HomPINNs.

**Remark 2**.: _In certain situations, particularly when the number of solutions of (1) is large or the parameter space of \(\lambda\) is in high-dimension, HomPINNs may output values of \(\lambda^{*}\) that are not accurate enough. Our strategy is to regard the value of \(\lambda^{*}\) as the hyper-parameter in the homotopy physics-informed neural network for the forward problems (we will call it forward HomPINNs in the following) [30] and train it to obtain the optimized network parameters \(\mathbf{\theta}^{*}\). Once the training is completed, the optimized parameters \(\lambda^{*}\) and \(\mathbf{\theta}^{*}\) will initialize the NN for inverse problems (\(\mathbf{\lambda}_{0}=\lambda^{*}\) and \(\mathbf{\theta}_{0}=\mathbf{\theta}^{*}\)). Then we can run the proposed HomPINNs again to identify the correct \(\lambda\)._

## 3 Numerical examples

In this section, we conduct numerical experiments to verify the proposed methodology. We begin with one-dimensional problems to study the effect of hyper-parameters on the algorithm's performance and to provide general rules of thumb for implementing the method. Building on the above experiment, we proceed to solve a two-dimensional elliptic DE problem. All the tests are run on a desktop with the following specifications: Intel Core i9-10920X CPU, RTX 3090 GPU, and 64 GB DDR4 RAM memory. The number of collocation points is chosen such that the model performance is not influenced.

### 1D example with two solutions

The first one-dimensional DE considered is from [25] and can be expressed as:

\[\begin{cases}\dfrac{\partial^{2}u(x)}{\partial x^{2}}=-\lambda- \lambda u^{4},\ \ x\in(0,1)\\ \dfrac{\partial u(x)}{\partial x}\bigg{|}_{x=0}=u(x)|_{x=1}=0.\end{cases} \tag{5}\]
With \(\lambda=1.20\), (5) has two solutions, and 80 observations are randomly collected from these two solutions. As there are two solutions, \(\{u_{1}(x_{i}),u_{2}(x_{i})\}\), at a specific observation location \(x_{i}\), \(u_{i}\), the observation at \(x_{i}\), is randomly selected from \(\{u_{1}(x_{i}),u_{2}(x_{i})\}\). Also, all the observation locations \(\{x_{i}\}_{i=1}^{N_{o}}\) are randomly selected in the domain of interest. Such a strategy of generating observations is applied in all the tests in the present study.

We first discuss how to determine \(M\) in HomPINN. The NN used in this experiment is a fully connected network, with 1, 30, 30, 30, and \(M\) neurons in each layer, respectively. The network is initialized following the strategy in [31] to start the training in the first homotopy step. Table 1 summarizes the results with different \(M\). We observe sudden changes in the identified values of \(\lambda\), training loss, and testing loss when \(M\) changes from 1 to 2. This behavior suggests that the exact number of solutions should not be one. Additionally, when \(M\) increases from 2 to 3 or 4, the values of the identified \(\lambda\) from HomPINN do not change, which implies that there should be two solutions, i.e., \(M=2\). Moreover, the results in Table 1 also indicate that even when \(M\) is set to a value greater than the exact number of solutions, HomPINN can still accurately identify \(\lambda\). In light of the above results, the correct number of \(M\) can be estimated by trial and error.

Moreover, the predicted results made by the optimized HomPINNs are illustrated in Fig. 3 with \(M=2\), along with the observations from the two different solutions of (5).

The change of training loss during the homotopy processes for different \(M\) is depicted in Fig. 4. When using \(M=1\), the training loss maintains a value of \(2\cdot 10^{-2}\) consistently across the homotopy process. This behavior also illustrates that \(M=1\) is not a wise choice for this example. On the other hand, when \(M\) is set as 2, 3, or 4, the training losses of HomPINNs decrease as the homotopy steps proceed and converge at nearly a sub-linear rate.

Our next task is to examine the effects of the number of observations (\(N_{o}\)). Based on the prior results, we use \(M=2\) and decrease \(N_{o}\) from 80 to 60, 40, and 20, and the results are shown in Table 2. A smaller amount of observations produce a larger error of the identified \(\lambda\), training loss, and testing loss, yet they remain on a similar level. These results indicate that the proposed method can still be effective even when the amount of observations is limited. This property is favorable when computing resources are constrained. Boosting the number of observations helps to achieve a higher level of accuracy from HomPINNs.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(M\) & PARAM value & PARAM error & Train loss & Test loss \\ \hline
1 & 1.3258 & 1.26E-1 & 2.04E-2 & 1.18E-1 \\ \hline
2 & 1.2000 & 3.90E-6 & 1.89E-9 & 2.13E-6 \\ \hline
3 & 1.2000 & 4.60E-6 & 4.63E-9 & 2.07E-6 \\ \hline
4 & 1.2000 & 1.50E-6 & 2.98E-9 & 1.98E-6 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Values of the identified \(\lambda\) from HomPINN, errors of the identified \(\lambda\) from HomPINN, training losses, and testing losses, using different input \(M\) to HomPINN. The correct value of \(\lambda\) is 1.20.

Figure 3: Observations sampled from the two solutions of (5) and predictions made from HomPINNs. Here the x-axis is \(x\), and the y-axis is \(u\). The black circles are the sampled observations. The red solid line with squares and the blue solid line with stars are the predicted \(u\) values from HomPINNs. The correct parameter value \(\lambda=1.20\), and with the correct estimate of \(M=2\), the identified parameter given by HomPINNs \(\lambda\) is 1.200.
### 1D example with multiple solutions

We delve into a more intricate example having 7 solutions [26]:

\[\left\{\begin{aligned} &\frac{\partial^{2}u(x)}{\partial x^{2}}=u^{4}-pu^{2}, \ \ x\in(0\,1),\\ &\left.\frac{\partial u(x)}{\partial x}\right|_{x=0}=\left.u(x) \right|_{x=1}=0.\end{aligned}\right. \tag{6}\]

The proposed HomPINN attempts to specify the value of \(p\) given 210 observations (\(N_{o}=210\)). The observations \(\{x_{i},u_{i}\}_{i=1}^{210}\) are randomly selected from 7 solutions of (6) with the exact value of \(p=18.00\). However, we discover that using HomPINN once is unable to give an accurate enough output of \(p\) because the current example is more complex. To alleviate this inaccuracy, we implement the strategy in Remark 2. Our test shows that two processes of HomPINN are sufficient in this example.

Table 3 shows the outputs of the two processes of HomPINN with different \(M\). We observe the same jump of errors when we move from \(M=6\) to \(M=7\), implying that (6) should have 7 solutions. More insights can be found in Fig.

5. Fig. 5(a) displays predictions from one complete process of HomPINN, and the corresponding output of \(p\) is 17.3516. Then, we use \(p=17.3516\) as an input to train the forward HomPINNs, whose output of \(u\) are depicted in Fig. 5(b). Then, we use \(p=17.3516\) and the NN parameters of the trained forward HomPINNs to start the second

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(N_{o}\) & PARAM value & PARAM error & Train loss & Test loss \\ \hline
20 & 1.2000 & 3.94E-5 & 6.35E-9 & 9.49E-6 \\ \hline
40 & 1.2000 & 2.99E-5 & 5.81E-9 & 8.60E-6 \\ \hline
60 & 1.2000 & 1.13E-5 & 2.26E-9 & 7.66E-6 \\ \hline
80 & 1.2000 & 0.39E-5 & 1.89E-9 & 2.13E-6 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Values of the identified parameter from HomPINNs, errors of the identified parameter from HomPINN, training loss, and testing loss when different numbers of observations are used in HomPINNs. The correct value of \(\lambda\) is 1.20.

Figure 4: Training loss during the homotopy steps with \(M=1,2,3,4\). Here the x-axis indicates the number of epochs during the homotopy process, where \(\alpha_{k}\) (\(k=1,2,\cdots,11\)) is the index of each homotopy step. And the y-axis is the training loss in the logarithmic scales.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(M\) & 1st PARAM & 2rd PARAM & PARAM error & Train loss \\ \hline
5 & 17.3050 & 17.2868 & 7.13E-1 & 1.51E-1 \\ \hline
6 & 17.4538 & 17.6210 & 3.79E-1 & 3.65E-2 \\ \hline
7 & 17.3516 & 17.9982 & 1.84E-3 & 2.15E-5 \\ \hline
8 & 17.2918 & 17.9987 & 1.26E-3 & 1.79E-5 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Values of the identified parameter in the first and second homotopy process of HomPINN, and training losses. The correct value of \(p\) is 18.00.
homotopy process of HomPINN, and its results are shown in Fig. 5(c). As shown in Fig. 5(b), even though ruining HomPINN one homotopy process only outputs a value of \(p\) that are not very accurate, HomPINN provides a better prediction of \(u\) with the given value of \(p\) from HomPINN. In other words, the NN parameters (\(\mathbf{\theta}\)) in the trained forward HomPINN are closer to their global optimum. As a result, using these NN parameters to restart the second homotopy process of HomPINNs helps to find a more accurate value of \(p\). Once \(p\) is specified from HomPINNs, we use another forward HomPINNs to obtain all the solutions of (6), as shown in Fig. 5(c).

In the next test, the observations \(\{x_{i},u_{i}\}_{i=1}^{180}\) are randomly selected from three of the seven solutions (see black circles in Fig. 6(a)). Under this configuration, we found that the trained HomPINNs after one homotopy process can recover \(p\) well. We subsequently apply forward HomPINNs with the recovered value of \(p\) to predict \(u\), as shown in Fig. 6(b). The value and error of \(p\) from the forward HomPINNs and training losses are also shown in Table 4.

\begin{table}
\begin{tabular}{c c c c} \hline \hline \(M\) & PARAM value & PARAM error & Training loss \\ \hline
1 & 6.1780 & 7.13E-1 & 8.86E0 \\ \hline
2 & 15.9029 & 3.79E-1 & 1.46E0 \\ \hline
3 & 17.9983 & 1.84E-3 & 5.77E-5 \\ \hline
4 & 17.9989 & 1.26E-3 & 1.62E-5 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Values of the identified parameter (\(p\)) from HomPINN, errors of the identified parameters from HomPINN, and final losses. The observation are sampled from three of the seven solutions (only partial information is given). The correct value of \(p\) is 18.

Figure 5: Results from two processes of HomPINN. Here the x-axis is \(x\), and the y-axis is \(u\). The black circles are observations sampled from the 7 solutions, and seven solid lines represent the predicted solutions of (6). (a): the results after the first homotopy process of HomPINNs whose output of \(p\) is 17.3516. (b): the results from the forward HomPINNs given parameter \(p=17.3516\). (c): the second homotopy process of HomPINNs whose output of \(p\) is 17.9982.
It is encouraging to note that the proposed HomPINNs are able to specify the unknown parameters accurately even when the observations do not contain all the solutions of a DE. In fact, observing fewer solutions can simplify HomPINN training as overlapping or intersecting among different solutions can make the training more challenging, in particular when the DEs is defined in a high-dimension domain. Once the unknown parameter is specified (\(p=17.9983\)), the forward HomPINN are utilized to find all the solutions of (6) (see Fig. 6(b)). The training loss in HomPINN to specify \(p\) is shown in Fig. 7.

### Gray-Scott simulations

The Gray-Scott model mathematically represents the reaction-diffusion process of two chemical species interacting with each other and their environment. Its origins date back to the 1980s when Gray and Scott [32] first introduced it. Since then, it has been used to simulate a wide range of phenomena, such as pattern formation [28], chemical waves [33], and chemical turbulence [34]. The model describes the changes in concentration of the two species, \(A(x,y,t)\) and \(S(x,y,t)\). The present study considers the steady state of the Gray-Scott model with the zero-flux boundary condition:

\[\begin{cases}D_{A}\Delta A+SA^{2}-(\mu+\rho)A=0,\\ D_{S}\Delta S-SA^{2}+\rho(1-S)=0,\\ \mathbf{n}\cdot\nabla A|_{\partial\Omega}=\mathbf{n}\cdot\nabla S|_{\partial \Omega}=0.\end{cases} \tag{7}\]

Here \(\mathbf{n}\) denotes the outward normal at the domain boundary, the two species have diffusion coefficients represented by \(D_{A}\) and \(D_{S}\), respectively. Parameter \(\rho\) denotes the constant supply of species \(S\) into the system, and \(\mu\) represents the rate of decay of species \(A\).

Figure 6: The results after one homotopy process of HomPINNs and the predictions made by the forward HomPINNs after the parameter is recovered. Here x-axis is \(x\), and the y-axis is \(u\). The black circles are observations sampled from three of the seven solutions, and solid lines represent the predicted solutions of (6). Left: the trained HomPINNs after one homotopy process. The recovered parameter: \(p=17.9983\). Right: all the potential solutions of (6) predicted by the forward HomPINNs given \(p=17.9983\) from the proposed HomPINN.

Figure 7: The training loss in HomPINN of the example in 3.2. Here \(\alpha_{k}\) in the x-axis is the indices of the homotopy step, and the y-axis is the training loss (logarithmic scale). The observations for this training are sampled from three of the seven solutions, and the black solid line is the training loss in HomPINN when \(M=3\). The training loss decreases within each homotopy step, and the output value of \(p\) is17.9983.
To generate observations, we use \(D_{A}=2.5\times 10^{-4}\), \(D_{S}=5\times 10^{-4}\), \(\rho=4.0\times 10^{-2}\), and \(\mu=6.5\times 10^{-2}\). The observations do not contain all the solutions of (7) with the parameters chosen. Our goal is to utilize HomPINNs to identify the values of the parameters and retrieve all the solutions of (7) given the parameters from HomPINNs. To tackle this complicated task, we first employed trial and error to determine the number of observations (\(N_{o}\)) and NN structures (the number of neurons in each layer and the number of layers) that are most effective for HomPINNs.

We first randomly generate a total of 10,000 observations \(\{\mathbf{x}_{i},\mathbf{u}_{i}\}_{i=1}^{10,000}\) in the four known solutions (see black circles in Fig. 8). Then We randomly select 1,000 and 5,000 samples from \(\{\mathbf{x}_{i},\mathbf{u}_{i}\}_{i=1}^{10,000}\) as the observations input to HomPINNs.

We also investigate the network structure by altering the number of neurons in each layer (width) and the number of layers. Table 5 summarizes the values of the identified parameters from HomPINNs and training loss under different conditions of the number of observations, network width, and network layers. We discover that HomPINNs can fail for several reasons, such as limited observations, inadequate network width, or insufficient network layers. HomPINNs become more effective when the number of observations exceeds 5,000, the network width is larger than 128, and the number of layers is greater than 6. When the number of given observations is too small, HomPINNs may not be able to distinguish different patterns (solutions) in the observations accurately. This can lead to overfitting of the model to the observations, which negatively impacts its generalization capacity. The number of neurons and layers in a NN determines the width and depth of the network. If one of them is too small, HomPINNs may not be able to capture the complexity of the observations accurately and unable to generalize the patterns of solutions well. The performance of HomPINNs can be enhanced by increasing the number of observations, network width, and the number of network layers to capture complex patterns in observations and provide more accurate predictions.

Based on the above analysis, we opted to select 5,000 observations randomly from the four known solutions (see Fig. 9) and configure the NN with 6 hidden layers, each of which has 128 neurons.

After the first process of HomPINN, the values of the identified parameters are \(D_{A}=2.020\times 10^{-4}\), \(D_{S}=4.849\times 10^{-4}\), \(\rho=3.678\times 10^{-2}\), and \(\mu=6.882\times 10^{-2}\) (see Fig. 10 and Table 5), which indicates that HomPINNs recover the unknown parameters to some extent. Nevertheless, it is unable to discern different solutions associated with the observations. From Fig. 10, the results produced by HomPINN may mix observations from multiple solutions, yet the given solutions are still satisfying the governing equation (7).

We implement the strategy in Remark 2 to improve the accuracy. Using the parameters specified by the first homotopy process HomPINN as inputs, we perform the forward HomPINNs following the process outlined in [30]. Fig. 11 depicts the forward HomPINN predictions to \(A(\cdot)\) and \(S(\cdot)\), which are closer to the exact solution illustrated in Fig. 8, although the input parameters from the forward HomPINNs are not very accurate. Subsequently, the

Figure 8: The four known solutions of the Gray–Scott model with the chosen parameters to obtain observation input to HomPINN. These contour plots show the distribution of the component \(A(x,y)\) (top four figures) and \(S(x,y)\) (bottom four figures) on a domain of \([0,1]\times[0,1]\), with the horizontal and vertical axes representing x and y, respectively. The parameters used here: \(D_{A}=2.5\times 10^{-4}\), \(D_{S}=5\times 10^{-4}\), \(\rho=4.0\times 10^{-2}\), and \(\mu=6.5\times 10^{-2}\).
Figure 10: Predictions of the Gray–Scott model after the first homotopy process of HomPINNs. The contour plots are for component \(A(x,y)\). The horizontal and vertical axes represent x and y, respectively. The domain is \([0,1]\times[0,1]\). Learned parameters from HomPINN: \(D_{A}=2.020\times 10^{-4}\), \(D_{S}=4.849\times 10^{-4}\), \(\rho=3.678\times 10^{-2}\), and \(\mu=6.882\times 10^{-2}\).

Figure 9: Sampled observations of the Gray–Scott model in 2D. The contour plots are for the component \(A(x,y)\). The horizontal and vertical axes represent x and y, respectively, and the domain is \([0,1]\times[0,1]\). The black circles are sampled observations.
NN parameters in the trained forward HomPINNs along with the DE parameters specified from the first homotopy process of HomPINNs are used to initialize the second homotopy process of HomPINNs. With DE parameters specified in the second homotopy process of HomPINNs being closer to the correct ones, our predictions of \(A(\cdot)\) and \(S(\cdot)\) are improved in the second homotopy process of HomPINNs (see Fig. 12). Tables 6 and 7 provide a summary of the training outcomes for the two homotopy processes of HomPINNs, including the DE parameters obtained in each homotopy step, the parameter errors, and the training losses. The results demonstrate that the DE parameters can be correctly determined by HomPINNs. Once we identify the correct DE parameters from HomPINN, these parameters are input to the forward HomPINNs to discover all the solutions of the governing equation (7), as depicted in Fig. 13.
Figure 13: All the solutions of (7) from the forward HomPINNs given the parameters after the second homotopy process of HomPINN. The contour plots are for component \(A(x,y)\). The horizontal and vertical axes represent x and y, respectively. The parameters used here for the HomPINN are: \(D_{A}=2.409\times 10^{-4}\), \(D_{S}=5.064\times 10^{-4}\), \(\rho=4.076\times 10^{-2}\), and \(\mu=6.432\times 10^{-2}\).
## 4 Conclusion and Discussion

To tackle the complex behavior arising from non-uniqueness, symmetry, and bifurcations in the solution space, we proposed a novel HomPINN method for solving inverse problems in nonlinear DEs with multiple solutions. The proposed method utilizes a combination of homotopy continuation and NNs in the training process to estimate unknown parameters in the DEs and is able to produce precise predictions. We evaluated the effectiveness of our proposed method by conducting experiments on different tasks. For DEs in one dimension, we explore the effect of various factors on the model's performance, such as estimates of the number of solutions, the number of observations, and observations that do not include all the DE solutions. Next, the two-dimensional Gray-Scott problem is investigated to determine the influence of the network structure. Through the tests conducted, we determine the optimal number of observations and network structure. With the determined hyper-parameters, HomPINN successfully specifies the values of the four parameters in the Gray-Scott problem, although the input observations only have 4 of the solutions. All the solutions to the Grey-Scott problem are finally recovered with HomPINN, through the use of the parameters from the proposed HomPINN.

Our experimental results indicate the following advantages of the proposed model: 1. HomPINN offers a flexible and adaptable solution for inverse problems with large-scale datasets; 2. HomPINN is able to tackle both one and two-dimensional problems; 3. the optimal number of observations and network structure has been identified in our tests; 4. HomPINN is able to solve inverse problems even if we do not know all the solutions of a DE; 5. It is possible to discover other solutions of a DE that are not available in observations.

In the future, we will focus on creating techniques to automatically identify the number of solutions or investigating more efficient homotopy continuation methods. In addition, we aspire to evaluate our method on a wider range of problems, including those with more intricate geometries and boundary conditions, to thoroughly assess its performance and versatility.

## 5 Acknowledgements

G.L. and H.Z. gratefully acknowledge the support of the National Science Foundation (DMS-1555072, DMS-2053746, and DMS-2134209), Brookhaven National Laboratory Subcontract 382247, and U.S. Department of Energy (DOE) Office of Science Advanced Scientific Computing Research program DE-SC0021142 and DE-SC0023161.

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline \multirow{2}{*}{\(k\)} & \multicolumn{6}{c}{PARAM value} & \multicolumn{6}{c}{PARAM error} & \multirow{2}{*}{Train Loss} \\ \cline{2-3} \cline{7-10}  & \(D_{A}\) & \(D_{S}\) & \(\rho\) & \(\mu\) & \(D_{A}\) & \(D_{S}\) & \(\rho\) & \(\mu\) & \\ \hline  & 1.000E-5 & 2.000E-5 & 1.000E-3 & 1.000E-3 & 2.400E-4 & 4.800E-4 & 3.900E-2 & 6.400E-2 & 6.157E-3 \\
1 & 1.620E-4 & 4.709E-4 & 3.820E-2 & 8.931E-2 & 8.800E-5 & 2.911E-5 & 1.800E-3 & 2.431E-2 & 5.051E-4 \\
2 & 1.875E-4 & 4.816E-4 & 3.710E-2 & 7.155E-2 & 6.252E-5 & 1.840E-5 & 2.904E-3 & 6.552E-3 & 2.296E-5 \\
3 & 1.948E-4 & 4.821E-4 & 3.688E-2 & 6.954E-2 & 5.521E-5 & 1.795E-5 & 3.121E-3 & 4.537E-3 & 2.643E-5 \\
4 & 1.970E-4 & 4.847E-4 & 3.688E-2 & 6.909E-2 & 5.297E-5 & 1.526E-5 & 3.120E-3 & 4.088E-3 & 7.933E-6 \\
5 & 1.989E-4 & 4.836E-4 & 3.678E-2 & 6.897E-2 & 5.110E-5 & 1.639E-5 & 3.221E-3 & 3.965E-3 & 4.712E-6 \\
6 & 1.999E-4 & 4.840E-4 & 3.675E-2 & 6.884E-2 & 5.010E-5 & 1.602E-5 & 3.248E-3 & 3.842E-3 & 2.091E-6 \\
7 & 2.008E-4 & 4.841E-4 & 3.677E-2 & 6.884E-2 & 4.917E-5 & 1.586E-5 & 3.234E-3 & 3.835E-3 & 2.057E-6 \\
8 & 2.013E-4 & 4.842E-4 & 3.676E-2 & 6.885E-2 & 4.872E-5 & 1.579E-5 & 3.236E-3 & 3.848E-3 & 1.307E-6 \\
9 & 2.016E-4 & 4.848E-4 & 3.677E-2 & 6.883E-2 & 4.837E-5 & 1.522E-5 & 3.230E-3 & 3.829E-3 & 7.243E-7 \\
10 & 2.020E-4 & 4.849E-4 & 3.678E-2 & 6.882E-2 & 4.798E-5 & 1.507E-5 & 3.224E-3 & 3.816E-3 & 6.376E-7 \\ \hline \hline \end{tabular}
\end{table}
Table 6: Values and errors of the identified parameters from the first homotopy process of HomPINNs, and training losses among the first homotopy process of HomPINNs. The correct parameters are \(D_{A}=2.5\times 10^{-4}\), \(D_{S}=5\times 10^{-4}\), \(\rho=0.04\), and \(\mu=0.065\).

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline \multirow{2}{*}{\(k\)} & \multicolumn{6}{c}{PARAM value} & \multicolumn{6}{c}{PARAM error} & \multirow{2}{*}{Train Loss} \\ \cline{2-3} \cline{7-10}  & \(D_{A}\) & \(D_{S}\) & \(\rho\) & \(\mu\) & \(D_{A}\) & \(D_{S}\) & \(\rho\) & \(\mu\) & \\ \hline  & 2.020E-4 & 4.849E-4 & 3.678E-2 & 6.882E-2 & 4.798E-5 & 1.507E-5 & 3.224E-3 & 3.816E-3 & 1.339E-5 \\
1 & 2.384E-4 & 5.057E-4 & 4.072E-2 & 6.424E-2 & 1.159E-5 & 5.699E-6 & 7.212E-4 & 7.581E-4 & 9.054E-6 \\
2 & 2.397E-4 & 5.059E-4 & 4.073E-2 & 6.428E-2 & 1.028E-5 & 5.889E-6 & 7.286E-4 & 7.199E-4 & 5.097E-6 \\
3 & 2.404E-4 & 5.061E-4 & 4.074E-2 & 6.429E-2 & 9.638E-6 & 6.098E-6 & 7.407E-4 & 7.072E-4 & 3.194E-6 \\
4 & 2.406E-4 & 5.063E-4 & 4.075E-2 & 6.430E-2 & 9.430E-6 & 6.298E-6 & 7.539E-4 & 7.009E-4 & 1.891E-6 \\
5 & 2.407E-4 & 5.064E-4 & 4.076E-2 & 6.430E-2 & 9.268E-6 & 6.360E-6 & 7.585E-4 & 7.009E-4 & 1.129E-6 \\
6 & 2.408E-4 & 5.063E-4 & 4.076E-2 & 6.430E-2 & 9.222E-6 & 6.339E-6 & 7.610E-4 & 7.009E-4 & 7.102E-7 \\
7 & 2.409E-4 & 5.064E-4 & 4.076E-2 & 6.431E-2 & 9.129E-6 & 6.360E-6 & 7.589E-4 & 6.945E-4 & 4.132E-7 \\
8 & 2.409E-4 & 5.064E-4 & 4.076E-2 & 6.431E-2 & 9.106E-6 & 6.396E-6 & 7.610E-4 & 6.945E-4 & 2.435E-7 \\
9 & 2.409E-4 & 5.063E-4 & 4.076E-2 & 6.432E-2 & 9.083E-6 & 6.350E-6 & 7.585E-4 & 6.818E-4 & 1.642E-7 \\
10 & 2.409E-4 & 5.064E-4 & 4.076E-2 & 6.432E-2 & 9.106E-6 & 6.360E-6 & 7.585E-4 & 6.818E-4 & 1.249E-7 \\ \hline \hline \end{tabular}
\end{table}
Table 7: Values and errors of the identified parameters from the second homotopy process of HomPINN, and training losses among the second homotopy process of HomPINNs. The correct parameter are \(D_{A}=2.5\times 10^{-4}\), \(D_{S}=5\times 10^{-4}\), \(\rho=0.04\), and \(\mu=0.065\).