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

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

376 Chapter 8 Other problems

To approximately solve the boundary value inverse problem (8.91)–(8.93), (8.95),

(8.96), we will use the iteration method related with reﬁnement of a boundary con-

dition of type (8.94).

8.3.2 Iteration method

We write the boundary value inverse problem as an operator equation of the ﬁrst kind.

To the boundary condition (8.94) we put in correspondence equation (8.96), i. e.,

Aμ = ϕ. (8.97)

The linear operator A is deﬁned for functions given on some portion of the boundary

γ , and its values are functions given on the other portion of the boundary (on γ

∗

). To

calculate the values of A, we solve the boundary value problem (8.91)–(8.95).

To approximately solve the ill-posed problem (8.97), we employ an iteration method

based on the passage to a problem with a symmetric operator, where, instead of (8.97),

we solve the equation

∗

Aμ = A

∗

ϕ. (8.98)

In using the explicit iteration method for (8.98), we have:

k+1

− μ

k+1

− A

∗

Aμ

= A

∗

ϕ, k = 0, 1,.... (8.99)

For non-stationary functions given on γ and γ

∗

, we deﬁne the Hilbert spaces

L(γ, [0, T ]) and L(γ

∗

, [0, T ]) with the scalar products and norms

(u,v) =



u(x)v(x) dx dt, u=



(u, u),

(u,v)

∗



∗

u(x)v(x) dx dt, u

∗



(u, u)

∗

respectively. In the steepest descend method, the iteration parameters can be calculated

by the rule

k+1

r



Ar



∗

, r

= A

∗

Aμ

− A

∗

ϕ. (8.100)

Let us dwell on a point of primary concern in the practical use of (8.99) related with

the necessity to calculate the values of the operator conjugate in A. For functions μ

and ν given respectively on γ and γ

∗

we have:

(Aμ, ν)

∗

= (μ, A

∗

ν). (8.101)

We obtain the values of A as

Aμ = u(x, t), x ∈ γ

∗

Section 8.3 Identiﬁcation of the boundary condition in two-dimensional problems 377

Here, u(x, t) is the solution of the boundary value problem (8.91)–(8.95) (the function

u(x, t) deﬁnes the ground state of the system). The values of the conjugate operator

can be found from the function ψ(x, t), to be found from the solution of an auxiliary

boundary value problem (conjugate state). For this boundary value problem to be for-

mulated, we multiply equation (8.91) through by ψ(x, t), and then perform integration

over  and time:





∂u

∂t

+ Lu



ψ dx dt = 0, (8.102)

where

Lu =−



β=1

∂

∂x



k(x)

∂u

∂x



For the ﬁrst term in (8.102), we have



∂u

∂t

ψ dx dt =−



∂ψ

∂t

dx dt, (8.103)

provided that, in view of (8.95),

ψ(x, T ) = 0, x ∈ . (8.104)

Similar manipulations with the second term in (8.102), performed with due regard

for the boundary conditions (8.92)–(8.94), yield:



Luψ dx dt

=−



μψ dx dt +



∂

∂ψ

∂n

dx dt +



uLψ dx dt. (8.105)

We choose the boundary conditions for the conjugate state in the form

k(x)

∂ψ

∂n

(x, t) = ν(x

, t), x ∈ ∂γ

∗

, (8.106)

k(x)

∂ψ

∂n

(x, t) = 0, x ∈ , (8.107)

k(x)

∂u

∂ψ

(x, t) = 0, x ∈ γ. (8.108)

With such boundary conditions, substitution of (8.103), (8.105) yields the equality





−

∂ψ

∂t

+Lψ



dx dt −



μψ dx dt +



∗

uν dx dt = 0. (8.109)

We assume that the function ψ(x, t) satisﬁes the equation

−

∂ψ

∂t

+ Lψ = 0, x ∈ , 0 < t < T. (8.110)

378 Chapter 8 Other problems

Figure 8.11 Calculation grid over space

Thereby, the conjugate state is to be found as the solution of the direct problem (8.104),

(8.106)–(8.108), (8.110).

Under the present conditions, from (8.109) we obtain:



μψ dx dt =



∗

uν dx dt.

The latter equality can be brought to the form (8.101), and for the values of A

∗

obtain the representation

∗

ν = ψ(x, t), x ∈ γ,

where ψ(x, t) is the solution of the boundary value problem (8.104), (8.106)–(8.108),

(8.110).

8.3.3 Difference problem

Let us start the present consideration with the formulation of a difference problem for

the direct problem (8.91)–(8.95). Taking the fact into account that on the calculation-

domain boundary second-kind boundary conditions are given, we use a grid whose

nodes are displaced by half the step size from the domain boundary (see Figure 8.11).

Section 8.3 Identiﬁcation of the boundary condition in two-dimensional problems 379

Along each of the directions x

,β= 1, 2, the grid is uniform, and let



| x

= (i

− 1/2)h

, i

= 1, 2,...,N

, N

= l

For the grid in the rectangle  we use the settings ω = ω

× ω

With regard to the boundary conditions (8.92)–(8.94), to the differential operator L

we put in correspondence a two-dimensional difference operator deﬁned as the sum of

one-dimensional operators:

 = 

+ 

. (8.111)

With the grid displaced by half the grid size, it seems reasonable to deﬁne the differ-

ence operator 

as follows:



y =

⎧

⎨

⎩

−a

+ h

, x

/ h

, x

= h

/2,

−(a

¯x

)

, h

/2 < x

< l

− h

/2,

¯x

/ h

, x

= l

− h

/2,

where, for instance, a

, x

) = k(x

− 0.5h

, x

). In a similar way, we deﬁne the

difference operator 



y =

⎧

⎨

⎩

−a

, x

+ h

/ h

, x

= h

/2,

−(a

¯x

)

, h

/2 < x

< l

− h

/2,

¯x

/ h

, x

= l

− h

with a

, x

) = k(x

, x

− 0.5h

In the Hilbert space L

(ω) of mesh functions set on ω, we deﬁne the scalar product

(y,v)



x∈ω

y(x)v(x)h

It is an easy matter to check that the difference operator  deﬁned by (8.111) is a

self-adjoint operator non-negative in L

(ω), i.e.  = 

∗

≥ 0.

We denote as y

the difference solution at the time t

= nτ , where τ>0 is the time

step size (N

τ = T ). An approximate solution of the direct problem (8.91)–(8.95) can

be found as the solution of the difference problem

n+1

− y

+ (σ y

n+1

+ (1 − σ)y

) = f

n = 0, 1,...,N

− 1,

(8.112)

= 0, x ∈ ω. (8.113)

The weighted difference scheme (8.112), (8.113) is absolutely stable in the case

of σ ≥ 1/2. For simplicity, we restrict ourselves to the case of the purely implicit

scheme (σ = 1), in which for the solution of the difference problem (8.112), (8.113)

there holds the following stability estimate:

y

n+1

≤y

+τ  f

, n = 0, 1,...,N

− 1.

380 Chapter 8 Other problems

The inhomogeneity of the right-hand side of (8.112) results from the boundary con-

ditions (8.94), so that in the case of σ = 1 we can put:

(x) =

⎧

⎨

⎩

0, x

∈ ω

, x

= h

/2,

0, x

∈ ω

, h

/2 < x

< l

− h

/2,

n+1

)/ h

, x

∈ ω

, x

= l

− h

/2.

(8.114)

Thus, the boundary conditions are included on the operator level into the difference

equation. With the form of the inhomogeneous right-hand side taken into account, the

above estimate proves that the solution of the difference problem is stable with respect

to the initial data and boundary conditions.

Let us formulate now the boundary value inverse problem on the mesh level. In

the case under consideration (see (8.91)–(8.93), (8.95), (8.96)) the function μ

n = 1, 2,...,N

in the right-hand side of (8.114) is not given. Instead of this function,

the solution at the nodes adjacent to the boundary γ

∗

is assumed to be known:

, 0.5h

) = ϕ

), x

∈ ω

, n = 1, 2,...,N

. (8.115)

Hence, on the mesh level we can speak about a problem in which it is required to iden-

tify the right-hand side of special form (8.114) from boundary observations (8.115).

8.3.4 Iterative reﬁnement of the boundary condition

Similarly to the continuous case, here we use the formulation of the boundary value

inverse problem in the form of a ﬁrst-kind operator equation. From the given boundary

condition (right-hand side set in the form (8.114) with known mesh function μ

n = 1, 2,...,N

), we put in correspondence the function ϕ

), n = 1, 2,...,N

(see (8.115)):

Aμ = ϕ. (8.116)

To calculate ϕ

), n = 1, 2,...,N

, we have to solve the difference boundary value

problem (8.112)–(8.114). The operator equation (8.116) is the difference analogue

to (8.97).

The iteration method is to be applied to the symmetrized equation (see (8.98))

∗

Aμ = A

∗

ϕ. (8.117)

The explicit iteration method for the approximate solution of (8.117) can be written in

the form

k+1

− μ

k+1

+ A

∗

Aμ

= A

∗

ϕ, k = 0, 1,... . (8.118)

The matter of calculating the values of A

∗

deserves special consideration. To this

end, ﬁrst of all we have to give mesh analogues of the Hilbert spaces L(γ , [0, T ]) and

Section 8.3 Identiﬁcation of the boundary condition in two-dimensional problems 381

L(γ

∗

, [0, T ]), to be subsequently referred to as H and H

∗

, respectively. In H and H

∗

the scalar product and the norm are deﬁned as

(u,v) =



n=1



∈ω

, l

τ, u=



(u, u),

(u,v)

∗



n=1



∈ω

, 0)v

, 0)h

τ, u

∗



(u, u)

∗

Then, in the steepest descend method, for the iteration parameters in (8.118) we have:

k+1

r



Ar



∗

, r

= A

∗

Aμ

− A

∗

ϕ. (8.119)

The operator adjoint to A can be found from the equality

(Aμ, ν)

∗

= (μ, A

∗

ν). (8.120)

Here, the mesh functions μ

), n = 1, 2,...,N

and ν

), n = 1, 2,...,N

are

given in the calculation domain ω at the nodes adjacent to the boundaries γ and γ

∗

respectively.

To formulate the difference problem for the conjugate state, we multiply the differ-

ence equation for the ground state

n+1

− y

+ y

n+1

= f

, (8.121)

by the mesh function ψ

n+1

τ and sum it up over all nodes x ∈ ω and all n =

0, 1,...,N

− 1:



n=1



x∈ω



− y

n−1

+ y



τ =



n=1



x∈ω

n−1

τ. (8.122)

With regard to (8.114) and with regard to the introduced notation, for the right-hand

side of (8.122) we immediately obtain:



n=1



x∈ω

n−1

τ = (μ, ψ ). (8.123)

As a result, we can relate this right-hand side with the right-hand side in (8.120).

Using the initial condition (8.113), we can rearrange the ﬁrst term in the left-hand

side in (8.122) as



n=1



x∈ω

− y

n−1

τ =−



n=1



x∈ω

n+1

− ψ

τ, (8.124)

382 Chapter 8 Other problems

provided that the following additional condition holds:

= 0, x ∈ ω. (8.125)

With regard to the self-adjointness of  and L

(ω), for the second term in the right-

hand side of (8.122) we have:



n=1



x∈ω

y

τ =



n=1



x∈ω

ψ

τ. (8.126)

Substitution of (8.123), (8.124), and (8.126) into (8.122) yields:



n=1



x∈ω



−

n+1

− ψ

+ ψ



τ = (μ, ψ ). (8.127)

We choose the mesh function ψ

(x) so that the left side of (8.127) is related with the

left side in (8.120).

We assume that the function ψ

(x) satisﬁes the difference equation

−

n+1

− ψ

+ ψ

= g

, n = 1, 2,...,N

(8.128)

supplemented with the initial condition (8.125). Hence, with the given right-hand side

(x) the conjugate state is to be found using the purely implicit scheme (8.125),

(8.128) on the grid displaced by the time step size.

From (8.127) and (8.128), we obtain



n=1



x∈ω

τ = (μ, ψ ). (8.129)

In view of (8.115) and (8.116), the left-hand side of (8.129) coincides with the left-

hand side in (8.120) provided that the right-hand side in (8.128) is given as

(x) =

⎧

⎨

⎩

)/ h

, x

∈ ω

, x

= h

/2,

0, x

∈ ω

, h

/2 < x

< l

− h

/2,

0, x

∈ ω

, x

= l

− h

/2.

(8.130)

This choice of the right-hand side corresponds to the solution of the difference bound-

ary value problem with boundary conditions of the second kind (see (8.106)).

From the equality between the right-hand sides of (8.120) and (8.129), for the values

of A

∗

we obtain:

∗

= ψ

, l

), x

∈ ω

, n = 1, 2,...,N

. (8.131)

Here, the mesh function ψ

(x) is to be found as the solution of the difference problem

(8.125), (8.128), (8.130).

Section 8.3 Identiﬁcation of the boundary condition in two-dimensional problems 383

8.3.5 Program realization

The program presented below embodies the above iteration solution method for the

boundary value inverse problem (8.91)–(8.93), (8.95), (8.96) in the case of k(x) = 1.

The iteration method is terminated based on discrepancy.

Program PROBLEM17

C PROBLEM17 - BOUNDARY-VALUE INVERSE PROBLEM

C TWO-DIMENSIONAL PROBLEM

C ITERATIVE REFINEMENT OF THE BOUNDARY CONDITION

IMPLICIT REAL

8 ( A-H, O-Z )

PARAMETER ( DELTA = 0.005D0, N1 = 50, N2 = 50, M = 101 )

DIMENSION A(9

N2), X1(N1), X2(N2), T(M)

+ ,Y(N1,N2), F(N1,N2), G(N1)

+ ,FI(N1,M), FID(N1,M), FS(N1,M), FIK(N1,M)

+ ,U(N1,M),UA(N1,M), UK(N1,M), R(N1,M), AR(N1,M)

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 - MESH SIZE OVER SPACE;

C TAU - TIME STEP;

C DELTA - INPUT-DATA INACCURACY LEVEL;

C FI(N1,M) - EXACT DIFFERENCE BOUNDARY CONDITION;

C FID(N1,M) - DISTURBED DIFFERENCE BOUNDARY CONDITION;

C U(N1,M) - EXACT SOLUTION OF THE INVERSE PROBLEM

C (BOUNDARY CONDITION);

C UA(N1,M) - APPROXIMATE SOLUTION OF THE INVERSE PROBLEM;

C EPSR - RELATIVE SOLUTION INACCURACY IN THE DIFFERENCE

C PROBLEM;

C EPSA - ABSOLUTE SOLUTION INACCURACY IN THE DIFFERENCE

C PROBLEM;

C EQUIVALENCE ( A(1), A0 ),

( A(N+1), A1 ),

( A(2

N+1), A2 ),

X1L = 0.D0

X1R = 1.D0

X2L = 0.D0

X2R = 0.5D0

TMAX = 1.D0

PI = 3.1415926D0

EPSR = 1.D-5

EPSA = 1.D-8

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

C GRID

H1 = (X1R-X1L) / N1

384 Chapter 8 Other problems

H2 = (X2R-X2L) / N2

TAU = TMAX / (M-1)

DOI=1,N1

X1(I) = X1L + (I-0.5D0)

END DO

DOJ=1,N2

X2(J) = X2L + (J-0.5D0)

END DO

DOK=1,M

T(K) = (K-1)

TAU

END DO

N=N1

DOI=1,9

A(I) = 0.0

END DO

C DIRECT PROBLEM

C PURELY IMPLICIT DIFFERENCE SCHEME

C INITIAL CONDITION

DOI=1,N1

DOJ=1,N2

Y(I,J) = 0.D0

END DO

C BOUNDARY CONDITION

CALL FLUX (U, X1, T, N1, M)

DOK=2,M

C DIFFERENCE-SCHEME COEFFICIENTS IN THE DIRECT PROBLEM

DOI=1,N1

G(I) = U(I,K)

END DO

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

N+1), F, Y,

+ H1, H2, N1, N2, TAU, G)

C SOLUTION OF THE DIFFERENCE PROBLEM

IDEFAULT(1) = 0

IREPT = 0

CALL SBAND5 (N, N1, A(1), Y, F, EPSR, EPSA)

C INPUT DATA FOR THE INVERSE PROBLEM

DOI=1,N1

FI(I,K) = Y(I,1)

END DO

C DISTURBING OF MEASURED QUANTITIES

DOI=1,N1

DOK=2,M

FID(I,K) = FI(I,K) + 2.

DELTA

(RAND(0)-0.5)

Section 8.3 Identiﬁcation of the boundary condition in two-dimensional problems 385

END DO

C INVERSE PROBLEM

C RIGHT SIDE OF THE SYMMETRIZED EQUATION

C DIRECT PROBLEM

C PURELY IMPLICIT DIFFERENCE SCHEME

C INITIAL CONDITION (ATt=T)

DOI=1,N1

DOJ=1,N2

Y(I,J) = 0.D0

END DO

DO K = M,2,-1

C DIFFERENCE-SCHEME COEFFICIENTS IN THE DIRECT PROBLEM

DOI=1,N1

G(I) = FID(I,K)

END DO

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

N+1), F, Y,

+ H1, H2, N1, N2, TAU, G)

C SOLUTION OF THE DIFFERENCE PROBLEM

IDEFAULT(1) = 0

IREPT = 0

CALL SBAND5 (N, N1, A(1), Y, F, EPSR, EPSA)

C RIGHT-HAND SIDE

DOI=1,N1

FS(I,K) = Y(I,N2)

END DO

C ITERATION METHOD

IT=0

C INITIAL APPROXIMATION

DOI=1,N1

DOK=2,M

UK(I,K) = 0.D0

END DO

100IT=IT+1

C GROUND STATE

C INITIAL CONDITION

DOI=1,N1

DOJ=1,N2