# Numerical Inversion for the Multiple Fractional Orders in the Multiterm TFDE

Authors: Chunlong Sun; Gongsheng Li; Xianzheng Jia
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

The fractional order in a fractional diffusion model is a key parameter which characterizes the anomalous diffusion behaviors. This paper deals with an inverse problem of determining the multiple fractional orders in the multiterm time-fractional diffusion equation (TFDE for short) from numerics. The homotopy regularization algorithm is applied to solve the inversion problem using the finite data at one interior point in the space domain. The inversion fractional orders with random noisy data give good approximations to the exact order demonstrating the efficiency of the inversion algorithm and numerical stability of the inversion problem.


---

## Body

## 1. Introduction

The partial differential equations of fractional order have played an important role in modeling of the anomalous phenomena and in the theory of the complex systems during the last two decades; see, for example, [1–8]. The so-called time-fractional diffusion equation (TFDE) that is obtained from the classical diffusion equation by replacing the first-order time derivative by a fractional derivative of order α with 0<α<1 has to be especially mentioned. On the other hand, by the attempts to describe some real processes with the equations of the fractional order, several researches were confronted with the situation that the order α of the time-fractional derivative from the corresponding models did not remain constant and changed, say, in the interval from 0 to 1, from 1 to 2, or even from 0 to 2. To manage these phenomena, several approaches were suggested. One of them introduces the fractional derivatives of the variable order, that is, the derivatives with the order that can change with the time or/and depending on the spatial coordinates [9–11], and the other way is to employ the distributed order TFDE, or the multiterm TFDE in discretization. Let Ω be a bounded domain in Rd (d≥1) with smooth boundary ∂Ω, and let T>0; the multiterm homogeneous TFDE with variable coefficient in Ω is given as (1)∂αu∂tα+∑s=1Srs∂βsu∂tβs=∇·D∇u,x∈Ω,0<t≤T,where u=u(x,t) denotes the state variable at space point x and time t, α denotes the principal fractional order, and β1,β2,…,βS are the multiterm fractional orders of the time derivatives, which satisfy the condition (2)0<βS<βS-1<⋯<β1<α<1,and r1,r2,…,rS are positive constants, and D(x)>0 is the smooth diffusion coefficient tensor. All of the above time-fractional derivatives are defined in the sense of Caputo; for example, the fractional derivative of the order β∈(0,1) is given by (3)∂βu∂tβ=1Γ1-β∫0t∂ux,s∂sdst-sβ.See, for example, Podlubny [12] and Kilbas et al. [13] for the definition and properties of Caputo’s derivative.There are still a few research works reported on the multiterm TFDE like (1). On theoretical analysis and analytical methods for the forward problem, we refer to Daftardar-Gejji and Bhalekar [14], Luchko [15, 16], Jiang et al. [17], Ding et al. [18, 19], and Li et al. [20], and for numerical methods and simulations we refer to [21–23], and so on.However, for real problems, the fractional orders, the initial distribution, the diffusion coefficient, or the source term cannot be obtained directly and we have to determine them by some additional measurements, which contributes to inverse problems arising in the fractional diffusion models. There are still some researches on inverse problems for the one-term TFDE; see, for example, Murio [24], Liu et al. [25, 26], Sakamoto and Yamamoto [27], Tuan [28], Chi et al. [29], Yamamoto and Zhang [30], Luchko et al. [31], Wei et al. [32, 33], and Liu et al. [34]; also see Jin and Rundell [35] for a tutorial review on inverse problems for anomalous diffusion processes.It is noted that the research works stated above are almost related to coefficient identification problems in the one-term time/space fractional diffusion equations. However, it is also important to deal with inverse problems of determining the fractional orders in the fractional differential equations since the fractional order is an essential index characterizing the anomalous diffusion. As for inverse problems of determining fractional orders in the single-term time/space fractional diffusion models, we refer to [36–42], and so on. On the other hand, there are few literatures concerned with the inverse problems in the multiterm TFDEs to our knowledge. Li and Yamamoto [39] studied an inverse problem of identifying the multiple fractional orders in the multiterm TFDE, and they gave the uniqueness result using Laplace transform and analytical method, and later they considered the similar model [42], and also the uniqueness of determining the fractional orders, the number of the fractional terms, and the spatially varying coefficient simultaneously is proved. Recently, Sun et al. [43] considered a simultaneous inversion problem for determining the space-dependent diffusion and source coefficients in the multiterm TFDE using the optimal perturbation regularization algorithm, and quite a few numerical inversions are presented.Based on the above analysis, we are to deal with the inverse problem of determining the multiple fractional orders in the multiterm TFDE with the additional measurements at the interior point from numerics. The uniqueness results for such kind of inverse problems have been obtained (c.f. [39, 42], e.g.), but numerical inversions are still open to be implemented. Based on the difference solution to the forward problem, we perform numerical inversions by utilizing the homotopy regularization algorithm not only with the accurate data but also with random noisy data. The inversion fractional orders approximate to the exact orders as the noise level gets smaller demonstrating a numerical stability of the inverse problem here.The rest of the paper is organized as follows. In Section2, an implicit finite difference solution to the forward problem is given and the inverse problem of determining the fractional orders is formulated. In Section 3, the homotopy regularization algorithm is introduced to solve the inversion problem and numerical inversions are presented, and concluding remarks are given in Section 4.


## 2. The Forward Problem and the Inverse Problem

Consider the forward problem given by (1) with the initial condition (4)ux,0=u0x,x∈Ω,and the homogeneous Dirichlet boundary condition (5)u∂Ω=0,0<t≤T,where the initial function is smooth enough and satisfies the consistency condition with the boundary condition. Let the diffusion coefficient D(x) be C1-class and take positive values on Ω¯; then the forward problem (1), (4)-(5) has a unique solution for suitable initial functions (c.f. [20], e.g.). Here we focus our attention on the finite difference solution to the forward problem. For completeness of the paper, we give an implicit finite difference scheme in 1D case for solving the forward problem. For further details, see [22, 23], and so on.


### 2.1. The Difference Scheme to the Forward Problem

LetΩ=(0,l) for l>0. For given integer numbers M and N, discretizing the space domain by xi=ih(i=0,1,…,M) and the time domain by tn=nτ(n=0,1,…,N), we have by definition (3) (6)∂αu∂tαxi,tn+1=τ-αΓ2-α∑k=0nuxi,tn+1-k-uxi,tn-kk+11-α-k1-α+Oτ,∂βsu∂tβsxi,tn+1=τ-βsΓ2-βs∑k=0nuxi,tn+1-k-uxi,tn-kk+11-βs-k1-βs+Oτ,for s=1,2,…,S, respectively; here h=l/M is the space mesh step and τ=T/N is the time mesh step.LetD(x)∈C1(Ω¯) and D(x)>0 for x∈Ω¯. By discretizing the term (∂/∂x)(D(x)ux) using the ordinary integer-order difference method and denoting uin≈u(xi,tn) and Di=D(xi), we get(7)uin+1-uin1+∑s=1Sriτα-βsΓ2-αΓ2-βs+∑k=1nuin+1-k-uin-kk+11-α-k1-α+∑s=1Srsk+11-βs-k1-βsτα-βsΓ2-αΓ2-βs=DiταΓ2-αh2ui-1n+1-2uin+1+ui+1n+1+Di-Di-1ταΓ2-αh2uin+1-ui-1n+1+Rn+1,where Rn=O(ταh2+τα+1) is the truncated term. Denoting v=1+∑s=1Srsτα-βs(Γ(2-α)/Γ(2-βs)), dividing by v on two sides of (7), and ignoring the truncated term, we get the following: (8)-pi+qiui-1n+1+1+2pi-qiuin+1-piui+1n+1=uin-∑k=1ndkuin+1-k-uin-k,where (9)pi=DiταΓ2-αh2v,qi=Di-Di-1ταΓ2-αh2v,dk=k+11-α-k1-αv+∑s=1Sk+11-βs-k1-βsvrsτα-βsΓ2-αΓ2-βsfor k=1,2,…,n. The initial boundary value conditions are discretized as (10)ui0=u0xi;u0n=0,uMn=0. Let (11)Un=u1n,u2n,…,uM-1n;U0=u10,u20,…,uM-10,ck=1-d1,k=1,dk-1-dk,k=2,…,n;and B=(bij)(M-1)×(M-1), where bij is defined by (12)bij=0,j>i+1,-pi,j=i+1,1+2pi-qi,j=i,-pi-qi,j=i-1,0,j<i-1,for i=1,2,…,M-1,j=1,2,…,M-1.Thus, we have the implicit finite difference scheme in the matrix form given as(13)BU1=U0;BUn+1=c1Un+c2Un-1+⋯+cnU1+dnU0.Theorem 1.
The implicit difference scheme (13) has only one solution, and it is of unconditional stability and convergence for any finite time T>0.Proof.
By the assumptions for the diffusion coefficientD(x), we have pi>0 and pi-qi>0 for i=1,…,M-1. So the coefficient matrix B of (13) is strictly diagonally dominant; then the difference equation (13) has only one solution. Moreover, thanks to the equality ∑k=1nck+dn=1 and with a similar method as used in [22], we get the unconditional stability and convergence of the difference scheme for any given finite time T>0.


### 2.2. The Inverse Problem

For the forward problem (1), (4)-(5), if the fractional orders α and βs (s=1,…,S) in (1) are unknown, we encounter the inverse problem of determining these multiple fractional orders. Suppose that there are some measured points in the space domain; for example, let x0∈Ω be the measured point, and we have the additional information given as (14)ux0,t=θt,0<t≤T,and the inverse problem is to determine α∈(0,1) and βs∈(0,1)(s=1,…,S) using the overposed condition (14) based on (1) and the initial boundary conditions (4) and (5).As stated in Section1, the general uniqueness results have been proved in [39, 42] by Dr. Li and Professor Yamamoto. We give the following lemma for the completeness of the paper.Lemma 2 (see [39]).
Suppose that the fractional ordersα and βs(s=1,…,S), the diffusion coefficient tensor D(x), and the coefficients rs(s=1,…,S) in (1) satisfy the conditions given in Section 1, and the initial function u0(x) is smooth enough, and u0(x)≥0 and u0(x)≢0 for x∈Ω. Then all the fractional orders can be determined uniquely by the additional data {u(x0,t)=θ(t)} for t∈(0,T).The uniqueness result is very important for inverse problems in theory. By Lemma2 we need to utilize the additional data θ(ti) measured at x0 for i=1,2,…; however, in concrete computations, we find that the numerical inversions can also be performed only employing a few of the additional data. So it is still meaningful to study inverse problems from numerics.


## 2.1. The Difference Scheme to the Forward Problem

LetΩ=(0,l) for l>0. For given integer numbers M and N, discretizing the space domain by xi=ih(i=0,1,…,M) and the time domain by tn=nτ(n=0,1,…,N), we have by definition (3) (6)∂αu∂tαxi,tn+1=τ-αΓ2-α∑k=0nuxi,tn+1-k-uxi,tn-kk+11-α-k1-α+Oτ,∂βsu∂tβsxi,tn+1=τ-βsΓ2-βs∑k=0nuxi,tn+1-k-uxi,tn-kk+11-βs-k1-βs+Oτ,for s=1,2,…,S, respectively; here h=l/M is the space mesh step and τ=T/N is the time mesh step.LetD(x)∈C1(Ω¯) and D(x)>0 for x∈Ω¯. By discretizing the term (∂/∂x)(D(x)ux) using the ordinary integer-order difference method and denoting uin≈u(xi,tn) and Di=D(xi), we get(7)uin+1-uin1+∑s=1Sriτα-βsΓ2-αΓ2-βs+∑k=1nuin+1-k-uin-kk+11-α-k1-α+∑s=1Srsk+11-βs-k1-βsτα-βsΓ2-αΓ2-βs=DiταΓ2-αh2ui-1n+1-2uin+1+ui+1n+1+Di-Di-1ταΓ2-αh2uin+1-ui-1n+1+Rn+1,where Rn=O(ταh2+τα+1) is the truncated term. Denoting v=1+∑s=1Srsτα-βs(Γ(2-α)/Γ(2-βs)), dividing by v on two sides of (7), and ignoring the truncated term, we get the following: (8)-pi+qiui-1n+1+1+2pi-qiuin+1-piui+1n+1=uin-∑k=1ndkuin+1-k-uin-k,where (9)pi=DiταΓ2-αh2v,qi=Di-Di-1ταΓ2-αh2v,dk=k+11-α-k1-αv+∑s=1Sk+11-βs-k1-βsvrsτα-βsΓ2-αΓ2-βsfor k=1,2,…,n. The initial boundary value conditions are discretized as (10)ui0=u0xi;u0n=0,uMn=0. Let (11)Un=u1n,u2n,…,uM-1n;U0=u10,u20,…,uM-10,ck=1-d1,k=1,dk-1-dk,k=2,…,n;and B=(bij)(M-1)×(M-1), where bij is defined by (12)bij=0,j>i+1,-pi,j=i+1,1+2pi-qi,j=i,-pi-qi,j=i-1,0,j<i-1,for i=1,2,…,M-1,j=1,2,…,M-1.Thus, we have the implicit finite difference scheme in the matrix form given as(13)BU1=U0;BUn+1=c1Un+c2Un-1+⋯+cnU1+dnU0.Theorem 1.
The implicit difference scheme (13) has only one solution, and it is of unconditional stability and convergence for any finite time T>0.Proof.
By the assumptions for the diffusion coefficientD(x), we have pi>0 and pi-qi>0 for i=1,…,M-1. So the coefficient matrix B of (13) is strictly diagonally dominant; then the difference equation (13) has only one solution. Moreover, thanks to the equality ∑k=1nck+dn=1 and with a similar method as used in [22], we get the unconditional stability and convergence of the difference scheme for any given finite time T>0.


## 2.2. The Inverse Problem

For the forward problem (1), (4)-(5), if the fractional orders α and βs (s=1,…,S) in (1) are unknown, we encounter the inverse problem of determining these multiple fractional orders. Suppose that there are some measured points in the space domain; for example, let x0∈Ω be the measured point, and we have the additional information given as (14)ux0,t=θt,0<t≤T,and the inverse problem is to determine α∈(0,1) and βs∈(0,1)(s=1,…,S) using the overposed condition (14) based on (1) and the initial boundary conditions (4) and (5).As stated in Section1, the general uniqueness results have been proved in [39, 42] by Dr. Li and Professor Yamamoto. We give the following lemma for the completeness of the paper.Lemma 2 (see [39]).
Suppose that the fractional ordersα and βs(s=1,…,S), the diffusion coefficient tensor D(x), and the coefficients rs(s=1,…,S) in (1) satisfy the conditions given in Section 1, and the initial function u0(x) is smooth enough, and u0(x)≥0 and u0(x)≢0 for x∈Ω. Then all the fractional orders can be determined uniquely by the additional data {u(x0,t)=θ(t)} for t∈(0,T).The uniqueness result is very important for inverse problems in theory. By Lemma2 we need to utilize the additional data θ(ti) measured at x0 for i=1,2,…; however, in concrete computations, we find that the numerical inversions can also be performed only employing a few of the additional data. So it is still meaningful to study inverse problems from numerics.


## 3. The Inversion Algorithm and Numerical Inversions

In this section, we present numerical inversions for the inverse problem of (1), (4)-(5) with (14). The inversion algorithm we utilize is the homotopy regularization algorithm (see [43, 44], e.g.), which is a combination of the homotopy method with the optimal perturbation algorithm. We give a sketch for the inversion algorithm in the following.


### 3.1. The Homotopy Regularization Algorithm

For the fractional ordersα and βs(s=1,2,…,S) satisfying the order condition (2), let I=(0,1)×(0,1)×⋯×(0,1)⊂RS+1. We denote a vector (15)a=α,β1,β2,…,βS∈I and equip the Euclidean norm (16)a=α2+β12+⋯+βS2.For any givena∈I, denote u(a)(x,t) as the unique solution to the forward problem. Combined with the additional information (14), solving the inverse problem is equivalent to a minimization problem (17)mina∈I⁡uax0,t-θtL20,T2,where u(a)(x0,t)-θ(t)L2(0,T)2=∫0T[u(a)(x0,t)-θ(t)]2dt.Following the homotopy regularization idea, it turns out that the following minimization problem is solved instead of (17): (18)mina∈I⁡1-μuax0,t-θtL20,T2+μa2,where μ∈(0,1) is the homotopy parameter which takes values from near 1 decreasingly approximating to 0. By linearization and numerical differentiation approximations, solving (18) is transformed to solve a normal equation combined with an iteration process via (19)μE+1-μGTGδaj=1-μGTη-ξ,aj+δaj⟶aj+1,j=0,1,…,where δaj is a perturbation vector for any given aj∈I, j denotes the iterative number, and a0 is the initial iteration, and (20)G=gqiQ∗S+1,gqi=uaj+τeix0,tq-uajx0,tqτ;here τ is the numerical differential step and ei(i=1,2,…,S+1) is the unit basis vector in RS+1, and (21)ξ=uajx0,t1,…,uajx0,tQT;η=θt1,…,θtQT.The algorithm can be terminated as long as an optimal perturbation satisfies the conditiona¯j2≤eps; here eps is the given convergent precision. Like that done in [43, 44], we employed the homotopy parameter by (22)μ=μj=11+exp⁡σj-j0,where j is the number of iterations, j0 is a preestimated number, and σ>0 is an adjusted parameter.In the following, we will take two-term case (S=1) and three-term case (S=2) as example to perform the inversion algorithm in 1D case. We set the model parameters l=π, T=1, and D=1 and u0(x)=sin⁡(x) if there is no specification. In addition, we choose j0=5 and σ=0.5 in (22) to determine the homotopy parameter and utilize the grid steps h=l/100 and τ=T/100 in the computation of the forward problem by the difference scheme (13).


### 3.2. Inversion forS=1

In the case ofS=1, let r1=0.5, and the fractional orders we have to determine are α and β1, and a=(α,β1) is the exact solution of the fractional orders. By the exact fractional orders, we compute the forward problem to get the additional data at x0=l/2 with which we reconstruct the fractional orders. It is noticeable that the dimension of the additional data is related to the discretization of the time domain; that is, there are u(l/2,ti) for i=0,1,…,N; here ti is the time at which the measurement is made, and we have a N+1-dimensional vector representing the additional data given as (23)ul2,t0,ul2,t1,…,ul2,tN.However, the inversion could be performed using a few of the additional data, and it is meaningful to investigate the influence of the number (dimension) of the additional data on the inversion algorithm.


#### 3.2.1. Inversion with Accurate Data

Taking two additional datasets at the time oft=0.3 and t=0.4 and choosing the initial iteration as a0=0, the inversion results are listed in Table 1, where ainv is the inversion solution, Err=a-ainv/a denotes the solutions error, and j is the number of iterations.Table 1
The inversion results with different fractional orders forS=1.
a
a
i
n
v
Err
j
(0.9, 0.1)
(0.900000, 0.100000)
1.7583
e
-
9
30
(0.8, 0.2)
(0.800000, 0.200000)
1.2048
e
-
9
30
(0.7, 0.3)
(0.700000, 0.300000)
1.5710
e
-
8
30
(0.6, 0.4)
(0.600000, 0.400000)
2.0909
e
-
8
32
(0.4, 0.3)
(0.400000, 0.300000)
5.9200
e
-
8
34
(0.3, 0.2)
(0.300000, 0.200000)
1.9477
e
-
8
34
(0.2, 0.1)
(0.200000, 0.100000)
1.5804
e
-
7
33Without loss of generality, choosing the exact fractional orders asα=0.8 and β1=0.4, we implement the inversion algorithm using different number of the additional data. The inversion results are listed in Table 2, where L denotes the number (dimension) of the additional data and [t] denotes the measured time to get the additional data.Table 2
The inversion results with number of the additional data forS=1.
L
[
t
]
a
i
n
v
E
r
r
j
2
[
0.3,0.4
]
(0.800000, 0.400000)
2.569
e
-
8
30
3
[
0.3,0.4,0.5
]
(0.800000, 0.400000)
3.288
e
-
9
28
4
[
0.3,0.4,0.5,0.6
]
(0.800000, 0.400000)
1.042
e
-
8
26
5
[
0.3,0.4,0.5,0.6,0.7
]
(0.800000, 0.400000)
3.712
e
-
9
25
6
[
0.3,0.4,0.5,0.6,0.7,0.8
]
(0.800000, 0.400000)
3.448
e
-
9
24
7
[
0.2,0.3,0.4,0.5,0.6,0.7,0.8
]
(0.800000, 0.400000)
6.181
e
-
9
23
8
[
0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9
]
(0.800000, 0.400000)
1.885
e
-
8
22
9
[
0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9
]
(0.800000, 0.400000)
9.076
e
-
9
22From Table1, it can be seen that the inversion can be realized only utilizing two additional datasets, and the inversion solutions are good approximations to the exact solutions. From Table 2, we can see that there are some influences for choosing the number of the additional data on the inversion algorithm, and the inversion results become a little better as the number of the additional data increases and the number of the iterations decreases. Moreover, if using much more additional data in the above, for example, taking all of the measured data given by (23), the inversion solution is still ainv=(0.800000,0.400000), and the solutions error is Err=2.4636e-8, but the number of iterations decreases to j=17.


#### 3.2.2. 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 inversions is often given as(24)θεt=θt+εζ,where ε>0 is the noise level and ζ is a random vector ranging within [-1,1].Also we take the exact fractional orders asα=0.8 and β=0.4; that is, the exact solution of the inverse problem here is given as a=(0.8,0.4). It is noticeable that the inversion with data noises could always fail if still using the above inversion parameters, and we perform the inversion using the completely additional data at ti for i=0,1,…,N and choosing the initial iteration as nonzero vector of a0=(0.1,0.1). The average inversion results with continuous ten-time inversions are listed in Table 3, where a¯inv is the average inversion solution of the ten-time inversions and E¯rr=a¯inv-a2/a2 is the relative average error in the solutions.Table 3
The inversion results using noisy data forS=1.
ε
a
-
i
n
v
E
-rr
j-
5%
(0.789840, 0.422128)
2.7223
e
-
2
14.6
1%
(0.799757, 0.401140)
1.3032
e
-
3
14.3
0.1%
(0.799990, 0.400094)
1.0588
e
-
4
13.9
0.01%
(0.799999, 0.400009)
1.0086
e
-
5
15.0
0.001%
(0.800000, 0.400001)
9.6878
e
-
7
16.0From Table3, we find that the inversion is satisfactory in the case of using random noisy data, and the inversion errors become small when reducing the noise level.


### 3.3. Inversion forS=2

In the case ofS=2, the fractional orders we are to determine are α,β1 and β2, and a=(α,β1,β2) is the exact solution in this case. Let r1=1,r2=0.2, and also choose x0=l/2 as the measured point. We utilize nine additional datasets at the time of ti=i/10 for i=1,2,…,9 and choose the initial iteration as a0=(0.1,0.05,0.01) here. The inversion results with accurate data are listed in Table 4.Table 4
The inversion results with different fractional orders forS=2.
a
a
i
n
v
Err
j
(0.9, 0.7, 0.5)
(0.900000, 0.700000, 0.500000)
1.0216
e
-
9
41
(0.8, 0.6, 0.4)
(0.800000, 0.600000, 0.400000)
4.5193
e
-
9
40
(0.7, 0.5, 0.3)
(0.700000, 0.500000, 0.300000)
4.3115
e
-
8
39
(0.6, 0.5, 0.4)
(0.600000, 0.500000, 0.400000)
2.9329
e
-
9
45
(0.5, 0.4, 0.3)
(0.500000, 0.400000, 0.300000)
2.4011
e
-
9
45
(0.4, 0.3, 0.2)
(0.400000, 0.300000, 0.200000)
8.8878
e
-
9
45
(0.3, 0.2, 0.1)
(0.300000, 0.200000, 0.100000)
9.9418
e
-
8
45As done forS=1, we can also realize the inversion in the case of using noisy data for S=2 although the ill-posedness of the inversion becomes more serious than that of S=1. Taking a=(0.9,0.7,0.5) as an example and using a0=(0.3,0.2,0.1) as the initial iteration, the average results also with continuous ten-time inversions are listed in Table 5, where ε, a¯inv, and E¯rr are all the same as those used in Table 3. It is important to note that we perform the inversion algorithm with convergent precision as eps=1e-2 for ε=1% and ε=0.1% and eps=1e-4 for ε=0.01% and ε=0.001%, respectively.Table 5
The inversion results using noisy data forS=2.
ε
a
-
i
n
v
E
-rr
1%
(0.853161, 0.762050, 0.411797)
9.4438
e
-
2
0.1%
(0.878925, 0.716188, 0.520835)
2.7124
e
-
2
0.01%
(0.900086, 0.699732, 0.500985)
8.2268
e
-
4
0.001%
(0.900007, 0.699989, 0.500042)
3.4895
e
-
5From Tables3 and 5, we can see that the inversion solutions with random noisy data approximate to the exact solutions as the noise level gets smaller. Although the inversion results for S=2 with large noises are not so good as that for S=1, they are still satisfactory and show a numerical stability for the inverse problem.


## 3.1. The Homotopy Regularization Algorithm

For the fractional ordersα and βs(s=1,2,…,S) satisfying the order condition (2), let I=(0,1)×(0,1)×⋯×(0,1)⊂RS+1. We denote a vector (15)a=α,β1,β2,…,βS∈I and equip the Euclidean norm (16)a=α2+β12+⋯+βS2.For any givena∈I, denote u(a)(x,t) as the unique solution to the forward problem. Combined with the additional information (14), solving the inverse problem is equivalent to a minimization problem (17)mina∈I⁡uax0,t-θtL20,T2,where u(a)(x0,t)-θ(t)L2(0,T)2=∫0T[u(a)(x0,t)-θ(t)]2dt.Following the homotopy regularization idea, it turns out that the following minimization problem is solved instead of (17): (18)mina∈I⁡1-μuax0,t-θtL20,T2+μa2,where μ∈(0,1) is the homotopy parameter which takes values from near 1 decreasingly approximating to 0. By linearization and numerical differentiation approximations, solving (18) is transformed to solve a normal equation combined with an iteration process via (19)μE+1-μGTGδaj=1-μGTη-ξ,aj+δaj⟶aj+1,j=0,1,…,where δaj is a perturbation vector for any given aj∈I, j denotes the iterative number, and a0 is the initial iteration, and (20)G=gqiQ∗S+1,gqi=uaj+τeix0,tq-uajx0,tqτ;here τ is the numerical differential step and ei(i=1,2,…,S+1) is the unit basis vector in RS+1, and (21)ξ=uajx0,t1,…,uajx0,tQT;η=θt1,…,θtQT.The algorithm can be terminated as long as an optimal perturbation satisfies the conditiona¯j2≤eps; here eps is the given convergent precision. Like that done in [43, 44], we employed the homotopy parameter by (22)μ=μj=11+exp⁡σj-j0,where j is the number of iterations, j0 is a preestimated number, and σ>0 is an adjusted parameter.In the following, we will take two-term case (S=1) and three-term case (S=2) as example to perform the inversion algorithm in 1D case. We set the model parameters l=π, T=1, and D=1 and u0(x)=sin⁡(x) if there is no specification. In addition, we choose j0=5 and σ=0.5 in (22) to determine the homotopy parameter and utilize the grid steps h=l/100 and τ=T/100 in the computation of the forward problem by the difference scheme (13).


## 3.2. Inversion forS=1

In the case ofS=1, let r1=0.5, and the fractional orders we have to determine are α and β1, and a=(α,β1) is the exact solution of the fractional orders. By the exact fractional orders, we compute the forward problem to get the additional data at x0=l/2 with which we reconstruct the fractional orders. It is noticeable that the dimension of the additional data is related to the discretization of the time domain; that is, there are u(l/2,ti) for i=0,1,…,N; here ti is the time at which the measurement is made, and we have a N+1-dimensional vector representing the additional data given as (23)ul2,t0,ul2,t1,…,ul2,tN.However, the inversion could be performed using a few of the additional data, and it is meaningful to investigate the influence of the number (dimension) of the additional data on the inversion algorithm.


### 3.2.1. Inversion with Accurate Data

Taking two additional datasets at the time oft=0.3 and t=0.4 and choosing the initial iteration as a0=0, the inversion results are listed in Table 1, where ainv is the inversion solution, Err=a-ainv/a denotes the solutions error, and j is the number of iterations.Table 1
The inversion results with different fractional orders forS=1.
a
a
i
n
v
Err
j
(0.9, 0.1)
(0.900000, 0.100000)
1.7583
e
-
9
30
(0.8, 0.2)
(0.800000, 0.200000)
1.2048
e
-
9
30
(0.7, 0.3)
(0.700000, 0.300000)
1.5710
e
-
8
30
(0.6, 0.4)
(0.600000, 0.400000)
2.0909
e
-
8
32
(0.4, 0.3)
(0.400000, 0.300000)
5.9200
e
-
8
34
(0.3, 0.2)
(0.300000, 0.200000)
1.9477
e
-
8
34
(0.2, 0.1)
(0.200000, 0.100000)
1.5804
e
-
7
33Without loss of generality, choosing the exact fractional orders asα=0.8 and β1=0.4, we implement the inversion algorithm using different number of the additional data. The inversion results are listed in Table 2, where L denotes the number (dimension) of the additional data and [t] denotes the measured time to get the additional data.Table 2
The inversion results with number of the additional data forS=1.
L
[
t
]
a
i
n
v
E
r
r
j
2
[
0.3,0.4
]
(0.800000, 0.400000)
2.569
e
-
8
30
3
[
0.3,0.4,0.5
]
(0.800000, 0.400000)
3.288
e
-
9
28
4
[
0.3,0.4,0.5,0.6
]
(0.800000, 0.400000)
1.042
e
-
8
26
5
[
0.3,0.4,0.5,0.6,0.7
]
(0.800000, 0.400000)
3.712
e
-
9
25
6
[
0.3,0.4,0.5,0.6,0.7,0.8
]
(0.800000, 0.400000)
3.448
e
-
9
24
7
[
0.2,0.3,0.4,0.5,0.6,0.7,0.8
]
(0.800000, 0.400000)
6.181
e
-
9
23
8
[
0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9
]
(0.800000, 0.400000)
1.885
e
-
8
22
9
[
0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9
]
(0.800000, 0.400000)
9.076
e
-
9
22From Table1, it can be seen that the inversion can be realized only utilizing two additional datasets, and the inversion solutions are good approximations to the exact solutions. From Table 2, we can see that there are some influences for choosing the number of the additional data on the inversion algorithm, and the inversion results become a little better as the number of the additional data increases and the number of the iterations decreases. Moreover, if using much more additional data in the above, for example, taking all of the measured data given by (23), the inversion solution is still ainv=(0.800000,0.400000), and the solutions error is Err=2.4636e-8, but the number of iterations decreases to j=17.


### 3.2.2. 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 inversions is often given as(24)θεt=θt+εζ,where ε>0 is the noise level and ζ is a random vector ranging within [-1,1].Also we take the exact fractional orders asα=0.8 and β=0.4; that is, the exact solution of the inverse problem here is given as a=(0.8,0.4). It is noticeable that the inversion with data noises could always fail if still using the above inversion parameters, and we perform the inversion using the completely additional data at ti for i=0,1,…,N and choosing the initial iteration as nonzero vector of a0=(0.1,0.1). The average inversion results with continuous ten-time inversions are listed in Table 3, where a¯inv is the average inversion solution of the ten-time inversions and E¯rr=a¯inv-a2/a2 is the relative average error in the solutions.Table 3
The inversion results using noisy data forS=1.
ε
a
-
i
n
v
E
-rr
j-
5%
(0.789840, 0.422128)
2.7223
e
-
2
14.6
1%
(0.799757, 0.401140)
1.3032
e
-
3
14.3
0.1%
(0.799990, 0.400094)
1.0588
e
-
4
13.9
0.01%
(0.799999, 0.400009)
1.0086
e
-
5
15.0
0.001%
(0.800000, 0.400001)
9.6878
e
-
7
16.0From Table3, we find that the inversion is satisfactory in the case of using random noisy data, and the inversion errors become small when reducing the noise level.


## 3.2.1. Inversion with Accurate Data

Taking two additional datasets at the time oft=0.3 and t=0.4 and choosing the initial iteration as a0=0, the inversion results are listed in Table 1, where ainv is the inversion solution, Err=a-ainv/a denotes the solutions error, and j is the number of iterations.Table 1
The inversion results with different fractional orders forS=1.
a
a
i
n
v
Err
j
(0.9, 0.1)
(0.900000, 0.100000)
1.7583
e
-
9
30
(0.8, 0.2)
(0.800000, 0.200000)
1.2048
e
-
9
30
(0.7, 0.3)
(0.700000, 0.300000)
1.5710
e
-
8
30
(0.6, 0.4)
(0.600000, 0.400000)
2.0909
e
-
8
32
(0.4, 0.3)
(0.400000, 0.300000)
5.9200
e
-
8
34
(0.3, 0.2)
(0.300000, 0.200000)
1.9477
e
-
8
34
(0.2, 0.1)
(0.200000, 0.100000)
1.5804
e
-
7
33Without loss of generality, choosing the exact fractional orders asα=0.8 and β1=0.4, we implement the inversion algorithm using different number of the additional data. The inversion results are listed in Table 2, where L denotes the number (dimension) of the additional data and [t] denotes the measured time to get the additional data.Table 2
The inversion results with number of the additional data forS=1.
L
[
t
]
a
i
n
v
E
r
r
j
2
[
0.3,0.4
]
(0.800000, 0.400000)
2.569
e
-
8
30
3
[
0.3,0.4,0.5
]
(0.800000, 0.400000)
3.288
e
-
9
28
4
[
0.3,0.4,0.5,0.6
]
(0.800000, 0.400000)
1.042
e
-
8
26
5
[
0.3,0.4,0.5,0.6,0.7
]
(0.800000, 0.400000)
3.712
e
-
9
25
6
[
0.3,0.4,0.5,0.6,0.7,0.8
]
(0.800000, 0.400000)
3.448
e
-
9
24
7
[
0.2,0.3,0.4,0.5,0.6,0.7,0.8
]
(0.800000, 0.400000)
6.181
e
-
9
23
8
[
0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9
]
(0.800000, 0.400000)
1.885
e
-
8
22
9
[
0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9
]
(0.800000, 0.400000)
9.076
e
-
9
22From Table1, it can be seen that the inversion can be realized only utilizing two additional datasets, and the inversion solutions are good approximations to the exact solutions. From Table 2, we can see that there are some influences for choosing the number of the additional data on the inversion algorithm, and the inversion results become a little better as the number of the additional data increases and the number of the iterations decreases. Moreover, if using much more additional data in the above, for example, taking all of the measured data given by (23), the inversion solution is still ainv=(0.800000,0.400000), and the solutions error is Err=2.4636e-8, but the number of iterations decreases to j=17.


## 3.2.2. 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 inversions is often given as(24)θεt=θt+εζ,where ε>0 is the noise level and ζ is a random vector ranging within [-1,1].Also we take the exact fractional orders asα=0.8 and β=0.4; that is, the exact solution of the inverse problem here is given as a=(0.8,0.4). It is noticeable that the inversion with data noises could always fail if still using the above inversion parameters, and we perform the inversion using the completely additional data at ti for i=0,1,…,N and choosing the initial iteration as nonzero vector of a0=(0.1,0.1). The average inversion results with continuous ten-time inversions are listed in Table 3, where a¯inv is the average inversion solution of the ten-time inversions and E¯rr=a¯inv-a2/a2 is the relative average error in the solutions.Table 3
The inversion results using noisy data forS=1.
ε
a
-
i
n
v
E
-rr
j-
5%
(0.789840, 0.422128)
2.7223
e
-
2
14.6
1%
(0.799757, 0.401140)
1.3032
e
-
3
14.3
0.1%
(0.799990, 0.400094)
1.0588
e
-
4
13.9
0.01%
(0.799999, 0.400009)
1.0086
e
-
5
15.0
0.001%
(0.800000, 0.400001)
9.6878
e
-
7
16.0From Table3, we find that the inversion is satisfactory in the case of using random noisy data, and the inversion errors become small when reducing the noise level.


## 3.3. Inversion forS=2

In the case ofS=2, the fractional orders we are to determine are α,β1 and β2, and a=(α,β1,β2) is the exact solution in this case. Let r1=1,r2=0.2, and also choose x0=l/2 as the measured point. We utilize nine additional datasets at the time of ti=i/10 for i=1,2,…,9 and choose the initial iteration as a0=(0.1,0.05,0.01) here. The inversion results with accurate data are listed in Table 4.Table 4
The inversion results with different fractional orders forS=2.
a
a
i
n
v
Err
j
(0.9, 0.7, 0.5)
(0.900000, 0.700000, 0.500000)
1.0216
e
-
9
41
(0.8, 0.6, 0.4)
(0.800000, 0.600000, 0.400000)
4.5193
e
-
9
40
(0.7, 0.5, 0.3)
(0.700000, 0.500000, 0.300000)
4.3115
e
-
8
39
(0.6, 0.5, 0.4)
(0.600000, 0.500000, 0.400000)
2.9329
e
-
9
45
(0.5, 0.4, 0.3)
(0.500000, 0.400000, 0.300000)
2.4011
e
-
9
45
(0.4, 0.3, 0.2)
(0.400000, 0.300000, 0.200000)
8.8878
e
-
9
45
(0.3, 0.2, 0.1)
(0.300000, 0.200000, 0.100000)
9.9418
e
-
8
45As done forS=1, we can also realize the inversion in the case of using noisy data for S=2 although the ill-posedness of the inversion becomes more serious than that of S=1. Taking a=(0.9,0.7,0.5) as an example and using a0=(0.3,0.2,0.1) as the initial iteration, the average results also with continuous ten-time inversions are listed in Table 5, where ε, a¯inv, and E¯rr are all the same as those used in Table 3. It is important to note that we perform the inversion algorithm with convergent precision as eps=1e-2 for ε=1% and ε=0.1% and eps=1e-4 for ε=0.01% and ε=0.001%, respectively.Table 5
The inversion results using noisy data forS=2.
ε
a
-
i
n
v
E
-rr
1%
(0.853161, 0.762050, 0.411797)
9.4438
e
-
2
0.1%
(0.878925, 0.716188, 0.520835)
2.7124
e
-
2
0.01%
(0.900086, 0.699732, 0.500985)
8.2268
e
-
4
0.001%
(0.900007, 0.699989, 0.500042)
3.4895
e
-
5From Tables3 and 5, we can see that the inversion solutions with random noisy data approximate to the exact solutions as the noise level gets smaller. Although the inversion results for S=2 with large noises are not so good as that for S=1, they are still satisfactory and show a numerical stability for the inverse problem.


## 4. Conclusions

The numerical determination problem for the fractional orders in the multiterm TFDE is investigated using some measurements at the interior point of the space domain. The inversion problem is unique, and numerical inversions with random noisy data are performed successfully by using the homotopy regularization algorithm. The homotopy regularization algorithm can also be utilized to determine the multiple fractional orders in the multidimensional case as long as a solution to the forward problem can be worked out.It is noted that numerical inversions can be implemented smoothly in the case that the fractional orders satisfy the order condition (2) given in Section 1, which is just a necessary condition for the uniqueness of the inverse problem. The inversion results are very satisfactory if coping with accurate data; however, they become a little bad in the case of S=2 when using noisy data with noises greater than 1%. Therefore, we can say that the stability of the inverse problem here with more fractional derivatives, that is, S≥2, could be severely ill-posed in spite of the fact that uniqueness is valid. We will deal with the inverse problem in the multidimensional case and we have to seek more effective inversion algorithms for the multiterm TFDE with S>2 in the future work.



---

*Source: 3204959-2017-08-17.xml*