# Recovery of the time-dependent source term in the stochastic fractional diffusion equation with heterogeneous medium

In this work, an inverse problem in the fractional diffusion equation with random source is considered. The measurements used are the statistical moments of the realizations of single point data u(x_0,t,ω). We build the representation of the solution u in integral sense, then prove that the unknowns can be bounded by the moments theoretically. For the numerical reconstruction, we establish an iterative algorithm with regularized Levenberg-Marquardt type and some numerical results generated from this algorithm are displayed. For the case of highly heterogeneous media, the Generalized Multiscale finite element method (GMsFEM) will be employed.

## Authors

• 13 publications
• 3 publications
01/06/2021

### Numerical analysis for stochastic time-space fractional diffusion equation driven by fractional Gaussion noise

In this paper, we consider the strong convergence of the time-space frac...
12/13/2019

### Inverse source in two-parameter anomalous diffusion, numerical algorithms and simulations over graded time-meshes

We consider an inverse source two-parameter sub-diffusion model subject ...
12/20/2020

### Numerical solution of an inverse random source problem for the time fractional diffusion equation via PhaseLift

This paper is concerned with the inverse random source problem for a sto...
09/07/2021

### Numerical approximations for the fractional Fokker-Planck equation with two-scale diffusion

Fractional Fokker-Planck equation plays an important role in describing ...
06/03/2019

### A new nonlocal forward model for diffuse optical tomography

The forward model in diffuse optical tomography (DOT) describes how ligh...
09/07/2021

### L1 scheme on graded mesh for subdiffusion equation with nonlocal diffusion term

The solution of time fractional partial differential equations in genera...
04/07/2021

### Upscaling errors in Heterogeneous Multiscale Methods for the Landau-Lifshitz equation

In this paper, we consider several possible ways to set up Heterogeneous...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

### 1.1 Mathematical statement

The mathematical model in this work is stated as follows:

 ⎧⎪⎨⎪⎩∂αtu+Au=f(x)[g1(t)+g2(t)˙W(t)]=:f(x)g(t,ω),(x,t)∈D×(0,T],u(x,t)=0,(x,t)∈∂D×(0,T],u(x,0)=0,x∈D. (1)

The domain has sufficiently smooth boundary, and with denotes the Djrbashyan-Caputo fractional derivative, defined as

 ∂αtψ(t)=1Γ(1−α)∫t0(t−τ)−αψ′(τ) dτ,

where is the gamma function. The lower bound is set to ensure the well definedness of the Ito integral and this can be seen in the next section. The operator is an elliptic operator defined as , and may be highly heterogeneous. The source term contains the targeted unknowns , while the spatial component is given.

is the white noise derived from the Brownian motion and then

constitutes an Ito process, see [41] for details.

Our data is the moments of the realizations of on a single point with the restriction , which is different from [31]. This condition will make the inverse problem more challenging in mathematics, but is meaningful in practical application. For instance, regarding equation (1) as the contaminant diffusion system, solution will be the concentration of pollutant, and is the location of pollution source, in which it is severely polluted. If the pollutant is harmful for human body, it is not allowed to observe inside the support of considering staff’s health. Reflected on mathematics, we should set the restriction , even though it will increase the difficulty. Actually, given , this inverse problem will be reduced to Volterra integral equations, and the stability result and iterative algorithm follow straightforwardly. See Remark (1) for details.

The precise mathematical description of this inverse problem is given as follows: observe , then use the statistical moments of the measurements to reconstruct simultaneously.

### 1.2 Physical background and literature

In microscopic level, the random motion of a single particle can be viewed as a diffusion process. The classical diffusion equation can be deduced to describe the motion of particles, if we assume the key condition, the mean squared displacement of jumps after a long time is proportional to time, i.e. . However, recently, people found some anomalous diffusion phenomena [5, 8, 17, 26], in which the assumption is violated. Sometimes it may possess the asymptotic behavior of , i.e. The different rate will lead to a reformulation to the diffusion equation, introducing the time fractional derivative in it, and the corresponding equations are called fractional differential equations (FDEs). We list some applications of FDEs, to name a few, the thermal diffusion in media with fractal geometry [39], ion transport in column experiments [19], dispersion in a heterogeneous aquifer [1], non-Fickian diffusion in geological formations [6], the analysis on viscoelasticity in material science [35, 49]. [37] provides an extensive list.

If uncertainty is added in the source term, the FDE system will become more complicated and meaningful. Since it is common to meet a diffusion source, which is defined as a stochastic process to describe the uncertain character imposed by nature. As a consequence, it is worth to investigate the diffusion system with a random source. In such situation the solution will be written as a stochastic process, which makes the analysis more challenging.

In addition, to deal with the case of highly heterogeneous medium , the Generalized Multiscale Finite Element Method (GMsFEM [13]) will be used to simulate the forward problem of equation (1). The introduction of GMsFEM will be given in section (4.3.1).

For a comprehensive understanding of fractional calculus and FDEs, see [25, 45, 7] and the references therein. For inverse problems in FDEs, [24] is an extensive review. See [14, 50, 40] for inverse source and coefficient problems; see [22] for unique continuation principle; see [21]

for Carleman estimate in FDEs; see

[18, 27, 42] for fractional Calderon problem. Furthermore, if we extend the assumption to a more general case , the multi-term fractional diffusion equations and even the distributed-order differential equations will be generated, [43, 47, 30, 29]. For numerical methods for inverse problems, see [4, 3] and the references therein. Literature about the GMsFEM and its applications can be found in [13, 12, 9, 11, 10].

### 1.3 Main result and outline

Throughout this paper, the following restrictions on spatial component , observation point , and the unknowns are supposed to be valid.

###### Assumption 1.
• , and set such that ;

• changes its sign times on and ;

• and for ;

• , namely, .

Now we can state the main result, which says the unknowns can be limited by some statistical moments of observations .

###### Theorem 1.

Under Assumption (1), let satisfy equation (6) and define

 Cα=1/Γ(2−α), Bη=∥v(x0,⋅)∥−1L1(0,η), η>0 be small.

Then the following estimates for and are valid.

• If ,

• If ,

 ∥g1∥L1(0,T−η)≤(BηCfT+1)N+1−1BηCfT(CαBηT1−α E[∥u(x0,⋅,ω)∥L1(0,T)]+2Mη).
•  ∥g2∥L2(0,T−η)≤η1/2Bη ∥∥V[I1−αtu(x0,t,ω)]∥∥1/2L1(0,T).

In this theorem, means the fractional integral operator, and

are the notations for expectation and variance, respectively. These knowledge can be seen in section (

2). Noting that the stochastic process is independent of the sign of by the properties of in section (2), sequentially we consider instead of . That’s why the norm of is estimated.

The remaining part of this manuscript is structured as follows. Section (2

) includes the preliminaries, such as the probability space

and the stochastic solution for equation (1). Also, some auxiliary results like the reverse convolution inequality and the maximum principles in fractional diffusion equations are collected. In section (3), we prove Theorem (1). After that the numerical reconstruction for the unknowns is investigated in section (4). We construct the regularized Levenberg-Marquardt iteration (15), and prove its convergence–Proposition (2). Some numerical results generated by iteration (15) are also displayed. Furthermore, some brief knowledge of GMsFEM is provided in this section.

## 2 Preliminaries

### 2.1 Brownian motion and Ito isometry formula

To state the Ito formula, firstly we need to give the setting of probability space.

###### Definition 1.

We call a probability space if denotes the nonempty sample space, is the algebra of and is the probability measure.

With the above definition, the expectation and variance

of a random variable

can be given as

 E[X]=∫ΩX(ω) dP(ω), V[X]=E[(X−E[X])2].

The Brownian motion , which is also called Wiener process in mathematics, has the following properties,

• has continuous paths;

• has independent increments and satisfies

 W(t)−W(s)∼N(0,t−s), 0≤s≤t,

where

is the normal distribution.

Now the essential tool, Ito isometry formula can be stated.

###### Lemma 1.

([41]). Let be a probability space and satisfy the following properties.

• is -measurable, where denotes the Borel -algebra on

• for some .

Then the Ito integral , where denotes the random measure derived from , is well defined, and it follows that

 E[(∫S0ψ(t,ω) dW(t))2]=E[∫S0ψ2(t,ω) dt].

### 2.2 Stochastic weak solution

The randomness from means that we can not differentiate in for each . As a consequence, we will define the weak solution of equation (1) in the integral sense.

Firstly, the fractional integral operator and the corresponding Ito integral are given.

###### Definition 2.

The fractional integral operator is defined as

 Iαtψ(t)=Γ(α)−1∫t0(t−τ)α−1ψ(τ) dτ,t>0.

Then we define as

 Iαtg(t,ω)=Iαtg1(t)+Γ(α)−1∫t0(t−τ)α−1g2(τ) dW(τ).

Now we explain the necessity of the restriction . For , from the conditions and , we have is square-integrable on Then Lemma (1) yields that the Ito integral is well defined.

In addition, the direct calculation gives that

 Iαt∂αtψ(t)=ψ(t)−ψ(0),

which implies the next definition of the weak solution for equation (1).

###### Definition 3 (Stochastic weak solution).

The stochastic process is called as a stochastic weak solution of equation (1) if for each and , it holds that

 ⟨u(⋅,t,ω),ψ(⋅)⟩L2(D)+⟨IαtAu(⋅,t,ω),ψ(⋅)⟩L2(D)=Iαtg(t,ω) ⟨f(⋅),ψ(⋅)⟩L2(D), t∈(0,T].

### 2.3 Auxiliary lemmas

Here we list some auxiliary lemmas which will be used later. First the reverse convolution inequality is given.

###### Lemma 2.

([33, Lemma 3.1]). Let and be arbitrarily given. Suppose that , and on .

• If keeps its sign on , then

 (2)
• If only keeps its sign on , then

 ∥φ1∥L1(T1,T2)∥φ2∥L1(0,η)≤ (3) +2∥φ1∥L1(T2,T2+η)∥φ2∥L1(0,η).

The next lemmas are the maximum principles in FDEs.

###### Lemma 3.

(Maximum principle, [34, Theorem 2]). Fix let satisfy the following fractional diffusion equation

 ∂αtψ+Aψ=F(x,t), (x,t)∈D×(0,T), (4)

and define . If , then

 ψ(x,t)≤max{0,max{ψ(x,t):(x,t)∈λT}}, (x,t)∈¯¯¯¯¯D×[0,T].
###### Lemma 4.

(Strong maximum principle, [32, Theorem 1.1]). We set as the solution of equation (4) and let on . Then for any , the set is at most a finite set.

In addition, we gives a representation lemma for the weak solution .

###### Lemma 5.

([31, Lemma 5]) The weak solution of equation (1) can be written as

 u(x,t,ω)=Iαtg(t,ω)f(x)+∫t0Iαtg(τ,ω)vt(x,t−τ) dτ,t∈(0,T], (5)

where is the solution of the following deterministic fractional diffusion equation

 ⎧⎪⎨⎪⎩∂αtv+Av=0,(x,t)∈D×(0,T],v(x,t)=0,(x,t)∈∂D×(0,T],v(x,0)=f(x),x∈D. (6)

## 3 Main result

### 3.1 Moments representation

With Lemmas (1) and (5), the moments we used can be represented in terms of the unknowns . See the lemma below.

###### Lemma 6.
 E[I1−αtu(x0,t,ω)] =∫t0g1(τ)v(x0,t−τ) dτ, (7) V[I1−αtu(x0,t,ω)] =∫t0g22(τ)[v(x0,t−τ)]2 dτ,t∈(0,T].
###### Proof.

With Lemma (5) and the fact that

 ∫ts(t−τ)−α(τ−s)α−1 dτ=B(1−α,α)=Γ(1−α)Γ(α)/Γ(1),

here is the Beta function, the next result can be deduced,

 I1−αtu(x,t,ω)= f(x)Γ(α)Γ(1−α)∫t0(t−τ)−α∫τ0(τ−s)α−1g(s,ω) ds dτ +1Γ(1−α)∫t0(t−τ)−α∫τ0Iαtg(τ−s,ω)vt(x,s) ds dτ = 1Γ(α)Γ(1−α)[f(x)∫t0g(s,ω)∫ts(t−τ)−α(τ−s)α−1 dτ ds +∫t0vt(x,s)∫t−s0g(r,ω)∫tr+s(t−τ)−α(τ−s−r)1−α dτ dr ds] = f(x)∫t0g(s,ω) ds+∫t0vt(x,s)∫t−s0g(r,ω) dr ds = f(x)∫t0g(s,ω) ds+∫t0g(r,ω)[v(x,t−r)−v(x,0)] dr = ∫t0g(τ,ω)v(x,t−τ) dτ.

Then we have

 I1−αtu(x0,t,ω)=∫t0g1(τ)v(x0,t−τ) dτ+∫t0g2(τ)v(x0,t−τ) dW(τ).

Applying Ito formula in Lemma (1) to the above equality leads to (7). ∎

###### Remark 1.

In [31], the authors use integration by parts on the right side of (7) to deduce the following second kind Volterra equations,

 G1(t) =f−1(x0)E[I1−αtu(x0,t,ω)]−f−1(x0)∫t0G1(τ)vt(x0,t−τ) dτ, G2(t) =f−2(x0)V[I1−αtu(x0,t,ω)]−2f−2(x0)∫t0G2(τ)v(x0,t−τ)vt(x0,t−τ) dτ,

where

 G1(t)=∫t0g1(τ) dτ,G2(t)=∫t0g22(τ) dτ.

However, since in this work, we can only start the analysis from (7). Due to the convolution structure, the estimates of the unknowns on the partial interval are attained. See the next subsection for details.

### 3.2 Proof of Theorem (1)

From (7), we build the proof of Theorem (1).

###### Proof of Theorem (1) (a).

Let , , then inserting (2) to (7) straightforwardly yields that

 ∥∥E[I1−αtu(x0,t,ω)]∥∥L1(0,T)≥∥g1∥L1(0,T−η) ∥v(x0,⋅)∥L1(0,η).

For the left side, we have

 ∥∥E[I1−αtu(x0,t,ω)]∥∥L1(0,T) ≤1Γ(1−α)E[∫T0∫t0(t−τ)−α|u(x0,τ,ω)| dτ dt] =1Γ(2−α)E[∫T0(T−τ)1−α|u(x0,τ,ω)| dτ] ≤CαT1−αE[∥u(x0,⋅,ω)∥L1(0,T)],

then

###### Proof of Theorem (1) (b).

Let be small and assume that changes sign on , for convenience, we set and . By (7), for , we can write

 ∫ttkg1(τ)v(x0,t−τ) dτ=E[I1−αtu(x0,t,ω)]−k∑j=1Sj,

where

 Sj=∫tjtj−1g1(τ)v(x0,t−τ) dτ.

Using (3) to the above equality with , we can obtain that for ,

 ∥g1∥L1(tk,tk+1)≤ Bη∥∥∫ttkg1(τ)v(x0,t−τ) dτ∥∥L1(tk,tk+1+η)+2∥g1∥L1(tk+1,tk+1+η) (8) ≤ Bη(∥∥E[I1−αtu(x0,t,ω)]∥∥L1(tk,tk+1+η)+k∑j=1∥Sj∥L1(tk,tk+1+η)) +2∥g1∥L1(tk+1,tk+1+η).

For , from the condition that we have

 ∥g1∥L1(tk+1,tk+1+η)= ∫tk+1+ηtk+1|g1(τ)| dτ≤Mη. (9)

For , it holds that

 ∥Sj∥L1(tk,tk+1+η)≤ ∫tk+1+ηtk∫tjtj−1|g1(τ)|v(x0,t−τ) dτ dt = ∫tjtj−1|g1(τ)|∫tk+1+ηtkv(x0,t−τ) dt dτ = ∫tjtj−1|g1(τ)| ∥v(x0,⋅)∥L1(tk−τ,tk+1+η−τ) dτ.

Assumption (1) and Lemma (3) give that . Consequently,

 ∥Sj∥L1(tk,tk+1+η)≤ CfT∥g1∥L1(tj−1,tj),j=1,⋯,k. (10)

Inserting (9) and (10) into (8) yields that

 ∥g1∥L1(tk,tk+1)≤ Bη∥∥E[I1−αtu(x0,t,ω)]∥∥L1(tk,tk+1+η)+BηCfT∥g1∥L1(0,tk)+2Mη. (11)

Fix , we have

 ∥g1∥L1(0,t1)≤Bη∥∥E[I1−αtu(x0,t,ω)]∥∥L1(0,t1+η)+2Mη. (12)

Now we claim that for ,

 ∥g1∥L1(0,tk)≤(BηCfT+1)k−1BηCfT(Bη∥∥E[I1−αtu(x0,t,ω)]∥∥L1(0,tk+η)+2Mη),

and prove it by induction. The case of is valid by (12). Now assume that the claim holds for , then for , the estimate (11) gives that

 ∥g1∥L1(0,tl+1)≤ ∥g1∥L1(0,tl)+Bη∥∥E[I1−αtu(x0,t,ω)]∥∥L1(tl,tl+1+η)+BηCfT∥g1∥L1(0,tl)+2Mη ≤ (BηCfT+1)(BηCfT+1)l−1BηCfT(Bη∥∥E[I1−αtu(x0,t,ω)]∥∥L1(0,tl+η)+2Mη) +Bη∥∥E[I1−αtu(x0,t,ω)]∥∥L1(tl,tl+1+η)+2Mη ≤ (BηCfT+1)l+1−1BηCfT(Bη∥∥E[I1−αtu(x0,t,ω)]∥∥L1(0,tl+1+η)+2Mη).

So the claim is valid, and recalling that , we have

 ∥g1∥L1(0,T−η)≤ (BηCfT+1)N+1−1BηCfT(Bη∥∥E[I1−αtu(x0,t,ω)]∥∥L1(0,T)+2Mη) ≤ (BηCfT+1)N+1−1BηCfT(CαBηT1−α E[∥u(x0,⋅,ω)∥L1(0,T)]+2Mη).

The proof is complete. ∎

###### Proof of Theorem (1) (c).

Note that keeps its sign on . Analogous to the proof of Theorem (1) , setting , in (2), then (7) gives that

 ∥∥V[I1−αtu(x0,t,ω)]∥∥L1(0,T) ≥∥g22∥L1(0,T−η) ∥[v(x0,⋅)]2∥L1(0,η) =∥g2∥2L2(0,T−η) ∥v(x0,⋅)∥2L2(0,η).

Holder inequality yields that

 ∥v(x0,⋅)∥L2(0,η)≥∥1∥−1L2(0,η)∥v(x0,⋅)∥L1(0,η)=η−1/2B−1η.

Consequently,

 ∥g2∥L2(0,T−η)≤η1/2Bη ∥∥V[I1−αtu(x0,t,ω)]∥∥1/2L1(0,T).

The proof of Theorem (1) is complete. ∎

## 4 Numerical reconstruction

### 4.1 Regularized Levenberg-Marquardt iteration

The discretized formulation of integral equation (7) is derived as follows. Denote the uniform mesh on the interval as and set . From Lemmas (3) and (4), we have is at most a finite set. Thus we can set

 v(x0,t1)>0, (13)

if the mesh size is chosen appropriately.

Define

 E(tn)=E[I1−αtu(x0,tn,ω)],V(tn)=V[I1−αtu(x0,tn,ω)].

Then from (7) we have

 E(tn) =∫tn0g1(τ)v(x0,tn−τ) dτ=n∑k=1∫tktk−1g1(τ)v(x0,tn−τ) dτ ≈Δtn∑k=1[g1(tk−1)v(x0,tn−tk−1)+g1(tk)v(x0,tn−tk)]/2 =Δt[g1(0)v(x0,tn)/2+g1(tn)v(x0,0)/2+n−1∑k=1g