Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics
Подождите немного. Документ загружается.
216 Chapter 6 Right-hand side identification
The problem is being solved on a uniform grid with h = 0.02 and τ = 0.02.
Correctness of the problem is illustrated by calculation data obtained for the prob-
lem with perturbed boundary conditions at t = T (function ϕ(x) in (6.113)). In the
framework of the quasi-real experiment, the obtained solution was perturbed at t = T
by the law
ϕ
δ
(x) = ϕ(x) + 2δ(σ(x) − 1/2), x ∈ ω,
where σ(x) is a random function normally distributed over the interval [0, 1], and the
parameter δ defines the inaccuracy level.
Figure 6.13 shows the exact and reconstructed dependence of the right-hand side on
the spatial variable, and also the exact and approximate difference end-time solutions
of the direct problem obtained with δ = 0.0005. Similar data obtained with an inac-
curacy increased to δ = 0.001 and δ = 0.0002 are shown respectively in Figures 6.14
and 6.15. Considerable sensitiveness of the reconstructed right-hand side to the so-
lution inaccuracy at t = T is worth noting: at a relative inaccuracy of 0.1 % in the
given ϕ(x), the right-hand side could be determined accurate to approximately 25 %.
Figure 6.13 Solution of the problem obtained with δ = 0.0005
Section 6.4 Identification of time-independent right-hand side: parabolic equations 217
Figure 6.14 Solution of the problem obtained with δ = 0.0005
Figure 6.15 Solution of the problem perturbed with δ = 0.0002
218 Chapter 6 Right-hand side identification
6.5 Reconstruction of the right-hand side of an ellipti-
cal equation from observation data obtained at the
boundary
Below, we consider the classical inverse problem in the potential theory in which it is
required to determine the unknown right-hand side of the elliptic equation in the case
in which additional data are given on the boundary of the calculation domain.
6.5.1 Statement of the inverse problem
Consider a model inverse problem in which it is required to determine the unknown
right-hand side from observation data obtained at the domain boundary. To simplify
the consideration, restrict ourselves to the two-dimensional Poisson equation. Con-
sider first the formulation of the direct problem.
In a bounded domain the function u(x), x = (x
1
, x
2
) satisfies the equation
−u ≡−
2
α=1
∂
2
u
∂x
2
α
= f (x), x ∈ . (6.144)
Consider a Dirichlet problem in which equation (6.144) is supplemented with the fol-
lowing first-kind homogeneous boundary conditions:
u(x) = 0, x ∈ ∂. (6.145)
The direct problem is formulated in the form (6.144), (6.145), with known right-hand
side f (x) in (6.144).
Among the inverse problems for elliptic equations, consider the right-hand side
identification problem. We assume that additional measurements are feasible only
on the domain boundary. In addition to (6.145), the following second-kind boundary
conditions are also considered:
∂u
∂n
(x) = μ(x), x ∈ ∂, (6.146)
where n is the external normal to .
In this general formulation the solution of the inverse problem in which it is required
to determine the pair of functions {u(x), f (x)} from conditions (6.144)–(6.146) is not
unique. The latter statement requires no special comments: it suffices to consider the
inverse problem in a circle with the right-hand side dependent on the distance from
the center of the circle. The non-uniqueness stems from the fact that we are trying to
reconstruct a two-dimensional function (the right-hand side f (x) from a function with
lower dimensionality (μ(x), x ∈ ∂).
Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 219
6.5.2 Uniqueness of the inverse-problem solution
Unique determination of the right-hand side is possible in the case in which the un-
known right-hand side is independent of one of the variables. In fact, the latter was
the case of the right-hand side identification problem for a parabolic equation, with the
unknown dependence of the right-hand side on time or spatial variables to be recon-
structed.
Not trying to consider the general case, turn to a typical example. We assume that
the right-hand side (6.144) can be represented as
f (x) = ϕ
1
(x
2
) + x
1
ϕ
2
(x
2
). (6.147)
We pose a problem in which it is required to determine two functions ϕ
α
(x
2
), α = 1, 2,
independent of one of the variables (namely, of the variable x
1
), from (6.144)–(6.146).
We reformulate the inverse problem (6.144)–(6.147) by eliminating the unknown
functions ϕ
α
(x
2
), α = 1, 2. Double differentiation of (6.144) with respect to x
1
with
allowance for (6.147) gives:
∂
2
∂x
2
1
u = 0, x ∈ . (6.148)
In this way, we arrive at a boundary value problem for the composite equation (6.145),
(6.146), (6.148).
Let us show that the solution of problem (6.145), (6.146), (6.148) is unique. For this
to be shown, it suffices to prove that the solution of the problem with homogeneous
boundary conditions
∂u
∂n
(x) = 0, x ∈ ∂, (6.149)
is u(x) ≡ 0, x ∈ .
We multiply equation (6.148) by u(x) and perform integration over the whole do-
main ); this yields
∂
2
∂x
2
1
uudx = 0.
Taking into account the homogeneous boundary conditions (6.145) and permutability
of the operators ∂/∂x
1
and , we obtain
vvdx = 0,v=
∂u
∂x
1
.
The pair of homogeneous boundary conditions (6.145), (6.149) guarantees that
v(x) = 0, x ∈ ∂.
Under these conditions, we have
vvdx =
2
α=1
∂v
∂x
α
2
dx = 0
220 Chapter 6 Right-hand side identification
Figure 6.16 Calculation domain
and, hence, v(x) = 0 throughout the whole domain . From
∂u
∂x
1
= 0, x ∈
and boundary conditions (6.145) it follows that the only solution of problem (6.145),
(6.148), (6.149) is u(x) ≡ 0, x ∈ .
More informative a priori estimates for the solution of the boundary value problem
(6.145), (6.146), (6.148) can also be obtained. This matter will be considered on the
difference level below.
6.5.3 Difference problem
We assume the calculation domain to be a rectangle:
={x | x = (x
1
, x
2
), 0 < x
α
< l
α
,α= 1, 2}.
For the sides of we use conditions indicated in Figure 6.16, so that
∂ =
1
∪
2
∪
3
∪
4
.
We seek the right-hand side of (6.144) in class (6.147) under the following addi-
tional conditions posed on the sides
2
and
4
of the rectangle:
∂u
∂x
1
(0, x
2
) = μ
1
(x
2
),
∂u
∂x
1
(l
1
, x
2
) = μ
2
(x
2
). (6.150)
With boundary conditions given on
1
and/or
3
(see (6.146)), the problem becomes
overspecified.
Along both direction x
α
, α = 1, 2, we introduce a uniform grid
ω
α
={x
α
| x
α
= i
α
h
α
, i
α
= 0, 1,...,N
α
, N
α
h
α
= l
α
},
Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 221
so that
ω
α
={x
α
| x
α
= i
α
h
α
, i
α
= 1, 2,...,N
α
− 1, N
α
h
α
= l
α
},
∂ω
α
={x
α
| x
α
= 0, l
α
}.
For the grid in the rectangle we use the notations
ω = ω
1
× ω
2
={x | x = (x
1
, x
2
), x
α
∈ ω
α
,α= 1, 2},
ω = ω
1
× ω
2
.
In the standard notation adopted in the theory of difference schemes, at internal
nodes we define the difference Laplace operator
y = y
¯x
1
x
1
+ y
¯x
2
x
2
, x ∈ ω.
We put into correspondence to the direct problem (6.144), (6.145) the difference prob-
lem
−y = f (x), x ∈ ω, (6.151)
y(x) = 0, x ∈ ∂ω. (6.152)
In the inverse problem, the right-hand side is sought in the class (6.147) from the
additional conditions (6.150). To pass to a difference analogue of (6.145), (6.147),
(6.150), we define the mesh function v =−y not only at internal nodes (see (6.151)),
but also on the set of boundary nodes.
We can conveniently introduce fictitious nodes with i
1
=−1 and i
1
= N
1
+ 1to
extend the grid over the variable x
1
by one node from either side. We approximate
boundary conditions (6.150) on the extended grid. Accurate to O(h
2
1
),wehave:
y(h
1
, x
2
) − y(−h
1
, x
2
)
2h
1
= μ
1
(x
2
), (6.153)
y(l
1
+ h
1
, x
2
) − y(l
1
− h
1
, x
2
)
2h
1
= μ
2
(x
2
). (6.154)
Taking into account the boundary conditions (6.145) at the left boundary, we obtain
v(0, x
2
) =−y(0.x
2
) =−
y(h
1
, x
2
) − 2y(0, x
2
) + y(−h
1
, x
2
)
h
2
1
.
With (6.153) taken into account, we arrive at the expression
v(0, x
2
) =−
2
h
2
1
y(h
1
, x
2
) +
2
h
1
μ
1
(x
2
). (6.155)
In a similar way, on
4
we obtain
v(l
1
, x
2
) =−
2
h
2
1
y(l
1
− h
1
, x
2
) −
2
h
1
μ
1
(x
2
). (6.156)
222 Chapter 6 Right-hand side identification
Double difference differentiation (6.151) yields the equation
v
¯x
1
x
1
= 0, x ∈ ω. (6.157)
The boundary conditions for this equation have the form (6.154), (6.155). With
known v, the solution is to be determined (see (6.151), (6.152)) from
−y = v(x), x ∈ ω, (6.158)
y(x) = 0, x ∈ ∂ω. (6.159)
In this manner, we arrive at a system of two difference Poisson equations for the pair
{y,v}. These two equations are interrelated via boundary conditions (6.155), (6.156).
We can conveniently reformulate the boundary value problem with non-
homogeneous boundary conditions (6.155)–(6.157) as a problem with homoge-
neous boundary conditions for a non-homogeneous equation at the internal nodes.
With (6.155), at near-boundary nodes we have:
2v(h
1
, x
2
) − v(2h
1
, x
2
)
h
2
1
=
v(0, x
2
)
h
2
1
=−
2
h
4
1
y(h
1
, x
2
) +
2
h
3
1
μ
1
(x
2
).
From (6.156), we obtain
2v(l
1
− h
1
, x
2
) − v(l
1
− 2h
1
, x
2
)
h
2
1
=−
2
h
4
1
y(l
1
− h
1
, x
2
) −
2
h
3
1
μ
2
(x
2
).
Let us define difference operators A
α
, α = 1, 2 on the set of mesh functions vanish-
ing at the boundary nodes:
A
α
y =−y
¯x
α
x
α
, x ∈ ω.
Then, the boundary value problem (6.155)–(6.157) can be written as
A
1
v =−A
0
y + φ, x ∈ ω. (6.160)
Here, the difference operator A
0
is defined by
A
0
y =
⎧
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎩
2
h
4
1
y(h
1
, x
2
), x
1
= h
1
,
0, h
1
< x
1
< l
1
− h
1
,
2
h
4
1
y(l
1
− h
1
, x
2
), x
1
= l
1
− h
1
.
The right-hand side φ in nonzero only at near-boundary nodes:
φ =
⎧
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎩
−
2
h
3
1
μ
1
(x
2
), x
1
= h
1
,
0, h
1
< x
1
< l
1
− h
1
,
2
h
3
1
μ
2
(x
2
), x
1
= l
1
− h
1
.
Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 223
Figure 6.17 The mesh pattern for the difference scheme
In the introduced notation, the boundary value problem (6.158), (6.159) takes the
form
(A
1
+ A
2
)y = v, x ∈ ω. (6.161)
In this way, we pass from (6.155)–(6.159) to the system (6.160), (6.161). The latter
system can be conveniently written as the following single operator equation:
Ay = φ. (6.162)
Here,
A = ( A
1
+ A
2
)A
1
+ A
0
. (6.163)
The mesh pattern used in this difference scheme is shown in Figure 6.17.
In the ordinary way, in the Hilbert space H = L
2
(ω) we introduce the scalar product
and the norm:
(y,w) ≡
x∈ω
y(x)w(x)h
1
h
2
, y≡(y, y)
1/2
.
In H ,
A
α
= A
∗
α
> 0, A
0
= A
∗
0
≥ 0
and, hence, in (6.162) we have A = A
∗
> 0. By virtue of this, the difference problem
(6.162) has a unique solution.
We gave a sufficiently general scheme for constructing a discrete analogue to the
non-classical boundary value problem (6.145), (6.148), (6.150) suitable for solving
even more complex problems. In the case under consideration, we can substantially
simplify the problem by explicitly writing the solution of problem (6.155)–(6.157).
It should be noted that such transformation allows for specific features of the inverse
problem in the greatest possible extent and most clearly on the difference level.
The general solution of the difference equation (6.157) is a linear function of x
1
:
v(x
1
, x
2
) =
1 −
x
1
l
1
v(0, x
2
) +
x
1
l
1
v(l
1
, x
2
).
224 Chapter 6 Right-hand side identification
With boundary conditions (6.155), (6.156) taken into account, we obtain
v(x
1
, x
2
) =−
1 −
x
1
l
1
2
h
2
1
y(h
1
, x
2
) −
x
1
l
1
2
h
2
1
y(l
1
−h
1
, x
2
) +ψ(x
1
, x
2
), (6.164)
where
ψ(x
1
, x
2
) =
1 −
x
1
l
1
2
h
1
μ
1
(x
2
) −
x
1
l
1
2
h
1
μ
2
(x
2
).
Substitution of (6.164) into (6.151) leads us to the difference equation
−y +
1 −
x
1
l
1
2
h
2
1
y(h
1
, x
2
) +
x
1
l
1
2
h
2
1
y(l
1
−h
1
, x
2
) = ψ(x), x ∈ ω. (6.165)
In this way, the solution of the inverse right-hand side identification problem for
equation (6.151) in class (6.147) has reduced to the solution of the boundary value
problem (6.152), (6.165). Equation (6.165) is a loaded difference equation.
6.5.4 Solution of the difference problem
To find the solution of the difference problem, we use the variable separation method.
This approach can be applied to the difference problem (6.152) with the composite
difference operator A defined by (6.153). Let us dwell on a second possibility, when
to be sought is the solution of the boundary-layer problem for the loaded difference
elliptic equation (6.152), (6.165).
We write problem (6.152), (6.165) as the equation
(A
1
+ A
2
)y +q
1
(x
1
)y(h
1
, x
2
) +q
2
(x
1
)y(l
1
−h
1
, x
2
) = ψ(x), x ∈ ω, (6.166)
in which q
α
≥ 0, α = 1, 2.
We denote as λ
k
, v
k
(x
2
), k = 1, 2,...,N
2
− 1 the eigenvalues and eigenfunctions
of A
2
:
A
2
v = λv.
The solution of this difference spectral problem is well known:
λ
k
=
4
h
2
2
sin
2
kπ h
2
2l
2
,v
k
(x
2
) =
2
l
2
sin
kπ x
2
l
2
.
The eigenfunctions are orthonormal functions:
(v
k
,v
m
)
(2)
= δ
km
,δ
km
=
1, k = m,
0, k = m,
where
(v, w)
(2)
=
N
2
−1
i
2
=1
v(x
2
)w(x
2
)h
2
Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 225
is the scalar product in L
2
(ω
2
).
We seek the solution of problem (6.166) as an expansion in the eigenfunctions of
A
2
:
y(x
1
, x
2
) =
N
2
−1
k=1
c
k
(x
1
)v
k
(x
2
).
Substitution into (6.166) leads us to the necessity to solve the difference problems
(A
1
+ λ
k
)c
k
(x
1
) + q
1
(x
1
)c
k
(h
1
) + q
2
(x
1
)c
k
(l
1
− h
1
) = ψ
k
(x
1
), (6.167)
where
ψ
k
(x
1
) = (ψ, v
k
)
(2)
, k = 1, 2,...,N
2
− 1.
The matter of solution of these N
2
−1 one-dimensional difference problems should
be given particular attention. In the case of (6.167) it is required to find the solution of
the difference boundary value problem
−w
¯x
1
x
1
+ λw + q
1
(x
1
)w(h
1
) + q
2
(x
1
)w(l
1
− h
1
) = r(x
1
), (6.168)
w(0) = 0,w(l
1
) = 0. (6.169)
We represent the solution of problem (6.168), (6.169) as
w(x
1
) = s(x
1
) + s
(1)
(x
1
)w(h
1
) + s
(2)
(x
1
)w(l
1
− h
1
). (6.170)
We insert this expression into (6.168) and isolate the functions s
α
, α = 1, 2, collecting
the terms with w(h
1
), w(l
1
−h
1
) and equating them to zero. This gives us three three-
point difference equations for the auxiliary functions s(x
1
), s
α
(x
1
), α = 1, 2. With
allowance for (6.169), we assume the boundary conditions to be homogeneous, so that
−s
¯x
1
x
1
+ λs = r (x
1
), (6.171)
s(0) = 0, s(l
1
) = 0, (6.172)
−s
(1)
¯x
1
x
1
+ λs
(1)
=−q
1
(x
1
), (6.173)
s
(1)
(0) = 0, s
(1)
(l
1
) = 0, (6.174)
−s
(2)
¯x
1
x
1
+ λs
(2)
=−q
2
(x
1
), (6.175)
s
(2)
(0) = 0, s
(2)
(l
1
) = 0. (6.176)
After solving the three standard problems (6.171)–(6.176), the functions w(h
1
) and
w(l
1
− h
1
) can be found. From representation (6.170), it readily follows that
w(h
1
) = s(h
1
) + s
(1)
(h
1
)w(h
1
) + s
(2)
(h
1
)w(l
1
− h
1
), (6.177)
w(l
1
− h
1
) = s(l
1
− h
1
) + s
(1)
(l
1
− h
1
)w(h
1
) + s
(2)
(l
1
− h
1
)w(l
1
− h
1
). (6.178)