# An Iterative Regularization Method for Identifying the Source Term in a Second Order Differential Equation

Authors: Fairouz Zouyed; Sebti Djemoui
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper discusses the inverse problem of determining an unknown source in a second order differential equation from measured final data. This problem is ill-posed; that is, the solution (if it exists) does not depend continuously on the data. In order to solve the considered problem, an iterative method is proposed. Using this method a regularized solution is constructed and an a priori error estimate between the exact solution and its regularized approximation is obtained. Moreover, numerical results are presented to illustrate the accuracy and efficiency of this method.


---

## Body

## 1. Introduction

LetH be a separable Hilbert space with the inner product (
·
,
·
) and the norm ·. Consider the problem of finding the source term f
∈
H in the following system: (1)
u
′
′
t
+
2
A
u
′
t
+
A
2
u
t
=
f
,
0
<
t
<
T
,
u
0
=
0
,
u
′
0
=
0
,with the additional data(2)
u
T
=
g
,where A
:
D
(
A
)
⊂
H
→
H is a positive self-adjoint linear operator with a compact resolvent; we denote by σ
(
A
) the spectrum of the operator A.The problem (1) is an abstract version of the system(3)
u
t
t
x
,
t
-
2
Δ
u
t
x
,
t
+
Δ
2
u
x
,
t
=
f
x
,
0
<
t
<
T
,
x
∈
Ω
,
u
x
,
t
=
Δ
u
x
,
t
=
0
,
0
≤
t
≤
T
,
x
∈
∂
Ω
,
u
x
,
0
=
u
t
x
,
0
=
0
,
x
∈
Ω
,which arises in the mathematical study of structural damped vibrations of string or a beam [1–3]. Also this problem can be considered as a biparabolic problem in the abstract setting. For physical motivation we cite the biparabolic model proposed in [4] for more adequate mathematical description of heat and diffusion processes than the classical heat equation. For other models we refer the reader to [5–7].For most classical partial differential equations, the reconstruction of source functions from the final data or a partial boundary data is an inverse problem with many applications in several branches of sciences and engineering, such as geophysical prospecting and pollutant detection [8–12].The main difficulty of inverse source identification problems is that they are ill-posed, that is, even if a solution exists, it does not depend continuously on the data; in other words, small error in the data measurement can induce enormous error to the solution. Thus, special regularization methods that restore the stability with respect to measurements errors are needed. In the present work, we focus on an iterative method proposed by Kozlov and Maz’ya [13, 14] for solving the problem; it is based on solving a sequence of well-posed boundary value problems such that the sequence of solutions converges to the solution for the original problem. It has been successfully used for solving various classes of ill-posed elliptic, parabolic, and hyperbolic problems [5, 15–21].We note that although the interest in inverse problem has rapidly increased during this decade, the literature devoted to the class of problems (1) is quite scarce.The paper is organized as follows. Section2 gives some tools which are useful for this study; in Section 3 we introduce some basic results and we show the ill-posedness of the inverse problem; Section 4 gives a regularization solution and error estimation between the approximate solution and the exact one; the numerical implementation is described in Section 5 to illustrate the accuracy and efficiency of this method.


## 2. Preliminaries

Let(
φ
n
)
n
≥
1
⊂
H be an orthonormal eigenbasis corresponding to the eigenvalues (
λ
n
)
n
≥
1 such that (4)
A
φ
n
=
λ
n
φ
n
,
n
∈
N
∗
,
0
<
λ
1
≤
λ
2
⋯
≤
⋯
,
lim
n
→
∞
⁡
λ
n
=
+
∞
.
ξ
=
∑
n
=
1
∞
E
n
ξ
,
E
n
ξ
=
ξ
,
φ
n
φ
n
,
∀
ξ
∈
H
.We denote by {
T
(
t
)
=
e
-
t
A
}
t
≥
0 the analytic semigroup generated by -
A on H, (5)
T
t
ξ
=
∑
n
=
1
∞
e
-
λ
n
t
E
n
ξ
,
∀
ξ
∈
H
.For α
>
0, the space H
α is given by (6)
H
α
=
ξ
∈
H
:
∑
n
=
1
∞
1
+
λ
n
2
α
E
n
ξ
2
<
∞
,with the norm (7)
ξ
H
α
=
∑
n
=
1
∞
1
+
λ
n
2
α
E
n
ξ
2
1
/
2
,
ξ
∈
H
α
.We achieve this section by a result concerning nonexpansive operators.Definition 1.
A linear bounded operatorL
:
H
→
H is called nonexpansive if L
≤
1.LetL be an nonexpansive operator; to solve the equation(8)
I
-
L
φ
=
ψ
,we state a convergence theorem for a successive approximation method.Theorem 2 (see [22], p. 66).
LetL be a nonexpansive, self-adjoint positive operator on H. Let ψ
∈
H be such that (8) has a solution. If 1 is not eigenvalue of L, then the successive approximations (9)
φ
n
+
1
=
L
φ
n
+
ψ
,
n
=
0,1
,
2
,
…converge to a solution to (8) for any initial data φ
0
∈
H. Moreover, L
n
φ
→
0 for every φ
∈
H, as n
→
∞.


## 3. Basic Results

### 3.1. The Direct Problem

LetZ
=
D
(
A
)
×
H with the norm U
Z
2
=
A
ξ
1
2
+
ξ
2
2, U
=
ξ
1
ξ
2
∈
Z.For a givenf
∈
H, consider the direct problem(10)
w
′
′
t
+
2
A
w
′
t
+
A
2
w
t
=
f
,
0
<
t
<
T
,
w
0
=
0
,
w
′
0
=
0
.Making the change of variable w
′
=
v, we can write the second order equation in (10) as a first order system in the space Z as follows:(11)
z
′
t
=
A
z
t
+
F
,
0
<
t
<
T
,
z
0
=
0
,where z
=
w
v, F
=
0
f, and A
=
0
I
-
A
2
-
2
A.The linear operatorA is unbounded with the domain D
(
A
)
=
D
(
A
2
)
×
D
(
A
) and it is the infinitesimal generator of strongly continuous semigroup {
S
(
t
)
=
e
t
A
}
t
≥
0. Moreover {
S
(
t
)
}
t
≥
0 is analytic (see [1]) and it admits the following explicit form: (12)
S
t
U
=
∑
n
=
1
∞
e
t
B
n
P
n
U
,
U
=
ξ
1
ξ
2
∈
Z
,where B
n
=
0
1
-
λ
n
2
-
2
λ
n and {
P
n
}
n
≥
1 is a complete family of orthogonal projections in Z given by P
n
=
d
i
a
g
(
E
n
,
E
n
).Using matrix algebra, we obtain(13)
e
t
B
n
=
e
-
λ
n
t
+
λ
n
t
e
-
λ
n
t
t
e
-
λ
n
t
-
λ
n
2
t
e
-
λ
n
t
-
λ
n
t
e
-
λ
n
t
+
e
-
λ
n
t
.From the semigroup theory (see [23]), the problem (11) admits a unique solution z
∈
C
(
[
0
,
T
)
,
Z
) given by (14)
z
=
∫
0
t
S
t
-
s
F
d
s
.Hence, (15)
z
=
∫
0
t
∑
n
=
1
∞
e
t
-
s
B
n
P
n
F
d
s
=
∫
0
t
∑
n
=
1
∞
σ
n
1
t
,
s
σ
n
2
t
,
s
σ
n
3
t
,
s
σ
n
4
t
,
s
·
0
f
,
φ
n
φ
n
d
s
,such that (16)
σ
n
1
t
,
s
=
e
-
λ
n
t
-
s
+
λ
n
t
-
s
e
-
λ
n
t
-
s
,
σ
n
2
t
,
s
=
t
-
s
e
-
λ
n
t
-
s
,
σ
n
3
t
,
s
=
-
λ
n
2
t
-
s
e
-
λ
n
t
-
s
,
σ
n
4
t
,
s
=
-
λ
n
t
-
s
e
-
λ
n
t
-
s
+
e
-
λ
n
t
-
s
.As a consequence, we obtain the following theorem.Theorem 3.
The problem (10) admits a unique solution w
∈
C
(
[
0
,
T
)
,
D
(
A
)
)
∩
C
1
(
[
0
,
T
)
,
H
) given by(17)
w
t
=
K
t
f
=
A
-
2
I
-
I
+
t
A
e
-
t
A
f
=
∑
n
=
1
∞
1
-
1
+
t
λ
n
e
-
t
λ
n
λ
n
2
f
,
φ
n
φ
n
.


### 3.2. Ill-Posedness of the Inverse Problem

Now, we wish to solve the inverse problem, that is, find the source termf in the system (1). Making use of the supplementary condition (2) and defining the operator K
(
T
)
:
f
→
g, we have (18)
g
=
u
T
=
K
T
f
=
∑
n
=
1
∞
σ
n
E
n
f
,where σ
n
=
(
1
-
(
1
+
T
λ
n
)
e
-
T
λ
n
)
/
λ
n
2.It is easy to see thatK
(
T
) is a self-adjoint compact linear operator. On the other hand, (19)
g
=
∑
n
=
1
∞
E
n
g
=
∑
n
=
1
∞
σ
n
E
n
f
,so(20)
σ
n
E
n
f
=
E
n
g
,which implies (21)
E
n
f
=
1
σ
n
E
n
g
,and therefore (22)
f
=
K
T
-
1
g
=
∑
n
=
1
∞
1
σ
n
E
n
g
.Note that 1
/
σ
n
→
∞ as n
→
∞, so the inverse problem is ill-posed; that is, the solution does not depend continuously on the given data. Hence this problem cannot be solved by using classical numerical methods.Remark 4.
As many boundary inverse value problems for partial differential equations which are ill-posed, the study of the problem (1) is reduced to the study of the equation K
(
T
)
f
=
g, where K
(
T
) is a compact self-adjoint operator in the Hilbert space H. This equation can be rewritten in the following way: (23)
f
=
I
-
γ
K
T
f
+
γ
g
=
L
f
+
γ
g
,where γ is a positive number satisfying γ
<
1
/
K
(
T
).
In the next section, we will show that the operatorL is nonexpansive and 1 is not eigenvalue of L, so it follows from Theorem 2 that (
f
n
)
n
∈
N
∗ converges and (
I
-
γ
K
(
T
)
)
n
f
→
0, for every f
∈
H, as n
→
∞.


## 3.1. The Direct Problem

LetZ
=
D
(
A
)
×
H with the norm U
Z
2
=
A
ξ
1
2
+
ξ
2
2, U
=
ξ
1
ξ
2
∈
Z.For a givenf
∈
H, consider the direct problem(10)
w
′
′
t
+
2
A
w
′
t
+
A
2
w
t
=
f
,
0
<
t
<
T
,
w
0
=
0
,
w
′
0
=
0
.Making the change of variable w
′
=
v, we can write the second order equation in (10) as a first order system in the space Z as follows:(11)
z
′
t
=
A
z
t
+
F
,
0
<
t
<
T
,
z
0
=
0
,where z
=
w
v, F
=
0
f, and A
=
0
I
-
A
2
-
2
A.The linear operatorA is unbounded with the domain D
(
A
)
=
D
(
A
2
)
×
D
(
A
) and it is the infinitesimal generator of strongly continuous semigroup {
S
(
t
)
=
e
t
A
}
t
≥
0. Moreover {
S
(
t
)
}
t
≥
0 is analytic (see [1]) and it admits the following explicit form: (12)
S
t
U
=
∑
n
=
1
∞
e
t
B
n
P
n
U
,
U
=
ξ
1
ξ
2
∈
Z
,where B
n
=
0
1
-
λ
n
2
-
2
λ
n and {
P
n
}
n
≥
1 is a complete family of orthogonal projections in Z given by P
n
=
d
i
a
g
(
E
n
,
E
n
).Using matrix algebra, we obtain(13)
e
t
B
n
=
e
-
λ
n
t
+
λ
n
t
e
-
λ
n
t
t
e
-
λ
n
t
-
λ
n
2
t
e
-
λ
n
t
-
λ
n
t
e
-
λ
n
t
+
e
-
λ
n
t
.From the semigroup theory (see [23]), the problem (11) admits a unique solution z
∈
C
(
[
0
,
T
)
,
Z
) given by (14)
z
=
∫
0
t
S
t
-
s
F
d
s
.Hence, (15)
z
=
∫
0
t
∑
n
=
1
∞
e
t
-
s
B
n
P
n
F
d
s
=
∫
0
t
∑
n
=
1
∞
σ
n
1
t
,
s
σ
n
2
t
,
s
σ
n
3
t
,
s
σ
n
4
t
,
s
·
0
f
,
φ
n
φ
n
d
s
,such that (16)
σ
n
1
t
,
s
=
e
-
λ
n
t
-
s
+
λ
n
t
-
s
e
-
λ
n
t
-
s
,
σ
n
2
t
,
s
=
t
-
s
e
-
λ
n
t
-
s
,
σ
n
3
t
,
s
=
-
λ
n
2
t
-
s
e
-
λ
n
t
-
s
,
σ
n
4
t
,
s
=
-
λ
n
t
-
s
e
-
λ
n
t
-
s
+
e
-
λ
n
t
-
s
.As a consequence, we obtain the following theorem.Theorem 3.
The problem (10) admits a unique solution w
∈
C
(
[
0
,
T
)
,
D
(
A
)
)
∩
C
1
(
[
0
,
T
)
,
H
) given by(17)
w
t
=
K
t
f
=
A
-
2
I
-
I
+
t
A
e
-
t
A
f
=
∑
n
=
1
∞
1
-
1
+
t
λ
n
e
-
t
λ
n
λ
n
2
f
,
φ
n
φ
n
.


## 3.2. Ill-Posedness of the Inverse Problem

Now, we wish to solve the inverse problem, that is, find the source termf in the system (1). Making use of the supplementary condition (2) and defining the operator K
(
T
)
:
f
→
g, we have (18)
g
=
u
T
=
K
T
f
=
∑
n
=
1
∞
σ
n
E
n
f
,where σ
n
=
(
1
-
(
1
+
T
λ
n
)
e
-
T
λ
n
)
/
λ
n
2.It is easy to see thatK
(
T
) is a self-adjoint compact linear operator. On the other hand, (19)
g
=
∑
n
=
1
∞
E
n
g
=
∑
n
=
1
∞
σ
n
E
n
f
,so(20)
σ
n
E
n
f
=
E
n
g
,which implies (21)
E
n
f
=
1
σ
n
E
n
g
,and therefore (22)
f
=
K
T
-
1
g
=
∑
n
=
1
∞
1
σ
n
E
n
g
.Note that 1
/
σ
n
→
∞ as n
→
∞, so the inverse problem is ill-posed; that is, the solution does not depend continuously on the given data. Hence this problem cannot be solved by using classical numerical methods.Remark 4.
As many boundary inverse value problems for partial differential equations which are ill-posed, the study of the problem (1) is reduced to the study of the equation K
(
T
)
f
=
g, where K
(
T
) is a compact self-adjoint operator in the Hilbert space H. This equation can be rewritten in the following way: (23)
f
=
I
-
γ
K
T
f
+
γ
g
=
L
f
+
γ
g
,where γ is a positive number satisfying γ
<
1
/
K
(
T
).
In the next section, we will show that the operatorL is nonexpansive and 1 is not eigenvalue of L, so it follows from Theorem 2 that (
f
n
)
n
∈
N
∗ converges and (
I
-
γ
K
(
T
)
)
n
f
→
0, for every f
∈
H, as n
→
∞.


## 4. Iterative Procedure and Convergence Results

The alternating iterative method is based on reducing the ill-posed problem (1) to a sequence of well-posed boundary value problems and consists of the following steps.First, we start by lettingf
0
∈
H be arbitrary; the initial approximation u
0 is the solution to the direct problem (24)
u
0
′
′
+
2
A
u
0
′
+
A
2
u
0
=
f
0
,
0
<
t
<
T
,
u
0
0
=
0
,
u
0
′
0
=
0
.Then, if the pair (
f
k
,
u
k
) has been constructed, let(25)
f
k
+
1
=
f
k
-
γ
u
k
T
-
g
,where γ is such that (26)
0
<
γ
<
1
K
T
,and K
(
T
)
=
sup
n
∈
N
∗
⁡
(
1
-
(
1
+
T
λ
n
)
e
-
λ
n
T
)
/
λ
n
2.Finally, we getu
k
+
1 by solving the problem (27)
u
k
+
1
′
′
+
2
A
u
k
+
1
′
+
A
2
u
k
+
1
=
f
k
+
1
,
0
<
t
<
T
,
u
k
+
1
0
=
0
,
u
k
+
1
′
0
=
0
.Let us iterate backwards in (25) to obtain(28)
f
k
+
1
=
f
k
-
γ
K
T
f
k
+
γ
g
=
I
-
γ
K
T
f
k
+
γ
g
=
I
-
γ
K
T
k
+
1
f
0
+
γ
∑
j
=
0
k
I
-
γ
K
T
j
g
.Now, we introduce some properties and tools which are useful for our main theorems.Lemma 5.
The norm of the operatorK
(
t
) is given by (29)
K
t
=
sup
n
∈
N
∗
⁡
1
-
1
+
t
λ
n
e
-
λ
n
t
λ
n
2
=
1
-
1
+
t
λ
1
e
-
λ
1
t
λ
1
2
.Proof.
We aim to find the supremum of the function(
1
-
(
1
+
t
λ
n
)
e
-
λ
n
t
)
/
λ
n
2, n
∈
N
∗, and for this purpose, fix t, let μ
=
λ
t, and define the function (30)
G
1
μ
=
1
-
1
+
μ
e
-
μ
μ
2
,
for
μ
≥
μ
1
=
λ
1
t
.We compute (31)
G
1
′
μ
=
μ
2
+
2
μ
+
2
e
-
μ
-
2
μ
3
.Put (32)
h
μ
=
μ
2
+
2
μ
+
2
e
-
μ
-
2
.Hence,(33)
G
1
′
μ
=
h
μ
μ
3
.To study the monotony of G
1, it suffices to determine the sign of h. We have (34)
h
′
μ
=
-
μ
2
e
-
μ
<
0
,
∀
μ
>
0
,and then h is decreasing; moreover h
(
μ
)
⊂
]
-
2,0
[, ∀
μ
>
0. Hence G
1
′
(
μ
)
<
0, ∀
μ
≥
μ
1, which implies that G
1 is decreasing and (35)
sup
μ
≥
μ
1
⁡
G
1
μ
=
G
1
μ
1
.Therefore, (36)
sup
n
≥
1
⁡
1
-
1
+
λ
n
t
e
-
λ
n
t
λ
n
2
=
1
-
1
+
λ
1
t
e
-
λ
1
t
λ
1
2
.Proposition 6.
For the linear operatorL
=
I
-
γ
K
(
T
), one has the following properties: (1)
Lis positive and self-adjoint,
(2)
L is nonexpansive,
(3)
1 is not an eigenvalue of L.Proof.
Form properties of operatorA and the definition of L it follows that L is self-adjoint and nonexpansive positive operator and from the inequality (37)
0
<
1
-
γ
1
-
1
+
T
λ
e
-
λ
T
λ
2
<
1
,
for
λ
∈
σ
A
,it follows that the point spectrum of L
,
σ
p
(
L
)
⊂
]
0,1
[. Then 1 is not eigenvalue of the operator L.Lemma 7.
Ifλ
>
0, one has the estimates(38)
1
1
+
λ
2
≤
max
⁡
3
T
2
,
1
1
-
1
+
T
λ
e
-
λ
T
λ
2
,
(39)
0
<
1
-
1
+
t
λ
e
-
λ
t
λ
2
<
T
2
,
∀
t
∈
0
,
T
.Proof.
To establish (38), let us first prove that(40)
1
3
+
μ
2
≤
1
-
1
+
μ
e
-
μ
μ
2
,
∀
μ
>
0
,which is equivalent to prove that (41)
G
2
μ
=
3
-
3
+
μ
2
1
+
μ
e
-
μ
≥
0
,
∀
μ
>
0
.We have (42)
G
2
′
μ
=
μ
μ
-
1
2
e
-
μ
≥
0
,
∀
μ
>
0
.Then, G
2 is nondecreasing and it follows that G
2
(
μ
)
⊂
]
0,3
[. SoG
2
(
μ
)
≥
0, ∀
μ
>
0.
Choosingμ
=
T
λ in (40), we obtain (43)
1
3
+
T
λ
2
≤
1
-
T
λ
+
1
e
-
T
λ
T
λ
2
.So,(44)
T
2
max
⁡
3
,
T
2
1
+
λ
2
≤
1
-
1
+
T
λ
e
-
T
λ
λ
2
.From (44), we deduce (38).
Now, we prove the estimate (39). It is easy to verify that (45)
G
3
μ
=
1
-
1
+
μ
e
-
μ
-
μ
2
<
0
,
∀
μ
>
0
.Then, if we choose μ
=
t
λ, we get(46)
1
-
1
+
t
λ
e
-
t
λ
<
t
2
λ
2
,
∀
λ
>
0
,
∀
t
∈
0
,
T
.Hence, from (46), (39) follows.Theorem 8.
Letu be a solution to the inverse problem (1). Let f
0
∈
H be an arbitrary initial data element for the iterative procedure proposed above and let u
k be the kth approximate solution. Then (i)
The method converges; that is,(47)
sup
t
∈
0
,
T
⁡
u
k
t
-
u
t
⟶
0
,
as
k
⟶
∞
.
(ii)
Moreover, if, for someα
=
1
+
θ, θ
>
0, f
0
-
f
∈
H
α, that is, f
0
-
f
H
α
≤
E, then the rate of convergence of the method is given by (48)
sup
t
∈
0
,
T
⁡
u
k
t
-
u
t
≤
T
2
C
E
k
-
α
/
2
,where C is a positive constant independent of k.Proof.
(i) From (28), we get (49)
f
k
=
I
-
γ
K
T
k
f
0
+
I
-
I
-
γ
K
T
k
K
T
-
1
g
,and then(50)
f
k
=
I
-
γ
K
T
k
f
0
-
f
+
f
,which implies that (51)
u
k
t
-
u
t
=
K
t
f
k
-
f
=
K
t
I
-
γ
K
T
k
f
0
-
f
.Hence,(52)
u
k
t
-
u
t
≤
K
t
I
-
γ
K
T
k
f
0
-
f
.From Lemma 5 and (39) we have(53)
sup
t
∈
0
,
T
⁡
K
t
=
sup
t
∈
0
,
T
⁡
1
-
1
+
t
λ
1
e
-
t
λ
1
λ
1
2
<
T
2
.Combining (52) and (53) and passing to the supremum with respect to t
∈
[
0
,
T
], we obtain (54)
sup
t
∈
0
,
T
⁡
u
k
t
-
u
t
≤
T
2
I
-
γ
K
T
k
f
0
-
f
⟶
0
,
as
k
⟶
∞
.
(ii) By part (i), we have(55)
u
k
t
-
u
t
2
≤
T
4
∑
n
=
1
∞
1
-
γ
1
-
1
+
λ
n
T
e
-
λ
n
T
λ
n
2
2
k
f
0
-
f
,
φ
n
2
,and hence (56)
u
k
t
-
u
t
2
≤
T
4
∑
n
=
1
∞
1
-
γ
1
-
1
+
λ
n
T
e
-
λ
n
T
λ
n
2
2
k
1
+
λ
n
2
-
α
1
+
λ
n
2
α
f
0
-
f
,
φ
n
2
.Using the inequality (38), we obtain (57)
u
k
t
-
u
t
2
≤
T
4
max
⁡
3
T
2
,
1
α
∑
n
=
1
∞
1
-
γ
β
n
2
k
β
n
α
1
+
λ
n
2
α
f
0
-
f
,
φ
n
2
,where β
n
=
(
(
1
-
(
1
+
λ
n
T
)
e
-
λ
n
T
)
/
λ
n
2
)
.
So, it follows that(58)
u
k
t
-
u
t
2
≤
T
4
max
⁡
3
T
2
,
1
α
sup
0
≤
β
n
≤
T
2
⁡
1
-
γ
β
n
2
k
β
n
α
f
0
-
f
H
α
2
.Put (59)
ϕ
β
=
1
-
γ
β
2
k
β
α
,
0
≤
β
≤
T
2
.We compute (60)
ϕ
′
β
=
1
-
γ
β
2
k
-
1
β
α
-
1
-
γ
2
k
+
α
β
+
α
.Setting ϕ
′
(
β
)
=
0, it follows that β
∗
=
α
/
(
2
k
+
α
)
γ is the critical point of ϕ. It is easy to see that the maximum of ϕ is attained at β
∗. So (61)
sup
0
≤
β
≤
T
2
⁡
ϕ
β
≤
ϕ
β
∗
=
1
-
γ
β
∗
2
k
β
∗
α
≤
β
∗
α
=
α
2
k
+
α
γ
α
,and hence(62)
sup
0
≤
β
≤
T
2
⁡
ϕ
β
≤
α
2
γ
α
k
-
α
.Combining (58) and (62), we obtain (63)
sup
t
∈
0
,
T
⁡
u
k
t
-
u
t
2
≤
T
4
α
2
γ
max
⁡
3
T
2
,
1
α
1
k
α
E
2
.
Since in practice the measured datag is never known exactly but only up to an error of, say, δ
>
0, it is our aim to solve the equation K
(
T
)
f
=
g from the knowledge of a perturbed right-hand side g
δ satisfying(64)
g
-
g
δ
<
δ
,where δ
>
0 denotes a noise level. In the following theorem, we consider the case of inexact data.Theorem 9.
Letα
=
1
+
θ, (
θ
>
0
), f
0 be an arbitrary initial data element for the iterative procedure proposed above such that (
f
0
-
f
)
∈
H
α, let u
k be the kth approximations solution for the exact data g, and let u
k
δ be the kth approximations solution corresponding to the perturbed data g
δ such that (64) holds. Then one has the following estimate:(65)
sup
t
∈
0
,
T
⁡
u
k
t
-
u
t
≤
T
2
δ
γ
k
+
C
E
1
k
α
/
2
.Proof.
Let(66)
f
k
=
I
-
γ
K
T
k
f
0
+
γ
∑
j
=
0
k
-
1
I
-
γ
K
T
j
g
,
u
k
t
=
K
t
f
k
,
f
k
δ
=
I
-
γ
K
T
k
f
0
+
γ
∑
j
=
0
k
-
1
I
-
γ
K
T
j
g
δ
,
u
k
δ
t
=
K
t
f
k
δ
.Using the triangle inequality, we obtain (67)
u
k
δ
-
u
≤
u
k
δ
-
u
k
+
u
k
-
u
.From Theorem 8, we have(68)
sup
t
∈
0
,
T
⁡
u
k
t
-
u
t
≤
T
2
C
E
1
k
α
/
2
.On the other hand, (69)
u
k
δ
t
-
u
k
t
=
K
t
f
k
δ
-
f
k
≤
T
2
γ
∑
j
=
0
k
-
1
I
-
γ
K
T
j
g
δ
-
g
≤
T
2
δ
γ
∑
j
=
0
k
-
1
I
-
γ
K
T
j
≤
T
2
δ
γ
∑
j
=
0
k
-
1
I
-
γ
K
T
j
.Since (70)
I
-
γ
K
T
≤
1
,it follows that(71)
sup
t
∈
0
,
T
⁡
u
k
δ
t
-
u
k
t
≤
T
2
δ
γ
k
.Combining (68) and (71) and passing to the supremum with respect to t
∈
[
0
,
T
], we obtain the estimate (65).Remark 10.
If we choose the number of the iterationsk
(
δ
) so that k
(
δ
)
→
0 as δ
→
0, we obtain (72)
sup
t
∈
0
,
T
⁡
u
k
δ
t
-
u
t
⟶
0
,
as
k
⟶
+
∞
.


## 5. Numerical Implementation

In this section, an example is devised for verifying the effectiveness of the proposed method. Consider the problem of finding a pair of functions(
u
(
x
,
t
)
,
f
(
x
)
), in the system(73)
∂
2
∂
t
2
u
x
,
t
-
2
∂
2
∂
x
2
∂
∂
t
u
x
,
t
+
∂
4
∂
x
4
u
x
,
t
=
f
x
,
t
,
x
∈
0,1
×
0,1
,
u
0
,
t
=
u
1
,
t
=
0
,
t
∈
0,1
,
u
x
,
0
=
u
t
x
,
0
=
0
,
x
∈
0,1
,
u
x
,
1
=
g
x
,
x
∈
0,1
.Denote(74)
A
=
-
∂
2
∂
x
2
,
with
D
A
=
H
0
1
0,1
∩
H
2
0,1
⊂
H
=
L
2
0,1
.
λ
n
=
n
2
π
2
,
φ
n
=
2
sin
⁡
n
π
x
,
n
=
1,2
,
…are eigenvalues and orthonormal eigenfunctions, which form a basis for H.The solution of the above problem is given by(75)
u
x
,
t
=
∑
n
=
1
∞
1
-
1
+
n
π
2
t
e
-
n
π
2
t
n
π
4
f
n
φ
n
,where f
n
=
(
f
,
φ
n
)
=
2
∫
0
1
‍
f
(
s
)
s
i
n
(
n
π
s
)
d
s, n
=
1,2
,
…Now, to solve the inverse problem, making use of the supplementary condition and defining the operatorK
:
f
→
g, we have(76)
g
x
=
u
x
,
1
=
K
f
x
=
2
∑
n
=
1
∞
1
-
1
+
n
π
2
e
-
n
π
2
n
π
4
∫
0
1
f
s
sin
⁡
n
π
s
d
s
sin
⁡
n
π
x
.Example 11.
In the following, we first selected the exact solutionf
(
x
) and obtained the exact data function g
(
x
) through solving the forward problem. Then we added a normally distributed perturbation to each data function and obtained vectors g
δ
(
x
). Finally we obtained the regularization solutions through solving the inverse problem with noisy data g
δ
(
x
) satisfying(77)
g
-
g
δ
L
2
0,1
2
≤
δ
.It is easy to see that if f
(
x
)
=
s
i
n
π
x, then (78)
u
x
,
t
=
1
-
1
+
π
2
t
e
-
π
2
t
π
4
sin
⁡
π
xis the exact solution of the problem (73). Consequently, g
(
x
)
=
(
(
1
-
(
1
+
π
2
)
e
-
π
2
)
/
π
4
)
s
i
n
(
π
x
)
.
Now, we propose to approximate the first and second space derivatives by using central difference and we consider an equidistant grid points to a spatial step sizex
0
=
0
<
x
1
<
⋯
<
x
N
+
1
=
1, (h
=
1
/
(
N
+
1
)), where N is a positive integer. We get the following semidiscrete problem:(79)
u
′
′
x
i
,
t
+
2
A
h
u
′
x
i
,
t
+
A
h
2
u
x
i
,
t
=
f
x
i
,
x
i
=
i
h
,
i
=
1
,
…
,
N
,
0
<
t
<
1
,
u
0
,
t
=
u
1
,
t
=
0
,
0
<
t
<
1
,
u
x
i
,
0
=
u
′
x
i
,
0
=
0
,
x
i
=
i
h
,
i
=
1
,
…
,
N
,
u
x
i
,
1
=
g
x
i
,
x
i
=
i
h
,
i
=
1
,
…
,
N
,where A
h is the discretisation matrix stemming from the operator A
=
-
d
2
/
d
x
2, and(80)
A
h
=
1
h
2
Tridiag
⁡
-
1,2
,
-
1is a symmetric, positive definite matrix, with eigenvalues (81)
μ
j
=
4
N
+
1
2
sin
2
⁡
j
π
2
N
+
1
,
j
=
1
,
…
,
N
,and orthonormal eigenvalues (82)
v
j
=
sin
⁡
m
j
π
N
+
1
1
≤
m
≤
N
,
j
=
1
,
…
,
N
.We assume that it is fine enough so that the discretization errors are small compared to the uncertainty δ of the data; this means that A
h is a good approximation of the differential operator A whose unboundedness is reflected in a large norm of A
h (see [24]).
Adding a random distributed perturbation to each data function, we obtain(83)
g
δ
=
g
+
ε
randn
⁡
size
g
,where ε indicates the noise level of the measurements data and the function randn(·) generates arrays of random numbers whose elements are normally distributed with mean 0, variance σ
2
=
1, and standard deviation σ
=
1. r
a
n
d
n(size(g)) returns an array of random entries that is of the same size as g. The noise level δ can be measured in the sense of root mean square error (RMSE) according to (84)
δ
=
g
δ
-
g
l
2
=
1
N
+
1
∑
i
=
0
N
g
x
i
-
g
δ
x
i
2
1
/
2
.The relative error is given as follows: (85)
RE
⁡
f
=
f
approximate
-
f
exact
l
2
f
exact
l
2
.The discrete iterative approximation of (66) is given by (86)
f
k
δ
x
i
=
I
-
γ
K
h
k
f
0
x
i
+
γ
∑
j
=
0
k
-
1
I
-
γ
K
h
j
g
δ
x
i
,
i
=
1
,
…
,
N
,where K
h
=
A
h
-
2
(
I
N
-
(
I
N
+
A
h
)
e
-
A
h
) and γ
<
1
/
K
h
=
(
μ
1
2
/
(
1
-
(
1
+
μ
1
)
e
-
μ
1
)
).Figures1–4 display that as the amount of noise ε decreases, the regularized solutions approximate better the exact solution.Figure 1
The comparison between the exact solutionf
e and its computed approximations f
a for N
=
60  k
=
4 and noisy level ε
=
1
0
-
3.Figure 2
The comparison between the exact solutionf
e and its computed approximations f
a for N
=
60  k
=
4 and noisy level ε
=
1
0
-
4.Figure 3
The comparison between the exact solutionf
e and its computed approximations f
a forN
=
60  k
=
5 and noisy level ε
=
1
0
-
3.Figure 4
The comparison between the exact solutionf
e and its computed approximations f
a for N
=
60  k
=
5 and noisy level ε
=
1
0
-
4
.Table1 shows that for k
=
4 or k
=
5 the relative error decreases with the decease of epsilon which is consistent with our regularization.Table 1
Relative error RE(f).
N
k
ε
RE(f)
60
4
10−3
0.2039
60
4
10−4
0.0945
60
5
10−3
0.3032
60
5
10−4
0.0305


## 6. Conclusion

In this paper, we have extended the iterative method to identify the unknown source term in a second order differential equation, convergence results were established, and error estimates have been obtained under an a priori bound of the exact solution. Some numerical tests have been given to verify the validity of the method.



---

*Source: 713403-2015-10-12.xml*