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

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

346 Chapter 8 Other problems

8.1.1 Statement of the problem

Among inverse mathematical physics problems, of primary signiﬁcance for practical

applications is the boundary value inverse problem. This problem is often encountered

in diagnostics, in the cases in which it is required to reconstruct, from additional mea-

surements made inside the calculation domain, the thermal boundary condition at the

domain boundary, where direct measurements are unfeasible.

This problem belongs to the class of conditionally well-posed problems and, for its

approximate solution, development of special regularization methods is under way. A

general approach for the solution of unstable problems for partial equations is the gen-

eralized inverse method. This method uses some perturbation of the initial equation,

the problem for the perturbed equation being a well-posed one. Here, the perturbation

parameter serves as regularization parameter.

In the consideration of the boundary value inverse problem for the one-dimensional

parabolic equation of the second order, the generalized inverse method can be devel-

oped considering the initial problem as a problem for the evolutionary equation of the

ﬁrst order. A second possibility is related with the consideration of the boundary value

inverse problem as a problem with initial data for the evolutionary equation of the sec-

ond order. Here, as the evolutionary variable, the spatial variable is used. That is why

here we speak of the continuation over the spatial variable in a boundary value inverse

problem.

Consider a heat conduction boundary value inverse problem in which it is required

to continue the solution over the spatial variable, which serves as time variable. The

problem is to be constructed as follows. The solution v(x, t) is to be found from the

equation

∂v

∂t

−

∂

∂x

= 0, 0 < x < l, 0 < t < T, (8.1)

supplemented with initial conditions, written in terms of the variables x and t, of the

form

v(0, t) = ϕ(t), 0 < t < T, (8.2)

∂v

∂x

(0, t) = 0, 0 < t < T, (8.3)

v(x, 0) = 0, 0 < x < l. (8.4)

Let us apply in the boundary value inverse problem (8.1)–(8.4) the following change

of variables: the variable x is changed for t, and t, for x (l → T , T → l). Next,

we denote the solution to be found as u(x, t) ( = v(t, x)).Foru(x, t), we obtain the

Section 8.1 Continuation over spatial variable in boundary value inverse problems 347

problem

∂

∂t

−

∂u

∂x

= 0, 0 < x < l, 0 < t < T, (8.5)

u(0, t) = 0, 0 < t < T, (8.6)

u(x, 0) = u

(x), 0 < x < l, (8.7)

∂u

∂t

(x, 0) = 0, 0 < x < l. (8.8)

In the new settings, the function ϕ(t) refers to u

(x).

The problem (8.5)–(8.8) is written as the operator equation

− Au = 0 (8.9)

with the initial conditions (8.7), (8.8). The operator A is deﬁned by

Au =

∂u

∂x

(8.10)

with the domain of deﬁnition

D(A) ={u | u = u(x, t), x ∈ [0, l], u(0, t) = 0}.

The introduced operator A is not a self-adjoint, sign-deﬁnite operator in H = L

(0, l).

8.1.2 Generalized inverse method

To approximately solve the inverse problem (8.7)–(8.9), we use the generalized inverse

method. We use this method in its version in which the approximate solution u

(x, t)

is to be determined from the perturbed equation

− Au

+ αA

∗

= 0, (8.11)

supplemented with the initial conditions

(x, 0) = u

(x), 0 < x < l, (8.12)

∂u

∂t

(x, 0) = 0, 0 < x < l. (8.13)

By virtue of (8.10), the operator conjugate in H = L

(0, l) to A is given by

∗

u =−

∂u

∂x

, (8.14)

and, in addition,

D(A

∗

) ={u | u = u(x, t), x ∈ [0, l], u(l, t) = 0}.

348 Chapter 8 Other problems

With regard to (8.10), (8.14), in the variant (8.11)–(8.13) of the generalized inverse

method it is required to solve the problem

∂

∂t

−

∂u

∂x

+ α

∂

∂x

= 0, 0 < x < l, 0 < t < T, (8.15)

(0, t) = 0, 0 < t < T, (8.16)

∂u

∂x

(l, t) = 0, 0 < t < T, (8.17)

(x, 0) = u

(x), 0 < x < l, (8.18)

∂u

∂t

(x, 0) = 0, 0 < x < l. (8.19)

Thus, the generalized inverse method leads us to a hyperbolic perturbation (see (8.15))

of the initial parabolic equation (8.5), this perturbation being a well-known one in the

computational practice.

Theorem 8.1 For the solution of the boundary value problem (8.15)–(8.19), there

holds the following a priori estimate:



∂u

∂t

(x, t)



+ α



∂u

∂x

(x, t)



≤ α exp



√





∂u

∂x

(x)



. (8.20)

Proof. Let us prove the above statement under more general conditions (namely, for

problem (8.11)–(8.13)). We multiply equation (8.11) scalarwise by du/dt; this yields:





+ αAu









. (8.21)

For the right-hand side of (8.21), we have:





≤

√





+ αAu





. (8.22)

Substitution of (8.22) into (8.21) yields the estimate



+ αAu



≤ α exp



√



Au



. (8.23)

Inequality (8.23) yields the following simpler estimate:

Au

≤α exp



1/4



Au

.

The desired estimate (8.20) for the perturbed problem (8.15)–(8.19) is a simple corol-

lary to the estimate (8.23).

Section 8.1 Continuation over spatial variable in boundary value inverse problems 349

Some other variants of the generalized inverse method can also be used. For in-

stance, we can determine the approximate solution from the equation

− Au

+ αAA

∗

= 0, (8.24)

supplemented with conditions (8.12), (8.13).

Hence, with (8.10) and (8.14), the approximate solution u

(x, t) is to be found from

the equation

∂

∂t

−

∂u

∂x

+ α

∂

∂x

∂t

= 0, 0 < x < l, 0 < t < T (8.25)

and conditions (8.16)–(8.19). Similarly to Theorem 8.1, we formulate the following

statement:

Theorem 8.2 For the solution of the boundary value problem (8.16)–(8.19), (8.25),

there holds the a priori estimate

u

(x, t)



∂u

∂t

(x, t)



≤ exp



1 + 2α

2α



u

(x)

. (8.26)

Proof. Here again, the consideration will be performed for the general perturbed evo-

lutionary equation. We can conveniently reformulate equation (8.24) as a system of

ﬁrst-order equations.

We deﬁne the vector U ={u

, u

} and the space H

as the direct sum of spaces H:

= H ⊕H. The addition in H

is performed coordinatewise, and the scalar product

is deﬁned as

(U, V ) = (u

) + (u

We deﬁne u

= u

, u

= du

/dt and write equation (8.24) as a system of ﬁrst-order

equations (U

={u

, u

}):

+ PU

= 0, (8.27)

where

P =

0 −E

−A αAA

∗

. (8.28)

Equation (8.27) is supplemented (see (8.12), (8.13)) with the initial conditions

(0) = U

={u

, 0}. (8.29)

In the notation introduced, we have:

U



=u



+u



=u





350 Chapter 8 Other problems

For problem (8.27)–(8.29), there holds the a priori estimate

U

(t)

≤ exp



1 + 2α

2α



U

(0)

. (8.30)

With (8.28), we have:

PU ={−u

, −Au

+ αAA

∗

We multiply equation (8.27) scalarwise by U

and obtain:

(u



+u



) + αA

∗



= (Au

, u

) + (u

, u

). (8.31)

For the right-hand side, we use the estimate

(Au

, u

) ≤ αA

∗



4α

u



, u

) ≤ εu



4ε

u



(8.32)

Choosing ε = 1/2 and substituting (8.32) into (8.31), we arrive at the estimate (8.30).

From (8.30), the estimate (8.26) for the solution of problem (8.16)–(8.19), (8.25)

follows.

8.1.3 Difference schemes for the generalized inverse method

Consider now the proposed variants of the generalized inverse method on the mesh

level. For simplicity, we restrict ourselves to the case of a one-dimensional boundary

value problem. To pass to a difference problem, we introduce the uniform grid

¯ω ={x | x = x

= ih, i = 0, 1,...,N, Nh = l},

where ω is the set of internal nodes, and ∂ω is the set of boundary nodes. The operators

A and A

∗

are introduced in the traditional manner. Let in the subscriptless notation we

have

Ay =



0, i = 0,

¯x

, i = 1, 2,...,N

with the domain of deﬁnition

D(A) ={y | y(x), x ∈¯ω, y

= 0}.

We deﬁne the scalar product and the norm in the mesh Hilbert space H = L

(ω) as

(z, y) =



x∈¯ω

z(x)y(x)h, y=



(y, y).

Section 8.1 Continuation over spatial variable in boundary value inverse problems 351

Then, for the adjoint operator A

∗

we have

∗

y =



−y

, i = 0, 1,...,N − 1,

0, i = N ,

and

D(A

∗

) ={y | y(x), x ∈¯ω, y

= 0}.

For the operator A

∗

A we have:

∗

Ay =

⎧

⎨

⎩

−y

/ h, i = 0,

−y

¯xx

, i = 1, 2,...,N − 1,

0, i = N .

Let us dwell ﬁrst on the variant (8.11) of the generalized inverse method. We determine

the difference solution of problem (8.11)–(8.13) from the equation

− Ay + α A

∗

Ay = 0 (8.33)

for y(x, t) ∈ H , supplemented with the initial conditions

y(x, 0) = u

(x), x ∈¯ω, (8.34)

(x, 0) = 0, x ∈¯ω. (8.35)

To approximately solve problem (8.33)–(8.35), we use a time-uniform grid with the

time step size τ . Consider the difference scheme

n+1

− 2y

+ y

n−1

−

A(y

+ y

n−1

) +

∗

A(y

n+1

+ 2y

+ y

n−1

) = 0,

n = 1, 2,.... (8.36)

To perform a stability study for the scheme (8.36), we introduce the settings

+ y

n−1

), w

− y

n−1

). (8.37)

Then, the difference scheme (8.36) can be written as

n+1

− w

− Av

∗

A(v

n+1

+ v

) = 0. (8.38)

We multiply the difference equation (8.38) scalarwise by

2(v

n+1

− v

) = y

n+1

− y

n−1

= τ(w

n+1

+ w

);

then we obtain

w

n+1



−w



+ αAv

n+1



− αAv



= τ(A v

n+1

+ w

). (8.39)

352 Chapter 8 Other problems

The right-hand side in (8.39) can be estimated as follows:

(Av

n+1

+ w

) ≤ βAv



4β

w

n+1

+ w



≤ βAv



1 + ε

4β

w

n+1



1 + ε

4βε

w



We choose ε = 1/2; then, we substitute the latter inequality into (8.39). This yields



1 −



w

n+1



+αAv

n+1



≤



1 +



w



+α



1 +



Av



. (8.40)

Using in the case of τ ≤ 4β/3 the estimate



1 −



−1

≤ exp



4β



and choosing β = (3α/2)

1/2

, from (8.40) we obtain the a priori estimate



1 −



w

n+1



+ αAv

n+1



≤ 



1 −



w



+ αAv





(8.41)

with

 = exp





2α



. (8.42)

Thus, the following statement is proved:

Theorem 8.3 For the difference scheme (8.36), in the case of τ ≤ 2(2α/3)

1/2

there

holds the estimate (8.37), (8.40).

In the variant (8.24) of the generalized inverse method, it is required to solve the equa-

tion

− Ay + α AA

∗

= 0

with the initial conditions (8.34), (8.35). In the numerical realization, we use the

scheme

n+1

− 2y

+ y

n−1

− Ay

+ α AA

∗

n+1

− y

n−1

2τ

= 0,

n = 1, 2,....

(8.43)

Theorem 8.4 The difference scheme (8.43) is -stable with

 = exp



1 + 4α

4α



and, for this scheme, there holds the a priori estimate

y





n+1

− y



≤ 



y

n−1





− y

n−1





. (8.44)

Section 8.1 Continuation over spatial variable in boundary value inverse problems 353

Proof. To obtain the desired a priori estimate, we scalarwise multiply the difference

equation (8.43) by

◦

n+1

− y

n−1

2τ

and, using the standard subscriptless notation adopted in the theory of difference

schemes, arrive at the equality

, y

◦

) − (Ay , y

◦

) + αA

∗

◦



= 0. (8.45)

For the ﬁrst two terms in the left side of (8.45) we have

, y

◦

) =

, y

)

(Ay , y

◦

) ≤ αA

∗

◦



4α

y

(8.46)

We insert (8.46) into (8.45); then, we arrive at the inequality

, y

)

≤

2α

y

. (8.47)

We add to the both sides of (8.47) the term

(y

)

(y



−y

n−1



)



+ y

n−1

− y

n−1



= 2(y

n−1

, y

) + τ y



Using the fact that

y

=y

n−1

+ τ y



≤ (1 + ε)y

n−1





1 +

4ε



y



and taking the inequality

n−1

, y

) ≤ βy

n−1



4β

y



into account, from (8.47) we obtain the inequality

(y

n−1



+y



)

≤



1 + ε

2α

+ 2β



y

n−1





2β

+ τ +

2α

1 + 4ε

4ε



y



, (8.48)

that holds for all positive ε and β.

354 Chapter 8 Other problems

We choose ε = τ/2 and β = 1/2; then, we obtain

1 + ε

2α

+ 2β = 1 +

2α

4α

(

− 1),

2β

+ τ +

2α

1 + 4ε

4ε

= 1 + τ +

4α

2α

(

− 1).

Here, the value of  is deﬁned in the conditions of the theorem. With regard to these

inequalities, inequality (8.48) yields the desired estimate (8.44).

8.1.4 Program

We numerically solve the boundary value inverse problem (8.1)–(8.4). In the frame-

work of the quasi-real computational experiment, to formulate the boundary conditions

(8.2), we consider the direct problem in which it is required to determine the function

v(x, t) from equation (8.1), initial condition (8.4), boundary condition (8.3) at the left

boundary, and the condition

v(l, t) = ψ(t), 0 < t < T

at the right boundary. To ﬁnd the approximate solution, we use the purely implicit

difference scheme.

The approximate solution of the inverse problem can be found using the variant

(8.24) of the generalized inverse method. In terms of the initial variables, this corre-

sponds to the solution of the problem (see (8.16)–(8.19), (8.25))

∂

∂x

−

∂v

∂t

+ α

∂

∂t

∂x

= 0, 0 < x < l, 0 < t < T,

(0, t) = ϕ(t)), 0 < t < T,

∂v

∂x

(0, t) = 0, 0 < t < T,

(x, 0) = 0, 0 < x < l,

∂v

∂t

(x, T ) = 0, 0 < x < l.

This boundary value problem is solved by the difference scheme (8.43).

The value of α is determined using the discrepancy criterion. We assume that the

input-data inaccuracy (the inaccuracy in setting ϕ(t))) plays a decisive role and, hence,

the approximation inaccuracy in choosing the regularization parameter can be ignored.

Program PROBLEM15

C PROBLEM15 - IDENTIFICATION OF THE BOUNDARY CONDITION

C NON-STATIONARY ONE-DIMENSIONAL PROBLEM

Section 8.1 Continuation over spatial variable in boundary value inverse problems 355

IMPLICIT REAL

8 ( A-H, O-Z )

PARAMETER ( DELTA = 0.01D0, N = 101, M = 101 )

DIMENSION X(N), Y(N,M), F(N,M), YY(M)

+ ,FI(M), FID(M), FIY(M), Q(M), QA(M)

+ ,A(M), B(M), C(M), FF(M)!M>=N

C PARAMETERS:

C XL, XR - LEFT AND RIGHT ENDS OF THE SEGMENT;

C N - NUMBER OF GRID NODES OVER SPACE;

C TMAX - MAXIMAL TIME;

C M - NUMBER OF GRID NODES OVER TIME;

C DELTA - INPUT-DATA INACCURACY LEVEL;

C FI(M) - EXACT DIFFERENCE BOUNDARY CONDITION;

C FID(M) - DISTURBED DIFFERENCE BOUNDARY CONDITION;

C Q(M) - EXACT SOLUTION OF THE INVERSE PROBLEM

C (BOUNDARY CONDITION);

C QA(M) - APPROXIMATE SOLUTION OF THE INVERSE PROBLEM;

XL = 0.D0

XR = 1.D0

TMAX = 1.D0

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

C GRID

H=(XR-XL)/(N-1)

DOI=1,N

X(I) = XL + (I-1)

END DO

TAU = TMAX / (M-1)

C DIRECT PROBLEM

C BOUNDARY REGIME

DOK=1,M

T = (K-1)

TAU

Q(K) = AF(T, TMAX)

END DO

C INITIAL CONDITION

T = 0.D0

DOI=1,N

Y(I,1) = 0.D0

END DO

C NEXT TIME LAYER

DOK=2,M

T=T+TAU

C DIFFERENCE-SCHEME COEFFICIENTS

C PURELY IMPLICIT SCHEME

DOI=2,N-1

A(I) = 1.D0 / (H