A hierarchical age-structured model of optimal vermin management by contraception

1.   Introduction

In recent years, fractional partial differential equations (FPDEs) have been widely used in natural science and engineering technology ^[1,2,3,4]. The advantage of FPDEs lies in their ability to better describe materials and processes that exhibit memory and genetic properties ^[5,6]. However, the solutions of FPDEs are much more complex. Many researchers have exploited diverse techniques for the investigation of FPDEs such as the finite difference method (FDM) ^[7], finite element method ^[8], spectral method ^[9], virtual element method ^[10], etc. The development of effective numerical methods to approximate FPDEs has been the goal of some researchers.

In recent years, neural networks (NNs) have been successfully applied to solve problems in various fields ^[11,12,13]. Due to the high expressiveness of NNs in functional approximation ^[14,15,16], using NNs to solve differential and integral equations has become an active and important research field. Physics-informed neural networks (PINNs) ^{[17,18,19,20]} are machine learning models that combine deep learning with physical knowledge. PINNs embed PDEs into the loss function of the NNs, enabling the NNs to learn solutions to PDEs. The PINNs algorithm is meshless and simple, and can be applied to various types of PDEs, including integral differential equations, FPDEs, and random partial differential equations. Moreover, PINNs solved the inverse problem of PDEs just as easily as they solved the forward problem ^[17]. PINNs have been successfully applied to solve various problems in scientific computing ^[21,22,23]. Pang et al. ^[24] used the FDM to approximate the fractional derivatives that cannot be automatically differentiated, thus extending the PINNs to fPINNs for solving FPDEs.

Despite the success of deep learning in the past, solving a wide range of PDEs is theoretically and practically challenging as complexity increases. Therefore, many aspects of PINNs need to be further improved to achieve more accurate predictions, higher computational efficiency, and robustness of training. Lu et al. ^[25] proposed DeepXDE, a deep learning library for solving PDEs, introduced a new residual-based adaptive refinement method to improve the training efficiency of PINNs, and new residual points were added at the position where the residuals of the PDEs were large, so that the discontinuities of PDEs could be captured well. Zhang et al. ^[26] combined fPINNs with the spectral method to solve the time-fractional phase field models. It had the characteristics of reducing the approximate number of discrete fractional operators, thus improving the training efficiency and obtaining higher error accuracy. Wu et al. ^[27] conducted a comprehensive study on two types of sampling of PINNs, including non-adaptive uniform sampling and adaptive non-uniform sampling, and the research results could also be used as a practical guide for selecting sampling methods. Zhang et al. ^[28] removed the soft constraints of PDEs in the loss function, and used the Lie symmetry group to generate the labeled data of PDEs to build a supervised learning model, thus effectively predicting the large amplitude and high frequency solutions of the Klein-Gordon equation. Zhang et al. ^[29] introduced the symmetry-enhanced physics-informed neural network (SPINN), which incorporated the invariant surface conditions derived from Lie symmetries or non-classical symmetries of PDEs into the loss function of PINNs, aiming to improve accuracy of PINNs. Lu et al. ^[30] and Xie et al. ^[31] introduced gradient-enhanced physics-informed neural networks (gPINNs) to solve PDEs and the idea of embedding the gradient information from the residuals of PDEs into the loss functions has also proven to be effective in other methods such as Gaussian process regression ^[32].

In this paper, inspired by the above works, gfPINNs are applied to solve the forward and inverse problems of the multiterm time-fractional Burger-type equation. The integer order derivatives are handled using the automatic differentiation capability of the NNs, while the fractional derivatives of the equation are approximated using finite difference discretization ^[33,34]. Subsequently, the residual information of the equation is then incorporated into the loss function of NNs and optimized to yield optimal parameters. For the inverse problems of the multiterm time-fractional Burger-type equation, their overall form are known but the coefficient and the orders of time-fractional derivatives are unknown. The gfPINNs explicitly incorporate information from the equation by including the differential operators of the equation directly into the optimization loss function. The parameters to be identified appear in the differential operators, which are then optimized by minimizing the loss function associated with those parameters. A numerical comparison between fPINNs and gfPINNs is conducted using numerical examples. The numerical results demonstrate the effectiveness of gfPINNs in solving the multiterm time-fractional Burger-type equation.

The structure of this paper is as follows. In Section 2, we define forward and inverse problems for the multiterm time-fractional Burger-type equation. In Section 3, we introduce fPINNs and gfPINNs and give the finite difference discretization to approximate the time-fractional derivatives. In Section 4, we demonstrate the effectiveness of gfPINNs in solving the forward and inverse problems of the multiterm time-fractional Burger-type equation by numerical examples, and compare the experimental results of fPINNs and gfPINNs. Finally, we give the conclusions of this paper in Section 5.

2.   Multiterm time-fractional Burger-type equation

We consider the following multiterm time-fractional Burger-type equation defined on the bounded domain $\Omega$:

where $(x, t)\in\Omega\times[0, T]$ and the initial and boundary conditions are given as

where $u(x, t)$ is the solution of the equation, $f(x, t)$ is the forcing term whose values are only known at scattered spatio-temporal coordinates, $v$ is the kinematic viscosity of fluid, $g(x)$ is a sufficiently smooth function, the fractional orders $\alpha$ and $\gamma$ have been restricted to (0, 1) and (1, 2), respectively, $_{0}^{C}D^{\theta}_{t}u(x, t)$ is the Caputo time-fractional derivative of order $\theta$ $(\theta > 0, n-1\le \theta < n)$ of $u(x, t)$ with respect to $t$ ^[35,36]:

where $\Gamma(\cdot)$ is the gamma function.

The forward and inverse problems of solving the multiterm time-fractional Burger-type equation are described as follows. For the forward problem, under the given preconditions of the fractional orders $\alpha$ and $\gamma$, the forcing term $f$, and the initial and boundary conditions, the solution $u(x, t)$ is solved. For the inverse problem, under the given preconditions of the initial and boundary conditions, the forcing term $f$, and additional concentration measurements at the final time $u(x, t) = h(x, t)$, the fractional orders $\alpha$ and $\gamma$, the flow velocity $v$, and the solution $u(x, t)$ are solved.

3.   Methodology

3.1.   fPINNs

This subsection introduces the idea of fPINNs and we consider both the forward and inverse problems, along with their corresponding NNs. We first consider the forward problem of the multiterm time-fractional Burger-type equation in the following form:

where $\mathscr{L}\{\cdot\}$ is a nonlinear operator and $\mathscr{L}\{u(x, t)\} = c_{1}{}_{0}^{C}D^{\alpha}_{t}u(x, t)+c_{2}{}_{0}^{C}D^{\gamma}_{t}u(x, t)+u(x, t)\frac{\partial u(x, t)}{\partial x}-v\frac{\partial^{2} u(x, t)}{\partial x^{2}}$. We divide the nonlinear operator $\mathscr{L}\{\cdot\}$ into two parts, $\mathscr{L} = \mathscr{L}_{AD} +\mathscr{L}_{nonAD}$. The first part is an integer derivative operator, which can be automatically differentiated (AD) using the chain rule. We have

and the second category consists of operators that lack automatic differentiation capabilities:

For $\mathscr{L}_{nonAD}$, we can discretize it using FDM and denote by $\mathscr{L}_{FDM}$ the discretization version of $\mathscr{L}_{nonAD}$.

During the NNs training process, our goal is to optimize its parameters in order to ensure that the approximate solution of the equation closely satisfies the initial and boundary conditions. The approximate solution is chosen as

where $u_{NN}$ represents the output of the NNs. The NNs acts as a surrogate model, approximating the relationship between spatio-temporal coordinates and the solution of the equation. It is defined by its weights and biases, forming the parameter vector $\mu$; see Figure 1 for a simple NN. This is fully connected with a single hidden layer consisting of three neurons. In this network, $x$ and $t$ are two inputs, which go through a linear transformation to obtain $x_{1} = w_{1}x+w_{4}t+b_{1}$, $x_{2} = w_{2}x+w_{5}t+b_{2}$, and $x_{3} = w_{3}x+w_{6}t+b_{3}$ in the hidden layer, and then, they go through a nonlinear transformation to get $Y_{i} = f(x_{i})$ for $i = 1, 2, 3$. We choose the hyperbolic tangent function $\tanh(\cdot)$. $Y_{i}$ to go through a linear transformation to obtain the output of the NNs, $u_{NN}(x, t;\mu) = w_{7}Y_{1}+w_{8}Y_{2}+w_{9}Y_{3}+b_{4}$. The vector of parameters $\mu$ is comprised of the weights $w_{i}$ and biases $b_{i}$. $\rho(0) = \rho(1) = 0$ and the auxiliary function $\rho(x)$ is preselected. $g(x)$ is the initial condition function such that it satisfies the initial and boundary conditions automatically.

The loss function of fPINNs for the forward problem with the approximate solution is defined as the mean-squared error of the equation residual

where $S_{F}\subset\Omega \times [0, T]$ and $\left | S_{F} \right |$ represents the number of training points. Then, we train the NNs to optimize the loss function of the forward problem with respect to the NNs parameters $\mu$, thus obtaining the optimal parameters $\mu_{best}$. Finally, we specify a set of arbitrary test points to test the trained NNs and observe the training performance.

The codes for solving the forward and inverse problems of the equation using NNs is similar. We only need to incorporate the parameters to be identified in the inverse problem into the loss function to be optimized in the forward problem, and no other changes are necessary. Next, we consider the following form of the inverse problem:

where $\xi$ is the parameter of the equation, so the loss function $L_{IV}$ for the inverse problem under consideration is

where $\alpha\in (0, 1)$ and $\gamma\in (1, 2)$, $S_{I_{1}}\subset\Omega \times [0, T]$ and $S_{I_{2}} \subset\Omega \times [0, T]$ are two sets of different training points, and $W_{I_{1}}$ and $W_{I_{2}}$ are preselected weight coefficients. We train the NNs to minimize the loss function, thereby obtaining $\alpha_{best}$ and $\gamma_{best}$, the flow velocity $v_{best}$, and the optimal parameters $\mu_{best}$ of the NNs.

3.2.   gfPINNs

We incorporate the residual information of the equation into the loss function of NNs and train the NNs to minimize this loss function, thus obtaining the optimal parameters of NNs. If the residuals in the PDEs are zero, then the gradient of the residuals in the PDEs should also be zero. Therefore, adding gradient information to the loss function is a necessary condition for training NNs. One motivation behind gfPINNs is that the residual in the loss function often fluctuates near zero. Penalizing the slope of the residual can reduce these fluctuations, making the residual closer to zero. In this section, we continue to consider the formulation of the forward and inverse problems of the equation discussed in the previous section.

We first consider the forward problem in the form of (3.1) and provide the loss function of gfPINNs for this form:

where

and the approximate solution of the equation is the same as Eq (3.4): $\tilde{u}(x, t) = \rho(x)u_{NN}(x, t)+g(x)$. The expression $L_{FW}$ as shown in Eq (3.5), where $W_{F}$, $W_{g_{1}F}$, and $W_{g_{2}F}$ are preselected weighting coefficients, $S_{g_{1}F}\subset\Omega \times [0, T]$ and $S_{g_{2}F}\subset\Omega \times [0, T]$ are two sets of different training points.

Next, we consider the inverse problem in the form of (3.6) and provide the loss function of gfPINNs for this form. The approach for the inverse problem of gfPINNs is similar to that of fPINNs. We provide the loss function for the inverse problem of gfPINNs.

where

and the expression $L_{IV}\{\mu, \xi = \{\alpha, \gamma, v\}\}$ as shown in Eq (3.7), where $W_{I}$, $W_{g_{1}I}$, $W_{g_{2}I}$, $W_{g_{1}I_{1}}$, $W_{g_{1}I_{2}}$, $W_{g_{2}I_{1}}$, and $W_{g_{2}I_{2}}$ are preselected weighting coefficients, $S_{g_{1}I_{1}}, S_{g_{2}I_{1}}\subset\Omega \times [0, T]$, $S_{g_{1}I_{2}}, and S_{g_{2}I_{2}} \subset\Omega \times [0, T]$ are four sets of different training points.

This defines the loss function of gfPINNs, which is exactly the same as discussed above for fPINNs. We train the NNs to obtain the optimal parameters of the NNs.

3.3.   Finite difference discretization for time-fractional derivatives

In the $x$ direction $[0, M]$, we take the mesh points $x_{p} = ih_{x}, i = 0, 1, 2, ..., M_{1}$, and in the $t$ direction $[0, T]$, we take the mesh points $t_{n} = n\tau, n = 0, 1, ..., N$, where $h_{x} = \frac{M}{M_{1}}$ and $\tau = \frac{T}{N}$ are the uniform spatial step size and temporal step size, respectively. Denote $\Omega _{h}\equiv \left \{ 0\le i \le M_{1} \right \}$, $\Omega _{\tau}\equiv \left \{ 0\le n \le N \right \}$. Suppose $u_{i}^{n}$ = $u(x_{i}, t_{n})$ is a grid function on $\Omega_{h} \times \Omega_{\tau}$.

We approximate the fractional derivatives of the equation using the finite difference discretization ^[33,34].

For $\alpha \in (0, 1)$, we have ${}_{0}^{C}D_{t}^{\alpha}u(x, t)\mid _{(x_{i}, t_{n})} = D_{\tau}^{\alpha} \tilde{u} _{i}^{n}+R_{1}(\tilde{u} _{i}^{n})$,

where $\tilde{u} _{i}^{n} = \tilde{u}(x_{i}, t_{n})$, $R_{1}\le C(\tau^{2-\alpha})$, and $a_{k}^{\alpha} = (k+1)^{1-\alpha}-k^{1-\alpha}$.

Lemma 3.1. ^[33] $\alpha \in (0, 1)$, $a_{l}^{\alpha} = (l+1)^{1-\alpha}-l^{1-\alpha}$, $l = 0, 1, 2, \dots$,

(1) $1 = a_{0}^{\alpha} > a_{1}^{\alpha} > a_{2}^{\alpha} > \dots > a_{l}^{\alpha} > 0, \; \lim_{l \to \infty} a_{l}^{\alpha}\to 0,$

(2) $(1-\alpha)^{l^{-\alpha}} < a_{l-1}^{(\alpha)} < (1-\alpha)(l-1)^{-\alpha}, \; l\ge 1$.

For $\gamma \in (1, 2)$, ${}_{0}^{C}D_{t}^{\gamma}u(x, t)\mid _{(x_{i}, t_{n})} = D_{\tau}^{\gamma} \tilde{u} _{i}^{n}+R_{2}(\tilde{u} _{i}^{n})$,

where $\delta_{t}u(x, t) = \frac{\partial u(x, t)}{\partial t}$, $R_{2}\le C(\tau^{3-\gamma})$, and $b_{k}^{\gamma} = (k+1)^{2-\gamma}-k^{2-\gamma}$.

Given the spatial position $x$, it can be seen from the finite difference discretization that the time-fractional derivative of $\tilde{u}(x, t)$ evaluated at time $t$ depends on the value of $\tilde{u}(x, t)$ calculated at all previous times $0$, $\tau$, $2\tau$, $\cdots$, $t$. We call the current time and the previous time the training points and the auxiliary points, respectively.

4.   Numerical examples

In this section, we demonstrate the effectiveness of gfPINNs in solving forward and inverse problems of the multiterm time-fractional Burger-type equation and we compared fPINNs with gfPINNs. We solve the forward problems of the equation and present the experimental results in Section 4.1. We solve the inverse problems and present the experimental results in Section 4.2.

We give a fabricated solution to the problem $u(x, t) = t^{p}\sin(\pi x)$. In the given approximate solution (3.4), the auxiliary function $\rho (\cdot)$ is defined as $\rho (\cdot) = 1 - \begin{Vmatrix}x \end{Vmatrix}_{2}^{2}$. We use the following form of $L^{2}$ relative error:

to measure the performance of the NNs, where $\tilde{u}$ denotes the approximated solution, $u$ is the exact solution, and $(x_{k}, t_{k})$ denotes the $k$-th test point.

We wrote the code in Python and took advantage of the automatic differentiation capability of TensorFlow ^[37]. The stochastic gradient descent Adam algorithm ^[38] was used to optimize the loss function. We initialized the NNs parameters using normalized Glorot initialization ^[39]. Otherwise, when training a neural network, we set the learning rate, the number of neurons, the number of hidden layers, and the activation function as $1\times10^{-3}$, 20, 4, and $\tanh(x)$, respectively.

4.1.   Forward problems

In this section, we consider the the multiterm time-fractional Burger-type equation of the form (2.1) with initial and boundary conditions (2.2). We let $v = 1$, $(x, t)\in [0, 1]\times[0, 1]$, and $g(x) = 0$, considering the smooth fabricated solution $u(x, t) = t^{p}\sin(\pi x)$ and the forcing term

Case 1: We choose $c_{1} = 1$, $c_{2} = 0$, and $\alpha = 0.5$, considering the smooth fabricated solution $u(x, t) = t^{4}\sin(\pi x)$ and the forcing term $f(x, t) = 3.5\frac{\Gamma(5)}{\Gamma(4.5)}t^{3.5} \sin$$(\pi x)+ t^{8} \sin(\pi x) \cos(\pi x) + \pi^{2} t^{4} \sin(\pi x)$. We consider $M_{1}-1$ training points of the spatial domain: $x_{i} = ih_{x}$ for $i = 1, 2, \cdots, M_{1}-1$ and $N$ training points of the time domain: $t_{n} = n\tau$ for $n = 1, 2, \cdots, N$. We do not need to place training points on the initial and boundary since the approximate solution $\tilde{u}(x, t) = tx(1-x)u_{NN}(x, t;\mu)$ satisfies the initial and boundary conditions automatically. For fPINNs, the loss function can be written as

The loss function of gfPINNs can be given as

By substituting Eqs (4.3)–(4.5) into Eq (3.8), we get the gfPINNs loss function $L_{gFW}$ with $W_{F} = 1$, $W_{g_{1}F} = 1$, and $W_{g_{2}F} = 1$. Next, we selected 2000 training points to train fPINNs and gfPINNs and other parameters of the NNs are set to those described at the beginning of this section. Figures 2–4 present a comparison between the predicted solutions from the fPINNs and gfPINNs models and the exact solution of the equation, demonstrating that gfPINNs can effectively solve the equation. Figure 5 shows the absolute errors between the exact solution and the solutions predicted by fPINNs and gfPINNs, and it can be seen that the prediction performance of gfPINNs is better than that of fPINNs. Figure 6 illustrates the $L^{2}$ relative errors of both fPINNs and gfPINNs models for a single experiment as the iteration count varies, showing that while both can achieve errors as low as $10^{-4}$, gfPINNs exhibits comparatively lower error and reduced oscillation.

Case 2: We choose $c_{1} = 0$, $c_{2} = 1$, and $\gamma = 1.5$, considering the smooth fabricated solution $u(x, t) = t^{4}\sin(\pi x)$ and the forcing term $f(x, t) = 2.5\frac{\Gamma(5)}{\Gamma(3.5)}t^{2.5} \sin$$(\pi x)+ t^{8} \sin(\pi x) \cos(\pi x) + \pi^{2} t^{4} \sin(\pi x)$. Similarly, we give the loss function of fPINNs as

For gfPINNs, the loss function can be written as

By substituting Eqs (4.6)–(4.8) into Eq (3.8), we get the gfPINNs loss function $L_{gFW}$ with $W_{F} = 1$, $W_{g_{1}F} = 0.16$, and $W_{g_{2}F} = 0.16$. Next, we selected 2000 training points to train fPINNs and gfPINNs and other parameters of the NNs are set to those described at the beginning of this section. Figures 7–9 present a comparison between the predicted solutions from the fPINNs and gfPINNs models and the exact solution of the equation, demonstrating that gfPINNs can effectively solve the equation. Figure 10 illustrates the absolute errors between the exact solution and the solutions predicted by both fPINNs and gfPINNs, revealing that the gfPINNs exhibit a relatively smaller absolute error. Figure 11 presents the iteration convergence curves for both the fPINNs and gfPINNs models for a single experiment, revealing that while both can achieve $L^{2}$ relative errors of $10^{−4}$ with increasing iterations, the prediction errors of gfPINNs are relatively low and more stable, resulting in superior prediction performance compared to fPINNs.

4.2.   Inverse problems

We use the code that solves the forward problem to solve the inverse problem. We simply add the parameters to be identified in the inverse problem to the list of parameters to be optimized in the forward problem, without changing anything else. In this section, gfPINNs are applied to solve the inverse problems of the multiterm time-fractional Burger-type equation of the form (3.6). We let $v = 1$, $(x, t)\in [0, 1]\times[0, 1]$, $g(x) = 0$, and considering additional concentration measurements at the final time $u(x, 1) = h(x, 1)$. Here, we still consider the smooth fabricated solution $u(x, t) = t^{p}\sin(\pi x)$ and the forcing term of formula (4.2).

Case 1: We choose $c_{1} = 1$ and $c_{2} = 0$. Similarly, we get the gfPINNs loss function $L_{gFW}$ with $W_{I} = 1$, $W_{g_{1}I} = 0.25$, and $W_{g_{2}I} = 0.25$. We set the fractional derivative to be 0.6. We selected 470 training points to train fPINNs and gfPINNs and other parameters of the NNs are set to those described at the beginning of this section. Figures 12–14 display a comparison between the predicted solutions from the fPINNs and gfPINNs models and the exact solution of the equation, demonstrating that gfPINNs can effectively solve the problem. Figure 15 illustrates the absolute errors between the exact solution and the solutions predicted by both fPINNs and gfPINNs, revealing that the gfPINNs exhibit a relatively smaller and more stable absolute error. Figure 16 illustrates the iteration convergence curves for the fPINNs and gfPINNs for a single experiment, indicating that although gfPINNs incur a higher computational cost for solving the inverse problem due to an additional loss term, both models can achieve $L^{2}$ relative errors of $10^{-4}$ as iterations progress, with gfPINNs showing a lower and more stable error curve compared to fPINNs.

Case 2: We choose $c_{1} = 0$ and $c_{2} = 1$. Similarly, we get the gfPINNs loss function $L_{gFW}$ with $W_{I} = 1$, $W_{g_{1}I} = 0.16$, and $W_{g_{2}I} = 0.0001$. We set the fractional derivative to be 1.6. We selected 400 training points to train fPINNs and gfPINNs and other parameters of the NNs are set to those described at the beginning of this section. For fPINNs and gfPINNs, we get the similar conclusion as Case 1 by training the NNs and observing the experimental results. Figures 17–19 display a comparison between the predicted solutions from the fPINNs and gfPINNs models and the exact solution of the equation, demonstrating that gfPINNs can effectively solve the problem. Figure 20 illustrates the absolute errors between the exact solution and the solutions predicted by both fPINNs and gfPINNs, revealing that the gfPINNs exhibit a relatively smaller absolute error. Figure 21 compares the $L^{2}$ relative errors of fPINNs and gfPINNs for a single experiment as iterations progress, revealing that while gfPINNs incurs a higher computational cost due to an additional loss term, both models can achieve an $L^{2}$ relative error of $10^{-3}$, with gfPINNs demonstrating a lower and more stable error curve than fPINNs.

5.   Conclusions

In this paper, the effectiveness of gfPINNs in solving the forward and inverse problems of the multiterm time-fractional Burger-type equation is verified through numerical examples. The $L^{2}$ relative errors for solutions predicted by both fPINNs and gfPINNs can achieve $10^{-4}$ for forward problems and $10^{-3}$ or even $10^{-4}$ for inverse problems. The experimental results indicate that gfPINNs demonstrate relatively lower and more stable errors with the increase of training iterations, thereby enhancing prediction performance. Nonetheless, the inclusion of an additional loss term in gfPINNs may result in a higher computational cost, such as when solving inverse problems, fPINNs exhibit faster convergence compared to gfPINNs.

Author contributions

Shanhao Yuan, Yanqin Liu, Qiuping Li and Chao Guo: Conceptualization, Methodology; Yibin Xu, Shanhao Yuan and Yanfeng Shen: Software, Visualization, Validation; Shanhao Yuan: Writing–Original draft preparation; Yanqin Liu: Writing–Reviewing & editing. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

We appreciated the support by the Natural Science Foundation of Shandong Province (ZR2023MA062), the National Science Foundation of China (62103079), the Belt and Road Special Foundation of The National Key Laboratory of Water Disaster Prevention (2023491911), and the Open Research Fund Program of the Data Recovery Key Laboratory of Sichuan Province (DRN19020).

Conflict of interest

The authors declare that they have no conflicts of interest.