# A Mollification Regularization Method for a Fractional-Diffusion Inverse Heat Conduction Problem

Authors: Zhi-Liang Deng; Xiao-Mei Yang; Xiao-Li Feng
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The ill-posed problem of attempting to recover the temperature functions from one measured transient data temperature at some interior point of a one-dimensional semi-infinite conductor when
the governing linear diffusion equation is of fractional type is discussed. A simple regularization method based on Dirichlet kernel mollification techniques is introduced. We also proposea priori and a posteriori parameter choice rules and get the corresponding error estimate between the exact solution and its regularized approximation. Moreover, a numerical example is provided to verify our theoretical results.


---

## Body

## 1. Introduction

It is well known that the classical diffusion equation requires second spatial derivatives and first time derivatives. However, people are shifting their partial focus to fractional-order differential equations with the realization that the use of fractional-order derivatives and integrals leads to formulas of certain physical processes which is more economical and useful than the classical approach in terms of Fick's laws of diffusion [1–4]. Some time fractional diffusion equations involving only a first-order spatial derivative and a half-order time derivative are successfully used for modeling some anomalous diffusion physical phenomena (see, e.g., [5–7] and the references therein). One can also find some applications of such fractional diffusion equations in several scenarios in [2, 7], such as relaxation to equilibrium in systems (such as polymers chains and membranes) with long temporal memory, anomalous transport in disordered systems, diffusion on fractals, and to model non-Markovian dynamical processes in protein folding. Fractional diffusion equations have been extensively investigated both in theory itself and numerical computation based on the broad applications in many application fields, especially in describing phenomena related to anomalous diffusion processes. The following is a partial list of articles which contain theoretical results and numerical tests. Some fundamental solutions and Green functions of fractional differential equations are given in [8, 9]. Finite difference techniques [10–12] and finite element methods [13, 14] provide some efficiency numerical results for several kinds of fractional differential equations. In [15], the decomposition method is used to construct analytical approximate solutions of time-fractional wave equation subject to specified boundary conditions. An analytical solution of a fractional diffusion equation by Adomian decomposition method is presented in [16]. One can also refer [17–19] to the regularity of the solution of fractional diffusion equations and some a priori estimates.In this paper, we consider the following one-dimensional fractional diffusion problem on a semi-infinite slab: suppose the temperaturef(t) at some interior point x=x0>0 is approximately measurable. For convenience, we set x0=1. The temperature u(x,t) at the 0<x<1 is desired and unknown. The mathematical description of the fractional diffusion problem is listed in the following. The unknown temperature u(x,t) satisfies
(1)∂∂xu(x,t)=-1κ∂1/2∂t1/2u(x,t)+u0πκt,x>0,t≥0,u(1,t)=f(t),
where κ is the constant diffusivity coefficient, u0=u(x,0), a constant. The half-time differentiation (∂1/2/∂t1/2)u(x,t) indicates the Riemann-Liouville fractional derivative with order α=1/2 which is given by the convolution integral
(2)(∂αu)(t)=1Γ(1-α)ddt∫0t‍u(s)(t-s)-αds,0≤t≤T,0<α<1,
where Γ(·) is the Gamma function. More detailed information on fractional derivatives and a general historical perspective may be found in [7, 20, 21]. The situations we will treat are those in which the system is initially at equilibrium, so that
(3)u(x,t)=u0,t<0,x≥0.
At t=0 a perturbation of the system commences by some unspecified process occurring at the boundary. During times of interest T, this perturbation does not affect regions remote from the x=0 boundary, so that the relationship
(4)limx→+∞u(x,t)=u0,t≤T
applies. In [22], Murio discussed the similar problems using mollification method with Gauss kernel. The idea used in this current work is a development of the ideas in [22]. However, there is only formal stability and discrete error discussion in the cited paper. Here, we give some analysis of error estimates under a priori and a posteriori regularization parameter, and a comparison of those two choice methods in our regularization method. In the following, the regularization method is outlined.In order to simplify the Fourier analysis of fractional diffusion problem (1), and in the rest of the paper, we assume without loss of generality, κ=1, u0=0. We also assume that all the functions involved are L2(ℝ) and use the corresponding L2 norm, as defined in the following, to measure errors:
(5)∥y∥=(∫ℝ‍|y(t)|2dt)1/2.
If the Fourier transform of a function y(t) is written as
(6)y^(ξ)=12π∫-∞+∞‍y(t)e-iξtdt,i=-1,
then the Fourier inversion formula reads
(7)y(t)=12π∫-∞+∞‍y^(ξ)eiξtdξ.
Fourier transforming (1) with respect to variable t, it follows that
(8)∂∂xu^(x,ξ)=-iξu^(x,ξ),x>0,ξ∈ℝ
according to the fact that [7]
(9)∂1/2∂t1/2h(t)^=iξh^(ξ).
The first-order ordinary differential equation (8) has the general solution
(10)u^(x,ξ)=u^(1,ξ)eiξ(1-x)=f^(ξ)eiξ(1-x),
or equivalently,
(11)e-iξ(1-x)u^(x,ξ)=f^(ξ).
Now, we can write the problem (1) as
(12)A^(ξ,x)u^(x,ξ)=f^(ξ),
where A^(ξ,x):=e-iξ(1-x) is a multiplication operator.Since we requireu^(x,·)∈L2(ℝ) and u0=0, we see from (4) and (10) that the real part of iξ is positive. Thereby, the multiplication operator A^(ξ,x) is a bounded linear operator for 0≤x<1. Then its inverse operator is unbounded. Therefore problem (12) is linear ill posed. For the general theory of linear ill posed problems, we can refer to [23]. It is worth pointing out that ill posed problems of a large number of diffusion equations, both fractional-order as well as integral order, have been discussed by many authors. Yang et al. [24–28] discuss the identification of source terms for some integral-order diffusion equations using some regularization strategies. Hon et al. [29, 30] apply some meshless methods to the ill posed problems of heat conduction equations. In [17, 22, 29, 31, 32], some uniqueness results and numerical methods are given for some fractional diffusion ill posed problems. Here, we apply a simple stabilizing method, namely, the mollification method with Dirichlet kernel [33], to stabilize the problem (1). Suppose that the measured data function fδ(t) satisfies
(13)∥f-fδ∥≤δ,
where δ is noise level. Take the Dirichlet function
(14)Dν(t):=2πsinνtt
as the mollifier kernel, where ν is a positive constant. Define operator Mν as
(15)Mνu(x,t):=∫-∞+∞‍Dν(τ)u(x,t-τ)dτ.
Then, we have the following associated problem: for some δ>0, find Mνuδ(x,t) that satisfies
(16)∂∂xMνuδ(x,t)=-∂1/2∂t1/2Mνuδ(x,t),x>0,t≥0,Mνuδ(1,t)=Mνfδ(t).
In fact, ν plays the role of regularization parameter.In Section2 an a priori parameter choice rule and the corresponding error estimate are discussed. In Section 3, we propose an a posteriori parameter choice rule and get the error bound. Numerical tests are given in Section 4.


## 2. The Error Estimate witha Priori Parameter Choice

In this section, the error estimate of the mollification regularization method will be derived under thea priori parameter choice rule. Suppose that the following source condition holds:
(17)∥u(0,·)∥≤E.
The source condition (17) is indispensable, otherwise, there can be no uniform convergence rates for any regularization method, that is, the convergence is arbitrarily slow (see Proposition 3.11 in [23]). For 0<x<1, by Parseval formula and triangle inequality, we know that
(18)∥Mνuδ(x,·)-u(x,·)∥=∥Mνuδ^(x,·)-u^(x,·)∥=∥e(1-x)iξDν^fδ^(ξ)-e(1-x)iξf^(ξ)∥≤∥e(1-x)iξDν^fδ^(ξ)-e(1-x)iξDν^f^(ξ)∥+∥e(1-x)iξDν^f^(ξ)-e(1-x)iξf^(ξ)∥:=I1+I2.
Since |eiξ|=e|ξ|/2, for I1, we get by (13)
(19)I1=∥e(1-x)iξDν^fδ^(ξ)-e(1-x)iξDν^f^(ξ)∥=(∫-νν‍|e(1-x)iξ(fδ^-f^)|2dξ)1/2≤e(1-x)ν/2δ.
For I2, we use the source condition (17) and obtain
(20)I2=∥e(1-x)iξDν^f^(ξ)-e(1-x)iξf^(ξ)∥=(∫|ξ|≥ν‍|e(1-x)iξf^(ξ)|2dξ)1/2=(∫|ξ|≥ν‍|e(1-x)iξe-iξu^(0,ξ)|2dξ)1/2≤e-ν/2xE.
Therefore,
(21)∥Mνuδ(x,·)-u(x,·)∥≤e(1-x)ν/2δ+e-ν/2xE.
Minimizing the right-hand side of (21), we set e(1-x)ν/2δ=e-ν/2xE and get
(22)ν=2(logEδ)2,
which also implies
(23)∥Mνuδ(x,·)-u(x,·)∥≤2E1-xδx,for0<x<1.
Therefore, we get the following theorem.Theorem 1.
Assume conditions (13), (17) hold. If the regularization parameter ν is taken by (22), then for 0<x<1, there holds the error estimate (23).Remark 2.
The error estimate (23) is order optimal in the sense of Tautenhahn [34]. In our application ∥u(0,·)∥ is usually not known, therefore we have no exact a priori bound E and cannot choose the parameter ν according to (22). However, if selecting ν=2(log(1/δ))2, we can obtain the convergence rate
(24)∥Mνuδ(x,·)-u(x,·)∥≤2δx,for0<x<1.
Theorem1 provides no information about the convergence and convergence rates of Mνuδ(x,·) at x=0. The question is settled by our next result. We now give the error estimate at x=0 under a stronger a priori assumption
(25)∥u(0,·)∥p≤E,p>0,
where ∥·∥p denotes the norm on Sobolev space Hp(ℝ) defined by
(26)∥f∥p2:=∫-∞∞‍(1+ξ2)p|f^(ξ)|2dξ
for f∈Hp(ℝ). We only need to reestimate I2 for x=0. Under the stronger a priori bound (25), it is not hard to get
(27)I2=(∫|ξ|≥ν‍|u^(0,ξ)|2dξ)1/2=(∫|ξ|≥ν‍(1+ξ2)-p(1+ξ2)p|u^(0,ξ)|2dξ)1/2≤1(1+ν2)pE≤Eν2p.
Then, (18), (19), and (27) lead to the following error bound:
(28)∥Mνuδ(0,·)-u(0,·)∥≤eν/2δ+Eν2p.
The error bound (28) does not provide the convergence as δ→0 obviously. Hence, we need to choose a proper parameter ν. For this reason, some proper lower bound of the right-hand side of (28) should be given. In order to minimize the right-hand side of (28), it is necessary to introduce the following lemma and its proof can be found in [34].Lemma 3 (see [34]).
Let the functionf(λ):(0,a]→ℝ be given by
(29)f(λ)=λb[dlog1λ]-c
with a constant c∈ℝ and positive constants a<1, b and d, then for the inverse function f-1(λ) one has
(30)f-1(λ)=λ1/b[dblog1λ]c/b(1+o(1))forλ→0.Minimizing the right-hand side of (28), we let eν/2δ=E/ν2p. Denote e-ν/2:=λ. Simple computation shows λ(log(1/λ))-4p=δ22p/E. By using (30) in Lemma 3, we obtain
(31)λ=δ22pE(logEδ22p)4p(1+o(1)),forδ→0,
which also implies that
(32)ν=2log2(δ22pE(logEδ22p)4p)(1+o(1)),forδ→0.
Therefore, the following estimate holds
(33)∥Mνuδ(0,·)-u(0,·)∥≤CE(2log2(δ22pE(logEδ22p)4p))(-2p),
where C is a constant, which also shows that the convergence rate at x=0 is logarithmic.Theorem 4.
Assume conditions (13), (25) hold. If the regularization parameter ν is taken as given by (32), then the error estimate (33) holds.Remark 5.
From the error estimate (23), as x→0, we see that the accuracy of regularization solution becomes progressively lower and even cannot get convergence for x=0. This is common in the theory of ill posed problems. Nevertheless, if a stronger a priori assumption (25) is imposed, the regularization solution converges to the exact solution at x=0, but only in a slower way.


## 3. The Error Estimate withA Posteriori Parameter Choice

In this section, we consider thea posteriori regularization parameter choice rule. Choose the regularization parameter ν as the solution of the equation
(34)d(ν):=∥Dν^fδ^(·)-fδ^(·)∥=τδ,
where τ>1 is a constant. To establish existence and uniqueness of solution of (34), we need the following lemma.Lemma 6.
Ifδ>0, then there hold the following:(a)
d(ν) is a continuous function;(b)
limν→0d(ν)=∥fδ^∥;(c)
limν→+∞d(ν)=0;(d)
d(ν) is a strictly decreasing function.The proof is very easy and we omit it here.Denote(35)ζ(x,t):=Mνuδ(x,t)-u(x,t).
We give the main result of this section as follows.Theorem 7.
Assume the conditions (11) and (17) hold and τ>1. Take the solution ν of (34) as the regularization parameter, then there holds the error estimate for 0<x<1:
(36)∥Mνuδ(x,t)-u(x,t)∥≤C(E,τ)E1-xδx,
where C(E,τ) is a constant depending on E,τ.Proof.
By (35), the Parseval formula, triangle inequality, and (17), we have
(37)∥ζ(0,·)∥=∥Mνuδ(0,·)-u(0,·)∥=∥Mνuδ^(0,·)-u^(0,·)∥=∥Dν^uδ^(0,·)-u^(0,·)∥=∥Dν^eiξfδ^-eiξf^∥=∥Dν^eiξ(fδ^-f^)+(Dν^-1)eiξf^∥≤∥Dν^eiξ(fδ^-f^)∥+E≤eν/2δ+E.
By virtue of (34), we know that
(38)τδ=∥Dν^fδ^-fδ^∥=∥(Dν^-1)(fδ^-f^+f^)∥≤∥(Dν^-1)(fδ^-f^)∥+∥(Dν^-1)f^∥≤δ+∥(Dν^-1)e-iξu^(0,ξ)∥≤δ+e-ν/2E,
which means
(39)(τ-1)δ≤e-ν/2E.
Thus it is obvious that
(40)eν/2δ≤Eτ-1.
Moreover, inserting (40) into (37), we get the following inequality:
(41)∥ζ(0,·)∥≤ττ-1E.
In addition,
(42)∥ζ(1,·)∥=∥Mνuδ(1,·)-u(1,·)∥=∥Mνuδ^(1,·)-u^(1,·)∥≤∥Mνuδ^(1,·)-uδ^(1,·)∥+∥uδ^(1,·)-u^(1,·)∥≤τδ+δ=(τ+1)δ.
It is easy to see
(43)∥ζ(x,·)∥2=∫-∞+∞‍|Dν^uδ^(x,ξ)-u^(x,ξ)|2dξ=∫-∞+∞‍|Dν^eiξ(1-x)fδ^(ξ)-eiξ(1-x)f^(ξ)|2dξ=∫-∞+∞‍|eiξ(1-x)(Dν^fδ^(ξ)-f^(ξ))|2dξ.
Denote θf=:Dν^fδ^(ξ)-f^(ξ). By using the Hölder inequality, we know that
(44)∥ζ(x,·)∥2=∫-∞+∞‍|eiξ(1-x)|2|θf|2x|θf|2(1-x)dξ≤[∫-∞+∞‍|(eiξθf)2(1-x)|1/(1-x)dξ]1-x[∫-∞+∞‍(θf2x)1/xdξ]x=[∫-∞+∞‍|(eiξθf)2|dξ]1-x[∫-∞+∞‍(θf2)dξ]x.
Therefore, we obtain
(45)∥ζ(x,·)∥≤∥ζ(0,·)∥1-x∥ζ(1,·)∥x,
which also implies that
(46)∥Mνuδ(x,·)-u(x,·)∥≤∥ζ(0,·)∥1-x∥ζ(1,·)∥x≤(ττ-1E)1-x((τ+1)δ)x.
Therefore, we complete the proof.


## 4. Numerical Examples

For linear heat diffusion, analytic solutions for the temperature distributionu(x,t),  x>0, t>0 in a semi-infinite solid with zero initial temperature, and u(0,t):=H(t) at the surface, are obtained using the integral equation
(47)u(x,t)=(H*∂∂tK)(x,t)=∫0t‍H(s)∂∂tK(x,t-s)ds,
where the kernel function
(48)K(x,t)=erfc(x2t)
is the temperature distribution corresponding to a unit step boundary temperature, H(t)=1,t>0. Here, erfc(·) denotes the complementary error function defined by
(49)erfc(z):=2π∫z∞‍e-t2dt.
We take the example from [22]. As an interesting and challenging test for the numerical method, we proposed the surface temperature function H(t), which is 1 between 0.2 and 0.6 and zero otherwise. The exact solutions for the FICHP, in this example, are the functions
(50)u(x,t)={0,0<t≤0.2,erfc(x2t-0.2),0.2<t≤0.6,erfc(x2t-0.2)-erfc(x2t-0.6),0.6<t<∞.
Consequently, the exact interior data temperature is given by
(51)f(t)={0,0<t≤0.2,erfc(12t-0.2),0.2<t≤0.6,erfc(12t-0.2)-erfc(12t-0.6),0.6<t<∞.
Suppose the vector F represents samples from the function f(t). The noisy discrete data function is generated by adding random errors to the exact data function, at every grid point, that is,
(52)Fδ(n)=F(n)+ϵn,|ϵn|≤δ,n=0,1,…,Nt,
where the (ϵn)'s are independent random Gaussian variables with variance σ2=ϵ2. The absolute and relative weighted l2 errors for the recovered interior temperatures are calculated as
(53)ei≔[1Nt+1∑n=0Nt‍|Mνuδ(x,n)-u(x,n)|]1/2,(54)er:=[(1/(Nt+1))∑n=0Nt|Mνuδ(x,n)-u(x,n)|]1/2[(1/(Nt+1))∑n=0Nt|u(x,n)|2]1/2,
respectively.It is easy to implement the algorithm described in (16). Using an available fast Fourier transform (FFT) subroutine, a simple program was written to test the algorithm for the above example in Matlab. Firstly, some tests are implemented to verify the effectiveness of the regularization method under a priori parameter choice rule at several interior points. In addition, we would like to compare the a posteriori parameter choice rule (34) with the a priori parameter choice rule (22). The a priori bound E in (22) can be calculated easily in this example. The Newton's bisection is used to solve (34), where we choose τ=1.1.Figures1 and 2 provide the comparisons between the exact solutions and regularization solutions with error bound δ=0.01 and δ=0.001 at interior point x=0.1,0.5,0.9 using a priori parameter choice rule (22), respectively. Figure 3 gives the comparison between the exact solutions and regularization solutions with error bound δ=0.001 at boundary point x=0 using a priori parameter choice rule (32). And intuitively, it seems that there are better numerical effects for closer distance from x=1. We list the error in Tables 1 and 2 to verify our result. Figure 4 demonstrates the comparisons between a priori and a posteriori parameter choice rules and Table 3 tells us that it has better effect using a posteriori parameter choice rule than a priori parameter choice rule. We also see that the difference of the numerical results between a priori choice rule and a posteriori choice rule is slight, which agrees with our theoretical results.Table 1
Errors witha priori rule.
x
δ
=
0.01
δ
=
0.001
e
r
e
i
e
r
e
i
0.1
0.1590
1.2154
0.1000
0.7643
0.5
0.0596
0.2729
0.0517
0.2370
0.9
0.0173
0.0497
0.0154
0.0443Table 2
Errors witha priori rule at x=0.
x
δ
=
0.001
e
r
e
i
0
0.1873
1.6748Table 3
The comparison betweena priori and a posteriori rules with δ=0.01.
x
A  priori
A  posteriori
0.8
e
r
=
0.0321
e
i
=
0.1063
e
r
=
0.0313
e
i
=
0.1010
0.2
e
r
=
0.1175
e
i
=
0.7817
e
r
=
0.1085
e
i
=
0.7214The comparisons of exact solution and regularization solution at (a)x=0.9, (b) x=0.5, (c) x=0.1, ν=34.3956, δ=0.01.
(a)
(b)
(c)The comparisons of exact solution and regularization solution at (a)x=0.9, (b) x=0.5, (c) x=0.1, ν=83.1949, δ=0.001.
(a)
(b)
(c)Figure 3
The comparisons of exact solution and regularization solution atx=0, ν=83.1949, δ=0.001.The comparisons of exact solution and regularization solution at (a)x=0.8, (b) x=0.2 with δ=0.01.
(a)
(b)



---

*Source: 109340-2013-02-03.xml*