Regularization methods for unknown source in ... - Semantic Scholar
Available online at www.sciencedirect.com
Mathematics and Computers in Simulation 85 (2012) 45–56
Regularization methods for unknown source in space fractional diffusion equation WenYi Tian a , Can Li a,b , Weihua Deng a,∗ , Yujiang Wu a b
a School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, People’s Republic of China Department of Applied Mathematics, School of Sciences, Xi’san University of Technology, Xi’san, Shaanxi 710054, People’s Republic of China
Received 5 June 2012; accepted 26 August 2012 Available online 13 October 2012
Abstract We discuss determining the unknown steady source in a space fractional diffusion equation and show that both the Fourier and wavelet dual least squares regularization methods work well for the ill-posed problem. The detailed error estimates are also strictly established for both of the methods. Moreover, the algorithm implementation and the corresponding numerical results are presented for the Fourier regularization method. © 2012 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Space fractional diffusion equation; Ill-posed problems; Regularization methods
1. Introduction During the past decades, the interest in fractional calculus and fractional differential equations has rapidly increased among engineers and scientists due to their vast potential of applications, including physical, chemical, mechanical engineering, signal processing and systems identification, biology systems, control theory, and finance etc. For a review the reader can refer to [16,17,21,27] and references therein. In fact, it can be noticed that one of the most successful and concrete applications of fractional calculus and fractional differential equations is to effectively characterize the anomalous diffusion. It is well known that the ordinary diffusion process is intimately related to the validity of the central limit theorem, which is characterized by the linear dependence of the mean square displacement x2 (t) ∼ κt with the diffusion coefficient κ. However, some diffusion processes, especially in various complex systems, no longer follow Gaussian behavior. This phenomenon is named anomalous diffusion which is described by the nonlinear growth of the mean square displacement x(t) of a diffusion particle over time t: x2 (t) ∼ κα tα , where κα is the diffusion coefficient, and α is the anomalous diffusion exponent. For different α, the anomalous diffusion is classified into subdiffusion (0 < α < 1), normal diffusion (α = 1), superdiffusion (α > 1), and ballistic diffusion (α = 2) [16,17]. And the Fick’s law is inevitable to be modified in order to precisely describe the anomalous diffusion behavior [14].
∗
Corresponding author. E-mail addresses: [email protected], [email protected] (W. Deng).
0378-4754/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.matcom.2012.08.011
46
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
The space fractional diffusion equation (will be given in the next section) we consider in this paper can describe the probability distribution of the particles having superdiffusion. Much progress has been made for numerically solving space fractional partial differential equations [5,9,12,14,15]. Here instead of further pursuing research in this direction, we discuss the space fractional inverse diffusion equation, i.e., to determine an unknown source, which depends only on the spatial variable, in the one dimensional space fractional diffusion equation. Determination of unknown source is to obtain the information about a physical object or system by observed datum, and it is one of the most important and well-studied problems in many branches of engineering sciences. Identifying unknown source is an inverse and severely ill-posed problem. We use Fourier and wavelet dual least squares regularization methods to solve the ill-posed problem. The error estimates are presented and for Fourier regularization method the numerical experiments are also performed. In fact, there are some works on the inverse problem of time fractional diffusion equations [2,3,13,18]. However, there is limited work on space fractional diffusion equation [28], which use Tikhonov regularization algorithm. The outline of the work is as follows. In Section 2 we briefly introduce the mathematical model. Section 3 presents the Fourier method to the Cauchy problem of one dimensional space fractional diffusion equation and establishes the error estimates. The wavelet dual least squares regularization method and its error estimates are given in Section 4. The numerical experiments are performed in Section 5, and we conclude this paper with some discussions in Section 6. 2. Mathematical model In this section, we briefly present the physical background of space fractional diffusion equation and introduce its concrete form. Following [14], denote u(x, t) as the probability distribution of the particles (or concentration of solute) at point x and time t. For the pure diffusion process, the conservation of mass equation can be expressed as ∂u(x, t) ∂J(x, t) =− + f (x), ∂t ∂x where J is the mass flux and f(x) a stable source. Usually, the mass flux J refers to the Fick’s first law ∂u(x, t) . ∂x However, the anomalous diffusion is given by the generalized Fick’s law
(2.1)
J(x, t) = −κ
(2.2)
α−1 }u, J(x, t) = −κα {p−∞ Dxα−1 − qx D∞
(2.3)
α are left and right where κα > 0, and p + q = 1 for p, q ≥ 0 (the case in [14] is p = q = 1/2), −∞ Dxα and x D∞ Riemann–Liouville fractional derivatives of order α (1 < α < 2), respectively, defined by x 1 ∂2 α (x − ζ)1−α u(ζ, t)dζ, −∞ Dx u(x, t) = (2 − α) ∂x2 −∞
and α x D∞ u(x, t)
(−1)α d2 = (2 − α) dx2
x
∞
(ζ − x)1−α u(ζ, t)dζ,
α is the smallest integer no less than α, and (·) denotes the Gamma function. Substituting (2.3) into (2.1) leads to ∂u α = κα ∇p,q u + f (x), (2.4) ∂t α = {p − ∞Dα + q xDα }, still being meaningful when α tends to 2, and corresponding to second order where ∇p,q x ∞ derivative. The existence and uniqueness of a variational solution to (2.4) with boundary and initial conditions is discussed in [25]. Our goal in this paper is to determine an unknown source f(x) in the space fractional diffusion equation (2.4) with the following conditions u(x, 0) = 0,
u(x, T ) = g(x),
(2.5)
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
47
that is to identify the source f(x) from the data g(x). We assume that the solution u(x, t) and f(x) exist for the exact data g(x), moreover f(x) is bounded in H s (R), i.e., there exists a constant M such that f s ≤ M,
(2.6)
where · s is the norm of Sobolev space 1 f s = ( (1 + |ξ|s )2 |fˆ (ξ)|2 dξ) 2 .
H s (R) (s
≥ 0), and it is defined by (2.7)
R
From (2.4), (2.5), and (2.6), we can define the operator K on L2 such that Kf (x) = g(x).
(2.8)
In fact, it can be easily shown that K is an operator from Hs to Hs . Remark 1. If we take p = q = α ∇p,q
−1 2 cos(πα/2)
in (2.4), then we obtain that the one dimensional fractional derivative α
is equivalent to the fractional Laplacian operator −(− ) 2 or the fractional Riesz derivative operator operator (more details one can see [10,26]), i.e., α −1 [−∞ Dxα + x D∞ α ]u(x, t) = −(− ) 2 u := −F−1 (|ξ|α F(u(ξ, t))). 2 cos(πα/2)
(2.9)
3. Fourier transform method Based on the Fourier transform method, there are some discussions [20,16,17] on the analytical solution of (2.4) with different initial and boundary conditions. There are also some results [1,7,8] on the inverse problem of the special case of (2.4), namely α → 2, which corresponds to the sideways heat equation. Here we further discuss the fractional case of (2.4), i.e., 1 < α < 2. If u( · , t) ∈ L2 (R), then the Fourier transform operator F : L2 (R) → L2 (R) about x is given by ∞ 1 u(ξ, ˆ t) = F(u(x, t)) = √ (3.1) u(x, t)e−iξx dx, ξ ∈ R, 2π −∞ α are obtained as [20,26] and the Fourier transforms of −∞ Dxα and x D∞
F{−∞ Dxα u(x, t)} = (iξ)α u(ξ, ˆ t), ξ ∈ R, and α F{x D∞ u(x, t)} = (−iξ)α u(ξ, ˆ t), ξ ∈ R.
Then we have α F{∇p,q u(x, t)} = (p(iξ)α + q(−iξ)α )u(ξ, ˆ t), ξ ∈ R.
(3.2)
According to the definition of complex power function, the principal value of the Fourier coefficient in (3.2) is πα πα p(iξ)α + q(−iξ)α = [cos( ) + i(p − q) sin(± )]|ξ|α . (3.3) 2 2 Applying the Fourier transform with respect to x in both sides of (2.4) and using (3.2) and (3.3), we have the following equation in Fourier space ⎧ ˆ t) ⎨ ∂u(ξ, = (ξ, α)u(ξ, ˆ t) + fˆ (ξ), ξ ∈ R, 0 < t < T, ∂t (3.4) ⎩ u(ξ, ˆ 0) = 0, u(ξ, ˆ T ) = gˆ (ξ), ξ ∈ R, πα α α where (ξ, α) = κα cos( πα 2 )|ξ| + i[κα (p − q) sin(± 2 )|ξ| ]. It can be easily verified that the solution of (3.4) is
u(ξ, ˆ t) =
(e (ξ,α)t − 1)fˆ (ξ) .
(ξ, α)
(3.5)
48
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
Since u(ξ, ˆ T ) = gˆ (ξ), we get gˆ (ξ) =
(e (ξ,α)T − 1)fˆ (ξ) .
(ξ, α)
(3.6)
: H s (R) → H s (R) is a multiplication operator, From (3.6) we conclude that K gˆ (ξ) = Kf (x) =
e (ξ,α)T − 1 ˆ fˆ (ξ), f (ξ) = K
(ξ, α)
(3.7)
is self-adjoint. By an elementary calculation, we have and K
(ξ, α)ˆg(ξ) fˆ (ξ) = (ξ,α)T . e −1
(3.8)
In the practical applications, the data g(x) should be replaced by the measured one gδ (x), and it satisfies gδ (·) − g(·) ≤ δ, where · denotes L2 (R) norm. Then (3.8) becomes as δ δ (ξ) = (ξ, α)g (ξ) . f (3.9) e (ξ,α)T − 1 Obviously, the noisy data (gδ (ξ) − g(ξ)) is magnified by O(|ξ|α ) in the high frequency component of f (ξ), so the problem is ill-posed. For solving this problem, a simple way is to apply Fourier regularization method which has been used for handling ill-posed problem with classical derivatives [7,8]. The main idea of the Fourier regularization method is to eliminate all the high frequencies in the Fourier space and consider (3.9) only for |ξ| ≤ ξ max . Hence we obtain a regularized solution ∞ 1
(ξ, α)gδ (ξ) δ fξmax (x) = √ eiξx (ξ,α)T (3.10) χmax dξ, e −1 2π −∞ where χmax denotes a characteristic function of the interval [− ξ max , ξ max ], the choice of the optimal regularization parameter ξ max will be given in the sequel. Lemma 2. For 1 < α < 2, there exists |
κα |ξ|α
(ξ, α) , |≤ πα α −1 1 − e(κα cos( 2 )|ξ| T )
(3.11)
e (ξ,α)T
πα α α where (ξ, α) = κα cos( πα 2 )|ξ| + i[κα (p − q) sin(± 2 )|ξ| ].
Proof. For 1 < α < 2, and p + q = 1 for p, q ≥ 0, it is evident that πα πα | (ξ, α)| ≤ (κα cos( )|ξ|α )2 + (κα (p − q) sin(± )|ξ|α )2 2 2 ≤ κα |ξ|α . With the fact |ez | = eRe(z) and cos( πα 2 ) < 0, we have |
(ξ, α) | (ξ, α)| | (ξ, α)| | ≤ (ξ,α)T = Re( (ξ,α)T ) ||e |e −1 | − 1| − 1| κα |ξ|α κα |ξ|α ≤ = πα α πα α . )|ξ| T ) )|ξ| T ) (κα cos( (κα cos( 2 2 − 1| |e 1−e
e (ξ,α)T
1
Theorem 3. Under the assumptions (2.6) and gδ (·) − g(·) ≤ δ, if ξmax = ( Mδ ) α+s , then the following error estimate holds α
s
f − fξδmax ≤ CM α+s δ α+s (1 + o(1)),
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
where C = (1 + 2κα ) is a positive constant independent of M and δ. Proof. According to Parseval’s Relation, and using (3.8) and (3.10), we obtain δ δ ˆ ˆ ˆ f − fξδmax = fˆ − f ξmax ≤ f − f χmax + f χmax − fξmax
(ξ, α)
(ξ, α)
(ξ, α) g − (ξ,α)T gχmax + (ξ,α)T ( g − gδ )χmax −1 e −1 e −1 1 1
(ξ, α)
(ξ, α) =( | (ξ,α)T | (ξ,α)T g|2 dξ) 2 + ( ( g − gδ )|2 dξ) 2 . −1 −1 |ξ|>ξmax e |ξ|≤ξmax e
≤
e (ξ,α)T
Due to (2.6), we have
1
(ξ, α) 2 2 | (ξ,α)T g| dξ) ( −1 |ξ|>ξmax e
=(
=( ≤
When |ξ| ≥ |
1 1/α , (κα | cos( πα 2 )|T )
|ξ|>ξmax
1 2 2 |f (ξ)| dξ)
|ξ|>ξmax −s M. ξmax
s −2
(1 + |ξ| )
1 (1 + |ξ| ) |f (ξ)| dξ) 2 s 2
it obtains from Lemma 2 that
κα |ξ|α κα |ξ|α = πα α πα (κα cos( −(κα | cos( )|ξ| T ) )||ξ|α T ) 2 2 1−eα 1−e κα |ξ| ≤ ≤ 2κα |ξ|α ; 1 − e−1
(ξ, α) | ≤
(ξ,α)T e −1
and for |ξ|
a. We first get the datum g = {g(xk )}N k=1 representing values of g(x) on an equidistant grid, and then add a random uniformly distributed perturbation over (− 1, 1) to each data of g to obtain perturbation datum gδ = g + λrand(size(g)),
(5.1)
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
where λ is the noise level, and xk = (k − 1)h − a, h =
N 1
δ (gkδ − gk )2 . δ = g − g l2 = N
2a N−1 .
53
So the noise δ is measured by (5.2)
k=1
Usually it is difficult to get the analytic solution of the fractional partial differential equation (2.4) with some boundary conditions, so we have to obtain the datum g by numerical methods for a given f(x) at first. The numerical procedure of the Fourier regularization method is described as follows: Step 1. Obtaining the datum g by numerical method; Step 2. Calculating the vector gδ by Fast Fourier Transform (FFT); Step 3. Computing the vector N −1 2 Θ(ξ, α) )χ(ξ ) , gδ(ξ k k eΘ(ξ,α)T − 1 k=− N
(5.3)
2
with ξ k = πk/a and χ(ξ k ) = 1 for |ξ k | ≤ ξ max or else χ(ξ k ) = 0; Step 4. Acting the Inverse Fourier Transform (IFFT) on (5.3), and producing the vector {fξδmax (xi )}N i=1 . The accuracy of regularization result is given in the l2 norm defined by
N 1
δ Ef = f − fξmax l2 = (fk − (fξδmax )k )2 . N k=1
For the difficulty to find the analytic solution of (2.4), here we take finite difference method to get the datum g for a given f(x) [14,15]. And we discretize the Eq. (2.4) by a fractional Crank–Nicholson scheme within a finite interval [− a, a] as follows, un+1 − uni i = τ
i+1 N−i+1
1 1 α n+1 ωkα un+1 κα α [p ωk ui−k+1 + q i+k−1 ] 2 h k=0
k=0
i+1 N−i+1
1 1
ωkα uni+k−1 ] + κα α [p ωkα uni−k+1 + q 2 h k=0
k=0
(5.4)
1 1 + fin+1 + fin + O(τ 2 + h), 2 2 where uni = u(xi , tn ), τ and h are time and space step size respectively. The left and right Riemann–Liouville fractional derivatives in (2.4) can be respectively described as ([14,15]) x+a ]+1 h
1 α ωkα u(x − (k − 1)h, t) + O(h), −a Dx u(x, t) = α h k=0 a−x [ ]+1 h
1 α = α ωkα u(x − (k − 1)h, t) + O(h), x Da u(x, t) h k=0 α where ωkα = (−1)k . k [
54
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56 1.5
1.5 Numerical Exact
1
1.5 Numerical Exact
1
0.5
0.5
0.5
0
0
0
−0.5
−0.5
−0.5
−1
−1
−1
−1.5 −10
−5
0
5
10
−1.5 −10
−5
0
5
Numerical Exact
1
10
−1.5 −10
−5
0
5
10
Fig. 1. Comparison of the exact solution f with the regularization solution fξδmax with different noise for α = 1.1. 1.5
1.5 Numerical Exact
1
1.5 Numerical Exact
1
0.5
0.5
0.5
0
0
0
−0.5
−0.5
−0.5
−1
−1
−1
−1.5 −10
−5
0
5
10
−1.5 −10
−5
0
5
Numerical Exact
1
10
−1.5 −10
−5
0
5
10
Fig. 2. Comparison of the exact solution f with the regularization solution fξδmax with different noise for α = 1.5.
x2
−4 Example. Take f (x) = ( x4 − 3x and T = 1 in (2.4). It is in fact that f(x) approaches to 0 as |x| > 6, then the 2 ) e numerical experiment is implemented in the interval [− 10, 10]. In the numerical computation, we choose κα = 1, p = q = 0.5, s = 2, and N = 64. And f 2 can be computed by Maple, 1 i.e., M = 4.452844182. From Theorem 3, we can take ξmax = ( Mδ ) α+2 . Fig. 1–3 shows the pictures of the exact solution f and the regularization solution fξδmax for α = 1.1, 1.5, 1.9, respectively. And Table 1 lists the regularization parameter ξ max and the error Ef for different noise level λ, the error decreases when λ becomes small. From the numerical experiments, it is noted that the error decreases with δ becoming small for α ∈ (1, 2). 3
1.5
1.5 Numerical Exact
1
1.5 Numerical Exact
1
0.5
0.5
0.5
0
0
0
−0.5
−0.5
−0.5
−1
−1
−1
−1.5 −10
−5
0
5
10
−1.5 −10
−5
0
5
Numerical Exact
1
10
−1.5 −10
−5
0
Fig. 3. Comparison of the exact solution f with the regularization solution fξδmax with different noise for α = 1.9.
5
10
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
55
Table 1 The regularization parameter ξ max and the error Ef for different noise level λ. λ
10−1
10−2
10−3
10−4
10−5
α = 1.1
ξ max Ef
4.0675 0.0406
8.6406 0.0398
17.9457 0.0397
37.3143 0.0397
79.0675 0.0397
α = 1.5
ξ max Ef
3.4583 0.0337
6.6125 0.0198
12.7839 0.0106
24.4357 0.0102
48.6976 0.0102
α = 1.9
ξ max Ef
3.1491 0.0375
5.4188 0.0197
9.7296 0.0127
17.9282 0.0099
32.0046 0.0098
6. Conclusions We have extended two regularization methods, being ever used to regularize the classical pure diffusion equation, to solve the inverse space fractional diffusion equation which is ill-posed. Rigorous error estimates are established for both of the methods, and extensive numerical experiments are performed for Fourier regularization method. We hope this work will be interesting and/or helpful for further discussing the inverse problem of space fractional partial differential equations. Acknowledgements This work was supported by the Program for New Century Excellent Talents in University under Grant No. NCET-090438, the National Natural Science Foundation of China under Grant No. 10801067 and No. 11271173, the Fundamental Research Funds for the Central Universities under Grant No. lzujbky-2010-63 and No. lzujbky-2012-k26, the National Basic Research Program of China 973 Program No. 2011CB706900, and the Starting Research Fund from the Xi’an University of Technology under Grant No. 108-211206. References [1] F. Berntsson, A spectral method for solving the sideways heat equation, Inverse Problems 15 (1999) 891–906. [2] A.N. Bondarenko, D.S. Ivaschenko, Numerical methods for solving inverse problems for time fractional diffusion equation with variable coefficient, Journal of Inverse and III-posed Problems 17 (2009) 419–440. [3] J. Cheng, J. Nakagawa, M. Yamamoto, T. Yamazaki, Uniqueness in an inverse problem for a one-dimensional fractional diffusion equation, Inverse Problems 25 (2009) 115002. [14] A.S. Chaves, A fractional diffusion equation to describe Lévy flights, Physics Letters A 239 (1998) 13–16. [5] W.H. Deng, Finite element method for the space and time fractional Fokker–Planck equation, SIAM Journal on Numerical Analysis 47 (2008) 204–226. [6] F.F. Dou, C.L. Fu, Determining an unknown source in the heat equation by a wavelet dual least squares method, Applied Mathematics Letters 22 (2009) 661–667. [7] F.F. Dou, C.L. Fu, F.-L. Yang, Optimal error bound and Fourier regularization for identifying an unknown source in the heat equation, Journal of Computational and Applied Mathematics 230 (2009) 728–737. [8] L. Eldén, F. Berntsson, T. Regi´nska, Wavelet and Fourier method for solving the sideways heat equation, SIAM Journal on Scientific Computing 21 (2000) 2187–2205. [9] V.J. Ervin, N. Heuer, J.P. Roop, Numerical approximation of a time dependent, nonlinear, space-fractional diffusion equation, SIAM Journal on Numerical Analysis 45 (2007) 572–591. [10] H.C. Fogedby, Lévy flights in random environments, Physical Review Letters 73 (1994) 2517–2520. [11] C.L. Fu, C.Y. Qiu, Wavelet and error estimation of surface heat flux, Journal of Computational and Applied Mathematics 150 (2003) 143–155. [12] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker–Planck equation, Journal of Computational and Applied Mathematics 166 (2004) 209–219. [13] J.J. Liu, M. Yamamoto, A backward problem for the time-fractional diffusion equation, Applicable Analysis 89 (2010) 1769–1788. [14] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection-dispersion flow equations, Journal of Computational and Applied Mathematics 172 (2003) 65–77. [15] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Applied Numerical Mathematics 56 (2006) 80–90. [16] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Physics Report 339 (2000) 1–77.
56
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
[17] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, Journal of Physics A: Mathematical and General 37 (2004) R161–R208. [18] D.A. Murio, Stable numerical solution of a fractional-diffusion inverse heat conduction problem, Computers and Mathematics with Applications 53 (2007) 1492–1501. [19] I. Daubechies, Ten lectures on wavelets, in: CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, 1992. [20] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [21] I. Podlubny, Geometric and Physical interpretation of fractional integration and fractional differentiation, Fractional Calculus and Applied Analisys 5 (2002) 367–386. [22] A.L. Qian, A new wavelet method for solving an ill-posed problem, Applied Mathematics and Computation 203 (2008) 635–640. [23] T. Regi´nska, Sideways heat equation and wavelets, Journal of Computational and Applied Mathematics 63 (1995) 209–214. [24] T. Regi´nska, L. Eldén, Solving the sideways heat equation by a wavelet-Galerkin method, Inverse Problems 13 (1997) 1093–1106. [25] J.P. Roop, Variational Solution of the Fractional Advection-Dispersion Equation, The Graduate School of Clemson University, 2004. [26] S. Samko, A. Kilbas, O. Marichev, Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach, London, 1993. [27] E. Scalas, R. Gorenflo, F. Mainardi, Fractional calculus and continuous-time finance, Physica A 284 (2000) 376–384. [28] H. Wei, W. Chen, H.G. Sun, X.C. Li, A coupled method for inverse source problem of spatial fractional anomalous diffusion equations, Inverse Problems in Science and Engineering 18 (2010) 945–956. [29] X.T. Xiong, C.L. Fu, Determining surface temperature and heat flux by a wavelet dual least squares method, Journal of Computational and Applied Mathematics 201 (2007) 189–207.
Mathematics and Computers in Simulation 85 (2012) 45–56
Regularization methods for unknown source in space fractional diffusion equation WenYi Tian a , Can Li a,b , Weihua Deng a,∗ , Yujiang Wu a b
a School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, People’s Republic of China Department of Applied Mathematics, School of Sciences, Xi’san University of Technology, Xi’san, Shaanxi 710054, People’s Republic of China
Received 5 June 2012; accepted 26 August 2012 Available online 13 October 2012
Abstract We discuss determining the unknown steady source in a space fractional diffusion equation and show that both the Fourier and wavelet dual least squares regularization methods work well for the ill-posed problem. The detailed error estimates are also strictly established for both of the methods. Moreover, the algorithm implementation and the corresponding numerical results are presented for the Fourier regularization method. © 2012 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Space fractional diffusion equation; Ill-posed problems; Regularization methods
1. Introduction During the past decades, the interest in fractional calculus and fractional differential equations has rapidly increased among engineers and scientists due to their vast potential of applications, including physical, chemical, mechanical engineering, signal processing and systems identification, biology systems, control theory, and finance etc. For a review the reader can refer to [16,17,21,27] and references therein. In fact, it can be noticed that one of the most successful and concrete applications of fractional calculus and fractional differential equations is to effectively characterize the anomalous diffusion. It is well known that the ordinary diffusion process is intimately related to the validity of the central limit theorem, which is characterized by the linear dependence of the mean square displacement x2 (t) ∼ κt with the diffusion coefficient κ. However, some diffusion processes, especially in various complex systems, no longer follow Gaussian behavior. This phenomenon is named anomalous diffusion which is described by the nonlinear growth of the mean square displacement x(t) of a diffusion particle over time t: x2 (t) ∼ κα tα , where κα is the diffusion coefficient, and α is the anomalous diffusion exponent. For different α, the anomalous diffusion is classified into subdiffusion (0 < α < 1), normal diffusion (α = 1), superdiffusion (α > 1), and ballistic diffusion (α = 2) [16,17]. And the Fick’s law is inevitable to be modified in order to precisely describe the anomalous diffusion behavior [14].
∗
Corresponding author. E-mail addresses: [email protected], [email protected] (W. Deng).
0378-4754/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.matcom.2012.08.011
46
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
The space fractional diffusion equation (will be given in the next section) we consider in this paper can describe the probability distribution of the particles having superdiffusion. Much progress has been made for numerically solving space fractional partial differential equations [5,9,12,14,15]. Here instead of further pursuing research in this direction, we discuss the space fractional inverse diffusion equation, i.e., to determine an unknown source, which depends only on the spatial variable, in the one dimensional space fractional diffusion equation. Determination of unknown source is to obtain the information about a physical object or system by observed datum, and it is one of the most important and well-studied problems in many branches of engineering sciences. Identifying unknown source is an inverse and severely ill-posed problem. We use Fourier and wavelet dual least squares regularization methods to solve the ill-posed problem. The error estimates are presented and for Fourier regularization method the numerical experiments are also performed. In fact, there are some works on the inverse problem of time fractional diffusion equations [2,3,13,18]. However, there is limited work on space fractional diffusion equation [28], which use Tikhonov regularization algorithm. The outline of the work is as follows. In Section 2 we briefly introduce the mathematical model. Section 3 presents the Fourier method to the Cauchy problem of one dimensional space fractional diffusion equation and establishes the error estimates. The wavelet dual least squares regularization method and its error estimates are given in Section 4. The numerical experiments are performed in Section 5, and we conclude this paper with some discussions in Section 6. 2. Mathematical model In this section, we briefly present the physical background of space fractional diffusion equation and introduce its concrete form. Following [14], denote u(x, t) as the probability distribution of the particles (or concentration of solute) at point x and time t. For the pure diffusion process, the conservation of mass equation can be expressed as ∂u(x, t) ∂J(x, t) =− + f (x), ∂t ∂x where J is the mass flux and f(x) a stable source. Usually, the mass flux J refers to the Fick’s first law ∂u(x, t) . ∂x However, the anomalous diffusion is given by the generalized Fick’s law
(2.1)
J(x, t) = −κ
(2.2)
α−1 }u, J(x, t) = −κα {p−∞ Dxα−1 − qx D∞
(2.3)
α are left and right where κα > 0, and p + q = 1 for p, q ≥ 0 (the case in [14] is p = q = 1/2), −∞ Dxα and x D∞ Riemann–Liouville fractional derivatives of order α (1 < α < 2), respectively, defined by x 1 ∂2 α (x − ζ)1−α u(ζ, t)dζ, −∞ Dx u(x, t) = (2 − α) ∂x2 −∞
and α x D∞ u(x, t)
(−1)α d2 = (2 − α) dx2
x
∞
(ζ − x)1−α u(ζ, t)dζ,
α is the smallest integer no less than α, and (·) denotes the Gamma function. Substituting (2.3) into (2.1) leads to ∂u α = κα ∇p,q u + f (x), (2.4) ∂t α = {p − ∞Dα + q xDα }, still being meaningful when α tends to 2, and corresponding to second order where ∇p,q x ∞ derivative. The existence and uniqueness of a variational solution to (2.4) with boundary and initial conditions is discussed in [25]. Our goal in this paper is to determine an unknown source f(x) in the space fractional diffusion equation (2.4) with the following conditions u(x, 0) = 0,
u(x, T ) = g(x),
(2.5)
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
47
that is to identify the source f(x) from the data g(x). We assume that the solution u(x, t) and f(x) exist for the exact data g(x), moreover f(x) is bounded in H s (R), i.e., there exists a constant M such that f s ≤ M,
(2.6)
where · s is the norm of Sobolev space 1 f s = ( (1 + |ξ|s )2 |fˆ (ξ)|2 dξ) 2 .
H s (R) (s
≥ 0), and it is defined by (2.7)
R
From (2.4), (2.5), and (2.6), we can define the operator K on L2 such that Kf (x) = g(x).
(2.8)
In fact, it can be easily shown that K is an operator from Hs to Hs . Remark 1. If we take p = q = α ∇p,q
−1 2 cos(πα/2)
in (2.4), then we obtain that the one dimensional fractional derivative α
is equivalent to the fractional Laplacian operator −(− ) 2 or the fractional Riesz derivative operator operator (more details one can see [10,26]), i.e., α −1 [−∞ Dxα + x D∞ α ]u(x, t) = −(− ) 2 u := −F−1 (|ξ|α F(u(ξ, t))). 2 cos(πα/2)
(2.9)
3. Fourier transform method Based on the Fourier transform method, there are some discussions [20,16,17] on the analytical solution of (2.4) with different initial and boundary conditions. There are also some results [1,7,8] on the inverse problem of the special case of (2.4), namely α → 2, which corresponds to the sideways heat equation. Here we further discuss the fractional case of (2.4), i.e., 1 < α < 2. If u( · , t) ∈ L2 (R), then the Fourier transform operator F : L2 (R) → L2 (R) about x is given by ∞ 1 u(ξ, ˆ t) = F(u(x, t)) = √ (3.1) u(x, t)e−iξx dx, ξ ∈ R, 2π −∞ α are obtained as [20,26] and the Fourier transforms of −∞ Dxα and x D∞
F{−∞ Dxα u(x, t)} = (iξ)α u(ξ, ˆ t), ξ ∈ R, and α F{x D∞ u(x, t)} = (−iξ)α u(ξ, ˆ t), ξ ∈ R.
Then we have α F{∇p,q u(x, t)} = (p(iξ)α + q(−iξ)α )u(ξ, ˆ t), ξ ∈ R.
(3.2)
According to the definition of complex power function, the principal value of the Fourier coefficient in (3.2) is πα πα p(iξ)α + q(−iξ)α = [cos( ) + i(p − q) sin(± )]|ξ|α . (3.3) 2 2 Applying the Fourier transform with respect to x in both sides of (2.4) and using (3.2) and (3.3), we have the following equation in Fourier space ⎧ ˆ t) ⎨ ∂u(ξ, = (ξ, α)u(ξ, ˆ t) + fˆ (ξ), ξ ∈ R, 0 < t < T, ∂t (3.4) ⎩ u(ξ, ˆ 0) = 0, u(ξ, ˆ T ) = gˆ (ξ), ξ ∈ R, πα α α where (ξ, α) = κα cos( πα 2 )|ξ| + i[κα (p − q) sin(± 2 )|ξ| ]. It can be easily verified that the solution of (3.4) is
u(ξ, ˆ t) =
(e (ξ,α)t − 1)fˆ (ξ) .
(ξ, α)
(3.5)
48
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
Since u(ξ, ˆ T ) = gˆ (ξ), we get gˆ (ξ) =
(e (ξ,α)T − 1)fˆ (ξ) .
(ξ, α)
(3.6)
: H s (R) → H s (R) is a multiplication operator, From (3.6) we conclude that K gˆ (ξ) = Kf (x) =
e (ξ,α)T − 1 ˆ fˆ (ξ), f (ξ) = K
(ξ, α)
(3.7)
is self-adjoint. By an elementary calculation, we have and K
(ξ, α)ˆg(ξ) fˆ (ξ) = (ξ,α)T . e −1
(3.8)
In the practical applications, the data g(x) should be replaced by the measured one gδ (x), and it satisfies gδ (·) − g(·) ≤ δ, where · denotes L2 (R) norm. Then (3.8) becomes as δ δ (ξ) = (ξ, α)g (ξ) . f (3.9) e (ξ,α)T − 1 Obviously, the noisy data (gδ (ξ) − g(ξ)) is magnified by O(|ξ|α ) in the high frequency component of f (ξ), so the problem is ill-posed. For solving this problem, a simple way is to apply Fourier regularization method which has been used for handling ill-posed problem with classical derivatives [7,8]. The main idea of the Fourier regularization method is to eliminate all the high frequencies in the Fourier space and consider (3.9) only for |ξ| ≤ ξ max . Hence we obtain a regularized solution ∞ 1
(ξ, α)gδ (ξ) δ fξmax (x) = √ eiξx (ξ,α)T (3.10) χmax dξ, e −1 2π −∞ where χmax denotes a characteristic function of the interval [− ξ max , ξ max ], the choice of the optimal regularization parameter ξ max will be given in the sequel. Lemma 2. For 1 < α < 2, there exists |
κα |ξ|α
(ξ, α) , |≤ πα α −1 1 − e(κα cos( 2 )|ξ| T )
(3.11)
e (ξ,α)T
πα α α where (ξ, α) = κα cos( πα 2 )|ξ| + i[κα (p − q) sin(± 2 )|ξ| ].
Proof. For 1 < α < 2, and p + q = 1 for p, q ≥ 0, it is evident that πα πα | (ξ, α)| ≤ (κα cos( )|ξ|α )2 + (κα (p − q) sin(± )|ξ|α )2 2 2 ≤ κα |ξ|α . With the fact |ez | = eRe(z) and cos( πα 2 ) < 0, we have |
(ξ, α) | (ξ, α)| | (ξ, α)| | ≤ (ξ,α)T = Re( (ξ,α)T ) ||e |e −1 | − 1| − 1| κα |ξ|α κα |ξ|α ≤ = πα α πα α . )|ξ| T ) )|ξ| T ) (κα cos( (κα cos( 2 2 − 1| |e 1−e
e (ξ,α)T
1
Theorem 3. Under the assumptions (2.6) and gδ (·) − g(·) ≤ δ, if ξmax = ( Mδ ) α+s , then the following error estimate holds α
s
f − fξδmax ≤ CM α+s δ α+s (1 + o(1)),
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
where C = (1 + 2κα ) is a positive constant independent of M and δ. Proof. According to Parseval’s Relation, and using (3.8) and (3.10), we obtain δ δ ˆ ˆ ˆ f − fξδmax = fˆ − f ξmax ≤ f − f χmax + f χmax − fξmax
(ξ, α)
(ξ, α)
(ξ, α) g − (ξ,α)T gχmax + (ξ,α)T ( g − gδ )χmax −1 e −1 e −1 1 1
(ξ, α)
(ξ, α) =( | (ξ,α)T | (ξ,α)T g|2 dξ) 2 + ( ( g − gδ )|2 dξ) 2 . −1 −1 |ξ|>ξmax e |ξ|≤ξmax e
≤
e (ξ,α)T
Due to (2.6), we have
1
(ξ, α) 2 2 | (ξ,α)T g| dξ) ( −1 |ξ|>ξmax e
=(
=( ≤
When |ξ| ≥ |
1 1/α , (κα | cos( πα 2 )|T )
|ξ|>ξmax
1 2 2 |f (ξ)| dξ)
|ξ|>ξmax −s M. ξmax
s −2
(1 + |ξ| )
1 (1 + |ξ| ) |f (ξ)| dξ) 2 s 2
it obtains from Lemma 2 that
κα |ξ|α κα |ξ|α = πα α πα (κα cos( −(κα | cos( )|ξ| T ) )||ξ|α T ) 2 2 1−eα 1−e κα |ξ| ≤ ≤ 2κα |ξ|α ; 1 − e−1
(ξ, α) | ≤
(ξ,α)T e −1
and for |ξ|
a. We first get the datum g = {g(xk )}N k=1 representing values of g(x) on an equidistant grid, and then add a random uniformly distributed perturbation over (− 1, 1) to each data of g to obtain perturbation datum gδ = g + λrand(size(g)),
(5.1)
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
where λ is the noise level, and xk = (k − 1)h − a, h =
N 1
δ (gkδ − gk )2 . δ = g − g l2 = N
2a N−1 .
53
So the noise δ is measured by (5.2)
k=1
Usually it is difficult to get the analytic solution of the fractional partial differential equation (2.4) with some boundary conditions, so we have to obtain the datum g by numerical methods for a given f(x) at first. The numerical procedure of the Fourier regularization method is described as follows: Step 1. Obtaining the datum g by numerical method; Step 2. Calculating the vector gδ by Fast Fourier Transform (FFT); Step 3. Computing the vector N −1 2 Θ(ξ, α) )χ(ξ ) , gδ(ξ k k eΘ(ξ,α)T − 1 k=− N
(5.3)
2
with ξ k = πk/a and χ(ξ k ) = 1 for |ξ k | ≤ ξ max or else χ(ξ k ) = 0; Step 4. Acting the Inverse Fourier Transform (IFFT) on (5.3), and producing the vector {fξδmax (xi )}N i=1 . The accuracy of regularization result is given in the l2 norm defined by
N 1
δ Ef = f − fξmax l2 = (fk − (fξδmax )k )2 . N k=1
For the difficulty to find the analytic solution of (2.4), here we take finite difference method to get the datum g for a given f(x) [14,15]. And we discretize the Eq. (2.4) by a fractional Crank–Nicholson scheme within a finite interval [− a, a] as follows, un+1 − uni i = τ
i+1 N−i+1
1 1 α n+1 ωkα un+1 κα α [p ωk ui−k+1 + q i+k−1 ] 2 h k=0
k=0
i+1 N−i+1
1 1
ωkα uni+k−1 ] + κα α [p ωkα uni−k+1 + q 2 h k=0
k=0
(5.4)
1 1 + fin+1 + fin + O(τ 2 + h), 2 2 where uni = u(xi , tn ), τ and h are time and space step size respectively. The left and right Riemann–Liouville fractional derivatives in (2.4) can be respectively described as ([14,15]) x+a ]+1 h
1 α ωkα u(x − (k − 1)h, t) + O(h), −a Dx u(x, t) = α h k=0 a−x [ ]+1 h
1 α = α ωkα u(x − (k − 1)h, t) + O(h), x Da u(x, t) h k=0 α where ωkα = (−1)k . k [
54
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56 1.5
1.5 Numerical Exact
1
1.5 Numerical Exact
1
0.5
0.5
0.5
0
0
0
−0.5
−0.5
−0.5
−1
−1
−1
−1.5 −10
−5
0
5
10
−1.5 −10
−5
0
5
Numerical Exact
1
10
−1.5 −10
−5
0
5
10
Fig. 1. Comparison of the exact solution f with the regularization solution fξδmax with different noise for α = 1.1. 1.5
1.5 Numerical Exact
1
1.5 Numerical Exact
1
0.5
0.5
0.5
0
0
0
−0.5
−0.5
−0.5
−1
−1
−1
−1.5 −10
−5
0
5
10
−1.5 −10
−5
0
5
Numerical Exact
1
10
−1.5 −10
−5
0
5
10
Fig. 2. Comparison of the exact solution f with the regularization solution fξδmax with different noise for α = 1.5.
x2
−4 Example. Take f (x) = ( x4 − 3x and T = 1 in (2.4). It is in fact that f(x) approaches to 0 as |x| > 6, then the 2 ) e numerical experiment is implemented in the interval [− 10, 10]. In the numerical computation, we choose κα = 1, p = q = 0.5, s = 2, and N = 64. And f 2 can be computed by Maple, 1 i.e., M = 4.452844182. From Theorem 3, we can take ξmax = ( Mδ ) α+2 . Fig. 1–3 shows the pictures of the exact solution f and the regularization solution fξδmax for α = 1.1, 1.5, 1.9, respectively. And Table 1 lists the regularization parameter ξ max and the error Ef for different noise level λ, the error decreases when λ becomes small. From the numerical experiments, it is noted that the error decreases with δ becoming small for α ∈ (1, 2). 3
1.5
1.5 Numerical Exact
1
1.5 Numerical Exact
1
0.5
0.5
0.5
0
0
0
−0.5
−0.5
−0.5
−1
−1
−1
−1.5 −10
−5
0
5
10
−1.5 −10
−5
0
5
Numerical Exact
1
10
−1.5 −10
−5
0
Fig. 3. Comparison of the exact solution f with the regularization solution fξδmax with different noise for α = 1.9.
5
10
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
55
Table 1 The regularization parameter ξ max and the error Ef for different noise level λ. λ
10−1
10−2
10−3
10−4
10−5
α = 1.1
ξ max Ef
4.0675 0.0406
8.6406 0.0398
17.9457 0.0397
37.3143 0.0397
79.0675 0.0397
α = 1.5
ξ max Ef
3.4583 0.0337
6.6125 0.0198
12.7839 0.0106
24.4357 0.0102
48.6976 0.0102
α = 1.9
ξ max Ef
3.1491 0.0375
5.4188 0.0197
9.7296 0.0127
17.9282 0.0099
32.0046 0.0098
6. Conclusions We have extended two regularization methods, being ever used to regularize the classical pure diffusion equation, to solve the inverse space fractional diffusion equation which is ill-posed. Rigorous error estimates are established for both of the methods, and extensive numerical experiments are performed for Fourier regularization method. We hope this work will be interesting and/or helpful for further discussing the inverse problem of space fractional partial differential equations. Acknowledgements This work was supported by the Program for New Century Excellent Talents in University under Grant No. NCET-090438, the National Natural Science Foundation of China under Grant No. 10801067 and No. 11271173, the Fundamental Research Funds for the Central Universities under Grant No. lzujbky-2010-63 and No. lzujbky-2012-k26, the National Basic Research Program of China 973 Program No. 2011CB706900, and the Starting Research Fund from the Xi’an University of Technology under Grant No. 108-211206. References [1] F. Berntsson, A spectral method for solving the sideways heat equation, Inverse Problems 15 (1999) 891–906. [2] A.N. Bondarenko, D.S. Ivaschenko, Numerical methods for solving inverse problems for time fractional diffusion equation with variable coefficient, Journal of Inverse and III-posed Problems 17 (2009) 419–440. [3] J. Cheng, J. Nakagawa, M. Yamamoto, T. Yamazaki, Uniqueness in an inverse problem for a one-dimensional fractional diffusion equation, Inverse Problems 25 (2009) 115002. [14] A.S. Chaves, A fractional diffusion equation to describe Lévy flights, Physics Letters A 239 (1998) 13–16. [5] W.H. Deng, Finite element method for the space and time fractional Fokker–Planck equation, SIAM Journal on Numerical Analysis 47 (2008) 204–226. [6] F.F. Dou, C.L. Fu, Determining an unknown source in the heat equation by a wavelet dual least squares method, Applied Mathematics Letters 22 (2009) 661–667. [7] F.F. Dou, C.L. Fu, F.-L. Yang, Optimal error bound and Fourier regularization for identifying an unknown source in the heat equation, Journal of Computational and Applied Mathematics 230 (2009) 728–737. [8] L. Eldén, F. Berntsson, T. Regi´nska, Wavelet and Fourier method for solving the sideways heat equation, SIAM Journal on Scientific Computing 21 (2000) 2187–2205. [9] V.J. Ervin, N. Heuer, J.P. Roop, Numerical approximation of a time dependent, nonlinear, space-fractional diffusion equation, SIAM Journal on Numerical Analysis 45 (2007) 572–591. [10] H.C. Fogedby, Lévy flights in random environments, Physical Review Letters 73 (1994) 2517–2520. [11] C.L. Fu, C.Y. Qiu, Wavelet and error estimation of surface heat flux, Journal of Computational and Applied Mathematics 150 (2003) 143–155. [12] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker–Planck equation, Journal of Computational and Applied Mathematics 166 (2004) 209–219. [13] J.J. Liu, M. Yamamoto, A backward problem for the time-fractional diffusion equation, Applicable Analysis 89 (2010) 1769–1788. [14] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection-dispersion flow equations, Journal of Computational and Applied Mathematics 172 (2003) 65–77. [15] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Applied Numerical Mathematics 56 (2006) 80–90. [16] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Physics Report 339 (2000) 1–77.
56
W. Tian et al. / Mathematics and Computers in Simulation 85 (2012) 45–56
[17] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, Journal of Physics A: Mathematical and General 37 (2004) R161–R208. [18] D.A. Murio, Stable numerical solution of a fractional-diffusion inverse heat conduction problem, Computers and Mathematics with Applications 53 (2007) 1492–1501. [19] I. Daubechies, Ten lectures on wavelets, in: CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, 1992. [20] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [21] I. Podlubny, Geometric and Physical interpretation of fractional integration and fractional differentiation, Fractional Calculus and Applied Analisys 5 (2002) 367–386. [22] A.L. Qian, A new wavelet method for solving an ill-posed problem, Applied Mathematics and Computation 203 (2008) 635–640. [23] T. Regi´nska, Sideways heat equation and wavelets, Journal of Computational and Applied Mathematics 63 (1995) 209–214. [24] T. Regi´nska, L. Eldén, Solving the sideways heat equation by a wavelet-Galerkin method, Inverse Problems 13 (1997) 1093–1106. [25] J.P. Roop, Variational Solution of the Fractional Advection-Dispersion Equation, The Graduate School of Clemson University, 2004. [26] S. Samko, A. Kilbas, O. Marichev, Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach, London, 1993. [27] E. Scalas, R. Gorenflo, F. Mainardi, Fractional calculus and continuous-time finance, Physica A 284 (2000) 376–384. [28] H. Wei, W. Chen, H.G. Sun, X.C. Li, A coupled method for inverse source problem of spatial fractional anomalous diffusion equations, Inverse Problems in Science and Engineering 18 (2010) 945–956. [29] X.T. Xiong, C.L. Fu, Determining surface temperature and heat flux by a wavelet dual least squares method, Journal of Computational and Applied Mathematics 201 (2007) 189–207.
