Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics

Подождите немного. Документ загружается.

416 Chapter 8 Other problems

from (8.168)–(8.170) we obtain:

−v + c

)v + θ(x

(x) = 0, x ∈ , (8.171)

v(x) = 0, x ∈ ∂, (8.172)

∂v

∂n

(x) = 0, x ∈ ∂. (8.173)

The solution uniqueness for the inverse problem (8.165)–(8.167) will be proved if we

show that equalities (8.171)–(8.173) are only valid with v(x) = 0, θ(x

) = 0, x ∈ .

Problem (8.171)–(8.173) is considered with given c

) and u

(x). Here, we deal

with an inverse (but linear) problem in which it is required to determine the pair v(x),

θ(x

), the identiﬁcation problem for the right-hand side. Let us formulate some sufﬁ-

cient conditions that guarantee that v(x) = 0, θ(x

) = 0, x ∈ . Apart from the usual

assumptions that the solution and the coefﬁcient are smooth functions, we assume that

the solution u(x) in (8.165)–(8.167) is a constant-sign function, for instance,

u(x)>0, x ∈

.

This condition is guaranteed (maximum principle) by the assumption that

c(x

) ≥ 0, x ∈ , ψ(x)>0, x ∈ ∂.

Under such constraints, from (8.171) we obtain the following composite equation:

∂

∂x



(−v + c

)v)



= 0, x ∈ . (8.174)

It is required to prove that the solution of the boundary value problem (8.172)–(8.174)

is v(x) = 0. Then, from (8.171) it immediately follows (u

> 0) that θ(x

) = 0as

well.

We multiply equation (8.174) by some function η(x) such that

η(x) = 0, x ∈ ∂, (8.175)

and integrate the resulting equation over the domain . In view of (8.175), we obtain



w(−v + c

)v) dx = 0, (8.176)

where the following notation is used:

w(x) =

(x)

∂η(x)

∂x

With regard to the homogeneous boundary conditions (8.172) and (8.173), from

(8.176) we obtain:



(−w + c

)w)v dx = 0.

Section 8.5 Coefﬁcient inverse problem for elliptic equation 417

This equality yields v(x) = 0, x ∈  provided that we can choose η(x) so that

−w + c

)w = v(x), x ∈  (8.177)

holds. For the calculation domain under consideration (see Figure 8.26), consider the

following boundary value problem for equation (8.177) with mixed boundary condi-

tions

w(x) = 0, x ∈ 

∪ 

, (8.178)

∂w

∂n

(x) = 0, x ∈ 

∪ 

. (8.179)

There is no doubt that the solution of the standard boundary value problem (8.177)–

(8.179) does exist. In more general problems, we can use the result about solution

existence for equation (8.177), in the case in which the boundary conditions are given

in the form

∂w

∂x

(x) = 0, x ∈ ∂;

here we have a boundary value problem for the second-order elliptic equation with

given angled derivatives on the boundary.

It only remains to identify the form of η(x). To this end, using the function w(x)

found from (8.177)–(8.179), we solve the boundary value problem

∂

∂x



(x)

∂η(x)

∂x



∂w

∂x

(x), x ∈ 

with boundary conditions (8.175).

8.5.3 Difference inverse problem

Let us formulate the difference analogue of the coefﬁcient inverse problem (8.165)–

(8.167). In , we introduce a grid with step sizes h

, α = 1, 2 uniform along both

directions. First, we deﬁne the set of internal grid nodes

ω ={x | x = (x

, x

), x

= i

, i

= 1, 2,...,N

−1, N

= l

,β= 1, 2}

and let ∂ω be the set of boundary nodes. We designate as ∂ω

∗

the set of near-boundary

nodes, i.e.,

∂ω

∗

={x | x ∈ ω, x

= h

, x

= l

− h

,β= 1, 2}.

First, we approximate the direct problem (8.165), (8.166). For internal nodes, we

deﬁne the standard two-dimensional difference Laplace operator  on a ﬁve-point

mesh pattern:

y = y

¯x

+ y

¯x

, x ∈ ω.

418 Chapter 8 Other problems

We put in correspondence to the boundary value problem (8.165), (8.166) the differ-

ence Dirichlet problem

−y + c(x

)y = 0, x ∈ ω, (8.180)

y(x) = ϕ(x), x ∈ ∂ω. (8.181)

In the inverse problem, the mesh function c(x

), x

∈ ω

is unknown, but additional

conditions are given which correspond to (8.167). On the passage to the discrete prob-

lem, additional conditions can be put in correspondence with setting the mesh solution

at the nodes of ∂ω

∗

. The latter is the case in using the simplest approximation with

directed differences (with ﬁrst-order approximation accuracy). A similar situation is

observed in the case in which second-order approximations of the boundary condition

(8.167) are used. As an illustration, consider approximations of boundary conditions

on 

For the near-boundary node we have

u(h

, x

) − u(0, x

)

∂u

∂x

(0, x

) +

∂

∂x

(0, x

) + O(h

). (8.182)

With regard to the boundary conditions (8.166), (8.167), we arrive at the expression

y(h

, x

) = ϕ(h

, x

) − h

ψ(h

, x

The use of this expression corresponds to the approximation of boundary conditions

(8.167) of ﬁrst order. The higher (second) approximation order is achieved for the

solution of the boundary value problem (8.165), (8.166). In this case, it follows from

(8.165) that

∂

∂x

(0, x

) = u(0, x

) −

∂

∂x

(0, x

) = c(x

)ϕ(0, x

) −

∂

∂x

(0, x

With regard to (8.182), for the near-boundary nodes we can put

y(h

, x

) = ϕ(h

, x

) − h

ψ(h

, x

) + c(x

)ϕ(0, x

) − ϕ

¯x

(0, x

Yet, such approximations are not too much suitable for us in solving the inverse prob-

lem since, here, the value at the near-boundary node depends on the unknown (sought)

coefﬁcient c(x

). That is why we will use the simplest ﬁrst-kind approximations; in

this case, the additional conditions can be formulated in the form

y(x) = φ(x), x ∈ ∂ω

∗

. (8.183)

In the solution of applied problems the input data are normally given with some

inaccuracy. We assume that, compared to these inaccuracies, the approximation inac-

curacies can be disregarded with sufﬁciently ﬁne calculation grids used. In the exam-

ple under consideration (8.165)–(8.167), we restrict ourselves to the case in which the

Section 8.5 Coefﬁcient inverse problem for elliptic equation 419

largest inaccuracies are the measurement inaccuracies of the normal derivative at the

boundary. Hence, in approximating the solution of the inverse problem with input-data

inaccuracies, we can leave the boundary condition (8.181) unchanged, and instead of

(8.183), put

y(x) ≈ φ

(x), x ∈ ∂ω

∗

, (8.184)

where the parameter δ deﬁnes the inaccuracy level.

8.5.4 Iterative solution of the inverse problem

In the approximate solution of the inverse problem (8.180), (8.181), (8.184), we will

use gradient iteration methods. Here, the sought mesh function c(x

) is to be reﬁned

at each iteration step using the criterion for minimization of the discrepancy functional

(the discrepancy here is the accuracy with which the condition (8.184) is fulﬁlled).

We deﬁne the discrepancy as

J (c) =



x∈∂ω

∗

(y(x) − φ

(x))

h(x), (8.185)

where

h(x) =

⎧

⎨

⎩

, x

= h

, l

− h

, x

= h

, l

− h

, x

= h

, l

− h

, x

= h

, l

− h

+ h

)/2, x

= h

, l

− h

, x

= h

, l

− h

We assume that the sought function c(x

) belongs to the space of mesh functions

(ω

), in which for the scalar product and norm we use the settings

(c, d) =



∈ω

c(x

)d(x

, c=



(c, c).

To derive an expression for the gradient of J (c), we assume that, to the increment δc,

some increment δ J of the functional (8.185) and some increment δy of the solution of

problem (8.180), (8.181) corresponds.

Accurate to terms of second-order smallness, from (8.180) and (8.181) we obtain:

−δy + c(x

)δy + δcy= 0, x ∈ ω, (8.186)

δy(x) = 0, x ∈ ∂ω. (8.187)

For the discrepancy functional gradient, we immediately obtain:

δ J (c) = 2



x∈∂ω

∗

(y(x) − φ

(x)) δyh(x). (8.188)

The gradient of J



δ J (c) = ( J



(c), c).

420 Chapter 8 Other problems

To rearrange the right-hand side of (8.188), we multiply the equation (8.186) by

some mesh function ξ(x) h

, x ∈ ω and sum up the obtained equation over all

nodes in ω:



x∈ω

ξ(x)(−δy + c(x

)δy + δcy)h

= 0. (8.189)

We assume that

ξ(x) = 0, x ∈ ∂ω. (8.190)

Under such constraints, equation (8.189) yields:



x∈ω

δy (−ξ + c(x

)ξ)h



x∈ω

δcyξh

= 0. (8.191)

To derive the desired representation for J



(c), we relate the ﬁrst term in (8.191) with

the right-hand side of (8.188).

Let the function ξ(x) be deﬁned as the solution of the equation

−ξ + c(x

)ξ =−F (x), x ∈ ω. (8.192)

Here, the right-hand side F(x) is

F(x) =

⎧

⎪

⎨

⎪

⎩

2(y(x) − φ

(x))/ h

, x

= h

, l

− h

, x

= h

, l

− h

2(y(x) − φ

(x))/ h

, x

= h

, l

− h

, x

= h

, l

− h

+ h

(y(x) − φ

(x)), x

= h

, l

− h

, x

= h

, l

− h

0, x ∈ ω, x ∈ ∂ω

∗

For this right-hand side, from (8.188), (8.191), and (8.192) we obtain:

δ J (c) =



x∈ω

δcyξh



∈ω

δc





∈ω

y ξ h



With this equality, for the functional gradient we have the representation





∈ω

y ξ h

, x

∈ ω. (8.193)

To calculate the gradient, we have to solve the boundary value problem (8.180), (8.181)

for the ground state (the mesh function y(x)) and the boundary value problem (8.190),

(8.192) for the conjugate state (the mesh function ξ(x)).

In the case of the two-layer gradient iteration method the quantity c

) (k is the

iteration number) is to be reﬁned using the scheme

k+1

− c

k+1

+ J



) = 0, x

∈ ω

, k = 0, 1,... . (8.194)

Here, it should be taken into account that the equation J



nonlinear equation. The iterative process (8.194) can be terminated by the discrepancy

criterion.

Section 8.5 Coefﬁcient inverse problem for elliptic equation 421

8.5.5 Program

The gradient method described above was realized in its simplest version with constant

iteration parameter s

k+1

(simple iteration method). In the perturbation procedure for

input data, a situation with perturbed boundary conditions of the second kind (8.167)

is modeled.

Program PROBLEM19

C PROBLEM19 - IDENTIFICATION OF THE LOWEST COEFFICIENT

C IN THE ELLIPTIC EQUATION OF SECOND ORDER

C TWO-DIMENSIONAL PROBLEM

IMPLICIT REAL

8 ( A-H, O-Z )

PARAMETER ( DELTA = 0.02D0, N1 = 51, N2 = 51 )

DIMENSION A(12

N2), X1(N1), X2(N2), CK(N2), GR(N2)

+ ,YG1(N1), YG2(N2), YG3(N1), YG4(N2)

+ ,YD1(N1), YD2(N2), YD3(N1), YD4(N2)

+ ,YY1(N1), YY2(N2), YY3(N1), YY4(N2)

COMMON / SB5 / IDEFAULT(4)

COMMON / CONTROL / IREPT, NITER

C PARAMETERS:

C X1L, X2L - COORDINATES OF THE LEFT CORNER;

C X1R, X2R - COORDINATES OF THE RIGHT CORNER;

C N1, N2 - NUMBER OF NODES IN THE SPATIAL GRID;

C H1, H2 - SPATIAL STEP OF THE GRID;

C TAU - TIME STEP;

C DELTA - INPUT-DATA INACCURACY LEVEL;

C CK(N2) - COEFFICIENT TO BE RECONSTRUCTED;

C YG1(N1),

C YG2(N2),

C YG3(N1),

C YG4(N2) - MESH FUNCTION AT THE NEAR-BOUNDARY NODES;

C YD1(N1),

C YD2(N2),

C YD3(N1),

C YD4(N2) - DISTURBED MESH FUNCTION AT THE NEAR-BOUNDARY NODES;

C EPSR - RELATIVE ACCURACY FOR THE DIFFERENCE-PROBLEM SOLUTION;

C EPSA - ABSOLUTE ACCURACY FOR THE DIFFERENCE-PROBLEM SOLUTION;

C EQUIVALENCE ( A(1), A0 ),

( A(N+1), A1 ),

( A(2

N+1), A2 ),

( A(9

N+1), F ),

( A(10

N+1), U );

( A(11

N+1), US )

X1L = 0.D0

X1R = 1.D0

X2L = 0.D0

X2R = 1.D0

EPSR = 1.D-6

EPSA = 1.D-9

SS = - 40.D0

422 Chapter 8 Other problems

OPEN (01, FILE=’RESULT.DAT’) ! FILE TO STORE THE CALCULATED DATA

C GRID

H1 = (X1R-X1L) / (N1-1)

H2 = (X2R-X2L) / (N2-1)

DOI=1,N1

X1(I) = X1L + (I-1)

END DO

DOJ=1,N2

X2(J) = X2L + (J-1)

END DO

N=N1

DOI=1,12

A(I) = 0.D0

END DO

C DIRECT PROBLEM

C BOUNDARY CONDITIONS

CALL BNG (A(10

N+1), X1, X2, N1, N2)

C LOWEST COEFFICIENT

DOJ=1,N2

CK(J) = AC(X2(J))

END DO

C DIFFERENCE-SCHEME COEFFICIENT IN THE DIRECT PROBLEM

CALL FDS (A(1), A(N+1), A(2

N+1), CK, A(9

N+1), A(10

N+1),

+ H1, H2, N1, N2)

C SOLUTION OF THE DIFFERENCE PROBLEM

IDEFAULT(1) = 0

IREPT = 0

CALL SBAND5 (N, N1, A(1), A(10

N+1), A(9

N+1), EPSR, EPSA)

C OBSERVATIONAL DATA

CALL BNGDER (A(10

N+1), H1, H2, N1, N2, YG1, YG2, YG3, YG4)

C DISTURBING OF MEASURED VALUES

DOI=2,N1-1

H=H2

IF (I.EQ.2 .OR. I.EQ.N1-1) H = (H1+H2) / 2.D0

YD1(I) = YG1(I) + 2.

DELTA

(RAND(0)-0.5)

YD3(I) = YG3(I) + 2.

DELTA

(RAND(0)-0.5)

END DO

DOJ=2,N2-1

H=H1

IF (J.EQ.2 .OR. J.EQ.N2-1) H = (H1+H2) / 2.D0

YD2(J) = YG2(J) + 2.

DELTA

(RAND(0)-0.5)

YD4(J) = YG4(J) + 2.

DELTA

(RAND(0)-0.5)

END DO

Section 8.5 Coefﬁcient inverse problem for elliptic equation 423

C INVERSE PROBLEM

C ITERATION METHOD

IT=0

C INITIAL APPROXIMATION

DOJ=1,N2

CK(J) = 0.D0

END DO

100IT=IT+1

C GROUND STATE

C DIFFERENCE-SCHEME COEFFICIENTS IN THE DIRECT PROBLEM

CALL FDS (A(1), A(N+1), A(2

N+1), CK, A(9

N+1), A(10

N+1),

+ H1, H2, N1, N2)

C SOLUTION OF THE DIFFERENCE PROBLEM

IDEFAULT(1) = 0

IREPT = 0

CALL SBAND5 (N, N1, A(1), A(10

N+1), A(9

N+1), EPSR, EPSA)

C SOLUTION AT THE OBSERVATION POINTS

CALL BNGDER (A(10

N+1), H1, H2, N1, N2, YY1, YY2, YY3, YY4)

C CONJUGATE STATE

C DIFFERENCE-SCHEME COEFFICIENTS

CALL FDS (A(1), A(N+1), A(2

N+1), CK, A(9

N+1), A(10

N+1),

+ H1, H2, N1, N2)

C RIGHT-HAND SIDE

CALL RHS (A(9

N+1), N1, N2, H1, H2, YD1, YD2, YD3, YD4,

+ YY1, YY2, YY3, YY4)

C SOLUTION OF THE DIFFERENCE PROBLEM

IDEFAULT(1) = 0

IREPT = 0

CALL SBAND5 (N, N1, A(1), A(11

N+1), A(9

N+1), EPSR, EPSA)

C SIMPLE-ITERATION METHOD

C FUNCTIONAL GRADIENT

CALL GRAD (A(10

N+1), A(11

N+1), H1, H2, N1, N2, GR)

C NEXT APPROXIMATION

DOJ=2,N2-1

CK(J) = CK(J) + SS

GR(J)

END DO

424 Chapter 8 Other problems

C EXIT FROM THE ITERATION PROCEDURE BY THE DISCREPANCY CRITERION

SUM = 0.D0

SUMD = 0.D0

DOI=3,N1-2

SUM = SUM + (YD1(I) - YY1(I))

SUM = SUM + (YD3(I) - YY3(I))

END DO

SUMD = SUMD + 2.D0

(N1-4)

DELTA

DOJ=3,N2-2

SUM = SUM + (YD2(J) - YY2(J))

SUM = SUM + (YD4(J) - YY4(J))

END DO

SUMD = SUMD + 2.D0

(N2-4)

DELTA

SUM = SUM + (YD1(2) - YY1(2))

(H1+H2) / 2.D0

SUM = SUM + (YD4(2) - YY4(2))

(H1+H2) / 2.D0

SUM = SUM + (YD3(2) - YY3(2))

(H1+H2) / 2.D0

SUM = SUM + (YD3(N1) - YY3(N1))

(H1+H2) / 2.D0

SUMD = SUMD + 2.D0

(H1+H2)

DELTA

IF ( SUM.GT.SUMD ) GO TO 100

C SOLUTION

WRITE ( 01,

) (A(10

N+I), I=1,N)

WRITE ( 01,

) (CK(J), J=1,N2)

CLOSE (01)

STOP

END

DOUBLE PRECISION FUNCTION AC ( X2 )

IMPLICIT REAL

8 ( A-H, O-Z )

C LOWEST COEFFICIENT IN THE ELLIPTIC EQUATION

AC = 10.D0

RETURN

END

SUBROUTINE FDS (A0, A1, A2, CK, F, U, H1, H2, N1, N2)

C GENERATION OF THE COEFFICIENT ARRAY FOR THE DIFFERENCE SCHEME

C FOR THE ELLIPTIC EQUATION

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION A0(N1,N2), A1(N1,N2), A2(N1,N2)

+ ,CK(N2), F(N1,N2), U(N1,N2)

DOJ=2,N2-1

DOI=2,N1-1

A1(I-1,J) = 1.D0/(H1

H1)

A1(I,J) = 1.D0/(H1

H1)

A2(I,J-1) = 1.D0/(H2

H2)

A2(I,J) = 1.D0/(H2

H2)

A0(I,J) = A1(I,J) + A1(I-1,J) + A2(I,J) + A2(I,J-1) + CK(J)

F(I,J) = 0.D0

END DO

Section 8.5 Coefﬁcient inverse problem for elliptic equation 425

C FIRST-KIND HOMOGENEOUS BOUNDARY CONDITIONS

DOJ=2,N2-1

A0(1,J) = 1.D0

A1(1,J) = 0.D0

A2(1,J) = 0.D0

F(1,J) = U(1,J)

F(2,J) = F(2,J) + U(1,J) / (H1

H1)

END DO

DOJ=2,N2-1

A0(N1,J) = 1.D0

A1(N1-1,J) = 0.D0

A1(N1,J) = 0.D0

A2(N1,J) = 0.D0

F(N1,J) = U(N1,J)

F(N1-1,J) = F(N1-1,J) + U(N1,J) / (H1

H1)

END DO

DOI=2,N1-1

A0(I,1) = 1.D0

A1(I,1) = 0.D0

A2(I,1) = 0.D0

F(I,1) = U(I,1)

F(I,2) = F(I,2) + U(I,1) / (H2

H2)

END DO

DOI=2,N1-1

A0(I,N2) = 1.D0

A1(I,N2) = 0.D0

A2(I,N2) = 0.D0

A2(I,N2-1) = 0.D0

F(I,N2) = U(I,N2)

F(I,N2-1) = F(I,N2-1) + U(I,N2) / (H2

H2)

END DO

A0(1,1) = 1.D0

A1(1,1) = 0.D0

A2(1,1) = 0.D0

F(1,1) = U(1,1)

A0(N1,1) = 1.D0

A2(N1,1) = 0.D0

F(N1,1) = U(N1,1)

A0(1,N2) = 1.D0

A1(1,N2) = 0.D0

F(1,N2) = U(1,N2)

A0(N1,N2) = 1.D0

F(N1,N2) = U(N1,N2)

RETURN

END

SUBROUTINE RHS (F, N1, N2, H1, H2, YD1, YD2, YD3, YD4,

+ YY1, YY2, YY3, YY4)

C RIGHT-HAND SIDE IN THE EQUATION FOR THE CONGUGATE STATE