A tutorial on inverse problems for anomalous diffusion processes

We shall use extensively the two-parameter Mittag-Leffler function ${{E}_{\alpha ,\beta }}(z)$ (with $\alpha \gt 0$ and $\beta \in \mathbb{R}$) defined by [53, equation (1.8.17), p 40]

$$\begin{eqnarray}&&{{E}_{\alpha ,\beta }}(z)=\mathop{\sum }\limits_{k=0}^{\infty }\frac{{{z}^{k}}}{\Gamma (k\alpha +\beta )}\ \ z\in \mathbb{C}.\end{eqnarray} \tag{ 2.1 }$$

This function with $\beta =1$ was first introduced by Gösta Mittag-Leffler in 1903 [74] and then generalized by others [1, 36]. It can be verified directly that

$$\begin{eqnarray*}&&{{E}_{1,1}}(z)={{{\rm e}}^{z}},\quad {{E}_{2,1}}(z)={\rm cosh} \sqrt{z},\quad {{E}_{2,2}}(z)=\frac{{\rm sinh} \sqrt{z}}{\sqrt{z}}.\end{eqnarray*}$$

Hence it represents a generalization of the exponential function in that ${{E}_{1,1}}(z)={{{\rm e}}^{z}}$. The Mittag-Leffler function ${{E}_{\alpha ,\beta }}(z)$ is an entire function of z with order ${{\alpha }^{-1}}$ and type 1 [53, p 40]. Further, the function ${{E}_{\alpha ,1}}(-t)$ is completely monotone on the positive real axis ${{\mathbb{R}}^{+}}$ [82], and thus it is positive on ${{\mathbb{R}}^{+}}$; see also [90] for extension to the two-parameter Mittag-Leffler function ${{E}_{\alpha ,\beta }}(z)$. It appears in the solution representation for FDEs: the functions ${{E}_{\alpha ,1}}(-\lambda {{t}^{\alpha }})$ and ${{t}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-\lambda {{t}^{\alpha }})$ appear in the kernel of the time fractional diffusion problem with initial data and the right-hand side, respectively, and $x{{E}_{\alpha ,2}}(-\lambda {{x}^{\alpha }})$ and ${{x}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-\lambda {{x}^{\alpha }})$ are eigenfunctions to the fractional Sturm–Liouville problem with a zero potential, cf section 4.1.

In our discussions, the asymptotic behavior of the function ${{E}_{\alpha ,\beta }}(z)$ will play a crucial role. It satisfies the following exponential asymptotics [53, p 43, equations (2.8.17) and (2.8.18)], which was first derived by Djrbashian [21], and refined by many researchers [80, 105].

Lemma 2.1. Let $\alpha \in (0,2)$, $\beta \in \mathbb{R}$, and $\mu \in (\alpha \pi /2,{\rm min} (\pi ,\alpha \pi ))$. Then for $N\in \mathbb{N}$

$$\begin{eqnarray*}\begin{array}{rcl} {{E}_{\alpha ,\beta }}(z) & = & \frac{1}{\alpha }{{z}^{(1-\beta )/\alpha }}{{{\rm e}}^{{{z}^{1/\alpha }}}}-\mathop{\sum }\limits_{k=1}^{N}\frac{1}{\Gamma (\beta -\alpha k)}\frac{1}{{{z}^{k}}}+O(\frac{1}{{{z}^{N+1}}})\ {\rm with}|z|\to \infty ,|{\rm arg} (z)|\leqslant \mu , \\ {{E}_{\alpha ,\beta }}(z) & = & -\mathop{\sum }\limits_{k=1}^{N}\frac{1}{\Gamma (\beta -\alpha k)}\frac{1}{{{z}^{k}}}+O(\frac{1}{{{z}^{N+1}}})\ {\rm with}|z|\to \infty ,\ \mu \leqslant |{\rm arg} (z)|\leqslant \pi . \\ \end{array}\end{eqnarray*}$$

From these asymptotics, the Mittag-Leffler function ${{E}_{\alpha ,\beta }}(z)$ decays only linearly on the negative real axis ${{\mathbb{R}}^{-}}$, which is much slower than the exponential decay for the exponential function e^z. However, on the positive real axis ${{\mathbb{R}}^{+}}$, it grows exponentially, and the growth rate increases with the fractional order $0\lt \alpha \lt 2$. To illustrate the distinct feature, we plot the functions ${{E}_{\alpha ,1}}(-{{\pi }^{2}}{{t}^{\alpha }})$ and ${{t}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-{{\pi }^{2}}{{t}^{\alpha }})$ in figure 1 for several different α values, where $\lambda ={{\pi }^{2}}$ is the first Dirichlet eigenvalue of the negative Laplacian on the unit interval $\Omega =(0,1)$; see appendix A.1 for further details on the computation of the Mittag-Leffler function. Figure 1(a) can be viewed as the time evolution of $u(1/2,t)$, where $\partial _{t}^{\alpha }u-{{u}_{xx}}=0$ with $u(0,t)=u(1,t)=0$, and initial data ${{u}_{0}}(x)={\rm sin} \pi x$ (the lowest Fourier eigenmode). The slow decay behavior at large time is clearly observed. In particular, at t = 1, the function ${{E}_{\alpha ,1}}(-{{\pi }^{2}}t)$ still takes values distinctly away from zero for any $0\lt \alpha \lt 1$, whereas the exponential function ${{{\rm e}}^{-{{\pi }^{2}}t}}$ almost vanishes identically. In contrast, for t close to zero, the picture is reversed: the Mittag-Leffler function ${{E}_{\alpha ,1}}(-{{\pi }^{2}}t)$ decays much faster than the exponential function ${{{\rm e}}^{-{{\pi }^{2}}t}}$. The drastically different behavior of the function ${{E}_{\alpha ,1}}(-z)$, in comparison with the exponential function ${{{\rm e}}^{-z}}$, explains many unusual phenomena with inverse problems for FDEs to be described below. According to the exponential asymptotics, the function ${{E}_{\alpha ,\alpha }}(z)$ decays faster on the negative real axis ${{\mathbb{R}}^{-}}$, since $1/\Gamma (0)=0$, i.e., the first term in the expansion vanishes. This is confirmed numerically in figure 1(b). Even though not shown, it is noted that the function ${{E}_{\alpha ,\alpha }}(z)$ decays only quadratically on the negative real axis ${{\mathbb{R}}^{-}}$ for $\alpha \in (0,1)$ or $\alpha \in (1,2)$, which is asymptotically much slower than the exponential decay.

Zoom In Zoom Out Reset image size

The distribution of zeros of the Mittag-Leffer function ${{E}_{\alpha ,\beta }}(z)$ is of immense interest, especially in the related Sturm–Liouville problem; see section 4.1 below. The case of $\beta =1$ was first studied by Wiman [104]. It was revisited by Djrbashian [21], and many deep results were derived, especially for the case of $\alpha =2$. There are many further refinements [93]; see [83] for an updated account.

The Wright function ${{W}_{\rho ,\mu }}(z)$ is defined by [106, 107]

$$\begin{eqnarray*}&&{{W}_{\rho ,\mu }}(z)=\mathop{\sum }\limits_{k=0}^{\infty }\frac{{{z}^{k}}}{k!\Gamma (\rho k+\mu )},\qquad \mu ,\ \rho \in \mathbb{R},\;\rho \gt -1,\quad z\in \mathbb{C}.\end{eqnarray*}$$

This is an entire function of order $1/(1+\rho )$ [30, theorem 2.4.1]. It was first introduced in connection with a problem in number theory by Edward M Wright, and revived in recent years since it appears as the fundamental solution for FDEs [69]. The Wright function ${{W}_{\rho ,\mu }}(z)$ has the following the asymptotic expansion in one sector containing the negative real axis ${{\mathbb{R}}^{-}}$ [67, theorem 3.2]. Like before, the exponential asymptotics can be used to deduce the distribution of its zeros [66].

Lemma 2.2. Let $-1\lt \rho \lt 0$, $y=-z$, ${\rm arg} (z)\leqslant \pi$, $-\pi \lt {\rm arg} (y)\leqslant \pi$, and for all small $\epsilon \gt 0$, $|{\rm arg} (y)|\leqslant {\rm min} (3\pi (1+\rho )/2,\pi )-\epsilon$. Then

$$\begin{eqnarray*}&&{{W}_{\rho ,\mu }}(z)={{Y}^{1/2-\mu }}{{{\rm e}}^{-Y}}\left\{ \mathop{\sum }\limits_{m=0}^{M-1}{{A}_{m}}{{Y}^{-m}}+O({{Y}^{-M}}) \right\},\ \ Y\to \infty ,\end{eqnarray*}$$

where $Y=(1+\rho ){{({{(-\rho )}^{-\rho }}y)}^{1/(1+\rho )}}$ and the coefficients A_m, $m=0,1,\ldots$ are defined by the asymptotic expansion

$$\begin{eqnarray*}\begin{array}{rcl} \frac{\Gamma (1-\mu -\rho t)}{2\pi {{(-\rho )}^{-\rho t}}{{(1+\rho )}^{(1+\rho )(t+1)}}\Gamma (t+1)} & = & \mathop{\sum }\limits_{m=0}^{M-1}\frac{{{(-1)}^{m}}{{A}_{m}}}{\Gamma ((1+\rho )t+\mu +\frac{1}{2}+m)} \\ {} & {} & +O\left( \frac{1}{\Gamma ((1+\rho )t+\beta +\frac{1}{2}+M)} \right), \\ \end{array}\end{eqnarray*}$$

valid for ${\rm arg} (t),{\rm arg} (-\rho t)$, and ${\rm arg} (1-\mu -\rho t)$ all lying between $-\pi$ and π and t tending to infinity.

The Wright function ${{W}_{\rho ,\mu }}(z)$, $-1\lt \rho \lt 0$, decays exponentially on the negative real axis ${{\mathbb{R}}^{-}}$, in a manner similar to the exponential function e^z, but at a different decay rate. Its special role in fractional calculus is underscored by the fact that it forms the free-space fundamental solution ${{K}_{\alpha }}(x,t)$ to the one-dimensional time fractional diffusion equation [69] by

$$\begin{eqnarray}&&{{K}_{\alpha }}(x,t)=\frac{1}{2{{t}^{\frac{\alpha }{2}}}}{{W}_{-\frac{\alpha }{2},1-\frac{\alpha }{2}}}\left( -|x|/{{t}^{\alpha /2}} \right).\end{eqnarray} \tag{ 2.2 }$$

The multidimensional case is more complex and involves further special functions, in particular, the Fox H function [54, 91]. For $\alpha =1$, the formula (2.2) recovers the familiar free-space fundamental solution for the one-dimensional heat equation, i.e.

$$\begin{eqnarray*}&&K(x,t)=\frac{1}{2\sqrt{\pi t}}{{{\rm e}}^{-\frac{{{x}^{2}}}{4t}}},\end{eqnarray*}$$

which is a Gaussian distribution in x for any $t\gt 0$. In the fractional case, the fundamental solution ${{K}_{\alpha }}(x,t)$ exhibits quite different behavior than the heat kernel. To see this, we show the profile of ${{K}_{\alpha }}(x,t)$ in figure 2 for several α values; see appendix A.1 for a brief description of the computational details. For any $0\lt \alpha \lt 1$, ${{K}_{\alpha }}(x,t)$ decays slower at a polynomial rate as the argument $|x|/{{t}^{\alpha /2}}$ tends to infinity, i.e., having a long tail, when compared with the Gaussian density. The long tail profile is one of distinct features of slow diffusion [5]. Further, for any $\alpha \lt 1$, the profile is only continuous but not differentiable at x = 0. The kink at the origin implies that the solution operator to time fractional diffusion may only have a limited smoothing property.

Zoom In Zoom Out Reset image size

We shall follow the well-established practice in the inverse problem community, i.e., using the singular value decomposition, as the main tool for numerically analyzing the problem behavior [32]. Specifically, we shall numerically compute the forward map F, and analyze its behavior to gain insights into the inverse problem. Given a matrix ${\bf A}\in {{\mathbb{R}}^{n\times m}}$, its singular value decomposition is given by

$$\begin{eqnarray*}&&{\bf A}={\bf U}\Sigma {{{\bf V}}^{{\rm t}}},\end{eqnarray*}$$

where ${\bf U}=[{{{\bf u}}_{1}}\ \ldots \ {{{\bf u}}_{n}}]\in {{\mathbb{R}}^{n\times n}}$ and ${\bf V}=[{{{\bf v}}_{1}}\ \ldots \ {{{\bf v}}_{m}}]\in {{\mathbb{R}}^{m\times m}}$ are column orthonormal matrices, and $\Sigma \in {{\mathbb{R}}^{n\times m}}={\rm diag}({{\sigma }_{1}},\ldots ,{{\sigma }_{r}},0,\ldots ,0)$ is a diagonal matrix, with the diagonal elements $\{{{\sigma }_{i}}\}$ being nonnegative and listed in a descending order ${{\sigma }_{1}}\gt \ldots \gt {{\sigma }_{r}}\gt 0$, and r being the (numerical) rank of the matrix ${\bf A}$. The diagonal element ${{\sigma }_{i}}$ is known as the ith singular value, and the corresponding columns of ${\bf U}$ and ${\bf V}$, i.e., ${{{\bf u}}_{i}}$ and ${{{\bf v}}_{i}}$, are called the left and right singular vectors, respectively.

One simple measure of the conditioning of a linear inverse problem ${\bf Ax}={\bf b}$ is the condition number ${\rm cond}({\bf A})$, which is defined as the ratio of the largest to the smallest nonzero singular value, i.e.

$$\begin{eqnarray*}&&{\rm cond}({\bf A})={{\sigma }_{1}}/{{\sigma }_{r}}.\end{eqnarray*}$$

In particular, if the condition number is small, then the data error will not be amplified much. In the case of a large condition number, the issue is more delicate: it may or may not amplify the data perturbation greatly. A more complete picture is provided by the singular value spectrum $({{\sigma }_{1}},{{\sigma }_{2}},\ldots ,{{\sigma }_{r}})$. Especially, a singular value spectrum gradually decaying to zero without a clear gap is characteristic of many discrete ill-posed problems, which is reminiscent of the spectral behavior of compact operators. We shall adopt these simple tools to analyze related inverse problems below.

In addition, using singular value decomposition and regularization techniques, e.g. Tikhonov regularization or truncated singular value decomposition, one can conveniently obtain numerical reconstructions, even though this might not be the most efficient way to do so. However, we shall not delve into the extremely important question of practical reconstructions, since it relies heavily on a priori knowledge on the sought for solution and the statistical nature (Gaussian, Poisson, Laplace ...) of the contaminating noise in the data, which will depend very much on the specific application. We refer interested readers to the monographs [26, 41, 92] and the survey [49] for updated accounts on regularization methods for constructing stable reconstructing procedures and efficient computational techniques. We will also briefly mention below related works on the application of regularization techniques to inverse problems for FDEs.

First we consider the time fractional backward diffusion. By the linearity of the inverse problem, we may assume that equation (3.1) is prescribed with a homogeneous Dirichlet boundary condition, i.e., u = 0 at x = 0, 1, and the initial condition $u(0)=v$. Then the inverse problem reads: given the final time data $g=u(T)$, find the initial condition v. It arises in, for example, the determination of a stationary contaminant source in underground fluid flow.

To gain insight, we apply the separation of variables. Let $\{({{\lambda }_{j}},{{\phi }_{j}})\}$, with ${{\lambda }_{j}}={{(j\pi )}^{2}}$ and ${{\phi }_{j}}=\sqrt{2}{\rm sin} j\pi x$, be the Dirichlet eigenpairs of the negative Laplacian on the interval Ω. The eigenfunctions $\{{{\phi }_{j}}\}$ form an orthonormal basis of the ${{L}^{2}}(\Omega )$ space. Then using the Mittag-Leffler function ${{E}_{\alpha ,\beta }}(z)$ defined in (2.1), the solution u to equation (3.1) can be expressed as

$$\begin{eqnarray*}&&u(x,t)=\mathop{\sum }\limits_{j=1}^{\infty }{{E}_{\alpha ,1}}(-{{\lambda }_{j}}{{t}^{\alpha }})(v,{{\phi }_{j}}){{\phi }_{j}}(x).\end{eqnarray*}$$

Therefore, the final time data $g=u(T)$ is given by

$$\begin{eqnarray*}&&g(x)=\mathop{\sum }\limits_{j=1}^{\infty }{{E}_{\alpha ,1}}(-{{\lambda }_{j}}{{T}^{\alpha }})(v,{{\phi }_{j}}){{\phi }_{j}}(x).\end{eqnarray*}$$

It follows directly that the initial data v is formally given by

$$\begin{eqnarray*}&&v=\mathop{\sum }\limits_{j=1}^{\infty }\frac{(g,{{\phi }_{j}})}{{{E}_{\alpha ,1}}(-{{\lambda }_{j}}{{T}^{\alpha }})}{{\phi }_{j}}.\end{eqnarray*}$$

Since the function ${{E}_{\alpha ,1}}(-t)$ is completely monotone on the positive real axis ${{\mathbb{R}}^{+}}$ [82] for any $\alpha \in (0,1]$, the denominator in the representation does not vanish. In case of $\alpha =1$, the formula reduces to the familiar expression

$$\begin{eqnarray*}&&v=\mathop{\sum }\limits_{j=1}^{\infty }{{{\rm e}}^{{{\lambda }_{j}}T}}(g,{{\phi }_{j}}){{\phi }_{j}}.\end{eqnarray*}$$

This formula shows clearly the well-known, severely ill-posed nature of the backward diffusion problem: the perturbation in the jth Fourier mode $(g,{{\phi }_{j}})$ of the (noisy) data g is amplified by an exponentially growing factor ${{{\rm e}}^{{{\lambda }_{j}}T}}$, which can be astronomically large, even for a very small index j, if the terminal time T is not very small. Hence it is always severely ill-conditioned and we must multiply the jth Fourier mode of the data g by a factor ${{{\rm e}}^{{{\lambda }_{j}}T}}$ in order to recover the corresponding mode of the initial data v.

In the fractional case, by lemma 2.1, the Mittag-Leffler function ${{E}_{\alpha ,1}}(z)$ decays only linearly on the negative real axis ${{\mathbb{R}}^{-}}$, and thus the multiplier $1/{{E}_{\alpha ,1}}(-{{\lambda }_{j}}{{T}^{\alpha }})$ grows only linearly in ${{\lambda }_{j}}$, i.e., $1/{{E}_{\alpha ,1}}(-{{\lambda }_{j}}{{T}^{\alpha }})\sim {{\lambda }_{j}}$, which is very mild compared to the exponential growth ${{{\rm e}}^{{{\lambda }_{j}}T}}$ for the case $\alpha =1$, and thus the fractional case is only mildly ill-posed. Roughly, the jth Fourier mode of the initial data v now equals the jth mode of the data g multiplied by ${{\lambda }_{j}}$. More precisely, it amounts to the loss of two spatial derivatives [88, theorem 4.1]

$$\begin{eqnarray}&&\parallel v{{\parallel }_{{{L}^{2}}(\Omega )}}\leqslant c\parallel u(T){{\parallel }_{{{H}^{2}}(\Omega )}}.\end{eqnarray} \tag{ 3.2 }$$

Intuitively, the history mechanism of the anomalous diffusion process retains the complete dynamics of the physical process, including the initial data, and thus it is much easier to go backwards to the initial state v. This is in sharp contrast to classical diffusion, which has only a short memory and loses track of the preceding states quickly. This result has become quite well-known in the inverse problems community and has contributed to a belief that 'inverse problems for FDEs are less ill-conditioned than their classical counterparts'—throughout this paper we will see that this conclusion as a general statement can be quite far from the truth.

Does this mean that for all terminal time T the fractional case is always less ill-posed than the classical one? The answer is yes, in the sense of the norm on the data space in which the data g lies. Does this mean that from a computational stability standpoint that one can always solve the backward fractional problem more effectively than for the classical case? The answer is no, and the difference can be substantial. To illustrate the point, let J be the highest frequency mode required of the initial data v and assume that we believe we are able to multiply the first J modes ${{g}_{j}}=(g,{{\phi }_{j}})$, $j=1,2,\ldots ,J$, by a factor no larger than M (which roughly assumes that the noise levels in both cases are comparable). By the monotonicity of the function ${{E}_{\alpha ,1}}(-t)$ in t, it suffices to examine the Jth mode. For the heat equation ${{v}_{J}}:=(v,{{\phi }_{J}})={{{\rm e}}^{{{\lambda }_{J}}T}}{{g}_{J}}$ and provided that $T={{T}_{J}}\lt {{\lambda }_{J}}/{\rm log} M$ this is feasible. For a fixed J, if $T_{\alpha }^{\star }$ denotes the point where

$$\begin{eqnarray*}&&{{{\rm e}}^{-{{\lambda }_{J}}T_{\alpha }^{\star }}}={{E}_{\alpha ,1}}(-{{\lambda }_{J}}T_{\alpha }^{\star }),\end{eqnarray*}$$

then in the fractional case for $T\lt T_{\alpha }^{\star }$ the growth factor on g_J will exceed M for any $T\lt T_{\alpha }^{\star }$. In table 1, we present the critical value $T_{\alpha }^{\star }$ for several values of the fractional order α and the maximum number of modes J. The numbers in the table are very telling. For example, for the case J = 5, $\alpha =1/4$ and T = 0.02 (which is approximately one half the value of $T_{\alpha }^{\star }$), the growth factor is about 1.6 for the heat equation but about 113 for the fractional case. With J = 10 and $\alpha =1/4$ and $T=T_{\alpha }^{\star }$ the growth factor is around 336. If $T=T_{\alpha }^{\star }/10$ then it has again dropped to less than 2 for the heat equation but about 190 for the fractional case. Of course, for $T\gt T_{\alpha }^{\star }$ the situation completely reverses. With J = 10, $\alpha =1/4$ and $T=10\;T_{\alpha }^{\star }$ the growth factor is a possibly workable value of around 600; while for the heat equation it is greater than 10²⁵. We reiterate that the apparent contradiction between the theoretical ill-conditioning and numerical stability is due to the spectral cutoff present in any practical reconstruction procedure.

Table 1.  The critical values $T_{\alpha }^{*}$ for fractional backward diffusion.

$\alpha \backslash J$	3	5	10
1/4	0.0442	0.0197	0.0059
1/2	0.0387	0.0163	0.0049
3/4	0.0351	0.0142	0.0040

Next we examine the influence of the fractional order α on the inversion step more closely. To this end, we expand the initial condition v in the piecewise linear finite element basis functions defined on a uniform partition of the domain $\Omega =(0,1)$ with 100 grid points. Then we compute the discrete forward map F from the initial condition to the final time data $g=u(T)$, defined on the same mesh. Numerically, this can be achieved by a fully discrete scheme based on the L1 approximation in time and the finite difference method in space; see appendix A.2 for a description of the numerical method. The ill-posed behavior of the discrete inverse problem is then analyzed using singular value decomposition. A similar experimental setup will be adopted for other examples below.

The numerical results are shown in figure 3. The condition number of the (discrete) forward map F stays mostly around $O({{10}^{4}})$ for a fairly broad range of α values, which holds for all three different terminal times T. This can be attributed to the fact that for any $\alpha \in (0,1),$ backward fractional diffusion amounts to a two spacial derivative loss, cf (3.2). Unsurprisingly, as the fractional order α approaches unity, the condition number eventually blows up, recovering the severely ill-posed nature of the classical backward diffusion problem, cf figure 3(a). Further, we observe that the smaller is the terminal time T, the quicker is the blowup. The precise mechanism for this observation remains unclear. Interestingly, the condition number is not monotone with respect to the fractional order α, for a fixed T. This might imply potential nonuniqueness in the simultaneous recovery of the fractional order α and the initial data v. The singular value spectra at T = 0.01 are shown in figure 3(b). Even though the condition numbers for $\alpha =1/4$ and $\alpha =1/2$ are quite close, their singular value spectra actually differ by a multiplicative constant, but their decay rates are almost identical, thereby showing comparable condition numbers. This shift in singular value spectra can be explained by the local decay behavior of the Mittag-Leffler function, cf figure 1(a): the smaller is the fractional order α, the faster is the decay around t = 0.

Zoom In Zoom Out Reset image size

Even though the condition number is very informative about the (discretized) problem, it does not provide a full picture, especially when the condition number is large. In this case the singular value spectrum can be far more revealing. The spectra for two different α values are given in figure 4. At $\alpha =1/2$, the singular values decay at almost the same algebraic rate, irrespective of the terminal time T. This is expected from the two-derivative loss for any $\alpha \lt 1$. However, for $\alpha =1$, the singular values decay exponentially, and the decay rate increases dramatically with the increase of the terminal time T. For T = 0.001, there are a handful of 'significant' singular values, say above 10⁻³, but when the time T increases to T = 1, there is only one meaningful singular value remaining. The distribution of the singular values has important practical consequences. For a small time T, the first few singular values for the classical diffusion case actually might be much larger than that for the fractional case, which indicates that the classical case is actually numerically much easier to recover in this regime, concurring with the observations drawn from table 1. For example, at T = 0.001, the first twenty singular values are larger than the fractional counterpart, cf figure 3(b), and hence, the first twenty modes, i.e., left singular vectors, are more stable in the reconstruction procedure.

Zoom In Zoom Out Reset image size

The mathematical model (3.1) is rescaled with a unit diffusion coefficient. In practice, there is always a diffusion coefficient σ in the elliptic operator, i.e.

$$\begin{eqnarray*}&&\partial _{t}^{\alpha }u-\nabla \cdot (\sigma \nabla u)+qu=f.\end{eqnarray*}$$

For example, the thermal conductivity σ of the gun steel at moderate temperature is about $1.8\times {{10}^{-5}}\;{{{\rm m}}^{2}}\;{{{\rm s}}^{-1}}$ [11] and the diffusion coefficient of the oxygen in water at 25 °C is $2.10\times {{10}^{-5}}\;{\rm c}{{{\rm m}}^{2}}\;{{{\rm s}}^{-1}}$ [18]. Mathematically, this does not change the ill-posed nature of the inverse problem. However, the presence of a diffusion coefficient σ has important consequence: it enables the practical feasibility of the classical backward diffusion problem (and likely for many other inverse problems for the diffusion equation). Physically, a constant conductivity σ amounts to rescaling the final time T by $T^{\prime} =T/\sigma$. In the fractional case, a similar but nonlinear scaling law $T^{\prime} =T/{{\sigma }^{1/\alpha }}$ remains valid.

Numerically, time fractional backward diffusion has been extensively studied. Liu and Yamamoto [65] proposed a numerical scheme for the one-dimensional fractional backward problem based on the quasi-reversibility method [57], and derived error estimates for the approximation, under a priori smoothness assumption on the initial condition. This represents one of the first works on inverse problems in anomalous diffusion. Later, Wang and Liu [99] studied total variation regularization for two-dimensional fractional backward diffusion, and analyzed its well-posedness of the optimization problem and the convergence of an iterative scheme of Bregman type. Wei and Wang [102] developed a modified quasi-boundary value method for the problem in a general domain, and established error estimates for both a priori and a posteriori parameter choice rules. In view of better stability results in the fractional case, one naturally expects better error estimates than the classical diffusion equation, which is confirmed by these studies.

Next we consider the sideways problem for time fractional diffusion. There are several possible formulations, e.g., the quarter plane and the finite space domain. The quarter plane sideways fractional diffusion problem is as follows. Let $u(x,t)$ be defined in $(0,\infty )\times (0,\infty )$ by

$$\begin{eqnarray*}&&\partial _{t}^{\alpha }u-{{u}_{xx}}=0,\qquad x\gt 0,\ t\gt 0,\end{eqnarray*}$$

and the boundary and initial conditions

$$\begin{eqnarray*}&&u(x,0)=0\quad {\rm and}\quad u(0,t)=f(t),\end{eqnarray*}$$

where we assume $|u(x,t)|\leqslant {{c}_{1}}{{{\rm e}}^{{{c}_{2}}{{x}^{2}}}}$. We do not know the left boundary condition f, but are able to measure u at an intermediate point $x=L\gt 0$, $h(t)=u(L,t)$. The inverse problem is: given the (noisy) data h, find the boundary condition f. The solution u of the forward problem is given by a convolution integral with the kernel being the spatial derivative ${{K}_{\alpha ,x}}(x,s)$ of the fundamental solution ${{K}_{\alpha }}(x,s)$, cf (2.2), by

$$\begin{eqnarray*}&&u(x,t)=\int _{0}^{t}{{K}_{\alpha ,x}}(x,t-s)f(s){\rm d}s.\end{eqnarray*}$$

This representation is well known for the case $\alpha =1$, and it was first derived by Carasso [11]; see also [8] for related discussions. It leads to a convolution integral equation for the unknown f in terms of the given data h

$$\begin{eqnarray*}&&h(t)=\int _{0}^{t}{{R}_{\alpha }}(t-s)f(s){\rm d}s,\end{eqnarray*}$$

where the convolution kernel ${{R}_{\alpha }}(s)$ is given by a Wright function in the form

$$\begin{eqnarray*}\begin{array}{rcl} {{R}_{\alpha }}(s) & = & \frac{1}{2{{s}^{\alpha }}}{{W}_{-\frac{\alpha }{2},2-\frac{\alpha }{2}}}(-L{{s}^{-\alpha /2}})=\mathop{\sum }\limits_{k=0}^{\infty }\frac{{{(-L)}^{k}}}{k!\;\Gamma (-\frac{\alpha }{2}k+2-\frac{\alpha }{2})}\;{{s}^{-k\frac{\alpha }{2}-\alpha }}. \\ \end{array}\end{eqnarray*}$$

In case of $\alpha =1$, i.e., classical diffusion, the kernel R(s) is given explicitly by

$$\begin{eqnarray*}&&R(s)=\frac{L}{2\sqrt{\pi }}{{s}^{-\frac{3}{2}}}{{{\rm e}}^{-\frac{{{L}^{2}}}{4s}}}\in {{C}^{\infty }}(0,\infty ).\end{eqnarray*}$$

Since all its derivatives vanish at s = 0, the classical theory of Volterra integral equations of the first kind [56] implies the extreme ill-conditioning of the problem. This is not surprising: we are, after all, mapping a function $f\in {{C}^{0}}(0,\infty )$ to an element in ${{C}^{\infty }}(0,\infty )$. The conditioning of the time fractional sideways problem again depends on the convolution kernel ${{R}_{\alpha }}$ and its derivatives at s = 0 and in this case is the value of the Wright function ${{W}_{-\alpha /2,2-\alpha /2}}(-z)$ and its derivatives as $z\to \infty$. These are again zero (in fact the Wright function ${{W}_{-\alpha /2,2-\alpha /2}}(z)$ also decays exponentially to zero for large negative arguments, cf lemma 2.2), and thus the fractional sideways problem is also severely ill-posed. However, this analysis does not show their difference in the degree of ill-posedness: even though both are severely ill-posed, their practical computational behavior can still be quite different, as we shall see below.

To see their difference in the degree of ill-posedness, we examine another variant of the sideways problem on a finite interval $\Omega =(0,1)$, with Cauchy data prescribed on the axis x = 0, i.e. given zero initial condition ${{u}_{0}}=u(x,0)=0$, recovering $h=u(1,t)$ from the lateral Cauchy data at x = 0:

$$\begin{eqnarray*}&&u(0,t)=f(t),\;\;{{u}_{x}}(0,t)=g(t),\;\;t\geqslant 0.\end{eqnarray*}$$

This problem is also known as the lateral Cauchy problem in the literature. In the case $\alpha =1$, it is known that the inverse problem is severely ill-posed [8, 33]. To gain insight into the fractional case, we apply the Laplace transform. With $\mathop{\;}\limits^{\widehat{{}}}$ being the Laplace transform in time, and noting the Laplace transform of the Caputo derivative $\widehat{\partial _{t}^{\alpha }u}={{z}^{\alpha }}\hat{u}(z)-{{z}^{\alpha -1}}{{u}_{0}}$ [53, lemma 2.24], we deduce

$$\begin{eqnarray*}&&{{z}^{\alpha }}\hat{u}(x,z)-{{\hat{u}}_{xx}}(x,z)=0,\;\;\hat{u}(0)=\hat{f},\;\;{{\hat{u}}_{x}}(0)=\hat{g}.\end{eqnarray*}$$

The general solution $\hat{u}$ is given by $\hat{u}(x,z)=\hat{f}{\rm cosh} {{z}^{\alpha /2}}x+\hat{g}\frac{{\rm sinh} {{z}^{\alpha /2}}x}{{{z}^{\alpha /2}}x}$ and thus the solution $\hat{h}(z)=\hat{u}(1,z)$ at x = 1 is given by

$$\begin{eqnarray*}&&\hat{h}=\hat{f}{\rm cosh} {{z}^{\alpha /2}}+\hat{g}\frac{{\rm sinh} {{z}^{\alpha /2}}}{{{z}^{\alpha /2}}}.\end{eqnarray*}$$

The solution h(t) can then be recovered by an inverse Laplace transform

$$\begin{eqnarray}&&h(t)=\frac{1}{2\pi {\rm i}}{{\int }_{{\rm Br}}}\hat{h}{{{\rm e}}^{zt}}{\rm d}z,\end{eqnarray} \tag{ 3.3 }$$

where ${\rm Br}=\{z\in \mathbb{C}:\Re z=\sigma ,\sigma \gt 0\}$ is the Bromwich path. Upon deforming the contour suitably, this formula will allow the development of an efficient numerical scheme for the sideways problem via quadrature rules [103], provided that the lateral Cauchy data is available for all $t\gt 0$. The expression (3.3) indicates that, in the fractional case, the sideways problem still suffers from severe ill-posedness in theory, since the high frequency modes of the data perturbation are amplified by an exponentially growing multiplier ${{{\rm e}}^{{{z}^{\alpha /2}}}}$. However, numerically, the degree of ill-posedness decreases dramatically as the fractional order α decreases from unity to zero, since as α → 0⁺, the multipliers are growing at a much slower rate, and thus we have a better chance of recovering many more modes of the boundary data. In other words, both the classical and fractional sideways problems are severely ill-posed in the sense of error estimates between the norms in the data and unknowns; but with a fixed frequency range, the behavior of the time fractional sideways problem can be much less ill-posed. Hence, anomalous diffusion mechanism does help substantially since much more effective reconstructions are possible in the fractional case.

Next we illustrate the point numerically. The numerical results for the sideways problem are given in figure 5. It is observed that the degree of ill-posedness of the finite-dimensional discretized version of the inverse problem indeed decreases dramatically with the decrease of the fractional order α, cf figure 5(a), which agrees well with the preceding discussions. Surprisingly, for T = 1 there is a sudden transition around $\alpha =1/2$, below which the sideways problem behaves as if nearly well-posed, but above which the conditioning deteriorates dramatically with the increase of the fractional order α and eventually it recovers the properties of the classical sideways problem. Similar transitions are observed for other terminal times. This might be related to the discrete setting, for which there is an inherent frequency cutoff. Further, as the fractional order α approaches zero, the problem reaches a quasi-steady state much quicker and thus the forward map F can have only fairly localized elements along the main diagonal. To give a more complete picture, we examine the singular value spectrum in figure 5(b). Unlike the backward diffusion problem discussed earlier, the singular values are actually decaying only algebraically, even for $\alpha =1$, and then there might be a few tiny singular values contributing to the large condition number. The larger is the fractional order α, the more tiny singular values are in the spectrum. Hence, in the discrete setting, even for $\alpha =3/4$, the problem is still nearly well-posed, despite the large apparent condition number, since a few tiny singular values with a distinct gap from the rest of the spectrum are harmless in most regularization techniques.

Zoom In Zoom Out Reset image size

Physically this can also be observed in figure 6, where the forward map F is from the Dirichlet boundary condition x = 1 to the flux boundary condition at x = 0, in a piecewise linear finite element basis. Pictorially, the forward map F is only located in the upper left corner and has a triangular structure, which reflects the casual or Volterra nature of the sideways problem for the fractional diffusion equation. We note that the causal structure should be utilized in developing reconstruction techniques, via, e.g., Lavrentiev regularization [56]. For small α values, e.g., $\alpha =1/4$, the finite element basis at the right end point x = 1 is almost instantly transported to the left end point x = 0, whose magnitude is slightly decreased, but with little diffusive effect, resulting a diagonally dominant forward map. However, as the fractional order α increases towards unity, the diffusive effect eventually kicks in, and the information spreads over the whole interval. Further, for large α values, it takes much longer time to reach the other side and there is a lag of information arrival, which explains the presence of tiny singular values. The larger is the fractional order α, the smaller is the magnitude, i.e., the less is the amount of the information reached the other side. Hence, one feasible approach is to recover only the boundary condition over a smaller subinterval of the measurement time interval. This idea underlies one popular engineering approach, the sequential function specification method [4, 64].

Zoom In Zoom Out Reset image size

The sideways problem for the classical diffusion has been extensively studied, and many efficient numerical methods have been developed and analyzed [8, 11, 24, 25]. In the fractional case, however, there are only a few works on numerical schemes, mostly for one-dimensional problems, and there seems no theoretical study on stability etc. Murio [76, 77] developed several numerical schemes, e.g., based on space marching and finite difference method, for the sideways problem, but without any analysis. Qian [85] discussed about the ill-posedness of the quarter plane formulation of the sideways problem using the Fourier analysis, based on which a mollifier method was proposed, with error estimates provided. In [87], the recovery of a nonlinear boundary condition from the lateral Cauchy data was studied using an integral equation approach, and a convergent fixed point iteration method was suggested. The influence of the imprecise specification of the fractional order α on the reconstruction was examined. Zheng and Wei [113] proposed a mollification method for the quarter plane formulation of the sideways problem, by convoluting the fractional derivative with a smooth kernel, and derived error estimates for the approximation, under a prior bounds on the solution. The Cauchy problem of the time fractional diffusion has been numerically studied in [114]. In particular, with the separation of variables, a Volterra integral equation reformulation of the problem was derived, from which the ill-posedness of the Cauchy problem follows directly. All these works are concerned with the one-dimensional case, and the high dimensional case has not been studied.

A third classical linear inverse problem for the diffusion equation is the inverse source problem, i.e., the recovery of the source term f from lateral boundary data or final time data. Clearly, one piece of boundary data or final time data alone is insufficient to uniquely determine a general source term, due to dimensional disparity. To restore the possible uniqueness, as usual, we look for only a space- or time-dependent component of the source term f. With different combinations of the data and source term, we get several different (and not equivalent) formulations of the inverse source problems. Below we examine several of them briefly. By the linearity of the forward problem, we without loss of generality, assume a zero initial data v = 0 and a zero potential q = 0 throughout this part.

First, suppose we can measure the solution u at the final time t = T, and aim at recovering either a space dependent or time dependent component of the source term f. Like before, we resort to the separation of variables. For the case of a space dependent only source term f(x), the solution u to the forward problem is given by

$$\begin{eqnarray*}\begin{array}{rcl} u(t) & = & \mathop{\sum }\limits_{j=1}^{\infty }\int _{0}^{t}{{(t-\tau )}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-{{\lambda }_{j}}{{(t-\tau )}^{\alpha }})(f,{{\phi }_{j}}){{\phi }_{j}}{\rm d}\tau \\ {} & = & \mathop{\sum }\limits_{j=1}^{\infty }\frac{1}{{{\lambda }_{j}}}(1-{{E}_{\alpha ,1}}(-{{\lambda }_{j}}{{t}^{\alpha }}))(f,{{\phi }_{j}}){{\phi }_{j}}. \\ \end{array}\end{eqnarray*}$$

Hence the measured data $g=u(T)$ is given by

$$\begin{eqnarray*}&&g=\mathop{\sum }\limits_{j=1}^{\infty }\frac{1}{{{\lambda }_{j}}}(1-{{E}_{\alpha ,1}}(-{{\lambda }_{j}}{{T}^{\alpha }}))(f,{{\phi }_{j}}){{\phi }_{j}}.\end{eqnarray*}$$

By taking inner product with ${{\phi }_{j}}$ on both sides, we arrive at the following representation of the source term f in terms of the measured data g

$$\begin{eqnarray}&&f=\mathop{\sum }\limits_{j=1}^{\infty }{{\lambda }_{j}}\frac{(g,{{\phi }_{j}})}{1-{{E}_{\alpha ,1}}(-{{\lambda }_{j}}{{T}^{\alpha }})}{{\phi }_{j}}.\end{eqnarray} \tag{ 3.4 }$$

By the complete monotonicity of the Mittag-Leffler function ${{E}_{\alpha ,1}}(-t)$ on the positive real axis ${{\mathbb{R}}^{+}}$ [82], we deduce

$$\begin{eqnarray*}&&1\gt {{E}_{\alpha ,1}}(-{{\lambda }_{1}}{{T}^{\alpha }})\gt {{E}_{\alpha ,1}}(-{{\lambda }_{2}}{{T}^{\alpha }}),\end{eqnarray*}$$

and thus the formula (3.4) is well defined for any $T\gt 0$, and gives the precise condition for the existence of a source term. Even with a modest value of the terminal time T, the factor $1-{{E}_{\alpha ,1}}(-{{\lambda }_{1}}{{T}^{\alpha }})$ is close to unity for all small α values, especially for those close to zero. Each frequency component $(f,{{\phi }_{j}})$ differs from $(g,{{\phi }_{j}})$ essentially by a factor ${{\lambda }_{j}}$, which amounts to two derivative loss in space. Actually one can show

$$\begin{eqnarray*}&&\parallel f{{\parallel }_{{{L}^{2}}(\Omega )}}\leqslant c\parallel g{{\parallel }_{{{H}^{2}}(\Omega )}}.\end{eqnarray*}$$

This behavior is identical with that for the backward fractional diffusion. The statement holds also for the inverse source problem for the classical diffusion case. This is not surprising, since with a space dependent source term f, the solution u to the forward problem can be split into the steady solution u_s and the decaying solution u_d, i.e., $u={{u}_{s}}+{{u}_{d}}$, where u_s and u_d solve

$$\begin{eqnarray*}&&-u_{s}^{\prime\prime }=f,\;\;{{u}_{s}}(0)={{u}_{s}}(1)=0,\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&\partial _{t}^{\alpha }{{u}_{d}}-{{u}_{d,xx}}=0,\;\;{{u}_{d}}(0,x)=f(x),\;\;{{u}_{d}}(0,t)={{u}_{d}}(1,t)=0,\end{eqnarray*}$$

respectively. By the decay behavior of the solution u_d, the steady state component u_s is dominating, which amounts to a two spatial derivative loss. This is fully confirmed by the numerical experiments, cf figure 7. It is observed that the condition number is almost independent of the fractional order α, and it is of order $O({{10}^{3}})$, reflecting the mildly ill-posed nature of the inverse problem. In particular, for large terminal time T, the singular value spectra are almost identical for all fractional orders, decaying to zero at an algebraic rate, cf figure 7(b).

Zoom In Zoom Out Reset image size

Next we turn to the time dependent case, i.e., seeking a source term f of the form $f(x,t)=p(t)q(x)$, with a known spacial component q(x), from the final time data $g=u(T)$. Mathematically, the inverse problem even for the classical diffusion equation has not been completely analyzed. The inclusion of a nontrivial term q(x) is important since without this there is nonuniqueness. To see this, we take u to satisfy ${{u}_{t}}-{{u}_{xx}}=f(t)$ on $(0,1)\times (0,T)$ with initial data $u(x,0)=1$ and a homogeneous Neumann boundary condition $-{{u}_{x}}(0,t)={{u}_{x}}(1,t)=0$. Then one solution satisfying $u(x,T)=g(x)=1$ is given by $u(x,t)=1$ and $f\equiv 0$, but another is $u(x,t)={\rm cos} (2\pi t/T)$ and $f=(-2\pi /T){\rm sin} (2\pi t/T)$. Likewise, in the fractional case, we can take $u={\rm cos} (2\pi t/T)$ for the second solution and define f to be its αth order Djrbashian–Caputo fractional derivative in time.

Like previously, the solution u to (3.1) is given by

$$\begin{eqnarray*}&&u(t)=\mathop{\sum }\limits_{j=1}^{\infty }\int _{0}^{t}{{(t-\tau )}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-{{\lambda }_{j}}{{(t-\tau )}^{\alpha }})p(\tau ){\rm d}\tau (q,{{\phi }_{j}}){{\phi }_{j}}.\end{eqnarray*}$$

Hence the measured data $g(x)=u(x,T)$ is given by

$$\begin{eqnarray*}&&g(x)=\mathop{\sum }\limits_{j=1}^{\infty }\int _{0}^{T}{{(T-\tau )}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-{{\lambda }_{j}}{{(T-\tau )}^{\alpha }})p(\tau ){\rm d}\tau (q,{{\phi }_{j}}){{\phi }_{j}}(x).\end{eqnarray*}$$

By taking inner product with ${{\phi }_{j}}$ on both sides, we deduce

$$\begin{eqnarray*}&&(g,{{\phi }_{j}})=(q,{{\phi }_{j}})\int _{0}^{T}{{(T-\tau )}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-{{\lambda }_{j}}{{(T-\tau )}^{\alpha }})p(\tau ){\rm d}\tau .\end{eqnarray*}$$

In the case of $\alpha =1$, the formula recovers the relation

$$\begin{eqnarray*}&&(g,{{\phi }_{j}})=(q,{{\phi }_{j}})\int _{0}^{T}{{{\rm e}}^{-{{\lambda }_{j}}(T-\tau )}}p(\tau ){\rm d}\tau ,\end{eqnarray*}$$

which resembles a finite-time Laplace transform or moment problem, and thus severely smoothing, which renders the inverse source problem severely ill-posed. Intuitively, the term ${{{\rm e}}^{-{{\lambda }_{j}}(T-t)}}$ can only pick up the information for t close to the terminal time T, and for t away from T, the information is severely damped, especially for high frequency modes, which leads to the severely ill-posed nature of the inverse problem. In the fractional case, the forward map F from the unknown to the data is clearly compact, and thus the problem is still ill-posed. However, the kernel ${{t}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-{{\lambda }_{j}}{{t}^{\alpha }})$ is less smooth and decays much slower, and one might expect that the problem is less ill-posed than the canonical diffusion counterpart. To examine the point, we present the numerical results for the inverse problem in figure 8. It is severely ill-posed irrespective of the fractional order α: the singular values decay exponentially to zero without a distinct gap in the spectrum. In particular, for the terminal time T = 1, the spectrum is almost identical for all fractional orders α. For small T, the singular values still decay exponentially, but the rate is different: the smaller is the fractional order α, the faster is the decay, cf figure 8(a). Consequently, a few more modes of the source term $p(\tau )$ might be recovered. In other words, due to a slower local decay of the exponential function ${{{\rm e}}^{-\lambda t}}$, compared with the Mittag-Leffler function ${{t}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-\lambda {{t}^{\alpha }})$, cf figure 1(a), actually more frequency modes can be picked up by normal diffusion than the fractional counterpart, cf figure 8(a). This indicates that with sufficiently accurate data, at a small time instance, the sideways problem for normal diffusion may allow recovering more modes, i.e., anomalous diffusion does not help solve the inverse problem.

Zoom In Zoom Out Reset image size

In practice, the accessible data can also be the flux data at the end point, e.g., x = 0 or x = 1. We briefly discuss the case of recovering a time dependent component p(t) in the source term $f=q(x)p(t)$ from the flux data at x = 0. By repeating the preceding argument, the data $g:=-{{u}_{x}}(0,t)$ is related to the unknown p(t) by

$$\begin{eqnarray*}&&g(t)=-\mathop{\sum }\limits_{j=1}^{\infty }\int _{0}^{t}{{(t-\tau )}^{\alpha -1}}{{E}_{\alpha ,\alpha }}(-{{\lambda }_{j}}{{(t-\tau )}^{\alpha }})p(\tau ){\rm d}\tau (q(x),{{\phi }_{j}})\phi _{j}^{\prime }(0).\end{eqnarray*}$$

In [88, theorem 4.4], a stability result was established for the recovery of the time dependent component p(t). Along the same line of thought, under reasonable assumptions, one can deduce that

$$\begin{eqnarray*}&&\parallel p{{\parallel }_{C[0,T]}}\leqslant c\parallel \partial _{t}^{\alpha }g{{\parallel }_{C[0,T]}}.\end{eqnarray*}$$

The inverse problem roughly amounts to taking the αth order Djrbashian–Caputo fractional derivative in time. Hence as the fractional order α decreases from unity to zero, it becomes less and less ill-posed. For α close to zero, it is nearly well-posed, at least numerically. In other words, anomalous diffusion can mitigate the degree of ill-posedness for the inverse problem. To illustrate the discussion, we present in figure 9 some numerical results, where the forward map F is from the time dependent component p(t) to the flux data g(t) at x = 0, both defined over the interval $[0,T]$, discretized using a continuous piecewise linear finite element basis. The condition number of the discrete forward map F decreases monotonically as the fractional order α decreases from unity to zero, confirming the preceding discussions. Further, the terminal time T does not affect the condition number to a large extent.

Zoom In Zoom Out Reset image size

It is widely accepted in inverse heat conduction that an inverse problem will be severely ill-posed when the data and unknown are not aligned in the same space/time direction, and only mildly ill-posed when they do align with each other. Our discussions with the inverse source problems indicate that the observation remains valid in the time fractional diffusion case. In particular, although not presented, we note that the inverse source problem of recovering a space dependent component from the lateral Cauchy data is severely ill-posed for both fractional and normal diffusion. In the simplest case of a space dependent-only source term, it is mathematically equivalent to unique continuation, a well known example of severely ill-posed inverse problems.

The inverse source problems for the classical diffusion equation have been extensively studied; see e.g., [7, 9, 37]. Inverse source problems for FDEs have also been numerically studied. Zhang and Xu [111] established the unique recovery of a space dependent source term in (3.1) with pure Neumann boundary data and overspecified Dirichlet data at x = 0. This is achieved by an eigenfunction expansion and Laplace transform, and the uniqueness follows from a unique continuation principle of analytic functions. Sakamoto and Yamamoto [89] discussed the inverse problem of determining a spatially varying function of the source term by final overdetermined data in multi-dimensional space, and established its well-posedness in the Hadamard sense except for a discrete set of values of the diffusion constant, using an analytic Fredholm theory. Very recently, Luchko et al [68] showed the uniqueness of recovering a nonlinear source term from the boundary measurement, and developed a numerical scheme of fixed point iteration type. Aleroev et al [2] showed the uniqueness of recovering a space dependent source term from integral type observational data. Recently, there are many numerical studies on this class of inverse problems. In [101], the numerical recovery of a spatially varying function of the source term from the final time data in a general domain was studied using a quasi-boundary value problem method; see also [98, 112] related studies. Wang et al [100] proposed to determine the space-dependent source term from the final time data in multi-dimension using a reproducing kernel Hilbert space method.

Now we consider a nonlinear inverse coefficient problem for the time fractional diffusion equation: given the final time data $g=u(T)$, find the potential q in the model

$$\begin{eqnarray}&&\partial _{t}^{\alpha }u-{{u}_{xx}}+qu=0\quad \;{\rm in}\;\Omega ,\end{eqnarray} \tag{ 3.5 }$$

with a homogeneous Neumann boundary condition and initial data v. The parabolic counterpart has been extensively studied [15, 16, 38], where it was shown that the problem is nearly well-posed in the Hardamard sense in suitable Hölder space, under certain conditions, using the strong maximum principle. In [38], an elegant fixed point method was developed, and the monotone convergence of the method was established. It can be adapted straightforwardly to the fractional case: given an initial guess q⁰, compute the update q^k recursively by

$$\begin{eqnarray*}&&{{q}^{k+1}}=\frac{g^{\prime\prime} -\partial _{t}^{\alpha }u(x,T;{{q}^{k}})}{g},\end{eqnarray*}$$

where the notation $u(x,T;{{q}^{k}})$ denote the solution to problem (3.5) with the potential q^k at t = T. Since the strong maximum principle is still valid for the time fractional diffusion equation [110], the scheme is monotonically convergent, under suitable conditions.

As the terminal time $T\to \infty$, the problem recovers a steady-state problem, and the scheme amounts to twice numerical differentiation in space and converges within one iteration, provided that the data g is accurate enough. Hence, it is natural to expect that the convergence of the scheme will depends crucially on the time T: the larger is the time T, the closer is the solution u to the steady state solution; and thus the faster is the convergence of the fixed point scheme. By lemma 2.1, as the fractional order α approaches zero, the solution u decays much faster around t = 0 than the classical one, i.e., $\alpha =1$. In other words, the fractional diffusion problem can reach a 'quasi-steady state' much faster than the classical one, especially for α close to zero, and the scheme will then converge much faster.

To illustrate the point, we present in figure 10 some numerical results of reconstructing a discontinuous potential $q=1+2x{{\chi }_{[0,0.5]}}+2(1-{{\chi }_{(0.5,1]}}x)$ (with ${{\chi }_{S}}$ being the characteristic function of the set S). In order to illustrate the convergence behavior of the fixed point scheme we take exact data. In the figure, e denotes the relative ${{L}^{2}}(\Omega )$ error. The numerical results fully confirm the preceding discussions: at a fixed time T, the smaller is the fractional order α, the faster is the convergence; and at fixed α, the larger is the time T, the faster is the convergence. Numerically, one also observes the monotone convergence of the scheme.

Zoom In Zoom Out Reset image size

Generally, the recovery of a coefficient in FDEs has not been extensively studied. Cheng et al [14] established the unique recovery of the fractional order $\alpha \in (0,1)$ and the diffusion coefficient from the lateral boundary measurements. It represents one of the first mathematical works on invere problems for FDEs, and has inspired many follow-up works. Yamamoto and Zhang [109] established conditional stability in determining a zeroth-order coefficient in a one-dimensional FDE with one half order Caputo derivative by a Carleman estimate. Carleman estimates for time fractional diffusion were discussed in [13, 62, 108]. In [73], the unique determination of the spatial coefficient and/or the fractional order from the data on a subdomain was shown for a positive initial condition. Wang and Wu [97] studied the simultaneous recovery of two time varying coefficients, i.e., a kernel function and a source function, from the additional integral observation in multi-dimension, using a fixed point theorem. All these works are concerned with the theoretical analysis, ant there are even fewer works on the numerical analysis of related inverse problems. Li et al [58] suggested an optimal perturbation algorithm for the simultaneous numerical recovery of the diffusion coefficient and fractional order in a one-dimensional time fractional FDE. In [50], the authors considered the identification of a potential term from the lateral flux data at one fixed time instance corresponding to a complete set of source terms, and established the unique determination for 'small' potentials. Further, a Newton type method was proposed in [50], and its convergence was shown.

Even though our discussions have focused on time fractional diffusion, which involves one single fractional derivative in time, it is also possible to consider equations where the time derivative involves multiple factional orders, i.e., $\sum _{k=1}^{m}{{c}_{k}}\partial _{t}^{{{\alpha }_{k}}}$ for a sequence ${{\alpha }_{1}}\gt {{\alpha }_{2}}\gt \ldots \gt {{\alpha }_{m}}$ [44, 60]; see [59] for some first uniqueness results for inverse coefficient problems in the multi-dimensional case. Further extensions include the distributed-order, spatially and/or temporally variable-order and tempered fractional diffusion, to better capture certain physical processes, for which, however, related inverse problems have not been discussed at all.

One of the very first undetermined coefficient problems for PDEs was discussed in the paper by Jones [52] (see also [8, chapter 13]). This is to determine the coefficient a(t) from

$$\begin{eqnarray*}\begin{array}{rcl} {{u}_{t}} & = & a(t){{u}_{xx}},\quad 0\lt x\lt \infty ,\quad t\gt 0 \\ u(x,0) & = & 0,\quad -a(t){{u}_{x}}(0,t)=g(t),\quad 0\lt t\lt T \\ \end{array}\end{eqnarray*}$$

under the over-posed condition of measuring the 'temperature' at x = 0

$$\begin{eqnarray*}&&u(0,t)=\psi (t)\end{eqnarray*}$$

In [52], Jones provided a complete analysis of the problem, by giving necessary and sufficient conditions for a unique solution as well as determining the exact level of ill-conditioning. The key step in the analysis is a change of variables and conversion of the problem to an equivalent integral equation formulation. Perhaps surprisingly, this approach involves the use of a fractional derivative as we now show.

The assumptions are that g is continuous and positive and ψ is continuously differentiable with $\psi (0)=0$ and $\psi ^{\prime} \gt 0$ on $(0,T)$. In addition, the function h(t) defined by

$$\begin{eqnarray*}&&h(t)=\frac{\sqrt{\pi }g(t)}{\int _{0}^{t}{{(t-\tau )}^{-1/2}}\psi ^{\prime} (\tau ){\rm d}\tau }\end{eqnarray*}$$

satisfies ${{{\rm lim} }_{t\to 0}}h(t)={{h}_{0}}\gt 0$. Note that h is the ratio of the two data functions; the flux g and the Djrbashian–Caputo derivative of order $1/2$ of ψ. If we define ${{h}_{i}}={\rm inf} h$ and ${{h}_{s}}={\rm sup} h$ on $[0,T]$ and look at the space $\mathcal{G}\;:=\;\{a\in C[0,T):h_{i}^{2}\leqslant a(t)\leqslant h_{s}^{2}\}$, then it was shown that any $a\in \mathcal{G}$ satisfies the inverse problem must also solve the integral equation

$$\begin{eqnarray*}&&a(t)=\frac{\sqrt{\pi }g(t)}{\int _{0}^{t}\psi ^{\prime} (\tau ){{[\int _{\tau }^{t}a(s){\rm d}s]}^{-1/2}}{\rm d}\tau }=:\mathcal{T}a,\end{eqnarray*}$$

and vice-versa. The main result in [52] is that the operator $\mathcal{T}$ has a unique fixed point on $\mathcal{G}$ and indeed $\mathcal{T}$ is monotone in the sense of preserving the partial order on $\mathcal{G}$, i.e., if ${{a}_{1}}\geqslant {{a}_{2}}$ then $\mathcal{T}{{a}_{1}}\leqslant \mathcal{T}{{a}_{2}}$.

Given these developments, it might seem that a parallel construction for the time fractional diffusion counterpart, $\partial _{t}^{\alpha }u=a(t){{u}_{xx}}$, would be relatively straightforward but this seems not to be the case. The basic steps for the parabolic version require items that just are not true in the fractional case, such as the product rule, and without these the above structure cannot be replicated or at least not without some further ingenuity.

First we consider the following Sturm–Liouville problem on the unit interval $\Omega =(0,1)$: find $u\in H_{0}^{1}(\Omega ){\mkern 1mu} \cap {{H}^{\beta }}(\Omega )$ and $\lambda \in \mathbb{C}$ such that

$$\begin{eqnarray}&&-{}_{0}^{C}D_{x}^{\beta }u+qu=\lambda u\quad {\rm in}\;\Omega ,\end{eqnarray} \tag{ 4.1 }$$

with a homogeneous Dirichlet boundary condition $u(0)=u(1)=0$. A Sturm–Liouville problem of this form was considered by Djrbashian [19, 22] in 1960s to construct certain biorthogonal basis for spaces of analytic functions; see also [78]. Like before, with $\beta =2$, it recovers the classical Sturm–Liouville problem. In the case of a general potential q, in the fractional case, little is known about the analytical properties of the eigenvalues and eigenfunctions. For the case of a zero potential q = 0, there are countably many eigenvalues $\{{{\lambda }_{j}}\}$ to (4.1), which are zeros of the Mittag-Leffler function ${{E}_{\beta ,2}}(-\lambda )$. The corresponding eigenfunctions are given by $x{{E}_{\beta ,2}}(-{{\lambda }_{j}}{{x}^{\beta }})$. Using the exponential asymptotics on the Mittag-Leffler function in lemma 2.1, one [50, 93] can show that asymptotically, the eigenvalues ${{\lambda }_{j}}$ are distributed as

$$\begin{eqnarray*}&&|{{\lambda }_{j}}|\sim {{(2\pi j)}^{\beta }}\quad {\rm and}\quad {\rm arg} ({{\lambda }_{j}})\sim \frac{(2-\beta )\pi }{2}.\end{eqnarray*}$$

Hence, for any $\beta \in (1,2)$, there are only a finite number of real eigenvalues to (4.1), and the rest appears as complex conjugate pairs.

It is well known that eigenvalues contain valuable information about the boundary value problem. For example it is known that the sequence of Dirichlet eigenvalues can uniquely determine a potential q symmetric with respect to the point $x=1/2$, and together with additional spectral information, one can uniquely determine a general potential q; see [12, 86] for an overview of results on the classical inverse Sturm–Liouville problem. In the fractional case, the eigenvalues are generally genuinely complex, and a complex number may carry more information than a real one. Thus one naturally wonders whether these complex eigenvalues do contain more information about the potential. Numerically the answer is affirmative. To illustrate this, we show some numerical reconstructions in figure 11, obtained by using a frozen Newton method and representing the sought-for potential q in Fourier series [50]. The Dirichlet eigenvalues can be computed efficiently using a Galerkin finite element method [45]. One observes that one single Dirichlet spectrum can uniquely determine a general potential q. Unsurprisingly, as the fractional order β tends two, the reconstruction becomes less and less accurate, since in the limit $\beta =2$, the Dirichlet spectrum cannot uniquely determine a general potential q. Theoretically, the surprising uniqueness in the fractional case remains to be established.

Zoom In Zoom Out Reset image size

Naturally, one can also consider the Riemann–Liouville case:

$$\begin{eqnarray}&&-{}_{0}^{R}D_{x}^{\beta }u+qu=\lambda u\quad \;{\rm in}\;\Omega ,\end{eqnarray} \tag{ 4.2 }$$

with $u(0)=u(1)=0$. Like before, little is known about the analytical properties of the eigenvalues and eigenfunctions. In the case of a zero potential q = 0, there are countably many eigenvalues to (4.2), which are zeros of the Mittag-Leffler function ${{E}_{\beta ,\beta }}(-\lambda )$, and the corresponding eigenfunctions are given by ${{x}^{\beta -1}}{{E}_{\beta ,\beta }}(-{{\lambda }_{j}}{{x}^{\beta }})$. Further, the asymptotics of the eigenvalues are still valid. Hence, for any $\beta \in (1,2)$, there are only a finite number of real eigenvalues to (4.2), and the rest appears as complex conjugate pairs.

The numerical results from the Dirichlet spectrum in the Riemann–Liouville case are shown in figure 12. For a general potential q, the reconstruction represents only the symmetric part, which is drastically different from the Djrbashian–Caputo case, but identical with that for the classical Sturm–Liouville problem. Further, if we assume that the potential q is known on the left half interval, then the Dirichlet spectrum allows uniquely reconstructing the potential q on the remaining half interval, cf figure 12(b). These results indicate that in the Riemann–Liouville case the complex spectrum is not more informative than the classical Sturm–Liouville problem. The precise mechanism underlying the fundamental differences between the Djrbashian–Caputo and Riemann–Liouville cases awaits further study. However, as in the classical case $\beta =2$, one can show that the linearized derivative of the map $q\to u(1;\lambda ,q)$ around q = 0 cannot span more than the subspace of even functions in ${{L}^{2}}(\Omega )$.

Zoom In Zoom Out Reset image size

In general, the Sturm–Liouville problem with a fractional derivative remains completely elusive, and numerical methods such as finite element method [46] provide a valuable (and often the only) tool for studying its analytical properties. For a variant of the fractional Sturm–Liouville problem, which contains a fractional derivative in the lower-order term, Malamud [71] established the existence of a similarity transformation, analogous to the well-known Gel'fand–Levitan–Marchenko transformation, and also the unique recovery of the potential from multiple spectra. In the classical case, the Gel'fand–Levitan–Marchenko transformation lends itself to a constructive algorithm [86]; however, it is unclear whether this is true in the fractional case. In [50], the authors proposed a Newton type method for reconstructing the potential, which numerically exhibits very good convergence behavior. However, a rigorous convergence analysis of the scheme is still missing. Further, the uniqueness and nonuniqueness issues of related inverse Sturm–Liouville problems are outstanding. Last, as noted above, there are other possible choices of the space fractional derivative, e.g., fractional Laplacian and Riesz derivative. It is unknown whether the preceding observations are valid for these alternative derivatives.

One classical elliptic inverse problem is the Cauchy problem for the Laplace equation, which plays a fundamental role in the study of many elliptic inverse problems [40]. A first example was given by Jacques Hadamard [31] to illustrate the severe ill-posedness of the Cauchy problem, which motivated him to introduce the concept of well-posedness and ill-posedness for problems in mathematical physics. So a natural question is whether the Cauchy problem for the fractional elliptic equation is also as ill-posed? To illustrate this, we consider the following fractional elliptic problem on the rectangular domain $\Omega =\{(x,y)\in {{\mathbb{R}}^{2}}:0\lt x\lt 1,0\lt y\lt 1\}$

$$\begin{eqnarray}&&{}_{0}^{C}D_{x}^{\beta }u+{}_{0}^{C}D_{y}^{\beta }u=0\;{\rm in}\;\Omega ,\end{eqnarray} \tag{ 4.3 }$$

with the fractional order $\beta \in (1,2)$ and the Cauchy data

$$\begin{eqnarray*}&&u(x,0)=g(x)\quad {\rm and}\quad \frac{\partial u}{\partial \nu }(x,0)=h(x),\;\;0\lt x\lt 1,\end{eqnarray*}$$

where ν is the unit outward normal direction. With $\beta =2$, it recovers the Cauchy problem for the Laplace equation. By applying the separation of variables, we can assume that $u(x,y)=\phi (x)\psi (y)$, which directly gives for some scalar $\lambda \in \mathbb{C}$ that

$$\begin{eqnarray*}\begin{array}{rcl} {}_{0}^{C}D_{x}^{\beta }\phi (x) & = & -\lambda \phi (x), \\ {}_{0}^{C}D_{y}^{\beta }\psi (y) & = & \lambda \psi (y), \\ \end{array}\end{eqnarray*}$$

Let $({{\lambda }_{j}},{{\phi }_{j}})$ be a Dirichlet eigenpair of the Caputo derivative operator $-{}_{0}^{C}D_{x}^{\beta }$ on the unit interval $D=(0,1)$, i.e., ${{\phi }_{j}}(x)=x{{E}_{\beta ,2}}(-{{\lambda }_{j}}{{x}^{\beta }})$, and $|{{\lambda }_{j}}|\to \infty$ as $j\to \infty$ [51]; see section 4.1 for further details. With the choice $\phi ={{\phi }_{j}}$ and the Cauchy data pair $(g,{{h}_{j}})=(0,-x{{E}_{\beta ,2}}({{\lambda }_{j}}{{x}^{\beta }})/{{\lambda }_{j}})$, the component ${{\psi }_{j}}$ satisfies

$$\begin{eqnarray*}&&{}_{0}^{C}D_{y}^{\beta }{{\psi }_{j}}(y)={{\lambda }_{j}}{{\psi }_{j}}(y)\;{\rm in}\;y\in (0,\infty ),\end{eqnarray*}$$

with the initial condition ${{\psi }_{j}}(0)=0$ and $\frac{{\rm d}}{{\rm d}y}{{\psi }_{j}}(0)=1/{{\lambda }_{j}}$. Using the relation $\frac{{\rm d}}{{\rm d}x}{{x}^{\gamma -1}}{{E}_{\beta ,\gamma }}(\lambda {{x}^{\beta }})=\lambda {{x}^{\gamma -2}}{{E}_{\beta ,\gamma -1}}(\lambda {{x}^{\beta }})$ [53, p 46], we deduce that the solution ${{\psi }_{j}}$ to the fractional ordinary differential equation is given by

$$\begin{eqnarray*}&&{{\psi }_{j}}(y)=y{{E}_{\beta ,2}}({{\lambda }_{j}}{{y}^{\beta }})/{{\lambda }_{j}}.\end{eqnarray*}$$

Hence, ${{u}_{j}}(x,y)=x{{E}_{\beta ,2}}(-{{\lambda }_{j}}{{x}^{\beta }})y{{E}_{\beta ,2}}({{\lambda }_{j}}{{y}^{\beta }})/\lambda _{j}^{2}$ is a solution to the Cauchy problem with g = 0 and ${{h}_{j}}(x)=-x{{E}_{\beta ,2}}(-{{\lambda }_{j}}{{x}^{\beta }})/{{\lambda }_{j}}$. By the exponential asymptotics of the Mittag-Leffler function, cf lemma 2.1, we deduce that ${{h}_{j}}(x)\to 0$ as $j\to \infty$, whereas for any $y\gt 0$, the solution ${{u}_{j}}(x,y)\to \infty$ as $j\to \infty$, in view of the exponential growth of the Mittag-Leffler function ${{E}_{\beta ,2}}(z)$, cf lemma 2.1. This indicates that the Cauchy problem for the fractional elliptic equation is also exponentially ill-posed. However, the interesting question of the degree of ill-posedness, in comparison withthe classical case, is unclear and certainly worthy of further study. Further, we note that the numerical solution of the fractional elliptic equation (4.3) is highly nontrivial, and there seems to be no efficient yet rigorous solver in the literature and this seems to be due to a lack of theory about the solution to such problems.

Now we return to the backward diffusion problem with fractional derivatives in the space variable(s). Let $\Omega =(0,1)$ be the unit interval. Then the one-dimensional space fractional diffusion equation is given by

$$\begin{eqnarray*}&&{{u}_{t}}-{}_{0}^{C}D_{x}^{\beta }u=0,\quad (x,t)\in \Omega \times (0,\infty ),\end{eqnarray*}$$

where the fractional order $\beta \in (1,2)$. The equation is equipped with the following initial condition $u(x,0)=v$ and zero boundary condition $u(0,t)=u(1,t)=0$. The backward problem is: given the final time data $g(x)=u(x,T)$, find the initial data v. Since the Djrbashian–Caputo derivative operator ${}_{0}^{C}D_{x}^{\beta }$ with the zero Dirichlet boundary is sectorial on suitable spaces [42], the existence of a solution u follows from the analytic semigroup theory [43, 81], and formally it can be represented by

$$\begin{eqnarray*}&&u(t)={{{\rm e}}^{-At}}v,\end{eqnarray*}$$

where A is the representation of the Djrbashian–Caputo derivative operator $-{}_{0}^{C}D_{x}^{\beta }$ on its domain. Formally, the solution v to the space fractional backward problem is given by

$$\begin{eqnarray*}&&v={{{\rm e}}^{AT}}g.\end{eqnarray*}$$

In case of $\beta =2$, using the eigenpairs $\{({{\lambda }_{j}},{{\phi }_{j}})\}$ and the ${{L}^{2}}(\Omega )$ orthogonality of the eigenfunctions $\{{{\phi }_{j}}\}$, it recovers the well known formula

$$\begin{eqnarray*}&&v=\mathop{\sum }\limits_{j=1}^{\infty }{{{\rm e}}^{{{\lambda }_{j}}T}}(g,{{\phi }_{j}}){{\phi }_{j}}.\end{eqnarray*}$$

The growth factor ${{{\rm e}}^{{{\lambda }_{j}}T}}$ explains the severely ill-posed nature of the inverse problem. In the fractional case, such an explicit representation is no longer available since the corresponding eigenfunctions $\{{{\phi }_{j}}\}$ are not orthogonal in ${{L}^{2}}(\Omega )$ (actually they can be almost linearly dependent), due to the non self adjoint nature of the Djrbashian–Caputo derivative operator $-{}_{0}^{C}D_{x}^{\beta }$. Nonetheless, according to the discussions in section 4.1, the eigenvalues $\{{{\lambda }_{j}}\}$ increase to infinity with the index j, and asymptotically lies on two rays. Hence, one naturally expects that the backward problem is also exponentially ill-posed. However, the magnitudes (and the real parts) of the eigenvalues grow at a rate slower than that of the standard Sturm–Liouville problem, and thus the space fractional backward problem is less ill-posed than the classical one. To illustrate the point, we present the numerical results in figure 13. For all fractional orders β, the singular values decay exponentially, but the decay rate increases dramatically with the increase of the fractional order β and the terminal time T. Hence, anomalous superdiffusion does not change the exponentially ill-posed nature of the backward problem, but numerically it does enable recovering more Fourier modes of the initial data v.

Zoom In Zoom Out Reset image size

Last, we note that for other choices of the fractional derivative, e.g., the Riemann–Liouville fractional derivative and the fractional Laplacian [6, 55], the magnitude of eigenvalues of the operator also tends to infinity, and the growth rate increases with the fractional order β. Therefore, the preceding observations on the space fractional backward problem are expected to be valid for these choices as well.

Last we return to the classical sideways diffusion problem but now with a fractional derivative in space rather than in time. Let $\Omega =(0,1)$ be the unit interval. Then the one-dimensional space fractional diffusion equation is given by

$$\begin{eqnarray*}&&{{u}_{t}}-{}_{0}^{C}D_{x}^{\beta }u=0,\quad (x,t)\in \Omega \times (0,\infty ),\end{eqnarray*}$$

where the fractional order $\beta \in (1,2)$. The equation is equipped with an initial condition $u(x,0)=0$ and the following lateral Cauchy boundary conditions

$$\begin{eqnarray*}&&u(0,t)=f(t)\ \quad {\rm and}\quad \ {{u}_{x}}(0,t)=g(t),\ \ t\gt 0.\end{eqnarray*}$$

We wish to compute the solution at x = 1, i.e., $h(t)\;:=\;u(1,t)$. In the case $\beta =2$, the model recovers the standard diffusion equation, and we have already discussed the severe ill-conditioning of the classical case. Due to the nonlocal nature of the fractional derivative, one might expect that in the space fractional case, the sideways problem is less ill-posed. To see this, we take Laplace transform in time to arrive at (with $\mathop{\;}\limits^{\widehat{{}}}$ denoting the Laplace transform)

$$\begin{eqnarray*}&&z\hat{u}(x,z)-{}_{0}^{C}D_{x}^{\beta }\hat{u}(x,z)=0,\end{eqnarray*}$$

with the initial conditions (at x = 0)

$$\begin{eqnarray*}&&\hat{u}(0,z)=\hat{f}(z)\quad {\rm and}\quad {{\hat{u}}_{x}}(0,z)=\hat{g}(z).\end{eqnarray*}$$

The solution $\hat{u}(x,z)$ to the initial value problem is given by

$$\begin{eqnarray*}&&\hat{u}(x,z)=\hat{f}(z){{E}_{\beta ,1}}(z{{x}^{\beta }})+\hat{g}(z)x{{E}_{\beta ,2}}(z{{x}^{\beta }})\end{eqnarray*}$$

and thus

$$\begin{eqnarray*}&&\hat{h}(z)=\hat{f}(z){{E}_{\beta ,1}}(z)+\hat{g}(z){{E}_{\beta ,2}}(z).\end{eqnarray*}$$

Like before, the boundary condition h(t) at x = 1 can be found by an inverse Laplace transform

$$\begin{eqnarray*}&&h(t)={{\int }_{{\rm Br}}}{{{\rm e}}^{zt}}\hat{h}(z){\rm d}z.\end{eqnarray*}$$

In case of $\beta =2$, this gives ${\rm cosh} \sqrt{z}$ and ${\rm sinh} \sqrt{z}/\sqrt{z}$ multipliers to the data $\hat{f}(z)$ and $\hat{g}(z)$ resulting in the exponential ill-conditioning of the sideways heat problem. In the case of a general $\beta \in (1,2)$, the exponential asymptotics in lemma 2.1 indicates that the problem still suffers from exponentially growing multipliers to the data, and thus the problem is still severely ill-conditioned. Simple computation shows that the multiplier is asymptotically larger for the fractional order β closer to unity. In other words, anomalous diffusion in space does not mitigate the ill-conditioned nature of the sideways problem, but actually worsens the conditioning severely.

To further illustrate the point, we compute the forward map F from the Dirichlet boundary condition at x = 1 to the flux at x = 0 numerically with a finite element in space and finite difference in time scheme, cf appendix A.3 for the details. The numerical results are presented in figure 14. The singular value spectra clearly show the ill-posedness nature of the space fractional sideways problem: as the fractional order β increases from one to two, the majority of the singular values move upward, the decay of the singular values slows down, and thus the sideways problem becomes less and less ill-posed (but still severely so). Further, there are more tiny singular values kicking in as the fractional order β decreases to one, which indicates the inherent rank deficiency of the forward map F and might be relevant in the uniqueness of the inverse problem. This confirms the preceding analysis: the degree of ill-posedness worsens with the decrease of the fractional order β, and the fractional counterpart is more ill-posed than the classical one. In other words, anomalous diffusion actually severely worsens the conditioning of the already very ill-posed sideways problem. Further, the numerical results tend to indicate that the Djrbashian–Caputo derivative with an order $\beta \in (1,2)$ acts as an interpolation between the diffusion and convection, which results in a history mechanism in space: when the history piece runs from the left to the right, it is unlikely to propagate the information in the reverse direction; and the closer is the fractional order β to unity, the stronger is the directional effect. The latter is not counterintuitive, since in the limit of $\beta =1$, the Djrbashian–Caputo fractional derivative ${}_{0}^{C}D_{x}^{\beta }\;u$ recovers the first order derivative $\frac{\partial u}{\partial x}$, and the problem is of convection type, and surely no information can be convected backwards!

Zoom In Zoom Out Reset image size

In the case $\beta =2$, one may equally measure the lateral Cauchy data at x = 1, and aims at recovering the Dirichlet boundary condition at x = 0. Clearly, this does not change the nature of the inverse problem, and it is equally ill-posed. Due to the directional nature of the Djrbashian–Caputo derivative ${}_{0}^{C}D_{x}^{\beta }$, one naturally wonders whether this 'directional' feature does influence the ill-posed nature of the sideways problem. To illustrate the point, we repeat the preceding arguments and deduce

$$\begin{eqnarray*}&&z\hat{u}(x,z)-{}_{0}^{C}D_{x}^{\beta }\hat{u}(x,z)=0,\end{eqnarray*}$$

with the boundary conditions at x = 1

$$\begin{eqnarray*}&&\hat{u}(1,z)=\hat{f}(z)\quad {\rm and}\quad {{\hat{u}}_{x}}(1,z)=\hat{g}(z).\end{eqnarray*}$$

To derive the solution, denote the initial conditions at x = 0 by $\tilde{f}(z)=\hat{u}(0,z)$ and $\tilde{g}(z)={{\hat{u}}_{x}}(0,z)$. Then the solution $\hat{u}(x,z)$ to the initial value problem is given by

$$\begin{eqnarray*}&&\hat{u}(x,z)=\tilde{f}(z){{E}_{\beta ,1}}(z{{x}^{\beta }})+\tilde{g}x{{E}_{\beta ,2}}(z{{x}^{\beta }}).\end{eqnarray*}$$

Use the differentiation formula $\frac{{\rm d}}{{\rm d}x}{{x}^{\gamma -1}}{{E}_{\beta ,\gamma }}(z{{x}^{\beta }})=z{{x}^{\gamma -2}}{{E}_{\beta ,\gamma -1}}(z{{x}^{\beta }})$ [53, p 46] we deduce that at x = 1, there hold

$$\begin{eqnarray*}\begin{array}{rcl} \tilde{f}(z){{E}_{\beta ,1}}(z)+\tilde{g}{{E}_{\beta ,2}}(z) & = & \hat{f}(z), \\ \tilde{f}(z){{E}_{\beta ,0}}(z)+\tilde{g}{{E}_{\beta ,1}}(z) & = & {{z}^{-1}}\hat{g}(z). \\ \end{array}\end{eqnarray*}$$

Solving the linear system yields the solution to the sideways problem

$$\begin{eqnarray*}&&\tilde{f}(z)=\frac{{{E}_{\beta ,1}}(z)\hat{f}(z)-{{z}^{-1}}{{E}_{\beta ,2}}(z)\hat{g}(z)}{{{E}_{\beta ,1}}{{(z)}^{2}}-{{E}_{\beta ,0}}(z){{E}_{\beta ,2}}(z)},\end{eqnarray*}$$

and accordingly the solution $h(t)\equiv u(0,t)$ is given by an inverse Laplace transform. The growth factors of the data $\hat{f}$ and $\hat{g}$ are ${{E}_{\beta ,1}}(z)/({{E}_{\beta ,1}}{{(z)}^{2}}-{{E}_{\beta ,0}}(z){{E}_{\beta ,2}}(z))$ and ${{z}^{-1}}{{E}_{\beta ,2}}(z)/({{E}_{\beta ,1}}{{(z)}^{2}}-{{E}_{\beta ,0}}(z){{E}_{\beta ,2}}(z))$, respectively. The growth of these factors at large z argument determines the degree of ill-conditioning of the sideways problem. To this end, we appeal to the exponential asymptotic of the Mittag-Leffler function ${{E}_{\alpha ,\beta }}(z)$, cf lemma 2.1, and note that Bromwhich path lies in the sector $|{\rm arg} \;z|\leqslant \pi /2$ to deduce that for large $|z|$,there holds

$$\begin{eqnarray*}\begin{array}{rcl} \quad {{E}_{\beta ,1}}{{(z)}^{2}} & \sim & \frac{1}{{{\beta }^{2}}}{{{\rm e}}^{2{{z}^{1/\beta }}}}-\frac{2}{\beta \Gamma (1-\beta )z}{{{\rm e}}^{{{z}^{1/\beta }}}}, \\ {{E}_{\beta ,0}}(z){{E}_{\beta ,2}}(z) & \sim & \frac{1}{{{\beta }^{2}}}{{{\rm e}}^{2{{z}^{1/\beta }}}}-\frac{1}{\beta \Gamma (2-\beta )}{{z}^{1/\beta -1}}{{{\rm e}}^{{{z}^{1/\beta }}}}. \\ \end{array}\end{eqnarray*}$$

Hence, the numerator ${{E}_{\beta ,1}}{{(z)}^{2}}-{{E}_{\beta ,0}}(z){{E}_{\beta ,2}}(z)$ behaves like

$$\begin{eqnarray*}&&{{E}_{\beta ,1}}{{(z)}^{2}}-{{E}_{\beta ,0}}(z){{E}_{\beta ,2}}(z)\sim \frac{1}{\beta \Gamma (2-\beta )}{{z}^{1/\beta -1}}{{{\rm e}}^{{{z}^{1/\beta }}}}\qquad {\rm as}\;|z|\to \infty .\end{eqnarray*}$$

This together with the exponential asymptotic of ${{E}_{\beta ,1}}(z)$ and ${{E}_{\beta ,2}}(z)$ from lemma 2.1 indicates that the multipliers for $\hat{f}$ and $\hat{g}$ are growing at most at a very low-order polynomial rate, for large z. Hence, the high-frequency components of the data noise are not amplified much (at most polynomially instead of exponentially). The analysis indicates that the sideways problem with the lateral Cauchy data specified at the point x = 1 is nearly well-posed, as long as the fractional order β is away from two, for which it recovers the classical ill-posed sideways problem for the heat equation.

Next we illustrate the preceding discussions numerically. The behavior of the forward map F from the Dirichlet boundary at x = 0 to the flux data at x = 1 is shown in figure 15. For a wide range of values of the fractional order β, the condition number of the forward map F is of order 100, which is fairly mild, in view of the size of the linear system, i.e., 100 × 100. When the fractional order β increases towards two, the inverse problem recovers the classical sideways problem, and as expected, the condition number increases dramatically. However, the onset of the blowup depends on the terminal time T: the smaller is the time T, the smaller seems the onset value. The precise mechanism for this phenomenon is still unknown. For $\beta \leqslant 7/4$, the singular value spectrum only spans a narrow interval, resulting in a very small condition number. Physically, like before, this can be explained as the 'convective' nature of the Djrbashian–Caputo fractional derivative: as the fractional order β tends to unity, the information at x = 0 is transported to x = 1, free from distortion, and thus the inverse problem is almost well-posed. In summary, depending on the location of the over-specified data, anomalous superdiffusion can either help or aggravate the conditioning of the sideways problem.

Zoom In Zoom Out Reset image size

Last, we would like to note that the study of space fractional inverse problems, either theoretical or numerical, is fairly scarce. This is partly attributed to the relatively poor understanding of forward problems for FDEs with a space fractional derivative: there are only a few mathematical studies on one-dimensional space fractional diffusion, and no mathematical study on multi-dimensional problems involving space fractional derivatives (of either Riemann–Liouville or Caputo type). Nonetheless, our preliminary numerical experiments show distinct new features for related inverse problems, which motivate their analytical studies.

Like many special functions, the efficient and accurate numerical computation of the Mittag-Leffler function ${{E}_{\alpha ,\beta }}(z)$ is delicate [28, 29, 94]. An efficient algorithm relies on partitioning the complex plane $\mathbb{C}$ into different regions, where different approximations, i.e., power series, integral representation and exponential asymptotic for small values of the argument, intermediate values and large values, respectively, are used for efficient numerical computation; see [94] for the some partition and error estimates. The special case of the Mittag-Leffler function ${{E}_{\alpha ,\beta }}(z)$ with a real argument $z\in \mathbb{R}$, which plays a predominant role in time-fractional diffusion, can also be efficiently computed with the Laplace transform and suitable quadrature rules [27].

The computation of the Wright function ${{W}_{\rho ,\mu }}(z)$ is even more delicate. In theory, like before, it can be computed using power series for small values of the argument and a known asymptotic formula for large values, and for the intermediate case, values are obtained by using an integral representation [67]. The integral representation for the Wright function ${{W}_{\rho ,\mu }}(z)$ for intermediate values in the case of interest for the fundamental solution in one-dimension (where $\rho =-\alpha /2\lt 0$, $0\lt \mu =1+\rho \lt 1$ and $z=-x$, $x\gt 0$) is given by

$$\begin{eqnarray*}&&{{W}_{\rho ,\mu }}(-x)=\int _{0}^{\infty }K(x,\rho ,\mu ,r){\rm d}r\end{eqnarray*}$$

where the kernel $K(x,\rho ,\mu ,r)$ is given by

$$\begin{eqnarray*}&&K(x,\rho ,\mu ,r)={{r}^{-\mu }}{{{\rm e}}^{-r+x{\rm cos} (\pi \rho ){{r}^{-\rho }}}}{\rm sin} (x{\rm sin} (\pi \rho ){{r}^{-\rho }}+\pi \mu ).\end{eqnarray*}$$

This is a singular kernel with a leading order ${{r}^{-\mu }}={{r}^{-1-\rho }}$, with successive singular kernels of the form ${{r}^{-1-2\rho }}$, ${{r}^{-1-3\rho }}$ etc, upon expanding the terms. Hence, a direct treatment via numerical quadrature is inefficient. A more efficient approach is to use the change of variable $s={{r}^{-\rho }}$, i.e., $r={{s}^{-1/\rho }}$, and the transformed kernel is

$$\begin{eqnarray*}&&\tilde{K}(x,\rho ,\mu ,s)={{(-\rho )}^{-1}}{{s}^{{{(-\rho )}^{-1}}(-\mu +1)-1}}{{{\rm e}}^{-{{s}^{-\frac{1}{\rho }}}+x{\rm cos} (\pi \rho )s}}{\rm sin} (x{\rm sin} (\pi \rho )s+\pi \mu ).\end{eqnarray*}$$

The fundamental solution of the one-dimensional time-fractional diffusion equation is expressed in terms of a Wright function ${{W}_{\rho ,\mu }}(-x)$ with the choice $\rho =-\alpha /2$ and $\mu =1+\rho$, cf (2.2). In this case the resulting kernel $\tilde{K}$ simplifies to

$$\begin{eqnarray*}&&\tilde{K}(x,\rho ,1+\rho ,s)={{(-\rho )}^{-1}}{{{\rm e}}^{-{{s}^{-\frac{1}{\rho }}}+x{\rm cos} (\pi \rho )s}}{\rm sin} (x{\rm sin} (\pi \rho )s+(1+\rho )\pi ).\end{eqnarray*}$$

This kernel is free from the grave singularity, and thus the quadrature method is quite effective. In general, the integral can be computed efficiently via the Gauss–Jacobi quadrature, with the weight function ${{s}^{{{(-\rho )}^{-1}}(-\mu +1)-1}}$. We note that an algorithm for the Wright function ${{W}_{\rho ,\mu }}(z)$ over the whole complex plane $\mathbb{C}$ with rigorous error analysis is still missing. The endeavor in this direction would almost certainly involve dividing the complex domain $\mathbb{C}$ into different regions, and using different approximations on each region separately.

We describe a finite difference method for the initial boundary value problem for the one-dimensional time-fractional diffusion equation

$$\begin{eqnarray*}&&\partial _{t}^{\alpha }u-{{u}_{xx}}+qu=f(x,t)\quad (x,t)\in \Omega \times (0,T),\end{eqnarray*}$$

with the initial condition $u(x,0)=v$ and boundary conditions

$$\begin{eqnarray*}&&u(0,t)=g(t)\quad {\rm and}\quad u(1,t)=h(t),\quad t\gt 0.\end{eqnarray*}$$

There are many efficient numerical schemes for discretizing the problem. The discretization in space can be achieved by the standard central difference scheme, Galerkin finite element method [47] or spectral method, and the discretization in time can be achieved with the L1 approximation [63, 96] and convolution quadrature [48]. We shall adopt the L1 approximation in time and the central difference in space. Specifically, we divide the interval $[0,T]$ into uniform subintervals, with nodes ${{t}_{k}}=k\tau$, k = 0,..., K, and a time step size $\tau =T/K$. Similarly, we partition the spatial domain Ω into uniform subintervals, with grid points ${{x}_{i}}={\rm i}h$, i = 0,...,N, and mesh size $h=1/N$. Then the L1-approximation of the Djrbashian–Caputo fractional derivative $\partial _{t}^{\alpha }u(x,{{t}_{k}})$ developed in [63, 96] is given by:

$$\begin{eqnarray}&&\partial _{t}^{\alpha }u(x,{{t}_{k}})\approx {{\tau }^{-\alpha }}\left[ {{b}_{0}}u(x,{{t}_{k}})-{{b}_{k-1}}u(x,{{t}_{0}})+\mathop{\sum }\limits_{j=1}^{k-1}({{b}_{j}}-{{b}_{j-1}})u(x,{{t}_{k-j}}) \right],\end{eqnarray} \tag{ A.1 }$$

where the weights b_j are given by

$$\begin{eqnarray*}&&{{b}_{j}}=({{(j+1)}^{1-\alpha }}-{{j}^{1-\alpha }})/\Gamma (2-\alpha ),\ j=0,1,\ldots ,K-1.\end{eqnarray*}$$

If the solution $u(x,t)$ is C² continuous in time, the local truncation error of the L1 approximation is bounded by $c{{\tau }^{2-\alpha }}$ for some c depending only on u [63, equation (3.3)]. In general, one can show that the scheme is only first-order accurate. Next with the central difference scheme in space and the notation $u_{i}^{k}\approx u({{x}_{i}},{{t}_{k}})$, we arrive at the following fully discrete scheme

$$\begin{eqnarray*}&&\begin{array}{l} \left[ {{b}_{0}}u_{i}^{k}-{{b}_{k-1}}u_{i}^{0}+\mathop{\sum }\limits_{j=1}^{k-1}({{b}_{j}}-{{b}_{j-1}})u_{i}^{k-j} \right] \\ \quad +\frac{u_{i-1}^{k}-2u_{i}^{k}+u_{i+1}^{k}}{{{h}^{2}}}+{{q}_{i}}u_{i}^{k}=f_{i}^{k},\;\;i=1,\ldots ,N-1, \\ \end{array}\end{eqnarray*}$$

with ${{q}_{i}}=q({{x}_{i}})$ and $f_{i}^{k}=f({{x}_{i}},{{t}_{k}})$. We note that at each time step, one needs to solve a tridiagonal linear system. However, the right-hand side at the current step involves all previous steps, which can be quite expensive for a small step size, and this will most likely be required due to the first order in time convergence. This history piece represents one of the main computational challenges for time fractional differential equations. There are high-order schemes, e.g., convolution quadrature generated by the second-order backward difference formula [48]. Further, we note that the finite difference scheme in space can be replaced with the Galerkin finite element method, which is especially suitable for high dimensional problems on a general domain and elliptic operator involving variable coefficients [47].

Now we describe a fully discrete scheme based on the backward Euler method in time and a Galerkin finite element method in space for the space fractional diffusion problem on the unit interval $\Omega =(0,1)$

$$\begin{eqnarray*}&&{{u}_{t}}-{}_{0}^{C}D_{x}^{\beta }\;u+qu=f\quad \;{\rm in}\;\Omega \times (0,T],\end{eqnarray*}$$

with the initial condition u = v and the Dirichlet boundary condition

$$\begin{eqnarray*}&&u(0,t)=g(t)\quad {\rm and}\quad u(1,t)=h(t),\;\;t\gt 0.\end{eqnarray*}$$

The Galerkin finite element method relies on the variational formulation for the fractional elliptic problem

$$\begin{eqnarray*}&&-{}_{0}^{C}D_{x}^{\beta }u+qu=f\quad \;{\rm in}\;\Omega ,\end{eqnarray*}$$

with a homogeneous Dirichlet boundary condition $u(0)=u(1)=0$, recently developed in [45]. The variational formulation of the problem is given by: find $u\in U\equiv H_{0}^{\beta /2}(\Omega )$ such that

$$\begin{eqnarray*}&&-({}_{0}^{R}D_{x}^{\beta /2}u,\;{}_{x}^{R}D_{1}^{\beta /2}v)+(qu,v)=(f,v)\quad \forall v\in V,\end{eqnarray*}$$

where ${}_{0}^{R}D_{x}^{\gamma }v$ and ${}_{x}^{R}D_{1}^{\gamma }v$ are the left-sided and right-sided Riemann–Liouville derivative of order $\gamma \in (0,1)$ defined by (with ${{c}_{\gamma }}=1/\Gamma (1-\gamma )$)

$$\begin{eqnarray*}&&{}_{0}^{R}D_{x}^{\gamma }v(x)={{c}_{\gamma }}\frac{{\rm d}}{{\rm d}x}\int _{0}^{x}{{(x-s)}^{-\gamma }}v(s){\rm d}s\quad {\rm and}\quad {}_{x}^{R}D_{1}^{\gamma }v(x)=-{{c}_{\gamma }}\frac{{\rm d}}{{\rm d}x}\int _{x}^{1}{{(x-s)}^{-\gamma }}v(s){\rm d}s,\end{eqnarray*}$$

respectively, and the test space V is given by

$$\begin{eqnarray*}&&V=\left\{ v\in {{H}^{\beta /2}}(\Omega ):\;v(1)=0,\;({{x}^{1-\beta }},v)=0 \right\}.\end{eqnarray*}$$

For the finite element discretization we first divide the unit interval Ω into a uniform mesh, with the grid points ${{x}_{i}}={\rm i}h$, i = 0,...,N and mesh size $h=1/N$. Then for ${{U}_{h}}\subset U$ we take the continuous piecewise linear finite element space, and for ${{V}_{h}}\subset V$ we construct it from U_h. Specifically, with the finite element basis ${{\phi }_{i}}(x)$, $i=1,\ldots ,N-1$, we take ${{\tilde{\phi }}_{i}}(x)={{\phi }_{i}}(x)-{{\gamma }_{i}}(1-x)\in {{V}_{h}}$, where the constant ${{\gamma }_{i}}$ is determined by the integral condition $({{x}^{1-\beta }},{{\tilde{\phi }}_{i}})=0$, i.e., ${{\gamma }_{i}}={{h}^{2-\beta }}({{({\rm i}-1)}^{3-\beta }}+{{({\rm i}+1)}^{3-\beta }}-2{{{\rm i}}^{3-\beta }})$. The computation of the leading term in the stiffness matrix and mass matrix can be carried out analytically, and the part involving the potential q can be computed efficiently using quadrature rules; see [46] for details.

Now for the time-dependent problem, like before, we divide the time interval $[0,T]$ into uniform subintervals, with ${{t}_{k}}=k\tau$, k = 0,...,K, and the time step size $\tau =T/K$. Then with the backward Euler method in time, and the finite element method in space, the approximate solution u_h^k at time t_k can be split into $u_{h}^{k}=\tilde{u}_{h}^{k}+{{s}^{k}}$ with the particular solution ${{s}^{k}}=g({{t}_{k}})+(h({{t}_{k}})-g({{t}_{k}}))x$ with the homogeneous solution $\tilde{u}_{h}^{k}\in {{U}_{h}}$ satisfying

$$\begin{eqnarray*}&&\begin{array}{l} {{\tau }^{-1}}(\tilde{u}_{h}^{k},{{v}_{h}})-({}_{0}^{R}D_{x}^{\beta /2}\tilde{u}_{h}^{k},\;{}_{x}^{R}D_{1}^{\beta /2}{{v}_{h}})+(q\tilde{u}_{h}^{k},{{v}_{h}}) \\ \qquad \quad =({{f}^{k}},{{v}_{h}})+{{\tau }^{-1}}(u_{h}^{k-1}-{{s}^{k}},{{v}_{h}})-(q{{s}^{k}},{{v}_{h}})\quad \forall {{v}_{h}}\in {{V}_{h}}. \\ \end{array}\end{eqnarray*}$$

We note the resulting linear system is of lower Hessenberg form, due to the nonlocality of the fractional derivative operator. However, the coefficient matrix does not change during the time stepping procedure, and thus an LU factorization might be applied to speedup the computation.