# Numerical Identification of Multiparameters in the Space Fractional Advection Dispersion Equation by Final Observations

Authors: Dali Zhang; Gongsheng Li; Guangsheng Chi; Xianzheng Jia; Huiling Li
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

This paper deals with an inverse problem for identifying multiparameters in 1D space fractional advection dispersion equation (FADE) on a finite domain with final observations. The parameters to be identified are the fractional order, the diffusion coefficient, and the average velocity in the FADE. The forward problem is solved by a finite difference scheme, and then an optimal perturbation regularization algorithm is introduced to determine the three parameters simultaneously. Numerical inversions are performed both with the accurate data and noisy data, and several factors having influences on realization of the algorithm are discussed. The inversion solutions are in good approximations to the exact solutions demonstrating the efficiency of the proposed algorithm.


---

## Body

## 1. Introduction

Many diffusion processes in nature and engineering, such as contaminants transport in the soil, oil flow in porous media, long distance transport of pollutants in the groundwater, and so forth, are referred to as anomalous diffusion, where the particle plume spreads faster or slower than predicted by the classical diffusion model. In recent two decades, the space FADE has been found to be an efficient mathematical model to describe such anomalous diffusion phenomena (see [1–8], e.g.). If the usual second-order derivative ∂2/∂x2 in space is replaced by a fractional derivative of order 1<α<2, denoted as ∂α/∂xα, we get the space FADE given as
(1.1)∂c∂t=D∂αc∂xα-v∂c∂x+S(x,t),
where c=c(x,t) denotes concentration of the solute in space point x and time t, 1<α<2 denotes the order of the fractional derivative, D>0 is the dispersion coefficient, v>0 denotes the average velocity, and S(x,t) is a source term. The derivative ∂αc/∂xα here denotes a fractional derivative at space point x from the left hand-side of the domain in the meaning of Grünwald-Letnikov or Riemann-Liouville (see [9], e.g.).There are quite a few research works on the fractional differential equations, refer to [10–16] for numerical methods and model simulations on (1.1), and there are still a few studies in the time fractional diffusion equations recently (see [17–22], e.g.). However, research on inverse problems for the space FADE has not been paid much attention to our knowledge. Wei et al. [23] studied an inverse source problem in the spatial fractional diffusion equation and solved the inverse problem numerically based on the best perturbation method with Tikhonov regularization. Chi et al. [24] considered an inverse problem of determining the space-dependent source function f(x) in (1.1) in the case of S(x,t)=e-tf(x) with final observations and presented numerical inversions by applying the optimal perturbation regularization algorithm.In this paper, we will deal with an inverse problem of simultaneously determining the fractional order, the diffusion coefficient, and the average velocity in (1.1) on a finite domain also with the final observations. Rodrigues et al. [25] considered an inverse problem of simultaneously determining the diffusion coefficient and the source term in 1D integer-order diffusion equation by the conjugate gradient method with aids of an adjoint problem, and there seem to have few similar works reported in the known literatures especially for the space fractional diffusion equations. It is noticeable that the inversion problem here for simultaneously identifying these three kinds of model parameters becomes difficult as compared with those problems of determining a single unknown, and regularization strategies have to be utilized to overcome the ill-posedness. The forward problem is solved numerically by utilizing an implicit difference scheme based on the shifted Grünwall formula which was discussed in [10, 12, 24]. In regard of simultaneous inversion for the three kinds of parameters, we will apply the optimal perturbation regularization algorithm which was developed in [26–29]. However, since the inverse problem is different, the inversion algorithm has to be changed in detail so as to get effective solutions. The inversion algorithm is performed not only in the case of using the accurate data, but also with the random noisy data. The inversion solutions give good approximations to the exact solutions showing that the inversion algorithm is stable and suitable for the inverse problem of determining multiparameters in the FADE. The paper is arranged as follows.In Section2, an implicit finite difference scheme for solving the forward problem is given, and a numerical example is presented to support the numerical method. Section 3 gives the inverse problem of determining the parameters of the fractional order, the diffusion coefficient, and the flow velocity simultaneously by final observations, and the optimal perturbation regularization algorithm is introduced. In Section 4, numerical inversions both with the accurate data and noisy data are implemented, and the factors influencing realization of the inversion algorithm are discussed. Finally, several concluding remarks are given in Section 5.


## 2. Numerical Solution to the Forward Problem

### 2.1. The Implicit Difference Scheme

For0<x<l and 0<t<T, consider an initial boundary value problem for the FADE given by (1.1), where the initial boundary value conditions are given as
(2.1)c(x,0)=g(x),0≤x≤l,c(0,t)=c(l,t)=0,0≤t≤T,
where the function g(x) is assumed to be continuous on [0,l] and g(0)=g(l)=0.Let the space step beh=l/M and the time step be r=T/N, with xi=ih and tn=nr for i=0,1,…,M and n=0,1,…,N. According to standard Grünwald formula (see [9], e.g.), a discrete approximation to the fractional derivative is given as follows:
(2.2)∂αc(x,t)∂xα=limh→0h=x/M1hα∑k=0M(-1)k(αk)c(x-kh,t),
with which a modified approximation which is called the shifted Grünwald formula is obtained with the help of Gamma function (see [10, 12], e.g.)
(2.3)∂αc(x,t)∂xα=1hα∑k=0MΓ(k-α)Γ(-α)Γ(k+1)c(x-(k-1)h,t)+O(h).
Then we have a difference scheme by (1.1) with the initial boundary value conditions (2.1)
(2.4)cin+1-cinr=-vcin+1-ci-1n+1h+D1hα∑k=0i+1Γ(k-α)Γ(-α)Γ(k+1)ci-k+1n+1+Sin+1,c0n=0,cMn=0,ci0=g(xi),
where cin≈c(xi,tn), Sin≈S(xi,tn) for i=0,1,…,M and n=0,1,…,N-1.Denote(2.5)gk=Γ(k-α)Γ(-α)Γ(k+1),
and p=vr/h and q=Dr/hα, then this can be rearranged to yield
(2.6)-qg0ci+1n+1+(1+p-qg1)cin+1-(p+qg2)ci-1n+1-q∑k=3i+1gkci-k+1n+1=cin+rSin+1,
for i=1,2,…,M-1 and n=0,1,…,N-1. Let cn=(c1n,c2n,…,cM-1n)T, Sn=(S1n,S2n,…,SM-1n)T, and A=(aij) is the coefficient matrix given by (2.6), where the coefficients aij are defined by
(2.7)aij={0,j>i+1,-qg0,j=i+1,1+p-qg1,j=i,-p-qg2,j=i-1,-qgi-j+1,j<i-1,
for i=1,2,…,M-2 and j=1,2,…,M-1, and
(2.8)aM-1,j={-qgM-j,j=1,2,…,M-3,-p-qg2,j=M-2,1+p-qg1,j=M-1,
for i=M-1, respectively. Thus, we get the implicit finite difference scheme in matrix form given as follows:
(2.9)Acn+1=cn+rSn+1,n=0,1,…,N-1.For the above difference scheme, similarly as done in [10, 12], it is not difficult to get its convergence and stability with the help of spectrum analysis to the coefficient matrix.Lemma 2.1 (see [10, 12]).
The implicit difference scheme defined by (2.9) is unconditionally stable, and the difference solution is convergent to the exact solution of the forward problem as h,r→0 for finite time domain.


### 2.2. Numerical Testification

In this subsection, an example is presented to show effectiveness of the finite difference scheme (2.9) for solving the forward problem numerically. For the forward problem (1.1)-(2.1), set l=π, T=1, and take the initial function as g(x)=x(π-x), the source term as
(2.10)S(x,t)=e-t[x(x-π)+v(π-2x)-D(πΓ(2)Γ(2-α)x1-α-Γ(3)Γ(3-α)x2-α)],
and utilize the definition of Grünwald-Letnikov fractional derivative; the α-order fractional derivatives of the functions x and x2 are given as
(2.11)∂αx∂xα=Γ(2)Γ(2-α)x1-α,∂α(x2)∂xα=Γ(3)Γ(3-α)x2-α,
respectively, then the exact solution of the froward problem can be obtained which is given as c(x,t)=x(π-x)e-t.In the numerical simulation, we will takeα=1.6 as example, and set the diffusion coefficient D=0.1, the velocity v=1; numerical results are listed in Tables 1 and 2, respectively, where M=N denotes the number of space and time grids, respectively, and Err denotes relative error in the exact and numerical solutions at the given time. Moreover, the exact and numerical solutions at t=0.4 and at t=0.8 for M=N=40 are plotted in Figure 1, respectively.Table 1
Errors in the solutions att=0.4 with number of grids.
M
(
N
)
20
40
80
100
200
400
Err
0.0088
0.0055
0.0031
0.0026
0.0014
6.97
e
-
4Table 2
Errors in the solutions att=0.8 with number of grids.
M
(
N
)
20
40
80
100
200
400
Err
0.0165
0.0105
0.0060
0.0049
0.0026
0.0014The exact and numerical solutions withM=N=40.
(a)
t
=
0.4
(b)
t
=
0.8From Tables1 and 2 and Figure 1, we can see that the difference scheme given by (2.9) is of numerical convergence; and the numerical solutions are in good agreement with the exact solution. In what follows, numerical inversions can be performed based on the the above numerical method given by (2.9) by applying the optimal perturbation regularization algorithm.


## 2.1. The Implicit Difference Scheme

For0<x<l and 0<t<T, consider an initial boundary value problem for the FADE given by (1.1), where the initial boundary value conditions are given as
(2.1)c(x,0)=g(x),0≤x≤l,c(0,t)=c(l,t)=0,0≤t≤T,
where the function g(x) is assumed to be continuous on [0,l] and g(0)=g(l)=0.Let the space step beh=l/M and the time step be r=T/N, with xi=ih and tn=nr for i=0,1,…,M and n=0,1,…,N. According to standard Grünwald formula (see [9], e.g.), a discrete approximation to the fractional derivative is given as follows:
(2.2)∂αc(x,t)∂xα=limh→0h=x/M1hα∑k=0M(-1)k(αk)c(x-kh,t),
with which a modified approximation which is called the shifted Grünwald formula is obtained with the help of Gamma function (see [10, 12], e.g.)
(2.3)∂αc(x,t)∂xα=1hα∑k=0MΓ(k-α)Γ(-α)Γ(k+1)c(x-(k-1)h,t)+O(h).
Then we have a difference scheme by (1.1) with the initial boundary value conditions (2.1)
(2.4)cin+1-cinr=-vcin+1-ci-1n+1h+D1hα∑k=0i+1Γ(k-α)Γ(-α)Γ(k+1)ci-k+1n+1+Sin+1,c0n=0,cMn=0,ci0=g(xi),
where cin≈c(xi,tn), Sin≈S(xi,tn) for i=0,1,…,M and n=0,1,…,N-1.Denote(2.5)gk=Γ(k-α)Γ(-α)Γ(k+1),
and p=vr/h and q=Dr/hα, then this can be rearranged to yield
(2.6)-qg0ci+1n+1+(1+p-qg1)cin+1-(p+qg2)ci-1n+1-q∑k=3i+1gkci-k+1n+1=cin+rSin+1,
for i=1,2,…,M-1 and n=0,1,…,N-1. Let cn=(c1n,c2n,…,cM-1n)T, Sn=(S1n,S2n,…,SM-1n)T, and A=(aij) is the coefficient matrix given by (2.6), where the coefficients aij are defined by
(2.7)aij={0,j>i+1,-qg0,j=i+1,1+p-qg1,j=i,-p-qg2,j=i-1,-qgi-j+1,j<i-1,
for i=1,2,…,M-2 and j=1,2,…,M-1, and
(2.8)aM-1,j={-qgM-j,j=1,2,…,M-3,-p-qg2,j=M-2,1+p-qg1,j=M-1,
for i=M-1, respectively. Thus, we get the implicit finite difference scheme in matrix form given as follows:
(2.9)Acn+1=cn+rSn+1,n=0,1,…,N-1.For the above difference scheme, similarly as done in [10, 12], it is not difficult to get its convergence and stability with the help of spectrum analysis to the coefficient matrix.Lemma 2.1 (see [10, 12]).
The implicit difference scheme defined by (2.9) is unconditionally stable, and the difference solution is convergent to the exact solution of the forward problem as h,r→0 for finite time domain.


## 2.2. Numerical Testification

In this subsection, an example is presented to show effectiveness of the finite difference scheme (2.9) for solving the forward problem numerically. For the forward problem (1.1)-(2.1), set l=π, T=1, and take the initial function as g(x)=x(π-x), the source term as
(2.10)S(x,t)=e-t[x(x-π)+v(π-2x)-D(πΓ(2)Γ(2-α)x1-α-Γ(3)Γ(3-α)x2-α)],
and utilize the definition of Grünwald-Letnikov fractional derivative; the α-order fractional derivatives of the functions x and x2 are given as
(2.11)∂αx∂xα=Γ(2)Γ(2-α)x1-α,∂α(x2)∂xα=Γ(3)Γ(3-α)x2-α,
respectively, then the exact solution of the froward problem can be obtained which is given as c(x,t)=x(π-x)e-t.In the numerical simulation, we will takeα=1.6 as example, and set the diffusion coefficient D=0.1, the velocity v=1; numerical results are listed in Tables 1 and 2, respectively, where M=N denotes the number of space and time grids, respectively, and Err denotes relative error in the exact and numerical solutions at the given time. Moreover, the exact and numerical solutions at t=0.4 and at t=0.8 for M=N=40 are plotted in Figure 1, respectively.Table 1
Errors in the solutions att=0.4 with number of grids.
M
(
N
)
20
40
80
100
200
400
Err
0.0088
0.0055
0.0031
0.0026
0.0014
6.97
e
-
4Table 2
Errors in the solutions att=0.8 with number of grids.
M
(
N
)
20
40
80
100
200
400
Err
0.0165
0.0105
0.0060
0.0049
0.0026
0.0014The exact and numerical solutions withM=N=40.
(a)
t
=
0.4
(b)
t
=
0.8From Tables1 and 2 and Figure 1, we can see that the difference scheme given by (2.9) is of numerical convergence; and the numerical solutions are in good agreement with the exact solution. In what follows, numerical inversions can be performed based on the the above numerical method given by (2.9) by applying the optimal perturbation regularization algorithm.


## 3. The Inverse Problem and the Inversion Algorithm

In many cases for the anomalous diffusion model given by (1.1), the fractional order, the diffusion coefficient, and other physical quantities are often unknown and cannot be measured easily. However, these physical quantities can be identified from some known data which can be observed or measured easily based on the forward problem. The method of obtaining these quantities is considered as an inverse problem. Suppose that the fractional order α, the diffusion coefficient D, and the flow velocity v in (1.1) are all unknown, and we have final observation as the additional information given as
(3.1)c(x,T)=θ(x),0≤x≤l.
Thus, an inverse problem of simultaneously determining these three parameters α, D, and v is formulated by (1.1) and the initial boundary value conditions (2.1) with the additional condition (3.1). In what follows, we are to deal with the above multiparameters identification problem from numerics by utilizing an optimal perturbation algorithm.As we know, most of inversion algorithms are based on regularization strategies so as to overcome ill-posedness of the inverse problem, and different kinds of inverse problems could need different approximate methods on the basis of conditional well-posedness analysis (see [25, 30–33], e.g.). It is noted that the optimal perturbation algorithm has been attested to be effective for identifying model parameters at least for the ordinary integer-order diffusion equations (see [26–29, 34], e.g.), and it is also efficient for determining source functions in the FADE (see [23, 24], e.g.). In this section, we will also employ it to determine the three parameters α, D, and v simultaneously in the FADE. The inversion algorithm is given as follows.Denotea=(α,D,v)=(a1,a2,a3), and suppose that c(x,t;a) is the unique solution of the forward problem for any prescribed a∈R3; then a feasible way to solve the inverse problem here is to solve the minimization problem
(3.2)mina∈R3{∥c(x,T;a)-θ(x)∥22+μ∥a∥22},
where μ>0 is the regularization parameter, c(x,T;a) is the computational output obtained by the solution taking values at t=T, and θ(x) is the final observation referred to (3.1).Suppose that(3.3)aj+1=aj+δaj,j=0,1,…,
here δaj is called a perturbation for given aj. Thus, in order to get aj+1 from the given aj, we only need to get an optimal perturbation δaj. In what follows, for convenience of writing, aj and δaj are abbreviated as a and δa, respectively. Then, we only need to work out an perturbation vector
(3.4)δa=(δa1,δa2,δa3).Taking Taylor’s expansion forc(x,T;a+δa) at a and ignoring higher order terms, we can get
(3.5)c(x,T;a+δa)≈c(x,T;a)+∇Tc(x,T;a)·δa.
Noting to (3.2), and with the help of (3.5), we define an error functional given as follows:
(3.6)F(δa)=∥∇Tc(x,T;a)·δa-[θ(x)-c(x,T;a)]∥22+μ∥δa∥22.
Now, discretizing the space domain [0,l] with 0=x1<x2<⋯<xK=l, where K denotes the number of grids, the above L2 norm can then be reduced to the discrete Euclidean norm given as
(3.7)F(δa)=∥Bδa-(η-ξ)∥22+μ∥δa∥22,
where
(3.8)B=(bks)K×3,bks=c(xk,T;a+τes)-c(xk,T;a)τ,s=1,2,3,
and τ is the numerical differentiation step, and e1=(1,0,0), e2=(0,1,0), and e3=(0,0,1), and
(3.9)ξ=(c(x1,T;a),c(x2,T;a),…,c(xK,T;a)),η=(θ(x1),θ(x2),…,θ(xK)).It is not difficult to testify that minimizing functional (3.7) is equivalent to the solving of the following normal equation (see [35], e.g.)
(3.10)(μI+BTB)δa=BT(η-ξ).
So, an optimal perturbation can be solved via the formula:
(3.11)δa=(μI+BTB)-1BT(η-ξ),
with which an optimal inversion solution can be approximated by the iteration procedure (3.3) as long as the perturbation is satisfying a given precision. Obviously, it lies in suitable choices of the regularization parameter, the numerical differential step, the initial iteration, and the convergent precision to perform the above inversion algorithm. In the next section, numerical inversions will be performed by employing the above inversion algorithm.


## 4. Numerical Inversion

Taking the space domain as[0,π], the final time as T=1, the source term as S(x,t)=xe-t, and the initial function as g(x)=x2(π-x), we will carry out numerical simulations in this section. For the given exact parameters α, D, and v in (1.1), the forward problem is solved by the difference scheme (2.9) with M=100, N=20, and then the additional data are obtained with which the inversion algorithm is performed to reconstruct the exact parameters. All computations are finished in a PC of Tsinghua Tongfang.


### 4.1. Influence of Regularization Parameter on the Algorithm

It is important to choose an optimal regularization parameter by theoretical analysis in solving inverse problems, where regularization methods are needed. However, it is still a feasible approach to the choice of regularization parameter by trial and error especially for moderate ill-posed inverse problems. The regularization parameter here is selected by numerical testification, and its influence on the inversion algorithm for given differential step is discussed in this subsection.Suppose that the fractional order isα=1.6, the diffusion coefficient is D=1, and the average velocity is v=1, that is, a=(1.6,1,1) is regarded as the exact solution in this subsection. In the concrete implementation of the inversion algorithm, set initial iteration as a0=(1.1,0,0), numerical differential step as τ=0.1, and convergent criterion as ∥δa∥2≤1e-4. The computational results are listed in Table 3, where μ denotes the regularization parameter, ainv denotes the inversion solution, Tcpu/I is the CPU time of each iteration, Tcpu is the total CPU time for each inversion, j denotes the number of iterations, and Err=∥ainv-a∥2/∥a∥2 is the relative error in the solutions.Table 3
Influence of regularization parameter on the algorithm withτ=0.1.
μ
a
inv
Err
T
cpu
/
j
1
(1.5998012, 1.0004052, 1.0003096)
2.5630155
e
-
4
10.190/59
8
e
-
1
(1.5998300, 1.0003477, 1.0002653)
2.1975589
e
-
4
9.516/52
6
e
-
1
(1.5998732, 1.0002609, 1.0001987)
1.6463489
e
-
4
8.281/45
5
e
-
1
(1.5998862, 1.0002350, 1.0001787)
1.4814645
e
-
4
7.156/41
3
e
-
1
(1.5999263, 1.0001542, 1.0001167)
9.6913294
e
-
5
6.000/33
1
e
-
1
(1.5999617, 1.0000825, 1.0000618)
5.1475320
e
-
5
4.250/24
1
e
-
2
(1.5999665, 1.0000747, 1.0000552)
4.6246166
e
-
5
3.391/19
1
e
-
4
(1.5999779, 1.0000495, 1.0000365)
3.0598630
e
-
5
3.375/19
1
e
-
6
(1.5999780, 1.0000493, 1.0000364)
3.0467546
e
-
5
3.406/19
1
e
-
8
(1.5999780, 1.0000492, 1.0000364)
3.0466238
e
-
5
3.344/19
1
e
-
10
(1.5999780, 1.0000492, 1.0000364)
3.0466225
e
-
5
3.157/19
0
DivergentIt is noticeable that the inversion results have little changes if utilizing any positive regularization parameters smaller than1e-10, and the solutions error still remain Err=3.0466225e-5 even with very small regularization parameter μ=1e-36. The inversion errors varying with regularization parameters in this case are plotted in Figure 2(a). Furthermore, if taking numerical differential step as τ=0.01, and other inversion parameters unchanged, the computational results are listed in Table 4, where μ, ainv, Err, and Tcpu/j all denote the same meanings as used in Table 3.Table 4
Influence of regularization parameter on the algorithm withτ=0.01.
μ
a
inv
Err
T
cpu
/
j
5
(1.6015151, 0.99697413, 0.99767406)
1.9229322
e
-
3
32.657/176
4.5
(1.6013808, 0.99724171, 0.99787964)
1.7528679
e
-
3
29.687/163
4
(1.6012283, 0.99754585, 0.99811333)
1.5595733
e
-
3
27.704/150
3.5
(1.6010527, 0.99789612, 0.99838251)
1.3369560
e
-
3
24.703/137
3
(1.6009217, 0.99815749, 0.99858340)
1.1708440
e
-
3
22.500/122
2.5
(1.6007547, 0.99849085, 0.99883965)
9.5898118
e
-
4
20.094/108
2.4
(1.6007389, 0.99852247, 0.99886397)
9.3887868
e
-
4
19.313/105
2.3
(1.6006934, 0.99861329, 0.99893378)
8.8116322
e
-
4
18.985/103
2.2
(1.6006824, 0.99863531, 0.99895072)
8.6716578
e
-
4
18.469/100
2.1
(1.6006235, 0.99875285, 0.99904107)
7.9247157
e
-
4
17.640/96
≤2
DivergentThe solutions errors with regularization parameter and differential step.
(a)
τ
=
0.1
(b)
μ
=
0.1From Tables3 and 4 and Figure 2(a), we can find that regularization parameter has important impact on the inversion algorithm for given numerical differential step, and the regularization parameter must be chosen large for small numerical differential step. However, the inversion algorithm can be performed smoothly with any small regularization parameters strictly greater than zero if choosing large numerical differential steps.


### 4.2. Influence of Numerical Differential Step on the Algorithm

By the above computations together with (3.8), we find that numerical differentiation is another key point in the realization of the inversion algorithm, and suitable differential steps are necessary in order to realize a successful inversion. As done for the regularization parameter, we will investigate the influence of the differential step on the algorithm by testification.Also set the exact solution to bea=(1.6,1,1), and take the regularization parameter as μ=0.1, and other inversion parameters also unchanged as used in the above. The inversion results are listed in Table 5, where τ denotes the numerical differential step, and ainv, Err, and Tcpu/j all denote the same meanings as referred in the above. Moreover, the errors in the solutions varying with the numerical differential steps for μ=0.1 are plotted in Figure 2(b).Table 5
Influence of numerical differential step on the algorithm withμ=0.1.
τ
a
inv
Err
T
cpu
/
j
≥0.9
Divergent
0.85
(1.6001245, 0.99973957, 0.99978970)
1.6724819
e
-
4
6.563/36
0.8
(1.6001340, 0.99969424, 0.99977646)
1.8814739
e
-
4
6.453/35
0.7
(1.6001260, 0.99971463, 0.99979024)
1.7604399
e
-
4
6.172/34
0.6
(1.6001232, 0.99972209, 0.99979564)
1.7153187
e
-
4
5.937/32
0.5
(1.6001146, 0.99974254, 0.99981050)
1.5903046
e
-
4
5.359/29
0.4
(1.6000995, 0.99977749, 0.99983595)
1.3759274
e
-
4
3.672/20
0.3
(1.5999046, 1.0002119, 1.0001567)
1.3124502
e
-
4
4.984/27
0.2
(1.5999303, 1.0001532, 1.0001138)
9.5146183
e
-
5
5.203/28
0.1
(1.5999617, 1.0000825, 1.0000618)
5.1475320
e
-
5
4.375/24
0.05
(1.5999773, 1.0000477, 1.0000361)
2.9959973
e
-
5
4.484/24
0.04
(1.5999772, 1.0000475, 1.0000360)
2.9877074
e
-
5
5.344/29
≤0.03
DivergentFrom Table5 and Figure 2(b), we can see that the numerical differential step plays a similar role as to that of the regularization parameter in the realization of the algorithm. The inversion algorithm could be a failure in case of using large or small differential steps. The differential steps should fall in the interval of (0.03,0.9) in this example.


### 4.3. Influence of Fractional Order on the Algorithm

In this subsection, we will investigate influence of the fractional order on the inversion algorithm with the convergent criterion given as∥δa∥2≤10-6. By the above computations, the regularization parameter here is μ=0.1, the differential step is τ=0.1, and the initial iteration is a0=(1.1,0,0) as before. The computational results are listed in Table 6, where a denotes the exact solution, α denotes the fractional order, and ainv, Err, and Tcpu/j also refer to the same notations as used in the above.Table 6
Influence of fractional order on the inversion algorithm.
α
a
a
inv
Err
T
cpu
/
j
1.9
(1.9, 1, 1)
(1.90000, 1.00000, 1.00000)
1.7848
e
-
7
3.547/21
1.8
(1.8, 1, 1)
(1.80000, 1.00000, 1.00000)
1.6868
e
-
7
3.859/23
1.7
(1.7, 1, 1)
(1.70000, 1.00000, 1.00000)
3.4783
e
-
7
4.344/26
1.6
(1.6, 1, 1)
(1.60000, 1.00000, 1.00000)
4.9312
e
-
7
5.500/33
1.5
(1.5, 1, 1)
(1.50000, 1.00000, 1.00000)
9.7401
e
-
7
7.797/47
1.4
(1.4, 1, 1)
(1.40000, 1.00000, 1.00000)
2.6188
e
-
6
13.031/79
1.3
(1.3, 1, 1)
(1.30000, 1.00001, 1.00001)
9.5267
e
-
6
24.437/148
1.2
(1.2, 1, 1)
(1.24150, 0.838272, 0.845473)
0.1227
87.734/531From Table6, we find that the inversion results are not so good when the fractional order is smaller than 1.2, which demonstrates that ill-posedness of the inversion problem could become severe if coping with small fractional orders.


### 4.4. Inversion with Noisy Data

It is difficult to perform an inversion algorithm in the case of using random noisy data, especially for inverse problems arising from the fractional diffusion. Noting computational errors and data noises, the additional information utilized for real inverse problems is often given as(4.1)θε(x)=θ(x)+εζ,
where θ(x) is the accurate additional information given by (3.1), ε>0 is noise level, and ζ is a random vector ranged in [-1,1]. In the following, we will take the fractional order α=1.5 as example, that is, the exact solution in this subsection is given as a=(1.5,1,1).Generally speaking, in case of using regularization strategy to damp the noises, an optimal regularization parameter should be chosen according to the noise level. However, the situation becomes simple for the inverse problem considered here. Without loss of generality, we set the noise levels beε=1% and ε=5% and we choose regularization parameter as μ=0.5, μ=0.1 and μ=0.01, respectively, and other inversion parameters also unchanged as in the above. The average inversion results by ten-time computations are listed in Tables 7 and 8 respectively, where ε denotes the noise level, a-inv denotes the average inversion solution of the ten-time computations, E-rr=∥a-inv-a∥2/∥a∥2 denotes the average relative error in the solutions, and T-cpu with j- denote the average CPU time and iteration number of the inversion, respectively.Table 7
Inversion results with regularization parameters forε=1%.
μ
a
-
inv
E
-
rr
T
-
cpu
j
-
0.01
(1.4992660, 1.0017420, 1.0015202)
1.1766681
e
-
3
7.2
39.9
0.1
(1.4992650, 1.0017420, 1.0015203)
1.1768453
e
-
3
8.7
47
0.5
(1.4992640, 1.0017427, 1.0015217)
1.1776617
e
-
3
13.8
75.1Table 8
Inversion results with regularization parameters forε=5%.
μ
a
-
inv
E
-
rr
T
-
cpu
j
-
0.01
(1.4947065, 1.0158355, 1.0141558)
1.0618188
e
-
2
8.4
42.0
0.1
(1.4947063, 1.0158360, 1.0141562)
1.0618514
e
-
2
8.3
49.1
0.5
(1.4947055, 1.0158378, 1.0141577)
1.0619698
e
-
2
13.1
79.0From Tables7 and 8, we can see clearly that for the given noise level, the inversion solutions are almost unchanged, although the number of iterations and the corresponding CPU time increase to some extent when using large regularization parameters. Furthermore, Table 9 gives inversion results by using regularization parameter as μ=0.1 with different noise levels, where ε also denotes the noise level and a-inv, E-rr, and T-cpu with j- denote the same meanings as used in the above.Table 9
Inversion results with random noises usingμ=0.1.
ε
a
-
inv
E
-
rr
T
-
cpu
j
-
10
%
(1.4841, 1.0568, 1.0515)
0.0380
11.9
58.3
5
%
(1.4947, 1.0158, 1.0142)
0.0106
8.3
49.1
3
%
(1.4973, 1.0073, 1.0065)
0.0049
8.8
47.4
1
%
(1.4993, 1.0017, 1.0015)
0.0012
8.7
47.0
0.5
%
(1.4997, 1.0008, 1.0007)
5.33
e
-
4
8.4
46.8
0.1
%
(1.4999, 1.0001, 1.0001)
9.87
e
-
5
8.1
46.9
0.01
%
(1.5000, 1.0000, 1.0000)
1.03
e
-
5
8.1
47.0From Table9, it can be seen that the inversion results are satisfactory in the case of using noisy data, and the error in the solutions becomes small and approaches to zero with the noise level goes to zero demonstrating that the inversion algorithm is of numerical stability.


## 4.1. Influence of Regularization Parameter on the Algorithm

It is important to choose an optimal regularization parameter by theoretical analysis in solving inverse problems, where regularization methods are needed. However, it is still a feasible approach to the choice of regularization parameter by trial and error especially for moderate ill-posed inverse problems. The regularization parameter here is selected by numerical testification, and its influence on the inversion algorithm for given differential step is discussed in this subsection.Suppose that the fractional order isα=1.6, the diffusion coefficient is D=1, and the average velocity is v=1, that is, a=(1.6,1,1) is regarded as the exact solution in this subsection. In the concrete implementation of the inversion algorithm, set initial iteration as a0=(1.1,0,0), numerical differential step as τ=0.1, and convergent criterion as ∥δa∥2≤1e-4. The computational results are listed in Table 3, where μ denotes the regularization parameter, ainv denotes the inversion solution, Tcpu/I is the CPU time of each iteration, Tcpu is the total CPU time for each inversion, j denotes the number of iterations, and Err=∥ainv-a∥2/∥a∥2 is the relative error in the solutions.Table 3
Influence of regularization parameter on the algorithm withτ=0.1.
μ
a
inv
Err
T
cpu
/
j
1
(1.5998012, 1.0004052, 1.0003096)
2.5630155
e
-
4
10.190/59
8
e
-
1
(1.5998300, 1.0003477, 1.0002653)
2.1975589
e
-
4
9.516/52
6
e
-
1
(1.5998732, 1.0002609, 1.0001987)
1.6463489
e
-
4
8.281/45
5
e
-
1
(1.5998862, 1.0002350, 1.0001787)
1.4814645
e
-
4
7.156/41
3
e
-
1
(1.5999263, 1.0001542, 1.0001167)
9.6913294
e
-
5
6.000/33
1
e
-
1
(1.5999617, 1.0000825, 1.0000618)
5.1475320
e
-
5
4.250/24
1
e
-
2
(1.5999665, 1.0000747, 1.0000552)
4.6246166
e
-
5
3.391/19
1
e
-
4
(1.5999779, 1.0000495, 1.0000365)
3.0598630
e
-
5
3.375/19
1
e
-
6
(1.5999780, 1.0000493, 1.0000364)
3.0467546
e
-
5
3.406/19
1
e
-
8
(1.5999780, 1.0000492, 1.0000364)
3.0466238
e
-
5
3.344/19
1
e
-
10
(1.5999780, 1.0000492, 1.0000364)
3.0466225
e
-
5
3.157/19
0
DivergentIt is noticeable that the inversion results have little changes if utilizing any positive regularization parameters smaller than1e-10, and the solutions error still remain Err=3.0466225e-5 even with very small regularization parameter μ=1e-36. The inversion errors varying with regularization parameters in this case are plotted in Figure 2(a). Furthermore, if taking numerical differential step as τ=0.01, and other inversion parameters unchanged, the computational results are listed in Table 4, where μ, ainv, Err, and Tcpu/j all denote the same meanings as used in Table 3.Table 4
Influence of regularization parameter on the algorithm withτ=0.01.
μ
a
inv
Err
T
cpu
/
j
5
(1.6015151, 0.99697413, 0.99767406)
1.9229322
e
-
3
32.657/176
4.5
(1.6013808, 0.99724171, 0.99787964)
1.7528679
e
-
3
29.687/163
4
(1.6012283, 0.99754585, 0.99811333)
1.5595733
e
-
3
27.704/150
3.5
(1.6010527, 0.99789612, 0.99838251)
1.3369560
e
-
3
24.703/137
3
(1.6009217, 0.99815749, 0.99858340)
1.1708440
e
-
3
22.500/122
2.5
(1.6007547, 0.99849085, 0.99883965)
9.5898118
e
-
4
20.094/108
2.4
(1.6007389, 0.99852247, 0.99886397)
9.3887868
e
-
4
19.313/105
2.3
(1.6006934, 0.99861329, 0.99893378)
8.8116322
e
-
4
18.985/103
2.2
(1.6006824, 0.99863531, 0.99895072)
8.6716578
e
-
4
18.469/100
2.1
(1.6006235, 0.99875285, 0.99904107)
7.9247157
e
-
4
17.640/96
≤2
DivergentThe solutions errors with regularization parameter and differential step.
(a)
τ
=
0.1
(b)
μ
=
0.1From Tables3 and 4 and Figure 2(a), we can find that regularization parameter has important impact on the inversion algorithm for given numerical differential step, and the regularization parameter must be chosen large for small numerical differential step. However, the inversion algorithm can be performed smoothly with any small regularization parameters strictly greater than zero if choosing large numerical differential steps.


## 4.2. Influence of Numerical Differential Step on the Algorithm

By the above computations together with (3.8), we find that numerical differentiation is another key point in the realization of the inversion algorithm, and suitable differential steps are necessary in order to realize a successful inversion. As done for the regularization parameter, we will investigate the influence of the differential step on the algorithm by testification.Also set the exact solution to bea=(1.6,1,1), and take the regularization parameter as μ=0.1, and other inversion parameters also unchanged as used in the above. The inversion results are listed in Table 5, where τ denotes the numerical differential step, and ainv, Err, and Tcpu/j all denote the same meanings as referred in the above. Moreover, the errors in the solutions varying with the numerical differential steps for μ=0.1 are plotted in Figure 2(b).Table 5
Influence of numerical differential step on the algorithm withμ=0.1.
τ
a
inv
Err
T
cpu
/
j
≥0.9
Divergent
0.85
(1.6001245, 0.99973957, 0.99978970)
1.6724819
e
-
4
6.563/36
0.8
(1.6001340, 0.99969424, 0.99977646)
1.8814739
e
-
4
6.453/35
0.7
(1.6001260, 0.99971463, 0.99979024)
1.7604399
e
-
4
6.172/34
0.6
(1.6001232, 0.99972209, 0.99979564)
1.7153187
e
-
4
5.937/32
0.5
(1.6001146, 0.99974254, 0.99981050)
1.5903046
e
-
4
5.359/29
0.4
(1.6000995, 0.99977749, 0.99983595)
1.3759274
e
-
4
3.672/20
0.3
(1.5999046, 1.0002119, 1.0001567)
1.3124502
e
-
4
4.984/27
0.2
(1.5999303, 1.0001532, 1.0001138)
9.5146183
e
-
5
5.203/28
0.1
(1.5999617, 1.0000825, 1.0000618)
5.1475320
e
-
5
4.375/24
0.05
(1.5999773, 1.0000477, 1.0000361)
2.9959973
e
-
5
4.484/24
0.04
(1.5999772, 1.0000475, 1.0000360)
2.9877074
e
-
5
5.344/29
≤0.03
DivergentFrom Table5 and Figure 2(b), we can see that the numerical differential step plays a similar role as to that of the regularization parameter in the realization of the algorithm. The inversion algorithm could be a failure in case of using large or small differential steps. The differential steps should fall in the interval of (0.03,0.9) in this example.


## 4.3. Influence of Fractional Order on the Algorithm

In this subsection, we will investigate influence of the fractional order on the inversion algorithm with the convergent criterion given as∥δa∥2≤10-6. By the above computations, the regularization parameter here is μ=0.1, the differential step is τ=0.1, and the initial iteration is a0=(1.1,0,0) as before. The computational results are listed in Table 6, where a denotes the exact solution, α denotes the fractional order, and ainv, Err, and Tcpu/j also refer to the same notations as used in the above.Table 6
Influence of fractional order on the inversion algorithm.
α
a
a
inv
Err
T
cpu
/
j
1.9
(1.9, 1, 1)
(1.90000, 1.00000, 1.00000)
1.7848
e
-
7
3.547/21
1.8
(1.8, 1, 1)
(1.80000, 1.00000, 1.00000)
1.6868
e
-
7
3.859/23
1.7
(1.7, 1, 1)
(1.70000, 1.00000, 1.00000)
3.4783
e
-
7
4.344/26
1.6
(1.6, 1, 1)
(1.60000, 1.00000, 1.00000)
4.9312
e
-
7
5.500/33
1.5
(1.5, 1, 1)
(1.50000, 1.00000, 1.00000)
9.7401
e
-
7
7.797/47
1.4
(1.4, 1, 1)
(1.40000, 1.00000, 1.00000)
2.6188
e
-
6
13.031/79
1.3
(1.3, 1, 1)
(1.30000, 1.00001, 1.00001)
9.5267
e
-
6
24.437/148
1.2
(1.2, 1, 1)
(1.24150, 0.838272, 0.845473)
0.1227
87.734/531From Table6, we find that the inversion results are not so good when the fractional order is smaller than 1.2, which demonstrates that ill-posedness of the inversion problem could become severe if coping with small fractional orders.


## 4.4. Inversion with Noisy Data

It is difficult to perform an inversion algorithm in the case of using random noisy data, especially for inverse problems arising from the fractional diffusion. Noting computational errors and data noises, the additional information utilized for real inverse problems is often given as(4.1)θε(x)=θ(x)+εζ,
where θ(x) is the accurate additional information given by (3.1), ε>0 is noise level, and ζ is a random vector ranged in [-1,1]. In the following, we will take the fractional order α=1.5 as example, that is, the exact solution in this subsection is given as a=(1.5,1,1).Generally speaking, in case of using regularization strategy to damp the noises, an optimal regularization parameter should be chosen according to the noise level. However, the situation becomes simple for the inverse problem considered here. Without loss of generality, we set the noise levels beε=1% and ε=5% and we choose regularization parameter as μ=0.5, μ=0.1 and μ=0.01, respectively, and other inversion parameters also unchanged as in the above. The average inversion results by ten-time computations are listed in Tables 7 and 8 respectively, where ε denotes the noise level, a-inv denotes the average inversion solution of the ten-time computations, E-rr=∥a-inv-a∥2/∥a∥2 denotes the average relative error in the solutions, and T-cpu with j- denote the average CPU time and iteration number of the inversion, respectively.Table 7
Inversion results with regularization parameters forε=1%.
μ
a
-
inv
E
-
rr
T
-
cpu
j
-
0.01
(1.4992660, 1.0017420, 1.0015202)
1.1766681
e
-
3
7.2
39.9
0.1
(1.4992650, 1.0017420, 1.0015203)
1.1768453
e
-
3
8.7
47
0.5
(1.4992640, 1.0017427, 1.0015217)
1.1776617
e
-
3
13.8
75.1Table 8
Inversion results with regularization parameters forε=5%.
μ
a
-
inv
E
-
rr
T
-
cpu
j
-
0.01
(1.4947065, 1.0158355, 1.0141558)
1.0618188
e
-
2
8.4
42.0
0.1
(1.4947063, 1.0158360, 1.0141562)
1.0618514
e
-
2
8.3
49.1
0.5
(1.4947055, 1.0158378, 1.0141577)
1.0619698
e
-
2
13.1
79.0From Tables7 and 8, we can see clearly that for the given noise level, the inversion solutions are almost unchanged, although the number of iterations and the corresponding CPU time increase to some extent when using large regularization parameters. Furthermore, Table 9 gives inversion results by using regularization parameter as μ=0.1 with different noise levels, where ε also denotes the noise level and a-inv, E-rr, and T-cpu with j- denote the same meanings as used in the above.Table 9
Inversion results with random noises usingμ=0.1.
ε
a
-
inv
E
-
rr
T
-
cpu
j
-
10
%
(1.4841, 1.0568, 1.0515)
0.0380
11.9
58.3
5
%
(1.4947, 1.0158, 1.0142)
0.0106
8.3
49.1
3
%
(1.4973, 1.0073, 1.0065)
0.0049
8.8
47.4
1
%
(1.4993, 1.0017, 1.0015)
0.0012
8.7
47.0
0.5
%
(1.4997, 1.0008, 1.0007)
5.33
e
-
4
8.4
46.8
0.1
%
(1.4999, 1.0001, 1.0001)
9.87
e
-
5
8.1
46.9
0.01
%
(1.5000, 1.0000, 1.0000)
1.03
e
-
5
8.1
47.0From Table9, it can be seen that the inversion results are satisfactory in the case of using noisy data, and the error in the solutions becomes small and approaches to zero with the noise level goes to zero demonstrating that the inversion algorithm is of numerical stability.


## 5. Conclusion

(i) An inverse problem for identifying multiparameters of the fractional order, the diffusion coefficient, and the average velocity simultaneously in the FADE with final observations is investigated. An implicit finite difference scheme is employed to solve the forward problem, and an optimal perturbation regularization algorithm is applied to reconstruct the three parameters. By the numerical simulations, we conclude that the optimal perturbation regularization algorithm is of numerical stability, and it is efficient for multiparameters identification problem arising from the FADE.(ii) There are few factors that have influences on the inversion algorithm, but the regularization parameter and the numerical differential step here seem to play a much more important role in the realization of the algorithm. The inversion algorithm can be performed successfully for the regularization parameter and the numerical differential step belonging to suitable intervals, respectively. However, the inversion algorithm may be a failure if the numerical differential step is taking too small values. In addition, it seems to be insensitive to the choice of the regularization parameter in concrete realization of the algorithm showing that the inverse problem discussed here is of moderate ill-posedness.(iii) By the inversion computations, we find that the fractional order in the FADE has some influence on the forward problem and the inverse problem. If the fractional orderα goes to 2, the solutions error becomes small and the computational complexity is reduced; however, if α is less than 1.2, the ill-posedness of the inverse problem becomes severe, and the inversion results could be unacceptable. We will deal with the forward problem and the corresponding inverse problem with the fractional order smaller than 1.2 in our sequent works.



---

*Source: 740385-2012-12-24.xml*