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

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

336 Chapter 7 Evolutionary inverse problems

C SOLUTION AT THE OBSERVATION POINTS

CALL PU (IM, JM, A(10

N+1), N1, N2, FF, L)

DO IP = 1,L

FIK(IP,K) = FF(IP)

END DO

C CONJUGATE STATE

T = TMAX + TAU

C INITIAL CONDITION

DOI=1,N

A(10

N+I) = 0.D0

END DO

DOK=2,M+1

C DIFFERENCE-SCHEME COEFFICIENTS

CALL FDS (A(1), A(N+1), A(2

N+1), A(9

N+1), A(10

N+1),

+ H1, H2, N1, N2, TAU)

CALL RHS (A(9

N+1), N1, N2, H1, H2, FID, FIK, IM, JM, L, M, K)

C SOLUTION OF THE DIFFERENCE PROBLEM

IDEFAULT(1) = 0

IREPT = 0

CALL SBAND5 (N, N1, A(1), A(10

N+1), A(9

N+1), EPSR, EPSA)

END DO

C SIMPLE ITERATION METHOD

C NEXT APPROXIMATION

CALL BS (A(10

N+1), A(1), A(N+1), A(2

N+1), A(12

N+1), N1, N2 )

DOI=1,N

A(11

N+I) = A(11

N+I) - SS

A(12

N+I)

END DO

C EXIT FROM THE ITERATIVE PROCESS BY THE DISCREPANCY CRITERION

SUM = 0.D0

DOK=1,M

DOIP=1,L

SUM = SUM + (FID(IP,K) - FIK(IP,K))

TAU

END DO

SUM = SUM/(L

TMAX)

SL2 = DSQRT(SUM)

WRITE (

) IT, SL2

IF ( SL2.GT.DELTA ) GO TO 100

C SOLUTION

WRITE ( 01,

) (A(11

N+I), I=1,N)

WRITE ( 01,

) ((FI(IP,K), IP=1,L), K=1,M)

CLOSE (01)

STOP

Section 7.5 Continuation of non-stationary ﬁelds from point observation data 337

END

SUBROUTINE INIT (U, X1, X2, N1, N2)

C INITIAL CONDITION

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION U(N1,N2), X1(N1), X2(N2)

DOI=1,N1

DOJ=1,N2

U(I,J) = 0.D0

IF ((X1(I)-0.6D0)

2 + (X2(J)-0.6D0)

2.LE.0.04D0)

+ U(I,J) = 1.D0

END DO

RETURN

END

SUBROUTINE FDS (A0, A1, A2, F, U, H1, H2, N1, N2, TAU)

C GENERATION OF DIFFERENCE-SCHEME COEFFICIENTS

C FOR PARABOLIC EQUATION WITH CONSTANT COEFFICIENTS

C IN THE CASE OF PURELY IMPLICIT SCHEME

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION A0(N1,N2), A1(N1,N2), A2(N1,N2), F(N1,N2), U(N1,N2)

DOJ=2,N2-1

DOI=2,N1-1

A1(I-1,J) = 1.D0/(H1

H1)

A1(I,J) = 1.D0/(H1

H1)

A2(I,J-1) = 1.D0/(H2

H2)

A2(I,J) = 1.D0/(H2

H2)

A0(I,J) = A1(I,J) + A1(I-1,J) + A2(I,J) + A2(I,J-1)

+ + 1.D0/TAU

F(I,J) = U(I,J)/TAU

END DO

C FIRST-KIND HOMOGENEOUS BOUNDARY CONDITION

DOJ=2,N2-1

A0(1,J) = 1.D0

A1(1,J) = 0.D0

A2(1,J) = 0.D0

F(1,J) = 0.D0

END DO

DOJ=2,N2-1

A0(N1,J) = 1.D0

A1(N1-1,J) = 0.D0

A1(N1,J) = 0.D0

A2(N1,J) = 0.D0

F(N1,J) = 0.D0

END DO

DOI=2,N1-1

A0(I,1) = 1.D0

A1(I,1) = 0.D0

338 Chapter 7 Evolutionary inverse problems

A2(I,1) = 0.D0

F(I,1) = 0.D0

END DO

DOI=2,N1-1

A0(I,N2) = 1.D0

A1(I,N2) = 0.D0

A2(I,N2) = 0.D0

A2(I,N2-1) = 0.D0

F(I,N2) = 0.D0

END DO

A0(1,1) = 1.D0

A1(1,1) = 0.D0

A2(1,1) = 0.D0

F(1,1) = 0.D0

A0(N1,1) = 1.D0

A2(N1,1) = 0.D0

F(N1,1) = 0.D0

A0(1,N2) = 1.D0

A1(1,N2) = 0.D0

F(1,N2) = 0.D0

A0(N1,N2) = 1.D0

F(N1,N2) = 0.D0

RETURN

END

SUBROUTINE RHS (F, N1, N2, H1, H2, FI, FIK, IM, JM, L, M, K)

C RIGHT-HAND SIDE IN THE EQUATION FOR THE CONJUGATE STATE

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION F(N1,N2), FI(L,M), FIK(L,M), IM(L), JM(L)

DOJ=2,N2-1

DOI=2,N1-1

DOIP=1,L

IF (I.EQ.IM(IP).AND.J.EQ.JM(IP)) THEN

F(I,J) = F(I,J) + (FIK(IP,M+2-K)-FI(IP,M+2-K))/(H1

H2)

END IF

END DO

RETURN

END

SUBROUTINE MEGP (XP, YP, IM, JM, L, H1, H2)

C OBSERVATION POINTS

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION IM(L), JM(L), XP(L), YP(L)

XP(1) = 0.17D0

Section 7.5 Continuation of non-stationary ﬁelds from point observation data 339

YP(1) = 0.53D0

XP(2) = 0.36D0

YP(2) = 0.87D0

XP(3) = 0.42D0

YP(3) = 0.39D0

XP(4) = 0.58D0

YP(4) = 0.48D0

XP(5) = 0.83D0

YP(5) = 0.25D0

XP(6) = 0.11D0

YP(6) = 0.15D0

XP(7) = 0.76D0

YP(7) = 0.71D0

XP(8) = 0.28D0

YP(8) = 0.33D0

XP(9) = 0.35D0

YP(9) = 0.65D0

XP(10) = 0.49D0

YP(10) = 0.24D0

DO IP = 1,L

IM(IP) = XP(IP)/H1 + 1

IF ((XP(IP) - IM(IP)

H1).GT. 0.5D0

H1) IM(IP) = IM(IP)+1

JM(IP) = YP(IP)/H2 + 1

IF ((YP(IP) - JM(IP)

H2).GT. 0.5D0

H2) JM(IP) = JM(IP)+1

END DO

RETURN

END

SUBROUTINE PU (IM, JM, U, N1, N2, F, L)

C SOLUTION AT THE MEASUREMENT POINTS

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION IM(L), JM(L), F(L), U(N1,N2)

DO IP = 1,L

I = IM(IP)

J = JM(IP)

F(IP) = U(I,J)

END DO

RETURN

END

SUBROUTINE BS ( U, A0, A1, A2, B, N1, N2 )

C FUNCTIONAL GRADIENT

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION U(N1,N2), A0(N1,N2), A1(N1,N2), A2(N1,N2), B(N1,N2)

DOJ=2,N2-1

DOI=2,N1-1

B(I,J) = A0(I,J)

U(I,J)

- A1(I-1,J)

U(I-1,J)

- A1(I,J)

U(I+1,J)

- A2(I,J-1)

U(I,J-1)

- A2(I,J)

U(I,J+1)

END DO

340 Chapter 7 Evolutionary inverse problems

END DO

RETURN

END

We have restricted ourselves to the case in which in equation (7.257) we have k(x) = 1

and, in addition, purely implicit schemes are used (the difference-scheme coefﬁcients

for the ground state are generated in the subroutine FDS, and for the conjugate state,

in the subroutines FDS and RHS.

7.5.6 Computational experiments

The basic grid was the uniform grid with h

= 0.02 and h

= 0.02 for the problem in

unit square. In the realization of the quasi-real experiment, as input data, the solution

of the direct problem at some points of the calculation domain was taken. The problem

is solved till the time T = 0.025, the time step size of the grid being τ = 0.00025. In

the direct problem, the initial condition is given in the form

(x, 0) =



1,(x

− 0.6)

+ (x

− 0.6)

≤ 0.04,

0,(x

− 0.6)

+ (x

− 0.6)

> 0.04.

Figure 7.18 Observation points and initial conditions in the direct problem

Figure 7.18 shows the location of the observation points and indicates the portion

of the calculation domain in which the initial condition was localized. As input data,

the solution at the observation points is used (see Figure 7.19).

Section 7.5 Continuation of non-stationary ﬁelds from point observation data 341

Figure 7.19 Solution at the observation points

Figure 7.20 Inverse-problem solution obtained with δ = 0.02

342 Chapter 7 Evolutionary inverse problems

Consider calculated data that illustrate possibilities in reconstructing the initial con-

dition from point observation data. Figure 7.20 shows the solution of the inverse prob-

lem obtained with the inaccuracy level δ = 0.02 (plotted are contour lines following

with the grid size u = 0.1). The effect due to the inaccuracy level is illustrated in

Figures 7.21 and 7.22.

Figure 7.21 Inverse-problem solution obtained with δ = 0.04

Figure 7.22 Inverse-problem solution obtained with δ = 0.01

Section 7.6 Exercises 343

7.6 Exercises

Exercise 7.1 Examine convergence of the approximate solution found by formulas

(7.14), (7.36) to the exact solution of problem (7.6), (7.7) H

, D = S

Exercise 7.2 Formulate the non-local problem that arises in the optimal control prob-

lem (7.38)–(7.41) with a non-self-adjoint operator A.

Exercise 7.3 Using the program PROBLEM10, examine the time dependence of the

approximate-solution inaccuracy. Compare the experimental data with the theoretical

results (see estimate (7.35)).

Exercise 7.4 Examine the generalized inverse method (7.109), (7.138) to approxi-

mately solve the ill-posed problem (7.101), (7.102).

Exercise 7.5 Consider the additive scheme of component-by-component splitting (to-

tal approximation scheme)

n+1/2

− y

− 

+ α

(σ y

n+1/2

+ (1 − σ)y

) = 0,

n+1

− y

n+1/2

− 

n+1/2

+ α

(σ y

n+1

+ (1 − σ)y

n+1/2

) = 0

for the Cauchy problem posed for equation (7.158).

Exercise 7.6 Based on the data calculated by the program PROBLEM11, examine the

effect of exact-solution smoothness on the initial-condition reconstruction inaccuracy.

Exercise 7.7 Consider the possibility of using iteration methods constructed around

the minimization of the discrepancy functional in solving the inverted-time problem

(7.159)–(7.161).

Exercise 7.8 Construct an iteration method for reﬁning the initial condition in scheme

(7.166), (7.167) as applied to the approximate solution of the retrospective inverse

problem with a non-self-adjoint, positively deﬁned operator .

Exercise 7.9 Using the program PROBLEM12, experimentally examine how the

choice of B in (7.179) affects the possibility of identifying a solution with desired

smoothness.

Exercise 7.10 Construct a difference scheme for the non-local boundary value prob-

lem

−



k(x)



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

u(0) + αu(l) = μ

−k(0)

(0) = μ

Modify the sweep algorithm so that to make it appropriate for the solution of the related

non-local difference problem.

344 Chapter 7 Evolutionary inverse problems

Exercise 7.11 Construct the difference scheme for the generalized inverse method

(7.220)–(7.222) as applied to the approximate solution of the ill-posed problem

(7.187)–(7.189) and perform a stability study for this scheme.

Exercise 7.12 Perform computational experiments (program PROBLEM13) to inves-

tigate into the solution inaccuracy in the Cauchy problem for the Laplace equation as

dependent on the smoothness of exact initial conditions.

Exercise 7.13 Formulate optimality conditions for the minimization problem for

(7.274), (7.275) under constraints (7.270), (7.271).

Exercise 7.14 Derive optimality conditions for the difference problem (7.271)–

(7.274), (7.275) in the case of a non-self-adjoint operator >0.

Exercise 7.15 Modify the program PROBLEM14 so that the initial condition be re-

ﬁned by the method of minimal discrepancies instead of the simple iteration method.

8 Other problems

Above, two classes of inverse problems for mathematical physics equations have been

considered in which it was required to identify the right-hand side of an equation or the

initial condition for the equation. Among other problems important for applications,

boundary value inverse problems deserve mention in which to be reconstructed are the

boundary conditions. For approximate solution of the latter problems, methods using

some perturbation of the equation or methods using non-locally perturbed boundary

conditions can be applied. In the case in which the generalized inverse method is used

to solve some boundary value inverse problem, namely, in treating the spatial coordi-

nate as the evolutionary coordinate, special emphasis is to be placed on the hyperbolic

regularization method, used to pass from the hyperbolic to a parabolic equation. In the

present chapter, possibilities offered by the generalized inverse method as applied to

problems with perturbed boundary conditions are discussed with the example of the

boundary value inverse problem for the one-dimensional parabolic equation of sec-

ond order. For a more general two-dimensional problem, an algorithm with iteratively

reﬁned boundary condition is used. The problems most difﬁcult for examination are

coefﬁcient inverse problems for mathematical physics equations. Here, we have re-

stricted ourselves to the matter of numerical solution of two coefﬁcient problems. In

the ﬁrst problem it is required to determine the higher coefﬁcient as a function of the

solution for a one-dimensional parabolic equation. We describe a computational al-

gorithm that solves the coefﬁcient inverse problem for the two-dimensional elliptic

equation in the case in which the unknown coefﬁcient does not depend on one of the

two coordinates.

8.1 Continuation over the spatial variable in the boundary

value inverse problem

In this section, we consider the boundary value inverse problem for the one-

dimensional parabolic equation of second order (heat conduction equation). In this

problem, it is required to reconstruct the boundary condition from measurements per-

formed inside the calculation domain. This problem belongs to the class of condition-

ally well-posed problems and, to be solved stably, it requires the use of regularization

methods. Here, the generalized inverse method is to be applied under conditions in

which the problem is considered as an evolutionary one with respect to the spatial

variable. The use of the generalized inverse method leads, in particular, to the well-

known hyperbolic regularization of boundary value inverse problems.