For soliton solutions of discrete NLEs with $M$ solitons, we first construct $M$ traveling wave invariants to characterize the translational invariance of each soliton: where $v_m$ denotes the exact physical velocity of the $m$-th soliton (fixed as analytical values derived from soliton solution formulations), $n$ is the discrete spatial variable, and $t$ is the temporal variable.

For any spacetime point $(n,t)$, we define $M+1$ prediction components $\{u_0, u_1, \dots, u_M\}$ via a shared neural network:

We denote the parameters in the neural network, namely the weights and biases, by $\theta$. In the design of the loss function, to represent the residual of the discrete mKdV equation, we perform automatic differentiation on the time dimension of the discrete mKdV equation, and handle the spatial dimension using forward and backward difference methods. That is, for its spatial domain, we use $n \pm i$ to represent $u_{n \pm i}$, and then obtain the residual $\xi(u, n \pm i, \theta)$. The specific expression of the loss function is as follows: with where $w_1$, $w_2$, $w_3$, $w_4$, and $w_5$ are weights for the five loss components, corresponding to the physics equation loss, initial condition loss, boundary condition loss, symmetry loss, and translation loss, respectively. $\mathcal{L}_\mathrm{F}(\theta)$ denotes the physics equation loss, $N_\mathrm{F}$ is the number of sampling points for computing PDE residuals, $\xi(u, n\pm i, \theta)$ represents the nonlinear terms of NLEs; $\mathcal{L}_\mathrm{I}(\theta)$ denotes the initial condition loss, $N_\mathrm{I}$ is the number of sampling points for the initial condition, $n_i$ is the spatial sampling point at the initial time; $\mathcal{L}_\mathrm{B}(\theta)$ denotes the boundary condition loss, which constrains the neural network output in the boundary region to match the given boundary condition, $N_\mathrm{B}$ is the number of sampling points at the boundary, when the boundary condition is $u(t, n_\mathrm{B}) = u_\mathrm{B}(t)$ ($n_\mathrm{B}$ is the point on the boundary), $(t_i, n_{\mathrm{B},i})$ is the spatiotemporal sampling point on the boundary; $\mathcal{L}_\mathrm{S}(\theta)$ denotes the symmetry loss, which constrains the neural network output to satisfy the central symmetry property. $N_\mathrm{S}$ is the number of sampling points for computing the symmetry loss; $(t_i, n_i)$ is a selected spatiotemporal sampling point, and $(2c_t - t_i, 2c_n - n_i)$ is its symmetric point with respect to the learnable symmetry center $(c_t, c_n)$. $\mathcal{L}_\mathrm{T}(\theta)$ denotes the translation loss, $N_\mathrm{T}$ is the number of sampling points for computing the translation loss, $u_0$ denotes the baseline value of the neural network prediction, $u_m$ ($m = 1,2,\dots,M$) represents the $m$th component of the neural network prediction at the sampling points, and the loss term constrains the difference between the baseline value $u_0$ and the average of $u_0$ and the sum of $M$ prediction components $\sum_{m = 1}^M u_m$ at each of the $N_\mathrm{T}$ sampling points.

The complete architecture is shown in figure 1.

This paper adopts the Hammersley sampling method [39] to generate the sampling training set $ N_\mathrm{F} $. First, the Adam optimizer is used for training, and then the Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimizer is used for further refined training. In equation (7), we assign different weights to each loss term, and different weight configurations enable the neural network to correctly learn the physical laws of the equation.

In the neural network architecture adopted in this paper, it first includes an input layer for processing two-dimensional spatiotemporal data. For different soliton solutions, the number of hidden layers used is 4, and the number of neurons corresponding to each layer is 32/64. The activation function used is the hyperbolic tangent (tanh) function. In terms of optimizer settings, the learning rate of Adam is set to $ 1 \times 10^{-3} $, while the learning rate of L-BFGS is set to 1.

To effectively evaluate the performance of the model, we use the relative error as an indicator to evaluate the model’s effect. Its function is defined as follows:

The code runs in a Python 3.11.7 environment using PyTorch version 2.7.0+cu126 on a high-performance computer equipped with an NVIDIA GeForce RTX 4060 laptop GPU.

Based on the above descriptions, the core concepts of TC-PINNs are summarized as follows, and their implementation details are summarized in table 1:

First, we study one-soliton solution of equation (1). We use one-soliton solution in [35] as follows: where $\xi_1 = \left(2\ln\lambda_1\right)n + \left(\lambda_1^4 - \frac{1}{\lambda_1^4} + \frac{2}{\lambda_1^2} - 2\lambda_1^2\right)t$, where $\lambda_1 \gt 0$ is the spectral parameter with a value of $\frac{1}{2}$, the soliton velocity corresponding to this one-soliton solution is given by $v = \frac{\lambda_1^8 - 2\lambda_1^6 + 2\lambda_1^2 - 1}{2\lambda_1^4 \ln \lambda_1}$. The training domain is $n \in [-5, 5]$ and $t \in [-1, 1]$. The Hammersley sampling method is adopted, with $N_\mathrm{F} = 2000$ physical sampling points and $N_\mathrm{B} = 200$ boundary sampling points. The initial time $t_0$ is set to $0.0$, and the number of initial points $N_\mathrm{I}$ is $100$. First, we conduct experiments using PINNs without soft/hard constraints, where the loss function weights are $w_1 = 10$, $w_2 = 100$, $w_3 = 100$. The training process is divided into two stages: the first stage uses the Adam optimizer for $2000$ iterations, and the second stage uses the L-BFGS optimizer for $50$ iterations, and the $\text{RE}_{\text{Ins}}$ is $9.6054 \times 10^{-4}$. Next, we apply soft and hard constraints based on the soliton symmetry property, respectively. For the soft constraint, a symmetry loss term is added to the loss function with $w_4 = 10$, and the $\text{RE}_{\text{Ins}}$ is $7.5113 \times 10^{-4}$. For the hard constraint, the symmetry property is embedded into the neural network without adding a symmetry loss term, and the $\text{RE}_{\text{Ins}}$ is $5.5394 \times 10^{-4}$. Finally, we apply soft and hard constraints based on the soliton translational property, respectively. For the soft constraint, a translational loss term is added to the loss function with $w_5 = 10$, and the $\text{RE}_{\text{Ins}}$ is $2.9499 \times 10^{-4}$. For the hard constraint, the translational property is embedded into the neural network without adding a translational loss term, and the $\text{RE}_{\text{Ins}}$ is $1.83885 \times 10^{-4}$. The spatiotemporal heatmaps under the five methods are compared in figure 2.

To enhance the credibility of the results, each method is tested with five different random seeds, and the relative errors are presented in the form of mean $\pm$ standard deviation. We also perform extrapolation to compare the generalization ability of the models, with the extrapolation range set to $n \in [-5, 5]$ and $t \in [-1.2, -1] \cup [1, 1.2]$. The experimental results are also stored in table 2. The results show that after applying soft or hard constraints, the model can learn the symmetry or translation properties of solitons and improve the accuracy of numerical simulations. Moreover, the improvement in model accuracy is more significant under hard constraints, and the extrapolation error accuracy of hard constraints is higher than that of soft constraints, indicating that hard constraints have superior generalization ability. Furthermore, under the same conditions, the prediction performance of TC-PINNs is more excellent. In terms of computational efficiency, the proposed h-TC-PINNs achieve the best trade-off between accuracy and training speed, requiring the least training time while maintaining the lowest errors. We plot the above experimental data into a boxplot for a more intuitive visualization, as shown in figure 3.

To observe the spatiotemporal evolution of the one-soliton solution, we present the evolution diagrams obtained by h-TC-PINNs at three time instants. As can be seen from figure 4, the predicted solutions almost coincide with the true solutions, indicating that the model has successfully learned the nonlinear characteristics of the equation.

Subsequently, we compare the soft and hard constraints using the translation or symmetry properties of the two-soliton solution. We use the two-soliton solution from [35]:

where Subsequently, we take $\lambda_1 = \frac{5}{2}$, $\lambda_2 = \frac{1}{2}$ and perform numerical simulations on its two bell-shaped soliton solution (11). The training domain is the same as that for the single soliton, with $N_\mathrm{F} = 3000$ Hammersley sampling points generated and $N_\mathrm{B} = 200$ boundary sampling points. The initial time $t_0$ is set to $0.0$ and $N_\mathrm{I} = 100$ initial points are used. The soliton velocities corresponding to this two-soliton solution are given by $v_1 = \frac{\lambda_1^8 - 2\lambda_1^6 + 2\lambda_1^2 - 1}{2\lambda_1^4 \ln \lambda_1}$ and $v_2 = \frac{\lambda_2^8 - 2\lambda_2^6 + 2\lambda_2^2 - 1}{2\lambda_2^4 \ln \lambda_2}$, respectively. First, we conduct experiments without any soft or hard constraints. In the loss function, $w_1 = 10$, $w_2 = 100$, $w_3 = 100$. The model is first trained with the Adam optimizer for 2000 iterations and then with the L-BFGS optimizer for 50 iterations. After training, the $\text{RE}_{\text{Ins}}$ is $1.7925 \times 10^{-3}$. Next, we apply soft and hard constraints respectively using the soliton symmetry property. For the soft constraint, a symmetry loss term is added to the loss function with $w_4 = 10$, and the $\text{RE}_{\text{Ins}}$ is $9.1458 \times 10^{-4}$. For the hard constraint, instead of adding a symmetry loss term, the symmetry property is embedded into the neural network, and the $\text{RE}_{\text{Ins}}$ is $8.8923 \times 10^{-4}$. Finally, we apply soft and hard constraints respectively using the soliton translation property. For the soft constraint, a translation loss term is added to the loss function with $w_5 = 10$, and the $\text{RE}_{\text{Ins}}$ is $7.3507 \times 10^{-4}$. For the hard constraint, the translational property is embedded into the neural network without adding a translational loss term, and the $\text{RE}_{\text{Ins}}$ is $5.0489 \times 10^{-4}$. The spatiotemporal heatmaps under the five methods are compared in figure 5.

To enhance the credibility of the results, each method is tested with five different random seeds, and the relative errors are presented in the form of mean $\pm$ standard deviation. We also perform extrapolation to compare the generalization ability of the models, with the extrapolation range set to $n \in [-5, 5]$ and $t \in [-1.2, -1] \cup [1, 1.2]$. The experimental results are also stored in table 3. The results show that after applying soft or hard constraints, the model can learn the symmetry or translation properties of solitons and improve the accuracy of numerical simulations. Moreover, the improvement in model accuracy is more significant under hard constraints, and the extrapolation error accuracy of hard constraints is higher than that of soft constraints, indicating that hard constraints have superior generalization ability. In terms of computational efficiency, although the TC-PINNs require slightly longer training time than SDE-PINNs, they deliver substantially better prediction accuracy and generalization performance, representing a favorable trade-off between computational cost and solution quality for two-soliton problems. Furthermore, under the same conditions, the prediction performance of TC-PINNs is more excellent. We plot the above experimental data into a boxplot for a more intuitive visualization, as shown in figure 6.

Finally, we further extend the proposed method to RS solutions of the discrete mKdV equation. We take the first-order RS solution in [35] as an example:

For the numerical simulation of first-order RS solution, the spatiotemporal domain for training is the same as that in section 3.1, and the number of sampling points and neural network settings are also the same as those in section 3.1, the soliton velocity corresponding to this first-order RS solution is given by $v = -\frac{675}{64}$. First, we conduct experiments using PINNs without any soft or hard constraints. In this case, the loss function has $w_1 = 10$, $w_2 = 100$, $w_3 = 100$. The training process is divided into two stages: the first stage uses the Adam optimizer for $2000$ iterations, and the second stage uses the L-BFGS optimizer for $50$ iterations. The final $\text{RE}_{\text{Ins}}$ is $3.2045 \times 10^{-3}$. Next, we apply soft and hard constraints respectively using the soliton symmetry property. For the soft constraint, a symmetry loss term is added to the loss function with $w_4 = 10$, and the $\text{RE}_{\text{Ins}}$ is $2.3498 \times 10^{-3}$. For the hard constraint, instead of adding a symmetry loss term, the symmetry property is embedded into the neural network, and the $\text{RE}_{\text{Ins}}$ is $1.0941 \times 10^{-3}$. Finally, we apply soft and hard constraints based on the soliton translational property, respectively. For the soft constraint, a translational loss term is added to the loss function with $w_5 = 10$, resulting in a $\text{RE}_{\text{Ins}}$ of $5.3870 \times 10^{-4}$. For the hard constraint, the translational property is embedded into the neural network without adding a translational loss term, and the $\text{RE}_{\text{Ins}}$ is $3.7200 \times 10^{-4}$. The spatiotemporal heatmaps under the five methods are compared in figure 7.

To enhance the credibility of the results, each method is tested with five different random seeds, and the relative errors are presented in the form of mean $\pm$ standard deviation. We also perform extrapolation to compare the generalization ability of the models, with the extrapolation range set to $n \in [-5, 5]$ and $t \in [-1.2, -1] \cup [1, 1.2]$. The experimental results are also stored in table 4. The results show that after applying soft or hard constraints, the model can learn the symmetry or translation properties of solitons and improve the accuracy of numerical simulations. Moreover, the improvement in model accuracy is more significant under hard constraints, and the extrapolation error accuracy of hard constraints is higher than that of soft constraints, indicating that hard constraints have superior generalization ability. In terms of computational efficiency, the proposed h-TC-PINNs achieve the best trade-off between accuracy and training speed, requiring the least training time while maintaining the lowest errors. Furthermore, under the same conditions, the prediction performance of TC-PINNs is more excellent. We plot the above experimental data into a boxplot for a more intuitive visualization, as shown in figure 8.

To explore the effect of activation functions on model performance, we conduct comparative experiments using five typical activation functions, namely tanh, sin, ReLU, SiLU and Sigmoid. The quantitative results of inside predicting and extrapolation relative errors are presented in table 8.

It can be observed that the model with ReLU activation function exhibits the worst performance, with both errors reaching the order of $10^{-1}$, indicating severe instability during training. In contrast, smooth activation functions including tanh, sin, SiLU and Sigmoid achieve superior accuracy. In this work, we adopt the commonly used tanh activation function.

On the basis of determining the activation function, we further investigate the influence of network depth and width on model performance. Experiments are conducted on networks with different depths (2, 4, 6 layers) and widths (16, 32, 64 neurons per layer), and the corresponding errors are listed in table 9.

The results show that all configurations achieve comparable performance with errors at the same order of magnitude. For h-TC-PINNs, $\text{RE}_{\text{Ins}}$ range from $1.1468\times 10^{-4}$ to $2.2858\times 10^{-4}$, and $\text{RE}_{\text{Ext}}$ range from $1.8964\times 10^{-3}$ to $4.9042\times 10^{-3}$, indicating that variations in network depth and width do not lead to significant performance differences. Meanwhile, h-TC-PINNs consistently outperform standard PINNs across all architectures, with both $\text{RE}_{\text{Ins}}$ and $\text{RE}_{\text{Ext}}$ being lower than those of standard PINNs.

Considering both computational efficiency and model stability, the 4-layer network with 32 neurons per layer is chosen as the default structure for its optimal trade-off. While deeper or wider networks can yield slightly lower errors in specific cases, they introduce higher computational overhead without substantial performance gains. The selected configuration achieves balanced performance ($\text{RE}_{\text{Ins}} = 2.3491\times 10^{-4}$, $\text{RE}_{\text{Ext}} = 4.9491\times 10^{-3}$) while maintaining a reasonable computational cost. Meanwhile, the network with 64 neurons per layer can be adopted as an alternative when dealing with more complex soliton solutions that require stronger feature fitting capability.

After determining the activation function and network structure, an ablation study on loss weights is conducted to evaluate their effects on the model’s inside predicting and extrapolation performance, as illustrated in figures 9 and 10. The boxplots clearly present the median errors and distribution dispersion of h-TC-PINNs and standard PINNs under different weight combinations.

Regarding the individual influence of each weight, when the PDE residual weight $w_1$ increases from 1 to 10, the median values of both inside predicting and extrapolation relative errors decrease significantly for both models. When $w_1$ is further raised to 50, the box widths become noticeably wider, indicating an increase in error dispersion. For the initial condition weight $w_2$, the configuration $w_2 = 10$ achieves lower median error and more concentrated distribution, while $w_2 = 1$ and $w_2 = 50$ lead to larger error fluctuations. As the boundary condition weight $w_3$ grows from 10 to 100, the median extrapolation errors of both models decline gradually with more concentrated error distributions.

By comparing the two models, h-TC-PINNs exhibits narrower box widths and more stable error distributions than standard PINNs under most weight settings. In contrast, standard PINNs suffers from more significant error fluctuations especially under $w_1 = 50$ and $w_3 = 10$. Based on all observations, the combination $(w_1,w_2,w_3) = (10,10,100)$ achieves the lowest median error and the narrowest box width for both models, which is therefore adopted as the optimal loss weight combination in all our experiments.

To explore the sensitivity of the proposed model to the time step $\Delta t$, one hundred time step samples are uniformly sampled within the range of $[0.01, 0.1]$. For each $\Delta t$, the inside predicting and extrapolation relative errors are recorded. A series of statistical tests are conducted to quantitatively determine whether there exists a significant correlation between the magnitude of the time step and the model prediction errors. These methods include Pearson correlation analysis, Spearman correlation analysis, and the coefficient of determination $R^2$. The distribution of prediction errors with respect to different time steps is visualized in figure 11, which shows no apparent trend between $\Delta t$ and model errors.

To comprehensively evaluate the relationship between time step and prediction errors, each statistical metric is used for a clear purpose, with specific judgment standards as follows: specifically, the Pearson correlation coefficient is used to check the linear correlation between time step and prediction errors—a value close to 0 means no obvious linear correlation, while values close to 1 or $-1$ mean a strong positive or negative linear correlation. The Spearman correlation coefficient is used to check the monotonic correlation between them—a value close to 0 means no obvious monotonic correlation, and values close to 1 or $-1$ mean a strong monotonic correlation. The corresponding $p$-value is used to judge significance: if $p \lt 0.05$, the correlation is statistically significant; if $p \ge 0.05$, there is no statistically significant correlation. In addition, the coefficient of determination $R^2$ is used to measure how much the change of prediction errors can be explained by the time step—a value close to 0 means the time step has almost no explanatory power for error changes, while a value close to 1 means the time step can well explain the error changes.

As presented in table 10, the Pearson and Spearman correlation coefficients between the time step and both prediction errors are all close to zero. All corresponding $p$-values are far greater than the significance level of $0.05$, indicating that no statistically significant linear or monotonic correlation exists between $\Delta t$ and model accuracy. Furthermore, the $R^2$ values are nearly zero, which confirms that the variation in prediction errors cannot be attributed to the change in the time step. These statistical results are consistent with the visual distribution in figure 11, where error points scatter randomly without any clear increasing or decreasing tendency as $\Delta t$ changes.

In summary, statistical analyses demonstrate that changing the time step within the tested interval imposes no significant impact on model performance. All our experiments use the minimum time step $\Delta t = 0.01$.

To verify the rationality of the traveling wave coordinate transformation in h-TC-PINNs, three comparative schemes are tested under the same network structure and training settings: the original PINNs, PINNs with fixed $v_m$ as input, and the proposed h-TC-PINNs with learnable $v_m$.

Table 11 shows that simply adding $v_m$ as a fixed input variable brings limited improvement to the model performance. Compared with the original PINNs, PINNs+$v_m$ achieves a slight reduction in $\text{RE}_{\text{Ins}}$, but its extrapolation performance does not show obvious improvement and even lags slightly behind, indicating that fixing $v_m$ as an input fails to effectively utilize the traveling wave physical priors. In sharp contrast, the proposed h-TC-PINNs with learnable $v_m$ traveling wave transformation achieves substantial improvements in both inside predicting and extrapolation accuracy, yielding much lower errors than the other two schemes. This fully demonstrates that the learnable traveling wave transformation can effectively embed inherent physical priors into the model, enhancing its fitting accuracy and generalization ability, and further verifying the rationality and effectiveness of the adopted traveling wave transformation.

To validate the model’s applicability under unknown $v_m$, 50 initial velocities $v_0$ are sampled within $[-40.0, 10.0]$, and five random seeds are tested for each $v_0$. The results are presented in figure 12.

As can be seen, the model achieves the lowest errors near the true velocity (marked by the red dashed line). Specifically, within the range $[-20, 0]$ that accounts for 40% of the entire sampling interval, both $\text{RE}_{\text{Ins}}$ and $\text{RE}_{\text{Ext}}$ remain consistently low, showing excellent performance. As the deviation from the true velocity increases, the error gradually rises, and individual random seeds show instability. Nevertheless, the model maintains stable and low errors for most cases, demonstrating strong robustness to unknown $v_m$ and effectiveness without precise physical prior information.