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

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

406 Chapter 8 Other problems

C NODE CLOSEST TO THE OBSERVATION POINT XD

ND = (XD-XL)/H+1

IF ( (XD-XL - (ND-1)

H).GT.0.5

H)ND=ND+1

C DIRECT PROBLEM

C BOUNDARY CONDITION

DOK=1,M

T=K

TAU

BR(K) = AG(T)

END DO

C INITIAL CONDITION

T = 0.D0

DOI=1,N+1

Y(I) = 0.D0

U(I,1) = Y(I)

END DO

C NEXT TIME LAYER

DOK=1,M

T=K

TAU

C DIFFERENCE-SCHEME COEFFICIENTS

C PURELY IMPLICIT LINEARIZED SCHEME

DOI=2,N

U1 = (Y(I) + Y(I-1)) / 2

U2 = (Y(I+1) + Y(I)) / 2

A(I) = AK(U1) / (H

B(I) = AK(U2) / (H

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

F(I) = Y(I) / TAU

END DO

C BOUNDARY CONDITION AT THE LEFT END POINT

B(1) = 0.D0

C(1) = 1.D0

F(1) = 0.D0

C BOUNDARY CONDITION AT THE RIGHT END POINT

A(N+1) = 0.D0

C(N+1) = 1.D0

F(N+1) = AG(T)

C SOLUTION OF THE PROBLEM ON THE NEXT TIME LAYER

CALL PROG ( N+1, A, C, B, F, ALPHA, BETA, Y )

DOI=1,N+1

U(I,K+1) = Y(I)

END DO

C SOLUTION AT THE OBSERVATION POINT

Section 8.4 Coefﬁcient inverse problem for the nonlinear parabolic equation 407

DOK=1,M+1

PHI(K) = U(ND,K)

PHID(K) = PHI(K)

END DO

C DISTURBING OF MEASURED VALUES

DOK=2,M+1

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

DELTA

(RAND(0)-0.5)

END DO

C INVERSE PROBLEM

C CHOICE OF THE STEP IN THE PIECEWISE-CONSTANT APPROXIMATION

L=1

100 CONTINUE

ML=M/L

DO LK = 1, L+1

T = (LK-1)

TAU

UL(LK) = AG(T)

END DO

C INITIAL CONDITION

T = 0.D0

SD = 0.D0

DOI=1,N+1

Y1(I) = 0.D0

END DO

DOLK=1,L

C DETERMINATION OF THE UNKNOWN COEFFICIENT OVER THE SUB-INTERVAL

C ITERATION GOLDEN-SECTION METHOD

AS = 0.1D0

BS = 10.D0

R1 = (DSQRT(5.0D0)-1.D0)/2

R2=R1

HS=BS-AS

C SOLUTION OF THE PROBLEM OVER THE SUB-INTERVAL

C AT A GIVEN COEFFICIENT

CALL STEPB ( N, ND, M, ML, LK, H, TAU, A, C, B, F, AL, BET

+ ,Y, Y1, AS, AKS, UL, PHID, YA )

CALL STEPB ( N, ND, M, ML, LK, H, TAU, A, C, B, F, AL, BET

+ ,Y, Y1, BS, AKS, UL, PHID, YB )

CS=AS+R2

DS=AS+R1

CALL STEPB ( N, ND, M, ML, LK, H, TAU, A, C, B, F, AL, BET

+ ,Y, Y1, CS, AKS, UL, PHID, YC )

CALL STEPB ( N, ND, M, ML, LK, H, TAU, A, C, B, F, AL, BET

+ ,Y, Y1, DS, AKS, UL, PHID, YD )

KS=0

200 KS=KS+1

IF (YC.LT.YD) THEN

BS=DS

YB=YD

408 Chapter 8 Other problems

DS=CS

YD=YC

HS=BS-AS

CS=AS+R2

CALL STEPB ( N, ND, M, ML, LK, H, TAU, A, C, B, F, AL, BET

+ ,Y, Y1, CS, AKS, UL, PHID, YC )

ELSE

AS=CS

YA=YC

CS=DS

YC=YD

HS=BS-AS

DS=AS+R1

CALL STEPB ( N, ND, M, ML, LK, H, TAU, A, C, B, F, AL, BET

+ ,Y, Y1, DS, AKS, UL, PHID, YD )

END IF

IF (DABS(YB-YA).GT.EPSA.OR.HS.GT.EPSH) GO TO 200

C FOUND SOLUTION

PS=AS

YP=YA

IF (YB.LT.YA) THEN

PS=BS

YP=YB

ENDIF

DOI=1,N+1

Y1(I) = Y(I)

END DO

C TOTAL INACCURACY

SD=SD+YP

END DO

WRITE(01,

) L, KS, SD

C EXIT FROM THE PROGRAM BY THE DISCREPANCY CRITERION

L=L

IF (L.EQ.2) THEN

SD0=SD

IF (SD.GT.TMAX

DELTA

2 .AND. L .LE. M) GO TO 100

ELSE

IF (SD.LT.SD0) THEN

IF (SD.GT.TMAX

DELTA

2 .AND. L .LE. M) GO TO 100

END IF

C RECORDING OF THE SOLUTION IN A FILE

WRITE ( 01,

) ((U(I,K),I=1,N+1), K=1,M+1)

WRITE ( 01,

) (PHI(K), K=1,M+1)

WRITE ( 01,

) (PHID(K), K=1,M+1)

WRITE ( 01,

) (AKS(K), K=1,M+1)

CLOSE ( 01 )

STOP

END

SUBROUTINE PROG ( N, A, C, B, F, AL, BET, Y )

Section 8.4 Coefﬁcient inverse problem for the nonlinear parabolic equation 409

IMPLICIT REAL

8 ( A-H, O-Z )

C SWEEP METHOD

DIMENSION A(N), C(N), B(N), F(N), Y(N), AL(N+1), BET(N+1)

AL(1) = B(1) / C(1)

BET(1) = F(1) / C(1)

DOI=2,N

SS = C(I) - AL(I-1)

A(I)

AL(I) = B(I) / SS

BET(I) = ( F(I) + BET(I-1)

A(I))/SS

END DO

Y(N) = BET(N)

DO I = N-1, 1, -1

Y(I) = AL(I)

Y(I+1) + BET(I)

END DO

RETURN

END

DOUBLE PRECISION FUNCTION AK(U)

IMPLICIT REAL

8 ( A-H, O-Z )

C COEFFICIENT AT HIGHER DERIVATIVES

AK = 0.5D0 + U

RETURN

END

DOUBLE PRECISION FUNCTION AG(T)

IMPLICIT REAL

8 ( A-H, O-Z )

C BOUNDARY CONDITION AT THE RIGHT END POINT

AG=T

RETURN

END

SUBROUTINE STEPB ( N, ND, M, ML, LK, H, TAU, A, C, B, F, AL, BET

+ ,Y, Y1, AK, AKS, UL, PHI, FD )

IMPLICIT REAL

8 ( A-H, O-Z )

C CALCULATION OF THE DISCREPANCY ON THE SUBINTERVAL

C AT A GIVEN CONSTANT COEFFICIENT

DIMENSION A(N+1), C(N+1), B(N+1), F(N+1), Y(N+1), Y1(N+1)

+ ,ALPHA(N+2), BETA(N+2), PHI(M+1), AKS(M+1), UL(M+1)

AKS(LK+1) = AK

FD = 0.D0

DOI=1,N+1

Y(I) = Y1(I)

END DO

T = (LK-1)

TAU

DOK=1,ML

T=T+TAU

C DIFFERENCE-SCHEME COEFFICIENTS

410 Chapter 8 Other problems

C PURELY IMPLICIT SCHEME

DOI=2,N

U1 = (Y(I) + Y(I-1)) / 2

U2 = (Y(I+1) + Y(I)) / 2

A(I) = AFK(U1, LK+1, UL, AKS, M) / (H

B(I) = AFK(U2, LK+1, UL, AKS, M) / (H

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

F(I) = Y(I) / TAU

END DO

C BOUNDARY CONDITION AT THE LEFT END POINT

B(1) = 0.D0

C(1) = 1.D0

F(1) = 0.D0

C BOUNDARY CONDITION AT THE RIGHT END POINT

A(N+1) = 0.D0

C(N+1) = 1.D0

F(N+1) = AG(T)

C SOLUTION OF THE PROBLEM ON THE NEXT TIME LAYER

CALL PROG ( N+1, A, C, B, F, ALPHA, BETA, Y )

C DISCREPANCY AT THE OBSERVATION POINT

FD = FD + (Y(ND) - PHI((LK-1)

ML+K+1))

TAU

END DO

RETURN

END

DOUBLE PRECISION FUNCTION AFK ( U, L, UL, AL, M )

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION UL(M+1), AL(M+1)

C PIECEWISE-CONSTANT COEFFICIENT

AFK = AL(2)

DOIL=2,L

IF (U.LE.UL(IL) .AND. U.GT.UL(IL-1)) AFK = AL(IL)

END DO

RETURN

END

In the subroutine-function AK, the nonlinear coefﬁcient for the direct problem is

speciﬁed; the data obtained by solving this program are used as input data for the

inverse problem. The subroutine-function AG calculates the boundary conditions at

the right boundary. The piecewise constant completion is effected in AFK.

Section 8.4 Coefﬁcient inverse problem for the nonlinear parabolic equation 411

8.4.6 Computational experiments

For input-data setting, we solve the direct problem for l = 1, T = 1 on the uniform

grid with h = 0.01 and τ = 2

−7

. Here, the simplest linearized difference scheme is

used in which the nonlinear coefﬁcient was taken from the previous time layer, i. e.,

n+1

− y

+ A(y

n+1

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

− 1.

Below, calculation data are given, obtained for the case in which the exact solution

was obtained for

k(u) = 0.5 + u,

so that the coefﬁcient to be found exhibits three-fold variation in the calculation do-

main. The direct-problem solution is shown in Figure 8.22.

In solving the inverse problem, we have restricted ourselves to one observation point

and to the case in which

= 0.3, M = 1,ϕ

) = y

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

where y

(x) is the mesh solution of the direct problem. These data were perturbed

with a normally distributed function:

) = ϕ

) + 2δ(σ(t

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

The solution of the inverse problem obtained at the inaccuracy level deﬁned by

δ = 0.02 is shown in Figure 8.23. Considering the fact that the boundary condition

was set in the form g(t) = t, in one and the same graph we have plotted the exact and

approximate input data for the inverse problem (the solution at the point x = z

= 0.6

and the sought coefﬁcient k(u)). At the chosen level of inaccuracy, for the coefﬁcient

to be reconstructed it sufﬁces to have two subintervals ( p = 2 in (8.147), (8.156)).

The effect due to inaccuracies can be ﬁgured out considering Figures 8.24 and 8.25.

412 Chapter 8 Other problems

Figure 8.22 Direct-problem solution

Figure 8.23 Inverse-problem solution obtained with δ = 0.02

Section 8.4 Coefﬁcient inverse problem for the nonlinear parabolic equation 413

Figure 8.24 Inverse-problem solution obtained with δ = 0.04

Figure 8.25 Inverse-problem solution obtained with δ = 0.01

414 Chapter 8 Other problems

Figure 8.26 Calculation domain

8.5 Coefﬁcient inverse problem for elliptic equation

In this section, we consider a problem in which it is required to ﬁnd the lowest coef-

ﬁcient in the second-order elliptic equation from data given at the calculation-domain

boundary. We assume that the unknown coefﬁcient is independent of one of the vari-

ables. With this simplest example, we illustrate the possibility to examine the solution

unicity for coefﬁcient inverse problems. The computational algorithm is constructed

for a model problem in a rectangle.

8.5.1 Statement of the problem

The problem is considered in its simplest form, with the calculation domain being a

rectangle:

 ={x | x = (x

, x

), 0 < x

< l

,α= 1, 2}.

For individual segments of the boundary of , we use the following settings (see

Figure 8.26):

∂ = 

∪ 

As usually, we start from the formulation of the direct problem. We assume that the

function u(x), x = (x

, x

) satisﬁes the equation

−u + c(x

)u = 0, x ∈ , (8.165)

where

u ≡



α=1

∂

∂x

Section 8.5 Coefﬁcient inverse problem for elliptic equation 415

Equation (8.165) is supplemented with boundary conditions of the ﬁrst kind:

u(x) = ϕ(x), x ∈ ∂. (8.166)

In the direct problem (8.165), (8.166), the lowest coefﬁcient is assumed to depend just

on x

In the inverse problem of interest, the coefﬁcient c(x

) is unknown. This coefﬁcient

is to be determined from some additional data. We consider the case in which the

additional data are obtained from measurements performed on the domain boundary

in the form

∂u

∂n

(x) = ψ(x), x ∈ ∂, (8.167)

where n is the external normal to .

Previously, we considered a linear inverse problem in which it was required to

determine the unknown right-hand side independent of one of the variables. The

posed inverse problem (8.165)–(8.167), in which it is required to ﬁnd a function pair

{u(x), c(x

)}, is a more complex problem. The difﬁculties in the consideration of this

problem largely result from its nonlinearity.

8.5.2 Solution uniqueness for the inverse problem

We can speak of a certain overdetermination of the coefﬁcient inverse problem under

consideration. In problem (8.165)–(8.167), identiﬁcation of the unknown coefﬁcient

c(x

) implies determination of a function of one variable over the interval [0, l

] from

two functions of the variable x

(the function ψ(x) on the sides 

, 

) and from

two functions of the variable x

(the function ψ(x) on the sides 

, 

). Hence, in

general it would be reasonable to pose the question about sufﬁcient additional data for

determining the coefﬁcient c(x

To examine solution uniqueness for the inverse problem (8.165)–(8.167), one can

follow the traditional approach that is normally used in the cases of direct nonlinear

boundary value mathematical physics problems and nonlinear inverse problems. Sup-

pose that there exist two solutions of the inverse problem (8.165)–(8.167); we denote

these solutions as {u

(x), c

), β = 1, 2}, i. e.,

−u

+ c

= 0, x ∈ , (8.168)

(x) = ϕ(x), x ∈ ∂, (8.169)

∂u

∂n

(x) = ψ(x), x ∈ ∂, β = 1, 2. (8.170)

For the differences

v(x) = u

(x) − u

(x), θ (x

) = c

) − c