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

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

176 Chapter 6 Right-hand side identiﬁcation

Here, special computational algorithms for stable numerical differentiation must

be applied. To approximate the solution, we can use ﬁnite-difference regularization

algorithms in which the discretization step size serves the regularization parameter.

This method was discussed at length previously, when we considered the stationary

problem for the second-order ordinary differential equation. In the solution of the

inverse problem (6.46)–(6.49), the discretization step sizes over space and time must

be matched with the inaccuracy in u(x, t).

Evolutionary problems possess some speciﬁc feature of utmost signiﬁcance. The

current solution of the problem depends only on the prehistory, i.e., on the solutions at

preceding times, and does not depends on the solutions at future times. This point can

(and often must) be taken into account when developing computational algorithms.

This requirement is perfectly justiﬁed in the case of direct problems of type (6.45)–

(6.47) and can also be important when one considers inverse problems similar to prob-

lem (6.46)–(6.49).

Discussing the solution of evolutionary problems in the general context, we can

speak of computational algorithms of two types. In the computational algorithms of

the ﬁrst type, the current solution is to be calculated from the solution at preceding

times. Here, we speak of local algorithms for solving evolutionary problems. In the

global algorithms the current solution is sought using a procedure that involves future

times.

6.2.2 Global regularization

To approximately solve the inverse problem (6.46)–(6.49), we use the general regu-

larization scheme proposed by A. N. Tikhonov. First of all, we introduce some new

settings.

In the Hilbert space L

(), in the ordinary way we introduce the norm and the

scalar product:

(v, w) =



v(x)w(x) dx, v

= (v, v) =



(x) dx.

For functions v(x, t), w(x, t) ∈ H, where H = L

), we set:

(v, w)

∗



(v, w) dt =



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

∗

= ((v, v)

∗

)

1/2

We deﬁne the operator

Au =−



k(x)



on the set of functions satisfying the condition (6.44). In L

(), for the operator A we

have:

A = A

∗

≥ mE, m = κ

Section 6.2 Right-hand side identiﬁcation in the case of parabolic equation 177

Then, the inverse problem (6.46)–(6.48) can be written as

f = D u, (6.50)

where the operator

Du =

+ Au (6.51)

is deﬁned on the set of functions satisfying the initial condition (6.47). For the input

data (see (6.49)), we have:

u

− u

∗

≤ δ. (6.52)

In the Tikhonov regularization method, the approximate solution f

of problem

(6.50)–(6.52) is sought as the solution of the following variational problem:

(E + αD

∗

D) f

= Du

. (6.53)

A speciﬁc feature of the regularization method as applied to the solution of the inverse

evolutionary problem is manifested in the problem operator D deﬁned by (6.51). In

particular, it becomes necessary to explicitly deﬁne the operator D

∗

Provided that v(x, 0) = 0 and w(x, T ) = 0, we have:

(Dv,w)

∗







dt +



(A,w)dt

= (v, w)



−







dt +



(A,w)dt = (v, Dw)

∗

Thus, we deﬁne the operator D

∗

w =−

+ Aw (6.54)

on the set of functions

w(x, T ) = 0, 0 ≤ x ≤ l. (6.55)

In view of (6.47), (6.51), (6.54) and (6.55), equation (6.53) for the approximate

solution f

is an elliptic equation involving second time derivatives and forth spatial

derivatives. For f

, the boundary conditions at t = 0, T are as follows:

(x, 0) = 0, 0 ≤ x ≤ l, (6.56)

D f

(x, T ) = 0, 0 ≤ x ≤ l. (6.57)

These points should be given particular attention in the numerical realization.

In the case of the identiﬁcation problem for the right-hand side of the second-order

ordinary differential equation, the following two basic possibilities are available. The

ﬁrst possibility implies using the Tikhonov regularization method in interpreting the

right-hand side identiﬁcation problem as a problem in which it is required to solve a

178 Chapter 6 Right-hand side identiﬁcation

ﬁrst-kind operator equation. In the second possibility, the right-hand side identiﬁcation

problem is considered as a problem in which it is required to calculate the values of a

bounded operator. The latter possibility was realized in the scheme (6.52), (6.53) for

problem (6.46)–(6.49). It also makes sense to dwell here on the standard variant of the

Tikhonov regularization.

With the given right-hand side f (x, t), the solution of the boundary value problem

(6.45)–(6.47) uniquely deﬁnes the solution u(x, t). To reﬂect this correspondence, we

introduce the operator G:

G f = u. (6.58)

Instead of u(x, t), the function u

(x, t) is given and, in addition, the estimate (6.52)

holds.

The approximate solution f

of problem (6.49), (6.58) can be found as the solution

of the following problem:

( f

) = min

v∈H

(v), (6.59)

where

(v) = (Gv − u



∗

)

+ α(v

∗

)

. (6.60)

The following important circumstance deserves mention. In using algorithm (6.52),

(6.53), constructed around the interpretation of the identiﬁcation problem as an un-

bounded operator value problem, additional constraints on the sought function (bound-

ary conditions of type (6.56), (6.57)) are to be posed. This is not quite justiﬁed a pro-

cedure (see the results of numerical experiments on the identiﬁcation of the right-hand

side of the ordinary differential equation in the previous section). No such problems

are encountered in using (6.59), (6.60).

6.2.3 Local regularization

In the approximate solution of the problem in which it is required to identify the right-

hand side of a non-stationary equation from known solution, one can more conve-

niently employ an algorithm that makes it possible to determine the right-hand side

at a given time using input information only at preceding times. Compared to the

global regularization algorithm, here we, speaking generally, loose in the approximate-

solution accuracy, but save time. Consider some basic possibilities in the construction

of local regularization algorithms for the approximate solution of the inverse problem

(6.46)–(6.49). We will dwell here on the local analogue of the regularization of type

(6.59), (6.60).

The basic idea uses the fact that the right-hand side is deﬁned by the solution by

each ﬁxed moment. In other words, the numerical differentiation procedure is to be

regularized only with respect to spatial variables. In fact, the input data are to be

smoothed considering only part of all variables. Such algorithms can be realized em-

ploying preliminary discretization over time.

Section 6.2 Right-hand side identiﬁcation in the case of parabolic equation 179

We introduce the uniform grid over time

¯ω

= ω

∪{T }={t

= nτ, n = 0, 1,...,N

,τN

= T }.

The following subscriptless notation generally accepted in the theory of difference

schemes will be used:

y = y

, ˆy = y

n+1

, ˇy = y

n−1

y −ˇy

, y

ˆy − y

Let us formulate the inverse right-hand side identiﬁcation problem for the parabolic

equation after the partial discretization. For the main quantities, the same settings as

in the continuous case will be used. For simplicity, we restrict the consideration to the

case of purely explicit time approximation, in which the right-hand side is determined

(compare (6.48)) from the differential-difference relation

− u

n−1

+ Au

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

. (6.61)

The input data (direct-problem solution u

) is set accurate to some inaccuracy. We

assume the inaccuracy level to be characterized by a constant δ, so that

u

− u

≤δ, n = 1, 2,...,N

. (6.62)

We denote the approximate solution of problem (6.61), (6.62) at the moment t = t

as f

. Consider the matter of reconstruction of the function f

from given u

, u

n−1

and f

n−1

. We assume that the approximate right-hand side f

corresponds to the

solution of the boundary value problem

− w

n−1

+ Aw

= f

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

, (6.63)

= 0. (6.64)

First, we determine f

as the solution of the problem

( f

) = min

v∈H

(v), (6.65)

where H = L

() and

(v) =u

− w(v)

+ αv

. (6.66)

Here, in view of (6.63) and (6.64), w(v) is the solution of the difference problem

w −ˇw

+ Aw = v. (6.67)

180 Chapter 6 Right-hand side identiﬁcation

This allows us to write equation (6.66) as

(v) =



−

ˇw − G



+ αv

where

= G

∗



E + A



−1

In this way, the approximate solution at each time t = t

can be determined from

the equation

∗

+ α f

= G

∗

−

∗

ˇw. (6.68)

With (6.67), using the introduced notation, we obtain the following problem for w

−1

w = f

ˇw. (6.69)

6.2.4 Iterative solution of the identiﬁcation problem

In solving inverse mathematical physics problems, primary attention should be paid to

iteration methods that most fully realize the idea of ﬁnding the solution of the inverse

problem through successive solution of a set of direct problems. Considering approxi-

mate solution of problem (6.46)–(6.49) by iteration methods, here we dwell on global

regularization in the variant (6.58)–(6.60).

After symmetrization of (6.58), we write the two-layer iteration method as

k+1

− f

k+1

+ G

∗

G f

= G

∗

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

In using the steepest descend method, the iteration parameters can be calculated by

k+1



r



∗

Gr



∗



, r

= G

∗

G f

− G

∗

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

The number of iterations in (6.70), (6.71) is to be matched with the inaccuracy δ (see

(6.52)).

This approach can be of use provided that we are able to calculate the operators G

and G

∗

. Recall that v = G f

is the solution of the direct problem

+ Av = f

, 0 < t ≤ T, (6.72)

v(0) = 0. (6.73)

The values of the conjugate operator w = G

∗

v can be found by solving the direct

problem

−

+ Aw = v, 0 ≤ t < T , (6.74)

w(T ) = 0. (6.75)

Section 6.2 Right-hand side identiﬁcation in the case of parabolic equation 181

Thus, for a given iteration parameter, the passage to the next iteration is related, in

line with (6.70), with the solution of two direct problems, namely, (6.72), (6.73) and

(6.74), (6.75).

Consider some speciﬁc features in the computational realization of the iteration

method related, ﬁrst of all, with time discretization. We retain for the mesh functions

the same settings as for the continuous-argument functions.

We use a uniform grid ¯ω over the interval

 = [0, l] with a grid size h:

¯ω ={x | x = x

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

Here, ω is the set of internal nodes, and ∂ω is the set of boundary nodes.

We approximate the differential operator A at internal nodes accurate to the second

order with the difference operator

Ay =−(ay

¯x

)

, x ∈ ω, (6.76)

where, for instance, a(x) = k(x − 0.5h).

In the mesh Hilbert space L

(ω), we introduce the norm by the relation y=

(y, y)

1/2

, where

(y,w) =



x∈ω

y(x)w(x)h.

Recall that on the set of functions vanishing on ∂ω, for the self-conjugate operator A

under the constraints k(x) ≥ κ>0, q(x) ≥ 0 there holds the estimate

A = A

∗

≥ κλ

E, (6.77)

where

sin

πh

≥

We put into correspondence to the direct problem (6.72), (6.73) the following sym-

metric difference problem:

− v

n−1

A(v

n+1

+ v

) =

( f

n+1

+ f

n = 1, 2,...,N

(6.78)

= 0, x ∈ ω. (6.79)

In the latter case, the approximation inaccuracy is of second order both over time and

space. In a similar manner, other two-layer difference schemes can be considered. We

bring the problem (6.78), (6.79) to the operator notation v = Gf

, which deﬁnes the

operator G.

For two-dimensional mesh functions, we deﬁne a Hilbert space H = L

) in

which the scalar product and the norm are deﬁned as

(v, w)

∗

−1



n=1

)τ +

) +

), v

∗



(v, v)

∗

182 Chapter 6 Right-hand side identiﬁcation

A problem conjugate in H to (6.78), (6.79) is the difference problem (see (6.74),

(6.75))

−

− w

n−1

A(w

+ w

n−1

) =

+ v

n−1

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

,(6.80)

= 0, x ∈ ω. (6.81)

We can check this by scalarwise multiplying equation (6.78) by w

. The problem

(6.80), (6.81) can be written in a compact form as w = G

∗

In line with (6.70), the iteration method can be written as

k+1

− f

k+1

+ G

∗

= G

∗

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

At the ﬁrst stage, we have to calculate the right-hand side G

∗

; to this end, a

boundary value problem of type (6.80), (6.81) is to be solved. For the difference

= G

∗

− G

∗

to be calculated, it is required at each step to solve two boundary

value problems, ((6.78), (6.79) and (6.80), (6.81)). The iteration parameters can be

calculated (see (6.71)) by the formula

k+1



r



∗

Gr



∗



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

Determination of Gr

requires an additional boundary value problem of type (6.78),

(6.79) to be solved.

The regularization parameter here is the number K (δ) of iterations in (6.82). The

criterion for the exit from the iterative process is

Gf

K (δ)

− u



∗

≤ δ. (6.84)

It should be noted that the algorithm of interest is to be used under certain restric-

tions imposed on the right-hand side. In accordance with the applied symmetrization

at the expense of the operator G, the following constraints on the right-hand side are

imposed:

= 0, x ∈ ω,

= 0, n = 0, 1,...,N , x ∈ ∂ω.

The algorithm (6.78)–(6.84) is embodied in the program PROBLEM6. Given below

is the complete text of this program.

Program PROBLEM6

C PROBLEM6 - RIGHT-HAND SIDE IDENTIFICATION

C ONE-DIMENSIONAL NON-STATIONARY PROBLEM

Section 6.2 Right-hand side identiﬁcation in the case of parabolic equation 183

IMPLICIT REAL

8 ( A-H, O-Z )

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

DIMENSION X(N), Y(N,M), YD(N,M), F(N,M), FA(N,M)

+ ,V(N,M), W(N,M), GU(N,M), RK(N,M), YY(N)

+ ,A(N), B(N), C(N), D(N), E(N), FF(N)

C PARAMETERS:

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

C N - NUMBER OF GRID NODES OVER THE SPATIAL VARIABLE;

C TMAX - MAXIMAL TIME;

C M - NUMBER OF GRID NODES OVER TIME;

C DELTA - INPUT-DATA INACCURACY;

C Y(N,M) - EXACT DIFFERENCE SOLUTION OF THE BOUNDARY-VALUE

C PROBLEM;

C YD(N,M) - DISTURBED DIFFERENCE SOLUTION OF THE BOUNDARY-VALUE

C PROBLEM;

C F(N,M) - EXACT RIGHT-HAND SIDE;

C FA(N,M) - CALCULATED RIGHT-HAND SIDE;

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 INITIAL CONDITION

T = 0.D0

DOI=1,N

Y(I,1) = 0.D0

YD(I,1) = 0.D0

END DO

C NEXT TIME LAYER

DOK=2,M

T=T+TAU

C DIFFERENCE-SCHEME COEFFICIENTS

DOI=2,N-1

X1 = (X(I) + X(I-1)) / 2

X2 = (X(I+1) + X(I)) / 2

A(I) = AK(X1) / (2

B(I) = AK(X2) / (2

C(I) = A(I) + B(I) + 1.D0 / TAU

FF(I) = A(I)

Y(I-1,K-1)

184 Chapter 6 Right-hand side identiﬁcation

+ + (1.D0 / TAU - A(I) - B(I) )

Y(I,K-1)

+ + B(I)

Y(I+1,K-1)

+ + (AF(X(I),T) + AF(X(I),T-TAU)) / 2

END DO

C BOUNDARY CONDITIONS AT THE LEFT AND RIGHT END POINTS

B(1) = 0.D0

C(1) = 1.D0

FF(1) = 0.D0

A(N) = 0.D0

C(N) = 1.D0

FF(N) = 0.D0

C SOLUTION AT THE NEXT TIME LAYER

ITASK = 1

CALL PROG3 ( N, A, C, B, FF, YY, ITASK )

DOI=1,N

Y(I,K) = YY(I)

END DO

C NOISE ADDITION

DOK=2,M

YD(1,K) = 0.D0

YD(N,K) = 0.D0

DO I = 2,N-1

YD(I,K) = Y(I,K)

+ + 2.D0

DELTA

(RAND(0)-0.5D0) / DSQRT(TMAX

(XR-XL))

END DO

C SYMMETRIZATION (SOLUTION OF THE CONJUGATE PROBLEM)

C INITIAL CONDITION

T = TMAX

DOI=1,N

GU(I,M) = 0.D0

END DO

C NEXT TIME LAYER

DO K = M-1, 1, -1

T=T-TAU

C DIFFERENCE-SCHEME COEFFICIENTS

DOI=2,N

X1 = (X(I) + X(I-1)) / 2

X2 = (X(I+1) + X(I)) / 2

A(I) = AK(X1) / (2

B(I) = AK(X2) / (2

Section 6.2 Right-hand side identiﬁcation in the case of parabolic equation 185

C(I) = A(I) + B(I) + 1.D0 / TAU

FF(I) = A(I)

GU(I-1,K+1)

+ + (1.D0 / TAU - A(I) - B(I) )

GU(I,K+1)

+ + B(I)

GU(I+1,K+1)

+ + (YD(I,K) + YD(I,K+1)) / 2

END DO

C BOUNDARY CONDITIONS AT THE LEFT AND RIGHT END POINTS

B(1) = 0.D0

C(1) = 1.D0

FF(1) = 0.D0

A(N) = 0.D0

C(N) = 1.D0

FF(N) = 0.D0

C SOLUTION OF THE PROBLEM AT THE NEXT TIME LAYER

ITASK = 1

CALL PROG3 ( N, A, C, B, FF, YY, ITASK )

DOI=1,N

GU(I,K) = YY(I)

END DO

C INVERSE RIGHT-HAND SIDE IDENTIFICATION PROBLEM

IT=0

ITMAX = 1000

C INITIAL APPROXIMATION

DOK=1,M

DOI=1,N

FA(I,K) = 0.D0

END DO

100 IT = IT+1

C SOLUTION OF THE DIRECT PROBLEM FOR THE GIVEN RIGHT-HAND SIDE

C INITIAL CONDITION

T = 0.D0

DOI=1,N

V(I,1) = 0.D0

END DO

C NEXT TIME LAYER

DOK=2,M

T=T+TAU

C DIFFERENCE-SCHEME COEFFICIENTS

DOI=2,N

X1 = (X(I) + X(I-1)) / 2

X2 = (X(I+1) + X(I)) / 2