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

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

76 Chapter 3 Boundary value problems for elliptic equations

C SETTING OF MATRIX DIAGONALS, ZERO INITIAL APPROXIMATION,

C AND RIGHT-HAND SIDE OF THE SYSTEM:

IDEFAULT(1) = 0

IREPT = 0

DO1I=1,10

A(I) = 6.D0

A(N+I) = 1.D0

A(2

N+I) = 2.D0

Y(I) = 0.D0

1 CONTINUE

F( 1) = -8.D0

F( 2) = -6.D0

F( 3) = -4.D0

F( 4) = -2.D0

F( 5) = 0.D0

F( 6) = 22.D0

F( 7) = 24.D0

F( 8) = 26.D0

F( 9) = 28.D0

F(10) = 41.D0

C SUBROUTINE CALL FOR SOLVING THE SYSTEM:

EPSR = 1.D-6

EPSA = 1.D-7

CALL SBAND5 ( N, L, A, Y, F, EPSR, EPSA )

C PRINTING THE SOLUTION VECTOR:

WRITE ( 06,

)’ SOLUTION:’

DO2I=1,10

WRITE ( 06,

) Y(I)

2 CONTINUE

END

Below, the full text of the subroutine SBAND5 is given.

Subroutine SBAND5

SUBROUTINE SBAND5 ( N, L, A, Y, F, EPSR, EPSA )

IMPLICIT REAL

8 ( A-H, O-Z )

DIMENSION A(1), Y(N), F(N)

COMMON / SB5 / IDEFAULT,

I1R, I2R,

IFREE

COMMON / CONTROL / IREPT, NITER

IF ( IDEFAULT .EQ. 0 ) THEN

I1R=N

I2R=2

IG = 3

Section 3.4 Program realization and numerical examples 77

ELSE IF ( IDEFAULT .NE. 0 ) THEN

IG = IFREE + L

END IF

IT = IG + N

ITL=IT +L

IQ = ITL + N

IF ( IREPT .EQ. 0 ) THEN

A(I1R+N) = 0.D0

DO9J=1,L

A(I2R+J+N-L) = 0.D0

A(IT+J) = 0.D0

9 CONTINUE

DO10J=1,N

G = A(J) - A(I1R+J-1)

A(ITL+J-1)

- A(I2R+J-L)

A(ITL+J-L)

G=1.D0/G

A(IG+J) = DSQRT(G)

A(ITL+J) = G

( A(I1R+J) + A(I2R+J) )

10 CONTINUE

DO20J=1,N

A(J) = A(IG+J)

A(J)

A(IG+J) - 2.D0

A(I1R+J) = A(IG+J)

A(I1R+J)

A(IG+J+1)

A(I2R+J) = A(IG+J)

A(I2R+J)

A(IG+J+L)

20 CONTINUE

IREPT = 1

END IF

DO30J=N,1,-1

A(ITL+J) = Y(J) / A(IG+J)

Y(J) = A(ITL+J) - A(I1R+J)

A(ITL+J+1) - A(I2R+J)

A(ITL+J+L)

A(ITL+J-L) = 0.D0

30 CONTINUE

RR0 = 0.D0

DO 40 J = 1,N

G = A(ITL+J)

A(ITL+J) = F(J)

A(IG+J) - A(J)

G - Y(J)

+ A(I1R+J-1)

A(ITL+J-1) + A(I2R+J-L)

A(ITL+J-L)

F(J) = A(ITL+J) - G

RR0 = RR0 + F(J)

F(J)

40 CONTINUE

NIT = 0

EPSNIT = 1.D0

RR = RR0

50 CONTINUE

IF ( EPSNIT .GE. EPSR .AND. RR .GE. EPSA ) THEN

IF ( NIT .EQ. 0 ) THEN

78 Chapter 3 Boundary value problems for elliptic equations

RRI = 1.D0 / RR

DO60J=N,1,-1

A(IQ+J) = F(J)

A(ITL+J) = A(IQ+J) + A(I1R+J)

A(ITL+J+1) + A(I2R+J)

A(ITL+J+L)

A(ITL+J-L) = 0.D0

60 CONTINUE

ELSE IF ( NIT .GE. 1 ) THEN

BK = RR

RRI

RRI = 1.D0 / RR

DO61J=N,1,-1

A(IQ+J) = F(J) + BK

A(IQ+J)

A(ITL+J) = A(IQ+J) + A(I1R+J)

A(ITL+J+1) + A(I2R+J)

A(ITL+J+L)

A(ITL+J-L) = 0.D0

61 CONTINUE

END IF

TQ = 0.D0

DO70J=1,N

G = A(ITL+J)

A(ITL+J) = A(IQ+J) + A(J)

G + A(I1R+J-1)

A(ITL+J-1)

+ A(I2R+J-L)

A(ITL+J-L)

A(IT+J) = G + A(ITL+J)

TQ = TQ + A(IT+J)

A(IQ+J)

70 CONTINUE

AK=RR/TQ

RR = 0.D0

DO 80 J = 1,N

Y(J) = Y(J) + AK

A(IQ+J)

F(J) = F(J) - AK

A(IT+J)

RR = RR + F(J)

F(J)

80 CONTINUE

NIT = NIT + 1

EPSNIT = DSQRT( RR/RR0 )

GO TO 50

END IF

DO90J=N,1,-1

A(ITL+J) = Y(J) + A(I1R+J)

A(ITL+J+1) + A(I2R+J)

A(ITL+J+L)

Y(J) = A(IG+J)

A(ITL+J)

90 CONTINUE

NITER = NIT

RETURN

END

Section 3.4 Program realization and numerical examples 79

3.4.3 Program

To approximately solve the boundary value problem (3.60), (3.61), the following pro-

gram was used.

Program PROBLEM2

C PROBLEM2 - DIRICHLET BOUNDARY-VALUE PROBLEM FOR THE ELLIPTIC

C EQUATION WITH VARIABLE COEFFICIENTS IN RECTANGLE

IMPLICIT REAL

8 ( A-H, O-Z )

PARAMETER ( N1 = 101, N2 = 101 )

DIMENSION A(7

N2), Y(N1,N2), F(N1,N2),

BL(N2), BR(N2), BB(N1), BT(N1)

COMMON / SB5 / IDEFAULT(4)

COMMON / CONTROL / IREPT, NITER

C THE ARRAY A MUST BE SUFFICIENTLY LONG TO STORE

C THE COEFFICIENTS OF THE SYMMETRIC MATRIX OF THE

C DIFFERENCE PROBLEM, THE SET PRINCIPAL DIAGONAL A0,

C THE SUPERDIAGONAL A1, THE UPPER REMOTED DIAGONAL A2,

C THE SOLUTION VECTOR Y, THE RIGHT-HAND SIDE VECTOR F,

C AND THE VECTORS USED IN THE ITERATION ALGORITHM

C FOR SOLVING THE DIFFERENCE EQUATION

C ( SEE SUBROUTINE SBAND5 ).

C DATA INPUT:

C X1L, X2L - COORDINATES OF THE LEFT BOTTOM CORNER

C OF THE RECTANGULAR CALCULATION DOMAIN;

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

C X1D, X2D - COORDINATES OF THE RIGHT UPPER CORNER

C OF SUBDOMAIN D2;

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

C DIRECTIONS

C EPSR - DESIRED RELATIVE ACCURACY OF THE ITERATIVE

C APPROXIMATION TO THE SOLUTION.

C EPSA - DESIRED ABSOLUTE ACCURACY OF THE ITERATIVE

C APPROXIMATION TO THE SOLUTION.

X1L = 0.D0

X1R = 1.D0

X2L = 0.D0

X2R = 1.D0

EPSR = 1.D-6

EPSA = 0.D0

N=N1

DOI=1,7

A(I) = 0.D0

END DO

C BOUNDARY CONDITIONS

80 Chapter 3 Boundary value problems for elliptic equations

DOJ=1,N2

BL(J) = 0.D0

BR(J) = 0.D0

END DO

DOI=1,N1

BB(I) = 0.D0

BT(I) = 0.D0

END DO

C TO BE FOUND IN THE SUBROUTINE FDS_EL ARE THE CENTRAL, RIGHT,

C AND UPPER COEFFICIENTS OF THE DIFFERENCE SCHEME

C ON THE FIVE-POINT MESH PATTERN (A0, A1, AND A2, RESPECTIVELY),

C AND THE RIGHT-HAND SIDE F.

CALL FDS_EL ( X1L, X1R, X2L, X2R, N1, N2, BL, BR, BB, BT,

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

N+1), F )

C THE SUBROUTINE SBAND5 SOLVES THE SOLUTION OF THE DIFFERENCE

C PROBLEM BY THE ALTERNATE-TRIANGULAR APPROXIMATE

C FACTORIZATION - CONJUGATE GRADIENT METHOD

IDEFAULT(1) = 0

IREPT = 0

C INITIAL APPROXIMATION

DOI=1,N1

DOJ=1,N2

Y(I,J) = 0.D0

END DO

CALL SBAND5 ( N, N1, A, Y, F, EPSR, EPSA )

OPEN ( 01, FILE = ’RESULT.DAT’ )

WRITE ( 01,

) NITER

WRITE ( 01,

) ((Y(I,J),I=1,N1),J=1,N2)

CLOSE ( 01 )

STOP

END

Among the main components of the program, the subroutine FDS EL, used to set

the difference-problem coefﬁcients, deserves mention.

SUBROUTINE FDS_EL ( X1L, X1R, X2L, X2R, N1, N2, BL, BR, BB, BT,

H1, H2, A0, A1, A2, F )

C SETTING OF DIFFERENCE-SCHEME COEFFICIENTS

Section 3.4 Program realization and numerical examples 81

C FOR SOLVING THE DIRICHLET BOUNDARY-VALUE PROBLEM

C FOR THE SECOND-ORDER ELLIPTIC EQUATION

C WITH VARIABLE COEFFICIENTS IN RECTANGLE

IMPLICIT REAL

8 ( A-H, O-Z )

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

BL(N2), BR(N2), BB(N1), BT(N1)

C PARAMETERS

C INPUT PARAMETERS:

C X1L, X1R - LEFT AND RIGHT END POINTS OF THE SEGMENT (THE FIRST

C VARIABLE);

C X2L, X2R - LEFT AND RIGHT END POINTS OF THE SEGMENT (THE SECOND

C VARIABLE);

C N1, N2 - NUMBER OF NODES OVER THE FIRST AND SECOND VARIABLES;

C BL(N2) - BOUNDARY CONDITION AT THE LEFT BOUNDARY (X1 = X1L);

C BR(N2) - BOUNDARY CONDITION AT THE RIGHT BOUNDARY (X1 = X1RL);

C BB(N1) - BOUNDARY CONDITION AT THE LOWER BOUNDARY (X2 = X2L);

C BT(N1) - BOUNDARY CONDITION AT THE UPPER BOUNDARY (X2 = X2RL).

C OUTPUT PARAMETERS:

C H1, H2 - STEPS OF THE UNIFORM RESTANGULAR GRID;

C A0, A1, A2 - DIFFERENCE-SCHEME COEFFICIENTS

C (ON THE FIVE-POINT MESH PATTERN WITH REGARD FOR

C SYMMETRY);

C F - RIGHT-HAND SIDE OF THE DIFFERENCE SCHEME.

C NOTES:

C THE COEFFICIENTS AND THE RIGHT-HAND SIDE OF THE DIFFERENCE

C EQUATION ARE SET IN THE SUBROUTINES-FUNCTIONS AK, AQ, AF

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

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

H12=H1/H2

H21=H2/H1

HH = H1

C INTERNAL NODES

DOJ=2,N2-1

X2 = X2L + (J-1)

DOI=2,N1-1

X1 = X1L + (I-1)

A1(I-1,J) = H21

AK(X1-0.5D0

H1,X2)

A1(I,J) = H21

AK(X1+0.5D0

H1,X2)

A2(I,J-1) = H12

AK(X1,X2-0.5D0

H2)

A2(I,J) = H12

AK(X1,X2+0.5D0

H2)

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

+HH

AQ(X1,X2)

F(I,J) = HH

AF(X1,X2)

END DO

82 Chapter 3 Boundary value problems for elliptic equations

C LEFT AND RIGHT BOUNDARIES: FIRST-KIND BOUNDARY CONDITION

DOJ=2,N2-1

A1(1,J) = 0.D0

A2(1,J) = 0.D0

A0(1,J) = 1.D0

F(1,J) = BL(J)

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

A2(N1,J) = 0.D0

A0(N1,J) = 1.D0

F(N1,J) = BR(J)

END DO

C BOTTOM AND UPPER BOUNDARY: FIRST-KIND BOUNDARY CONDITION

DOI=2,N1-1

A1(I,1) = 0.D0

A2(I,1) = 0.D0

A0(I,1) = 1.D0

F(I,1) = BB(I)

A1(I,N2) = 0.D0

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

A0(I,N2) = 1.D0

F(I,N2) = BT(I)

END DO

C LEFT BOTTOM CORNER

A1(1,1) = 0.D0

A2(1,1) = 0.D0

A0(1,1) = 1.D0

F(1,1) = 0.5D0

(BL(1) + BB(1))

C LEFT UPPER CORNER

A1(1,N2) = 0.D0

A0(1,N2) = 1.D0

F(1,N2) = 0.5D0

(BL(N2) + BT(1))

C RIGHT BOTTOM CORNER

A2(N1,1) = 0.D0

A0(N1,1) = 1.D0

F(N1,1) = 0.5D0

(BB(N1) + BR(1))

C RIGHT UPPER CORNER

A0(N1,N2) = 1.D0

F(N1,N2) = 0.5D0

(BT(N1) + BR(N2))

RETURN

END

DOUBLE PRECISION FUNCTION AK ( X1, X2 )

IMPLICIT REAL

8 ( A-H, O-Z )

C COEFFICIENTS AT THE HIGHER DERIVATIVES

Section 3.4 Program realization and numerical examples 83

AK = 1.D0

IF ( X1 .LE. 0.5D0.AND. X2 .LE. 0.5D0 ) AK = 10.D0

RETURN

END

DOUBLE PRECISION FUNCTION AQ ( X1, X2 )

IMPLICIT REAL

8 ( A-H, O-Z )

C COEFFICIENTS AT THE LOWER TERMS OF THE EQUATION

AQ = 0.D0

IF ( X1 .LE. 0.5D0.AND. X2 .LE. 0.5D0 ) AK = 0.D0

RETURN

END

DOUBLE PRECISION FUNCTION AF ( X1, X2 )

IMPLICIT REAL

8 ( A-H, O-Z )

C RIGHT-HAND SIDE OF THE EQUATION

AF = 1.D0

IF ( X1 .LE. 0.5D0.AND. X2 .LE. 0.5D0 ) AF = 0.D0

RETURN

END

In the example given below we consider a problem for the elliptic equation with

piecewise-constant coefﬁcients.

3.4.4 Computational experiments

Consider data obtained in experiments on numerical solution of a model bound-

ary value problem. The problem simulates processes in a piecewise-homogeneous

medium, with equation coefﬁcients assumed constant in separate subdomains. In the

calculation domain , a rectangle 

is singled out (see Figure 3.3). The boundary

value problem (3.60), (3.61) is solved in the case in which

k(x), q(x), f (x) =



, q

, f

, x ∈ 

, q

, f

, x ∈ 

The above program listing refers to the basic variant in which

= 1.25, l

= 1, s

= 0.5, s

= 0.5,

= 1,κ

= 10, q

= 0, q

= 0, f

= 1, f

= 0,

and the boundary conditions are assumed homogeneous (μ(x) = 0).

84 Chapter 3 Boundary value problems for elliptic equations

Figure 3.3 Calculation domain

The ﬁrst point of interest concerns the efﬁciency of the computational algorithm

used. We are speaking, ﬁrst of all, about the rate of convergence of the iteration

method, and about the convergence of approximate solution to the exact solution. Cal-

culated data obtained on a sequence of reﬁned (from N = N

= N

= 25 to N = 200)

grids are given in Table 3.1.

N 25 50 100 200 400

n 17 22 32 47 69

max

0.06161 0.06240 0.06266 0.06278 0.06284

Table 3.1 Computations performed on a sequence of consecutively reﬁned grids

Here n is the total number of performed iterations (the program is terminated on reach-

ing the relative accuracy ε = 10

−6

), and y

max

is the maximum magnitude of the ap-

proximate solution.

In the case of the problem with variable coefﬁcients, the dependence n = O(N

1/2

)

holds, typical of steepest iteration methods. Considering the problem with discon-

tinuous coefﬁcients, we cannot bargain for a difference-solution convergence rate of

O(h

), readily achieved in the case of problems with smooth solutions. The actual

rate of convergence can be ﬁgured out to an extent considering data on the maximum

magnitude of the approximate solution.

Figure 3.4 shows the approximate solution of the problem obtained for the basic

variant. It is of interest to trace the effect due to the coefﬁcient ratio κ

/κ

.Weﬁx

= 1 and vary κ

. Two cases of relatively high and low coefﬁcient κ

are illustrated

by Figures 3.5 and 3.6.

The use of high (low) coefﬁcients at the higher-order derivatives corresponds to the

method of ﬁctitious domains. With high coefﬁcients (see Figure 3.5), on the boundary

Section 3.5 Exercises 85

Figure 3.4 Solution obtained on the grid N = N

= N

= 100

Figure 3.5 Solution obtained with κ

= 100

of the irregular domain 

we have homogeneous Dirichlet conditions. The second

limiting case of low κ

gives the homogeneous Neumann approximate boundary con-

ditions (see Figure 3.6).

3.5 Exercises

Exercise 3.1 In an irregular domain , we treat the boundary value problem

−



α=1

∂

∂x



k(x)

∂u

∂x



+ q(x)u = f (x), x ∈ , (3.64)

u(x) = μ(x), x ∈ ∂ (3.65)

with

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