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

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

276 Chapter 7 Evolutionary inverse problems

In a similar manner, difference schemes for other variants of the generalized inverse

method can be constructed. Let us examine, for instance, two-layer difference schemes

for a pseudo-parabolic perturbation. To approximately solve equation (7.91) with the

initial condition (7.82), we use the difference scheme

n+1

− y

− (σ y

n+1

+ (1 − σ)y

) + α

n+1

− y

= 0. (7.134)

The scheme (7.134) can be written in the canonical form (7.85) with

A =−, B = E + (α − σ τ ). (7.135)

Theorem 7.16 The difference scheme (7.85), (7.135) is ρ-stable in H for all τ>0

with

ρ = 1 +

τ, (7.136)

provided that σ ≤ 0.

Proof. In the case under consideration, the condition B = B

∗

> 0 and the right

inequality of (7.96) obviously hold. With (7.135), the left inequality of (7.96) assumes

the form

c(E + (α − σ τ )) ≥ 

and, hence, we can put c = 1/α; then, for ρ we have (7.136).

Remark 7.17 The estimate of the difference solution of scheme (7.134) with ρ given

by (7.136) is perfectly consistent with the estimate for the solution of the differential

problem (see estimate (7.126)).

Remark 7.18 Accurate to the adopted notation, the scheme (7.134) coincides with the

ordinary weighted scheme immediately written for (7.82), (7.91):

n+1

− y

− (σ



n+1

+ (1 − σ



) = 0. (7.137)

In sufﬁces to put in (7.137) σ



= σ − α/τ. The only essential thing here is that the

weight σ



in (7.137) is negative.

The presented schemes, as well as all other known schemes, can be obtained by

regularization of unstable difference schemes. Moreover, based on an analysis of reg-

ularized difference schemes, other versions of the generalized inverse method can be

constructed.

The difference scheme (7.134) of the pseudo-parabolic generalized inverse method

is nothing else than a regularized scheme of type (7.85), (7.88), (7.97). The difference

scheme of the basic variant (7.130) of the generalized inverse method refers to the

use of regularized difference schemes with perturbed operators A and B (see (7.131)).

Section 7.2 Regularized difference schemes 277

The regularized difference scheme (7.85), (7.89), (7.97) can be used to construct a new

version of the generalized inverse method in which the approximate solution is sought

as the solution of a Cauchy problem for the equation

− Au

+ αA

= 0. (7.138)

In construction of regularized difference schemes, the perturbed differential prob-

lem itself becomes unnecessary. Note also that for difference problems the construc-

tion of a stable scheme around the regularization principle is facilitated substantially

because general results are available concerning the necessary and sufﬁcient conditions

for ρ-stability of the difference schemes.

7.2.4 Regularized additive schemes

Certain difﬁculties may arise in the computational realization of difference schemes

constructed within the framework of the generalized inverse method. The latter is

related, ﬁrst of all, with the fact that the perturbed problem is a singularly perturbed

problem (with small perturbation parameters at higher derivatives). It should be noted

that a perturbed equation is a higher-order equation. To discuss this matter, here we

imply that the standard variant of the generalized inverse method (7.103), (7.104) is

used.

In the approximate solution of problem (7.77)–(7.79), it seems reasonable to use the

difference scheme

n+1

− y

− y

+ α

(σ y

n+1

+ (1 − σ)y

) = 0. (7.139)

The latter implies that in the two-parametric scheme (7.130) we have σ

= 0 and

= σ . The scheme (7.139) is ρ-stable (see Theorem 7.15) if σ ≥ 1/2.

Above, we showed that the difference schemes constructed within the framework of

the generalized inverse method can be obtained in two ways. The ﬁrst approach implies

approximation of the perturbed differential problem. The second approach is based on

using the regularization principle for operator-difference schemes. In connection with

the present consideration of scheme (7.139), another interpretation of the generalized

inverse method scheme can prove useful, in which this scheme is considered as a

scheme with smoothed mesh solution.

We consider scheme (7.139) as a scheme of the predictor-corrector type. Consider

the simplest case of σ = 1. In this case, at the predictor stage we seek the mesh

solution ˜y

n+1

using the explicit scheme

˜y

n+1

− y

− y

= 0. (7.140)

In accordance with (7.139), in the case of σ = 1 the corrector stage is

n+1

−˜y

n+1

+ α

n+1

= 0. (7.141)

278 Chapter 7 Evolutionary inverse problems

The explicit scheme (7.140) is unstable; a stable solution can be obtained using

(7.141). Equation (7.141) can be considered as an Euler equation in the smoothing

problem for the mesh function ˜y

n+1

) = min

v∈H

(v), (7.142)

(v) =v −˜y

n+1



+ ταv

. (7.143)

In the latter interpretation, we can speak of the generalized inverse method in the form

(7.140), (7.141) (or (7.140), (7.142), (7.143)) as of a local regularization algorithm.

Such a relation can be traced especially distinctly on the difference level.

Realization of the difference scheme (7.139) requires solution of the following

forth-order difference elliptic equation:

ασ

v +

v = f.

At present, this task presents rather a difﬁcult computational problem. It would there-

fore be desirable to have regularization difference schemes simpler from the standpoint

of computational realization. In this connection, additive spatial-variable split schemes

deserve mention in which transition to the next time layer necessitates solution of one-

dimensional difference problems.

We construct additive operator-difference schemes starting from scheme (7.139).

We begin with the interpretation of scheme (7.139) as a regularization scheme. As a

generating scheme, we consider, as usual, an explicit scheme in which for the general-

ized inverse method we have:

n+1

− y

+ (α

− )y

= 0. (7.144)

To an additive regularization scheme, the following form of (7.139) corresponds:

(E + αστ

)

n+1

− y

+ (α

− )y

= 0.

In the latter case, the transition from (7.144) is interpreted as

B = E → B + αστ

We can change scheme (7.139) to the following somewhat different yet equivalent

form:

n+1

− y

+ (E + αστ

)

−1

(α

− )y

= 0. (7.145)

Scheme (7.145) corresponds to the use of the following version of multiplicative reg-

ularization:

A = α

−  → (E + αστ

)

−1

A. (7.146)

Section 7.2 Regularized difference schemes 279

It is the regularization (7.145), (7.146) that will be used as the basis for the desired

regularized additive operator-difference schemes.

In view of (7.80), we have the following additive representation for :

 =



β=1



,

y =−(a

¯x

)

,β= 1, 2. (7.147)

Either of these operator terms is positive,



= 

∗

> 0,β= 1, 2,

and both are generally non-permutable:



= 



In the case of three-dimensional problems, the operator  can be factorized to three

one-dimensional pairwise non-permutable positive operators.

We construct the additive scheme using constructions analogous to (7.145), and

used for each operator term in (7.147). The latter yields the regularized scheme

n+1

− y



β=1

(E + αστ

)

−1

(α

− 

= 0. (7.148)

To obtain the stability condition, we write scheme (7.148) in the form

n+1

− y



β=1

= 0, (7.149)

where

= (E + αστ

)

−1

(α

− 

), β = 1, 2. (7.150)

In the case of interest, for the operators we have A

= A

∗

, β = 1, 2; hence,

A =



β=1

= A

∗

. (7.151)

The scheme (7.149), (7.151) is ρ-stable provided that the two-sided inequality

(7.96) is fulﬁlled; in the case of interest, this inequality assumes the form

1 − ρ

E ≤ A ≤

1 + ρ

E. (7.152)

Inequality (7.152) will be fulﬁlled, for instance, in the case of

1 − ρ

2τ

E ≤ A

≤

1 + ρ

2τ

E,β= 1, 2.

280 Chapter 7 Evolutionary inverse problems

With allowance for (7.150), we rewrite this inequality as follows:

1 − ρ

2τ

(E + αστ

) ≤ α

− 

≤

1 + ρ

2τ

(E + αστ

). (7.153)

Taking into account the inequality

α

− 

≥−

4α

we obtain that for the left inequality in (7.153) to be fulﬁlled it sufﬁces for us to put

ρ = 1 +

2α

. (7.154)

The right inequality (7.153) is fulﬁlled in all cases if σ ≥ 1. This allows the following

statement to be formulated:

Theorem 7.19 The regularized additive scheme (7.148) is ρ-stable in H, with ρ de-

ﬁned by (7.154) with τ>0, provided that σ ≥ 1.

The realization of the additive scheme (7.148) requires the inversion of the one-

dimensional difference operators E +αστ

, β = 1, 2. The following computational

scheme can be proposed which demonstrates an intimate interrelation between the

regularized schemes under consideration and additive-averaged difference schemes.

We determine the auxiliary mesh functions y

(β)

n+1

, β = 1, 2 by solving the equations

(β)

n+1

− y

2τ

+ (E + αστ

)

−1

(α

− 

= 0,β= 1, 2. (7.155)

Afterwards, the next-layer solution is given by

n+1



β=1

(β)

n+1

. (7.156)

We change equation (7.155) to the form analogous to (7.139):

(β)

n+1

− y

2τ

− 

+ α

(σ y

(β)

n+1

+ (1 − σ)y

) = 0.

Thus, the determination of y

(β)

n+1

, β = 1, 2 corresponds to the case in which the gener-

alized inverse method is applied to each operator term in (7.147).

Above (see (7.140), (7.141)), a close relation between the generalized inverse

method and local regularization algorithms (smoothing at each time layer, see (7.140),

(7.142), (7.143)) was noted. In the latter interpretation, the regularized scheme (7.148)

corresponds to the case in which, as a local regularization algorithm, we use smoothing

along individual directions (see (7.155)) followed by averaging according to (7.156).

Section 7.2 Regularized difference schemes 281

We can relate with the regularized additive scheme (7.148) or, more precisely, with

the additive-averaged scheme (7.155), (7.156), the generalized inverse method used to

solve problem (7.101), (7.102) with

A =



β=1

, A

= A

∗

≥ 0,β= 1, 2,..., p. (7.157)

In the case of two-component splitting, we have p = 2. The approximate solution of

problem (7.101), (7.102), (7.157) is to be found as the solution of a Cauchy problem

for the equation

−



β=1

+ α



β=1

= 0. (7.158)

For the solution of problem (7.104), (7.158), the following estimate of stability with

respect to initial data holds:

u

(t)≤exp



4α



u

.

The latter estimate is perfectly consistent with the above ρ-stability estimate (see

(7.154)) for the regularized difference scheme (7.148).

7.2.5 Program

To approximately solve the inverted-time problem, we use the regularized additive

scheme (7.148). The computational realization here is based on the interpretation

of this scheme as the additive-averaged scheme (7.155), (7.156). The ﬁve-diagonal

difference problems are solved by the ﬁve-point sweep algorithm using the subroutine

PROG5. A detailed description of the algorithm is given in Section 6.1.

The input data for the inverted-time problem are taken from the solution of the direct

problem (7.77)–(7.79). To ﬁnd the approximate solution, the following purely implicit

additive-averaged scheme is used:

(β)

n+1

− y

2τ

+ 

n+1

= 0,β= 1, 2,

n+1



β=1

(β)

n+1

The end-time difference solution of the direct problem is perturbed to use this solu-

tion as input data in the inverse problem.

The value of the regularization parameter was chosen considering the discrepancy,

using the sequence

= α

, q > 0

282 Chapter 7 Evolutionary inverse problems

with given α

= 0. and q = 0.75. To calculate the discrepancy, we additionally solved

the direct problem in which initial data were taken equal to the solution of the inverse

problem.

For simplicity, in the program presented below we have restricted ourselves to the

simplest case of time-independent coefﬁcients, with k(x) = 1 in (7.77), and the

weighting parameter in the regularized additive scheme (7.155), (7.156) is σ = 1.

Program PROBLEM11

C PROBLEM11 INVERTED-TIME PROBLEM

C TWO-DIMENSIONAL PROBLEM

C ADDITIVE REGULARIZED SCHEME

IMPLICIT REAL

8 ( A-H, O-Z )

PARAMETER ( DELTA = 0.01, N1 = 101, N2 = 101, M = 101 )

DIMENSION U0(N1,N2), UT(N1,N2), UTD(N1,N2), U(N1,N2)

+ ,X1(N1), X2(N2), Y(N1,N2), Y1(N1,N2), Y2(N1,N2), YY(N1)

+ ,A(N1), B(N1), C(N1), D(N1), E(N1), F(N1) ! N1 >= N2

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 - STEP OVER SPACE;

C TAU - TIME STEP;

C DELTA - INPUT-DATA INACCURACY;

C Q - MULTIPLIER IN THE REGULARIZATION PARAMETER;

C U0(N1,N2) - RECONSTRUCTED INITIAL CONDITION;

C UT(N1,N2) - END-TIME SOLUTION OF THE DIRECT PROBLEM;

C UTD(N1,N2)- DISTURBED SOLUTION OF THE DIRECT PROBLEM;

C U(N1,N2) - APPROXIMATE SOLUTION OF THE INVERSE PROBLEM;

X1L = 0.D0

X1R = 1.D0

X2L = 0.D0

X2R = 1.D0

TMAX = 0.025D0

PI = 3.1415926D0

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

C GRID

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

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

TAU = TMAX / (M-1)

DOI=1,N1

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

END DO

DOJ=1,N2

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

END DO

C DIRECT PROBLEM

Section 7.2 Regularized difference schemes 283

C ADDITIVE-AVERAGE DIFFERENCE SCHEME

C INITIAL CONDITION

DOI=1,N1

DOJ=1,N2

U0(I,J) = AU(X1(I),X2(J))

Y(I,J) = U0(I,J)

END DO

DOK=2,M

C SWEEP OVER X1

DOJ=2,N2-1

C COEFFICIENTS OF THE DIFFERENCE SCHEME IN THE DIRECT PROBLEM

B(1) = 0.D0

C(1) = 1.D0

F(1) = 0.D0

A(N1) = 0.D0

C(N1) = 1.D0

F(N1) = 0.D0

DOI=2,N1-1

A(I) = 1.D0 / (H1

H1)

B(I) = 1.D0 / (H1

H1)

C(I) = A(I) + B(I) + 0.5D0 / TAU

F(I) = 0.5D0

Y(I,J) / TAU

END DO

C SOLUTION OF THE DIFFERENCE PROBLEM

ITASK1 = 1

CALL PROG3 ( N1, A, C, B, F, YY, ITASK1 )

DOI=1,N1

Y1(I,J) = YY(I)

END DO

C SWEEP OVER X2

DOI=2,N1-1

C COEFFICIENTS OF THE DIFFERENCE SCHEME IN THE DIRECT PROBLEM

B(1) = 0.D0

C(1) = 1.D0

F(1) = 0.D0

A(N2) = 0.D0

C(N2) = 1.D0

F(N2) = 0.D0

DOJ=2,N2-1

A(J) = 1.D0 / (H2

H2)

B(J) = 1.D0 / (H2

H2)

C(J) = A(J) + B(J) + 0.5D0 / TAU

284 Chapter 7 Evolutionary inverse problems

F(J) = 0.5D0

Y(I,J) / TAU

END DO

C SOLUTION OF THE DIFFERENCE PROBLEM

ITASK1 = 1

CALL PROG3 ( N2, A, C, B, F, YY, ITASK1 )

DOJ=1,N2

Y2(I,J) = YY(J)

END DO

C ADDITIVE AVERAGING

DOI=1,N1

DOJ=1,N2

Y(I,J) = 0.5D0

(Y1(I,J) + Y2(I,J))

END DO

C DISTURBING OF MEASURED QUANTITIES

DOI=1,N1

DOJ=1,N2

UT(I,J) = Y(I,J)

UTD(I,J) = Y(I,J)

END DO

DOI=2,N1-1

DOJ=2,N2-1

UTD(I,J) = UT(I,J) + 2.

DELTA

(RAND(0)-0.5)

END DO

C INVERSE PROBLEM

C ITERATIVE PROCESS TO ADJUST

C THE VALUE OF THE REGULARIZATION PARAMETER

IT=0

ITMAX = 100

ALPHA = 0.001D0

Q = 0.75D0

100IT=IT+1

C ADDITIVE-AVERAGED REGULARIZATION SCHEMES

C INITIAL CONDITION

DOI=1,N1

DOJ=1,N2

Y(I,J) = UTD(I,J)

END DO

DOK=2,M

Section 7.2 Regularized difference schemes 285

C SWEEP OVER X1

DOJ=2,N2-1

C COEFFICIENTS OF THE DIFFERENCE SCHEME IN THE INVERSE PROBLEM

DO I = 2,N1-1

A(I) = ALPHA / (H1

B(I) = 4.D0

ALPHA / (H1

C(I) = 6.D0

ALPHA / (H1

4) + 0.5D0 / TAU

D(I) = 4.D0

ALPHA / (H1

E(I) = ALPHA / (H1

F(I) = 0.5D0

Y(I,J) / TAU

+ - (Y(I+1,J) - 2.D0

Y(I,J) + Y(I-1,J)) / (H1

END DO

C(1) = 1.D0

D(1) = 0.D0

E(1) = 0.D0

F(1) = 0.D0

B(2) = 0.D0

C(2) = 5.D0

ALPHA / (H1

4) + 0.5D0 / TAU

C(N1-1) = 5.D0

ALPHA / (H1

4) + 0.5D0 / TAU

D(N1-1) = 0.D0

A(N1) = 0.D0

B(N1) = 0.D0

C(N1) = 1.D0

F(N1) = 0.D0

C SOLUTION OF THE DIFFERENCE PROBLEM

ITASK2 = 1

CALL PROG5 ( N1, A, B, C, D, E, F, YY, ITASK2 )

DOI=1,N1

Y1(I,J) = YY(I)

END DO

C SWEEP OVER X2

DOI=2,N1-1

C COEFFICIENTS IN THE DIFFERENCE SCHEME FOR THE INVERSE PROBLEM

DOJ=2,N2-1

A(J) = ALPHA / (H2

B(J) = 4.D0

ALPHA / (H2

C(J) = 6.D0

ALPHA / (H2

4) + 0.5D0 / TAU

D(J) = 4.D0

ALPHA / (H2

E(J) = ALPHA / (H2

F(J) = 0.5D0

Y(I,J) / TAU

+ - (Y(I,J+1) - 2.D0

Y(I,J) + Y(I,J-1)) / (H2

END DO

C(1) = 1.D0

D(1) = 0.D0

E(1) = 0.D0

F(1) = 0.D0

B(2) = 0.D0

C(2) = 5.D0

ALPHA / (H2

4) + 0.5D0 / TAU