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

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

156 Chapter 5 Solution methods for ill-posed problems

Exercise 5.11 Suppose that, in the solution of problem (5.1), (5.2), the total number

of iterations in the method (5.31)–(5.33) is chosen from the condition

Au

n(δ)

− f

≤δ<Au

n(δ)−1

− f



and, in addition, Au

− f

 >δ. Show that, in this case, u

n(δ)

− u→0asδ → 0.

Exercise 5.12 Prove the ill-posedness of the problem in which it is required to solve

the ﬁrst-kind Fredholm integral equation



K (x, s)μ(s) ds = ϕ(x), c ≤ x ≤ d

with a kernel K (x, s) continuously differentiable with respect to both variables in the

case of μ(x) ∈ C[a, b], ϕ(x) ∈ C[c, d].

Exercise 5.13 Consider the continuation problem for a gravitational ﬁeld disturbed

with local anomalies with non-negative anomalous (excessive with respect to sur-

rounding rocks) density. We seek the solution of the continuation problem (5.53),

(5.55) in the form (see (5.57))



∞

−∞

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

assuming, in addition, that the anomalies are localized at depths greater than H . Show

that μ(x ) ≥ 0, −∞ < x < ∞.

Exercise 5.14 In the program PROBLEM4, implement the conjugate gradient method.

Examine the efﬁciency of the method as applied to problems with small inaccuracies

in the gravitational ﬁeld at the earth surface (δ ≤ 0.01).

Exercise 5.15 Extend the program PROBLEM4 to the case in which the approximate

solution of the continuation problem is solved in the form of a double-bed potential.

Based on the solution of the model problem, perform an analysis of possibilities of-

fered by this approach.

6 Right-hand side identiﬁcation

An important class of inverse mathematical physics problem is the determination of

unknown right-hand sides of equations. In such problems, additional information

about the solution is provided either throughout the calculation domain or over some

part of the domain, in particular, on the domain boundary. Another class consists in

calculating a differential operator for an approximately given right-hand side and solu-

tion which is approximately known throughout the whole domain. Also worth noting

are speciﬁc features of right-hand side determination algorithms for non-stationary

problems. Here the values of the right-hand side are successively determined at ﬁxed

times. In the non-stationary case, additional information about the solution is required

on parts of the calculation domain. At the end of this chapter, computational algo-

rithms for solving such inverse boundary value problems involving stationary or non-

stationary mathematical physics equations are considered. Computational data are

provided that are obtained in several numerical experiments on the approximate solu-

tion of inverse problems.

6.1 Reconstruction of the right-hand side from known solu-

tion in the case of stationary problems

In this section, we consider a problem in which it is required to reconstruct the right-

hand side of a second-order ordinary differential equation whose solution is known

accurate to some given inaccuracy. A simplest computational algorithm uses stan-

dard difference approximations, the mesh size being the regularization parameter. The

possibility of using alternative approaches is noted.

6.1.1 Problem statement

As a simplest mathematical model, we use the second-order ordinary differential equa-

tion

−



k(x)



+ q(x)u = f (x), 0 < x < l (6.1)

under the following usual constraints imposed on the coefﬁcients:

k(x) ≥ κ>0, q(x) ≥ 0.

In the direct problem, it is required to determine the function u(x) from equation

(6.1) with given coefﬁcients k(x), q(x), right-hand side f (x), and additional condi-

tions given on the boundary. In a simplest case, ﬁrst-kind homogeneous conditions are

given:

u(0) = 0, u(l) = 0. (6.2)

158 Chapter 6 Right-hand side identiﬁcation

Formulate now the inverse problem. We assume that the right-hand side f (x) is

unknown. To ﬁnd it, we have to set some additional information about the solution.

Since, here, a function of x is being sought, it is therefore desirable that the additional

information be also given as a function of x. This allows us to assume that known is

the solution u(x), x ∈ [0, l].

With known exact u(x), x ∈ [0, l], the right-hand side can be uniquely found, and

it follows from (6.1) that

f (x) =−



k(x)



+ q(x)u, 0 < x < l. (6.3)

The use of formula (6.3) implies, for instance, that f (x) ∈ C(0, l) if u(x) ∈ C

[0, l]

in the case of k(x) ∈ C

[0, l] and q(x) ∈ C[0, l].

The point here is that the input data (in the case of interest, the solution u(x), x ∈

[0, l]) are given approximately and do not belong to the mentioned smoothness class.

A typical situation is such that, for instance, instead of the exact solution u(x), x ∈

[0, l] we know the function u

(x), x ∈ [0, l] and, in addition, u

(x) ∈ C[0, l] and

u

(x) − u(x)

C[0,l]

≤ δ, (6.4)

where

v(x)

C[0,l]

= max

x∈[0,l]

|v(x)|.

Problem (6.3), (6.4) is a problem in which it is required to calculate the values

of a differential operator. This problem belongs to the class of classically ill-posed

problems and, to be solved, needs some regularizing algorithm to be applied. Consider

some basic possibilities along this line.

6.1.2 Difference algorithms

Most naturally, difference methods can be used in calculating the right-hand side f (x),

0 < x < l. Over the interval

 = [0, l], we introduce a uniform grid 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 denote the approximate solution of problem (6.3), (6.4) at internal nodes as

(x). Using the standard notation adopted in the theory of difference schemes, we set

(x) =−(ay

¯x

)

+ cy, x ∈ ω, (6.5)

where for the mesh function y(x) we have:

y(x) = u

(x), x ∈¯ω.

Section 6.1 Right-hand side reconstruction from known solution: stationary problems 159

For the difference coefﬁcients, we set:

a(x) = k(x − 0.5h), x, x + h ∈ ω, c(x) = q(x), x ∈ ω.

The difference algorithm (6.5) is a regularizing one provided that a rule is given

by which the discretization grid size h = h(δ) can be chosen as a function of the

inaccuracy δ such that the norm of the inaccuracy z = f

− f vanishes as δ → 0.

We would also like to know the rate of convergence of the approximate solution to the

exact solution.

To optimize the discretization grid size and minimize the inaccuracy, we identify

two classes of exact solutions. We assume that

u(x)

[0,l]

≤ M, (6.6)

or, alternatively,

u(x)

[0,l]

≤ M, (6.7)

where

v(x)

[0,l]

= max

x∈[0,l]

max

0≤k≤n



(x)



In the case of (6.6) or (6.7), we impose on the coefﬁcients the constraints

k(x) ∈ C

[0, l], q(x) ∈ C

[0, l]

k(x) ∈ C

[0, l], q(x) ∈ C

[0, l].

Consider the problem for the inaccuracy z(x) = f

(x) − f (x), x ∈ ω. From (6.3),

(6.5), we easily obtain:

z(x) = z

(1)

(x) + z

(2)

(x), x ∈ ω, (6.8)

where

(1)

(x) =−(au

¯x

)

+ cu − f (x), x ∈ ω, (6.9)

(2)

(x) = f

(x) + (au

¯x

)

− cu , x ∈ ω. (6.10)

Let us estimate now the individual terms in (6.8) as dependent on the a priori as-

sumptions about the exact solution of problem (6.3). The ﬁrst term (z

(1)

(x)) is the

inaccuracy in approximating the differential operator with the difference operator. The

second term ( z

(2)

(x)) reﬂects the effect induced by the input-data inaccuracy.

For the approximation inaccuracy, we readily obtain:

(1)

(x) =−(au

¯x

)

+ cu +



k(x)



− q(x)u

=−

i+1

− a

) −

i+1

+ a

) −

i+1

− a

)

+cu + O(h

) +

+ k(x)

− q(x)u.

160 Chapter 6 Right-hand side identiﬁcation

Here β = 1 under the assumption of (6.7) and β = 2, provided that the inequality

(6.6) holds.

From similar representations, for the inaccuracy we obtain the estimate

z

(1)

(x)≤M

, (6.11)

with the used settings

z(x)=z(x)

C(ω)

= max

x∈ω

|z(x)|.

The constant M

in (6.11) depends not only on the solution but also on the coefﬁ-

cient k(x) and on the derivatives of k(x), i. e., M

= M

(u, k).

For the second inaccuracy component, from (6.5), (6.10) we obtain

z

(2)

(x)=−(a(y − u)

¯x

)

+ c(y − u)

≤



max



a(x), a(x + h)



+c(x)



max

x∈¯ω

y(x) − u(x).

With (6.4), we obtain

z

(2)

(x)≤M

δ. (6.12)

The constant

= max



a(x), a(x + h)



+c(x)

depends on the discretization grid size, i. e., M

= M

(h), and increases with decreas-

ing h.

With (6.11), (6.12), from (6.8)–(6.10) we obtain the desired estimate for the inaccu-

racy:

z(x)≤M

+ M

(h)δ. (6.13)

It follows from here that for the convergence with δ → 0 it is required that δh

−2

→ 0.

In the classes of a priori assumptions (6.6) and (6.7) about the solution, from the es-

timate (6.13) we obtain the optimal discretization grid size and the rate of convergence

opt

= O(δ

1/4

), z(x)=O(δ

1/2

)

in the case of (6.6) (β = 2) and

opt

= O(δ

1/3

), z(x)=O(δ

1/3

)

in the case of (6.7) (β = 1).

Thus, we see that the discretization grid size h in the difference calculation of the

right-hand side can be used as regularization parameter. Its value must be matched

with the input-data inaccuracy δ. The optimal value of h depends on the constants M

and M

, i. e., h

opt

= h

opt

, M

). Most frequently, no direct calculation of these con-

stants is possible; in such cases, the discretization grid size can be chosen considering

the discrepancy criterion.

Section 6.1 Right-hand side reconstruction from known solution: stationary problems 161

In numerical solution of unstable problems, discretization (passage to a ﬁnite-

difference problem) always brings about a regularization effect. This note applies both

to difference and projection methods. Here, as the discretization parameter, the dimen-

sion of the problem (discretization step) can be used. We can say that the numerical

methods possess the self-regularization property. Yet, if we try to raise the dimension

of the problem in order to reﬁne the solution, starting from a certain moment we ob-

tain progressively worsening results: here, ill-posedness of the problem is manifested.

That is why we have to terminate the calculations in due time and restrict ourselves to

an approximate solution obtained with some optimal discretization.

Discretization offers rather restricted opportunities in regularization of problems. In

computational practice, they often take another strategy. Methods are used in which the

continuous problem possesses regularizing properties. In the course of discretization,

the closeness of the discrete problem to the regularized differential problem due to the

use, for instance, of a sufﬁciently ﬁne grid cannot be explicitly traced. In this case,

additional regularizing properties resulting from discretization are ignored as exerting

only an insigniﬁcant inﬂuence. More accurate account of discretization inaccuracies is

also possible in the cases in which the effect due to the inaccuracies cannot be ignored.

In particular, the value of regularization parameter is matched not only with the right-

hand side inaccuracy, but also with the inaccuracy of the operator in the ﬁrst-kind

equation (the so-called general discrepancy principle).

6.1.3 Tikhonov regularization

Consider problem (6.2)–(6.4) on the operator level. On the set of functions vanishing

at the end points of the segment [0, l] (see (6.2)) we deﬁne the operator D as

Du =−



k(x)



+ q(x)u.

Then, with exact input data the problem (6.2), (6.3) can be written as

f = D u. (6.14)

In H = L

(0, l),wehave:

D = D

∗

≥ mE, m = κ

Hence, the operator D

−1

exists and

0 < D

−1

≤

The latter allows us to pass from equation (6.14) to the equation

A f = u, (6.15)

162 Chapter 6 Right-hand side identiﬁcation

where A = A

∗

= D

−1

Instead of the exact right-hand side of (6.15), the function u

is known to some

accuracy given by

u

− u≤δ. (6.16)

The problem (6.15), (6.16) was considered previously. This problem can be stably

solved by the Tikhonov regularization method. In this case, the approximate solu-

tion f

, α = α(δ) is to be found from

( f

) = min

v∈H

(v), (6.17)

where

(v) =Av − u



+ αv

. (6.18)

Using the adopted notation, we can write the functional (6.18) as

(v) =D

−1

v − u



+ αv

This functional is inconvenient for use because the values of D can be calculated easier

than the inversion of the operator. Instead of (6.18), we can use the more general

functional

(v) =D

−1

v − u



+ αv

with G = G

∗

> 0. We set G = D

∗

D; this yields

(v) =v − Du



+ αDv

. (6.19)

In this way, we arrive at a version of the Tikhonov regularization method for the

approximate solution of problem (6.14), (6.16) based on the solution of the variational

problem (6.17), (6.19). The situation around this algorithm was clariﬁed previously by

considering the intermediate problem (6.15), (6.16).

Method (6.17), (6.19) as applied to the approximate solution of problem (6.14),

(6.16) can be substantiated using the same scheme as in the approximate solution of

problem (6.15) (6.16) by method (6.17), (6.18). This makes further consideration of

this point unnecessary.

The Euler equation for (6.17), (6.19) has the form

(E + αD

∗

D) f

= Du

. (6.20)

In the case of (6.2)–(6.4), equation (6.20) yields a boundary value problem for the

fourth-order ordinary differential equation. For the solution of (6.20) the following

standard a priori estimate is valid:

 f

≤

√

u

. (6.21)

This estimate shows that the approximate solution is stable with respect to weak right-

hand side perturbations.

Section 6.1 Right-hand side reconstruction from known solution: stationary problems 163

6.1.4 Other algorithms

We present some other possibilities for the approximate solution of problem (6.14),

(6.16). First of all, considering the case in which one obtains the approximate solution

from equation (6.20), we can give another interpretation to the Tikhonov regularization

method.

At the ﬁrst stage, we deﬁne an auxiliary function

f = D u

. (6.22)

Thus,

f is the solution of the problem with noisy right-hand side. At the second stage,

the found solution is to be processed by applying a smoothing treatment. This can

most naturally be done by minimizing the functional

(v) =v −

f 

+v

. (6.23)

The minimum condition for the functional gives

(E + αG) f

f . (6.24)

In this interpretation, method (6.17), (6.19) refers to the case with G = D

∗

D chosen

in (6.17), (6.22), (6.23).

With

G ≥ D

∗

for the solution of (6.22), (6.24) there holds the estimate (6.21). Note one such possi-

bility that is of interest from the computational point of view. In the case of G = D

have G = D

∗

D > 0, and the method requires inversion of the operator E +αG. In the

model problem under consideration it is required to solve a boundary value problem

for a fourth-order equation. With

G = D

+ 2

√

equation (6.24) takes the form

(E +

√

α D)

f . (6.25)

Hence, we can restrict ourselves to double inversion of the operator E +

√

α D (i. e.,

to solving two boundary value problems for a second-order equation). Note that algo-

rithm (6.22), (6.25) can be considered as double smoothing performed with the opera-

tor G = D.

Iteration methods deserve a special detailed consideration. Here, however, just one

brief remark will be made. The approximate solution f

n(δ)

of problem (6.15), (6.16)

can be found using the simple-iteration method

k+1

− f

+ A f

= u

, k = 0, 1,...,

164 Chapter 6 Right-hand side identiﬁcation

which was considered at length previously. Turning back to problem (6.14), (6.16),

rewrite the iterative process as

k+1

− f

+ f

= Du

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

Thus, for the solution of the problem at the next iteration step to be obtained, one has

to invert the operator D. The regularization effect can be achieved with a properly

chosen reconditioner B = D.

6.1.5 Computational and program realization

We construct a computational algorithm for problem (6.14), (6.16) around the

Tikhonov regularization method (6.17), (6.19). First, pass to a discrete problem, as-

suming the discretization step to be sufﬁciently small.

On the set of mesh functions deﬁned on the grid ¯ω and vanishing on ∂ω, we deﬁne

the difference operator

Dy =−(ay

¯x

)

+ cy, x ∈ ω.

On discretization, we pass from problem (6.14) to the problem

ϕ = Dy. (6.27)

Instead of the mesh function y(x), x ∈ ω, the function y

(x), x ∈ ω is now given.

The inaccuracy level is deﬁned by the parameter δ:

y

− y≤δ, (6.28)

where ·is the norm in H = L

(ω):

y=(y, y)

1/2

,(y,v) =



x∈ω

y(x)v(x)h.

In H ,wehave:

D = D

∗

≥ κλ

E,λ

sin

πh

≥

In the Tikhonov method as applied to problem (6.27), (6.28) there arises (see (6.20))

the equation

(E + α D

)ϕ

= Dy

(6.29)

for the approximate solution ϕ

(x), x ∈¯ω. The value of α = α(δ) can be found using

the discrepancy principle:

y

− D

−1

=δ. (6.30)

Section 6.1 Right-hand side reconstruction from known solution: stationary problems 165

Note some speciﬁc features in the computational realization of the method.

Equation (6.29) is a difference equation with ﬁve-diagonal matrix. This equation

can be solved by the sweep method (Thomas algorithm). Below, we give the sweep-

method calculation formulas for the following ﬁve-diagonal system:

− D

+ E

= F

, (6.31)

−B

+ C

− D

+ E

= F

, (6.32)

i−2

− B

i−1

+ C

− D

i+1

+ E

i+2

= F

, i = 2,...,N − 2,(6.33)

N −1

N −3

− B

N −1

N −2

+ C

N −1

− D

N −1

= F

N −1

, (6.34)

N −2

− B

N −1

+ C

= F

. (6.35)

The sweep method for system (6.31)–(6.35) uses the following solution representa-

tion:

= α

i+1

− β

i+1

i+2

+ γ

i+1

, i = 0, 1,...,N − 2, (6.36)

N −1

= α

+ γ

. (6.37)

Like in the ordinary three-point sweep method, we ﬁnd ﬁrst the sweep coefﬁcients.

Using (6.36), we express y

i−1

and y

i−2

in terms of y

and y

i+1

i−1

= α

− β

i+1

+ γ

, i = 1, 2,...,N − 1, (6.38)

i−2

= (α

i−1

− β

i−1

− β

i−1

i+1

+ α

i−1

+ γ

i−1

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

(6.39)

Substitution of (6.38) and (6.39) into (6.33) yields

− A

i−1

+ α

i−1

− B

))y

= (D

+ β

i−1

− B

))y

i+1

− E

i+2

+ F

− A

i−1

− γ

i−1

− B

), i = 2, 3,...,N − 2.

In view of (6.36), we obtain

i+1

= S

−1

+ β

i−1

− B

)), β

i+1

= S

−1

i+1

= S

−1

− A

i−1

− γ

i−1

− B

)),

(6.40)

where

= (C

− A

i−1

+ α

i−1

− B

)), i = 2, 3,...,N − 2. (6.41)

The recurrent relations (6.40), (6.41) can be used provided that the coefﬁcients

,β

and γ

(i = 1, 2) are known. It follows from (6.31) and (6.36) that

= C

−1

,β

= C

−1

,γ

= C

−1

. (6.42)