Parameter Estimation for a Type of Fractional Diffusion Equation Based on Compact Difference Scheme

: Numerical solution and parameter estimation for a type of fractional diffusion equation are considered. Firstly, the symmetrical compact difference scheme is applied to solve the forward problem of the fractional diffusion equation. The stability and convergence of the symmetrical difference scheme are presented. Then, the Bayesian method is considered to estimate the unknown fractional order α of the fractional diffusion equation model. To validate the efﬁciency of the symmetrical numerical scheme and the estimation method, some simulation tests are considered. The simulation results demonstrate the accuracy of the compact difference scheme and show that the proposed estimation algorithm can provide effective statistical characteristics of the parameter.

Recently, the research about the inverse problem has attracted more and more interest, one can refer to [22][23][24] and the references therein. Most of the literature on the inverse problem of the fractional differential equations exploited deterministic techniques, such as exact matching, least squares optimizations, without considering measurement error and numerical error. In general, very roughly speaking, one may split the corresponding approaches of recovering the order into the following two categories: solving the corresponding inverse problems analytically or using numerical-analytical methods [25][26][27][28][29][30][31], or using some more soft, metaheuristic, or statistical optimization and regularization techniques [32][33][34][35][36][37]. However, the observations are usually contaminated with measurement error, and the forward problem of the models will bring the numerical error, so the uncertainties are non-ignorable. Some studies for the inverse problem under uncertainties have been provided, such as the extended maximum likelihood method [38], sensitivity analysis method [39], and the Bayesian method [40,41].
The work of this paper was mainly motivated by the study in [41]. However, in the paper, Fan et al. considered homogeneous fractional diffusion equation, and the center Box difference method was used to obtain the numerical solution of the equation. An inhomogeneous fractional diffusion equation was considered in this paper, and a highorder compact difference scheme instead of the second-order center Box difference method was used to obtain the corresponding numerical solution. We considered the Bayesian method to estimate the parameter of the fractional diffusion equation model based on a compact difference scheme. Firstly, a high-order compact difference scheme is employed to solve the forward problem of the fractional diffusion equation. The convergence and stability of the difference scheme are proved, where the difference scheme has 2 − α order in the temporal direction and 4 order in the spatial direction. Then, the Bayesian method is considered to estimate the unknown fractional order α of the model, where some additional observations contaminated with measurement errors must be provided. To validate the efficiency of the numerical scheme and the estimation method, some simulation tests are considered. The simulation results demonstrate the accuracy of the compact difference scheme and show that the proposed estimation algorithm can provide effective statistical characteristics of the parameters.
The outline of this paper is as follows. In Section 2, mathematical model of the problem is introduced, and a compact difference scheme is constructed to solve the forward problem. Section 3 introduces the Bayesian method for the inverse problem. In Section 4, numerical tests are presented to validate the performance of the estimation algorithm. Section 5 gives some discussion and remarks.

Mathematical Model of the Problem and the Compact Difference Scheme
This section described the mathematical model of the problem, and a compact difference scheme was constructed to solve the forward problem.

Mathematical Model of the Problem
Consider the following fractional parabolic equation, where 0 < α < 1 and Q > 0 are constants, φ, ψ 1 , and ψ 2 are given functions. The time fractional partial derivative ∂ α u ∂t α is defined in the Caputo sense as the following where Γ() is the gamma function. If α is a known parameter, Equations (1)-(3) represent a forward problem; extensive studies have been provided. However, in some practical situations, we cannot know the exact values of α. Under this situation, additional data must be provided to estimate the unknown parameter, which can be called the inverse problem. The additional data may be observed at some measured time for a fixed point x f ∈ [0, 1]. For simplicity, let x f = 1 in this paper, and we have the following observations, Let F be an operator, F : θ → u(1, t; θ), where θ = α; then, the inverse problem can be seen as the solution of the following equation, In fact, the observations may be contaminated with measurement errors inevitably. We will consider the Bayesian approach to deal with the inverse problem.

The Compact Difference Scheme for the Forward Problem
Assuming θ in the operator F to be known in this subsection, we construct the highorder compact difference scheme to solve (6).
be the grid function space defined on Ω hτ . Introducing the following notations One can use Lemma 1 to deal with (4).
And d j satisfies Lemma 2.

Lemma 2 ([43]
). Assume 0 < α < 1, then it holds that (1) d j decreases monotonically as j increases, and 0 < d j ≤ 1; The following Lemma 3 will play the key role to construct the high-order compact difference scheme.

Theorem 1. The difference scheme in Equations
The proof of Theorem 1 can be referred to the corresponding proof in [45].
|e i |; then, we can obtain the following error equations where R k i = O(τ 2−α + h 4 ). To prove the following theorems, introducing the following notations is the solution of (1)-(3), and {u k i |0 ≤ i ≤ M, 0 ≤ k ≤ N} is the solution of (14)- (16). Then, we have where C is a positive constant independent of h and τ.
Proof. Multiplying (17) with hδ α t e k i and summing up for i from 1 to M-1, one can obtain Then, we will estimate each term of (21). From the discrete Green formula, we have where inequality −(a + b) 2 ≥ −2(a 2 + b 2 ) is used. By the discrete Green formula, one has From Cauchy-Schwarz inequality, we have Inserting Equations (22)-(24) into (21), and take ε = 4/3, we obtain Noticing R k i = O(τ 2−α + h 4 ), from Lemma 6 in [45], we have where C 1 is a positive constant independent of h and τ. Then, we have The proof is completed.
Similar to the proof of Theorem 2, we have the following result about the stability of the difference scheme in Equations (14)- (16).

Parameter Estimation for the Inverse Problem
Considering the measurement error and numerical error in the observations, then (6) has the following form, where ε ∼ N(0, σ 2 ), and σ 2 is assumed to be a known parameter. θ = α is the parameter to be estimated from the observations governed by (26).
Assume α ∼ U(0, 1) is the prior probability distributions. We have the prior probability density function of θ, From (26), we have the following likelihood function, where Y is the observation vector, n is the length of Y. By the Bayes theorem, the posterior probability density function is The Bayesian method for estimating θ = α is as follows, (1) Draw an initial value of θ 0 ∼ U(0, 1), where N sim is the number of iterations.
In the Bayesian method, we will use the posterior mean as the estimation of θ. (29) can be obtained through (5), where a higher compact difference scheme is used. However, the traditional Crank-Nicolson difference scheme and other algorithms such as finite element method and finite volume method can also be used to solve (5). The advantage of the compact difference is that it is highly accurate, easy to program, and easily parallel. The advantage of the Bayesian estimation algorithm is that it can deal with contaminated observations and inhomogeneous fractional diffusion equation, and the most disadvantage of the Bayesian estimation algorithm is that it is time consuming.

Simulation
This section considers the following example. The configuration of the computer is Intel(R) Core(TM) i5-8265U CPU @ 1.60 GHz 1.80 GHz, Memory 4 GB (Lenovo, Beijing, China). All the codes are compiled in MATLAB 2009(a).

Simulation for the Forward Problem
As for the forward problem, we will validate the performance of the scheme (14)- (16) for solving (30). Introducing the following maximum error, and we define the following convergence order as For calculating R τ , h needs to be fixed and small enough. While for R h , τ needs to be fixed and small enough.
To compare the errors between the numerical solution and the exact solution of (30), Figure 1 was provided. Taking temporal and spatial step sizes as τ = h = 1/10, 1/20, 1/40, 1/80 respectively, let α = 0.2, and we can see that the value of the error is within [0, 1.6 × 10 −5 ]. The numerical solution can approximate the exact solution well enough.  Table 1 shows the maximum errors in the temporal directions for α = 0.2, 0.6, and 0.8, respectively, where we fixed the spatial step to be h = 1/2000. We can see that the temporal convergence order is 2 − α.  Table 2 shows the maximum errors in the spatial directions for α = 0.2, 0.6, and 0.8, respectively, where we fixed the temporal step τ = 1/10,000. In most of the cases, we can conclude that the spatial convergence order is 4. However, we can see that the convergence of the considered numerical scheme is not steady under the condition that α = 0.8. Table 2. Maximum errors and convergence order for solving (30) in the spatial direction, where Q = 1.8, τ = 1/10,000, and α = 0.2, 0.6, and 0. 8

Simulation for the Inverse Problem
As for the inverse problem, we will testify to the validity of the Bayesian method provided. To obtain the observations Y governed by (26), let θ = α = 0.2 in Equation (1) be the given value (note: θ = α = 0.6 and 0.8 will take the same conclusion, only θ = α = 0.2 is considered in this subsection for the sake of brevity). Firstly, we applied the compact difference scheme (14)- (16) to discrete (30), where M = N = 100. A series of points {u k M }, k = 0, 1, · · · , N were generated. Then, we picked up every fourth (two or one) point from the generated points, finally add the noise on each data, and we get the observations with length n = 25 (n = 50, n = 100).
However, the length of observations n and σ 2 may affect the estimation results of θ. In the following tests, we will evaluate such effects caused by different values of n and σ 2 . In the following tests, the length of the Markov chain is taken to be N sim = 5000, and the burn in number is taken to be 1000.
Let n = 100 be fixed, and assume σ ∈ [1 × 10 −8 , 5 × 10 −7 ]. For different σ, we will inspect the posterior mean, standard error, and quantile of α. Table 3 was obtained, and we can see that the measurement error with bigger σ may take more uncertainty of the estimation results. To inspect the effect of different n for α, let σ = 2 × 10 −7 be fixed, assume n = 25, n = 50, and n = 100. For different n, we also inspect the posterior mean, standard error, and quantile of α. Table 4 was obtained, when n = 50 and n = 100, the different values of n make little impact on the estimation results. While when n = 25, there is more uncertainty in the estimation results of α. We can draw the conclusion that there is little influence for estimating α under bigger values of different n. From the computational time, it can be seen that we need more estimation time when a bigger n is given. In the tests of simulation for the inverse problem, we set the seed as 'sum(100*clock)' in the 'rand' function to specify the settings of the random number generator. Then, the estimation results can always be the same.

Conclusions
This paper considered both the forward problem and inverse problem of a type of fractional diffusion equation. Firstly, a compact difference scheme was provided to solve the forward problem of the fractional diffusion equation. The convergence and stability of the difference scheme were proved, and the convergence order is O(τ 2−α + h 4 ). Secondly, as for the inverse problem, a Bayesian algorithm was provided to estimate the fractional order α of the diffusion equation under the condition that some additional observations contaminated with measurement errors were provided. Numerical experiments were provided to testify the convergence results both in the temporal direction and in the spatial direction and the estimation efficiency of the Bayesian algorithm.