Identification of the degradation coefficient for an anomalous diffusion process in hydrology

TL;DR

This paper develops a hybrid variational regularization and Laplace approximation method to accurately identify the degradation coefficient in anomalous diffusion processes in hydrology, addressing ill-posed inverse problems.

Contribution

It introduces a novel combined deterministic-stochastic approach for inverse degradation coefficient problems, ensuring uniqueness, stability, and quantifying uncertainty.

Findings

Proved unique determination of the degradation coefficient from average flux data.

Established existence, stability, and convergence of the solution.

Numerical examples demonstrate the method's efficiency and robustness.

Abstract

In hydrology, the degradation coefficient is one of the key parameters to describe the water quality change and to determine the water carrying capacity. This paper is devoted to identify the degradation coefficient in an anomalous diffusion process by using the average flux data at the accessible part of boundary. The main challenges in inverse degradation coefficient problems (IDCP) is the average flux measurement data only provide very limited information and cause the severe ill-posedness of IDCP. Firstly, we prove the average flux measurement data can uniquely determine the degradation coefficient. The existence and uniqueness of weak solution for the direct problem are established, and the Lipschitz continuity of the corresponding forward operator is also obtained. Secondly, to overcome the ill-posedness, we combine the variational regularization method with Laplace approximations (LA) to solve the IDCP. This hybrid method is essentially the combination of deterministic regularization method and stochastic method. Thus, it is able to calculate the minimizer (MAP point) more rapidly and accurately, but also enables captures the statistics information and quantifying the uncertainty of the solution. Furthermore, the existence, stability and convergence of the minimizer of the variational problem are proved. The convergence rate estimate between the LA posterior distribution and the actual posterior distribution in the sense of Hellinger distance is given, and the skewness are introduced for characterizing the symmetry or slope of LA solution, especially the relationship with the symmetry of the measurement data. Finally, the one-dimensional and two-dimensional numerical examples are presented to confirm the efficiency and robustness of the proposed method.

Tables

Table 1. Table 1: Numerical results for Example 1 with various N and ε 𝜀 \varepsilon ( α = 0.3 𝛼 0.3 \alpha=0.3 ).

$N\setminus \epsilon$	$0.0001$	$0.0005$	$0.001$
$1$	$0.0525$	$0.0528$	$0.0590$
$2$	$0.0444$	$0.0477$	$0.0506$
$3$	$0.0111$	$0.0287$	$0.0409$
$4$	$0.0095$	$0.0257$	$0.0344$
$5$	$0.0090$	$0.0197$	$0.0295$

Table 2. Table 2: Numerical results for Example 1 with various μ 𝜇 \mu and ε 𝜀 \varepsilon ( α = 0.3 𝛼 0.3 \alpha=0.3 ).

$\mu \setminus \epsilon$	$0.0001$	$0.0005$	$0.001$
$\delta$	$0.1086$	$0.2090$	$0.2970$
${\delta }^{\frac{3}{2}}$	$0.0090$	$0.0197$	$0.0295$
${\delta }^2$	$0.0386$	$0.0608$	$0.1437$

Table 3. Table 3: Numerical results for Example 1 with various α 𝛼 \alpha and ε 𝜀 \varepsilon .

$\alpha \setminus \epsilon$	$0.0001$	$0.0005$	$0.001$
$0.1$	$0.0086$	$0.0117$	$0.0252$
$0.3$	$0.0090$	$0.0197$	$0.0296$
$0.5$	$0.0079$	$0.0204$	$0.0398$
$0.7$	$0.0074$	$0.0351$	$0.0517$
$0.9$	$0.0071$	$0.0447$	$0.0642$

Table 4. Table 4: Numerical results for Example 1 for different basis functions and ε 𝜀 \varepsilon with α = 0.3 𝛼 0.3 \alpha=0.3 .

Type of basis function$\setminus \epsilon$	$0.0001$	$0.0005$	$0.001$
polynomial basis	$0.0129$	$0.0279$	$0.0436$
Trigonometric basis	$0.0090$	$0.0197$	$0.0295$

Table 5. Table 5: Numerical results for Example 1 for various type of data and ε 𝜀 \varepsilon with α = 0.3 𝛼 0.3 \alpha=0.3 .

Type of data$\setminus \epsilon$	$0.0001$	$0.0005$	$0.001$	$0.005$	$0.01$	$0.05$
Average flux data	$0.0090$	$0.0197$	$0.0295$	$0.0662$	$0.1371$	$0.2307$
Direct flux data	$0.0116$	$0.0125$	$0.0203$	$0.0445$	$0.0497$	$0.0539$

Table 6. Table 6: Numerical results for Example 3 for the skewness of exact solution is 0

Measurement data	skewness
Data on ${\Lambda }_0$	$0.3285$
Data on ${\Lambda }_1$	$-0.2088$
Data on ${\Lambda }_0\cup {\Lambda }_1$	$-0.0033$

Table 7. Table 7: Numerical results for Example 7 for the skewness of exact solution is 0

Measurement data	skewness
Data on ${\Lambda }_1\cup {\Lambda }_2$	$(-0.2905,0.2551)$
Data on ${\Lambda }_1\cup {\Lambda }_2\cup {\Lambda }_4$	$(-0.0170,0.1308)$
Data on $\partial \Omega$	$(-0.0150,-0.0150)$

Table 8. Table 8: Results for Example 7 for the 95% confidence interval at different ( x , y ) ∈ Ω 𝑥 𝑦 Ω (x,y)\in\Omega

y$\setminus x$	$0.2$	$0.5$	$0.8$
0.2	$[3.1721\times {10}^{-2},3.1814\times {10}^{-2}]$	$[4.5196\times {10}^{-2},4.5230\times {10}^{-2}]$	$[3.1716\times {10}^{-2},3.1810\times {10}^{-2}]$
0.5	$[4.4759\times {10}^{-2},4.4853\times {10}^{-2}]$	$[6.1178\times {10}^{-2},6.1279\times {10}^{-2}]$	$[4.4749\times {10}^{-2},4.4851\times {10}^{-2}]$
0.8	$[3.1713\times {10}^{-2},3.1812\times {10}^{-2}]$	$[4.5174\times {10}^{-2},4.5268\times {10}^{-2}]$	$[3.1711\times {10}^{-2},3.1801\times {10}^{-2}]$

Equations

$$q=-\kappa (x,u)\nabla u,$$

$$\frac{\partial u}{\partial t}+\nabla \cdot q=f(x,t,u),$$

$$\frac{\partial u}{\partial t}-\nabla \cdot (\kappa (x,u)\nabla u)=f(x,t,u).$$

$${}_0{D}_t^{\alpha }u+\nabla \cdot q=f(x,t,u),$$

$$\begin{array}{r}{}_0{D}_t^{\alpha }u=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t(t-s{)}^{-\alpha }\frac{\partial u(x,s)}{\partial s}ds.\end{array}$$

$${}_0{D}_t^{\alpha }u-\nabla \cdot (\kappa (x,u)\nabla u)=f(x,t,u).$$

$$\{\begin{array}{l}\begin{array}{rl}{}_0{D}_t^{\alpha }u-\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial u}{\partial x_l})+q(x)u& =\varphi (x)v(t),\text{in}D=\Omega \times (0,T),\\ u& =0,\text{on}\partial \Omega \times (0,T),\\ u& =0,\text{in}\Omega \times \{0\},\end{array}\end{array}$$

$$\begin{array}{r}\sum _{k,l=1}^d{a}_{kl}(x){\xi }_k{\xi }_l\ge {\lambda }_0\sum _{k=1}^d{\xi }_k^2,x\in \overline{\Omega },\xi \in R^d,{\lambda }_0>0.\end{array}$$

$$\begin{array}{r}{\int }_{\Lambda \times (0,T)}\frac{\partial u_j^i(x,t;q)}{\partial {\nu }_A}h(s,t)dsdt\end{array}$$

$$\begin{array}{r}\frac{\partial u_j^i}{\partial {\nu }_A}=\sum _{k,l=1}^d{a}_{kl}(x)\frac{\partial u_j^i}{\partial x_k}{\nu }_l(x),\end{array}$$

$$\begin{array}{rl}{\int }_{\Lambda \times (0,T)}\frac{\partial u_j^1(x,t;p)}{\partial {\nu }_A}h(x,t)dxdt& ={\int }_{\Lambda \times (0,T)}\frac{\partial u_j^1(x,t;q)}{\partial {\nu }_A}h(x,t)dxdt;\\ {\int }_{\Lambda \times (0,T)}\frac{\partial u_j^2(x,t;p)}{\partial {\nu }_A}h(x,t)dxdt& ={\int }_{\Lambda \times (0,T)}\frac{\partial u_j^2(x,t;q)}{\partial {\nu }_A}h(x,t)dxdt,\end{array}$$

$$\{\begin{array}{l}\begin{array}{rl}{}_t{D}_T^{\alpha }w-\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial w}{\partial x_l})+q(x)w& =0,\text{in}D=\Omega \times (0,T),\\ w& =h_0(x,t),\text{on}\partial \Omega \times (0,T),\\ w& =0,\text{in}\Omega \times \{T\}.\end{array}\end{array}$$

$$\{\begin{array}{l}\begin{array}{rl}{}_0{D}_t^{\alpha }\tilde{w}-\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial \tilde{w}}{\partial x_l})+q(x)\tilde{w}& =0,\text{in}D,\\ \tilde{w}& =h_0(x,T-t),\text{on}\partial \Omega \times (0,T),\\ \tilde{w}& =0,\text{in}\Omega \times \{0\}.\end{array}\end{array}$$

$$\begin{array}{rl}& {\int }_0^T{\int }_{\Omega }{\varphi }_j(x)v(t)w(x,t;q)dxdt\\ =& {\int }_0^T{\int }_{\Omega }{[}_0{D}_t^{\alpha }{u}_j^1(x,t;q)-\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial u_j^1(x,t;q)}{\partial x_l})+q(x)u_j^1(x,t;q)]w(x,t;q)dxdt\\ =& {\int }_0^T{\int }_{\Omega }({}_t{D}_T^{\alpha }w)u_j^{1}dxdt-{\int }_0^T{\int }_{\Omega }\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial u_j^1}{\partial x_l})wdxdt+{\int }_0^T{\int }_{\Omega }q(x)u_j^{1}wdxdt\\ =& {\int }_0^T{\int }_{\Omega }({}_t{D}_T^{\alpha }w)u_j^{1}dxdt-{\int }_0^T{\int }_{\Omega }\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial w}{\partial x_l})u_j^{1}dxdt+{\int }_0^T{\int }_{\Omega }q(x)u_j^{1}wdxdt\\ & +{\int }_0^T{\int }_{\partial \Omega }(\sum _{k,l=1}^d{a}_{kl}(x)\frac{\partial w}{\partial x_l}{\nu }^k)u_j^{1}dxdt-{\int }_0^T{\int }_{\partial \Omega }(\sum _{k,l=1}^d{a}_{kl}(x)\frac{\partial u_j^1}{\partial x_l}{\nu }^k)wdxdt\\ =& {\int }_0^T{\int }_{\Omega }{[}_t{D}_T^{\alpha }w-\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial w}{\partial x_l})+q(x)w]u_j^{1}dxdt-{\int }_0^T{\int }_{\partial \Omega }(\sum _{k,l=1}^d{a}_{kl}(x)\frac{\partial u_j^1}{\partial x_l}{\nu }^k)wdxdt\end{array}$$

$$\text{ =−∫0T∫Λ∂νA∂uj1(x,t;q)h(x,t)dxdt},$$

$$\begin{array}{r}{\int }_0^T{}_0{D}_t^{\alpha }u(t)w(t)dt={\int }_0^{T}u(t){}_t{D}_T^{\alpha }w(t)dt,\text{for}u,w\in AC[0,T],\text{and}w(T)=0,\end{array}$$

$$\begin{array}{r}{\int }_0^T{\int }_{\Omega }{\varphi }_j(x{)}_0{D}_t^{\alpha }v(t)w(x,t;q)dxdt=-{\int }_0^T{\int }_{\Lambda }\frac{\partial u_j^2(x,t;q)}{\partial {\nu }_A}h(x,t)dxdt.\end{array}$$

$$\begin{array}{rl}{\int }_0^T{\int }_{\Omega }{\varphi }_j(x)v(t)w(x,t;p)dxdt& =-{\int }_0^T{\int }_{\Lambda }\frac{\partial u_j^1(x,t;p)}{\partial {\nu }_A}h(x,t)dxdt;\\ {\int }_0^T{\int }_{\Omega }{\varphi }_j(x{)}_0{D}_t^{\alpha }v(t)w(x,t;p)dxdt& =-{\int }_0^T{\int }_{\Lambda }\frac{\partial u_j^2(x,t;p)}{\partial {\nu }_A}h(x,t)dxdt.\end{array}$$

$$\begin{array}{rl}{\int }_0^T{\int }_{\Omega }{\varphi }_j(x)v(t)w(x,t;q)dxdt=& {\int }_0^T{\int }_{\Omega }{\varphi }_j(x)v(t)w(x,t;p)dxdt;\\ {\int }_0^T{\int }_{\Omega }{\varphi }_j(x){}_0{D}_t^{\alpha }v(t)w(x,t;q)dxdt=& {\int }_0^T{\int }_{\Omega }{\varphi }_j(x){}_0{D}_t^{\alpha }v(t)w(x,t;p)dxdt.\end{array}$$

$$\begin{array}{rl}{\int }_0^{T}v(t)w(x,t;q)dt=& {\int }_0^{T}v(t)w(x,t;p)dt;\\ {\int }_0^T{}_0{D}_t^{\alpha }v(t)w(x,t;q)dt=& {\int }_0^T{}_0{D}_t^{\alpha }v(t)w(x,t;p)dt.\end{array}$$

$$\begin{array}{r}{\int }_0^T{}_t{D}_T^{\alpha }w(x,t;q)vdt-{\int }_0^T\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial w(x,t;q)}{\partial x_l})vdt+{\int }_0^{T}q(x)w(x,t;q)vdt=0;\\ {\int }_0^T{}_t{D}_T^{\alpha }w(x,t;p)vdt-{\int }_0^T\sum _{k,l=1}^d\frac{\partial }{\partial x_k}(a_{kl}(x)\frac{\partial w(x,t;p)}{\partial x_l})vdt+{\int }_0^{T}p(x)w(x,t;p)vdt=0.\end{array}$$

$$\begin{array}{r}{\int }_0^T{}_0{D}_t^{\alpha }v[w(x,t;q)-w(x,t;p)]dt=(q-p){\int }_0^{T}w(x,t;p)v(t)dt,\end{array}$$

$$\begin{array}{r}(q-p){\int }_0^{T}w(x,t;p)v(t)dt=0.\end{array}$$

$$\begin{array}{rl}\frac{\partial u_j^1(x,t;p)}{\partial {\nu }_A}h(x,t)dxdt& =\frac{\partial u_j^1(x,t;q)}{\partial {\nu }_A}h(x,t)dxdt;\\ \frac{\partial u_j^2(x,t;p)}{\partial {\nu }_A}h(x,t)dxdt& =\frac{\partial u_j^2(x,t;q)}{\partial {\nu }_A}h(x,t)dxdt,j=1,2,\cdots ,\end{array}$$

$$\begin{array}{r}{B}^s(D):={}_0{H}^s(0,T;L^2(\Omega ))\cap L^2(0,T;H_0^1(\Omega )),s\in (0,1)\end{array}$$

$$\begin{array}{r}\parallel u{\parallel }_{B^s(D)}={\left(\parallel u{\parallel }_{H^s(0,T;L^2(\Omega ))}^2+\parallel u{\parallel }_{L^2(0,T;H_0^1(\Omega ))}^2\right)}^{\frac{1}{2}}.\end{array}$$

$$\begin{array}{r}A(u,w)=F(w),\forall w\in B^{\frac{\alpha }{2}}(D),\end{array}$$

$$\begin{array}{rl}A(u,w):=& ({}_0{D}_t^{\frac{\alpha }{2}}u,{}_t{D}_T^{\frac{\alpha }{2}}w{)}_{L^2(D)}+{\int }_D\sum _{k,l=1}^d{a}_{kl}(x)\frac{\partial u_j}{\partial x_k}\frac{\partial w}{\partial x_l}dxdt\\ & +(qu,w{)}_{L^2(D)};\\ F(w):=& (\varphi v,w{)}_{L^2(D)}.\end{array}$$

$$\begin{array}{r}\parallel u{\parallel }_{B^{\frac{\alpha }{2}}(D)}\le C\parallel \varphi v{\parallel }_{L^2(D)},\end{array}$$

$$\begin{array}{r}\mathcal{F}:q(x)\in Q\mapsto \Phi \in R^{2\times N},\end{array}$$

Full text

Identification of the degradation coefficient for an anomalous diffusion process in hydrology

Guang-Hui Zheng Corresponding author. College of Mathematics and Econometrics, Hunan University, Changsha 410082, Hunan Province, China. Email: zhenggh2012@hnu.edu.cn   

Ming-Hui Ding College of Mathematics and Econometrics, Hunan University, Changsha 410082, Hunan Province, China. Email: minghuiding@hnu.edu.cn

ABSTRACT

In hydrology, the degradation coefficient is one of the key parameters to describe the water quality change and to determine the water carrying capacity. This paper is devoted to identify the degradation coefficient in an anomalous diffusion process by using the average flux data at the accessible part of boundary. The main challenges in inverse degradation coefficient problems (IDCP) is the average flux measurement data only provide very limited information and cause the severe ill-posedness of IDCP. Firstly, we prove the average flux measurement data can uniquely determine the degradation coefficient. The existence and uniqueness of weak solution for the direct problem are established, and the Lipschitz continuity of the corresponding forward operator is also obtained. Secondly, to overcome the ill-posedness, we combine the variational regularization method with Laplace approximations (LA) to solve the IDCP. This hybrid method is essentially the combination of deterministic regularization method and stochastic method. Thus, it is able to calculate the minimizer (MAP point) more rapidly and accurately, but also enables captures the statistics information and quantifying the uncertainty of the solution. Furthermore, the existence, stability and convergence of the minimizer of the variational problem are proved. The convergence rate estimate between the LA posterior distribution and the actual posterior distribution in the sense of Hellinger distance is given, and the skewness are introduced for characterizing the symmetry or slope of LA solution, especially the relationship with the symmetry of the measurement data. Finally, the one-dimensional and two-dimensional numerical examples are presented to confirm the efficiency and robustness of the proposed method.

keywords: fractional diffusion equation, average flux data, degradation coefficient, Laplace approximations, Hellinger distance, skewness.

1 Introduction

In recent years, fractional calculus and fractional differential equations have been more and more extensively used in many scientific fields. For example, physical, chemical, biology, engineering, medicine, hydrology, finance and so on, refer to [1, 2, 3, 4, 5, 6, 7, 8, 9, 10].

As is known to all, in hydrology, the normal solute diffusion obeys Darcy s law:

[TABLE]

and mass conservation law:

[TABLE]

where $q$ is diffusion flux, $\kappa$ is diffusion coefficient, $u$ is concentration of solute and $f$ denotes some source or sink. By substituting (1.1) into (1.2), the following classical diffusion equation can be derived

[TABLE]

It is well known that the classical diffusion equation can describe the normal diffusion quite well, and in probability theory, it corresponds to Brownian motion. However, there are an increasing number of so-called anomalous diffusion arises in real world, especially the diffusion phenomena occurred in some media with memory and hereditary properties [1, 2, 3, 4, 5]. The anomalous diffusion is not consistent with the classical mass conservation law, but satisfies what is called time fractional mass conservation law [11, 12], i.e.

[TABLE]

where $\alpha\in(0,1)$ ,and ${}_{0}D_{t}^{\alpha}u$ is the Caputo fractional derivative defined by [1, 2, 3, 4]

[TABLE]

Similarly, by combining Darcy s law (1.1), one can obtain the time fractional diffusion equation (TFDE)

[TABLE]

The time fractional diffusion equation, compared with the classical Brownian motion, is closely related to fractional Brownian motion and gradually accepted as an important tool for describing anomalous diffusion. Particularly in hydrology, the TFDE models sticking and trapping between mobile periods for contaminant particles in a porous medium [13] or a river flow [14]. About the direct problems for TFDE, i.e., initial value problem and initial boundary value problem, which have been studied extensively in the past few years [15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. However, in some practical problems, the boundary data on the whole boundary cannot be obtained. We only know the noisy data on a part of the boundary or at some interior points of the concerned domain, which will lead to some inverse problems, i.e., fractional inverse diffusion problems. Recently, there are also rapidly growing publications on the time fractional inverse diffusion problems, such as a tutorial review [25], inverse initial boundary value problems [25, 26, 27, 28, 29, 30], inverse source problems [31, 32, 33, 34, 35], and inverse coefficient problems [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48].

In hydrology and water resources management, the pollutant degradation coefficient is one of the key parameters to describe the water quality change and to determine the water carrying capacity [49, 50]. So, identification of degradation coefficient is very important for water quality evaluation, monitoring and protection. In anomalous diffusion, identification of degradation coefficient is often correspond to the inverse reaction coefficient problems (IRCP) for TFDE.

As for the inverse reaction (or degradation) coefficient problems for TFDE, there are only limited papers studying this topic. For example, Tuan [41] gave the uniqueness of IRCP by only using finite measurement data on the boundary. Jin and Rundell [39] established the uniqueness in determining the reaction coefficient from the direct flux measurements for one dimensional TFDE, and an algorithm of the quasi-Newton type is proposed for the numerical reconstruction. In [48], Yamamoto and Zhang obtained conditional stability in recovering the reaction coefficient in a one dimensional TFDE with one half order Caputo derivative by a Carleman estimate. By means of integral transform method, Miller and Yamamoto [40] proved the uniqueness for the IRCP for TFDE from the internal measurement data. Li et al. [37] suggested an optimal perturbation algorithm for the simultaneous numerical recovery of the diffusion coefficient and fractional order in a one dimensional TFDE. In [44], Sun and Wei proved the uniqueness in identifying the reaction coefficient for one dimensional TFDE, and used the conjugate gradient method to solve it numerically. The IRCP for TFDE from the final time data measurement was discussed by Jin and Rundell in [25].

In this paper, based on pollutant degradation model in hydrology and compared with the above one dimensional case which was studied more extensively, we consider the inverse degradation coefficient problem for TFDE in higher dimensions. Let $\Omega$ be a bounded domain in $\mathds{R}^{d}$ $(d\geq 1)$ with smooth boundary $\partial\Omega$. The pollutant degradation model can be described by TFDE as follows

[TABLE]

where the degradation coefficient $q\in C({\overline{\Omega}})$, $q\geq 0$, the diffusion coefficient Matrix $A=(a_{kl})_{d\times d}\in C^{1}({\overline{\Omega}})^{{d\times d}}$ and satisfy $a_{kl}=a_{lk}$, $1\leq k,l\leq d$,

[TABLE]

The input source formed with separated variables $\phi(x)v(t)$, where $v(t)$ is the time-varying strength of source, and $\phi(x)$ denotes the space-position information. For example, the usual modelling of point source has the form $\phi(x)=\delta(x-a_{0})$, here $\delta(x)$ is Dirac function and $a_{0}$ is the source location.

Throughout this paper, The solutions to system (1.6) will be denoted by $u_{j}^{i}(x,t;q)$ in order to indicate its dependence on the degradation coefficient $q$, and correspond to the input sources $\phi_{j}(x)v_{i}(t)$, $i=1,2$; $j=1,2,\cdots$. Hereafter, $C$ refers to a generic constant which may differ at different occurrences. Then, our aim is to solve the following inverse problem.

Inverse degradation coefficient problem (IDCP) for TFDE: Given the input source $\phi_{j}(x)v_{i}(t)$, $i=1,2$; $j=1,2,\cdots$, such that the measurement data set

[TABLE]

determine the degradation coefficient $q$ (see Figure 1.1 for a schematic illustration). Here $\Lambda\subseteq\partial\Omega$ is an accessible part of the boundary, $h$ is a nonzero nonnegative function, and

[TABLE]

where $\nu_{l}(x)$ is the $l$th component of the outward unit normal vector $\nu(x)$.

The measurement data set (1.7) is a weighted integral data on accessible part of the boundary, the weight function $h$ can be interpreted as a characterization of measure instrument. Compared with the direct flux measurement (or Neumann measurement data) on boundary, i.e. \frac{\partial u_{j}^{i}(x,t;q)}{\partial\nu_{A}}\big{|}_{\Lambda\times(0,T)} (see Jin and Rundell [39]), it is a average flux measurement, which is rather easier to measure as a practical matter [51, 52]. However, because of the average effect, the amounts of average flux measurement data is much less than the amounts of direct flux measurement data. The similar situation also arises in the popular topic called inverse scattering problem with phaseless data [53, 54, 55]. For IDCP, owing to the reduced measurement data, the average flux measurement reconstruction is much more severely ill-posed than the direct flux measurement reconstruction. Moreover, the high nonlinearity is also the inherent difficulty for IDCP. In this paper, we will prove that the limited average flux measurement is enough to identify the degradation coefficient uniquely. Then we combine the variational regularization method with Laplace approximations to solve the IDCP. This hybrid method is essentially the organic combination of deterministic method and stochastic method.

Our contributions of this work are fourfold:

• •

We establish the uniqueness of IDCP, i.e., the average flux measurement data (1.7) can uniquely determine the degradation coefficient $q$.

• •

Due to the severe ill-posedness of IDCP, the reconstruction accuracy of degradation coefficient are very sensitive to the data noise, we combine the $L^{r}$ ($r>1$) variational regularization method and Laplace approximations theory to overcome this difficulty. Especially, it can effectively capture the statistics information and quantifying the uncertainty of the solution.

• •

We analyze the Hellinger distance between the exact posterior measure and its approximation given by Laplace approximations and prove the second order convergence rate at the MAP point.

• •

Since the average flux measurement only provide very limited information, such as the symmetry of degradation coefficient. We first use a term called ”skewness” from statistics to characterize the symmetry of solution, and show the relation between the symmetry of solution and the symmetry of data. Furthermore, we recover numerically the skewness of degradation coefficient by improving the symmetry of measurement data.

The paper is organized as follows. In Section 2, we present the uniqueness for the IDCP. In Section 3, we consider the well-posedness of weak solution of the direct problem, and prove the Lipschitz continuity of the corresponding forward operator. Furthermore, the variational regularization scheme with $L^{r}$ ($r>1$) penalty term is introduced to overcome the ill-posedness of IDCP. The existence, stability, convergence theorem of the minimizer of the variational problem are established. Moreover, the first-order Fréchet derivative and second-order Fréchet derivative of forward operator are obtained by using adjoint method, and the conjugate gradient method is presented for calculating the minimizer (MAP point). In Section 4, the Bayesian theory and Laplace approximations method are recalled. The convergence rate estimate between the LA posterior distribution and the actual posterior distribution in the sense of Hellinger distance is given, and the confidence region and skewness are introduced for characterizing the accuracy and reliability of LA. Several one-dimensional and two-dimensional numerical examples are shown in Section 5. Finally, the conclusions are given in Section 6.

2 The uniqueness of the IDCP

In this section, the classical solutions of the problem (1.6) which belong to the space $C(\overline{D})\cap W^{1}_{t}(0,T)\cap C^{2}_{x}(\Omega)$ are considered [17], here $W^{1}_{t}(0,T)$ denotes the space of functions $f\in C^{1}(0,T]$ such that $f^{\prime}\in L^{1}(0,T)$. Motivated by the ideas in [57], we show that the average flux measurement data can determine the degradation coefficient uniquely.

Theorem 1.1. Let $\{\phi_{j}\}_{j=1}^{\infty}\in C(\Omega)$ be a complete set in $L^{2}(\Omega)$, $v\in C^{1}(0,T)$ and $h\in C_{0}(\Lambda\times(0,T))$ be given nonzero nonnegative functions, and $v$ satisfies $v(0)=0$. Assume $p(x)$, $q(x)\in C(\overline{\Omega})$, $p$, $q\geq 0$ on $\Omega$. Let $u_{j}^{i}(x,t;p)$, $u_{j}^{i}(x,t;q)$ be the classical solutions of problem (1.6) corresponding to the input sources $\phi_{j}(x)v_{i}(t)$ ($i=1,2$; $j=1,2,\cdots$) with the degradation coefficients $p$ and $q$ respectively. If we choose $v_{1}=v$, $v_{2}={{}_{0}D_{t}^{\alpha}v}$ such that

[TABLE]

then $q=p$ in $\Omega$.

Proof. Notice that $h\in C_{0}(\Lambda\times(0,T))$ be a given nonzero nonnegative function. Then we set $h_{0}=h$ on $(\Lambda\times(0,T))$ and $h_{0}=0$ on $(\partial\Omega\setminus\Lambda)\times(0,T)$, and introduce the function $w(x,t;q)$ as the solution of the following adjoint problem

[TABLE]

In fact, by using the transform formula $\tilde{w}(x,t)=w(x,T-t)$, (2.9) becomes [48]

[TABLE]

Since we choose $v_{1}=v$, from $(\ref{mix-eq1})$, $(\ref{mix-eq-weak1})$ and the Green’s formula we compute

[TABLE]

[TABLE]

where the second equality uses the following fractional integration by parts formula (see [46] Lemma 2.1):

[TABLE]

here $AC[0,T]$ is the space of functions which are absolutely continuous on $[0,T]$. Similarly, by setting $v_{2}(t)={{}_{0}D_{t}^{\alpha}v}$, it follows that

[TABLE]

For $p(x)$ and the corresponding function $w(x,t;p)$ given by $(\ref{mix-eq-weak1})$ ,we see that

[TABLE]

Then we see from $(\ref{condtion 1})$ that

[TABLE]

By the completeness of $\{\phi_{j}(x)\}_{j=1}^{\infty}$, we obtain

[TABLE]

Multiplying equation (2.9) by $v$, integrating by parts over $(0,T)$, we find

[TABLE]

The two expressions of $(\ref{eq2})$ are subtracted from each other, and using (2.11) we have

[TABLE]

we finally combine $(\ref{eq3})$, $(\ref{eq8})$ to discover

[TABLE]

The maximum principle for time factional diffusion equation (2.9) or (2.10) (see [17], [56]) can be applied to deduce that $w(x,t;p)>0$, then $q=p$ in $\Omega$.

Remark 1.2. Obviously, when the average flux measurement (2.8) is replaced by the direct flux measurement and the other assumptions of Theorem 1.1 are satisfied, the uniqueness result still holds, i.e., if

[TABLE]

then $q=p$ in $\Omega$. In [39], the direct flux measurement was used by Jin and Rundell for coefficient identification in one dimensional TFDE (see Theorem 3.1 (a) in [39]). However, they assume that the reaction coefficient must be known beforehand in the neighborhood of right boundary $x=1$.

Remark 1.3. Let $\epsilon\in(0,T)$, $v\in H^{1}(0,T)$ satisfying $v=0$ in $(0,\epsilon]$, $v>0$ in $(\epsilon,T)$, and the other assumptions of Theorem 1.1 are satisfied. By choosing $v_{1}=v$, $v_{2}={{}_{0}D_{t}^{\alpha}v}$, we see the uniqueness also holds. The proof is the same as in Theorem 1.1.

3 Variational regularization method

3.1 Weak solution of TFDE

In order to obtain the continuity of forward map and provide feasibility for numerical computation (such as finite element method), we study the weak solution of problem (1.6) with $q(x)\in L^{\infty}(\Omega)$. The existence, uniqueness and stability of weak solution are given. Let ${}_{0}C^{\infty}(0,T)$ denote the space of infinitely differentiable functions on $(0,T)$ with compact support in $(0,T]$. The Sobolev space ${}_{0}H^{s}(0,T)$ is the closure of ${}_{0}C^{\infty}(0,T)$ with respect to the norm $\parallel\cdot\parallel_{H^{s}(0,T)}$, where $\parallel\cdot\parallel_{H^{s}(0,T)}$ denotes the norm in the usual fractional Sobolev spaces $H^{s}(0,T)$ [58, 72].

We recall $D=(\Omega\times(0,T))$ and define a Hilbert space

[TABLE]

equipped with the norm

[TABLE]

Definition 3.1.1. (Weak solution) We call that $u\in B^{\frac{\alpha}{2}}(D)$ is a weak solution of the initial boundary value problem $(\ref{mix-eq1})$ provide

[TABLE]

where the bilinear form $A(\cdot,\cdot)$ and $F(\cdot)$ are defined by

[TABLE]

Theorem 3.1.2. If $q\in L^{\infty}(\Omega)$, $q>0$, and $a_{kl}(x)$, $\phi$, $v$ satisfy the same assumption in Theorem 1.1. Then the problem $(\ref{mix-eq1})$ exists a unique solution $u$ in $B^{\frac{\alpha}{2}}(D)$ and the solution $u$ satisfies

[TABLE]

where $C$ is a constant independent of $u$.

By using the Lax-Milgram theorem, the proof of Theorem 2.1 is standard (refer to [58, 19, 46]). Hence, we omit it.

3.2 The forward operator

In this section, we solve numerically the reaction coefficient $q(x)$ by problem $(\ref{mix-eq1})$. The inverse coefficient problem is formulated into a variational problem by using the Tikhonov regularization. Then the existence, stability and convergence of minimizer for the variational problem are provided.

Define a forward or solution operator

[TABLE]

where $Q=\{{q(x)\in L^{\infty}(\Omega)}\mid q_{\textrm{min}}\leq q(x)\leq q_{\textrm{max}}\}$, and

[TABLE]

and the norm of $R^{2\times N}$ is defined by $\parallel\Phi\parallel_{s}=(\sum\limits_{i=1}^{2}\sum\limits_{j=1}^{N}\mid\varphi_{ij}\mid^{s})^{\frac{1}{s}},\ s\geq 1$. To get the continuity of the forward operator, we restrict $h\in L^{2}(0,T;H_{0}^{\frac{1}{2}}(\Lambda))$.

Theorem 3.2.1. If $h\in L^{2}(0,T;H_{0}^{\frac{1}{2}}(\Lambda))$, and the other assumption is the same as Theorem 1.1. Then the nonlinear forward map $\mathcal{F}:Q\mapsto R^{2\times N}$ is the Lipschitz continuous.

Proof. Setting $u_{j}^{i}=u_{j}^{i}(x,t;q),\ \tilde{u}_{j}^{i}=u_{j}^{i}(x,t;\tilde{q})$, it follows that

[TABLE]

and

[TABLE]

Let $z_{j}^{i}=u_{j}^{i}-\tilde{u}_{j}^{i}$, then $z_{j}^{i}$ satisfies the following equation

[TABLE]

Since $v_{i}(t)\in H^{1}[0,T];\ \phi_{j}\in L^{2}(\Omega)$ and $v_{i}(t)\phi_{j}(x)\in H^{1}(0,T;L^{2}(\Omega))\subseteq L^{2}(\Omega\times(0,T))$, by using Theorem 3.1.2, we have

[TABLE]

Similarly from $(\ref{eq5})$

[TABLE]

Furthermore, by trace theorem, it implies

[TABLE]

and then

[TABLE]

Therefore, the Lipschitz continuity of $\mathcal{F}$ is proved.

3.3 Variational regularization method with $L^{r}$ ($r>1$) penalty term

In order to overcome the ill-posedness of the problem, we apply the variational regularization method with $L^{r}$ ($r>1$) penalty term to deal with it. Then the corresponding variational functional is defined as follows

[TABLE]

where $\Phi^{\delta}=(\varphi_{ij}^{\delta})_{2\times N}$ denotes the measure data matrix and satisfies $\parallel\Phi^{\delta}-\Phi\parallel_{s}\leq\delta$, here $\delta$ is the noise level. $\mu$ is the regularization parameter. Next, we will prove the existence, stability and convergence of minimizer of the variational functional (3.23). Although the proof is similar to [44, 59], for the completeness, we give some details.

Theorem 3.3.1. There exists a minimizer $q_{\mu}^{\delta}\in Q$ for variation functional $J(q)$.

Proof. Because of the nonnegativity of $J(q)$, there is a minimizing sequence $\{q_{k}\}$ in $Q$ such that

[TABLE]

Since $\{q_{k}\}\subseteq Q$ is bounded in $L^{\infty}(\Omega)$, there exists a subsequence, which is again denoted by $\{q_{k}\}$, such that $q_{k}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}q_{0}$ in $L^{\infty}(\Omega)$. Moreover, by using the reflexivity of $L^{r}(\Omega)$ and the density of $L^{r}(\Omega)$ $(r>1)$ in $L^{1}(\Omega)$ [60], it follows that $q_{k}\rightharpoonup q_{0}$ in $L^{r}(\Omega)$. Owing to the closed convexity of $Q$, by using Mazur Theorem, we find $q_{k}\rightarrow q_{0}$, and then $q_{0}\in Q$. From the weak lower semicontinuity of norm, and the Lipschitz continuity of $\mathcal{F}$, we have

[TABLE]

Therefore, $q_{0}$ is a minimizer of $J(q)$.

Theorem 3.3.2. Assume that $\{\Phi_{k}\}$ is a sequences which satisfy $\Phi_{k}\rightarrow\Phi^{\delta},\ k\rightarrow\infty$ in $R^{2\times N}$, and $\{q_{k}\}$ is a minimizer of $J(q)$ with $\Phi^{\delta}$ replaced by $\Phi_{k}$. Then the minimizers of $J(q)$ are stable with respect to the measurement data $\Phi^{\delta}$.

Proof. From the definition of $\{q_{k}\}$, we find

[TABLE]

Since $\{q_{k}\}\subseteq Q$ is bounded in $L^{\infty}(\Omega)$, there has a subsequence, still denoted by $\{q_{k}\}$, such that $q_{k}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}q_{0}$. Similar to the proof of Theorem 3.3.1, it implies $q_{k}\rightarrow q_{0}$. Based on continuity of $\mathcal{F}$ and weak lower semicontiunity of norm, by $(\ref{eq6})$ it follows that

[TABLE]

for all $q\in Q$. This deduce that $q_{0}$ is a minimizer of $J(q)$. Furthermore, we set $q=q_{0}$, and find that

[TABLE]

Thus, the stability result holds.

Definition 3.3.3. $q^{\dagger}\in Q$ is called an $L^{r}$-minimizing solution if

[TABLE]

Theorem 3.3.4. Suppose that the noise level sequence $\{\delta_{k}\}$ convergence monotonically to [math], and the corresponding measurement data $\Phi^{\delta_{k}}$ satisfy $\parallel\Phi^{\delta_{k}}-\Phi\parallel_{s}\leq\delta_{k}$. Moreover, assume that the regularization parameter $\mu(\delta)$ satisfies $\mu(\delta)\rightarrow 0$, and $\frac{\delta^{s}}{\mu(\delta)}\rightarrow 0$, (as $\delta\rightarrow 0$), and $\alpha(\delta)$ is also monotonically increasing. $\{q_{\mu(\delta_{k})}^{\delta_{k}}\}$ is a minimizer of $J(q)$ with $\Phi^{\delta}$ replaced by $\Phi^{\delta_{k}}$. Then $\{q_{\mu(\delta_{k})}^{\delta_{k}}\}$ has a subsequence which convergent to an $L^{r}$-minimizing solution, and satisfies that $\parallel q^{\dagger}\parallel_{L^{r}(\Omega)}=\lim\limits_{k\rightarrow\infty}\parallel q_{\mu(\delta_{k})}^{\delta_{k}}\parallel_{L^{r}(\Omega)}$.

Proof. From the definition of $q_{\mu(\delta_{k})}^{\delta_{k}}$, we get

[TABLE]

Then, take $k\rightarrow\infty$, we obtain

[TABLE]

Moreover, notice that

[TABLE]

and $\frac{\delta^{s}}{\mu(\delta)}\rightarrow 0,\ (k\rightarrow\infty)$, it implies that

[TABLE]

Since $\{q_{\mu(\delta_{k})}^{\delta_{k}}\}\subseteq Q$ is bounded, there has a subsequence, which denoted again by $\{q_{\mu(\delta_{k})}^{\delta_{k}}\}$ satisfying $q_{\mu(\delta_{k})}^{\delta_{k}}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}q_{0}\in L^{\infty}(\Omega)$. Likewise, by Mazur Theorem, we also get $q_{\mu(\delta_{k})}^{\delta_{k}}\rightarrow q_{0}\in L^{\infty}(\Omega)$. Hence, by (3.27), a similar argument implies that

[TABLE]

for all $q\in Q$ satisfying $\mathcal{F}(q)=\Phi$. Setting $q=q_{0}$ deduce that $\parallel q^{\dagger}\parallel_{L^{r}(\Omega)}=\parallel q_{0}\parallel_{L^{r}(\Omega)}$. That is, $q_{0}$ is an $L^{r}$-minimizing solution, and satisfies $\parallel q^{\dagger}\parallel_{L^{r}(\Omega)}=\parallel q_{0}\parallel_{L^{r}(\Omega)}=\lim\limits_{k\rightarrow\infty}\parallel q_{\mu(\delta_{k})}^{\delta_{k}}\parallel_{L^{r}(\Omega)}$.

3.4 The Fréchet derivative of variation functional

For simplicity, we only focus on the case $s=r=2$, then the variation functional becomes

[TABLE]

In order to find the minimizer of variation functional (3.29), the efficient evaluation of the Fréchet derivative is critical, here we adopt the adjoint method. First, we claim that the following asymptotic expansion formula of solution holds.

Theorem 3.4.1. The solution $u_{j}^{i}(q)$ is differentiable in the sense that: for any direction $\delta q\in L^{\infty}(\Omega)$, we have that

[TABLE]

where $\vartheta_{j}^{i}(q)[\delta q]$ satisfies the sensitive equation:

[TABLE]

Proof. Setting $\tilde{\vartheta}_{j}^{i}=u_{j}^{i}(q+\delta q)-u_{j}^{i}(q)-\vartheta_{j}^{i}(q)[\delta q]$. It is easy to verify that $\tilde{\vartheta}_{j}^{i}$ satisfies

[TABLE]

By Theorem 3.1.2, it deduces that $\|\tilde{\vartheta}_{j}^{i}\|_{B^{\frac{\alpha}{2}}(\Omega)}\leq\|\delta q\|_{L^{\infty}(\Omega)}\|\vartheta_{j}^{i}(q)[\delta q]\|_{L^{2}(\Omega)}$. Similarly, by using Theorem 3.1.2 to (3.31), we find $\|\vartheta_{j}^{i}(q)[\delta q]\|_{B^{\frac{\alpha}{2}}(\Omega)}\leq\|\delta q\|_{L^{\infty}(\Omega)}\|u_{j}^{i}(q)\|_{L^{2}(\Omega)}$, and then the claim now holds.

The following theorem shows that the solution $u_{j}^{i}(q)$ is twice Fréchet differentiable. Since the proof is analogous to Theorem 3.4.1 and is omitted.

Theorem 3.4.2. $\vartheta_{j}^{i}(q)[\delta q]$ is differentiable in the sense that: for the direction $\tilde{\delta}q\in L^{\infty}(\Omega)$, we have that

[TABLE]

where $\zeta_{j}^{i}(q)[\delta q,\tilde{\delta}q]$ satisfies the second-order sensitive equation:

[TABLE]

Theorem 3.4.3. For the direction $\delta q\in L^{\infty}(\Omega)$, the Fréchet derivative of variation functional $J(q)$ at $q\in Q$ is given by

[TABLE]

where $\varpi_{j}^{i}$ satisfies the adjoint equation:

[TABLE]

Proof. In fact, we only need to prove that the Fréchet derivative of $J_{1}(q)=\frac{1}{2}\parallel\mathcal{F}(q)-\Phi^{\delta}\parallel_{2}^{2}$ is given by

[TABLE]

Notice that $J_{1}(q)$ can be written as the Lagrangian multiplier formula,

[TABLE]

for any multiplier function $\lambda\in L^{2}(D)$. Then, from Theorem 3.4.2 and integration by parts, we have that

[TABLE]

By choosing $\lambda=\varpi_{j}^{i}$, and $\varpi_{j}^{i}$ satisfies (3.35), then we obtain (3.36).

3.5 The conjugate gradient method

The conjugate gradient method (CGM) combined with an appropriate stopping rule can serve as a regularization method, and it has been applied to various inverse problems [44, 46, 71]. However, the CGM is a deterministic regularization method which yield only a point estimate of the solution, without quantifying the associated uncertainties of measurement data noise. Here, we use CGM to calculate a maximum a posteriori (MAP) estimator of the solution, and then in conjunction with the Laplace approximations to sample the associated posterior distribution. The CGM we utilized is presented in Algorithm 1.

4 Bayesian theory and Laplace approximation

4.1 Convergence of Laplace approximation

From classical Bayesian theory, and notice that the observation error $\eta$ is an independent and identically distributed (i.i.d) Gauss random vector with mean zero and the covariance matrix $B=\delta^{2}I$, here $I$ is the unit matrix, we can write the minimization functional $(\ref{3.29})$ as [61]

[TABLE]

where $\parallel\cdot\parallel_{B}$ is a covariance weighted norm on $R^{2N\times 1}$ given by $\parallel\cdot\parallel_{B}=\parallel\delta^{-1}I\cdot\parallel_{2}$, and $\Phi^{\delta}=(\varphi_{11}^{\delta},\cdots,\varphi_{1N}^{\delta},\varphi_{21}^{\delta},\cdots,\varphi_{2N}^{\delta})^{T}\in R^{2N\times 1}$ is a stretched vector version of the one defined in $(\ref{3.29})$. Furthermore, the minimizer of $J_{B}$ defines the maximum a posteriori estimator

[TABLE]

The Laplace approximation (LA) in essence is a linearization around the MAP point $q_{\mathrm{MAP}}$ for sampling the posterior distribution of the solution (refer to [62]). It consists of approximating the posterior measure(or distribution) by $\omega\thickapprox N(q_{\mathrm{MAP}},C_{\mathrm{MAP}})$, here $C_{\mathrm{MAP}}=(J_{B}^{\prime\prime}(q_{\mathrm{MAP}}))^{-1}$ is the inverse of Hessian of $J_{B}(q_{\mathrm{MAP}})$. Next, motivated by [63], we analyze the Hellinger distance between the exact posterior measure and its approximation given by $N(q_{\mathrm{MAP}},C_{\mathrm{MAP}})$ and obtain the convergence as $q\rightarrow q_{\mathrm{MAP}}$. In [64], Wacher gave a bound of the Hellinger distance between the posterior and its LA, but it seems not to deduce the convergence. Moreover, the estimation of other distances (such as Kullback-Leibler divergence) between the posterior distribution and its approximation were applied in analysing the stochastic surrogate models (e.g., [65, 66]).

Lemma 5.1.1. For every $\varepsilon>0$ and $q\in Q$, there exists $M\in R$ such that, the forward map $\mathcal{F}:Q\mapsto R^{2N\times 1}$ satisfies,

[TABLE]

Proof. From theorem 3.1.2, and trace theorem, it follows that

[TABLE]

Thus

[TABLE]

and then

[TABLE]

here, $M=\log C_{N,\delta}$.

From the Laplace approximation, the Taylor expansion of $J_{B}(q)$ at MAP point $q_{\mathrm{MAP}}$ is given by

[TABLE]

In the infinite dimensional version of Beyesian theory. The posterior measure is absolutely continuous with respect to the prior, and the Rodon-Nikodym derivative between them is determined by the data likelihood, i.e.

[TABLE]

Similarly, the Laplace approximation version of Bayesian formula is defined as

[TABLE]

Lemma 5.1.2. The forward map $\mathcal{F}:Q\rightarrow R^{2N\times 1}$ satisfies:

(i) for every $\varepsilon>0$ and $r>0$, there exsits $M=M(\varepsilon,r)\in R$ such that, for all $q\in Q$ and $\parallel\Phi^{\delta}\parallel_{B}<r$,

[TABLE]

(ii) for every $r>0$, there exsits a $L=L(r)>0$ such that for all $q\in Q$ and $\Phi^{\delta}\in R^{2N\times 1}$ with $\max\{\parallel q\parallel_{L^{\infty}(\Omega)},\parallel\Phi^{\delta}\parallel_{B}\}<r$,

[TABLE]

The Lemma 5.1.2 is the direct results of Lemma 5.1.1, by using Lemma 2.1 in [63], so we omit the proof.

Theorem 5.1.3. Assume that $L(r)$ and $M$ are defined in Lemma 5.1.2, then the measure $w$ and its Laplace approximation measure $\tilde{w}$ are close with respect to the Hellinger distance, i.e., there is a constant $C$, such that

[TABLE]

and

[TABLE]

where the Hellinger distance is defined by

[TABLE]

Proof. By using Lemma 5.1.2 (ii), we obtain

[TABLE]

where $\mid\Omega\mid=\int_{\Omega}1dx$, $\mid Q_{r}\mid_{\infty}=\int_{L^{\infty}(\Omega)}1dw_{0}(q)$, and

[TABLE]

From the second-order necessary condition for minimizer, lemma 5.1.2 (i) and Fernique theorem (see [63]), we see that

[TABLE]

Moreover, notice $(\ref{eq16})$, it follows that

[TABLE]

From the definition of Hellinger distance, we get

[TABLE]

Now, using $(\ref{eq16})$ and $(\ref{eq17})$, it deduces

[TABLE]

and then

[TABLE]

Hence, we can find that

[TABLE]

Remark 5.1.4. From the proof of Theorem 5.1.3, we can see that the approximation error is derived from the Taylor expansion (4.40), and second-order convergence rate is given. Moreover, from the convergence estimation (4.45), we also find that when the noise level $\delta$ is decreased or the regularization parameter $\mu$ is increased, the Hellinger distance $d_{Hell}(w,\tilde{w})$ become small, and it provides some inspiration for selecting the regularization parameter.

4.2 Numerical algorithm of Laplace approximation

For effective numerical simulation, we take a finite dimensional approximation of the previous minimization problem as follows

[TABLE]

where $q^{M}=(q_{1},...,q_{M}),\ B=\delta^{2}I,\ B_{\mu}=\frac{\delta^{2}}{\mu}I$, and $\mathcal{\tilde{F}}$ is a finite-dimensional approximation of the continuous forward map $\mathcal{F}$. The Laplace approximation theory shows the posterior distribution $w\approx N(q_{\mathrm{MAP}},C_{\mathrm{MAP}})$. In the finite dimensional, the covariance matrix (see [67], Section 10.5)

[TABLE]

where $P$ is Jacobian matrix of forward operator $\tilde{\mathcal{F}}$ at $q^{M}$ point. Notice that the covariance formula (4.53) only uses the first order derivatives of $\tilde{\mathcal{F}}$. A standard implementation of Laplace approximation is presented in the following algorithm [67]:

Samples generated by (4.55) are drawn from $N(q_{\mathrm{MAP}},C_{\mathrm{MAP}})$, and so the ensemble of $N_{e}$ realizations $\{q^{j}\}_{j=1}^{N_{e}}$ provides an approximation to $N(q_{\mathrm{MAP}},C_{\mathrm{MAP}})$ and hence the posterior. Finally, we use the mean of the sample $\bar{q}_{n}=\frac{1}{n}\sum\limits_{j=1}^{n}q^{j}$ as an approximation of $q_{\mathrm{MAP}}$, where the convergence of $\bar{q}_{n}$ follows from the strong law of large numbers. From the classical Gaussian statistic theory, we find that $\bar{q}_{n}$ are consistent and best unbiased estimate of $q_{\mathrm{MAP}}$.

4.3 Confidence region and skewness

We can calculate the confidence region for the inferred parameters (such as degradation coefficient in IDCP) by the LA posterior samples. The confidence region is a set of points in an $d$-dimensional space, often represented as an ellipsoid around a point which is an estimated parameter. The confidence region shows that the real value of the identified parameter has a certain probability of falling around the numerical construction result and quantifies the level of confidence that the parameter lies in the region. Therefore, the confidence region gives the reliability of the construction result of the inferred parameters.

Definiton 4.3.1. (Confidence region) [68] The confidence region of $(1-\alpha)\times 100\%$ is defined as follows

[TABLE]

It is an ellipsoid centered on $q_{\mathrm{MAP}}$ point, and assume the eigenvalue of $C_{\mathrm{MAP}}$ are

[TABLE]

We find that the axes of the ellipsoid are $\sqrt{\frac{\lambda_{j}\mathrm{\chi^{2}_{\alpha}}(M)}{n}}$. Since $C_{\mathrm{MAP}}=\delta^{2}(\mu I+P^{T}P)^{-1}$, it deduces that $\lambda_{j}(C_{\mathrm{MAP}})=\frac{\delta^{2}}{\mu+\lambda_{j}(P^{T}P)}$, and the length of each axis is $\delta\sqrt{\frac{\chi_{\alpha}^{2}(M)}{n(\mu+\lambda_{j}(P^{T}P))}}$. Let $\bar{q}_{n,j}$ be the component of $\bar{q}_{n}$, then it satisfies $\mid\bar{q}_{n,j}-\bar{q}_{\mathrm{MAP},j}\mid\leq\delta\sqrt{\frac{\chi_{\alpha}^{2}(M)}{n(\mu+\lambda_{j}(P^{T}P))}}$. Thus, we can see that the confidence region is a multi-dimensional generalization of a confidence interval. It also can be seen that when the number of samples $n$, the regularization parameter $\mu$ and the eigenvalue of $P^{T}P$ increase or the measurement noise level $\delta$ decrease, the size of confidence region will decrease and give higher reliability.

Skewness is a statistic that studies the symmetry of data distribution. By measuring the skewness, we can determine the degree and direction of the asymmetry of the data distribution. The definition of skewness is given below

Definition 4.3.2. (Skewness) [69] Assume the third central moment of the random variable $X$ is exists, the following ratios:

[TABLE]

is called the skewness coefficient for X, i.e., skewness. Here $\varsigma$ is the mean, $\sigma$ is the standard deviation, $E$ is the expectation operator, $\gamma_{3}$ is the third central moment.

By normalization we can transform the identified parameter into a well defined probability density function for a hypothetical random variable, and calculate the skewness of random variable to estimate the symmetry of the inferred parameter (see Section 5).

5 Numerical experiments and discussions

In this section, we present some numerical examples to illustrate the feasibility of the Laplace approximation method for IDCP.

The noisy data are generated by adding random perturbations as follows [71]

[TABLE]

here $\zeta$ is a Gaussian random variable with zero mean and unit standard deviation. $\varepsilon$ indicates the noise level.

To show the accuracy of numerical solutions, we compute the relative $L^{2}$ error denoted by

[TABLE]

where $\bar{q}_{n}(x)$ is the approximate degradation coefficient reconstructed by the LA algorithm, and $q(x)$ is the exact solution.

In an iteration algorithm, aim to make $q_{k}$ to approximate $q_{\mathrm{MAP}}$ point, the important work is to find a suitable stopping rule. Using $E_{k}$ to represent the residuals of iterations of $k-1$ and $k$ steps, i.e.,

[TABLE]

when $E_{k}\leq eps$, we stop iterating as a stopping criterion.

5.1 Inversion for one-dimensional degradation coefficient

The domain $\Omega$ under consideration is a unite line $[0,1]$, and the boundary is indicated as $\Lambda_{0}=\{x=0\}$, $\Lambda_{1}=\{x=1\}$. The direct problem $(\ref{mix-eq1})$ is discretized by using $300$ uniform rectangular finite element. Set $T=1$, and the grid point on $[0,T]$ is $101$. The boundary data are obtained by solving the direct problem $(\ref{mix-eq1})$ with $a_{kl}(x)=1$, and setting positive function $h(x,t)=(1-t)$. We solve the direct problem, sensitive problem and adjoint problem by using the finite element method and construction of data $\varphi_{ij}^{\delta}$ using difference method.

The $q_{k}$ obtained by conjugate gradient method is an approximation to $q_{\mathrm{MAP}}$. Samples generated by algorithm 2 are drawn from $N(q_{\mathrm{MAP}},\ C_{\mathrm{MAP}})$, and so the ensemble $\{q^{(j)}\}_{j=1}^{Ne}$, the length $Ne$ is taken to be $10000$, provides an approximation to $N(q_{\mathrm{MAP}},\ C_{\mathrm{MAP}})$.

5.1.1 Smooth solution

$\mathbf{Example~{}~{}1.}$ In $(\ref{mix-eq1})$, the accessible boundary $\Lambda$ is taken to be $\Lambda_{1}$, and the degradation coefficient is given by

[TABLE]

First we investigate the effect of the amounts of basis functions. Trigonometric basis functions are used, i.e., $\phi_{j}\in span\{1,cos(2\pi x),sin(2\pi x),\cdots,cos(2N\pi x),sin(2N\pi x)\}$ in Example 1. In Table 1, for fractional order $\alpha=0.3$, we list the $L^{2}$ error $r_{e}$ of LA solution and exact solution under different amounts of basis functions and different noise levels, $N$ is the number of basis functions. From the Table 1, we can see that the error becomes smaller at the same noise level as the number of basis function increases. When the number of basis functions is the same, the result become worse with the increase of noise levels. In the following calculations, we chose $N=5$.

In Theorem 3.3.4, we introduce a slightly crude rule of regularization parameter selection, i.e.,

[TABLE]

The numerical results for Example 1 with various $\mu$ and $\varepsilon$ are shown in the Table 2. We find that at the same noise level, $\mu=\delta$ was selected to meet the rule $(\ref{eqq8})$, but the results are unsatisfactory. By choosing $\mu=\delta^{\frac{3}{2}}$ which satisfies $(\ref{eqq8})$, the error between the LA solution and the exact solution is small, and the result was desirable. However, the crude rule $(\ref{eqq8})$ for regularization parameter selection is not satisfied with $\mu=\delta^{2}$, and the error between LA solution and exact solution will gradually increase and result will be shock. Thus, the regularization parameter can be chosen neither too large nor too small, even if it satisfies $(\ref{eqq8})$.

The numerical result for Example 1 for various noise levels $\varepsilon=0.0001,0.0005,0.001$, and in the case of $\alpha=0.3,0.7$ are shown in Figure 5.2. The basis function is the trigonometric basis function in $L^{2}(\Omega)$. Compare the Figure $\ref{Example.1 for various noise le}(a)$ with Figure $\ref{Example.1 for various noise le}(b)$, the results in Figure $\ref{Example.1 for various noise le}(b)$ is worse than Figure $\ref{Example.1 for various noise le}(a)$. In Table 3, we further show the numerical errors $r_{e}$ of the Example 1 for different $\alpha$ and $\varepsilon$. It can be seen that the numerical results become worse as the noise levels increase. The error become larger as fractional order $\alpha$ increases. However, when the noise level is 0.0001, the numerical result is insensitive to fractional order $\alpha$.

Next, we change trigonometric basis function into polynomial basis function in $L^{2}(\Omega)$, i.e., $\phi_{j}\in span\{1,x,x^{2},\cdots,x^{N}\}$. We show the reconstruction results for Example 1 under various error levels and different types of basis function with $\alpha=0.3$ in the Table 4. The data indicate that the difference in types of basis function will affect the numerical results. The numerical results obtained by trigonometric basis functions are better than those obtained by polynomial basis functions. In the following example, we only consider the triangle basis functions.

In the Figure 5.3 show that the results by using the direct flux data is not sensitive to the fractional order $\alpha$. From Table 5, we find that the reconstruction results by using the direct flux data are better than the ones by using the average flux data. However, we also get satisfactory numerical results by using the average flux data when the noise levels are $\varepsilon=0.0001,\ 0.0005,\ 0.001$ respectively. The reconstruction error caused by using average flux data increases sharply when the noise level exceeds $0.001$, but using the direct flux data still give good results when the noise level exceeds $0.005$ and even reaches $\varepsilon=0.05$. This is because, compared to direct flux data, the amount of the average flux data is less and provides limited information. Moreover, the limited measurement data lead to higher sensitivity to noise and severally ill-posedness of IDCP. But, the average flux data is rather easier to measure as a practical matter, and it has been widely used recently. Hence, recovering the degradation coefficient accurately by using limited measurement data can be a real challenge.

For Example 1, we draw the 95% confidence interval in Figure 5.4 for the noise level $\varepsilon=0.0001,\ 0.0005$ with $\alpha=0.3$. The posterior mean $\bar{q}_{n}(x)$ is in excellent agreement with the exact solution, and the confidence interval quantifies its associated uncertainty. The confidence interval shrinks as the noise level $\varepsilon$ decreases. We observe that the confidence interval on one side near the observation data is relatively narrow, and the corresponding confidence interval is also relatively accurate. The confidence interval far from the observation data is wide and the corresponding confidence interval is relatively inaccurate. It’s a surprise that the position of the data affects the accuracy of our reconstruction results. In Example 2, Example 3, we will continue to verify this result.

$\mathbf{Example~{}~{}2}$. In $(\ref{mix-eq1})$, the accessible boundary $\Lambda$ is taken to be $\Lambda_{0}$, and the degradation coefficient is given by

[TABLE]

The numerical results for Example 2 for various levels of noise in the data are shown in Figure 5.5 with $\alpha=0.3,\ 0.7$. The results in Figure $\ref{tu-left}(b)$ is worse than Figure $\ref{tu-left}(a)$ as the fractional order $\alpha$ increase. Moreover, we also find that the position of the measurement data have a great influence on the reconstruction results. In the Figure 5.2, we let $\Lambda=\Lambda_{1}$, the average flux data are measured on the $\Lambda_{1}$ for Example 1, when the exact solution $q$ is biased to the right of the region, we obtain good numerical results. In the Figure 5.5, setting $\Lambda=\Lambda_{0}$, when the average flux data are measured on the $\Lambda_{0}$ for Example 2, the left biased $q$ give desired approximate result (see Example 3 for further discussion). We used LA algorithm to sample, calculated the corresponding posterior mean value $\bar{q}_{n}$, and drew 95% confidence interval in the Figure 5.6, we can get similar results of Example 1. The confidence interval shrinks as the noise level $\varepsilon$ decreases and near observation data is more accurate, far from the observation data is wide and inaccurate, see the Figure 5.6.

$\mathbf{Example~{}~{}3}$. In $(\ref{mix-eq1})$, the accessible boundary $\Lambda$ is taken to be $\Lambda=\Lambda_{0}\cup\Lambda_{1}$, and the degradation coefficient is given by

[TABLE]

The numerical results for Example 3 with various fractional order and noise levels are presented in Figure 5.7, it shows that the fractional order $\alpha$ increases but the result is not good. When we measure average flux data from both sides, and use LA method to recover the degradation coefficient, the error is small with noise levels are $\varepsilon=0.0001,\ 0.0005,\ 0.001$ respectively. Furthermore, the confidence interval has higher precision on both sides and lower precision in the middle of the region. The confidence interval are shown in Figure 5.8. This is consistent with our previous conclusions.

Now, we introduce the concept of skewness to describe the degree of LA solution’s deviation from the mean value, and further demonstrate the relation between the symmetry of solution and the symmetry of measurement data. The degradation coefficient is taken as Example 3, and the boundary average flux data are measured on different positions. Then the corresponding numerical results for Example 3 are displayed in Figure 5.9. As the error level increase to 0.001, when measuring the data on the boundary of $\Lambda_{1}$, the degradation coefficient recovered by LA method is left skewed, i.e., its skewness is positive. When the average flux data are measured on the boundary of $\Lambda_{0}$, the degradation coefficient is right skewed, i.e., its skewness is negative. When we measure average flux data on both sides of the boundary (symmetric data), that is to say $\Lambda=\Lambda_{1}\cup\Lambda_{0}$, the reconstructed degradation coefficient is also symmetric and give a accurate approximation, and its skewness is zero. Thus, we can use the symmetry of measurement data to capture the symmetry feature of the identified object. By using (4.58), the calculation of the corresponding skewness can be referred to Table 6.

5.1.2 Nonsmooth function

$\mathbf{Example~{}~{}4}.$ In this example, we consider the more challenging case of reconstructing a nonsmooth example with a cusp, and the degradation coefficient is prescribed as follows:

[TABLE]

$\mathbf{Example~{}~{}5}.$ We also consider an discontinuous example, and the degradation coefficient is:

[TABLE]

The numerical results for Example 4 and Example 5 for various noise levels in the case of $\alpha=0.3$ are shown in the Figure 5.10 with the regularization parameter $\mu=\delta^{\frac{3}{2}}$. It can be seen that the smaller the noise level, the better the numerical results. The error for nonsmooth Example 4 in the neighborhood of the cusp is large because of the smoothing nature of the prior. And the same result can be obtained from Example 5, the numerical results obtained in the neighborhood of segment point are not very desired. For the nonsmooth case, different regularization methods and different prior information are needed, such as TV prior [70], and we do not discuss the details here.

5.2 Inversion for two-dimensional degradation coefficient

The domain $\Omega$ under consideration is a unit square $[0,1]\times[0,1]$. Set $T=1$, and the grid point on $[0,T]$ is $101$. In the equation $(\ref{mix-eq1})$, we take $d=2$, the diffusion coefficient matrix is unitary. The forward problem is discretized using $2500$ uniform rectangular finite element. The number of basis functions is take N=5 and sample size $Ne=10000$. Set positive function $h(x,t)=(t-1)$. The $\partial\Omega=\Lambda_{1}\cup\Lambda_{2}\cup\Lambda_{3}\cup\Lambda_{4}$, where $\Lambda_{1}=(0,1]\times\{0\}$, $\Lambda_{2}=\{1\}\times(0,1]$, $\Lambda_{3}=[0,1)\times\{1\}$, $\Lambda_{4}=\{0\}\times(0,1)$. Next, we presents the numerical results for 2D cases.

$\mathbf{Example~{}~{}6}.$ We take the $\Lambda=\Lambda_{1}\cup\Lambda_{2}$, the degradation coefficient is:

[TABLE]

The LA solution of Example 6 tends to the front right of the region, we can see the Figure 5.11. In this case, we use the asymmetric boundary data to reconstruct the degradation coefficients, here the data is on the $\Lambda_{1}\cup\Lambda_{2}$, and it verifies the LA solution is close to the data location. The solution dependence on data location is more obvious in two dimensions.

$\mathbf{Example~{}~{}7}.$ We take the $\Lambda=\partial\Omega$, the degradation coefficient is:

[TABLE]

It is easy to see $q(x,y)$ is central symmetry, and the reconstruction results are displayed in Figure 5.12. Here, Figure $\ref{two-example 2}(a)$ is the exact solution and the Figure $\ref{two-example 2}(b)$, $\ref{two-example 2}(c)$, $\ref{two-example 2}(d)$ are correspond to average flux data measured on the different part of boundary. The skewness is calculated from the marginal density function. The data measured on boundary $\Lambda_{1}\cup\Lambda_{2}$, the graph is skewed to the front right of the region, that is, the skewness is negative in the x direction and positive in the y direction. When the average flux data measured on $\Lambda_{1}\cup\Lambda_{2}\cup\Lambda_{4}$, the graph slant toward the front of the area i.e., the skewness is zero in the x direction and positive in the y direction. We use the average flux data on the whole boundary $\partial\Omega$, the graph is accordance with the exact solution. The skewness correspond to the average flux data collected on different position of boundary can be seen in Table 7. Notice that here the corresponding marginal distributions of 2-D random variable are used to calculate the skewness by applying formula (4.58). In addition, the 95% confidence intervals of LA solutions at different $(x,y)\in\Omega$ are shown in the Table 8. It can be seen that, owing to the appropriate regularization parameter selection and average flux measurement on whole boundary, the size of confidence region is very small, i.e., the proposed LA algorithm gives a higher reliability.

6 Conclusions

In this paper, we study the identification of degeneracy coefficients in time-fractional diffusion equations (TFDE) by using the average flux data at the accessible part of boundary. We mainly prove that the average flux measurement data can uniquely determine the degradation coefficient. The Lipchitz continuity of the corresponding forward operator is obtained. Due to the average flux measurement data only provide very limited information, and lead to serious ill-posedness of IDCP. This paper combines Tikhonov regularization with Laplace method to overcome the ill-posedness. The existence, stability and convergence of solutions of variational problems are given. The paper introduce the sensitivity problem and the adjoint problem to find the minimizer of the variational problem by using the conjugate gradient method, and derive the mean and variance of the approximate posterior distribution by applying Bayesian theory and Laplace approximation. The Hellinger distance between the exact posterior measure and Laplace approximation is analyzed, and the second order convergence rate at MAP point is proved. The symmetry of the LA solution described by skewness is proposed, and find that the symmetry of solution is closely related to the symmetry of data. Finally, some numerical examples show that the method is not only accurate and flexible, but also can capture statistical information and quantify the uncertainty of the solution.

Acknowledgments

The work described in this paper was supported by the NSF of China (11301168).