Artificial intelligence based on deep learning has promoted the significant development of computer vision [1, 2] and especially large language model [3, 4] recently with the huge increase of the GPU computing capability. Physics-informed neural networks (PINNs) model depending on the deep neural network has been utilized to solve a lot of nonlinear partial differential equations (NPDEs), such as Burgers, nonlinear Schrödinger, and Navier–Stokes equations with the initial and boundary (IB) conditions [5–7]. Moreover, PINNs have been used to solve a complex nonlinear Schrödinger equation, enhancing optical fiber modeling and sensor/communication performance [8]. PINNs with a local time-updating discrete scheme are proposed for efficient and accurate forward and inverse soil consolidation analysis, demonstrating its superiority over continuous-time schemes, especially in high-dimensional problems [9]. In the scientific computation field, extensive studies have been carried out to improve, accelerate, and extend the original PINNs software. Convergence of the numerical simulation is expedited by introducing the global and local adaptive activation functions in the neural networks [10–12]. Improved PINNs are achieved by leveraging the conservative properties of the physical model to enhance the prediction accuracy [13]. Distributed GPU computing is made use of to reduce the training time in some pioneering research [14–16]. The adaptive sampling method significantly reduces convergence time and the required number of collocation points [17, 18]. A gradient-optimized model, achieved by analyzing the distribution of the back-propagated gradients of loss with respect to the neural network parameters, can balance the different loss terms to enhance the training efficiency [19, 20]. Incorporating physical causality into the loss function extends PINNs software capabilities to simulate the solutions of the dynamical systems exhibiting multi-scale, chaotic, or turbulent behaviors [21, 22]. The Fourier neural operator by parameterizing the integral kernel directly in Fourier space, is another method with a novel deep neural network to model turbulent flows with zero-shot super-resolution [23–25]. A deep learning framework that forcibly encodes a given physics structure in a recurrent convolutional neural network can facilitate the learning of the reaction-diffusion processes with high accuracy, robustness, and generalizability [26, 27]. PINNs with a self-adaptive macro model to accurately predict mixed traffic flow using sparse data, improving traffic simulation and interpretability [28]. The Deep Ritz Method uses deep learning to solve variational problems from PDEs, offering adaptability and scalability to high dimensions [29]. Deep Neural Networks are utilized to solve PDEs in computational mechanics by leveraging energy-based formulations as natural loss functions [30]. The adaptive deep learning method is introduced to analyze melting heat transfer in non-Newtonian fluids, efficiently solving boundary-layer problems using physics-informed constraints and transfer learning [31].

Here, we focus on the soliton wave solutions of the fourth-order Boussinesq equation using the PINNs model. The application of the Boussinesq equations includes the description of nonlinear phenomena in ocean engineering [32, 33], ion sound waves in plasmas [34, 35], and solitary waves in shallow water models [36, 37]. Additionally, the Boussinesq equation is explored for its ability to accurately simulate wave propagation, breaking, and interaction with nearshore topography, including dam-break scenarios [38]. Mathematically, the exact soliton wave solutions of the Boussinesq equation have been presented using the Hirota bilinear method, a comprehensive approach for constructing exact solutions of NPDEs [39, 40]. Boussinesq equations are in the form of ${u}_{tt}+{{ \mathcal N }}_{x}[u]=0$, in which ${{ \mathcal N }}_{x}[u]$ is a nonlinear operator acting on u with respect to x. Former research handling these kinds of NPDEs do not give a specific emphasis on the second-order operator in time with PINNs model [10, 19, 41]. Qian et al suggested that new forms for the PINN training loss function, including the gradient of the equation residual and the gradient of the initial-condition residual, are crucial to the error estimate but would be absent from the canonical PINNs formulation of the loss function in recent studies [42]. More attention should be devoted to solving the Boussinesq equation, which involves high-order spatial derivatives and second-order time derivatives, using the PINNs model. PINNs leverage automatic differentiation, which allows for the efficient calculation of gradients with respect to the inputs of the neural network. In contrast, traditional methods require discretizing differential operators, potentially introducing approximation errors, especially for high-order PDEs [6, 43]. Additionally, PINNs eliminate the need for mesh generation, offering significant advantages for problems with complex geometries or moving boundaries over traditional mesh-based approaches [16, 44]. Moreover, PINNs excel in solving inverse problems, enabling the inference of unknown parameters in PDEs from sparse data, thus outperforming traditional optimization-based techniques [6, 16].

In this paper, the improved PINNs (IPINNs) model is utilized to solve the Boussinesq equation with loss function containing a new initial first-order time derivative term by the aid of adaptive activation function and loss-balanced coefficients. The adaptive activation functions in the neural networks are employed to accelerate numerical convergence, as demonstrated in previous studies [10, 12]. In the current simulation, only the global adaptive activation function is used, achieved by adding a new scalable hyperparameter that dynamically alters the topology of loss function involved in the optimization process [10]. Loss-balanced coefficients are introduced based on the former research to enhance the prediction accuracy of PINNs [42, 45, 46]. Raissi previously suggested the potential extension of the original PINNs framework to solve the Boussinesq-type equations [5]. Recent studies have provided soliton [41, 47] and rogue [48] wave solutions of the Boussinesq equation with data-driven neural network model. In our current simulations, we revisit soliton wave solvers using the IPINNs model, achieving relative L₂-norm errors around 10⁻⁴ in all the simulations. This represents a significant improvement of two orders of magnitude in specific cases compared to previous results [41] using conventional PINNs. Parameter scans for the Boussinesq equation solvers are performed to demonstrate the versatility of the IPINNs model, revealing a close correlation between the width of two-soliton solutions and the training time. Notably, the increased computational demand for solving the fourth-order equation is mitigated by automatic differentiation, allowing the model to bypass the Courant–Friedrichs–Lewy (CFL) condition encountered in traditional methods. The diverse initial and boundary (IB) conditions, generated by the solvers with varying parameters, are utilized to validate the model within the same process. In some simulations, comparisons are drawn between the IPINNs model and a variant of the PINNs model, which omits the new initial first-order time derivative term.

The rest of the paper is organized as follows. In Section 2, we describe the updated loss function, which incorporates a new loss term and loss-balanced coefficients for the IPINNs model used in the current simulations. In Section 3, the one-soliton and two-soliton wave solvers for the Boussinesq equation are presented, showing good agreement between the exact and learned results. Finally, brief conclusions and discussions are provided in Section 5.

To describe the IPINNs model, the Boussinesq equation constrained by the IB conditions is firstly brought in, where X₁, X₂ and T₁, T₂ are spatial and temporal domain boundaries, respectively. I(x), Ip(x), B(X₁, t) and B(X₂, t) represent the IB conditions. Data-driven solutions for the Boussinesq equation are then simulated using the neural network setup of the IPINNs model. For the current Boussinesq equation, the loss function is defined as follows utilizing the mean squared error method: Here Θ denotes the weight matrices and bias vectors in the neural network architecture. λ_u, λ_I, λ_Ip, and λ_B are the loss-balanced coefficients and will be learned to balance the interplay among the different loss terms [42, 45, 46]. The new initial loss term would be neglected when the hyperparameter λ_Ip = 0. The newly added hyperparameter a adjusts the slope of the activation function, accelerating the training process compared to the fixed default value of 1.0 [10]. Mathematically, the adaptive activation function replaces the x variable with ax. Thus, the TANH (hyperbolic tangent) function, for example, changes from the common form (e^x − e^−x)/(e^x + e^−x) to (e^ax − e^−ax)/(e^ax + e^−ax). $({x}_{I}^{i},{t}_{I}^{i})$, $({x}_{B}^{i},{t}_{B}^{i})$, and $({x}_{u}^{i},{t}_{u}^{i})$ represent the initial, boundary, and collocation points for u(x, t) in the spatial and temporal domains(x ∈ [X₁, X₂], t ∈ [T₁, T₂]). u and $\hat{u}$(${\hat{u}}_{t}$) represent the learned and exact values with $f={u}_{tt}({\rm{\Theta }},a,t,x)\,+{{ \mathcal N }}_{x}[u({\rm{\Theta }},a,t,x)]$ on behalf of the predicted physical model. N_u, N_I(N_Ip) and N_B are the number of the collocation points, initial and boundary points. A new generalized loss term J_Ip (second term in J_I) is added in the initial condition due to the second-order derivative of u(x, t) with respect to t. Data-driven deep neural network is determined by the collocation points in the simulation domain and the IB points. The loss term J_u to be minimized at a finite set of collocation points can learn the physical model ${u}_{tt}+{{ \mathcal N }}_{x}[u]=0$. The other two loss terms J_I and J_B are trained to satisfy the constraints of the IB points. The automatic differentiation method and backpropagation algorithm are utilized to train the deep neural network by calculating the derivative of the loss function with respect to Θ, λ_u, λ_I, λ_Ip, λ_B, and a denoting the trainable variables.

To solve the Boussinesq equation, the IPINNs neural network model is configured with TANH activation functions and initialized Xavier methods. Dirichlet boundary conditions are applied to the spatial domain edges. The Adam optimizer is utilized to achieve a better set of initial learnable variables [49]. The L-BFGS algorithm [50], with second-order derivatives of the loss function, is then used to finetune the parameters to obtain higher accuracy. The combined use of the Adam and L-BFGS algorithms significantly accelerates numerical convergence and enhances robustness compared to the exclusive use of a single optimizer. The default maximum optimizer epochs for the Adam and L-BFGS algorithms are set as 10 000 and 60 000, respectively. The L-BFGS algorithm would always finish running the software earlier automatically by the convergence mechanism in all the current simulations. The initial learning rates for Adam optimizers and L-BFGS algorithm are 0.007 and 0.9 for the simulations. The simulations are conducted on Python 3.7 and TensorFlow_GPU-1.14 software, utilizing the NVIDIA K80 GPU computing platforms with the Linux system.

The nonlinear (1 + 1) dimensional fourth-order Boussinesq equation in the current simulation is as follows: The number of collocation points N_u is set as 10 000 by default generated by the Latin hypercube sampling method which are randomly parsed from the full high-resolution dataset in the IPINNs model. The default numbers of the IB training data are N_I = N_Ip = 256 and N_B = 100 in which the latter one is selected for both the lower and upper boundary conditions. In some simulations, the default parameters are modified, and these changes will be explicitly illustrated. The IB training data and collocation points are loaded into the IPINNs model to train the neural network parameters and predict the solutions of the Boussinesq equation. The analytic one-soliton and two-soliton solutions [39, 40] for the Boussinesq equation (3) to be compared with the learned results are shown by the following two equations, respectively: where In the colliding-soliton simulations, the parameter cc = 0 is used, while for the chasing-soliton cases, cc = 1. The parameter k in the one-soliton case, and k₁, k₂ in the two-soliton case, determine the different IB conditions. These parameters influence the solitons’velocity, width, and interaction behavior, thus playing a crucial role in shaping the dynamics of the soliton solutions within the simulation framework.

The IPINNs model has been employed to study one-soliton and two-soliton wave solutions for the fourth-order Boussinesq equation, providing insights into the dynamic behavior and interaction between solitons. The loss function, accounting for the Boussinesq equation’s second-order time dependence, includes a new initial first-order time derivative term, which is crucial for achieving high prediction accuracy. Adaptive activation function mechanism and loss-balanced coefficients are utilized to accelerate and enhance the PINNs model, as highlighted in previous research [10, 12, 42, 45]. The IPINNs model achieves high prediction accuracy, with relative L₂-norm errors around 10⁻⁴. This represents a significant improvement, with prediction precision increased by two orders of magnitude compared to previous results [41]. Comparisons among the TANH, SIN, GELU, SWISH, and MISH activation functions in the neural network are conducted using the colliding-soliton k = 1.1 case, confirming that the selection of TANH activation function is efficient for the current simulations. To validate the versatility and robustness of the IPINNs model, parameter scans are conducted for the two-soliton wave solvers of the Boussinesq equation. The training time is found to be closely related to the width of the two-soliton solutions when using the L-BFGS algorithm during colliding-soliton simulations. For instance, the width of the soliton with k = 1.7 is approximately 3.5 times narrower than that of the k = 0.7 case, with a corresponding training time that is about 11 times longer. However, the increased training time is not as pronounced as in traditional mesh-based algorithms, thanks to the use of automatic differentiation for the fourth-order Boussinesq equation. In the chasing-soliton scenario, when the condition k₂ = k₁ = 1.5 is met, the two solitons merge into one, resulting in a significant reduction in training time. The current model requires more initial data, which presents a challenge in obtaining the additional information.

The widespread applications of the PINNs software in solving NPDEs still require further understanding and optimization, as its prediction abilities are unstable in certain scenarios. To enhance the generalizability of deep neural network models, attention should be directed towards numerous NPDEs in physical fields that have yet to be solved using PINNs software. For instance, the modified Korteweg-de Vries–sine-Gordon equation, which describes few-cycle-pulse optical solitons and features rich structures such as soliton molecules, breather molecules, and breather-soliton molecules, presents a challenging problem for the PINNs model. Laser physics phenomena associated with the Boussinesq model will be studied in the future using the current model to assess its capability [51–53]. The soliton and lump solutions of the (2 + 1)-dimensional (3D) generalized variable-coefficient Hirota-Satsuma-Ito equation are worthy of study with the current model to verify its robustness and efficiency [54]. The computational efficiency of the current model when applied to 3D PDEs is a potential challenge and requires further study. For 3D PDEs with solutions exhibiting locally steep gradients or high-frequency components, adaptive sampling methods [16] and self-adaptive weight algorithms with scalable multi-GPU support [55] emerge as promising candidate approaches. The SAEE-PINN model using sparse data also provides valuable insights for our future research [28]. For PDEs with high-order or irregular domains, PINNs with a mesh-free method offer advantages over traditional methods that rely on mesh grid techniques [16, 56]. Future research will compare the current multilayer perceptron to the Kolmogorov–Arnold representation, which is also receiving attention in PINNs research. Combining domain decomposition techniques [15, 57, 58] with the current IPINNs model offers a promising approach for conducting large-domain simulations with higher prediction accuracy. Further advancements in heterogeneous parallel computing would significantly accelerate training through distributed data-parallel methods, which adjust the mini-batch sizes across workers. Incorporating more accelerated computing platforms in software development will also enhance the versatility and efficiency of software execution. Upgrading from TensorFlow 1 to TensorFlow 2 will be completed in our further work, improving support for distributed computing, automatic mixed precision, and other optimizations.

The authors declare that they have no conflict of interest.

This research does not involve human participants and/or animals.