A tutorial on inverse problems for anomalous diffusion processes

Over the last two decades, anomalous diffusion processes in which the mean squares variance grows slower or faster than that in a Gaussian process have found many applications. At a macroscopic level, these processes are adequately described by fractional differential equations, which involves fractional derivatives in time or/and space. The fractional derivatives describe either history mechanism or long range interactions of particle motions at a microscopic level. The new physics can change dramatically the behavior of the forward problems. For example, the solution operator of the time fractional diffusion diffusion equation has only limited smoothing property, whereas the solution for the space fractional diffusion equation may contain weak singularity. Naturally one expects that the new physics will impact related inverse problems in terms of uniqueness, stability, and degree of ill-posedness. The last aspect is especially important from a practical point of view, i.e., stably reconstructing the quantities of interest. In this paper, we employ a formal analytic and numerical way, especially the two-parameter Mittag-Leffler function and singular value decomposition, to examine the degree of ill-posedness of several ‘classical’ inverse problems for fractional differential equations involving a Djrbashian–Caputo fractional derivative in either time or space, which represent the fractional analogues of that for classical integral order differential equations. We discuss four inverse problems, i.e., backward fractional diffusion, sideways problem, inverse source problem and inverse potential problem for time fractional diffusion, and inverse Sturm–Liouville problem, Cauchy problem, backward fractional diffusion and sideways problem for space fractional diffusion. It is found that contrary to the wide belief, the influence of anomalous diffusion on the degree of ill-posedness is not definitive: it can either significantly improve or worsen the conditioning of related inverse problems, depending crucially on the specific type of given data and quantity of interest. Further, the study exhibits distinct new features of ‘fractional’ inverse problems, and a partial list of surprising observations is given below. (a) Classical backward diffusion is exponentially ill-posed, whereas time fractional backward diffusion is only mildly ill-posed in the sense of norms on the domain and range spaces. However, this does not imply that the latter always allows a more effective reconstruction. (b) Theoretically, the time fractional sideways problem is severely ill-posed like its classical counterpart, but numerically can be nearly well-posed. (c) The classical Sturm–Liouville problem requires two pieces of spectral data to uniquely determine a general potential, but in the fractional case, one single Dirichlet spectrum may suffice. (d) The space fractional sideways problem can be far more or far less ill-posed than the classical counterpart, depending on the location of the lateral Cauchy data. In many cases, the precise mechanism of these surprising observations is unclear, and awaits further analytical and numerical exploration, which requires new mathematical tools and ingenuities. Further, our findings indicate fractional diffusion inverse problems also provide an excellent case study in the differences between theoretical ill-conditioning involving domain and range norms and the numerical analysis of a finite-dimensional reconstruction procedure. Throughout we will also describe known analytical and numerical results in the literature.


Introduction
Diffusion is one of the most prominent transport mechanisms found in nature. At a microscopic level, it is related to the random motion of individual particles, and the use of the Laplace operator and the first-order derivative in the canonical diffusion model rests on a Gaussian process assumption on the particle motion, after Albert Einsteinʼs groundbreaking work [23]. Over the last two decades a large body of literature has shown that anomalous diffusion models in which the mean square variance grows faster (superdiffusion) or slower (subdiffusion) than that in a Gaussian process under certain circumstances can offer a superior fit to experimental data (see the comprehensive reviews [5,70,72,95] for physical background and practical applications). For example, anomalous diffusion is often observed in materials with memory, e.g., viscoelastic materials, and heterogeneous media, such as soil, heterogeneous aquifer, and underground fluid flow. At a microscopic level, the subdiffusion process can be described by a continuous time random walk [75], where the waiting time of particle jumps follows some heavy tailed distribution, whereas the superdiffusion process can be described by Lévy flights or Lévy walk, where the length of particle jumps follows some heavy tailed distribution, reflecting the long-range interactions among particles. Following the aforementioned micro-macro correspondence, the macroscopic counterpart of a continuous time random walk is a differential equation with a fractional derivative in time, and that for a Lévy flight is a differential equation with a fractional derivative in space. We will refer to these two cases as time fractional diffusion and space fractional diffusion, respectively, and it is generically called a fractional derivative equation (FDE) below. In general the fractional derivative can appear in both time and space variables.
Next we give the mathematical model in the simplest geometrical setting of one space dimension, taking the domain Ω = (0, 1). Then a general, linear FDE is given by where > T 0 is a fixed time, and it is equipped with suitable boundary and initial conditions. The fractional orders α ∈ (0, 1) and β ∈ (1, 2) are related to the parameters specifying the large-time behavior of the waiting-time distribution or long-range behavior of the particle jump distribution. For example, in hydrological studies, the parameter β is used to characterize the heterogeneity of porous medium [17]. In theory, these parameters can be determined from the underlying stochastic model, but often in practice, they are determined from experimental data [  2 , respectively, for which the model (1.1) recovers the standard one-dimensional diffusion equation, and thus generally the model inverse problems. Throughout the notation c, with or without a subscript, denote a generic constant, which may differ at different occurrences, but it is always independent of the unknown of interest.

Preliminaries
We recall two important special functions, Mittag-Leffler function and Wright function, and one useful tool for analyzing discrete ill-posed problems, singular value decomposition. years since it appears as the fundamental solution for FDEs [69]. The Wright function ρ μ W z ( ) , has the following the asymptotic expansion in one sector containing the negative real axis  − [67, theorem 3.2]. Like before, the exponential asymptotics can be used to deduce the distribution of its zeros [66].
, π π − < ⩽ y arg ( ) , and for all small ϵ > 0, π ρ π ϵ ⩽ + − y |arg( )| min (3 (1 ) 2, ) . Then (1 ) and the coefficients A m , = … m 0, 1, are defined by the asymptotic expansion which is a Gaussian distribution in x for any > t 0. In the fractional case, the fundamental solution α K x t ( , ) exhibits quite different behavior than the heat kernel. To see this, we show the profile of α K x t ( , ) in figure 2 for several α values; see appendix A.1 for a brief description of the computational details. For any α < < 0 1 , α K x t ( , ) decays slower at a polynomial rate as the argument α x t | | 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 α < 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.

Singular value decomposition
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  ∈ × A n m , its singular value decomposition is given by 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 σ σ σ … ( , , , ) . 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   [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.

Inverse problems for time fractional diffusion
In this section, we consider several model inverse problems for the following one-dimensional time fractional diffusion equation on the unit interval Ω = (0, 1): with the fractional order α ∈ (0, 1), the initial condition = u v (0) and suitable boundary conditions. Although we consider only the one-dimensional model, the analysis and computation in this part can be extended into the general multi-dimensional case, upon suitable modifications. Recall that ∂ α u t denotes the Djrbashian-Caputo fractional derivative of order α with respect to time t. For α = 1, the fractional derivative ∂ α u t coincides with the usual first-order derivative ′ u , and accordingly, the model (3.1) reduces to the classical diffusion equation. Hence it is natural to compare inverse problems for the model (3.1) with that for the standard diffusion equation. We shall discuss the following four inverse problems, i.e., the backward problem, sideways problem, inverse source problem and inverse potential problem. In the first three cases, we shall assume a zero potential q = 0. We will also discuss the solution of an inverse coefficient problem using fractional calculus.  Therefore, the final time data = g u T ( ) is given by j of the (noisy) data g is amplified by an exponentially growing factor λ e T j , 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 λ e T j in order to recover the corresponding mode of the initial data v.

Backward fractional diffusion
In the fractional case, by lemma 2. , which is very mild compared to the exponential growth λ e T j for the case α = 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 λ j . More precisely, it amounts to the loss of two spatial derivatives [88, theorem 4.1] 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 g ( , ) j j , = … j J 1, 2, , , 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 − in t, it suffices to examine the Jth mode. For the heat equation the growth factor is a possibly workable value of around 600; while for the heat equation it is greater than 10 25 . We reiterate that the apparent contradiction between the theoretical illconditioning and numerical stability is due to the spectral cutoff present in any practical reconstruction procedure.
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 Ω = (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 α ∈ (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 α = 1 4 and α = 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.
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 α = 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 α < 1. However, for α = 1, the singular values decay exponentially, and the decay rate   T 1 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.

Sideways fractional diffusion
This representation is well known for the case α = 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

where the convolution kernel α R s ( ) is given by a Wright function in the form
In case of α = 1, i.e., classical diffusion, the kernel R(s) is given explicitly by x This problem is also known as the lateral Cauchy problem in the literature. In the case α = 1, it is known that the inverse problem is severely ill-posed [ .

xx x
The general solution  u is given by and thus the solution The solution h(t) can then be recovered by an inverse Laplace transform 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 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 Inverse Problems 31 (2015) 035003 B Jin and W Rundell the data perturbation are amplified by an exponentially growing multiplier α e z 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 illposed. 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 α = 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 α = 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 α = 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.
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., α = 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]. 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 onedimensional 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- Hence the measured data = g u T ( ) is given by By taking inner product with ϕ j on both sides, we arrive at the following representation of the source term f in terms of the measured data g  ( ) 1 is given by = u x t ( , ) 1 and ≡ f 0, but another is π = u x t t T ( , ) cos (2 ) and π π = − f T tT ( 2 ) sin (2 ). Likewise, in the fractional case, we can take π = u tT cos (2 ) 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 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

By taking inner product with ϕ j on both sides, we deduce
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 ( ) might be recovered. In other words, due to a slower local decay of the exponential function λ − e t , compared with the Mittag-Leffler function , 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.
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 t : (0, ) x is related to the unknown p(t) by

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
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 T [0, ], 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.
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   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 , 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., α = 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.

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
To illustrate the point, we present in figure 10 some numerical results of reconstructing a discontinuous potential (with χ 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 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.
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 α ∈ (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 [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.

Fractional derivative as an inverse solution
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 and vice-versa. The main result in [52] is that the operator  has a unique fixed point on  and indeed  is monotone in the sense of preserving the partial order on , i.e., if ⩾ a a 1 2 then ⩽   a a 1 2 . Given these developments, it might seem that a parallel construction for the time fractional diffusion counterpart, ∂ = α u a t u ( ) t 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.

Inverse problems for space fractional diffusion
Now we turn to differential equations involving a fractional derivative in space. There are several possible choices of a fractional derivative in space, e.g., Djrbashian-Caputo fractional derivative, Riemann-Liouville fractional derivative, Riesz derivative, and fractional Laplacian [3], which all have received considerable attention. In recent years, the use of the fractional Laplacian is especially popular in high-dimensional spaces, and admits a welldeveloped analytical theory. We shall focus on the left-sided Djrbashian-Caputo fractional derivative ∈ (1, 2), and the one-dimensional case, and consider the following four inverse problems: inverse Sturm-Liouville problem, Cauchy problem for a fractional elliptic equation, backwards diffusion, and sideways problem.  ∈ (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]

Inverse Problems 31 (2015) 035003 B Jin and W Rundell
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 β = 2, the Dirichlet spectrum cannot uniquely determine a general potential q. Theoretically, the surprising uniqueness in the fractional case remains to be established. Naturally, one can also consider the Riemann-Liouville case:  Hence, for any β ∈ (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 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.

Cauchy problem for fractional elliptic equation
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 [53, p 46], we deduce that the solution ψ j to the fractional ordinary differential equation is given by . Nonetheless, according to the discussions in section 4.1, the eigenvalues λ { } 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.
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

Inverse Problems 31 (2015) 035003 B Jin and W Rundell
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.

Sideways problem
Last we return to the classical sideways diffusion problem but now with a fractional derivative in space rather than in time. Let Ω = (0, 1) be the unit interval. Then the onedimensional space fractional diffusion equation is given by 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.

Concluding remarks
Anomalous diffusion processes arise in many disciplines, and the physics behind is very different from normal diffusion. The unusual physics greatly influences the behavior of related forward problems. Further, it is well known that backward fractional diffusion is much less ill-posed than the classical backward diffusion, which has contributed to the belief that inverse problems for anomalous diffusion are always better behaved than that for the normal diffusion. In this work we have examined several exemplary inverse problems for anomalous diffusion processes in a numerical and semi-analytical manner. These include the sideways problem, backward problem, inverse source problem, inverse Sturm-Liouville problem and Cauchy problems. Our findings indicate that anomalous diffusion can give rise to very unusual new features, but they only partially confirm the belief: depending on the data and unknown, it may influence either positively or negatively the degree of ill-posedness of the inverse problem.
The mathematical study of inverse problems in anomalous diffusion is still in its infancy. There are only a few rigorous theoretical results on the uniqueness, existence and stability, which mostly focus on the one-dimensional case, and there are many more open problems 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 ( ) ( 1) 1 1 . We note that an algorithm for the Wright function ρ μ W z ( ) , over the whole complex plane  with rigorous error analysis is still missing. The endeavor in this direction would almost certainly involve dividing the complex domain  into different regions, and using different approximations on each region separately.
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]. . 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