An artificial neural network approach for a class of time-fractional diffusion and diffusion-wave equations

1.   Introduction

Fractional calculus is an important branch of mathematics. As its theory has been successfully applied to various fields in recent years, people have gradually found that fractional calculus can describe some "abnormal" phenomena in natural science and engineering applications. For example, mechanical behavior of viscoelastic materials ^[1], plasma motion under high temperature and high pressure ^[2], abnormal diffusion of extremely cold atoms or cell protoplasm ^[3], etc. Time fractional partial differential equation is an important model. When its fractional order $\alpha \in (0, 1)$, it is called time-fractional diffusion equation; when $\alpha \in (1, 2)$, it is called time-fractional diffusion-wave equation. The time-fractional diffusion equation can describe the abnormal slow diffusion phenomenon in porous media with fractal structure ^[4], while the time-fractional diffusion-wave equation can describe the propagation of mechanical scattered waves in viscoelastic media with power-law creep characteristics ^[5].

It is well known that the analytical solution of time-fractional diffusion and diffusion-wave equations is notoriously difficult to obtain. In order to promote the development and progress of science, a great quantity of researchers and scholars have studied the numerical approximation method of time-fractional diffusion and diffusion-wave equations. In recent years, Gu and his collaborators have proposed many effective methods to solve the time-fractional diffusion equations, such as quadratic spline collocation method ^[6], Lagrange interpolation method ^[7], parallel-in-time iterative algorithm ^[8] and the graded L1 scheme method ^[9]. For the time-fractional diffusion-wave equations, many classical numerical methods have been proposed, including finite difference method ^[10,11], element-free Galerkin method ^[12], wavelets method ^[13,14], meshless local collocation method ^[15], etc.

Moreover, artificial neural network method is also proposed to solve fractional differential equations. With the rapid development of machine learning technology, the demand for solving fractional differential equations by using artificial neural network method has increased by leaps and bounds, and some excellent results have been achieved. In 2017, Qu ^[16] applied cosine radial basis function neural network to solve fractional differential equations with initial value or boundary value problems for the first time. In 2018, Rostami et al. ^[17] adopted an iterative algorithm combining power series and neural network method to approximately seek the solution of high-order linear and ordinary differential equations. Rizaner et al. ^[18] proposed a numerical approch for approximately solving the first-order initial value problem using unsupervised radial basis function network, and verified the effectiveness with some numerical examples. In 2020, Hadian-Rasanan et al. ^[20] solved different types of Lane-Emden equation by constructing a single-layer orthogonal network with a fractional Legendre function as the activation function of the hidden layer, and used the Levenberg-Marquardt algorithm to train them. The above literature is to solve fractional ordinary differential equations by making use of artificial neural network method. However, there are few literatures on solving time-fractional partial equations by using artificial neural networks. Jafarian et al. regarded the artificial neural network method with a generalized sigmoid function as the cost function as the iterative scheme for solving linear fractional ordinary differential equations in ^[19]. Recently, Qu et al. ^[21] proposed a neural network method based on cosine basis and sine basis to solve a class of time-fractional diffusion and diffusion-wave equation with zero boundary conditions. Deep neural network methods for solving forward and inverse problems of time-fractional diffusion equations with conformable derivatives were introduced by Ye et al. ^[22]. In this paper, a new numerical method is proposed for a class of time-fractional diffusion and diffusion-wave equations with Dirichlet boundary conditions by combining Taylor series method with artificial neural network.

In this paper, We consider the time-fractional diffusion equation with the following form:

and the time-fractional diffusion-wave equation with the form as follows:

For the time-fractional diffusion equation (1.1) and the time-fractional diffusion-wave equation (1.2), we establish a single-layer neural network model with Taylor series of different power terms as neurons, and give the error iteration formula and neural network weight update formula. It is worth mentioning that we have adaptively adjusted the weights of different error terms, which will be shown in Section 4.

The main contributions of this paper are as below:

● The single layer functional link artificial neural network based on Taylor series is established for the first time to solve the time-fractional diffusion and diffusion-wave equation with more general form.

● The error function composed of initial error term, boundary error term and internal error term is constructed to update the unknown parameters in the neural network. In addition, the adaptive weight adjustment method is proposed to balance the different error terms.

● Compared with the traditional iterative numerical method, the proposed neural network method can calculate the solution at any position in the domain without repeating the whole iterative process.

The rest of this paper is arranged as follows. Section 2 briefly introduces the related definition and important properties of fractional calculus. In Section 3, we present the main idea of our neural network methods for solving the time-fractional diffusion and diffusion-wave equations (1.1), (1.2). In Section 4, we propose an adaptive weight adjustment method to adjust the weights of different error terms in the iteration. After that, four numerical examples are given to illustrate the effectiveness of the method in Section 5. Finally, we summarize the main work of our paper in Section 6.

2.   Preliminaries

This section gives some basic definition and related properties about fractional calculus; for more details see ^[23].

Definition 2.1. The Riemann-Liouville fractional integral of order $\alpha > 0$ for a function $f(x) \in L^1((0, Q])$ is defined as

where $\Gamma(\cdot)$ denotes the Gamma function.

Definition 2.2. The Caputo fractional derivative of order $\alpha > 0$ for a function $f(x) \in AC((0, Q])$ with $n \in \mathbb{N}^{+}$ is defined as

Since the derivative of a constant is 0, the following useful property holds in the Caputo sense:

Definition 2.3. The Mittag-Leffler function is defined by

where $z, \beta \in \mathbb{C}$, $\alpha > 0$, and $\mathbb{C}$ is the usual set of complex numbers. Particularly, $E_{1, 1}(z) = exp(z)$, that is, the common exponential function.

In order to facilitate the calculation of the Mittag-Leffler function later, the definition of the extended sine and cos functions and related properties are given ^[21,23]. The calculation code of the Mittag-Leffler function can be found in ^[24].

Definition 2.4. The extended sine function $Sin_{\alpha, \beta}(x)$ and cos function $Cos_{\alpha, \beta}(x)$ are defined as

where $x, \beta \in \mathbb{C}$, $\alpha > 0$. According to the Euler formulas, we have

In particular, we have $Sin_{1, 1}(x) = sin(x)$, $Cos_{1, 1}(x) = cos(x)$.

Lemma 2.1. For the Caputo fractional derivative of some special functions, we have three useful results as follows:

and

3.   Methodology

With regard to the time-fractional diffusion equation (1.1) and the time-fractional diffusion-wave equation (1.2), the proposed method is described in detail in this section. In order to solve the above two problems, we built a single layer neural network with a series of rapidly convergent components of Taylor series as neurons, as shown in Figure 1.

In simple terms, we construct the $k$-th test solution in the following form:

with the first few Taylor components:

and the higher order components of Taylor series may be expressed as the modalities as below:

Here, $N \in \mathbb{N}_{+}$, $\frac{(N+1)(N+2)}{2}$ and $\theta_{\frac{i(i+1)}{2}+j+1}^{k}$ represent the highest power of the series, the number of neurons and the unknown weights of the neural network in the $k$-th iteration, respectively.

Since the weights of the neural network are unknown, in the first iteration, the weights are randomly preseted by a Gaussian distribution. In the subsequent iterations, we renew the weights through the error function. For the initial boundary value problems like Eqs (1.1) and (1.2), the error function is composed of the following three terms:

where $\left\|\cdot\right\|_{F}$ is Frobenius matrix norm and $\left\|\cdot\right\|_{2}$ represents vector 2-norm. Let $x_{m} = \frac{(m-1)P}{N_{m}-1}$, $t_{n} = \frac{(n-1)Q}{N_{n}-1}$, $m = 1, 2, \cdots, N_{m}$, $n = 1, 2, \cdots, N_{n}$. The matrix $E_{i n}^{k}$ is

with

for $2 \leq m \leq N_{m}-1$, $2 \leq n \leq N_{n}$. The vector $E_{i n i}^{k}$ is

with

The matrix $E_{b o u}^{k}$ is

with

and

As for the internal points, by substituting the $k$-th test solution (3.1) into Eq (3.4) and making use of Definition 2.2, we further obtain that

Concerning the initial points, we gain the initial error by substituting the $k$-th test solution (3.1) into Eq (3.6)

With regard to the boundary points, substituting the $k$-th test solution (3.1) into Eqs (3.8) and (3.9) we acquire the boundary errors

and

Here, we use the gradient descent method to update the $k$-th unknown weights $\theta_{\frac{i(i+1)}{2}+j+1}^{k}$. For more details, please refer to ^[25].

The weight update formula is as follow:

with

Here, $\theta_{\frac{i(i+1)}{2}+j+1}^{0}$ are some numbers randomly generated by Gaussian distribution, $K$ is an integer representing the maximum number of iterations, and $\eta$ is the learning rate of neural network.

3.1.   The time-fractional diffusion-wave equation (1.2)

For the time-fractional wave equation (1.2), we just need to rewrite Eqs (3.10) and (3.14) into

and

with

respectively.

4.   Adaptive weight adjustment

It is well known that a system of partial differential equations can have infinite solutions without specifying proper bounds or initial conditions ^[26]. So is the time-fractional diffusion and diffusion-wave equation. Therefore, we believe that the gradient of internal error ($\nabla_{\bar{\theta}} E_{in}$), the gradient of initial error ($\nabla_{\bar{\theta}} E_{ini}$) and the gradient of boundary error ($\nabla_{\bar{\theta}} E_{bou}$) should match each other in magnitude, otherwise the network will eventually tend to learn any solution satisfying the equation due to the dominant effect of the gradient of internal error, and then, it is easy to fall back to the wrong prediction. In this section, let us talk over how to balance the interaction between internal error term ($E_{in}$), initial error term ($E_{ini}$) and boundary error term ($E_{bou}$). For simplicity, $\bar{\theta}$ denotes $\theta_{\frac{i(i+1)}{2}+j+1}$ in this section.

Recall the error function Eq (3.2) and weight update formula Eq (3.14) proposed in Section 3, we get that

Usually, the internal error term dominates the overall. To balance the interaction between different terms in this error, a simple method is to weight the initial error term and the boundary error term. Thus, Eq (3.2) and Eq (3.14) can be rewritten as

and

respectively, $\omega_{1}^{k}, \omega_{2}^{k}$ are the weighted coefficients. This method is similar to the weighted optimization method, which improves the optimization effect by penalizing the initial error term and boundary error term.

For the time-fractional diffusion-wave equation (1.2), Eq (3.16) can be rewritten as

Next, we will discuss how to choose the weighted coefficients $\omega_{1}^{k}, \omega_{2}^{k}$. For different models, the optimal constant can be very different, which means that we cannot find a fixed empirical formula that can be transferred between different time-fractional diffusion and diffusion-wave equations. Also, the error function always consists of different parts that are used to constrain the equation. It is impractical to manually assign different weights to different terms of the error function. The traditional parameter adjustment method can neither guarantee the efficiency nor ensure the optimization to the global optimal value. Therefore, inspired by the learning rate annealing algorithm in ^[27], we propose an adaptive weight coefficients adjustment method.

At the $(k+1)$th iterations, we calculate the estimated values of $\omega_{1}^{k}$ and $\omega_{2}^{k}$ by the formula as below:

where $\overline{\left|\nabla_{\bar{\theta}}E_{in}\right|}$, $\overline{\left|\nabla_{\bar{\theta}} E_{ini}\right|}$ and $\overline{\left|\nabla_{\bar{\theta}} E_{bou}\right|}$ denote the means of $\left|\nabla_{\bar{\theta}}E_{in}\right|$, $\left|\nabla_{\bar{\theta}} E_{ini}\right|$ and $\left|\nabla_{\bar{\theta}} E_{bou}\right|$.

Then, the weighted coefficients of $(k+1)$ iterations are updated by the weighted moving average method:

with $\xi$ is given as 0.1.

Finally, the algorithm of this paper is summarized as Algorithm 1.

5.   Numerical illustrations

In this section, we give four examples to illustrate the applicability of our proposed neural network method in solving time-fractional diffusion equation (1.1) and time-fractional diffusion-wave equation (1.2). All experiments were performed on a computer with the following configurations: Intel(R) Core(TM) i5-8265U CPU @ 1.60 GHz 1.80 GHz using MATLAB R2019a.

5.1.   The time-fractional diffusion equation (1.1)

Example 1. We consider Eq (1.1) with the conditions: $\lambda(x) = \mu(x) = 0$, $\varphi(x) = sin(\frac{\pi x}{P})$ and $f(x, t) = t^{1-\alpha} sin(\frac{\pi x}{P}) E_{1, 2-\alpha}(t) - (\frac{c \pi}{P})^{2} exp(t) sin(\frac{\pi x}{P})$. The exact solution of the equation may be calculated by means of Definition 2.4 and Lemma 2.1.

At first, the parameters of the model are set to the following: $P = 1$, $Q = 1$, $c = 1$ and $0 < \alpha < 1$. Then, the parameters of the proposed method we use are as follows: $N = 5$, $m = 5$, $n = 5$, learning rate $\mu = 4 \times 10^{-4}$, and accuracy threshold $\epsilon = 1 \times 10^{-7}$. When $\alpha = 0.9$, Figure 2 represent the comparisons between the exact solution and the neural network solution, and the number of iterations is 2000, 5000, 10000, 20000, respectively. We observe that the approximate solution gets better and better as the number of iterations increases. Then, we check the approximate solution on the grid of $N_{m} \times N_{n} = 21 \times 21$ in Figure 3a. The variation of the error value with the number of iterations is shown in Figure 4. Also, in Figure 3, a comparison of the exact and the approximate solution is given under $\alpha = 0.3, 0.5$ and $0.7$. In addition, we list the numerical solutions and absolute errors of some specific test points for $\alpha = 0.3, 0.5, 0.7$ and $0.9$ in Table 1. We can see that we get similar results for different fractional order $\alpha$. Eventually, the weights of the neural network are listed in Table 2.

Example 2. In this example, we consider Eq (1.1) with the conditions as follows:

in this situation, the exact solution is

Let $P = 1$, $Q = 1$, $c = 1$ and $0 < \alpha < 1$. In our proposed neural network method, the parameters are as follows: $N = 7$, $N_{m} = 5$, $N_{n} = 5$. We set the learning rate $\eta$ and accuracy threshold $\epsilon$ to $1 \times 10^{-4}$ and $1 \times 10^{-7}$, respectively. Through 20000 iterations, we obtain the neural network solution of Example 2 for different fractional orders $\alpha = 0.3, 0.5, 0.7$ and $0.9$. Figure 5 give a comparison of neural network solutions and exact solutions for different fractional orders. In Table 3, we list the numerical solutions and absolute errors of several test points under different $\alpha$. Next, we subdivided the computational grid into $N_{m} \times N_{n} = 21 \times 21$ and tested our neural network solution on Figure 6. In Figure 7, we plot the variation of the error for different fractional order. Finally, the weights of neurons after iteration are listed in Table 4.

Example 3. Let $P = 1$, $Q = 1$, $c = 1$ and $1 < \alpha < 2$, we consider Eq (1.2) with $\lambda(t) = 0$, $\mu(t) = P^{2}t^{3}$, $\varphi(x) = 0$, $\psi(x) = 0$ and $f(x, t) = \frac{\Gamma(4)}{\Gamma(4-\alpha)} t^{3-\alpha} x^{2} + 2c^{2}t^{3}$. Through calculation, we can gain the exact solution $u(x, t) = x^{2}t^{3}$.

In our proposrd method, we set $N = 7$ and sample gird $N_{m} \times N_{n} = 5 \times 5$. The learning rate $\eta$ and accuracy threshold $\epsilon$ are set to $1 \times 10^{-4}$ and $1 \times 10^{-7}$, respectively. With 20000 iterations, we acquire a good neural network solution. Figure 8 show a comparison of the neural network solution and the exact solution on sampling grid, their comparison on denser grid $N_{m} \times N_{n} = 21 \times 21$ and the change of the error value with the number of iterations, respectively. In Table 5, we list the numerical solutions and absolute errors of several test points for $\alpha = 1.9$. Also, Table 6 lists the weights of the proposed neural network in this example. The results illustrate that the proposed neural network method is also very effective for seeking the solution of the time-fractional diffusion-wave equations. (On account of too many graphs, we only list the situation of $\alpha = 1.9$ here.)

Example 4. In this example, we consider Eq (1.2) with the conditions as follows:

in this situation, the exact solution is

Let $P = 1$, $Q = 1$, $c = 1$ and $1 < \alpha < 2$, other parameters are listed as follows: $N = 5$, $N_{m} = N_{n} = 5$, $\eta = 6 \times 10^{-5}$ and $\epsilon = 1 \times 10^{-7}$. Similar to Example 3, we only show the case of $\alpha = 1.5$ by 40000 iterations. In Figure 9, we show the comparison of the neural network solution and the exact solution on sampling grid $N_{m} \times N_{n} = 5 \times 5$, their comparison on denser grid $N_{m} \times N_{n} = 21 \times 21$ and the change of error value with the number of iterations, respectively. In addition, we list the numerical solutions and absolute errors of several test points in Table 7. The weights of the neural network in this example are listed in Table 8.

6.   Conclusion

In this paper, different power terms of Taylor polynomial are used as neurons of the proposed artificial neural network to solve the time-fractional diffusion and diffusion-wave equations. The initial weight of the neural network is given randomly according to the Gaussian distribution, and what is more the weight in each iteration is updated by the error function. In order to balance the contribution of initial error term, boundary error term and internal error term in the error function, we weight different error terms, and give an adaptive weight adjustment algorithm. The trained network is found to be effective by testing other points within the solution region. In addition, the value of the error function decreases as the number of iterations increases. This illustrates that our method is very effective for solving time-fractional diffusion and diffusion-wave equations. In the future work, we will consider to extend our methods for solving other time-fractional partial differential equations, such as the time-fractional mixed diffusion and diffusion-wave equations ^[28], time-fractional diffusion equations with a time-invariant type variable order ^[29] and the generalized time space fractional diffusion equations with variable coefficients ^[30].

Acknowledgments

The authors express their sincere thanks to the anonymous reviewer for his/her careful reading of the paper, giving valuable comments and suggestions. It is their contributions that greatly improve the paper. The authors also thank the editors for their kind help.

This work was supported by the National Natural Science Foundation of China (Nos. 11801452, 11701456); Shaanxi Fundamental Science Research Project for Mathematics and Physics (Grant No. 22JSQ020); Fundamental Research Project of Natural Science in Shaanxi Province General Project (Youth) (Nos.2019JQ-415, 2019JQ-196); Chinese Universities Scientific Funds (Nos. 2452022104, 2452022105, 2452022372); and Innovation and Entrepreneurship Training Program for College Students of Northwest A & F University (202110712186, 202210712048, S202210712377).

Conflict of interest

The authors declare there is no conflict of interest.