Data-driven localized waves and parameter discovery in the massive
  Thirring model via extended physics-informed neural networks with interface
  zones
In this paper, we study data-driven localized wave solutions and parameter
discovery in the massive Thirring (MT) model via the deep learning in the
framework of physics-informed neural networks (PINNs) algorithm. Abundant
data-driven solutions including soliton of bright/dark type, breather and rogue
wave are simulated accurately and analyzed contrastively with relative and
absolute errors. For higher-order localized wave solutions, we employ the
extended PINNs (XPINNs) with domain decomposition to capture the complete
pictures of dynamic behaviors such as soliton collisions, breather oscillations
and rogue-wave superposition. In particular, we modify the interface line in
domain decomposition of XPINNs into a small interface zone and introduce the
pseudo initial, residual and gradient conditions as interface conditions linked
adjacently with individual neural networks. Then this modified approach is
applied successfully to various solutions ranging from bright-bright soliton,
dark-dark soliton, dark-antidark soliton, general breather, Kuznetsov-Ma
breather and second-order rogue wave. Experimental results show that this
improved version of XPINNs reduce the complexity of computation with faster
convergence rate and keep the quality of learned solutions with smoother
stitching performance as well. For the inverse problems, the unknown
coefficient parameters of linear and nonlinear terms in the MT model are
identified accurately with and without noise by using the classical PINNs
algorithm.
Junchao Chen, Jin Song, Zijian Zhou, Zhenya Yan
Data-driven localized waves and parameter discovery in the massive Thirring model via extended physics-informed neural networks with interface zones

###### Abstract

In this paper, we study data-driven localized wave solutions and parameter discovery in the massive Thirring (MT) model via the deep learning in the framework of physics-informed neural networks (PINNs) algorithm. Abundant data-driven solutions including soliton of bright/dark type, breather and rogue wave are simulated accurately and analyzed contrastively with relative and absolute errors. For higher-order localized wave solutions, we employ the extended PINNs (XPINNs) with domain decomposition to capture the complete pictures of dynamic behaviors such as soliton collisions, breather oscillations and rogue-wave superposition. In particular, we modify the interface line in domain decomposition of XPINNs into a small interface zone and introduce the pseudo initial, residual and gradient conditions as interface conditions linked adjacently with individual neural networks. Then this modified approach is applied successfully to various solutions ranging from bright-bright soliton, dark-dark soliton, dark-antidark soliton, general breather, Kuznetsov-Ma breather and second-order rogue wave. Experimental results show that this improved version of XPINNs reduce the complexity of computation with faster convergence rate and keep the quality of learned solutions with smoother stitching performance as well. For the inverse problems, the unknown coefficient parameters of linear and nonlinear terms in the MT model are identified accurately with and without noise by using the classical PINNs algorithm.

keywords: Deep learning, XPINNs algorithm with interface zones, Massive Thirring model, Data-driven localized waves, Parameter discovery +
Footnote †: journal: ****

## 1 Introduction

The application of deep learning in the form of neural networks (NNs) to solve partial differential equations (PDEs) has recently been explored in the field of scientific machine learning due to the universal approximation and powerful performance of NNs [1; 2; 3]. In particular, by adding a residual term from the given PDE to the loss function with the aid of automatic differentiation, the physics-informed neural networks (PINNs) approach has been proposed to accurately solve forward problems of predicting the solutions and inverse problems of identifying the model parameters from the measured data [4]. Unlike the traditional numerical techniques such as the finite element method and finite difference method, the PINNs method is a mesh-free algorithm and requires the relatively small amounts of data [5; 6]. Based on these advantages, the PINNs has also been applied extensively to different types of PDEs, including integro-differential equations [6], fractional PDEs [7], and stochastic PDEs [8; 9]. However, one of the main disadvantages of PINNs is the high computational cost associated with training NNs, which can lead to the poor performance and low accuracy. Especially, for modelling problems involving the long-term integration of PDEs, the large amount of data will lead to a
rapid increase in the cost of training. To reduce the computational cost, several improved methods based on PINNs have been introduced to accelerate the convergence of models without loss of performance. The conservative PINNs (cPINNs) method has been established on discrete subdomains obtained after dividing the computational domain for nonlinear conservation laws, where the subdomains are stitched together on the basis of flux continuity across subdomain interfaces [10]. The extended PINNs (XPINNs) approach with the generalized domain decomposition has been proposed to any type of PDEs, where the XPINNs formulation offers highly irregular, convex/non-convex space-time domain decomposition thereby reducing the training cost effectively [11]. The _hp_-variational PINNs (_hp_-VPINNs) framework based on the nonlinear approximation of NNs and _hp_-refinement via domain decomposition has has been formulated to solved PDEs using the variational principle [12].

Based on cPINNs and XPINNs, a parallel PINNs (PPINNs) via domain decomposition has recently been developed to effectively address the multi-scale and multi-physics problems [13; 14]. The augmented PINNs (APINNs) has been proposed with soft and trainable domain decomposition and flexible parameter [15]. The PINNs/XPINNs algorithm has been employed to approximate supersonic compressible flows [16] and high-speed aerodynamic flows [17], and to quantify the microstructure of polycrystalline nickel [18]. The highly ill-posed inverse water wave problems governed by Serre-Green-Naghdi equations have been solved by means of PINNs/XPINNs to choose optimal location of sensors [19]. The unified scalable framework for PINNs with temporal domain decomposition has been introduced for time-dependent PDEs[20]. The multi-stage training algorithm consisting of individual PINNs has recently been presented via domain decomposition only in the direction of time [21; 22]. Indeed, the theoretical analyses of PINNs/XPINNs including convergence and generalization properties have been proved for linear and nonlinear PDEs through the rigorous error estimates [23; 24; 25], and the conditions under which XPINNs improve generalization have been discussed by proving certain generalization bounds [26]. In addition, the activation function is one of the important hyperparameters in NNs, which usually acts on affine transformation. The activation function possesses the ability to introduce nonlinearity in NNs in order to capture nonlinear features, and popular activation functions include rectified linear units, maxout units, logistic sigmoid, hyperbolic tangent function and so on [27]. In order to improve the performance of NNs and the convergence speed of objective function, global and local adaptive activation functions have been proposed by introducing scalable parameters into the activation function and adding a slope recovery term into the loss function [28; 29]. Furthermore, a generalized Kronecker NNs for any adaptive activation functions has been established by Jagtap et al [30]. The comprehensive survey, the applications-based taxonomy of activation functions and the state-of-the-art adaptive activation functions as well as the systematic comparison of various fixed and adaptive activation functions have been performed in [27]. These adaptive activation functions have also been used extensively in XPINNs (for example, on the inverse supersonic flows [16]).

As a special class of nonlinear PDEs, integrable systems possess remarkable mathematical structures such as infinitely many conservation laws, Lax pairs, abundant symmetry, various exact solutions and so on. For such type of nonlinear system, the PINNs algorithm could benefit greatly from the underlying integrable properties and achieve a more accurate numerical prediction. Particularly, a rich variety of solutions with local features such as soliton, breather, rogue wave and periodic wave can be derived exactly by using the classical methods for studying integrable equations, including the inverse scattering transformation, the Darboux transformation, the Hirota bilinear method and the Backlund transformation. These solutions potentially provide us with a large number of training samples for numerical experiments in the framework of PINNs. For example, the classical and improved PINNs method based on explicit solutions has been used to study data-driven nonlinear wave and parameter discovery in many integrable equations, including KdV equation [4; 31], mKdV equation [31], nonlinear Schrodinger (NLS) equation [4; 32; 33; 34], derivative NLS equation [35; 36], Manakov system [37; 21; 38], Yajima-Oikawa system [39], Camassa-Holm equation [40] and etc. It is worth noting that a modified PINNs approach developed by introducing conserved quantities from integrable theory into NNs, has been shown to achieve higher prediction accuracy [41; 42]. Moreover, the PINNs scheme has been applied to discover Backlund transformation and Miura transform [43], which are two typical transformations in the field of integrable systems. The PINNs schemes based on Miura transformations have been established
as an implementation method of unsupervised learning for predicting new local wave solutions [44]. In addition, the PINN deep learning has been employed to solve forward and inverse problems of the defocusing NLS equation with a rational potential [45], the logarithmic NLS with \(\mathcal{PT}\)-symmetric harmonic potential [46], the \(\mathcal{PT}\)-symmetric saturable NLS equation [47] and other nearly integrable models [48, 49, 50, 51, 22]. More recently, the deep learning method with Fourier neural operator was used to study solitons of both integrable nonlinear wave equations [52] and integrable fractional nonlinear wave equations [53].

The massive Thirring (MT) model takes the following form in laboratory coordinates

\[\begin{split}\mathrm{i}u_{x}+v+u|v|^{2}&=0,\\ \mathrm{i}v_{t}+u+v|u|^{2}&=0.\end{split} \tag{1}\]

This two-component nonlinear wave evolution model is known to be completely integrable in the sense of possessing infinitely many conserved quantities [54] and by means of the inverse scattering transform method [55, 56]. It reduces to the massless Thirring model [57] if the linear terms \(u\) and \(v\) are ignored. In field theory, the MT model (1) corresponds to an exactly solvable example for one-dimensional nonlinear Dirac system arising as a relativistic extension of the NLS equation [57]. In nonlinear optics, both components in Eq.(1) represent envelopes of the forward and backward waves respectively [58]. Indeed, the MT model (1) is a particular case of the coupled mode equations with self-phase modulation [58]. This coupled mode system has been extensively used to describe pulse propagation in periodic or Bragg nonlinear optical media [59, 60, 61, 62, 63], deep water waves for a periodic bottom in ocean [64, 65] and superpositions of two hyperfine Zeeman sublevels in atomic Bose-Einstein condensates [66]. The Coleman correspondence between the sine-Gordon equation and the MT model has been found and used to generate solutions in [67]. Various solutions of the MT model like soliton, rogue wave and algebro-geometric solutions as well as other integrable properties have been widely studied by using many different techniques such as the Darboux transformation, Backlund transformation, Hirota bilinear method and dressing method [68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83]. With the aid of the bilinear Kadomtsev-Petviashvili (KP) hierarchy reduction approach, one of the authors has systematacially constructed tau-function solutions of the MT model (1) including bright, dark soliton and breather as well as rogue wave in compact form of determinants [84, 85].

In this paper, we aim to explore data-driven localized waves such as soliton, breather and rogue wave, as well as parameter discovery in the MT model (1) via the deep learning algorithm. In forward problems of predicting nonlinear wave, we will employ the traditional PINNs framework to learn the simple solutions including one-soliton of bright/dark type and first-order rogue wave. For the solutions of two-soliton, one-breather and second-order rogue wave with complicated structures, it is necessary to enlarge the computational domain to capture the complete pictures of dynamic behaviors such as soliton collisions, breather oscillations and rogue wave superpositions. These requirements force us to utilize the XPINNs approach with domain decomposition for the data-driven experiments correspondingly. More specifically, the computational domain is divided into a number of subdomains along the time variable \(t\) when we treat the MT model as a time evolution problem, hence the XPINNs architecture is also regarded as a multi-stage training algorithm in this situation. In particular, we slightly modify the interface line in the domain decomposition of XPINNs into a small interface zone shared by two adjacent subdomains. In the training process, the pseudo initial, residual and gradient points are randomly selected in the small interface zones and are subjected to three types of interface conditions for better stitching effect and faster convergence. In inverse problems of discovering parameter, the classical PINNs algorithm is applied to identify the coefficients of linear and nonlinear terms in the MT model (1) in the absence and presence of noise.

The remainder of the paper is organized as follows. In Section 2, we firstly introduce the classical PINNs method and develop the modified XPINNs approach with small interface zones, in which the domain decomposition with sampling points and the schematic of individual PINN with interface conditions are illustrated in detail. Then dynamic behaviors of data-driven localized waves in the MT model (1) ranging widely from one-soliton of bright/dark type, bright-bright soliton, dark-dark soliton, dark-antidark soliton, general breather, Kuznetsov-Ma breather and rogue wave of first- and second-order, are presented via the traditional PINNs
and modified XPINNs methods. In Section 3, based on three types of localized wave solution, we discuss data-driven parameter discovery in the MT model (1) with and without noise via the classical PINNs algorithm. Conclusions and discussions are given in the last section.

## 2 Data-driven localized waves of the MT model

In this section, we will use PINNs and XPINNs to solve forward problems of the MT model (1). More precisely, we will focus on the the MT model (1) with initial and Dirichlet boundary conditions as follows:

\[\left\{\begin{array}{l}\mathrm{i}u_{x}+v+u|v|^{2}=0,\\ \mathrm{i}v_{t}+u+v|u|^{2}=0,\\ u(x,T_{0})=u_{0}(x),\;\;v(x,T_{0})=v_{0}(x),\;\;x\in[L_{0},L_{1}],\\ u(L_{0},t)=u_{lb}(t),\;\;v(L_{0},t)=v_{lb}(t),\\ u(L_{1},t)=u_{ub}(t),\;\;v(L_{1},t)=v_{ub}(t),\end{array}\right.t\in[T_{0},T_ {1}]. \tag{2}\]

By taking the complex-valued solutions as \(u=p+\mathrm{i}q\) and \(v=r+\mathrm{i}\) is with the real-valued functions \((p,q)\) and \((r,s)\) being the real and imaginary parts of \(u\) and \(v\) respectively, the MT model (2) is decomposed into four real equations. According to the idea of PINNs, we need to construct a complex-valued neural network to approximate the solutions \(u\) and \(v\). The complex-valued PINNs for the MT model (2) can be defined as

\[\begin{array}{l}f_{u}:=\mathrm{i}\hat{u}_{x}+\hat{v}+\hat{u}|\vartheta|^{2},\\ f_{v}:=\mathrm{i}\hat{v}_{t}+\hat{u}+\hat{v}|\hat{u}|^{2}.\end{array} \tag{3}\]

By rewriting \(\hat{u}=\hat{p}+\mathrm{i}\hat{q}\) and \(\hat{v}=\hat{r}+\mathrm{i}\hat{s}\) with the real-valued functions \((\hat{p},\hat{q})\) and \((\hat{r},\hat{s})\), the above models are converted into the following real-valued PINNs:

\[\begin{array}{l}f_{p}:=\hat{p}_{x}+\hat{s}+\hat{q}(\hat{r}^{2}+\hat{s}^{2}),\;\;f_{q}:=-\hat{q}_{x}+\hat{r}+\hat{p}(\hat{r}^{2}+\hat{s}^{2}),\\ f_{r}:=f_{t}+\hat{q}+\hat{s}(\hat{p}^{2}+\hat{q}^{2}),\;\;f_{s}:=-\hat{s}_{t}+ \hat{p}+\hat{r}(\hat{p}^{2}+\hat{q}^{2}),\end{array} \tag{4}\]

which possess two inputs \((x,t)\) and four outputs \((\hat{p},\hat{q},\hat{r},\hat{s})\).

### Methodology

#### 2.1.1 Physics-informed neural networks

Based on the algorithm of PINNs, we firstly establish a fully connected NN of depth \(L\) with an input layer, \(L-1\) hidden-layers and an output layer. In the \(l\)th hidden-layer, we assume that it possesses \(N_{l}\) neurons, then the \(l\)th layer and the previous layer can be connected by using the affine transformation \(\mathcal{A}_{l}\) and the activation function \(\sigma\):

\[\mathbf{x}^{l}=\sigma(\mathcal{A}_{l}(\mathbf{x}^{l-1}))=\sigma(\mathbf{W}^{l }\mathbf{x}^{l-1}+\mathbf{b}^{l}), \tag{5}\]

where the weight matrix and bias vector are denoted by \(\mathbf{W}^{l}\in R^{N_{l}\times N_{l-1}}\), \(\mathbf{b}^{l}\in R^{N_{l}}\) respectively. Hence, the whole NN can be written as

\[\hat{\mathcal{U}}(\mathbf{x}\Theta)=(\mathcal{A}_{L}\circ\sigma\circ\mathcal{ A}_{L-1}\circ\cdots\circ\sigma\circ\mathcal{A}_{1})(\mathbf{x}), \tag{6}\]

where \(\hat{\mathcal{U}}(\mathbf{x};\Theta)\) represents four outputs \((\hat{p},\hat{q},\hat{r},\hat{s})\) for predicting the solutions \((p,q,r,s)\), and \(\mathbf{x}\) stands for two inputs \((x,t)\). The set \(\Theta=\{\mathbf{W}^{l},\mathbf{b}^{l}\}_{l=1}^{L}\in\hat{\mathcal{P}}\) contains the trainable parameters with \(\hat{\mathcal{P}}\) being the parameter space. Together with the initial and boundary value conditions, the loss function, which consists of different types of mean squared error (MSE), is defined as

\[\mathcal{L}(\Theta)=Loss=\mathbf{W}_{R}MSE_{R}+\mathbf{W}_{IB}MSE_{IB}, \tag{7}\]
where \(\mathbf{W}_{R}\), \(\mathbf{W}_{IB}\) are the adjustable weights for the residual and the initial-boundary terms respectively. The terms \(MSE_{R}\) and \(MSE_{IB}\) are given by

\[\begin{split} MSE_{R}=\frac{1}{N_{R}}\sum_{i=1}^{N_{R}}\left(\left| f_{p}(x_{R}^{i},t_{R}^{i})\right|^{2}+\left|f_{q}(x_{R}^{i},t_{R}^{i})\right|^{2}+ \left|f_{r}(x_{R}^{i},t_{R}^{i})\right|^{2}+\left|f_{s}(x_{R}^{i},t_{R}^{i}) \right|^{2}\right),\\ MSE_{IB}=\frac{1}{N_{IB}}\sum_{i=1}^{N_{B}}\left(\left|\hat{p}(x_{ IB}^{i},t_{IB}^{i})-p^{i}\right|^{2}+\left|\hat{q}(x_{IB}^{i},t_{IB}^{i})-q^{i} \right|^{2}+\left|f(x_{IB}^{i},t_{IB}^{i})-r^{i}\right|^{2}+\left|s(x_{IB}^{i},t_{IB}^{i})-s^{i}\right|^{2}\right),\end{split} \tag{8}\]

where \(\left\{x_{R}^{i},t_{R}^{i}\right\}_{i=1}^{N_{R}}\) is the set of random residual points, \(\left\{x_{IB}^{i},t_{IB}^{i}\right\}_{i=1}^{N_{IB}}\) is the set of training points randomly selected from the dataset of initial and boundary value conditions with \((p^{i},q^{i},r^{i},s^{i})\equiv\)\(\left[p(x_{IB}^{i},t_{IB}^{i}),q(x_{IB}^{i},t_{IB}^{i})\right]\), \(r(x_{IB}^{i},t_{IB}^{i})\), \(s(x_{IB}^{i},t_{IB}^{i})\)]. The \(N_{IB}\) and \(N_{R}\) represent the number of points in both terms of loss functions, the term \(MSE_{R}\) is used to measure the deviation from the differential equations of the MT model on the collocation points with the help of automatic differentiation, while the term \(MSE_{IB}\) enforces the given initial-boundary value conditions as a constraint for ensuring the well-posedness of the MT model. The PINN algorithm is designed to find the optimized parameters in the set \(\Theta\) by minimizing the total loss function. Consequently, the NN outputs generate the approximate solutions \(\hat{u}=\hat{p}+\mathrm{i}\hat{q}\) and \(\hat{v}=\hat{r}+\mathrm{i}\hat{s}\), which not only obey the differential equations in the MT model, but also satisfy the initial-boundary value conditions in Eq.(2).

#### 2.1.2 Extended physics-informed neural networks with interface zones

The basic idea of XPINNs is to divide the computational domain into a number of subdomains and to employ a separate PINN in each subdomain [11]. Then the training solutions that are subject to the proper interface conditions in the decomposed subdomains are stitched together to obtain a solution over the entire domain. In XPINNs and cPINNs [10; 11], a sufficient number of points need to be selected from a common interface line between two adjacent domains to substitute into the interface conditions. Here we slightly modify this interface line into a small interface zone that is jointly located in two adjacent domains. The main reason for this modification is that apart from the pseudo initial and residual conditions, especially the gradient condition will be introduced more reasonably as interface conditions for the smoothness of learned solutions. This simple improvement will result in better stitching performance and faster convergence as shown below.

To apply XPINNs with small interface zones into our model (2), the whole domain is firstly divided into \(N_{s}\) subdomains along the time interval \([T_{0},T_{1}]=[t_{0},t_{n}]=[t_{0},t_{1}]^{0}\bigcup[t_{1},t_{2}]^{1}\bigcup \cdots\bigcup[t_{n-1},t_{n}]^{N_{s}-1}\). The \(k\)-th subdomain (\(k\)=\(0,1,\cdots,N_{s}\)\(-1\)=\(n-1\)) is described as the set: \(\Omega_{k}\equiv\{(x,t)|x\in[x_{lb},x_{ub}]=[L_{0},L_{1}],t\in[t_{k},t_{k+1}]\}\). Then the small interface zones are introduced between two adjacent domains. As our division is along the time variable, the neighborhood of \(t_{k}\) is able to decide the interface zone: \(\Delta\Omega_{k}\equiv\{(x,t)|x\in[x_{lb},x_{ub}],t\in[t_{k}-e_{k}^{-},t_{k}+ e_{k}^{+}]\}\) between \(k-1\)-th subdomain and \(k\)-th one, where \(e_{k}^{+}\) and \(e_{k}^{-}\) are small parameters. Especially, if \(e_{k}^{\pm}=0\), the interface zones reduce to the interface lines of the original XPINNs. The schematic diagram of such domain decomposition with interface zones is displayed in Fig. 1.

Figure 1: The illustrated domain decomposition and distribution of sampling points at each stage.
This domain decomposition suggests that \(N_{s}\) stages of individual PINNs need to be employed gradually, where the \(k\)-th PINN is implemented in the \(k\)-th subdomain correspondingly. For the 0-th stage, the original PINN is used to train the model in the 0-th subdomain, where the loss function is the same as Eq.(7) without the additional terms. The dataset of the initial value is still generated from the initial value condition in the whole model (2), but one of the boundary value is obtained from the boundary value functions with the time subinterval \(t\in[t_{0},t_{1}]\).

For the \(k\)-th stage \((k>0)\), the loss function of PINN in the \(k\)-th subdomain is redefined by adding the interface conditions

\[\mathcal{L}(\Theta_{k})=\textit{Loss}_{k}=\textbf{W}_{R}^{(k)}MSE_{R}^{(k)}+ \textbf{W}_{PIB}^{(k)}MSE_{PIB}^{(k)}+\textbf{W}_{IR}^{(k)}MSE_{IR}^{(k)}+ \textbf{W}_{IG}^{(k)}MSE_{IG}^{(k)}, \tag{9}\]

where \(\textbf{W}_{R}^{(k)}\), \(\textbf{W}_{PIB}^{(k)}\), \(\textbf{W}_{IR}^{(k)}\) and \(\textbf{W}_{IG}^{(k)}\) are the alterable weights for the residual, the pseudo initial-boundary and interface terms where the last two weights correspond to residual and solution smoothness in the interface zone. The four terms of MSE are expressed as follows:

\[MSE_{R}^{(k)}=\frac{1}{N_{R}^{(k)}}\sum\limits_{i=1}^{N_{R}^{(k) }}\left(\left|f_{p}(x_{R_{k}^{\prime}}^{i},t_{R_{k}}^{i})\right|^{2}+\left|f_{ q}(x_{R_{k}^{\prime}}^{i},t_{R_{k}}^{i})\right|^{2}+\left|f_{r}(x_{R_{k}^{ \prime}}^{i},t_{R_{k}}^{i})\right|^{2}+\left|f_{s}(x_{R_{k}^{\prime}}^{i},t_{R _{k}}^{i})\right|^{2}\right), \tag{10}\] \[MSE_{PIB}^{(k)}=\frac{1}{N_{PIB}^{(k)}}\sum\limits_{i=1}^{N_{PIB }^{(k)}}\left(\left|\beta_{PIB_{k}^{\prime}}^{i},t_{PIB_{k}}^{i}\right)-p_{k}^{ i}\right|^{2}+\left|\vartheta_{PIB_{k}^{\prime}}^{i},t_{PIB_{k}}^{i} \right)-q_{k}^{i}\right|^{2}+\left|\hat{r}(x_{PIB_{k}^{\prime}}^{i},t_{PIB_{k} }^{i})-r_{k}^{i}\right|^{2}\] \[+\left|\delta(x_{PIB_{k}^{\prime}}^{i},t_{PIB_{k}}^{i})-s_{k}^{i }\right|^{2})\,,\] \[MSE_{IR}^{(k)}=\frac{1}{N_{IR}^{(k)}}\sum\limits_{i=1}^{N_{R}^{( k)}}\left(\left|f_{p}(x_{IR_{k}^{\prime}}^{i},t_{IR_{k}}^{i})-f_{p^{-}}^{i} \right|^{2}+\left|f_{q}(x_{IR_{k}^{\prime}}^{i},t_{IR_{k}}^{i})-f_{q^{-}}^{i} \right|^{2}+\left|f_{r}(x_{IR_{k}^{\prime}}^{i},t_{IR_{k}}^{i})-f_{r^{-}}^{i} \right|^{2}\right.\] \[\left.+\left|f_{s}(x_{IR_{k}^{\prime}}^{i},t_{IR_{k}}^{i})-f_{s} ^{i}\right|^{2}\right),\] \[MSE_{IG}^{(k)}=\frac{1}{N_{IG}^{(k)}}\sum\limits_{i=1}^{N_{IG}^{( k)}}\left(\left|\partial_{t}p(x_{IO_{k}^{\prime}}^{i},t_{IG_{k}}^{i})-\partial_{t}p^{i -}\right|^{2}+\left|\partial_{t}q(x_{IG_{k}^{\prime}}^{i},t_{IG_{k}}^{i})- \partial_{t}q^{i-}\right|^{2}+\left|\partial_{t}r(x_{IG_{k}^{\prime}}^{i},t_{IG _{k}}^{i})-\partial_{t}r^{i-}\right|^{2}\right.\] \[+\left|\partial_{t}s(x_{IG_{k}^{\prime}}^{i},t_{IG_{k}}^{i})- \partial_{t}s^{i-}\right|^{2})\,,\]

where \(\left\{x_{R_{k}^{\prime}}^{i},t_{R_{k}}^{i}\right\}_{i=1}^{N_{R}^{(k)}}\) is the set of random residual points in the \(k\)-th subdomain, \(\left\{x_{PIB_{k}^{\prime}}^{i},t_{PIB_{k}}^{i}\right\}_{i=1}^{N_{PIB}^{(k)}}\) is the set of training points randomly selected from the dataset of pseudo initial-boundary value conditions. The values of \((p_{k}^{i},q_{k}^{i},r_{k}^{i},s_{k}^{i})\) contain two parts: If the training points belong to the sub boundary, i.e., \((x_{PIB_{k}}^{i},t_{PIB_{k}}^{i})\in\left\{(x,t)|x=x_{lb},x_{ub},t\in[t_{k},t_{ k+1}]\right\}\), then \(p_{k}^{i}+\text{i}q_{k}^{i}=u_{lb}(t_{PIB_{k}}^{i}),u_{ub}(t_{PIB_{k}}^{i})\) and \(r_{k}^{i}+\text{i}s_{k}^{i}=v_{lb}(t_{PIB_{k}}^{i}),v_{ub}(t_{PIB_{k}}^{i})\): If

Fig. 2: The schematic of PINN for the MT model employed at a stage, where NN and physics-informed network share hyper-parameters with additional loss terms in the interface zone.
the training points belong to the interface zone, i.e., \((x^{i}_{PIB_{k}},t^{i}_{PIB_{k}})\in\Delta\Omega_{k}\), then \((p^{i}_{k^{\prime}},q^{i}_{k^{\prime}},t^{i}_{k^{\prime}},s^{i}_{k})\) are called pseudo initial values obtained from the predicting solution outputs \((p^{i}_{k-1},q^{i}_{k-1},p^{i}_{k-1},s^{i}_{k-1})\) at \((x^{i}_{PIB_{k}},t^{i}_{PIB_{k}})\) through the \((k-1)\)-th PINN with the last optimized parameters.

Furthermore, \(\left\{x^{i}_{IR_{k}},t^{i}_{IR_{k}}\right\}_{i=1}^{N_{IR}^{(k)}}\) and \(\left\{x^{i}_{IG_{k}},t^{i}_{IG_{k}}\right\}_{i=1}^{N_{IG}^{(k)}}\) denote the set of randomly selected residual points and gradient points respectively in the interface zone \(\Delta\Omega_{k}\). The \((f^{i}_{p^{\prime}},f^{i}_{q^{\prime}},f^{i}_{r^{\prime}},f^{i}_{s^{-}})\) and \((\partial_{t}p^{i-},\partial_{t}q^{i-},\partial_{t}r^{i-},\partial_{t}s^{i-})\) represent residual and gradient outputs at \((x^{i}_{IR_{k}},t^{i}_{IR_{k}})\) and \((x^{i}_{IG_{k}},t^{i}_{IG_{k}})\) respectively from the \(k\)\(-\)\(1\)-th PINN. Notice that due to the domain decomposition along the time direction, only the gradients with respect to the variable \(t\) need to be computed and then be imposed as one of interface conditions. The \(N_{R}^{(k)}\), \(N_{PIB}^{(k)}\), \(N_{IR}^{(k)}\) and \(N_{IG}^{(k)}\) represent the number of points in the corresponding sets. The distributions of these sampling points are illustrated in Fig. 1, where the (pseudo) initial-boundary points and residual points for each subdomain are marked with different color crosses and green dots respectively. In the interface zone \(\Delta\Omega_{k}\), residual and gradient points are distinguished by using red and blue dots. In fact, we can see clearly that the pseudo initial points are also taken in the interface zone \(\Delta\Omega_{k}\) and the relevant MSE term can be equivalently seen as one of interface conditions at \(k\)-th stage.

For four terms in the loss function (9), the basic roles of \(MSE_{R}^{(k)}\) and \(MSE_{PIB}^{(k)}\) are same as before described in the PINN algorithm on the \(k\)-th subdomain. The \(MSE_{IR}^{(k)}\) and \(MSE_{IG}^{(k)}\) correspond to the residual continuity and \(C^{1}\) solution continuity across the interface zone which is associated with two different sub-networks jointly. The last two terms and the relevant MSE for pseudo initial points acting as the interface conditions are responsible for transmitting the physical information from the \(k\)\(-\)\(1\)-th subdomain to the \(k\)-th one.

To better understand the XPINN algorithm, a schematic of modified PINN in a subdomain is shown in Fig. 2. Apart from the contributions of the NN part and the physics-informed part, and the loss function is changed by adding certain interface conditions that ensure the quality of stitching and improve the convergence rate. By minimizing the loss function below a given tolerance \(\epsilon\) or up to a prescribed maximum number of iterations, we can find the optimal values of weights \(\mathbf{W}\) and biases \(\mathbf{b}\). Then the outputs: \((I,f,\partial_{t})\circ(\hat{p},\hat{q},\hat{r},\hat{s})\) represent three parts: \(I\circ(\hat{p},\hat{q},\hat{r},\hat{s})=(\hat{p},\hat{q},\hat{r},\hat{s})\) are predicting solutions in the subdomain, \(f\circ(\hat{p},\hat{q},\hat{p},\hat{s})=(f_{p^{\prime}},f_{q^{\prime}},f_{q},f_ {q})\) and \(\partial_{t}\circ(\hat{p},\hat{q},\hat{r},\hat{s})=(\partial_{t}\hat{p}, \partial_{t}\hat{q},\partial_{t}\hat{r},\partial_{t}\hat{s})\) are the residuals and gradients on random points from the common interface zone. Finally, these predicting solutions in each subdomain are stitched together and the training errors are measured by the relative \(L^{2}\)-norm errors:

\[\mathrm{Err}r_{u}=\frac{\sqrt{\sum\limits_{i=1}^{N}|u_{i}-\partial_{i}|^{2}}}{ \sqrt{\sum\limits_{i=1}^{N}|u_{i}|^{2}}}\,,\qquad\mathrm{Err}r_{v}=\frac{ \sqrt{\sum\limits_{i=1}^{N}|v_{i}-\partial_{i}|^{2}}}{\sqrt{\sum\limits_{i=1}^ {N}|v_{i}|^{2}}}.\]

In the above PINN algorithm, the Adam and limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimizers are employed to train parameters for minimizing the loss function, where the former is one version of the stochastic gradient descent method and the latter is a full-batch gradient descent optimization algorithm with the aid of a quasi-Newton method. For smooth and regular solutions, the Adam optimizer is first-order accurate and robust, while the L-BFGS optimizer can seek a better solution due to its second-order accuracy. Therefore, one would first apply the Adam optimiser to achieve a small value of the loss function, and then switch to the L-BFGS optimiser to pursue a smaller loss function [4, 16]. In addition, weights are provided with Xavier initialization, biases are initialized to zeros and the activation function is selected as the hyperbolic tangent (tanh) function. Particularly, it is better to take the optimized weights and biases from the previous stage are as the initialization for the next stage in the XPINN algorithm when both networks have same hidden layers and neurons per layer. Therefore, in the following experiments for several types different data-driven solutions, we uniformly construct the NN with 7 hidden layers and 40 neurons per layer in each subdomain. The recommended weights in the loss function are taken as \(\mathbf{W}_{R}=\mathbf{W}_{IB}=1\) and \(\mathbf{W}_{R}^{(k)}=\mathbf{W}_{PIB}^{(k)}=1\), \(\mathbf{W}_{IR}^{(k)}=0.00001\), \(\mathbf{W}_{IG}^{(k)}=0.0001\) in these numerical simulations. At each stage, we first use the 5000 steps Adam optimization (1000 steps for first stage) with the default learning rate \(10^{-3}\) and then set the maximum iera
tions of the L-BFGS optimization to 50000. It is mentioned that in the L-BFGS optimization the iteration will stop when

\[\frac{L_{k}-L_{k+1}}{\text{max}\{|L_{k}|,|L_{k+1}|,1\}}\leq e_{m},\]

where \(L_{k}\) denotes loss in the \(n\)th step L-BFGS optimization and \(e_{m}\) stands for the machine epsilon. Here all the code is based on Python 3.9 and Tensorflow 2.7, and the default float type is always set to 'float64'.

In what follows, we will focus on the data-driven bright, dark soliton and breather as well as rogue wave solutions of the MT model (1). These localized solutions of arbitrary \(N\) order have been derived exactly in terms of determinants by using the bilinear KP hierarchy reduction technique [84; 85]. For these types of solutions, only the data-driven lower-order cases are considered here, and the higher-order cases can be performed similarly, but with more complicated computations. In particular, we will apply the classical PINN to the simple one-soliton and first-order rogue wave solutions, and the XPINN to two-soliton, one-breather and second-order rogue wave solutions.

### Data-driven bright one- and two-soliton solutions

The formulae of the bright one-soliton solutions of the MT model (1) are expressed in the form [84]

\[u=\frac{\text{i}\alpha_{1}^{*}e^{\xi_{1}}}{\omega_{1}\left[1-\frac{\text{i} \alpha_{1}|^{2}\omega_{0}\xi_{1}+\xi_{1}^{*}}{(\omega_{1}+\omega_{1})^{2}} \right]},\qquad v=\frac{\alpha_{1}^{*}e^{\xi_{1}}}{1+\frac{\text{i}\alpha_{1} |^{2}\omega_{0}\xi_{1}+\xi_{1}^{*}}{(\omega_{1}+\omega_{1})^{2}}}, \tag{11}\]

where \(\xi_{1}=\omega_{1}x-\frac{1}{\omega_{1}}t+\xi_{10}\) and \(\omega_{1}\), \(\alpha_{1}\), \(\xi_{10}\) are arbitrary complex constants. The symbol \({}^{*}\) stands for the complex conjugate hereafter.

To simulate the bright one-soliton, the parameters in exact solutions (11) are chosen as \(\alpha_{1}=1\), \(\omega_{1}=1+3\text{i}\) and \(\xi_{10}=0\). The intervals of the computational domain \([L_{0},L_{1}]\) and \([T_{0},T_{1}]\) are taken as \([-2,2]\) and \([-3,3]\), respectively. Then the initial and Dirichlet boundary conditions in (2) are presented explicitly. By taking the grid points for the space-time region as \(400\times 600\) with the equal step length, we get the data sets of the discretized initial-boundary value. For this simple one-soliton solutions, we employ the original PINN algorithm to conduct the numerical experiment. Here, \(N_{IB}=1000\) training points are randomly selected from the initial-boundary data sets and \(N_{R}=20000\) collocation points are generated in the whole domain by using the Latin Hypercube Sampling (LHS) method. By optimizing the learnable parameters continually in the NN, we arrive at the learned bright one-solitons \(\hat{u}\) and \(\hat{v}\). The relative \(L^{2}\)-norm errors for \((|\hat{u}|,|\hat{v}|)\) are \((1.183\text{e}{-03},7.085\text{e}{-04})\). Figs. 3(a) and 3(b) exhibit the three-dimensional (3D) structure and two-dimensional (2D) density profiles of the learned bright one-soliton for the solutions \((|\hat{u}|,|\hat{v}|)\), respectively. Fig. 3(c) shows the 2D density profiles of the point-wise absolute errors between exact and learned solutions. The maximum point-wise absolute errors are \((3.621\text{e}{-03},2.046\text{e}{-02})\), which appear at near the central peaks of soliton with the high gradients.

The expressions for the bright two-soliton solutions of the MT model (1) are given by

\[u=\frac{g}{f^{*}}\,,\quad v=\frac{h}{f}\,, \tag{12}\]

where tau functions \(f\), \(g\) and \(h\) are defined as

\[\begin{split}& f=1+c_{11}\text{e}^{\xi_{1}+\xi_{1}^{*}}+c_{21} \text{e}^{\xi_{2}+\xi_{1}^{*}}+c_{12}\text{e}^{\xi_{1}^{*}+\xi_{2}^{*}}+c_{22 \text{e}^{\xi_{2}^{*}+\xi_{2}^{*}}}+c_{121}\text{e}^{\xi_{1}+\xi_{2}+\xi_{1}^{ *}+\xi_{2}^{*}},\\ & g=\frac{\text{i}\alpha_{1}^{*}}{\omega_{1}}e^{\xi_{1}}+\frac{ \text{i}\alpha_{2}^{*}}{\omega_{2}}e^{\xi_{2}}-\frac{\text{i}\omega_{1}^{*}}{ \omega_{1}\omega_{2}}c_{121}\text{e}^{\xi_{1}+\xi_{2}+\xi_{1}^{*}}-\frac{\text {i}k_{2}^{*}}{\omega_{1}\omega_{2}}c_{122}\text{e}^{\xi_{1}+\xi_{2}+\xi_{2}^{ *}},\\ & h=\alpha_{1}^{*}e^{\xi_{1}}+\alpha_{2}^{*}e^{\xi_{2}}+c_{121} \text{e}^{\xi_{1}+\xi_{2}+\xi_{1}^{*}}+c_{122}\text{e}^{\xi_{1}+\xi_{2}+\xi_{2 }^{*}}\,,\end{split} \tag{13}\]
with

\[c_{ij^{*}}=\frac{\mathrm{i}\omega_{i}\kappa_{i}^{*}\kappa_{j}}{ \left(\omega_{i}+\omega_{j}^{*}\right)^{2}},\quad c_{12i^{*}}=\left(\omega_{1}- \omega_{2}\right)\omega_{i}^{*}\left[\frac{\alpha_{2}^{*}c_{1i^{*}}}{\omega_{1} \left(\omega_{2}+\omega_{i}^{*}\right)}-\frac{\alpha_{1}^{*}c_{2i^{*}}}{\omega_ {2}\left(\omega_{1}+\omega_{i}^{*}\right)}\right]\,,\] \[c_{121^{*}2^{*}}=\left|\omega_{1}-\omega_{2}\right|^{2}\left[ \frac{c_{11^{*}}c_{22^{*}}}{\left(\omega_{1}+\omega_{2}^{*}\right)\left( \omega_{2}+\omega_{1}^{*}\right)}-\frac{c_{12^{*}}c_{21^{*}}}{\left(\omega_{1 }+\omega_{1}^{*}\right)\left(\omega_{2}+\omega_{2}^{*}\right)}\right]\,.\]

Here \(\xi_{i}=\omega_{i}x-\frac{1}{\omega_{i}^{*}}t+\xi_{i0}\) and \(\omega_{i}\), \(\alpha_{i}\), \(\xi_{i0}\) (\(i\)=\(1,2\)) are arbitrary complex constants.

Fig. 4: The train loss of the data-driven bright two-soliton solution in each subdomain with Adam (blue) and L-BFGS (orange) optimizations.

Fig. 3: The data-driven bright one-soliton solution of the MT model (2): (a) and (b) The 3D and 2D density profiles of the learned one-soliton, respectively; (c) The 2D density profiles of the point-wise absolute errors between exact and learned solutions, \(\mathrm{Error}_{a}=\left|u-\hat{a}\right|\) and \(\mathrm{Error}_{b}=\left|v-\hat{v}\right|\).
For the bright two-soliton solutions, the parameters are taken as \(\alpha_{1}=\alpha_{2}=\frac{3}{2}\), \(\omega_{1}=\frac{3}{2}+\frac{3}{2}\)i, \(\omega_{2}=\frac{3}{5}+\frac{1}{2}\)i and \(\zeta_{10}=\zeta_{20}=0\). The computational domain \([L_{0},L_{1}]\) and \([T_{0},T_{1}]\) are defined as \([-5,5]\) and \([-3,3]\) respectively, then we have the explicit initial and Dirichlet boundary conditions in (2). The grid points for the space-time region are set to \(600\times 1200\) with the equal step length, the data sets of the discretized initial-boundary value can be obtained. Compared with the bright one-soliton solutions, the bright two-soliton solutions possess the relatively complicated structures. These features are not well captured by the traditional PINN algorithm, so we must apply the XPINN algorithm to perform this numerical experiment. The computational domain is divided to 5 subdomains as shown in Table 1, and \([e_{k}^{-},e_{k}^{+}]\) are chosen uniformly as \([5,3]\times\frac{3-(-3)}{1200}\) with \(k=1,2,3,4\) to decide four small interface zones. At each stage, the numbers of training points \(N_{PIB}^{(k)}\) are given in Table 1 respectively, and the numbers of collocation points \(N_{R}^{(k)}\) are set to \(N_{R}^{(k)}=20000\) uniformly. These points are generated by using the LHS method. To guarantee the smooth stitching at the line \(t=t_{k}\), it is necessary to take residual points and gradient points with the equal grids in the interface zone \(\Delta\Omega_{k}\). Hence, we simply take such two types of points as the corresponding part of the grid points in the process of discretising the exact solution. This implies that here the numbers of both types of sample points are \(N_{IR}^{(k)}=N_{IG}^{(k)}=8\times 600=4800\). The learned solutions in each subdomain are firstly obtained by training the learnable parameters in each NN. The loss function plots at each stage are given in Fig. 4. Furthermore, by stitching together these predicted solutions in sequence, we have the learned bright two-solitons solutions \(\hat{u}\) and \(\hat{v}\). The relative \(L^{2}\)-norm errors for \((|\hat{u}|,|\hat{v}|)\) in each subdomain and the whole domain are given in Table 1. The 3D structure and 2D density profiles of the learned bright two-soliton for the solutions \((|\hat{u}|,|\hat{v}|)\) in the whole domain are displayed in Figs. 5(a) and 5(b), respectively. The 2D density profiles of the point-wise absolute errors are depicted in Fig. 5(c), in which the maximum values (\(2.312\mathrm{e}{-}02,1.559\mathrm{e}{-}02\)) occur at near the central peaks of soliton. As we can see from Fig. 5(c), despite the mismatches appear in the high gradient regions, the slight mismatches along the interface line emerge due to the different sub-networks. To reduce such mismatch for smooth stitching, one can change the numbers of interface points, the weights of interface terms in the loss function adequately or introduce the higher-order gradient interface conditions.

For the purpose of comparison, we first apply the PINN algorithm directly to learn the bright two-soliton solution. The network structures and numbers of sample points are taken to be the same as in the numerical experiment for the bright one-soliton. But the Adam optimization and the maximum iterations of the L-BFGS optimization are adjusted to 20000 and 50000 respectively. The results show that the relative \(L^{2}\)-norm errors and the maximum point-wise absolute errors for \((|\hat{u}|,|\hat{v}|)\) reach at (\(4.533\mathrm{e}{-}01,5.376\mathrm{e}{-}01\)) and (\(2.951,3.141\)),

Fig. 5: The data-driven bright two-soliton solution of the MT model (2). (a) and (b) The 3D and 2D density profiles of the learned two-soliton, respectively; (c) The 2D density profiles of the point-wise absolute errors between exact and learned solutions, \(\mathrm{Error}_{u}=|u-\hat{u}|\) and \(\mathrm{Error}_{v}=|v-\hat{v}|\).
respectively. Thus, the original PINN seems to be very poor at predicting the bright two-soliton solution. On the other hand, we conduct the numerical experiment for the bright two-soliton by using the XPINN algorithm with interface lines, namely \([\epsilon_{k}^{-},\epsilon_{k}^{+}]=[0,0]\), where the domain decomposition, all network settings and numbers of sample points are the same as ones of Fig. 5. For training results, the relative \(L^{2}\)-norm errors for \((|\hat{u}|,|\hat{v}|)\) in each subdomain and the whole domain are listed in Table 2 and the maximum values of point-wise absolute errors reach at \((2.410\text{e}{-02},2.594\text{e}{-02})\). Hence, the XPINN algorithm with small interface zones performs slightly better in predicting effects than one with interface lines.

Furthermore, we discuss the influence of the size of small interface zone \(\Delta\Omega_{k}\), which is essentially governed by the thickness \([\epsilon_{k}^{-},\epsilon_{k}^{+}]\). Here we write \([\epsilon_{k}^{-},\epsilon_{k}^{+}]=[h_{k}^{-},h_{k}^{+}]\times\frac{3-(-3)}{ 1200}\) with integers \(h_{k}^{-}\) and \(h_{k}^{+}\). In order to preserve the smoothness of stitching at the line \(t=t_{k}\) thoroughly, we need to choose uniform sampling points with the same grid width. Hence, residual points and gradient points in the interface zone \(\Delta\Omega_{k}\) are taken the same as ones of the grid points in the discretisation of exact solution, which suggests that the numbers of both types of points can be calculated as \(N_{IR}^{(k)}=N_{IG}^{(k)}=(h_{k}^{-}+h_{k}^{+})\times\)600. Then, with respect to different size of interface zone \(\Delta\Omega_{k}\), we perform the training experiments for the bright two-soliton, where the domain decomposition, network architectures, the settings of training step and numbers of sample points are the same as ones in Fig. 5. According to the results shown in Table 3, one can see that relative \(L^{2}\) norm errors increase gradually as the thickness of interface zone enlarges.
### Data-driven dark one- and two-soliton solutions

The dark one-soliton solutions of the MT model (1) are written as [84]

\[u=\rho_{1}e^{\frac{i\phi}{\phi}}\frac{1-\frac{\mathrm{i}\omega_{1}^{*}K_{1}}{ \omega_{1}+\omega_{1}^{*}}e^{\xi_{1}+\xi_{1}^{*}}}{1-\frac{\mathrm{i}\omega_{1 }^{*}}{\omega_{1}+\omega_{1}^{*}}e^{\xi_{1}+\xi_{1}^{*}}},\qquad v=\rho_{2}e^{ \frac{1}{\phi}}\frac{1-\frac{\mathrm{i}\omega_{1}^{*}H_{1}}{\omega_{1}+\omega_ {1}^{*}}e^{\xi_{1}+\xi_{1}^{*}}}{1+\frac{\mathrm{i}\omega_{1}}{\omega_{1}+ \omega_{1}^{*}}e^{\xi_{1}+\xi_{1}^{*}}},\qquad\phi=\rho\left(\frac{\rho_{2}}{ \rho_{1}}x+\frac{\rho_{1}}{\rho_{2}}t\right), \tag{14}\]

with

\[K_{1}=-\frac{\omega_{1}-\mathrm{i}\alpha}{\omega_{1}^{*}+\mathrm{i}\alpha}, \qquad H_{1}=-\frac{\omega_{1}-\mathrm{i}\alpha(1+\rho_{1}\rho_{2})}{\omega_{ 1}^{*}+\mathrm{i}\alpha(1+\rho_{1}\rho_{2})},\qquad\xi_{1}=\frac{\rho_{2}}{ \alpha\rho_{1}}p_{1}x-\frac{\rho_{1}\alpha\rho}{\rho_{2}}\frac{t}{\omega_{1}} +\xi_{10}.\]

Here \(\alpha\), \(\rho_{1}\), \(\rho_{2}\) are real constants with the restriction that \(\rho=1+\rho_{1}\rho_{2}\neq 0\), \(\omega_{1}\), \(\xi_{10}\) are complex constants. These parameters have to satisfy the following algebraic constraint:

\[|\omega_{1}-\mathrm{i}\alpha\rho|^{2}=\alpha^{2}\rho_{1}\rho_{2}\rho\,.\]

For the dark one-soliton solutions, the parameter values \(\rho_{1}=\rho_{2}=\alpha=1\), \(\omega_{1}=-\sqrt{2}+2\mathrm{i}\) and \(\xi_{10}=0\) are chosen for the following numerical simulation. In the computational domain, the intervals \([L_{0},L_{1}]\) and \([T_{0},T_{1}]\) correspond to \([-3,3]\) and \([-3,3]\). The initial and Dirichlet boundary conditions in (2) are obtained explicitly from the above exact solutions. Similarly, the data sets of the discretized initial-boundary value are provided by taking the grid points \(400\times 600\) with the equal step length in the space-time region. In this simple case, the classical PINN algorithm is used to conduct the numerical experiment for dark one-soliton solutions. Furthermore, \(N_{IB}=1000\) training points are randomly chosen from the initial-boundary data sets and \(N_{R}=20000\) collocation points are generated in the whole domain by using the LHS method. After training the learnable parameters in the NN, the predicted dark one-soliton solutions \(\hat{u}\) and \(\hat{v}\) are successfully acquired. For this solutions, the relative \(L^{2}\)-norm errors for \((|\hat{u}|,|\hat{v}|)\) reach at \((3.998\mathrm{e}{-04},7.319\mathrm{e}{-04})\). In Figs. 6(a) and 6(b), we show the 3D structure and 2D density profiles of the learned dark one-soliton solutions \((|\hat{u}|,|\hat{v}|)\), respectively. The 2D density profiles of the point-wise absolute errors are displayed in Fig. 6(c), in which the maximum point-wise absolute errors \((5.333\mathrm{e}{-03},4.492\mathrm{e}{-03})\) are found in the high gradient regions near the central peaks of soliton.

The dark two-soliton solutions of the MT model (1) are expressed as

\[u=\rho_{1}e^{\frac{i\phi}{\phi}}\frac{\mathcal{S}}{f^{*}},\qquad v=\rho_{2}e^{ \frac{i\phi}{\phi}}\frac{h}{f},\qquad\phi=\rho\left(\frac{\rho_{2}}{\rho_{1}}x+ \frac{\rho_{1}}{\rho_{2}}t\right), \tag{15}\]

\begin{table}
\begin{tabular}{c|c c c c c c} \hline Layers-neurons & 5-40 & 5-60 & 7-40 & 7-60 & 9-40 & 9-60 \\ \hline \multirow{2}{*}{\(L^{2}\) error} & \(|\hat{u}|\) & 3.151e-03 & 3.660e-03 & 1.942e-03 & 2.922e-03 & 3.794e-03 & 2.448e-03 \\  & \(|\hat{v}|\) & 3.372e-03 & 3.443e-03 & 1.992e-03 & 3.418e-03 & 3.968e-03 & 2.131e-03 \\ \hline \end{tabular}
\end{table}
Table 4: Relative \(L^{2}\)-norm errors of the data-driven bright two-soliton solution in the whole domain with different numbers of hidden layers and neurons.

\begin{table}
\begin{tabular}{c|c c c} \hline Learning rate & \(8\times 10^{-4}\) & \(10^{-3}\) & \(1.2\times 10^{-3}\) \\ \hline \multirow{2}{*}{\(L^{2}\) error} & \(|\hat{u}|\) & 3.322e-03 & 1.942e-03 & 2.766e-03 \\  & \(|\hat{v}|\) & 3.265e-03 & 1.992e-03 & 3.132e-03 \\ \hline \end{tabular}
\end{table}
Table 5: Relative \(L^{2}\)-norm errors of the data-driven bright two-soliton solution in the whole domain with different learning rates in Adam optimization.
where tau functions \(f\), \(g\) and \(h\) are defined as

\[\begin{split}& f=1+d_{11^{*}}e^{\xi_{1}+\xi_{1}^{*}}+d_{22^{*}}e^{ \xi_{2}+\xi_{2}^{*}}+d_{11^{*}}d_{22^{*}}\Omega_{12}e^{\xi_{1}+\xi_{2}+\xi_{1}^ {*}+\xi_{2}^{*}},\\ & g=1+d_{11^{*}}^{*}K_{1}e^{\xi_{1}+\xi_{1}^{*}}+d_{22^{*}}^{*}K_ {2}e^{\xi_{2}+\xi_{2}^{*}}+d_{11^{*}}^{*}d_{22^{*}}^{*}K_{1}K_{2}\Omega_{12}e^{ \xi_{1}+\xi_{2}+\xi_{1}^{*}+\xi_{2}^{*}},\\ & h=1+d_{11^{*}}^{*}H_{1}e^{\xi_{1}+\xi_{1}^{*}}+d_{22}^{*}H_{2}e^ {\xi_{2}+\xi_{2}^{*}}+d_{11^{*}}^{*}d_{22}^{*}H_{1}H_{2}\Omega_{12}e^{\xi_{1}+ \xi_{2}+\xi_{1}^{*}+\xi_{2}^{*}},\end{split} \tag{16}\]

with

\[\begin{split}& d_{ii^{*}}=\frac{\mathrm{i}\omega_{i}}{\omega_{i}+ \omega_{i}^{*}},\;\;K_{i}=-\frac{\omega_{i}-\mathrm{i}\alpha}{\omega_{i}^{*}+ \mathrm{i}\alpha},\;\;H_{i}=-\frac{\omega_{i}-\mathrm{i}\alpha\rho}{\omega_{i}^ {*}+\mathrm{i}\alpha\rho},\\ &\Omega_{12}=\frac{|\omega_{1}-\omega_{2}|^{2}}{|\omega_{1}+ \omega_{2}^{*}|^{2}},\;\;\xi_{i}=\frac{\rho_{2}}{\alpha\rho_{1}}\omega_{i}x- \frac{\alpha\rho_{1}\rho}{\rho_{2}}\frac{t}{\omega_{i}}+\xi_{i0},\;\;i=1,2. \end{split}\]

Here \(\omega_{i}\), \(\xi_{i0}\) (\(i{=}1,2\)) are complex constants and \(\alpha\), \(\rho_{1}\), \(\rho_{2}\) are real constants with \(\rho{=}1+\rho_{1}\rho_{2}\neq 0\). These parameters need to satisfy the algebraic constraints:

\[|\omega_{i}-\mathrm{i}\alpha\rho|^{2}=\alpha^{2}\rho_{1}\rho_{2}\rho.\]

Following the detailed analysis of the maximum amplitudes of the dark one-soliton [84], we know that both components \(u\) and \(v\) allow dark and anti-dark solitons. Thus the dark two-soliton solutions could support three types of collisions: dark-dark solitons, dark-anti-dark solitons, and anti-dark-anti-dark ones. For the purpose of illustration, we will proceed here to learn only two cases: the dark-dark solitons and the dark-anti-dark solitons. The parameters in both cases are chosen as \(\omega_{2}=-\sqrt{2}+2\mathrm{i}\), \(\xi_{10}=\frac{1}{4}\) for the dark-dark solitons and \(\omega_{2}=\sqrt{2}+2\mathrm{i}\), \(\xi_{10}=-\frac{1}{2}\) for the dark-anti-dark solitons as well as the common values \(\rho_{1}=\rho_{2}=\alpha=1\), \(\omega_{1}=-1+\mathrm{i}\), \(\xi_{20}=0\).

\begin{table}
\begin{tabular}{c c c c c c c} \hline \multirow{2}{*}{\(L^{2}\) error} & \(t\)-domain & \([-4,-2.5]\) & \([-2.5,-1]\) & \([-1,1]\) & \([1,2.6]\) & \([2.6,4]\) & \([-4,4]\) \\ \cline{2-7}  & \(N_{lb}+(N_{Pl})\) & 1000 & 2000 & 2500 & 3000 & 3500 & \\ \hline \(|\tilde{u}|\) & & 1.270e-03 & 2.582e-03 & 2.637e-03 & 4.937e-03 & 7.284e-03 & 4.157e-03 \\ \(|\vartheta|\) & & 3.838e-04 & 8.310e-04 & 1.874e-03 & 3.844e-03 & 7.434e-03 & 3.655e-03 \\ \hline \end{tabular}
\end{table}
Table 6: Relative \(L^{2}\)-norm errors of the data-driven dark-dark soliton solution in each subdomain and the whole domain.

Figure 6: The data-driven dark one-soliton solution of the MT model (2): (a) and (b) The 3D and 2D density profiles of the learned one-soliton, respectively; (c) The 2D density profiles of the point-wise absolute errors between exact and learned solutions, \(\mathrm{Error}_{u}=|u-\tilde{u}|\) and \(\mathrm{Error}_{v}=|v-\phi|\).
In the following data-driven experiments, the computational domains of both cases \([L_{0},L_{1}]\times[T_{0},T_{1}]\) are taken as \([-6,6]\times[-4,4]\) and \([-5,5]\times[-4,4]\), respectively. From the definitions given in (2), both types of initial-boundary conditions can be written explicitly. The data sets of the discretized initial-boundary value are generated by inserting the grid points \(600\times 1200\) with the equal step length in the space-time regions. Similar to the bright two-soliton solutions, the dark two-soliton solutions with the relatively complicated structures force us to use the XPINN algorithm for pursuing the better numerical simulation results. For both cases, we divide the whole computational domain into 5 subdomains as given in Tables 6 and 7 respectively, and define \([e_{k}^{-},e_{k}^{+}]=[5,3]\times\frac{T_{1}-T_{0}}{1200}\) uniformly with \(k=1,2,3,4\) for the corresponding small interface zones. In Tables 6 and 7, we list the numbers of training points \(N_{PIB}^{(k)}\) at each stage. The number of collocation points in each subdomain is \(N_{R}^{(k)}=20000\) uniformly, which are produced by using the LHS method. To ensure smooth stitching at the line \(t=t_{k}\), we take residual points and gradient points in each interface zone as the corresponding part of the grid points in the discretisation of exact solution. Thus the numbers of two types of sample points are \(N_{IR}^{(k)}=N_{IG}^{(k)}=4800\). After training the learnable parameters in each NN and stitching together the predicted solutions in all subdomains, both types of the dark two-soliton solutions \(\hat{a}\) and \(\hat{o}\) are obtained in the whole computational domain. In each subdomain and the whole domain, the relative \(L^{2}\)-norm errors for \((|a|,|\hat{o}|)\) are calculated in Tables 6 and 7. For the dark-dark solitons and the dark-anti-dark solitons, the 3D structure and 2D density profiles of the learned solutions as well as the 2D density profiles of the point-wise absolute errors are shown in Figs. 7 and 8 respectively. In the comparisons between exact solutions and learned ones for two cases, the maximum absolute errors of the solutions \((|\hat{a}|,|\hat{o}|)\) are \((4.940\text{e}{-}02,4.699\text{e}{-}02)\) and \((6.169\text{e}{-}02,3.082\text{e}{-}02)\) near the central peaks and valleys of solitons that correspond to the high gradient regions. In addition, we can observe that apart from the large errors occurring near the wave peaks and valleys, the adjacent sub-networks with the different architectures lead to certain errors along the interface line. This type of mismatch can be reduced by adding more interface information, such as increasing the numbers of interface points and the weights of interface terms in the loss function, and imposing the additional interface conditions.
where tau functions \(f\), \(g\) and \(h\) are given by

\[\begin{split} f&=\frac{\epsilon^{\zeta_{1}^{\intercal} \zeta_{1}^{\intercal}\zeta_{1}^{\intercal}}}{\omega_{11}+\omega_{11}^{ \intercal}}\frac{\omega_{11}}{\omega_{21}}+\frac{\epsilon^{\zeta_{1}^{ \intercal}}}{\omega_{11}+\omega_{21}^{\intercal}}\frac{\omega_{11}}{\omega_{2 1}}+\frac{\epsilon^{\zeta_{1}^{\intercal}}}{\omega_{21}+\omega_{11}^{ \intercal}}+\frac{1}{\omega_{21}+\omega_{21}^{\intercal}},\\ g&=\frac{\Theta_{1}\Theta_{1}^{\intercal}-\epsilon^{ \zeta_{1}^{\intercal}\zeta_{1}^{\intercal}}}{\omega_{11}+\omega_{11}^{ \intercal}}\frac{\omega_{11}}{\omega_{21}^{\intercal}}+\frac{\Theta_{1} \Theta_{1}^{\zeta_{1}}}{\omega_{11}+\omega_{21}^{\intercal}}+\frac{\Theta_{1} ^{\intercal-1}\epsilon^{\zeta_{1}^{\intercal}}}{\omega_{21}^{\intercal}+ \omega_{11}^{\intercal}}\frac{\omega_{11}^{\intercal}}{\omega_{21}^{ \intercal}}+\frac{1}{\omega_{21}+\omega_{21}^{\intercal}},\\ h&=\frac{\Lambda_{1}\Lambda_{1}^{\intercal}- \epsilon^{\zeta_{1}^{\intercal}\zeta_{1}^{\intercal}}}{\omega_{11}+\omega_{11}^ {\intercal}}\frac{\omega_{11}^{\intercal}}{\omega_{21}^{\intercal}}+\frac{ \Lambda_{1}\epsilon^{\zeta_{1}}}{\omega_{11}+\omega_{21}^{\intercal}}+\frac{ \Lambda_{1}^{\intercal-1}\epsilon^{\zeta_{1}^{\intercal}}}{\omega_{21}^{ \intercal}+\omega_{11}^{\intercal}}\frac{\omega_{11}^{\intercal}}{\omega_{21} ^{\intercal}}+\frac{1}{\omega_{21}+\omega_{21}^{\intercal}},\end{split} \tag{18}\]

with

\[\Theta_{1}=\frac{\omega_{11}-\mathrm{i}\alpha}{\omega_{21}-\mathrm{i}\alpha}, \qquad\Lambda_{1}=\frac{\omega_{11}-\mathrm{i}\alpha\rho}{\omega_{21}-\mathrm{ i}\alpha\rho},\qquad\zeta_{1}=(\omega_{11}-\omega_{21})\left(\frac{\rho_{2}}{ \alpha\rho_{1}}x+\frac{\rho_{1}\mathrm{a}\rho}{\rho_{2}\omega_{11}\omega_{21}}t \right)+\zeta_{10}.\]

Here \(\omega_{11}\), \(\omega_{21}\), \(\zeta_{10}\) are complex constants and \(\alpha\), \(\rho_{1}\), \(\rho_{2}\) are real constants with \(\rho=1+\rho_{1}\rho_{2}\neq 0\), in which these parameters need to obey the algebraic constraint:

\[(\omega_{11}-\mathrm{i}\alpha\rho)(\omega_{21}-\mathrm{i}\alpha\rho)=-\alpha^{2 }\rho_{1}\rho_{2}\rho.\]

As reported in Ref.[84], the general breather of the MT model (1) contains two special cases: the Akhmediev breather and the Kuznetsov-Ma breather. For the sake of illustration, we will only consider the general one

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \multirow{2}{*}{\(L^{2}\) error} & \(t\)-domain & \([-7,-4.94]\) & \([-4.94,-2.88]\) & \([-2.88,-0.13]\) & \([-0.13,2.17]\) & \([2.17,4]\) & \([-7,4]\) \\ \cline{2-9}  & \(N_{lb}+(N_{Pl})\) & 1000 & 2000 & 2500 & 3000 & 3500 & \\ \hline \(|\hat{u}|\) & & 4.582e-04 & 2.405e-03 & 4.290e-03 & 3.358e-03 & 3.174e-03 & 3.112e-03 \\ \(|\phi|\) & & 5.341e-04 & 4.362e-03 & 1.304e-02 & 6.867e-03 & 3.140e-03 & 7.575e-03 \\ \hline \end{tabular}
\end{table}
Table 8: Relative \(L^{2}\)-norm errors of the data-driven one-breather solution in each subdomain and the whole domain.

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \multirow{2}{*}{\(L^{2}\) error} & \(t\)-domain & \([-4,-2.5]\) & \([-2.5,-1]\) & \([-1,1]\) & \([1,2.6]\) & \([2.6,4]\) & \([-4,4]\) \\ \cline{2-9}  & \(N_{lb}+(N_{Pl})\) & 1000 & 2000 & 2500 & 3000 & 3500 & \\ \hline \(|\hat{u}|\) & & 6.277e-04 & 1.468e-03 & 1.834e-03 & 2.622e-03 & 4.318e-03 & 2.449e-03 \\ \(|\phi|\) & & 4.381e-04 & 8.711e-04 & 1.135e-03 & 2.161e-03 & 3.759e-03 & 1.959e-03 \\ \hline \end{tabular}
\end{table}
Table 7: Relative \(L^{2}\)-norm errors of the data-driven dark-antidark soliton solution in each subdomains and the whole domain.

Fig. 8: The data-driven dark-antidark soliton solution of the MT model (2): (a) and (b) The 3D and 2D density profiles of the learned dark-antidark soliton, respectively; (c) The 2D density profiles of the point-wise absolute errors between exact and learned solutions, \(\text{Error}_{u}=|u-\hat{u}|\) and \(\text{Error}_{v}=|v-\hat{v}|\).
breather and the Kuznetsov-Ma breather here. In two cases, the parameters are taken as \(\rho_{1}=2\rho_{2}=\alpha=1\), \(\omega_{11}=4+2\mathrm{i}\), \(\omega_{21}=-\frac{12}{65}+\frac{99}{65}\); \(\zeta_{10}=0\) for the general one-breather and \(\rho_{1}=-2\rho_{2}=1\), \(\alpha=2\), \(\omega_{11}=5+\mathrm{i}\), \(\omega_{21}=\frac{1}{5}+\mathrm{i}\), \(\zeta_{10}=0\) for the Kuznetsov-Ma breather, respectively.

For the data-driven experiments of one-breather solutions, we set the computational domain \([L_{0},L_{1}]\times[T_{0},T_{1}]\) to \([-3,3]\times[-7,4]\) and \([-3,2]\times[-6,4]\) for two examples. The initial-boundary conditions of one-breather solutions can be obtained exactly from the formulae defined in (2). For both types of breather solutions, we still choose \(600\times 1200\) grid points with the same step size to generate the data sets of the discretized initial-boundary value. Similar to the bright and dark two-soliton solutions, we need to employ the XPINN algorithm to approximate the sophisticated breather solutions. More specifically, the whole domains for both types of breather are first divided into 5 subdomains as shown in Tables 8 and 9 respectively, and then the associated small interface zones are determined by defining \([e_{k}^{-},e_{k}^{+}]=[5,3]\times\frac{T_{1}-T_{0}}{1200}\) uniformly with \(k=1,2,3,4\). The numbers of training points \(N_{PIB}^{(k)}\) at each stage are given in Tables 8 and 9, and the number of collocation points \(N_{R}^{(k)}\) in each subdomain is set to 20000 equally, which are obtained by using the LHS method. In each interface zone, residual points and gradient points are taken as the corresponding part of the grid points in the discretisation of exact solution for the smooth stitching at the line \(t=t_{k}\). Hence the numbers of two types of sample points are given by \(N_{IR}^{(k)}=N_{IG}^{(k)}=4800\). By sequentially training the learnable parameters of each NN and combining the learned solutions throughout all subdomains, we have both types of the predicted one-breather solutions \(\hat{u}\) and \(\hat{v}\) in the whole computational domain. In Tables 8 and 9, we list their relative \(L^{2}\)-norm errors for \((|\hat{u}|,|\hat{v}|)\) in each subdomain and the whole domain. Figs. 9 and 10 exhibit the 3D structure and 2D density profiles of the learned solutions as well as the 2D density profiles of the point-wise absolute errors for the general one-breather and the Kuznetsov-Ma breather, respectively. The maximum values of the absolute error for two cases are \((8.944\mathrm{e}{-02},5.460\mathrm{e}{-02})\) and \((3.080\mathrm{e}{-02},6.705\mathrm{e}{-02})\), which occur at the peaks and valleys of the local wave with the high absolute gradients. In analogy to the bright and dark two-soliton solutions, one can see that the slight mismatches along the interface line also appear due to the different sub-networks. As mentioned above, we can introduce more interface information to achieve the smoother stitching. For example, these information can be obtained by adjusting the numbers of interface points, the weights of interface terms in the loss function and the sizes of interface zones adequately, or by imposing the additional smoothness conditions such as the higher-order gradients in the NN of each subdomain.

Figure 9: The data-driven one-breather solution of the MT model (2): (a) and (b) The 3D and 2D density profiles of the learned one-breather, respectively; (c) The 2D density profiles of the point-wise absolute errors between exact and learned solutions, \(\mathrm{Error}_{\alpha}=|u-\hat{u}|\) and \(\mathrm{Error}_{\epsilon}=|v-\hat{v}|\).
### Data-driven first- and second-order rogue wave solutions

In contrast to the soliton and breather solutions shown above, the rogue wave solutions are expressed in terms of rational polynomials. According to the derivation in [85], the first-order rogue wave solutions read

\[u=\rho_{1}\dot{e}^{\dot{\imath}\phi}\frac{(X_{1}^{*}-\tilde{\rho})(Y_{1}^{*}+ \tilde{\rho}^{*})}{X_{1}^{*}Y_{1}^{*}+\frac{1}{4}},\qquad v=\rho_{2}\dot{e}^{ \dot{\imath}\phi}\frac{(X_{1}^{*}-1)(Y_{1}^{*}+1)}{X_{1}Y_{1}+\frac{1}{4}}, \qquad\phi=\rho\left(\frac{\rho_{2}}{\rho_{1}}x+\frac{\rho_{1}}{\rho_{2}}t \right), \tag{19}\]

with

\[X_{1}=\frac{\rho_{2}\dot{\rho}}{\rho_{1}}x+\frac{\rho_{1}\rho \dot{\rho}}{\rho_{2}(\dot{\jmath}+\dot{\imath}\rho)^{2}}t+\frac{1}{2}\frac{ \dot{\rho}}{\dot{\rho}+\dot{\imath}\rho},\qquad Y_{1}=\frac{\rho_{2}\dot{\rho} }{\rho_{1}}x+\frac{\rho_{1}\rho\dot{\rho}}{\rho_{2}(\dot{\rho}-\dot{\imath} \rho)^{2}}t-\frac{1}{2}\frac{\dot{\rho}}{\dot{\rho}-\dot{\imath}\rho},\] \[\tilde{\rho}=\frac{\dot{\rho}}{\dot{\rho}-\dot{\imath}\rho_{1} \rho_{2}},\ \ \tilde{\rho}=\sqrt{-\rho_{1}\rho_{2}\rho},\ \ \rho=1+\rho_{1}\rho_{2}.\]

Here \(\rho_{1}\) and \(\rho_{2}\) are arbitrary real parameters which must satisfy the condition: \(-1<\rho_{1}\rho_{2}<0\).

In the first-order rogue wave solutions, we take the parameters \(\rho_{1}=1\) and \(\rho_{2}=-\frac{4}{5}\) to perform the data-driven experiment. For the computational domain, we take the intervals \([L_{0},L_{1}]\) and \([T_{0},T_{1}]\) as \([-2.8,2.8]\) and \([-1.8,1.8]\) respectively. Then the initial and Dirichlet boundary conditions in (2) are written explicitly, and the corresponding discretized data sets are generated by imposing \(600\times 400\) grid points with the equal step length in the space-time domain. Similar to the bright and dark one-soliton solutions, here we utilize the original PINN algorithm to simulate the simple first-order rogue wave solutions. In the initial-boundary data sets and the whole domain, we randomly select \(N_{IB}=1000\) training points and \(N_{R}=20000\) collocation points respectively by using the LHS method. By optimising the learnable parameters of the NN through the iterative training, the approximated first-order rogue wave solutions \(\tilde{u}\) and \(\tilde{v}\) are successfully obtained. In this case, the relative \(L^{2}\)-norm errors for \((|\hat{u}|,|\phi|)\) are (4.283e\(-04\), 5.764e\(-04\)). Fig. 11 exhibits the 3D structure and 2D density profiles [Figs. 11(a) and 11(b)] of the learned first-order rogue wave solutions \((|\hat{u}|,|\phi|)\), and the 2D density

\begin{table}
\begin{tabular}{c c c c c c c} \hline \(L^{2}\) error & \(t\)-domain & \([-6,-4.13]\) & \([-4.13,-2.25]\) & \([-2.25,0.25]\) & \([0.25,2.33]\) & \([2.33,4]\) & \([-6,4]\) \\ \hline \(N_{lb}\)\(+\)\((N_{Pl})\) & 1000 & 2000 & 2500 & 3000 & 3500 & \\ \hline \(|\hat{u}|\) & 4.724e-04 & 1.713e-03 & 4.448e-03 & 5.553e-03 & 4.872e-03 & 4.011e-03 \\ \(|\hat{\vartheta}|\) & 4.275e-04 & 7.680e-04 & 2.944e-03 & 5.393e-03 & 4.644e-03 & 3.458e-03 \\ \hline \end{tabular}
\end{table}
Table 9: Relative \(L^{2}\)-norm errors of the data-driven Kuznetsov-Ma breather solution in each subdomain and the whole domain.

Figure 10: The data-driven Kuznetsov-Ma breather solutions of the MT model (2): (a) and (b) The 3D and 2D density profiles of the learned Kuznetsov-Ma breather, respectively; (c) The 2D density profiles of the point-wise absolute errors between exact and learned solutions, \(\text{Error}_{\pi}=|u-\hat{u}|\) and \(\text{Error}_{\pi}=|v-\hat{v}|\).
profiles of the point-wise absolute errors [Fig. 11(c)], respectively. It can be clearly seen that the predicted solutions agree well with the exact ones, except for the maximum absolute errors (5.316e\(-\)03, 3.203e\(-\)03), which appear at the central peaks of rogue wave with the high gradients.

The second-order rogue wave solutions are given by [85]

\[u=\rho_{1}\mathrm{e}^{\mathrm{i}\phi}\frac{g}{f^{*}},\qquad v=\rho_{2}\mathrm{e }^{\mathrm{i}\phi}\frac{h}{f},\qquad\phi=\rho\left(\frac{\rho_{2}}{\rho_{1}}x+ \frac{\rho_{1}}{\rho_{2}}t\right), \tag{20}\]

where tau functions \(f\), \(g\) and \(h\) are determined by

\[f=\left|\begin{array}{cc}m_{11}^{(0,0,0)}&m_{13}^{(0,0,0)}\\ m_{31}^{(0,0,0)}&m_{33}^{(0,0,0)}\end{array}\right|,\quad g=\left|\begin{array} []{cc}m_{11}^{(-1,1,0)}&m_{13}^{(-1,1,0)}\\ m_{31}^{(-1,1,0)}&m_{33}^{(-1,1,0)}\end{array}\right|,\quad h=\left|\begin{array} []{cc}m_{11}^{(-1,0,1)}&m_{13}^{(-1,0,1)}\\ m_{31}^{(-1,0,1)}&m_{33}^{(-1,0,1)}\end{array}\right|, \tag{21}\]

with the elements

\[m_{11}^{(n,k,l)} =X_{1}Y_{1}+\frac{1}{4},\quad m_{13}^{(n,k,l)}=X_{1}(Y_{3}+\frac {1}{6}Y_{1}^{3})+\frac{1}{8}Y_{1}^{2}-\frac{1}{48},\quad m_{31}^{(n,k,l)}=Y_{ 1}(X_{3}+\frac{1}{6}X_{1}^{3})+\frac{1}{8}X_{1}^{2}-\frac{1}{48},\] \[m_{33}^{(n,k,l)} =(X_{3}+\frac{1}{6}X_{1}^{3})(Y_{3}+\frac{1}{6}Y_{1}^{3})+\frac{ 1}{16}(X_{1}^{2}-\frac{1}{6})(Y_{1}^{2}-\frac{1}{6})+\frac{1}{16}X_{1}Y_{1}+ \frac{1}{64}.\]

The variables \(X_{i}\) and \(Y_{i}\)\((i=1,3)\) are as follows:

\[X_{i}=\alpha_{i}x+\beta_{i}t+(n+\frac{1}{2})\theta_{i}+k\theta_{i}+l\zeta_{r}+ a_{i},\ \ Y_{i}=\alpha_{i}x+\beta_{i}^{*}t-(n+\frac{1}{2})\theta_{i}^{*}-k\theta_{i}^{*}- l\zeta_{r}+a_{i}^{*},\]

with the coefficients:

\[a_{1} =\zeta_{3}=0,\ \ \zeta_{1}=1,\ \ \alpha_{1}=\frac{\rho_{2} \hat{\rho}}{\rho_{1}},\ \ \alpha_{3}=\frac{\rho_{2}\hat{\rho}}{\delta\rho_{1}},\quad\beta_{1}=\frac{\rho_{ 1}\rho\hat{\rho}}{\rho_{2}(\hat{\rho}+\mathrm{i}\rho)^{2}},\] \[\beta_{3} =\frac{\rho_{1}\rho(4\mathrm{i}\rho_{1}\rho_{2}\rho-2\hat{\rho} \rho_{1}\rho_{2}-\hat{\rho})}{\delta\rho_{2}(\hat{\rho}+\mathrm{i}\rho)^{4}}, \ \ \theta_{1}=\frac{\hat{\rho}}{\hat{\rho}+\mathrm{i}\rho},\ \ \theta_{3}=\frac{\mathrm{i}\rho^{2}( \rho-1+\mathrm{i}\hat{\rho})}{6(\hat{\rho}+\mathrm{i}\rho)^{3}},\] \[\theta_{1} =\frac{\hat{\rho}}{\mathrm{i}(\rho-1)+\hat{\rho}},\ \ \theta_{3}=\frac{(\mathrm{i}\rho-\hat{\rho})(\rho-1)^{2}}{6[\mathrm{i}(\rho-1)+ \hat{\rho}]^{3}},\ \ \hat{\rho}=\sqrt{-\rho_{1}\rho_{2}\rho},\ \ \rho=1+\rho_{1}\rho_{2}.\]

Here \(a_{3}\) is an arbitrary complex parameter, \(\rho_{1}\rho_{2}\) are arbitrary real parameters which need to satisfy the condi

Figure 11: The data-driven first-order rogue wave solution of the MT model (2): (a) and (b) The 3D and 2D density profiles of the learned first-order rogue wave, respectively; (c) The 2D density profiles of the point-wise absolute errors between exact and learned solutions, \(\mathrm{Error}_{s}=|u-\hat{u}|\) and \(\mathrm{Error}_{r}=|v-\hat{v}|\).
tion: \(-1<\rho_{1}\rho_{2}<0\).

According to the Ref.[85], the second-order rogue wave solutions allow the wave patterns of the sole huge peak with \(a_{3}=0\) and three peaks with \(a_{3}\neq 0\). Here, we only perform the data-driven simulation for the second-order rogue wave with three peaks, where we take the parameters \(\rho_{1}=1\), \(\rho_{2}=-\frac{4}{5}\) and \(a_{3}=\frac{1}{2}\).

For the computational domain, the intervals \([L_{0},L_{1}]\) and \([T_{0},T_{1}]\) are chosen to be \([-3,2]\) and \([-6,4]\) respectively. As defined in (2), the initial and Dirichlet boundary conditions of the second-order rogue wave solutions can be written immediately. It is straightforward to obtain the data sets of the discretized initial-boundary value by imposing \(600\times 1200\) grid points in the space-time domain with the equal step length. In a similar way to the two-soliton and beather solutions, the complicated local features of second-order rogue wave solutions force us to use the XPINN algorithm for conducting the numerical experiments. First, the whole domain is divided into 5 subdomains as shown in Table 10, and \([\epsilon_{k}^{-},\epsilon_{k}^{+}]=[5,3]\times\frac{T_{1}-T_{0}}{1200}\) with \(k=1,2,3,4\) are defined uniformly for the small interface zones. At each stage, the numbers of training points \(N_{PIB}^{(k)}\) are provided in Table 10, and the number of collocation points \(N_{R}^{(k)}\) is equally fixed at 20000. Here we employ the LHS method to generate these random training points. In each interface zone, residual points and gradient points are selected as the corresponding part of the grid points in the discretisation of exact solution to ensure smooth stitching at the line \(t=t_{k}\). Hence the numbers of two types of sample points can be obtained as \(N_{IR}^{(k)}=N_{IG}^{(k)}=4800\). By training the learnable parameters of each NN in five subdomains and stitching these predicted solutions together, the learned second-order rogue wave solutions \(\hat{u}\) and \(\hat{v}\) can be constructed in the whole domain. The relative \(L^{2}\)-norm errors for \((|\hat{u}|,|\hat{v}|)\) in each subdomain and the whole domain are presented in Table 10. In Fig. 12, we display the 3D structure and 2D density profiles of the learned solutions as well as the 2D density profiles of the point-wise absolute errors for the second-order rogue wave, respectively. From Fig. 12(c), the maximum absolute errors are measured as \((4.988\mathrm{e}{-02},1.327\mathrm{e}{-01})\), which appear at the peaks and valleys of the rogue wave. Particularly, unlike the two-soliton and breather cases, the relatively larger errors at the peaks and valleys make the mismatches along the interface line unnoticeable. This is due to the fact that the small

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \(L^{2}\) error & \(t\)-domain & \([-3,-1.69]\) & \([-1.69,-0.38]\) & \([-0.38,1.38]\) & \([1.38,2.83]\) & \([2.83,4]\) & \([-3,4]\) \\ \cline{2-7}  & \(N_{Ib}{+}(N_{PI})\) & 1000 & 2000 & 2500 & 3000 & 3500 & \\ \hline \(|\hat{u}|\) & & 1.182e-03 & 3.587e-03 & 6.940e-03 & 4.241e-03 & 3.249e-03 & 4.739e-03 \\ \(|\hat{v}|\) & & 4.739e-04 & 2.813e-03 & 4.359e-03 & 4.114e-03 & 4.281e-03 & 3.573e-03 \\ \hline \hline \end{tabular}
\end{table}
Table 10: Relative \(L^{2}\)-norm errors of the data-driven second-order rogue wave solution in each subdomain and the whole domain.

Figure 12: The data-driven second-order rogue wave solution of the MT model (2): (a) and (b) The 3D and 2D density profiles of the learned second-order rogue wave, respectively; (c) The 2D density profiles of the point-wise absolute errors between exact and learned solutions, \(\mathrm{Error}_{u}=|u-\hat{u}|\) and \(\mathrm{Error}_{v}=|\bar{v}-\bar{v}|\).
regions with higher gradients, such as the sharp peaks in Fig. 12(a), make it difficult for NN learning to capture these features. This type of absolute error can be reduced by secondary or repetitive training by selecting more collocation points in these sharp regions.

## 3 Data-driven parameters discovery of the MT model

In this section, we consider the inverse problem of the MT model where both equations in Eq.(2) are replaced by

\[\begin{split}&\mathrm{i}u_{x}+\lambda_{1}v+\lambda_{2}u|v|^{2}=0, \\ &\mathrm{i}v_{t}+\lambda_{3}u+\lambda_{4}v|u|^{2}=0,\end{split} \tag{22}\]

with unknown real coefficients \(\lambda_{i}\) for \(i=1,2,3,4\). Correspondingly, the associated complex-valued PINNs for the MT model (2) are changed into the following form

\[\begin{split}& f_{u}:=\mathrm{i}\hat{u}_{x}+\lambda_{1}\hat{v}+ \lambda_{2}\hat{u}|\hat{v}|^{2},\\ & f_{v}:=\mathrm{i}\hat{v}_{t}+\lambda_{3}\hat{u}+\lambda_{4}\hat{ v}|\hat{u}|^{2}.\end{split} \tag{23}\]

which are rewritten as the real-valued PINNs as follows:

\[\begin{split}& f_{p}:=\hat{p}_{x}+\lambda_{1}\hat{s}+\lambda_{2} \hat{q}(\hat{r}^{2}+\hat{s}^{2}),\;\;f_{q}:=-\hat{q}_{x}+\lambda_{1}\hat{r}+ \lambda_{2}\hat{p}(\hat{r}^{2}+\hat{s}^{2}),\\ & f_{r}:=\hat{r}_{t}+\lambda_{3}\hat{q}+\lambda_{4}\hat{s}(\hat{ p}^{2}+\hat{q}^{2}),\;\;f_{s}:=-\hat{s}_{t}+\lambda_{3}\hat{p}+\lambda_{4}\hat{ r}(\hat{p}^{2}+\hat{q}^{2}),\end{split} \tag{24}\]

by taking the real-valued functions \((\hat{p},\hat{q})\) and \((\hat{r},\hat{s})\) as the real and imaginary parts of \((\hat{u},\hat{v})\), respectively.

For the inverse problem, the parameters \(\lambda_{i}\) with \(i=1,2,3,4\) in Eqs.(22) are unknown, some extra information of the internal region can be introduced into the PINNs of the MT model. This implementation leads to the additional loss term defined as

\[MSE_{IN}=\frac{1}{N_{IN}}\sum_{i=1}^{N_{IN}}\left(\left|\hat{p}(x^{i}_{IN},t^{ i}_{IN})-p^{i}_{IN}\right|^{2}+\left|\hat{q}(x^{i}_{IN},t^{i}_{IN})-q^{i}_{IN} \right|^{2}+\left|\hat{r}(x^{i}_{IN},t^{i}_{IN})-r^{i}_{IN}\right|^{2}+\left| \hat{s}(x^{i}_{IN},t^{i}_{IN})-s^{i}_{IN}\right|^{2}\right), \tag{25}\]

where \(\{x^{i}_{IN},t^{i}_{IN}\}_{i=1}^{N_{IN}}\) is the set of random internal points where the number of points is \(N_{IN}\), and the discretized solutions \((p^{i}_{IN},q^{i}_{IN},r^{i}_{IN},s^{i}_{IN})\) correspond to \([p(x^{i}_{IN},t^{i}_{IN}),q(x^{i}_{IN},t^{i}_{IN}),r(x^{i}_{IN},t^{i}_{IN}),s(x^ {i}_{IN},t^{i}_{IN})]\). Then the new loss function in (7) is modified as

\[\mathcal{L}(\Theta)=Loss=\mathbf{W}_{R}MSE_{R}+\mathbf{W}_{IB}MSE_{IB}+ \mathbf{W}_{IN}MSE_{IN}, \tag{26}\]

with the weight \(\mathbf{W}_{IN}\). Consequently, the redesigned PINNs of the MT model are used to learn the unknown parameters \(\lambda_{i}\)(\(i=1,2,3,4\)) simultaneously with the solutions \((\hat{p},\hat{q},\hat{r},\hat{s})\). It is noted that here the loss term of the initial-boundary value boundary conditions can be eliminated due to the introduction of the loss term for internal information [6; 7; 86].

Based on any given solutions of the MT model (1), the data-driven experiment for the parameter discovery

\begin{table}
\begin{tabular}{c c c} \hline \hline \multirow{2}{*}{MT model} & Item & \\ \cline{2-3}  & Parameters & Relative errors \\ \hline Correct & \(\lambda_{1}=1,\;\lambda_{2}=1\) & \(\lambda_{1}:0\%,\;\lambda_{2}:0\%\) \\  & \(\lambda_{3}=1,\;\lambda_{4}=1\) & \(\lambda_{3}:0\%,\;\lambda_{4}:0\%\) \\ \hline Identified & \(\lambda_{1}=0.9999868,\;\lambda_{2}=1.000026\) & \(\lambda_{1}:0.00132\%,\;\lambda_{2}:0.00262\%\) \\ (without noise) & \(\lambda_{3}=1.000012,\;\lambda_{4}=1.000082\) & \(\lambda_{3}:0.00118\%,\;\lambda_{4}:0.00819\%\) \\ \hline Identified & \(\lambda_{1}=0.9994288,\;\lambda_{2}=0.9972788\) & \(\lambda_{1}:0.05712\%,\;\lambda_{2}:0.27212\%\) \\ (5 \% noise) & \(\lambda_{3}=1.002426,\;\lambda_{4}=1.000761\) & \(\lambda_{3}:0.24258\%,\;\lambda_{4}:0.07610\%\) \\ \hline Identified & \(\lambda_{1}=0.9981983,\;\lambda_{2}=0.9947988\) & \(\lambda_{1}:0.18017\%,\;\lambda_{2}:0.52012\%\) \\ (10 \% noise) & \(\lambda_{3}=1.005139,\;\lambda_{4}=1.002173\) & \(\lambda_{3}:0.51393\%,\;\lambda_{4}:0.21727\%\) \\ \hline \hline \end{tabular}
\end{table}
Table 11: The correct MT model and the identified one with different noises driven by the bright one-soliton solutions.
can be performed theoretically. However, it is found from the last section for forward problems that although the complicated local solutions such as two-soliton, one-breather and second-order rogue wave can be learned with aid of the XPINN algorithm, the simple solutions such as one-soliton and first-order rogue wave predicted by means of the PINN framework show better accuracy with smaller absolute errors. Hence, here we only choose three types of simple solutions including the bright, dark one-soliton and first-order rogue wave to discover unknown parameters by using the classical PINN algorithm. For three cases, we establish the network architecture consisting of 7 hidden layers with 40 neurons per layer and fix the weights \(\mathbf{W}_{R}=\mathbf{W}_{IB}=\mathbf{W}_{IN}=1\) uniformly. In these trainings, we first perform the 5000 steps Adam optimization with the default learning rate \(10^{-3}\) and then use the L-BFGS optimization with the maximum iterations 50000.

In the numerical experiments driven by three types of solutions, the unknown parameter \(\lambda_{i}\) is initialized to \(\lambda_{i}=0\), and the relative error for these unknown parameters is defined as \(\frac{|\hat{\lambda}_{i}-\lambda_{i}|}{\lambda_{i}}\), where \(\hat{\lambda}_{i}\) and \(\lambda_{i}\) denote the predicted and true value respectively for \(i=1,2,3,4\). The parameters and the computational domains in the bright, dark one-soliton and first-order rogue wave solutions are still taken to be same as those for the forward problems in the last section. By imposing \(400\times 600\) grid points with the equal step length in each space-time region, we obtain the discretized data sets of the initial-boundary value and the internal true solution. Furthermore, \(N_{R}=20000\) collocation points, \(N_{IB}=1000\) initial-boundary points and \(N_{IN}=2000\) internal points are generated by using the LHS method. After optimizing the learnable parameters in the NN, the unknown parameters of the MT model in the final experimental results under three types of reference solutions are shown in Tables 11, 12 and 13. From any of the three cases, it can be seen that the unknown parameters are correctly identified by the network model and the relative errors increase observably with the noise intensity. In addition, whether or not the training data is affected by noise, these parameter errors are very small. This observation is an evidence for the robustness and efficiency of the PINN approach for the MT model.

\begin{table}
\begin{tabular}{c c c} \hline \hline \multirow{2}{*}{MT model} & \multicolumn{2}{c}{Item} \\ \cline{2-3}  & Parameters & Relative errors \\ \hline Correct & \(\lambda_{1}=1,\ \lambda_{2}=1\) & \(\lambda_{1}:0\%,\ \lambda_{2}:0\%\) \\  & \(\lambda_{3}=1,\ \lambda_{4}=1\) & \(\lambda_{3}:0\%,\ \lambda_{4}:0\%\) \\ \hline Identified & \(\lambda_{1}=1.000048,\ \lambda_{2}=1.000005\) & \(\lambda_{1}:0.00475\%,\ \lambda_{2}:0.00048\%\) \\ (without noise) & \(\lambda_{3}=0.9999870,\ \lambda_{4}=0.9999554\) & \(\lambda_{3}:0.00130\%,\ \lambda_{4}:0.00446\%\) \\ \hline Identified & \(\lambda_{1}=0.9991396,\ \lambda_{2}=0.9976323\) & \(\lambda_{1}:0.08604\%,\ \lambda_{2}:0.23677\%\) \\ (5 \% noise) & \(\lambda_{3}=0.9989458,\ \lambda_{4}=0.9981177\) & \(\lambda_{3}:0.10542\%,\ \lambda_{4}:0.18823\%\) \\ \hline Identified & \(\lambda_{1}=0.9976596,\ \lambda_{2}=0.9953481\) & \(\lambda_{1}:0.23404\%,\ \lambda_{2}:0.46519\%\) \\ (10 \% noise) & \(\lambda_{3}=0.9975830,\ \lambda_{4}=0.9962552\) & \(\lambda_{3}:0.24170\%,\ \lambda_{4}:0.37448\%\) \\ \hline \hline \end{tabular}
\end{table}
Table 13: The correct MT model and the identified one with different noises driven by the first-order rogue wave solutions

\begin{table}
\begin{tabular}{c c c} \hline \hline \multirow{2}{*}{MT model} & \multicolumn{2}{c}{Item} \\ \cline{2-3}  & Parameters & Relative errors \\ \hline Correct & \(\lambda_{1}=1,\ \lambda_{2}=1\) & \(\lambda_{1}:0\%,\ \lambda_{2}:0\%\) \\  & \(\lambda_{3}=1,\ \lambda_{4}=1\) & \(\lambda_{3}:0\%,\ \lambda_{4}:0\%\) \\ \hline Identified & \(\lambda_{1}=1.000278,\ \lambda_{2}=0.9995433\) & \(\lambda_{1}:0.02784\%,\ \lambda_{2}:0.04567\%\) \\ (without noise) & \(\lambda_{3}=1.000106,\ \lambda_{4}=0.9999900\) & \(\lambda_{3}:0.01057\%,\ \lambda_{4}:0.00100\%\) \\ \hline Identified & \(\lambda_{1}=1.003280,\ \lambda_{2}=0.9984017\) & \(\lambda_{1}:0.32800\%,\ \lambda_{2}:0.15983\%\) \\ (5 \% noise) & \(\lambda_{3}=1.002683,\ \lambda_{4}=0.9964039\) & \(\lambda_{3}:0.26831\%,\ \lambda_{4}:0.35961\%\) \\ \hline Identified & \(\lambda_{1}=1.006322,\ \lambda_{2}=0.9959229\) & \(\lambda_{1}:0.63217\%,\ \lambda_{2}:0.40771\%\) \\ (10 \% noise) & \(\lambda_{3}=1.005789,\ \lambda_{4}=0.9930986\) & \(\lambda_{3}:0.57889\%,\ \lambda_{4}:0.69014\%\) \\ \hline \hline \end{tabular}
\end{table}
Table 12: The correct MT model and the identified one with different noises driven by the dark one-soliton solutions.
## 4 Conclusions and discussions

In conclusion, we investigate data-driven localized wave solutions and parameter discovery in the MT model via the PINNs and modified XPINNs algorithm. For the forward problems, abundant localized wave solutions including soliton, breather and rogue wave are learned and compared to exact ones with relative \(L^{2}\)-norm and point-wise absolute errors. In particular, we have modified the interface line in the domain decomposition of XPINNs into the small interface zone shared by two adjacent subdomains. Then three types of interface condition including the pseudo initial, residual and gradient ones are imposed in the small interface zones. This modified XPINNs algorithm is applied to learn the data-driven solutions with relatively complicated structures such as bright-bright soliton, dark-dark soliton, dark-antidark soliton, general breather, Kuznetsov-Ma breather and second-order rogue wave. Numerical experiments show that the improved version of XPINNs can achieve better stitching effect and faster convergence for predicting dynamical behaviors of localized waves in the MT model. For the inverse problems, the coefficient parameters of linear and non-linear terms in the MT model are learned with clean and noise data by using the classical PINNs algorithm. The training results reveal that for three types of data-driven solutions, these unknown parameters can be recovered accurately even for the sample data with certain noise intensity.

## Acknowledgement

The work was supported by the National Natural Science Foundation of China (Grant Nos. 12226332, 119251080 and 12375003) and the Zhejiang Province Natural Science Foundation of China (Grant No. 2022SJ-GYZC01).
