# A Discretized Tikhonov Regularization Method for a Fractional Backward Heat Conduction Problem

Authors: Zhi-Liang Deng; Xiao-Mei Yang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We propose a numerical reconstruction method for solving a time-fractional backward heat conduction problem. Based on the idea of reproducing kernel approximation, we reconstruct the unknown initial heat distribution from a finite set of scattered measurements of transient temperature at a fixed final time. The standard Tikhonov regularization technique using the norm of reproducing the kernel Hilbert space as the penalty term is adopted to provide a stable solution when the measurement data contains noise. Numerical results indicate that the proposed method is efficient.


---

## Body

## 1. Introduction

LetΩ
⊂
ℝ
d, d
∈
ℕ, be a bounded domain with sufficiently smooth boundary ∂
Ω. Consider the following initial boundary value problem for time-fractional diffusion equation (TFDE):
(1)
0
D
t
γ
u
(
x
,
t
)
=
ℒ
u
,
x
∈
Ω
,
t
∈
[
0
,
T
]
,
u
(
x
,
0
)
=
f
(
x
)
,
x
∈
Ω
-
,
ℬ
u
(
x
,
t
)
:
=
μ
u
(
x
,
t
)
+
β
∂
u
∂
ν
=
0
,
x
∈
∂
Ω
,
t
∈
[
0
,
T
]
.
Here0
D
t
γ denotes the Caputo fractional derivative with respect to t and is defined by
(2)
0
D
t
γ
ψ
(
t
)
=
{
1
Γ
(
1
-
γ
)
∫
0
t
‍
ψ
′
(
τ
)
(
t
-
τ
)
γ
d
τ
,
0
<
γ
<
1
,
ψ
′
(
t
)
,
γ
=
1
.
Γ
(
·
) is the Gamma function and ℒ is a symmetric uniformly elliptic operator and T
>
0 is a fixed final time, μ and β are constants, and ν is the outward unit normal vector of the domain Ω. In what follows, let 0
<
γ
<
1 and ℒ be given by
(3)
ℒ
u
(
x
)
:
=
∑
i
=
1
d
‍
∂
∂
x
i
(
∑
j
=
1
d
‍
D
i
j
(
x
)
∂
∂
x
j
u
(
x
)
)
+
Q
(
x
)
u
(
x
)
,
x
∈
Ω
,
where Q
∈
C
(
Ω
-
), Q
(
x
)
≤
0 for x
∈
Ω
- and D
i
j
=
D
j
i
∈
C
1
(
Ω
-
), 1
≤
i, j
≤
d. Moreover, we assume that there exists a positive constant D
>
0 such that
(4)
D
∑
i
=
1
d
‍
η
i
2
≤
∑
i
,
j
=
1
d
‍
D
i
j
(
x
)
η
i
η
j
,
x
∈
Ω
-
,
η
∈
ℝ
d
.Recently, people are shifting their partial focus to fractional-order differential equations (FDEs) with the realization that the use of fractional-order derivatives and integrals leads to formulas of certain physical processes (for instance, some anomalous diffusion processes) which is more economical and useful than the classical approach in terms of Fick’s laws of diffusion. In this paper, we consider the fractional-order partial differential equation (FPDE) in (1), which is obtained from the standard diffusion equation by replacing the first-order time derivative with a fractional derivative of order γ, with 0
<
γ
<
1. Different models using this kind of FDEs have been proposed [1–3], and there has been significant interest in developing numerical schemes for their solution. Physically, the time-fractional partial differential equations describe the continuous time random walk problems (the non-Markovian process). The physical interpretation of the fractional derivative is that it represents a degree of memory in the diffusing material [4]. Actually, the convolution integral in the definition of the Caputo fractional-order derivatives for γ
∈
(
0,1
) at time t requires all the knowledge of classical derivative ψ
′
(
τ
) for τ
∈
(
0
,
t
), which reflects the “memory effect” of fractional derivatives. The utilization of the memory effect of fractional derivatives comes with a high cost regarding numerical solvability. Any algorithm using a discretization of a noninteger derivative has to take into account its nonlocal structure which means in general a high storage requirement and great overall complexity of the algorithm. Similar convolution model can be used to depict the tumor-immune system [5]. Fractional dimensional model is also used to characterize the binary images of DNA [6] and oscillators [7, 8]. In [9], a general approach is proposed to approximate ideal filters based on fractional calculus from the point of view of systems of fractional order. Recently, numerous attempts to solve TFDE can be found in the literature [4, 10–14]. We can refer to [15–18] for more specified knowledge about fractional calculus.The backward problem governed by time-fractional partial differential equation in (1) is to recover the heat distribution at any earlier time 0
≤
t
<
T from the measurement u
ϵ
(
x
,
T
), written as u
T
ϵ, which is noise-contaminated data for the exact temperature u
(
x
,
T
):
(5)
∥
u
ϵ
(
·
,
T
)
-
u
(
·
,
T
)
∥
L
2
(
Ω
)
≤
ϵ
for some known noise level ϵ
>
0. Such inverse problems have been considered by several authors. Based on the eigenfunction expansions, Sakamoto and Yamamoto [19] establish the unique existence of the weak solution and the asymptotic behavior as the time t goes to ∞ for the forward problem and prove the stability and uniqueness in the backward problem in time. For the one-dimensional case, Liu and Yamamoto [20] propose a regularizing scheme by the quasi-reversibility to restore the stability of the backward problem. In [21], a regularization by projection is applied to the same problem as in [20] and the corresponding convergence rates are obtained undera priori anda posteriori parameter choice strategies, respectively. Here, we pay our attention to the situation of stable reconstruction of the initial heat distribution f
(
x
) from some scattered noisy data of u
ϵ
(
·
,
T
). More specifically, the data u
ϵ
(
·
,
T
) are collected only at a finite set of points {
z
1
,
z
2
,
…
,
z
m
}
∈
Ω. We then reconstruct the initial temperature distribution u
(
x
,
0
) from the scattered noisy data at t
=
T. For solving the backward diffusion problem, we employ a discretized Tikhonov regularization by the Ritz approach coupled with the reproducing kernel Hilbert space (RKHS), which is proposed in [22].The Tikhonov regularization method has been widely studied and applied to all varieties of ill-posed problems [23, 24]. The discretized Tikhonov regularization method and its relative theories are also explored in detail [24]. We adopt the Tikhonov regularization method by a reproducing kernel Hilbert space into the backward problem (1). As we know, the theory and practice of reproducing kernel are a fast growing research area. The numerical methods by RKHS have been also rapidly developed in recent years [25, 26]. These developments are due to the increasing interest in the use of reproducing kernel for the solution of mathematical and engineering problems, for instance, machine learning [27], signal processing [28], stochastic processes [29], wavelet transforms [30], and so forth. For the details of RKHS, we are able to refer to [31]. However, to the authors’ knowledge, there are few applications of RKHS to inverse problems. We provide the partial list of the recent works. Takeuchi and Yamamoto [22] prove the convergence of the discretized Tikhonov regularization method by RKHS. Hon and Takeuchi [32] apply this method into a backward heat conduction problem for parabolic-type partial differential equations. In reproducing the kernel Hilbert space settings, an inverse source identification problem for parabolic equations is considered in [33]. In [34], Saitoh discusses comprehensively the corresponding applications of RKHS in inverse problems.The remainder of this paper is composed of five sections. In Section2, we discuss Green’s function for problem (1) and use it to construct the reproducing kernel. In Section 3, we state the reconstruction method of the fractional backward diffusion problem by using the reproducing kernel. In Section 4, some numerical examples are given to illustrate the effectiveness of our method. A summary is made in the Section 5. Finally, we list some existing knowledge about the reproducing kernel Hilbert space in the Appendix.


## 2. Green’s Function and the Reproducing Kernel

In this section, we explore Green’s function of system (1) and use it to construct the reproducing kernel. Firstly, let function G
(
x
,
t
;
ξ
) be Green’s function of the system (1); that is, G satisfies the following equations in distribution’s sense:
(6)
0
D
t
γ
G
(
x
,
t
;
ξ
)
=
ℒ
G
(
x
,
t
;
ξ
)
,
x
,
ξ
∈
Ω
,
t
∈
[
0
,
T
]
,
G
(
x
,
0
;
ξ
)
=
δ
(
x
-
ξ
)
,
x
,
ξ
∈
Ω
,
ℬ
G
(
x
,
t
;
ξ
)
:
=
μ
G
(
x
,
t
;
ξ
)
+
β
∂
G
∂
ν
=
0
,
x
∈
∂
Ω
,
t
∈
[
0
,
T
]
.
It is easy to know that problem (6) is equivalent to the following problem:
(7)
0
D
t
γ
G
(
x
,
t
;
ξ
)
-
ℒ
G
(
x
,
t
;
ξ
)
=
δ
(
x
-
ξ
)
δ
(
t
)
,
11111111111111111111
x
,
ξ
∈
Ω
,
t
∈
[
0
,
T
]
,
G
(
x
,
0
;
ξ
)
=
0
,
x
,
ξ
∈
Ω
,
ℬ
G
(
x
,
t
;
ξ
)
:
=
μ
G
(
x
,
t
;
ξ
)
+
β
∂
G
∂
ν
=
0
,
x
∈
∂
Ω
,
t
∈
[
0
,
T
]
.
Let us employ the Laplace transform to solve the system (6). The Laplace transform of a function φ on ℝ
+ is defined by
(8)
ϕ
(
s
)
:
=
∫
0
∞
e
-
s
t
φ
(
t
)
d
t
,
Re
(
s
)
>
0
.
The Laplace transform of the Caputo fractional derivative is given by [10]
(9)
∫
0
∞
e
-
s
t

D
t
γ
0
φ
(
t
)
d
t
=
s
γ
ϕ
(
s
)
-
∑
k
=
0
J
-
1
‍
s
γ
-
k
-
1
0
D
t
k
φ
(
0
)
,
where J
-
1
<
γ
≤
J. The Caputo fractional derivative appears more suitable to be treated by the Laplace transforming in that it requires the knowledge of initial values of the function and of its integer derivatives of order k
=
1,2
,
…
,
J
-
1. By the Laplace transforming about the time variable t, the system (6) becomes
(10)
s
γ
g
(
x
,
s
;
ξ
)
-
s
γ
-
1
δ
(
x
-
ξ
)
=
ℒ
g
(
x
,
s
;
ξ
)
,
x
,
ξ
∈
Ω
,
Re
(
s
)
>
0
,
μ
g
(
x
,
s
;
ξ
)
+
β
∂
g
∂
ν
(
x
,
s
;
ξ
)
=
0
,
where g denotes the Laplace transform of Green’s function G. It should be noted here that the system in (1) is only defined on t
∈
[
0
,
T
], not on [
0
,
+
∞
). When we make the Laplace transform, some necessary preprocess, for example, the function continuation technique, needs to be done on the solution of (1) to satisfy the condition of the Laplace transform. Because we do not use the value of u for t
>
T, the condition can be satisfied easily.Applying this technique of eigenfunction expansions to problems (10), we have that
(11)
g
(
x
,
s
;
ξ
)
=
∑
n
=
1
∞
‍
s
γ
-
1
s
γ
+
λ
n
2
φ
n
(
ξ
)
φ
n
(
x
)
,
where φ
n
(
x
) is the nth orthonormal eigenfunction and λ
n is the corresponding eigenvalue to the Sturm-Liouville problem
(12)
ℒ
φ
(
x
)
+
k
2
φ
(
x
)
=
0
,
subject to the boundary conditions
(13)
μ
φ
(
x
)
+
β
∂
φ
∂
ν
=
0
.
Taking the Laplace inverse of (11), we have that
(14)
G
(
x
,
t
;
ξ
)
=
∑
n
=
1
∞
‍
E
γ
(
-
λ
n
2
t
γ
)
φ
n
(
ξ
)
φ
n
(
x
)
,
where E
γ
(
z
) is the Mittag-Leffler function defined by
(15)
E
γ
(
z
)
=
∑
k
=
0
∞
‍
z
k
Γ
(
γ
k
+
1
)
.
For the details of the Mittag-Leffler function, one can refer to [3]. Subsequently, one can easily verify that the unique solution of system (1) with initial value u
(
x
,
0
)
=
f
∈
H can be written as
(16)
u
(
x
,
t
)
=
∫
Ω
‍
G
(
x
,
t
;
ξ
)
f
(
ξ
)
d
ξ
,
where G
(
x
,
t
;
ξ
) is Green’s function defined by (14). With the aid of (16), for each f
∈
L
2
(
Ω
), introducing the operator K
:
L
2
(
Ω
)
→
L
2
(
Ω
) by
(17)
K
f
(
x
)
:
=
∫
Ω
‍
G
(
x
,
T
;
ξ
)
f
(
ξ
)
d
ξ
,
we may formulate the inverse problem as an integral equation of the first kind
(18)
K
f
(
x
)
=
u
(
x
,
T
)
.The symmetry of Green’s function indicates that the operatorK is self-adjoint. We proceed by giving a brief account on the close connections of the ill posedness of operator equation (18) with the singular value system of operator K. According to the symmetry of operator K, we only need to discuss its eigensystem. It is easy to know that the eigensystem of K is given by
(19)
{
E
γ
(
-
T
γ
λ
n
2
)
,
φ
n
}
n
=
1
∞
,
where E
γ
(
-
T
γ
λ
n
2
) is the eigenvalue and φ
n is the corresponding eigenfunction. We can see the decay of the eigenvalues {
E
γ
(
-
T
γ
λ
n
2
)
}
n
=
1
∞ with the increase of λ
n from the following asymptotic behavior of the Mittag-Leffler function E
γ
(
·
) which can be found in [20, 21] or by the results in [3].Lemma 1 (see [20]).
Let0
<
γ
0
<
γ
1
<
1. Then there exist constants C
1
,
±
,
C
2
,
±
>
0 depending only on γ
0 and γ
1 such that
(20)
C
1
,
-
γ
e
x
1
/
γ
≤
E
γ
(
x
)
≤
C
1
,
+
γ
e
x
1
/
γ
,
∀
x
≥
0
,
(21)
C
2
,
-
Γ
(
1
-
γ
)
1
1
-
x
≤
E
γ
(
x
)
≤
C
2
,
+
Γ
(
1
-
γ
)
1
1
-
x
,
∀
x
≤
0
.
These estimates are uniform for all γ
∈
[
γ
0
,
γ
1
].Next, we consider the construction of reproducing kernel using Green’s functionG
(
x
,
t
;
ξ
). Define Φ
(
x
,
ξ
)
:
=
G
(
x
,
t
0
;
ξ
) for some t
0
∈
(
0
,
T
). The symmetry about the space variable x, ξ of Green’s function indicates that Φ
(
·
,
·
)
:
Ω
×
Ω
→
ℝ is symmetric. Now that Φ is a symmetric positive definite kernel, a unique RKHS in which the given kernel acts as the reproducing kernel can be constructed (see [26] for details). Henceforth we denote by H
t
0 the RHKS generated by the kernel Φ
(
x
,
ξ
)
=
G
(
x
,
t
0
;
ξ
). Actually, according to [26, Chapter 10], the inner product and norm on H
t
0 are defined by
(22)
(
f
,
g
)
H
t
0
:
=
∑
n
=
1
∞
‍
1
E
γ
(
-
λ
n
2
t
0
γ
)
(
f
,
φ
n
)
L
2
(
Ω
)
(
g
,
φ
n
)
L
2
(
Ω
)
,
∥
f
∥
H
t
0
2
:
=
∑
n
=
1
∞
‍
1
E
γ
(
-
λ
n
2
t
0
γ
)
|
(
f
,
φ
n
)
L
2
(
Ω
)
|
2
,
respectively. The space H
t
0 is actually given by
(23)
H
t
0
=
{
f
∈
L
2
(
Ω
)
:
∑
n
=
1
∞
‍
1
E
γ
(
-
λ
n
2
t
0
γ
)
|
(
f
,
φ
n
)
L
2
(
Ω
)
|
2
<
∞
}
.
The second inequalities in Lemma 1 show that, as the function of λ
n, 1
/
E
γ
(
-
λ
n
2
t
0
γ
)
∈
𝒮
ℐ, where 𝒮
ℐ denotes the collection of slowly increasing functions [26] defined by
(24)
𝒮
ℐ
:
=
{
f
:
ℝ
d
⟶
ℝ
∣
f
(
x
)
=
O
(
∥
x
∥
2
m
)
ℝ
d
as
∥
x
∥
2
⟶
∞
for
some
m
∈
ℕ
0
}
.
Hence, according to the theoretical results of [26, Chapter 10], we assert that the space H
t
0 is consistent with some Sobolev space H
k
(
Ω
) for k
>
d
/
2 and the norm on H
t
0 is equivalent to the norm ∥
f
∥
H
k
(
Ω
)
=
[
∑
|
α
|
≤
k
‍
∥
D
α
f
∥
L
2
(
Ω
)
2
]
1
/
2.


## 3. Formulation of the Inverse Problem and the Reconstruction Method

In order to find the initial temperature distributionf
(
x
), we would like to determine the solution of minimization problem
(25)
inf
⁡
f
∈
L
2
(
Ω
)
∥
K
f
-
u
ϵ
(
·
,
T
)
∥
L
2
(
Ω
)
.
However, the minimization element of (25) generally is a poor approximation of the desired initial function f due to the error ϵ in u
(
x
,
T
) and the ill posedness of operator equation (18). The Tikhonov regularization replaces the minimization problem (25) by the solution of a penalized least-squares problem
(26)
inf
⁡
f
∈
H
t
0
(
Ω
)
J
α
(
f
)
:
=
inf
⁡
f
∈
H
t
0
(
Ω
)
∥
K
f
-
u
ϵ
(
·
,
T
)
∥
L
2
(
Ω
)
2
+
α
∥
f
∥
H
t
0
2
with regularization parameter α
>
0. We can see from [23, Proposition 3.11] that the convergence of any regularization method can be arbitrarily slow in general. Actually, convergence rates can be given only on some subset of H
t
0, that is, undera priori assumptions on the exact data. Here, we assume the exact solution f belongs to the set of source conditions
(27)
ℳ
μ
:
=
{
f
∈
H
t
0
(
Ω
)
:
f
=
(
K
*
K
)
μ
ω
,
∥
ω
∥
L
2
(
Ω
)
≤
E
}
,
where E is thea priori bound and μ
>
0 is a constant. As in [23], we know that there exists a constant μ
0, named the “qualification” of the regularization method, such that 0
<
μ
≤
μ
0. For the Tikhonov regularization method, the qualification μ
0
=
1. However, according to (21), it is easy to know that ℳ
μ
2
⊂
ℳ
μ
1 as μ
1
>
μ
2 and ℳ
1
/
4
=
H
t
0. Therefore, here we only need consider the case of 0
<
μ
≤
1
/
4 for the Tikhonov method (26).In order to solve the minimization problem (26), some discretization scheme needs to be given. A natural way to obtain such discretization is to generate a finite dimensional approximation to the minimal element of the Tikhonov functional J
α. For this, we define a subspace V
N
:
=
span
⁡
{
Φ
(
·
,
ξ
)
:
ξ
∈
X
N
}
⊂
H
t
0, where X
N
:
=
{
ξ
1
,
ξ
2
,
…
,
ξ
N
}
⊂
Ω. This approximation by discretization is equivalent to finding the minimal norm least-squares solution of the equation
(28)
K
N
f
:
=
K
P
N
f
=
u
T
,
where P
N
:
H
t
0
→
V
N is the projection operator. Moreover, we produce a finite dimensional approximation f
α
,
N to K
†
u
T by minimizing the Tikhonov functional (26) over the finite dimensional space V
N. Denote by f
α
,
N
,
ϵ the minimizer of J
α for noise input data u
T
ϵ. It is well known that f
α
,
N
,
ϵ satisfy [23, 24]
(29)
f
α
,
N
,
ϵ
=
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
u
T
ϵ
.
Provided {
V
N
} is an expanding sequence, the convergence of the Tikhonov regularized solutions is proved [24]. Takeuchi and Yamamoto [22] show that the discretized Tikhonov regularized solutions converge to the exact solution without the monotonicity of V
N under ana priori choice strategy for N and α. However, one can see that, from the existing results, in both cases, the regularized solutions f
α
,
N converging to the exact solution depends on whether r
N
:
=
∥
K
(
I
-
P
N
)
∥
→
0,  (
N
→
∞
). Moreover, the convergence of r
N requires the knowledge of the fill distance h
X
N
,
Ω, which is defined by [26]
(30)
h
X
N
,
Ω
:
=
sup
⁡
x
∈
Ω
min
⁡
ξ
k
∈
X
N
∥
x
-
ξ
k
∥
.
The fill distance can be interpreted in various geometrical ways. For example, we can consider it as the radius of the largest ball which is completely contained in Ω and which does not contain a data site. In this sense h
X
N
,
Ω describes the largest data-site-free hole in Ω.Here, we utilize the same proof as in [32] to give the following lemma.Lemma 2.
Considerr
N
=
∥
K
(
I
-
P
N
)
∥
→
0 as N
→
∞.Proof.
Consider(31)
∥
K
(
I
-
P
N
)
∥
=
sup
⁡
f
∈
H
t
0
f
≠
θ
∥
K
(
I
-
P
N
)
f
∥
L
2
(
Ω
)
∥
f
∥
H
t
0
.
It is obvious that the kernel functionΦ
(
x
,
y
) is sufficiently smooth. Therefore, according to the error estimate in the Appendix, we can find a positive constant k such that the estimate
(32)
∥
f
-
P
N
f
∥
L
∞
(
Ω
)
≤
h
X
N
,
Ω
k
∥
f
∥
H
t
0
holds for all f
∈
H
t
0. Meanwhile, it is easy to know
(33)
∥
K
f
∥
L
2
(
Ω
)
=
[
∫
Ω
‍
|
K
f
|
2
d
x
]
1
/
2
=
[
∫
Ω
‍
|
∫
Ω
‍
G
(
x
,
T
;
ξ
)
f
(
ξ
)
d
ξ
|
2
d
x
]
1
/
2
≤
C
∥
f
∥
L
∞
(
Ω
)
,
∀
f
∈
H
t
0
,
where C
:
=
{
∫
Ω
‍
[
∫
Ω
‍
|
G
(
x
,
T
;
ξ
)
|
d
ξ
]
2
d
x
}
1
/
2 is a constant. Moreover, the property of RKHS, f
(
x
)
=
(
f
(
·
)
,
Φ
(
x
,
·
)
), leads to
(34)
∥
f
∥
L
∞
(
Ω
)
≤
sup
⁡
x
∈
Ω
Φ
(
x
,
x
)
∥
f
∥
H
t
0
.
Combining (32)–(34), we get
(35)
∥
K
(
I
-
P
N
)
f
∥
L
2
(
Ω
)
≤
M
∥
f
-
P
N
f
∥
L
∞
(
Ω
)
≤
M
h
X
N
,
Ω
k
∥
f
∥
H
t
0
,
where M is a constant. Consequently, substituting the estimate (35) into (31) and letting N
→
∞, we complete the proof.According to the classical results on the Tikhonov regularization for linear ill-posed problem (see, e.g., [23, 24]) and in view of (1), it holds that
(36)
∥
f
-
f
α
,
N
,
ϵ
∥
≤
∥
f
-
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
u
∥
+
∥
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
(
u
-
u
T
ϵ
)
∥
≤
∥
f
-
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
u
∥
+
ϵ
2
α
.
As in [35], we can estimate the noise-free term as follows:
(37)
∥
f
-
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
u
∥
=
∥
f
-
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
K
f
∥
≤
∥
[
I
-
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
K
N
]
f
∥
+
∥
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
(
K
N
-
K
)
f
∥
=
∥
[
I
-
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
K
N
]
(
K
*
K
)
μ
w
∥
+
∥
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
(
K
N
-
K
)
f
∥
≤
∥
[
I
-
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
K
N
]
×
[
(
K
*
K
)
μ
-
(
K
N
*
K
N
)
μ
]
w
∥
+
∥
[
I
-
(
K
N
*
K
N
+
α
I
)
-
1
K
N
*
K
N
]
(
K
N
*
K
N
)
μ
w
∥
+
1
α
∥
K
(
I
-
P
N
)
f
∥
≤
∥
[
(
K
*
K
)
μ
-
(
K
N
*
K
N
)
μ
]
w
∥
+
α
μ
E
+
1
α
∥
K
(
I
-
P
N
)
f
∥
≤
∥
K
*
K
-
K
N
*
K
N
∥
μ
E
+
α
μ
E
+
1
α
∥
K
(
I
-
P
N
)
f
∥
≤
C
∥
K
(
I
-
P
N
)
∥
2
μ
E
+
α
μ
E
+
1
α
∥
K
(
I
-
P
N
)
f
∥
≤
C
∥
K
(
I
-
P
N
)
∥
2
μ
E
+
α
μ
E
+
1
α
∥
K
(
I
-
P
N
)
∥
2
μ
+
1
,
where C is a constant. In view of the best possible error bound being ϵ
2
μ
/
(
2
μ
+
1
), the term ∥
K
(
I
-
P
N
)
∥ has to be chosen such that
(38)
∥
K
(
I
-
P
N
)
∥
≤
ϵ
1
/
(
2
μ
+
1
)
.
From the above discussions, we have the following theorem.Theorem 3.
Under assumptions (5) and (38), there holds that
(39)
∥
f
-
f
α
,
N
,
ϵ
∥
≤
C
ϵ
2
μ
/
(
2
μ
+
1
)
+
α
μ
E
+
ϵ
α
.
Moreover, if the regularization parameter α is chosen by α
=
O
(
ϵ
2
/
(
2
μ
+
1
)
), one then obtains the following estimate:
(40)
∥
f
-
f
α
,
N
,
ϵ
∥
≤
C
~
ϵ
2
μ
/
(
2
μ
+
1
)
,
where the constant C
~ does not depend on ϵ.


## 4. Numerical Tests

In this section, we present numerical results to illustrate the feasibility of the reconstruction method as described in the previous section.In practical situation, we only can get the scattered noisy data ofu
(
·
,
T
), that is, {
u
ϵ
(
z
1
,
T
)
,
u
ϵ
(
z
2
,
T
)
,
…
,
u
ϵ
(
z
m
,
T
)
}. As a result, instead of solving (25) we intend to deal with the following problem:
(41)
inf
⁡
f
∈
V
N
∥
K
f
-
u
ϵ
(
·
,
T
)
∥
ℝ
m
2
+
α
∥
f
∥
H
t
0
(
Ω
)
2
,
where
(42)
∥
K
f
-
u
ϵ
(
·
,
T
)
∥
ℝ
m
=
{
∑
j
=
1
m
‍
|
K
f
(
z
j
)
-
u
ϵ
(
z
j
,
T
)
|
2
}
1
/
2
.
Since V
N
=
span
⁡
{
Φ
(
·
,
ξ
)
∣
ξ
∈
{
ξ
1
,
ξ
2
,
…
,
ξ
N
}
}, the minimizer f
α
,
N
,
ϵ can be written as
(43)
f
α
,
N
,
ϵ
(
·
)
=
∑
k
=
1
N
‍
λ
~
k
Φ
(
·
,
ξ
k
)
=
∑
k
=
1
N
‍
λ
~
k
G
(
·
,
t
0
;
ξ
k
)
,
ξ
k
∈
X
N
.
From the definition of RKHS, it follows that
(44)
∥
f
α
,
N
,
ϵ
∥
H
t
0
2
=
∑
j
,
k
=
1
N
‍
λ
~
k
λ
~
j
G
(
ξ
k
,
t
0
;
ξ
j
)
.
In addition, we know that
(45)
K
G
(
z
j
,
t
0
;
ξ
k
)
=
∫
Ω
‍
G
(
z
j
,
T
;
ξ
)
G
(
ξ
,
t
0
;
ξ
k
)
d
ξ
=
G
(
z
j
,
t
0
+
T
;
ξ
k
)
for j
=
1,2
,
…
,
m and k
=
1,2
,
…
,
N. Now, it is easy to see that the coefficient vector λ
~
=
(
λ
~
k
)
k
=
1
N satisfies the following linear system:
(46)
(
A
*
A
+
α
B
)
λ
~
=
A
*
u
T
ϵ
,
where A is an m
×
N matrix, A
* is the transpose of A, and B is an N
×
N matrix defined by
(47)
A
j
,
k
=
G
(
z
j
,
t
0
+
T
;
ξ
k
)
,
j
=
1
,
…
,
m
,
k
=
1
,
…
,
N
,
B
j
,
k
=
G
(
ξ
j
,
t
0
;
ξ
k
)
,
j
,
k
=
1
,
…
,
N
.
If we truncate Green’s function to the former J terms, matrix B can be decomposed to the following product:
(48)
B
=
L
*
L
,
where L
* denotes the conjugate transpose of L and L is given by
(49)
L
=
(
L
j
,
k
)
=
(
[
E
γ
(
-
λ
j
2
t
0
γ
)
]
1
/
2
φ
j
(
ξ
k
)
)
,
j
=
1,2
,
…
,
J
,
k
=
1,2
,
…
,
N
.
Moreover, we turn to search for the minimizer λ
~ of the following functional:
(50)
F
α
(
λ
)
:
=
∥
A
λ
-
u
T
ϵ
∥
ℝ
m
2
+
∥
L
λ
∥
ℝ
J
2
.
After obtaining the vector λ
~, we substitute it into (43) and then get the regularized approximation f
α
,
N
,
ϵ.In our tests, the measurement points{
z
j
}
j
=
1
m, which are randomly generated by using the Matlab function rand(
·
), are scattered in the domain Ω. Now we generate the final measurement data at T with noise by
(51)
u
T
ϵ
=
u
(
z
,
T
)
+
2
π
ϵ
100
×
rand
(
z
)
×
norm
⁡
(
u
(
z
,
T
)
)
,
where z
=
{
z
j
}
j
=
1
m are the measurement points and rand
(
z
) generates a standard m
-dimensional random vector. To evaluate the proposed method, we compute the relative error of the reconstructed solutions denoted by R
(
f
):
(52)
R
(
f
)
=
∥
f
α
,
N
,
ϵ
(
·
)
-
f
(
·
)
∥
l
2
∥
f
(
·
)
∥
l
2
,
for
1
-dimensional
case
,
R
(
f
)
=
∥
f
α
,
N
,
ϵ
(
·
)
-
f
(
·
)
∥
∞
∥
f
(
·
)
∥
∞
,
for
2
-dimensional
case
,
where ∥
·
∥
l
2 denotes the l
2 norm and ∥
·
∥
∞ denotes the ∞ norm. Before we proceed, it is natural that we have to determine t
0, X
N, measurement points {
z
j
}
j
=
1
m, and observation time T to define A, B. According to the convergence theorem, smaller h
X
N
,
Ω yields better numerical solution, which implies that we need to choose as large N as possible. However, the ill posedness of the backward diffusion problem results in the ill condition of the matrix A
*
A, which causes us not to implement the numerical computation of the inverse of A
*
A
+
α
B when the regularizing term B is also ill conditioned. We use the following one-dimensional example to depict the change trend of the condition number of B, cond(
B
), with respect to N.Example 1.
Consider the following Dirichlet boundary value problem:(53)
0
D
t
γ
u
(
x
,
t
)
=
∂
2
u
(
x
,
t
)
∂
x
2
,
x
∈
(
0,1
)
,
t
∈
[
0
,
T
]
,
u
(
x
,
0
)
=
x
(
1
-
x
)
,
x
∈
[
0,1
]
,
u
(
0
,
t
)
=
u
(
1
,
t
)
=
0
,
t
∈
[
0,1
]
.
The forward problem has a unique solution
(54)
u
(
x
,
t
)
=
∑
n
=
1
∞
‍
d
n
E
γ
(
-
n
2
π
2
t
γ
)
sin
n
π
x
,
where the coefficient
(55)
d
n
=
2
∫
0
1
‍
u
(
x
,
0
)
sin
n
π
x
d
x
=
4
n
3
π
3
[
1
-
(
-
1
)
n
]
,
n
=
1,2
,
…
.
To clarify the numerical influence of some relative parameters butγ, we fix parameter γ
=
1
/
2 firstly.
In this test, we fixm
=
59 firstly. In Figure 1, for the cases of t
0
∈
{
1
e
-
8,1
e
-
9,1
e
-
10
}, we plot cond(
B
) versus the number of N running from 10 to 100, respectively. The displayed results in Figure 1 show that the condition number cond(
B
) increases exponentially as N increases. Nevertheless, we can remove such influence of the ill condition of B on numerical computation through modifying its small singular value as a fixed small constant ϵ
~. And in doing so, the numerical precision does not change significantly. Therefore, we can implement the proposed method without worrying so much about the size limitation of N. For the choice of t
0, we compare the computational results for some different t
0’s and T’s using exact final data u
(
x
,
T
) for Example 1. For this, a preset value for the regularization parameter α needs to be provided. Here, we simply choose α
=
1
0
-
5
(
max
⁡
⁡
(
A
*
A
)
/
max
⁡
⁡
(
B
)
) as in [32], where max
⁡
(
C
)
=
max
⁡
i
,
j
c
i
j for matrix C
=
(
c
i
j
). In the following computation, the Matlab code developed by Hansen [36, 37] is used to obtain the approximation solution for solving the discrete system (46). In addition, note that the Mittag-Leffler function is numerically realized by implementing the Matlab toolbox by Podlubny [38]. In Figure 2, we plot the relative errors for different t
0, t
0
=
1
e
-
5,1
e
-
7,1
e
-
10, versus the number of N running from 10 to 100 for ϵ
=
0 and T
=
25. The computational results show that, using smaller t
0, we have less relative error. In addition, we also see from Figure 2 that, with the increase of the number of N, the relative error becomes smaller firstly and then it remains steady after arriving to a certain extent.Figure 1
The effect of an increasing number ofN on the condition number of B for the cases when t
0
∈
{
1
e
-
8,1
e
-
9,1
e
-
10
}.Figure 2
The effect of an increasing number ofN on the relative error for the cases when t
0
∈
{
1
e
-
5,1
e
-
7,1
e
-
10
} with noise-free data.We also need to consider the effect of the numberJ of truncation term of Green’s function and the number of measurement points m on the numerical precision. With T
=
3, we display the numerical results for several m and J as the noise level ϵ
=
0 in Table 2. It can be seen that when J and m become sufficiently large, the relative errors almost remain at the level 1
0
-
4.Next, in the absence of thea priori information, we only evaluate the proposed algorithm in (46) for noisy data by using L-curve parameter choice method instead of that in Theorem 3. In Example 1, we fix N
=
59, m
=
59, and t
0
=
1
e
-
10. Table 1 reports the relative errors of f
α
,
N
,
ϵ for different noise levels ϵ and final measurement times T. These numerical results for all noisy cases are satisfactory. In general, it can been seen from Table 1 that, at the smaller T and ϵ, the numerical effects are better. In addition, when the measurement time T
=
25, the exact solution u
(
x
,
0
) and the numerical solution f
α
,
N
,
ϵ with the relative noise level ϵ
=
20
% are displayed in Figure 3. It can be observed that the method works even for the case of T
=
25 with noise level ϵ
=
20
% as well.Table 1
Example1: the relative errors with t
0
=
1
e
-
10, N
=
59, and m
=
59.
ϵ
T
1
5
10
0.1%
0.0143
0.0156
0.0177
1%
0.0295
0.0239
0.0332
2%
0.0388
0.0390
0.0387
5%
0.0483
0.0441
0.0505
10%
0.0528
0.0601
0.0703Table 2
Example1: the relative errors with t
0
=
1
e
-
10, N
=
59, and ϵ
=
0
%.
J
m
19
29
39
49
59
69
5
3.3721
e
-
3
3.4356
e
-
3
3.3769
e
-
3
3.3737
e
-
3
3.4008
e
-
3
3.3724
e
-
3
10
1.0725
e
-
3
1.5791
e
-
3
9.9642
e
-
4
9.9488
e
-
4
9.9616
e
-
4
9.9901
e
-
4
15
6.1811
e
-
4
3.4322
e
-
4
5.7376
e
-
4
3.2938
e
-
4
3.2952
e
-
4
3.3259
e
-
4
20
5.0304
e
-
3
3.5813
e
-
4
8.9305
e
-
4
4.7205
e
-
4
9.8808
e
-
4
5.4753
e
-
4
25
3.7310
e
-
3
4.3286
e
-
4
2.1742
e
-
4
1.5605
e
-
4
1.5743
e
-
4
1.5375
e
-
4
30
5.9706
e
-
3
2.6680
e
-
3
1.6015
e
-
4
1.5884
e
-
4
2.0451
e
-
4
1.3825
e
-
4Figure 3
Exact solutionu
(
x
,
0
) and numerical solution for ϵ
=
20
% with the final time T
=
25 and m
=
59, N
=
59, and t
0
=
1
e
-
10.Finally, we hope to use Example1 to show that the proposed algorithm is robust for order γ. For varying γ, we plot the relative error versus γ in Figure 4. The displayed results show that the numerical method is robust when γ is varying.Figure 4
The relative error versus orderγ for Example 1 with T
=
0.5, m
=
59, N
=
59, and t
0
=
1
e
-
10 for noise ϵ
=
1
%.Example 2.
Consider the following Neumann boundary value problem:(56)
0
D
t
1
/
2
u
(
x
,
t
)
=
∂
2
u
(
x
,
t
)
∂
x
2
,
x
∈
(
0
,
π
)
,
t
∈
(
0
,
T
)
,
u
(
x
,
0
)
=
cos
⁡
x
+
cos
⁡
2
x
,
x
∈
[
0
,
π
]
,
u
x
(
0
,
t
)
=
u
x
(
π
,
t
)
=
0
,
t
∈
[
0
,
T
]
.
The unique solution to (56) is given by
(57)
u
(
x
,
t
)
=
E
1
/
2
(
-
t
1
/
2
)
cos
⁡
x
+
E
1
/
2
(
-
4
t
1
/
2
)
cos
⁡
2
x
.
In this example, we show the numerical results for T
=
0.1 under the setting m
=
59, t
0
=
1
e
-
10, and N
=
55. Figure 5 shows that the proposed method is capable of giving satisfactory results for the case of the Neumann boundary condition.Figure 5
Exact solutionu
(
x
,
0
) and numerical solution for ϵ
∈
{
11
%
,
7
%
,
3
%
,
0.5
%
} with the final time T
=
0.1 and m
=
59, N
=
55, and t
0
=
1
e
-
10.Example 3.
We consider a two-dimensional fractional diffusion problem with the Dirichlet boundary value inΩ
=
(
0,1
)
×
(
0,1
):
(58)
0
D
t
γ
u
(
x
,
y
,
t
)
=
Δ
u
,
x
∈
Ω
,
t
∈
(
0
,
T
)
,
u
(
x
,
y
,
0
)
=
sin
π
x
sin
π
y
,
x
∈
Ω
-
,
u
(
x
,
y
,
t
)
=
0
,
x
∈
∂
Ω
,
t
∈
[
0
,
T
]
.
The initial distributionu
(
x
,
y
,
0
) is to be recovered by using the exact solution
(59)
u
(
x
,
y
,
t
)
=
E
γ
(
-
π
2
t
γ
)
sin
π
x
sin
π
y
.
We firstly display the choice of the set X
N in Figure 6. According to the analysis for 1-dimensional case about the number of J and m, we only deal with the case of t
0
=
1
e
-
10 and  m
=
100. The measurement points {
x
1
,
x
2
,
…
,
x
m
} are scattered in Ω, which are generated using the Matlab command rand. We report the relative error of f
α
,
N
,
ϵ for different final times T and noise levels ϵ in Table 3 for γ
=
1
/
2. The numerical comparison between exact solution and regularized solution is shown in Figure 7. The numerical results show that the proposed method is acceptable for the 2-dimensional example. We also consider the influence of varying γ on the numerical stability. The relative error versus γ is plotted in Figure 8, from which we can see that the proposed method is robust about parameter γ.Table 3
Example3: the relative errors with t
0
=
1
e
-
10 and m
=
100 for γ
=
1
/
2.
ϵ
T
1
2
5
10
1%
0.0131
0.0117
0.0138
0.0199
2%
0.0118
0.0141
0.0229
0.0201
5%
0.0272
0.0192
0.0304
0.0535Figure 6
The illustration of the setX
N.Figure 7
The comparison between the exact solution and regularized solution for Example3.Figure 8
The relative error versus orderγ for Example 3 with T
=
1,  m
=
100, and t
0
=
1
e
-
10 for noise ϵ
=
1
%.


## 5. Conclusion

In this paper, in a reproducing kernel Hilbert space setting, we propose a numerical reconstruction method, namely, the discretized Tikhonov regularization method, to recover the initial temperature distribution of the backward fractional diffusion problem. The implementation of the proposed method is simple and easy. Numerical tests show that the method is efficient.



---

*Source: 964373-2014-03-18.xml*