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

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

146 Chapter 5 Solution methods for ill-posed problems

Figure 5.2 Solution of the direct problem at various depths

we consider a model problem with two anomalies, circular in cross-section, with cen-

ters at (x

, z

) = (0.8, −0.3) and (x

, z

) = (1.1, −0.4). The exact solution of prob-

lem (5.53)–(5.55) is

u(x, z) = c

z − z

(x − x

)

+ (z − z

)

+ c

z − z

(x − x

)

+ (z − z

)

Suppose that c

= 0.3 and c

= 1.2, i. e., the deeper anomaly has fourfold greater

power.

The exact solution of the problem at various depths H is shown in Figure 5.2. From

the measured data, at H = 0 and H = 0.1 two sources of gravitational perturba-

tions can be suspected. At larger depths (as we approach the sources), two individual

anomalies can be identiﬁed (at H = 0.2 and especially at H = 0.25). It is this cir-

cumstance that necessitates the continuation of potential ﬁelds towards anomalies.

5.5.2 Integral equation

To ﬁnd an approximate solution of the continuation problem (5.53)–(5.55), one can

use many computational algorithms. We will dwell here on the use of the integral

equation method. We approximate the gravitational potential produced by anomalies

with a potential generated by an elemental bed with a carrier located at the segment

0 ≤ x ≤ L and at the depth z =−H . Accurate to a multiplier, we have

U (x, z) =



ln r μ(s) ds, r =



(x − s)

+ (z + H)

, z > −H.

Section 5.5 Program implementation and computational experiments 147

For the derivative with respect to the vertical variable we obtain:

u(x, z) =



z + H

(x − s)

+ (z + H)

μ(s) ds, z > −H. (5.56)

For the unknown density, from (5.54) and (5.56) we obtain the ﬁrst-kind integral equa-

tion



K (x, s)μ(s) ds = ϕ(x) (5.57)

with the symmetric kernel

K (x, s) =

(x − s)

+ H

We write equation (5.57) as the ﬁrst-kind equation

Aμ = ϕ (5.58)

in which

Aμ =



K (x, s)μ(s) ds.

The integral operator A : H → H , where H = L

(0, L), is a self-adjoint operator

(Aμ, ν) = (μ, Aν). Besides, this operator is a bounded operator because

(Aμ, μ) =



K (x, s)μ(s)μ(x) ds dx

≤





(x, s) ds dx



1/2





(s)μ

(x) ds dx



1/2

≤ γ

μ

where γ

< L/H.

The integral operator under consideration is a self-adjoint, bounded operator. Yet,

the iteration method cannot be directly applied to (5.58) because, generally speaking,

the operator A here is not a positive (non-negative) operator. Hence, symmetrization

of (5.58) is to be applied:

∗

Aμ = A

∗

ϕ, (5.59)

which we can ﬁnally use in the iteration method.

5.5.3 Computational realization

In the numerical solution of the integral equation (5.59), we use a uniform grid

¯ω ={x | x = (i − 0.5)h, i = 1, 2,...,N, Nh = L}.

We designate the approximate solution as y

= y(x

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

148 Chapter 5 Solution methods for ill-posed problems

We put into correspondence to the integral operator A the difference operator

y)(x

) =



j=1

K (x

, s

h. (5.60)

This approximation corresponds to the use of the rectangle quadrature formula. The

difference analogue of (5.59) is the equation

Ry = f, R = A

∗

, f = A

∗

ϕ. (5.61)

To approximately solve equation (5.61), we use the iteration method

k+1

− y

k+1

+ Ry

= f, k = 0, 1,... . (5.62)

Two choices of iteration parameters were considered:

1. Simple iteration method, in which τ

= τ = const;

2. Steepest descend method, in which

k+1

, r

)

(Rr

, r

)

, r

= Ry

− f.

To model the input-data inaccuracy, we perturbed the exact solution at the grid nodes

by the law

(x) = ϕ(x) + 2ζ(σ(x) − 1/2), x = x

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

where σ(x) is a random quantity normally distributed over the interval from 0 to 1.

The parameter ζ deﬁnes the inaccuracy level in the right-hand side of (5.58), and

ϕ

(x)) − ϕ(x)≤δ = ζ

√

The iterative process is to be terminated considering the discrepancy. Some other

speciﬁc features of the computational algorithm will be noted below.

5.5.4 Program

To approximately solve the model two-dimensional continuation problem for gravita-

tional ﬁelds, we used the program presented below.

Program PROBLEM4

C PROBLEM4 - CONTINUATION OF ANAMALOUS GRAVITATIONAL FIELD

PARAMETER ( HH = 0.225, DELTA = 0.1, ISCHEME = 1, N = 100 )

Section 5.5 Program implementation and computational experiments 149

DIMENSION Y(N), X(N), PHI(N), PHID(N), F(N), A(N,N)

+ ,V(N), R(N), AR(N), U(N), UY(N), UY1(N)

C PARAMETERS:

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

C HH - CARRIER DEPTH;

C H0 - CONTINUATION DEPTH;

C ISCHEME - SIMPLE-ITERATION METHOD (ISCHEME = 0),

C QUICKEST DESCEND METHOD (ISCHEME = 1);

C N - NUMBER OF GRID NODES;

C Y(N) - SOLUTION OF THE INTEGRAL EQUATION;

C U(N) - EXACT ANOMALOUS FIEKD AT THE DEPTHZ=H0;

C UY(N) - CALCULATED FIELD AT THE DEPTHZ=H0;

C UY1(N) - CALCULATED FIELD AT THE EARTH SURFACE;

C PHI(N) - EXACT ANOMALOUS FIELD AT THE EARTH SURFACE;

C PHID(N) - DISTURBED FIELD AT THE EARTH SURFACE;

C DELTA - INPUT-DATA INACCURACY LEVEL;

C A(N,N) - PROBLEM MATRIX;

XL=0.

XR=2.

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

C GRID

H = (XR - XL)/N

DOI=1,N

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

END DO

C DATA AT THE EARTH SURFACE

C1 = 0.3

X1 = 0.8

Z1 = - 0.3

C2 = 1.2

X2 = 1.1

Z2 = - 0.4

DO I = 1,N

PHI(I) = - C1/((X(I)-X1)

2+Z1

- C2/((X(I)-X2)

2+Z2

PHID(I) = PHI(I) + 2.

DELTA

(RAND(0)-0.5)/SQRT(XR-XL)

END DO

C EXACT IN-EARTH DATA

H0 = 0.2

DO I = 1,N

U(I) = - C1/((X(I)-X1)

2 + (H0+Z1)

(H0+Z1)

- C2/((X(I)-X2)

2 + (H0+Z2)

(H0+Z2)

END DO

C RIGHT-HAND SIDE

DO I = 1,N

SUM=0.

DO J = 1,N

SUM = SUM + AK(X(I),X(J),HH)

PHID(J)

150 Chapter 5 Solution methods for ill-posed problems

END DO

F(I) = SUM

END DO

C MATRIX ELEMENTS

DO I = 1,N

DO J = 1,N

SUM=0.

DO K = 1,N

SUM = SUM + AK(X(I),X(K),HH)

AK(X(K),X(J),HH)

END DO

A(I,J) = SUM

END DO

C ITERATIVE PROCESS

IT=0

ITMAX = 100

DO I = 1,N

Y(I) = 0.

END DO

100 IT = IT+1

C CACULATIONV=RY

DO I = 1,N

SUM=0.

DO J = 1,N

SUM = SUM + A(I,J)

Y(J)

END DO

V(I) = SUM

END DO

C ITERATION PARAMETER

DO I = 1,N

R(I) = V(I) - F(I)

END DO

IF ( ISCHEME.EQ.0 ) THEN

C SIMPLE-ITERATION METHOD

TAU = 0.1

END IF

IF ( ISCHEME.EQ.1 ) THEN

C QUICKEST DESCEND METHOD

DO I = 1,N

R(I) = V(I) - F(I)

END DO

DO I = 1,N

SUM=0.

DO J = 1,N

Section 5.5 Program implementation and computational experiments 151

SUM = SUM + A(I,J)

R(J)

END DO

AR(I) = SUM

END DO

SUM1 = 0.

SUM2 = 0.

DO I = 2,N-1

SUM1 = SUM1 + R(I)

R(I)

SUM2 = SUM2 + AR(I)

R(I)

END DO

TAU = SUM1/SUM2

END IF

C NEXT APPROXIMATION

DO I = 1,N

Y(I) = Y(I) - TAU

R(I)

END DO

C ITERATION LOOP EXIT

SUM0 = 0.

DO I = 1,N

SUM=0.

DO J = 1,N

SUM = SUM + AK(X(I),X(J),HH)

Y(J)

END DO

SUM0 = SUM0 + (SUM

H - PHID(I))

END DO

SL2 = SQRT(SUM0)

WRITE ( 01,

) IT, TAU, SL2

IF (SL2.GT.DELTA .AND. IT.LT.ITMAX) GO TO 100

CLOSE ( 01 )

C CALCULATION OF THE FOUND FIELD AT THE DEPTHZ=H0

DO I = 1,N

SUM=0.

DO J = 1,N

SUM = SUM + AK(X(I),X(J),HH-H0)

Y(J)

END DO

UY(I) = SUM

END DO

C CALCULATION OF THE FOUND FIELD AT THE EARTH SURFACE

DO I = 1,N

SUM=0.

DO J = 1,N

SUM = SUM + AK(X(I),X(J),HH)

Y(J)

END DO

UY1(I) = SUM

END DO

STOP

END

FUNCTION AK(X,S,HH)

152 Chapter 5 Solution methods for ill-posed problems

C INTEGRAL-EQUATION KERNEL

AK = HH/((X-S)

2+HH

RETURN

END

5.5.5 Computational experiments

Consider several examples for the solution of the continuation problem in the case of

H = 0.225. We use a computation grid whose total number of nodes is N = 100.

The most interesting dependence here is the accuracy in reconstructing the anomalous

gravitational ﬁeld versus the input-data inaccuracy δ. Figure 5.3 shows data obtained

by solving the problem with the input-data inaccuracy δ = 0.1 (in setting the ﬁeld at

the earth surface z = 0). Presented are the exact and found ﬁelds at the depth z =−2.

Figure 5.3 Solution of the continuation problem obtained with δ = 0.1

Similar data obtained at a larger inaccuracy level (δ = 0.2) are shown in Figure 5.4.

At the latter inaccuracy level, two individual anomalies are hard to identify (here, no

two-peak conﬁguration is observed). Even a decrease of δ down to δ = 0.01 little

saves the situation (see Figure 5.5).

Section 5.5 Program implementation and computational experiments 153

Figure 5.4 Solution of the continuation problem obtained with δ = 0.2

Figure 5.5 Solution of the continuation problem obtained with δ = 0.01

154 Chapter 5 Solution methods for ill-posed problems

We now provide a few statements concerning the use of these or those iteration

methods. In the program PROBLEM4, the possibility to use either the simple-iteration

method (constant iteration parameter) or the steepest descend method is provided.

A slight modiﬁcation allows one to implement the iterative conjugate gradient method,

to be used in the case of problems with small right-hand side inaccuracies, in which

the steepest descend method fails to provide a desired rate of convergence.

The number of iterations necessary for solving the problem with δ = 0.1bythe

steepest descend method is n = 6. The total number of iterations in the simple-

iteration method versus the iteration parameter is illustrated by the data in Table 5.1.

τ 0.05 0.1 0.15 0.2 0.25 0.3

n 43 21 14 11 9 32

Table 5.1 Total number of iterations in the simple-iteration method

In the iterative solution of (5.59), the optimal iteration parameter is τ = τ

= 2/ ¯γ

where

∗

A ≤¯γ

can be estimated invoking the estimate ¯γ

< L

5.6 Exercises

Exercise 5.1 Show that in the Tikhonov regularization method for the solution there

holds the a priori estimate

u

≤

√

 f



that expresses stability with respect to the right-hand side.

Exercise 5.2 Formulate conditions for convergence (analogue to Theorem 5.1, 5.2) of

the approximate solution in H

, D = D

∗

> 0 found by minimization of the functional

(v) =Av − f



+ αv

Exercise 5.3 Construct an example illustrating the necessity of condition δ

/α(δ) →

0 with δ → 0 (Theorem 5.1) for convergence of the approximate solutions to the exact

solution in the Tikhonov method.

Exercise 5.4 Suppose that the exact solution of problem (5.1) is

A

−2

u≤M,

where M = const > 0, and for the right-hand side inaccuracy the estimate (5.2) holds.

Then, for the approximate solution u

found as the solution of problem (5.4), (5.5)

there holds the a priori estimate

z≤

√

+ αM.

Section 5.6 Exercises 155

Exercise 5.5 Obtain an estimate for the rate of convergence in the Tikhonov method

under the a priori constraints

A

−p

u≤M, 0 < p < 2

posed on the exact solution of problem (5.1), (5.2).

Exercise 5.6 Let u

be the solution of problem (5.4), (5.5) and

m(α) =Au

− f



+ αu



ϕ(α) =Au

− f



, ψ(α) =u



Show that in the case of f

= 0 and 0 <α

<α

there hold the inequalities

m(α

)<m(α

), ϕ(α

)<ϕ(α

), ψ(α

)>ψ(α

Exercise 5.7 Show that for any δ ∈ (0,  f

) there exists a unique solution α = α(δ)

of the equation

ϕ(α) =Au

− f



= δ

such that α(δ) → 0 and u

α(δ)

− u→0asδ → 0 (substantiation of the choice of

regularization parameter based on the discrepancy).

Exercise 5.8 In solving problem (5.1), (5.2) with A = A

∗

> 0, the algorithm of

simpliﬁed regularization uses the equation

+ αu

= f

(5.63)

for ﬁnding the approximate solution. Formulate the related variational problem.

Exercise 5.9 Prove that in the case of

u

−1

≤ M

for the inaccuracy in the approximate solution of problem (5.1), (5.2) found from

(5.63) there holds the priori estimate

z

≤

Exercise 5.10 Suppose that for the approximate solution of problem (5.1), (5.2) one

uses the iteration method

k+1

− u

+ Au

= f

, k = 0, 1,...

with B = B

∗

> 0 and

B ≤ A ≤ γ

B,γ

> 0, 0 <τ <

Show that in the iteration method (5.31)–(5.33) u

n(δ)

− u

→ 0asδ → 0 provided

that n(δ) →∞and n(δ)δ → 0asδ → 0.