Inverse Problem of Finding the Coefficient of the Lowest Term in Two-dimensional Heat Equation with Ionkin-type Boundary Condition

0.84 MB·

by  Mansur I. Ismailov

Inverse Problem of Finding the Coefficient of the
Lowest Term in Two-dimensional Heat Equation
with Ionkin-type Boundary Condition
Mansur I. Ismailov*,**, Sait Erkovan*
*GebzeTechnicalUniversity,DepartmentofMathematics,41400,Gebze/Kocaeli,Turkey,
**InstituteofMathematicsandMechanics,NationalAcademyofSciencesofAzerbaijan,AZ1141,Baku,Azerbaijan,
[email protected],[email protected]
7
1
0 Abstract
2
We consider an inverse problem of determining the time-dependent lowest order coefficient of

b two-dimensional(2D)heatequationwithIonkinboundaryandtotalenergyintegraloverdetermi-
e nation condition. The well-posedness of the problem is obtained by generalized Fourier method
F
combined by the Banach fixed poind theorem. For obtaining a numerical solution of the inverse
3 problem, we propose the discretization method from a new combination. On the one hand, it
is known the traditional method of uniform finite difference combined with numerical integra-

] tion on a uniform grid (trapezoidal and Simpson’s), on the other hand, we give the method of
P
non-uniformfinitedifferenceiscombinedbyanumericalintegrationonanon-uniformgrid(with
A
Gauss-Lobatto nodes). Numerical examples illustrate how to implement the method.
.
h
Keywords: 2D heat equation; Ionkin-type boundary condition; Generalized Fourier method;
t
a Uniform finite diference method; Non-uniform finite difference method; Numerical integration
m
[ MSC 2010: Primary 35R30, Secondary 35K20, 65D30, 65M06

3
v 1. Introduction and Problem Formulation
4
3
0 Many researchers have been involved in the inverse coefficient problem for heat equation since this
9 topic is of importance in some engineering texts and many industrial applications. The 1D problems
0 havebeenwellinvestigated. However,therearesomanyworkswhichareneededtoinvestigateforthe
.
1 multidimentionalproblems(seethebooks[1,2,3]andreferencestherein), althoughmanyresearchers
0 have reported its difficulties.
7 Among these inverse problems, much attention is given to the determination of the lowest order
1
coefficient in heat equation, in particular, when this coefficient depends solely on time. Various
:
v methodsforfindingthelowestordercoefficientinamoregeneralmultidimensionalparabolicequation
i
X havebeenaddressedinnumerousworks,see[4,5,6,7,8]fortimedependentcoefficient,[9,10,11,12,
13, 14, 15] for space dependent coefficient, [16, 17, 18] for both time and space dependent coefficient.
r
a Theboundaryconditionsaremostfrequentlyclassical(Dirichlet,NeumannandRobin)andadditional
condition is most frequently specified as the solution at an interior point or an integral mean over the
entire domain.
In this paper, we consider the problem of determining the lowest coefficient that depends on time
only, for a two-dimensional parabolic equation with Ionkin type nonlocal boundary condition and the
totalenergymeasurement. Thistypeofnonlocalboundaryconditionsforpartialdifferentialequations
were considered by numerous authors starting from the classical work [19]. Nonlocal problems have
the specific feature that the corresponding spatial differential operator is nonself-adjoint and hence
the system of eigenfunctions is not complete and must be supplemented by associated functions. The
existenceanduniquenessofthesolutionofsuchaninverseproblemandwell-posednessofthisproblem
are examined by using the method of series expansion in terms of eigenfunctions and associated
functions of mentioned spatial differential operator.
1
Numerical methods for solving the problem of the identification of the lowest coefficient of 1D
parabolicequationsareconsideredinmanyworks. Theworks[20,21,22,23]dealingwiththenumer-
icalsolutionofinverseproblemsforparabolicequationsforanapplicationpointofview. Wehighlight
separately studies dealing with numerical solving inverse problems of finding time-dependent lowest
coefficientformultidimensionalparabolicequations[24,25,26]inviewofthepracticaluse. Discretiza-
tion in our problem is performed using uniform and non-uniform finite difference method combined
with suitable nodes for numerical integration on uniform and non-uniform grids.
We study the following problem for two-dimensional heat equation with Ionkin-type boundary
condition:
Let D ={(x,y):0<x,y <1}. In the domain D =D ×(0<t≤T), we consider
xy xy
∂u ∂2u ∂2u
= + −p(t)u(x,y,t)+f(x,y,t) (1.1)
∂t ∂x2 ∂y2
with the initial condition
u(x,y,0)=ϕ(x,y), 0≤x,y ≤1, (1.2)
nonlocal boundary conditions
u(0,y,t)=u(1,y,t), u (1,y,t)=0, 0≤y ≤1, 0≤t≤T,
x
(1.3)
u(x,0,t)=u(x,1,t)=0, 0≤x≤1, 0≤t≤T,
and overdetermination condition
1 1
(cid:90) (cid:90)
u(x,y,t)dxdy =E(t). (1.4)
0 0
The problem of finding a pair {p(t),u(x,y,t)} ∈ C[0,T]×C2,2,1(D) in (1.1)-(1.4) will be called an
inverse problem.
The paper is organized as follows. In Section 2, we recall some necessary results on basisness of
root functions concerning to two-dimensional spectral problem with Ionkin-type boundary condition,
on finite difference discretization on uniform and non-uniform grids and on some product rules of
numerical integration. The well-posedness of inverse problem (1.1)-(1.4) for small T is showed by
using generalized Fourier method combined with Banach fixed point theorem in Section 3. The
numerical methods for solving the inverse problem (1.1)-(1.4) by applying uniform and non-uniform
finite difference with suitable numerical integration rules are given in Section 4. We also present
several numerical examples intended to illustrate the behaviour of the proposed methods. The tests
were performed using MATLAB are discussed in Section 4. Finally, the concluding remark with the
comparasion of two different numerical methods on uniform and non-uniform grids are presented in
Section 5.
2. Preliminaries
In this section, we recall some necessary results on basisness of root functions concerning to two-
dimensional spectral problem with Ionkin-type boundary conditions (see [27]), on finite difference
discretizationonuniformandnon-uniformgrids(see[28,29])andonsomeproductrulesofnumerical
integration which are used for solving the problem numerically.
2.1. Spectral problem
Consider the following spectral problem:
∂2Z ∂2Z
+ +µZ =0, 0<x,y <1, (2.1)
∂x2 ∂y2
∂Z(1,y)
Z(0,y)=Z(1,y), =0, Z(x,0)=Z(x,1)=0, 0≤x,y ≤1, (2.2)
∂x
2
where µ is the separation parameter. We present the solution in the form
Z(x,y)=X(x)V(y). (2.3)
Substituting this expression into (2.1) and (2.2), we obtain the following problems
V(cid:48)(cid:48)(y)+λV(y)=0, 0<y <1, V(0)=V(1)=0, (2.4)
X(cid:48)(cid:48)(x)+γX(x)=0, 0<x<1, X(0)=X(1), X(cid:48)(1)=0, (2.5)
where γ =µ−λ. It is known the solutions of problem (2.4) have the form
√
λ =(πk)2, V (y)= 2sin(πky), k =1,2,....
k k
Here and in the following, we give the constants multiplying the eigenfunctions and associated func-
tions from normalization conditions.
The eigenvalues and the corresponding eigenfunctions of problem (2.5) were mentioned in [30]:
γ =(2πm)2, m=0,1,2,..., X =2, X (x)=4cos(2πmx), m=1,2,....
m 0 m
Consequently, the eigenvalues and eigenfunctions of problem (2.1), (2.2) with the representation (2.3)
have the form
µ =γ +λ =(2πm)2+(πk)2, Z (x,y)=X (x)V (y), m=0,1,2,..., k =1,2,....
m,k m k m,k m k
Here, the set of eigenfunctions Z (x,y) is not complete in the space L (D ); therefore, it is sup-
m,k 2 xy
∼
plemented the set of eigenfunctions by the following associated functions Z (x,y), m,k =1,2... in
m,k
[27] by following the lines of [31]
∼ √
Z =4(1−x)sin(2πmx) 2sin(πky).
m,k
We denote the system of root functions of problem (2.1), (2.2) as follows [27]:
∼
Z (x,y)=X (x)V (y), Z (x,y)=X (x)V (y), Z (x,y)=Z , k,m=1,2,...
0,k 0 k 2m−1,k m k 2m,k m,k
(2.6)
The adjoint problem of (2.1), (2.2) is
∂2W/∂x2+∂2W/∂y2+µW =0, 0<x,y <1,
(2.7)
W(0,y)=0, ∂W(0,y)/∂x=∂W(1,y)/∂x, W(x,0)=W(x,1)=0, 0≤x,y ≤1.
The root functions of the adjoint problem (2.7) have the form
W (x,y)=xV (y), W (x,y)=sin(2πmx)V (y),
0,k k 2m,k k
(2.8)
W (x,y)=xcos(2πmx)V (y), m,k =1,2,....
2m−1,k k
Here, the sequences (2.6), (2.8) form a biorthonormal system of functions on D ; i.e., for any ad-
xy
missible indices m,k,l, and p, we have (Z ,W ) = 1 if m = l and k = p and (Z ,W ) = 0
m,k lp m,k l,p
otherwise. And each of these sequence is a basis in the space L (D ) [27]. Here, the inner product
2 xy
is given by
1 1
(cid:90) (cid:90)
(ψ,ξ)= ψ(x,y)ξ(x,y)dxdy.
0 0
2.2. Finite difference discretization
We represent the derivatives in a differential equation with finite difference discretization both on a
uniform and on a non-uniform grid in order to solve the differential equation numerically.
We construct an explicit finite difference scheme for this problem. First, let us show this scheme
on a uniform grid.
3
2.2.1. Uniform grid
Let us consider a uniform grid with increments hx and hy with respect to x and y respectively and
time increment ht. Set
x =i·hx, i=0,1,...,Nx, hx·Nx=1,
i
y =j·hy, j =0,1,...,Ny, hy·Ny =1,
j
t =n·ht, n=0,1,...,Nt, ht·Nt=T.
n
Let un =u(x ,y ,t ), then from [29] we take the derivatives as follows:
ij i j n
∂u(cid:12)(cid:12)(cid:12) = unij+1−unij,
∂t(cid:12) ht
(xi,yj,tn)
∂2u(cid:12)(cid:12)(cid:12) = uni+1j −2unij +uni−1j, ∂2u(cid:12)(cid:12)(cid:12) = unij+1−2unij +unij−1.
∂x2(cid:12) hx2 ∂y2(cid:12) hy2
(xi,yj,tn) (xi,yj,tn)
Now let’s see the explicit finite difference scheme on a non-uniform grid.
2.2.2. Non-uniform grid
When we consider a non-uniform grid, even if we don’t have fixed increments for the spatial axes like
on a uniform grid but we have fixed increment ht for time. Set
x , i=0,1,...,Nx,
i
y , y =0,1,...,Ny,
j
t =n·ht, n=0,1,...,Nt, ht·Nt=T.
n
We take the first derivatives as
∂u(cid:12)(cid:12)(cid:12) = unij −uni−1j, for i=1,...,Nx, ∂u(cid:12)(cid:12)(cid:12) = unij+1−unij,
∂x(cid:12) x −x ∂t(cid:12) ht
(xi,yj,tn) i i−1 (xi,yj,tn)
and we take the second order derivatives on the spatial axes as follows [28]
∂2u(cid:12)(cid:12)(cid:12) = 2uni−1j + 2unij + 2uni+1j ,
∂x2(cid:12) (x −x )(x −x ) (x −x )(x −x ) (x −x )(x −x )
(xi,yj,tn) i−1 i i−1 i+1 i i−1 i i+1 i+1 i−1 i+1 i
∂2u(cid:12)(cid:12)(cid:12) = 2unij−1 + 2unij + 2unij+1 .
∂y2(cid:12) (y −y )(y −y ) (y −y )(y −y ) (y −y )(y −y )
(xi,yj,tn) j−1 j j−1 j+1 j j−1 j j+1 j+1 j−1 j+1 j
2.3. Numerical Integration
We consider several numerical integration formulas both for uniform and non-uniform grids.
2.3.1. On uniform grid
Let
(cid:90)(cid:90)
I = u(x,y)dxdy,
Dxy
x =ih, y =jk, h=1/n, k =1/m.
i j
The product trapezoidal rule is
n m
(cid:88)(cid:88)1
I ≈hk (u(x ,y )+u(x ,y )+u(x ,y )+u(x ,y )).
4 i−1 j−1 i−1 j i j−1 i j
i=1j=1
4
The product Simpson’s rule is
n/2m/2
hk (cid:88)(cid:88)
I ≈ [u(x ,y )+4u(x ,y )+u(x ,y )]
9 2i−2 2j−2 2i−1 2j−2 2i 2j−2
i=1 j=1
+4[u(x ,y )+4u(x ,y )+u(x ,y )]
2i−2 2j−1 2i−1 2j−1 2i 2j−1
+[u(x ,y )+4u(x ,y )+u(x ,y )],
2i−2 2j 2i−1 2j 2i 2j
where m and n must be even.
2.3.2. On non-uniform grid
We have product rules obtained from Gauss-Lobatto nodes and weights on the unit square. We get
desired nodes by moving these nodes to D from the unit square. Hence, if the increments of the
xy
nodes are not equal, then we can use the following product rule for numerical integration.
n m (cid:18) (cid:19)
1(cid:88)(cid:88) 1 1 1 1
I ≈ A u x + , y + ,
4 ij 2 i 2 2 j 2
i=1j=1
where x and y are nodes of Gauss-Lobatto nodes on unit square and A and A are corresponding
i j i j
weights to these nodes, respectively with A =A A [32].
ij i j
3. Well-posedness of inverse problem
We have the following assumptions on ϕ(x,y), f(x,y,t) and E(t).
Theorem 3.1 (Existence and uniqueness) Under the conditions
(A ) ϕ(x,y)∈C2,2(D ),
1 1 xy
(A ) ϕ(0,y)=ϕ(1,y), ϕ(x,0)=ϕ(x,1)=0, ϕ (1,y)=0,
1 2 x
(A ) ϕ ≤0, ϕ ≤0, m,k =1,2...,
1 3 0,2k−1 2m,2k−1
(A ) f(x,y,t)∈C(D), f(x,y,t)∈C2,2(D ), ∀t∈[0,T]
2 1 xy
(A ) f(0,y)=f(1,y), f(x,0)=f(x,1)=0,
2 2
(cid:104) (cid:105)
(A ) f (t)≥0, min f (t)≥ max f (t) 1−e−(π(2k−1))2T−(2πm)2T ,
2 3 2m,2k−1 2m,2k−1 2m,2k−1
0≤t≤T 0≤t≤T
m=0,1,..., k =1,2,...,
(A ) E(t)∈C1[0,T],
3 1
1 1
(cid:90) (cid:90)
(A ) E(0)= ϕ(x,y)dxdy,
3 2
0 0
(A ) E(t)>0, E(cid:48)(t)≤0, ∀t∈[0,T],
3 3
the inverse problem (1.1)-(1.4) has a unique solution for small T.
Proof Since (2.6) is a basis in L (D ), we present the solution of (1.1)-(1.3) in the following form
2 xy
for arbitrary p(t)∈C[0,T]:
∞ (cid:34) ∞ ∞ (cid:35)
(cid:88) (cid:88) (cid:88)
u(x,y,t)= α (t)Z (x,y)+ α (t)Z (x,y)+ α (t)Z (x,y) , (3.1)
0,k 0,k 2m−1,k 2m−1,k 2m,k 2m,k
k=1 m=1 m=1
5
where
(cid:90) t
α0,k(t)=ϕ0,k·e−(πk)2t−(cid:82)0tp(s)ds+ f0,k(τ)e−(πk)2(t−τ)−(cid:82)τtp(s)dsdτ,
0
(cid:90) t
α2m,k(t)=ϕ2m,ke−(2πm)2t−(πk)2t−(cid:82)0tp(s)ds+ f2m,k(τ)e−(2πm)2(t−τ)−(πk)2(t−τ)−(cid:82)τtp(s)dsdτ,
0
α2m−1,k(t)=[ϕ2m−1,k−4πm·ϕ2m,kt]e−(2πm)2t−(πk)2t−(cid:82)0tp(s)ds
(cid:90) t
+ [f2m−1,k(τ)−4πmf2m,k(τ)(t−τ)]e−(2πm)2(t−τ)−(πk)2(t−τ)−(cid:82)τtp(s)dsdτ.
0
Here
(cid:90)(cid:90)
ϕ = ϕ(x,y)W (x,y)dxdy,
m,k m,k
Dxy
(cid:90)(cid:90)
f (t)= f(x,y,t)W (x,y)dxdy, m=0,1,2,..., k =1,2,....
m,k m,k
Dxy
Under conditions A and A , the series (3.1), its t−partial derivative, the xx−second order and
1 2
yy−second order partial derivatives converge uniformly in D that their majorizing sums absolutely
convergence. Thus u(x,y,t) ∈ C2,2,1(D). Since the conditions that the series (3.1) can be termwise
differentiable by t and E(t)∈C1[0,T] with (A ) the overdetermination condition (1.4) is equivalent
3 2
to
1 1
(cid:90) (cid:90)
u (x,y,t)dxdy =E(cid:48)(t). (3.2)
t
0 0
Therefore we have
P(p(t))=p(t) (3.3)
from the equations (3.1)-(3.2) such that
1 (cid:34) (cid:88)∞ 4√2 √
P(p(t))= −E(cid:48)(t)+ f (t)−4 2π(2k−1)
E(t) π(2k−1) 0,2k−1
k=1
(cid:20) (cid:90) t (cid:21)
× ϕ0,2k−1e−(π(2k−1))2t−(cid:82)0tp(s)ds+ f0,2k−1(τ)e−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ
0
∞ ∞ √ (cid:32) √ √ (cid:33)
(cid:88) (cid:88) 4 2 16 2m 4 2(2k−1)
+ f (t)− +
π2(2k−1)m 2m,2k−1 2k−1 m
k=1m=1
×(cid:104)ϕ2m,2k−1e−(2πm)2t−(π(2k−1))2t−(cid:82)0tp(s)ds
(cid:90) t (cid:21)
+ f2m,2k−1(τ)e−(2πm)2(t−τ)−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ .
0
Now let us show that P is a contraction mapping in C+[0,T], for small T, where
C+[0,T]={p(t)∈C[0,T]:p(t)≥0}.
Moreover, it is easy to see that
P :C+[0,T]→C+[0,T]
under conditions (A ) , (A ) and (A ) . For all p (t),p (t)∈C+[0,T],
1 3 2 3 3 3 1 2
6
|P(p (t))−P(p (t))|
1 2
≤|E1(t)|(cid:88)∞ 4√2π(2k−1)(cid:104)|ϕ0,2k−1|(cid:12)(cid:12)(cid:12)e−(cid:82)0tp1(s)ds−e−(cid:82)0tp2(s)ds(cid:12)(cid:12)(cid:12)
k=1
+(cid:90) t|f0,2k−1(τ)|(cid:12)(cid:12)(cid:12)e−(cid:82)τtp1(s)ds−e−(cid:82)τtp2(s)ds(cid:12)(cid:12)(cid:12)dτ(cid:21)
0
∞ ∞ (cid:32) √ √ (cid:33)
1 (cid:88) (cid:88) 16 2m 4 2(2k−1)
+ +
|E(t)| 2k−1 m
k=1m=1
×(cid:20)|ϕ2m,2k−1|(cid:12)(cid:12)(cid:12)e−(cid:82)0tp1(s)ds−e−(cid:82)0tp2(s)ds(cid:12)(cid:12)(cid:12)+(cid:90) t|f2m,2k−1(τ)|(cid:12)(cid:12)(cid:12)e−(cid:82)τtp1(s)ds−e−(cid:82)τtp2(s)ds(cid:12)(cid:12)(cid:12)dτ(cid:21).
0
By using the mean value theorem we get
(cid:12)(cid:12)(cid:12)e−(cid:82)0tp1(s)ds−e−(cid:82)0tp2(s)ds(cid:12)(cid:12)(cid:12)≤T · max |p1(t)−p2(t)|=T(cid:107)p1−p2(cid:107)C[0,T].
0≤t≤T
The last inequality yields to the following existence and uniqueness results:
(cid:107)P(p )−P(p )(cid:107) ≤β(cid:107)p −p (cid:107) ,
1 2 C[0,T] 1 2 C[0,T]
where
T (cid:32)(cid:88)∞ √ (cid:90) T (cid:88)∞ √
β = 4 2π(2k−1)|ϕ |+ 4 2π(2k−1)|f (τ)|dτ
min |E(t)| 0,2k−1 0,2k−1
0≤t≤T k=1 0 k=1
∞ ∞ (cid:32) √ √ (cid:33)
(cid:88) (cid:88) 16 2m 4 2(2k−1)
+ + |ϕ |
2k−1 m 2m,2k−1
k=1m=1
(cid:90) T (cid:88)∞ (cid:88)∞ (cid:32)16√2m 4√2(2k−1)(cid:33) (cid:33)
+ + |f (τ)|dτ .
2k−1 m 2m,2k−1
0 k=1m=1
In the case β <1 the map P is contraction map in C+[0,T]. Obviously, this inequality is satisfied for
small T. Hence, P has unique fixed point by Banach fixed point theorem. (cid:3)
The following result on continuously dependence on the data of the solution of inverse problem
(1.1)-(1.4) holds.
Theorem 3.2 (Stability) Under assumptions (A )−(A ), the solution (p,u) depends continuously
1 3
upon the data.
Proof Let Φ = {ϕ,E,f} and Φ = {ϕ,E,f} be two sets of data, which satisfy the conditions (A )−
1
(A ). Denote (cid:107)Φ(cid:107) = (cid:107)ϕ(cid:107) + (cid:107)E(cid:107) + (cid:107)f(cid:107) . Suppose that there exist positive
3 C2,2(Dxy) C1[0,T] C2,2,0(D)
constants M and M such that
1 2
0<M ≤|E|, 0<M ≤|E|, (cid:107)Φ(cid:107)≤M , (cid:107)Φ(cid:107)≤M .
1 1 2 2
Let (p,u) and (p,u) be the solutions of inverse problem (1.1)-(1.4) corresponding the data Φ and Φ,
7
respectively. According to (3.3)
1 (cid:34) (cid:88)∞ 4√2 √
p(t)= −E(cid:48)(t)+ f (t)−4 2π(2k−1)
E(t) π(2k−1) 0,2k−1
k=1
(cid:20) (cid:90) t (cid:21)
× ϕ0,2k−1e−(π(2k−1))2t−(cid:82)0tp(s)ds+ f0,2k−1(τ)e−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ
0
∞ ∞ √ (cid:32) √ √ (cid:33)
(cid:88) (cid:88) 4 2 16 2m 4 2(2k−1)
+ f (t)− +
π2(2k−1)m 2m,2k−1 2k−1 m
k=1m=1
×(cid:104)ϕ2m,2k−1e−(2πm)2t−(π(2k−1))2t−(cid:82)0tp(s)ds
(cid:90) t (cid:21)
+ f2m,2k−1(τ)e−(2πm)2(t−τ)−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ ,
0
1 (cid:34) (cid:48) (cid:88)∞ 4√2 √
p(t)= −E (t)+ f (t)−4 2π(2k−1)
E(t) π(2k−1) 0,2k−1
k=1
(cid:20) (cid:90) t (cid:21)
× ϕ0,2k−1e−(π(2k−1))2t−(cid:82)0tp(s)ds+ f0,2k−1(τ)e−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ
0
∞ ∞ √ (cid:32) √ √ (cid:33)
(cid:88) (cid:88) 4 2 16 2m 4 2(2k−1)
+ f (t)− +
π2(2k−1)m 2m,2k−1 2k−1 m
k=1m=1
×(cid:104)ϕ2m,2k−1e−(2πm)2t−(π(2k−1))2t−(cid:82)0tp(s)ds
(cid:90) t (cid:21)
+ f2m,2k−1(τ)e−(2πm)2(t−τ)−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ .
0
Let us estimate the difference p−p at first. It is easy to compute the followings
(cid:12)(cid:12)E(cid:48)(t) E(cid:48)(t)(cid:12)(cid:12)
(cid:12) − (cid:12)≤M (cid:107)E−E(cid:107) ,
(cid:12)(cid:12)E(t) E(t)(cid:12)(cid:12) 3 C1[0,T]
(cid:12)(cid:12)(cid:12)(cid:88)∞ 1 (cid:32)f0,2k−1(t) − f0,2k−1(t)(cid:33)(cid:12)(cid:12)(cid:12)≤M (cid:107)f −f(cid:107) +M (cid:107)E−E(cid:107) ,
(cid:12)(cid:12) 2k−1 E(t) E(t) (cid:12)(cid:12) 4 C2,2,0(D) 5 C1[0,T]
k=1
(cid:12) (cid:12)
(cid:12)(cid:12)(cid:12)(cid:12)(cid:88)∞ (2k−1)(cid:18)E1(t)ϕ0,2k−1e−(π(2k−1))2t−(cid:82)0tp(s)ds− E1(t)ϕ0,2k−1e−(π(2k−1))2t(cid:82)0tp(s)ds(cid:19)(cid:12)(cid:12)(cid:12)(cid:12)
k=1
≤M (cid:107)ϕ−ϕ(cid:107) +M T(cid:107)p−p(cid:107) +M (cid:107)E−E(cid:107) ,
6 C2,2(Dxy) 7 C[0,T] 8 C1[0,T]
(cid:12)
(cid:12)(cid:12)(cid:12)(cid:88)∞ (2k−1)(cid:18)E1(t)(cid:90) tf0,2k−1(τ)e−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ
(cid:12)k=1 0
−E1(t)(cid:90) tf0,2k−1(τ)e−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ(cid:19)(cid:12)(cid:12)(cid:12)(cid:12)
0
≤M T(cid:107)f −f(cid:107) +M T2(cid:107)p−p(cid:107) +M T(cid:107)E−E(cid:107) ,
9 C2,2,0(D) 10 C[0,T] 11 C1[0,T]
(cid:12)(cid:12)(cid:12)(cid:88)∞ (cid:88)∞ 1 (cid:32)f2m,2k−1(t) − f2m,2k−1(t)(cid:33)(cid:12)(cid:12)(cid:12)≤M (cid:107)f −f(cid:107) +M (cid:107)E−E(cid:107) ,
(cid:12)(cid:12) (2k−1)m E(t) E(t) (cid:12)(cid:12) 12 C2,2,0(D) 13 C1[0,T]
k=1m=1
(cid:12)
(cid:12)(cid:12)(cid:12)(cid:88)∞ (cid:88)∞ (cid:18)2k4m−1 + 2km−1(cid:19)(cid:18)ϕ2Em,(2tk)−1e−(2πm)2t−(π(2k−1))2t−(cid:82)0tp(s)ds
(cid:12)
k=1m=1
ϕ (cid:19)(cid:12)
− 2m,2k−1e−(2πm)2t−(π(2k−1))2t−(cid:82)0tp(s)ds (cid:12)(cid:12)
E(t) (cid:12)
≤M (cid:107)ϕ−ϕ(cid:107) +M T(cid:107)p−p(cid:107) +M (cid:107)E−E(cid:107) ,
14 C2,2(Dxy) 15 C[0,T] 16 C1[0,T]
8
(cid:12)
(cid:12)(cid:12)(cid:12)(cid:88)∞ (cid:88)∞ (cid:18)2k4m−1 + 2km−1(cid:19)(cid:18)E1(t)(cid:90) tf2m,2k−1(τ)e−(2πm)2(t−τ)−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ
(cid:12)k=1m=1 0
+E1(t)(cid:90) tf2m,2k−1(τ)e−(2πm)2(t−τ)−(π(2k−1))2(t−τ)−(cid:82)τtp(s)dsdτ(cid:19)(cid:12)(cid:12)(cid:12)(cid:12)
0
≤M T(cid:107)f −f(cid:107) +M T2(cid:107)p−p(cid:107) +M T(cid:107)E−E(cid:107) ,
17 C2,2,0(D) 18 C[0,T] 19 C1[0,T]
where M , i=3,4,...,19 are constants that are determined from M and M . Then we get
i 1 2
(cid:16) (cid:17)
(1−M )(cid:107)p−p(cid:107) ≤M (cid:107)ϕ−ϕ(cid:107) +(cid:107)E−E(cid:107) +(cid:107)f −f(cid:107) ,
20 C[0,T] 21 C2,2(Dxy) C1[0,T] C2,2,0(D)
√ √ √ √
where M =4 2πT(M +TM )+4 2T(M +TM ) and M =max{4 2/πM +4 2πTM +
√ 20 √ 7 10√ √15 18 21 √ √4 9
4 2/π2M +4 2TM ,M +4 2/πM +4 2π(M +TM )+4 2/π2M +4 2(M +TM ),
√ 12 √ 17 3 5 8 11 13 16 19
4 2πM +4 2M }. The inequality M <1 holds for small T. Finally, we have
6 14 20
M
(cid:107)p−p(cid:107) ≤M (cid:107)Φ−Φ(cid:107), M = 21 .
C[0,T] 22 22 1−M
20
Similarly, we obtain the estimate the difference u−u from (3.1):
(cid:107)u−u(cid:107) ≤M (cid:107)Φ−Φ(cid:107).
C(D) 23
(cid:3)
Now, let’s see how we implement the numerical solution of the inverse problem in two different
methods.
4. Numerical methods for inverse problem
We will consider the examples of numerical solution of the inverse problem (1.1)-(1.4). For the
convenience of discussion of the numerical method, we will rewrite the equations as follows:
∂v/∂t=∂2v/∂x2+∂2v/∂y2+r(t)f(x,y,t), (4.1)
v(x,y,0)=ϕ(x,y), 0≤x,y ≤1, (4.2)
v(0,y,t)=v(1,y,t), v (1,y,t)=0, 0≤y ≤1, 0≤t≤T,
x
(4.3)
v(x,0,t)=v(x,1,t)=0, 0≤x≤1, 0≤t≤T,
1 1
(cid:90) (cid:90)
v(x,y,t)dxdy =E(t)r(t), (4.4)
0 0
by using the transformations
r(t)=e(cid:82)0tp(s)ds,
v(x,y,t)=u(x,y,t)r(t).
4.1. Uniform finite difference method
Weconsiderthefollowingexplicitfinitedifferencediscretizationoftheproblem(4.1)-(4.4)onauniform
grid.
vn+1−vn vn −2vn +vn vn −2vn +vn
ij ij = i+1j ij i−1j + ij+1 ij ij−1 +rnfn,
ht hx2 hy2 ij
(cid:18) (cid:19)
ht ht ht ht ht ht
⇒vn+1 = vn + vn + 1−2 −2 vn + vn + vn +ht·rnfn
ij hx2 i+1j hy2 ij+1 hx2 hy2 ij hx2 i−1j hy2 ij−1 ij
with the initial condition
ϕ =v0,
ij ij
9
nonlocal boundary conditions
vn −vn
vn =vn , Nxj Nx−1j =0⇒vn =vn ,
0j Nxj hx Nxj Nx−1j
vn =vn =0,
i0 iNy
and overdetermination condition
1 1
1 (cid:90) (cid:90)
rn = v(x,y,t )dxdy, (4.5)
En n
0 0
where
vn =v(x ,y ,t ), rn =r(t ), En =E(t ), fn =f(x ,y ,t ), n=0,...,Nt,
ij i j n n n ij i j n
i=0,...,Nx, j =0,...,Ny.
1 1
(cid:82) (cid:82)
When we approximate v(x,y,t)dxdy by the product trapezoidal rule
0 0
(cid:90)1(cid:90)1 (cid:88)Nx (cid:88)Ny 1
v(x,y,t)dxdy ≈hx·hy (v(x ,y ,t)+v(x ,y ,t)+v(x ,y ,t)+v(x ,y ,t)),
4 i−1 j−1 i−1 j i j−1 i j
0 0 i=1j=1
where t∈[0,T]. And similarly by the product Simpson’s rule
(cid:90)1(cid:90)1 hx·hy N(cid:88)x/2N(cid:88)y/2
v(x,y,t)dxdy ≈ [v(x ,y ,t)+4v(x ,y ,t)+v(x ,y ,t)]
9 2i−2 2j−2 2i−1 2j−2 2i 2j−2
0 0 i=1 j=1
+4[v(x ,y ,t)+4v(x ,y ,t)+v(x ,y ,t)]
2i−2 2j−1 2i−1 2j−1 2i 2j−1
+[v(x ,y ,t)+4v(x ,y ,t)+v(x ,y ,t)],
2i−2 2j 2i−1 2j 2i 2j
where Nx and Ny must be even.
4.2. Non-uniform finite difference method
Let us consider the same problem on a non-uniform grid.
vn+1−vn (cid:18) vn vn vn (cid:19)
ij ij =2 i−1j + ij + i+1j
ht (x −x )(x −x ) (x −x )(x −x ) (x −x )(x −x )
i−1 i i−1 i+1 i i−1 i i+1 i+1 i−1 i+1 i
(cid:18) vn vn vn (cid:19)
+2 ij−1 + ij + ij+1
(y −y )(y −y ) (y −y )(y −y ) (y −y )(y −y )
j−1 j j−1 j+1 j j−1 j j+1 j+1 j−1 j+1 j
+rnfn,
ij
2ht 2ht
⇒vn+1 = vn + vn
ij (x −x )(x −x ) i−1j (y −y )(y −y ) ij−1
i−1 i i−1 i+1 j−1 j j−1 j+1
(cid:18) (cid:19)
2ht 2ht
+ 1+ + vn
(x −x )(x −x ) (y −y )(y −y ) ij
i i−1 i i+1 j j−1 j j+1
2ht 2ht
+ vn + vn +ht·rnfn,
(x −x )(x −x ) i+1j (y −y )(y −y ) ij+1 ij
i+1 i−1 i+1 i j+1 j−1 j+1 j
with the initial condition
ϕ =v0,
ij ij
10