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

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

226 Chapter 6 Right-hand side identiﬁcation

Solvability of this system is controlled by the determinant

D = (1 − s

(1)

))(1 − s

(2)

− h

)) − s

(2)

(1)

− h

whose non-zero value can easily be guaranteed under certain constraints. Taking the

inequalities λ>0 and q

≥ 0, α = 1, 2 into account, we have: s

(α)

) ≥ 0,

0 ≤ x

≤ l

. Hence, we have D > 0 at sufﬁciently small l

, for instance.

As a matter of fact, the determinant of (6.177), (6.178) is always positive. To show

this, we have to recall (see (6.164)) the expressions for the mesh functions q

α = 1, 2:

) =



1 −



, q

) =

By virtue of this, for the solutions of the boundary value problems (6.173), (6.174) and

(6.175), (6.176) there hold the relations

(1)

) = s

(2)

− x

), 0 ≤ x

≤ l

The latter means, in particular, that it is unnecessary for us to solve (in the case of the

uniform computational grid used) one of the two difference boundary value problems,

(6.173), (6.174) or (6.175), (6.176). Hence,

D = 1 − 2s

(1)

) + (s

(1)

))

− (s

(1)

− h

))

and, with allowance for s

(1)

)>s

(1)

− h

), we obtain that D > 1.

With the mesh functions s(x

), s

(α)

,α= 1, 2 found from (6.171)–(6.176) and with

the mesh functions w(h

) and w(l

− h

) found from (6.177), (6.178), the solution of

problem (6.168), (6.169) can be found in the form (6.170). The program realization of

this algorithm is discussed below.

6.5.5 Program

The used computational algorithm employs the fast Fourier transform over the vari-

able x

. In the embodied algorithm, one can conveniently use the non-normalized

eigenfunctions

sin

kπ x

= sin

kπi

, x

= i

In the latter case, any mesh function f

deﬁned at i

= 1, 2,...,N

− 1 (on ω

) can

be expanded as

−1



k=1

sin

kπi

, i

= 1, 2,...,N

− 1,

where

−1



sin

kπi

, k = 1, 2,...,N

− 1.

Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 227

Let us use the subroutine SINT from the program library SLATEC. Here, we cite a

borrowed description of the subroutine, representing a component of the big program

package FFTPACK that embodies fast Fourier transform algorithms. You can use other

subroutines for fast Fourier transform available at hand.

Subroutine SINT

DECK SINT

SUBROUTINE SINT (N, X, WSAVE)

***

BEGIN PROLOGUE SINT

***

PURPOSE Compute the sine transform of a real, odd sequence.

***

LIBRARY SLATEC (FFTPACK)

***

CATEGORY J1A3

***

TYPE SINGLE PRECISION (SINT-S)

***

KEYWORDS FFTPACK, FOURIER TRANSFORM

***

AUTHOR Swarztrauber, P. N., (NCAR)

***

DESCRIPTION

C Subroutine SINT computes the discrete Fourier sine transform

C of an odd sequence X(I). The transform is defined below at

C output parameter X.

C SINT is the unnormalized inverse of itself since a call of SINT

C followed by another call of SINT will multiply the input sequence

C X by 2

(N+1).

C The array WSAVE which is used by subroutine SINT must be

C initialized by calling subroutine SINTI(N,WSAVE).

C Input Parameters

C N the length of the sequence to be transformed. The method

C is most efficient when N+1 is the product of small primes.

C X an array which contains the sequence to be transformed

C WSAVE a work array with dimension at least INT(3.5

N+16)

C in the program that calls SINT. The WSAVE array must be

C initialized by calling subroutine SINTI(N,WSAVE), and a

C different WSAVE array must be used for each different

C value of N. This initialization does not have to be

C repeated so long as N remains unchanged. Thus subsequent

C transforms can be obtained faster than the first.

C Output Parameters

C X For I=1,\dots,N

C X(I)= the sum from K=1 to K=N

X(K)

SIN(K

PI/(N+1))

C A call of SINT followed by another call of

228 Chapter 6 Right-hand side identiﬁcation

C SINT will multiply the sequence X by 2

(N+1).

C Hence SINT is the unnormalized inverse

C of itself.

C WSAVE contains initialization calculations which must not be

C destroyed between calls of SINT.

***

REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel

C Computations (G. Rodrigue, ed.), Academic Press,

C 1982, pp. 51-83.

***

ROUTINES CALLED RFFTF

***

REVISION HISTORY (YYMMDD)

C 790601 DATE WRITTEN

C 830401 Modified to use SLATEC library source file format.

C 860115 Modified by Ron Boisvert to adhere to Fortran 77 by

C (a) changing dummy array size declarations (1) to (

C (b) changing definition of variable SQRT3 by using

C FORTRAN intrinsic function SQRT instead of a DATA

C statement.

C 881128 Modified by Dick Valent to meet prologue standards.

C 891009 Removed unreferenced statement label. (WRB)

C 891009 REVISION DATE from Version 3.2

C 891214 Prologue converted to Version 4.0 format. (BAB)

C 920501 Reformatted the REFERENCES section. (WRB)

***

END PROLOGUE SINT

The auxiliary subroutine SINTI initializes the subroutine SINT.

Subroutine SINTI

DECK SINTI

SUBROUTINE SINTI (N, WSAVE)

***

BEGIN PROLOGUE SINTI

***

PURPOSE Initialize a work array for SINT.

***

LIBRARY SLATEC (FFTPACK)

***

CATEGORY J1A3

***

TYPE SINGLE PRECISION (SINTI-S)

***

KEYWORDS FFTPACK, FOURIER TRANSFORM

***

AUTHOR Swarztrauber, P. N., (NCAR)

***

DESCRIPTION

C Subroutine SINTI initializes the array WSAVE which is used in

C subroutine SINT. The prime factorization of N together with

C a tabulation of the trigonometric functions are computed and

C stored in WSAVE.

C Input Parameter

C N the length of the sequence to be transformed. The method

C is most efficient when N+1 is a product of small primes.

C Output Parameter

C WSAVE a work array with at least INT(3.5

N+16) locations.

C Different WSAVE arrays are required for different values

Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 229

C of N. The contents of WSAVE must not be changed between

C calls of SINT.

***

REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel

C Computations (G. Rodrigue, ed.), Academic Press,

C 1982, pp. 51-83.

***

ROUTINES CALLED RFFTI

***

REVISION HISTORY (YYMMDD)

C 790601 DATE WRITTEN

C 830401 Modified to use SLATEC library source file format.

C 860115 Modified by Ron Boisvert to adhere to Fortran 77 by

C (a) changing dummy array size declarations (1) to (

C (b) changing references to intrinsic function FLOAT

C to REAL, and

C (c) changing definition of variable PI by using

C FORTRAN intrinsic function ATAN instead of a DATA

C statement.

C 881128 Modified by Dick Valent to meet prologue standards.

C 890531 Changed all specific intrinsics to generic. (WRB)

C 890531 REVISION DATE from Version 3.2

C 891214 Prologue converted to Version 4.0 format. (BAB)

C 920501 Reformatted the REFERENCES section. (WRB)

***

END PROLOGUE SINTI

The program PROBLEM9 is used to obtain the approximate solution of the inverse

problem.

Program PROBLEM9

C PROBLEM9 IDENTIFICATION OF THE RIGHT-HAND SIDE

C OF THE POISSON EQUATION IN A RECTANGLE

PARAMETER ( DELTA = 0.0, N1 = 257, N2 = 257 )

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

+ ,FIT1(N2), FIT2(N2), FI1(N2), FI2(N2)

+ ,AM1(N2), AM2(N2), AMD1(N2), AMD2(N2)

+ ,A(N1), B(N1), C(N1), FF(N1), YY(N1), S(N1), S1(N1)

+ ,WSAVE(4

N2), WW(N2-2)

C PARAMETERS:

C X1L, X2L - COORDINATES OF THE LEFT BOTTOM CORNER

C OF THE RECTANGULAR CALCULATION DOMAIN;

C X1R, X2R - COORDINATES OF THE RIGHT TOP CORNER;

C N1, N2 - NUMBER OF DRID NODES OVER THE CORRESPONDING

C DIRECTIONS;

C U(N1,N2) - DIFFERENCE SOLUTION OF THE DIRECT PROBLEM;

C FI(N2),

C FI2(N2) - EXACT RIGHT-HAND SIDE COMPONENTS;

C AM1(N2),

230 Chapter 6 Right-hand side identiﬁcation

C AM2(N2) - BOUNDARY CONDITIONS;

C DELTA - INPUT-DATA INACCURACY LEVEL;

C AMD1(N2),

C AMD2(N2) - DISTURBED BOUNDARY CONDITIONS;

C FIS1(N2),

C FIS2(N2) - FOUND RIGHT-HAND SIDE COMPONENTS;

X1L = 0.

X1R = 1.

X2L = 0.

X2R = 1.

PI = 3.1415926

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

C GRID

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

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

DOI=1,N1

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

END DO

DOJ=1,N2

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

END DO

C DIRECT PROBLEM

C EXACT RIGHT-HAND SIDE

DOJ=2,N2-1

FIT1(J) = AF1(X2(J))

FIT2(J) = AF2(X2(J))

DOI=2,N1-1

F(I,J) = FIT1(J) + X1(I)

FIT2(J)

END DO

C FORWARD FOURIER TRANSFORM

CALL SINTI(N2-2,WSAVE)

DOI=2,N1-1

DOJ=2,N2-1

WW(J-1) = F(I,J)

END DO

CALL SINT(N2-2,WW,WSAVE)

DOJ=2,N2-1

F(I,J) = WW(J-1)

END DO

C SOLUTION OF THE ONE-DIMENSIONAL (OVER THE VARIABLE X1) PROBLEM

DOJ=2,N2-1

C DIFFERENCE-SCHEME COEFFICIENTS

ALAM = 4./H2

(SIN(PI

(J-1)/(2.

(N2-1))))

DOI=2,N1-1

Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 231

A(I) = 1. / (H1

H1)

B(I) = 1. / (H1

H1)

C(I) = A(I) + B(I) + ALAM

FF(I) = F(I,J)

END DO

C BOUNDARY CONDITIONS AT THE LEFT AND RIGHT END POINTS

B(1) = 0.

C(1) = 1.

FF(1) = 0.

A(N1) = 0.

C(N1) = 1.

FF(N1) = 0.

C SOLUTION OF THE ONE-DIMENSIONAL PROBLEM

ITASK = 1

CALL PROGS3 ( N1, A, C, B, FF, YY, ITASK )

DOI=1,N1

F(I,J) = YY(I)

END DO

C INVERSE FOURIER TRANSFORM

DOI=2,N1-1

DOJ=2,N2-1

WW(J-1) = F(I,J)

END DO

CALL SINT(N2-2,WW,WSAVE)

DOJ=2,N2-1

U(I,J) = 1./(2.

(N2-1))

WW(J-1)

END DO

C BOUNDARY CONDITIONS IN THE INVERSE PROBLEM

DOJ=1,N2

AM1(J) = U(2,J) / H1

AM2(J) = - U(N1-1,J) / H1

END DO

C DISTURBING OF MEASURED QUANTITIES

DOJ=1,N2

AMD1(J) = AM1(J) + 2.

DELTA

(RAND(0)-0.5)

AMD2(J) = AM2(J) + 2.

DELTA

(RAND(0)-0.5)

END DO

C INVERSE PROBLEM

C RIGHT-HAND SIDE OF THE LOADED EQUATION

DOI=2,N1-1

Q1=2.

(N1-I) / (X1R - X1L)

232 Chapter 6 Right-hand side identiﬁcation

Q2=2.

(I-1) / (X1R - X1L)

DOJ=2,N2-1

F(I,J) = Q1

AMD1(J) - Q2

AMD2(J)

END DO

C FORWARD FOURIER TRANSFORM

CALL SINTI(N2-2,WSAVE)

DOI=2,N1-1

DOJ=2,N2-1

WW(J-1) = F(I,J)

END DO

CALL SINT(N2-2,WW,WSAVE)

DOJ=2,N2-1

F(I,J) = WW(J-1)

END DO

C SOLUTION OF THE ONE-DIMENSIONAL (OVER THE VARIABLE X1) PROBLEM

DOJ=2,N2-1

C AUXILIARY FUNCTIONS

C DIFFERENCE-SCHEME COEFFICIENTS

ALAM = 4./H2

(SIN(PI

(J-1)/(2.

(N2-1))))

DOI=2,N1-1

A(I) = 1. / (H1

H1)

B(I) = 1. / (H1

H1)

C(I) = A(I) + B(I) + ALAM

FF(I) = F(I,J)

END DO

C BOUNDARY CONDITIONS AT THE LEFT AND RIGHT END POINTS

B(1) = 0.

C(1) = 1.

FF(1) = 0.

A(N1) = 0.

C(N1) = 1.

FF(N1) = 0.

C SOLUTION OF THE ONE-DIMENSIONAL PROBLEM

ITASK = 1

CALL PROGS3 ( N1, A, C, B, FF, S, ITASK )

DOI=2,N1-1

FF(I)=-2.

(N1-I) / ((N1-1.)

H1)

END DO

ITASK = 2

CALL PROGS3 ( N1, A, C, B, FF, S1, ITASK )

C SOLUTION OF THE SYSTEM OF TWO EQUATIONS

Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 233

DD = (1. - S1(2))

2 - (S1(N1-1))

W1 = ((1. - S1(2))

S(2) + S1(N1-1)

S(N1-1))/DD

W2 = ((1. - S1(2))

S(N1-1) + S1(N1-1)

S(2))/DD

DOI=2,N1-1

F(I,J) = S(I) + W1

S1(I) + W2

S1(N1+1-I)

END DO

C INVERSE FOURIER TRANSFORM

DOI=2,N1-1

DOJ=2,N2-1

WW(J-1) = F(I,J)

END DO

CALL SINT(N2-2,WW,WSAVE)

DOJ=2,N2-1

Y(I,J) = 1./(2.

(N2-1))

WW(J-1)

END DO

C RIGHT-HAND SIDE

DOI=2,N1-1

DOJ=2,N2-1

F(I,J) = - (Y(I+1,J) - 2.

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

H1)

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

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

H2)

END DO

DOJ=2,N2-1

FI1(J) = 2.

AMD1(J)/H1 - 2.

Y(2,J)/(H1

H1)

FI2(J) = - 2.

(AMD2(J)+AMD1(J))/((X1R-X1L)

H1)

+-2.

(Y(N1-1,J)-Y(2,J))/((X1R-X1L)

H1)

END DO

WRITE ( 01,

) (FI1(J),J=2,N2-1)

WRITE ( 01,

) (FIT1(J),J=2,N2-1)

WRITE ( 01,

) (FI2(J),J=2,N2-1)

WRITE ( 01,

) (FIT2(J),J=2,N2-1)

CLOSE ( 01 )

STOP

END

FUNCTION AF1 ( X2 )

C FIRST COMPONENT OF THE RIGHT-HAND SIDE OF THE EQUATION

AF1=1.

IF (X2.GT.0.5) AF1 = 0.

RETURN

END

FUNCTION AF2 ( X2 )

C SECOND COMPONENT OF THE RIGHT-HAND SIDE OF THE EQUATION

AF2=X2

234 Chapter 6 Right-hand side identiﬁcation

RETURN

END

SUBROUTINE PROGS3 ( N, A, C, B, F, Y, ITASK )

C SWEEP METHOD

C FOR TRIDIAGONAL MATRIX

C ITASK = 1: FACTORIZATION AND SOLUTION;

C ITASK = 2: ONLY SOLUTION

DIMENSION A(N), C(N), B(N), F(N), Y(N)

IF ( ITASK .EQ. 1 ) THEN

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

DOI=2,N

C(I) = C(I) - B(I-1)

A(I)

B(I) = B(I) / C(I)

END DO

ITASK = 2

END IF

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

DOI=2,N

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

A(I) ) / C(I)

END DO

Y(N) = F(N)

DO I = N-1, 1, -1

Y(I) = B(I)

Y(I+1) + F(I)

END DO

RETURN

END

The function of the additional programs requires no comments.

6.5.6 Computational experiments

The presented program PROBLEM9 solves the inverse problem (6.144), (6.145),

(6.150), in which to be found is the right-hand side (6.147) with

) =



1, 0 < x

< 0.5,

0, 0.5 < x

< 1,

) = x

The problem is solved in the unit square (l

= l

= 1). To obtain the input data

for the inverse problem, the direct problem (6.144), (6.145) is preliminarily solved at

a given right-hand side.

First of all, consider the solution data obtained for the inverse problem with un-

perturbed input data. Of interest here are numerical data obtained on a sequence of

progressively reﬁned grids (see Figures 6.18–6.20). The approximate solution is seen

Section 6.5 Right-hand side reconstruction from boundary data: elliptic equations 235

to converge to the exact solution. A sufﬁciently high accuracy can be obtained using

reﬁned grids.

More sensitive to input-data inaccuracies are calculation data obtained at different

levels of boundary-condition inaccuracies (6.150). These inaccuracies were modeled

in the ordinary way, for instance, as

˜μ

) = μ

) + 2δ(σ(x

) − 1/2), x

∈ ω

where σ(x

) is a random function normally distributed over the interval [0, 1], and the

parameter δ deﬁnes the inaccuracy level. Figure 6.21 shows data obtained by solving

the inverse problem with δ = 0.0003. This inaccuracy level corresponds to a relative

inaccuracy of 0.1%. The calculation grid N

= N

= 129 was used.

Figure 6.18 Solution of the problem obtained on the grid N

= N

= 65