Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics
Подождите немного. Документ загружается.
166 Chapter 6 Right-hand side identification
Substituting (6.38) with i = 1 into (6.32), we obtain
(C
1
− B
1
α
1
)y
1
= (D
1
− B
1
β
1
)y
2
− E
1
y
3
+ F
1
+ B
1
γ
1
.
Hence, we can set:
α
2
= S
−1
1
(D
1
− B
1
β
1
), S
1
= C
1
− B
1
α
1
,
β
2
= S
−1
1
E
1
,γ
2
= S
−1
1
(F
1
+ B
1
γ
1
).
(6.43)
To determine the coefficients α
N
,γ
N
in (6.37), we substitute (6.38), (6.39) with
i = N − 1 into (6.34). This yields the same formulas (6.40), (6.41) with i = N − 1.
Thus, all sweep coefficients can be found by the recurrent formulas (6.40), (6.41)
together with (6.42), (6.43).
For the solution to be found from (6.36), (6.37), we have to find y
N
. To this end, we
use equation (6.35). Substitution of (6.36) with i = N − 2 yields
y
N
= γ
N +1
, (6.44)
where γ
N +1
is to be found using formulas (6.40), (6.41) with i = N . Equality (6.44)
allows the solution of (6.31)–(6.35) to be found from (6.36), (6.37).
The above algorithm is realized in the subroutine PROG5:
Subroutine PROG5
SUBROUTINE PROG5 ( N, A, B, C, D, E, F, Y, ITASK )
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C SWEEP METHOD
C FOR FIVE-DIAGONAL MATRIX
C
C ITASK = 1: FACTORIZATION AND SOLUTION;
C
C ITASK = 2: SOLUTION ONLY
C
DIMENSION A(N), B(N), C(N), D(N), E(N), F(N), Y(N)
IF ( ITASK .EQ. 1 ) THEN
C
D(1) = D(1) / C(1)
E(1) = E(1) / C(1)
C(2) = C(2) - D(1)
*
B(2)
D(2) = ( D(2) - E(1)
*
B(2) ) / C(2)
E(2) = E(2) / C(2)
C
DOI=3,N
C(I) = C(I) - E(I-2)
*
A(I) + D(I-1)
*
( D(I-2)
*
A(I) - B(I) )
D(I) = ( D(I) + E(I-1)
*
( D(I-2)
*
A(I)-B(I)))/C(I)
E(I) = E(I) / C(I)
END DO
C
ITASK = 2
Section 6.1 Right-hand side reconstruction from known solution: stationary problems 167
END IF
C
F(1) = F(1) / C(1)
F(2) = ( F(2) + F(1)
*
B(2) ) / C(2)
DOI=3,N
F(I) = ( F(I) - F(I-2)
*
A(I)
*
- F(I-1)
*
( D(I-2)
*
A(I)-B(I)))/C(I)
END DO
C
Y(N) = F(N)
Y(N-1) = D(N-1)
*
Y(N) + F(N-1)
DO I = N-2, 1, -1
Y(I) = D(I)
*
Y(I+1) - E(I)
*
Y(I+2) + F(I)
END DO
RETURN
END
The identification problem was solved within the framework of a quasi-real experi-
ment. First, the direct problem (6.1), (6.2) with given coefficients and right-hand side
was solved. Afterwards, random perturbations were introduced into the mesh solution
and, then, the inverse problem was treated. The difference solution of the direct prob-
lem was found using the sweep method for tridiagonal matrices. Like in the case of
PROG5, in the subroutine PROG3 the matrix coefficients are not retained at the exit,
being replaced instead with the sweep coefficients.
Subroutine PROG3
SUBROUTINE PROG3 ( N, A, C, B, F, Y, ITASK )
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C SWEEP METHOD
C FOR TRIDIAGONAL MATRIX
C
C ITASK = 1: FACTORIZATION AND SOLUTION;
C
C ITASK = 2: SOLUTION ONLY
C
DIMENSION A(N), C(N), B(N), F(N), Y(N)
IF ( ITASK .EQ. 1 ) THEN
C
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
C
ITASK = 2
END IF
C
F(1) = F(1) / C(1)
DOI=2,N
168 Chapter 6 Right-hand side identification
F(I) = ( F(I) + F(I-1)
*
A(I) ) / C(I)
END DO
C
Y(N) = F(N)
DO I = N-1, 1, -1
Y(I) = B(I)
*
Y(I+1) + F(I)
END DO
RETURN
END
The regularization parameter was found from the discrepancy using the sequence
α
k
= α
0
q
k
, q > 0
with set α
0
and q.
Program PROBLEM5
C
C PROBLEM5 - RIGHT-HAND SIDE IDENTIFICATION
C ONE-DIMENSIONAL STATIONARY PROBLEM
C
IMPLICIT REAL
*
8 ( A-H, O-Z )
PARAMETER ( DELTA = 0.001D0, N = 201 )
DIMENSION Y(N), X(N), YD(N), F(N), FA(N), FT(N)
+ ,A(N), B(N), C(N), D(N), E(N), FF(N)
C
C PARAMETERS:
C
C XL, XR - LEFT AND RIGHT END POINTS OF THE SEGMENT;
C N - NUMBER OF GRID NODES;
C DELTA - INPUT-DATA INACCURACY;
C ALPHA - INITIAL REGULARIZATION PARAMETER;
C Q - MULTIPLIER IN THE REGULARIZATION PARAMETER;
C Y(N) - EXACT DIFFERENCE SOLUTION OF THE BOUNDARY-VALUE
C PROBLEM;
C YD(N) - DISTURBED DIFFERENCE SOLUTION OF THE BOUNDARY-VALUE
C PROBLEM;
C F(N) - EXACT RIGHT-HAND SIDE;
C FA(N) - CALCULATED RIGHT-HAND SIDE;
C
XL=0.
XR=1.
C
OPEN (01, FILE = ’RESULT.DAT’) ! FILE TO STORE THE CALCULATED DATA
C
C GRID
C
H=(XR-XL)/(N-1)
DOI=1,N
X(I) = XL + (I-1)
*
H
END DO
C
C DIFFERENCE-SCHEME COEFFICIENTS IN THE DIRECT PROBLEM
C
Section 6.1 Right-hand side reconstruction from known solution: stationary problems 169
B(1) = 0.D0
C(1) = 1.D0
FF(1) = 0.D0
A(N) = 0.D0
C(N) = 1.D0
FF(N) = 0.D0
DO I = 2,N-1
A(I) = AK(X(I)-0.5D0
*
H)/(H
*
H)
B(I) = AK(X(I)+0.5D0
*
H)/(H
*
H)
C(I) = A(I) + B(I) + AQ(X(I))
FF(I) = AF(X(I))
END DO
C
C SOLUTION OF THE DIFFERENCE PROBLEM
C
ITASK1 = 1
CALL PROG3 ( N, A, C, B, FF, Y, ITASK1 )
C
C NOISE ADDITION TO THE SOLUTION OF THE BOUNDARY-VALUE PROBLEM
C
YD(1) = 0.D0
YD(N) = 0.D0
DO I = 2,N-1
YD(I) = Y(I) + 2.
*
DELTA
*
(RAND(0)-0.5)
END DO
C
C RIGHT-HAND SIDE
C
FT(1) = 0.D0
FT(N) = 0.D0
DO I = 2,N-1
FT(I) = - AK(X(I)-0.5D0
*
H)
*
YD(I-1) / (H
*
H) +
+ ((AK(X(I)-0.5D0
*
H) + AK(X(I)+0.5D0
*
H))/(H
*
H) +
+ AQ(X(I)))
*
YD(I) -
+ AK(X(I)+0.5D0
*
H)
*
YD(I+1) / (H
*
H)
END DO
WRITE ( 01,
*
) (X(I), I = 1,N)
WRITE ( 01,
*
) (Y(I), I = 1,N)
WRITE ( 01,
*
) (YD(I), I = 1,N)
WRITE ( 01,
*
) (FT(I), I=1,N)
C
C ITERATIVE ADJUSTMENT OF THE REGULARIZATION PARAMETER
C
IT=0
ITMAX = 1000
ALPHA = 0.001D0
Q = 0.75D0
100IT=IT+1
C
C DIFFERENCE-SCHEME COEFFICIENTS IN THE INVERSE PROBLEM
C
C(1) = 1.D0
D(1) = 0.D0
E(1) = 0.D0
170 Chapter 6 Right-hand side identification
B(2) = 0.D0
C(2) = ALPHA
*
(AK(X(2)+0.5D0
*
H)/(H
**
2))
**
2+
+ ALPHA
*
( (AK(X(2)-0.5D0
*
H) +
+ AK(X(2)+0.5D0))/(H
**
2) + AQ(X(2)))
**
2 + 1.D0
D(2) = ALPHA
*
AK(X(2)+0.5D0
*
H)/(H
**
2)
*
+ ( (AK(X(2)+0.5D0
*
H) +
+ AK(X(2)+1.5D0
*
H))/(H
**
2) + AQ(X(3)) +
+ (AK(X(2)-0.5D0
*
H) +
+ AK(X(2)+0.5D0
*
H))/(H
**
2) + AQ(X(2)))
E(2) = ALPHA
*
AK(X(2)+0.5D0
*
H)
*
AK(X(2)+1.5D0
*
H)/(H
**
4)
A(N-1) = ALPHA
*
AK(X(N-1)-0.5D0
*
H)
*
AK(X(N-1)-1.5D0
*
H)/(H
**
4)
B(N-1) = ALPHA
*
AK(X(N-1)-0.5D0
*
H)/(H
**
2)
*
+ ( (AK(X(N-1)-0.5D0
*
H) + AK(X(N-1)-1.5D0
*
H))/(H
**
2) +
+ AQ(X(N-2) ) +
+ (AK(X(N-1)-0.5D0
*
H) + AK(X(N-1)+0.5D0
*
H))/(H
**
2) +
+ AQ(X(N-1)))
C(N-1) = ALPHA
*
(AK(X(N-1)-0.5D0
*
H)/(H
**
2))
**
2+
+ ALPHA
*
( (AK(X(N-1)-0.5D0
*
H) +
+ AK(X(N-1)+0.5D0
*
H))/(H
**
2) + AQ(X(N-1)))
**
2 + 1.D0
D(N-1) = 0.D0
A(N) = 0.D0
B(N) = 0.D0
C(N) = 1.D0
DO I = 2,N-2
A(I) = ALPHA
*
AK(X(I)-0.5D0
*
H)
*
AK(X(I)-1.5D0
*
H)/(H
**
4)
B(I) = ALPHA
*
AK(X(I)-0.5D0
*
H)/(H
**
2)
*
+ ( (AK(X(I)-0.5D0
*
H) + AK(X(I)-1.5D0
*
H))/(H
**
2) +
+ AQ(X(I-1)) +
+ (AK(X(I)-0.5D0
*
H) + AK(X(I)+0.5D0
*
H))/(H
**
2) + AQ(X(I)))
C(I) = ALPHA
*
(AK(X(I)-0.5D0
*
H)/(H
**
2))
**
2+
+ ALPHA
*
(AK(X(I)+0.5D0
*
H)/(H
**
2))
**
2+
+ ALPHA
*
( (AK(X(I)-0.5D0
*
H) +
+ AK(X(I)+0.5D0))/(H
**
2) + AQ(X(I)))
**
2 + 1.D0
D(I) = ALPHA
*
AK(X(I)+0.5D0
*
H)/(H
**
2)
*
+ ( (AK(X(I)+0.5D0
*
H) + AK(X(I)+1.5D0
*
H))/(H
**
2) +
+ AQ(X(I+1)) +
+ (AK(X(I)-0.5D0
*
H) + AK(X(I)+0.5D0
*
H))/(H
**
2) + AQ(X(I)))
E(I) = ALPHA
*
AK(X(I)+0.5D0
*
H)
*
AK(X(I)+1.5D0
*
H)/(H
**
4)
END DO
C
C SOLUTION OF THE DIFFERENCE PROBLEM
C
ITASK2 = 1
DO I = 1,N
FF(I) = FT(I)
END DO
CALL PROG5 ( N, A, B, C, D, E, FF, FA, ITASK2 )
C
C DISCREPANCY
C
B(1) = 0.D0
C(1) = 1.D0
FF(1) = 0.D0
A(N) = 0.D0
C(N) = 1.D0
FF(N) = 0.D0
Section 6.1 Right-hand side reconstruction from known solution: stationary problems 171
DO I = 2,N-1
A(I) = AK(X(I)-0.5D0
*
H)/(H
*
H)
B(I) = AK(X(I)+0.5D0
*
H)/(H
*
H)
C(I) = A(I) + B(I) + AQ(X(I))
FF(I) = FA(I)
END DO
C
ITASK1 = 1
CALL PROG3 ( N, A, C, B, FF, D, ITASK1 )
C
SUM = 0.D0
DO I = 2,N-1
SUM = SUM + (D(I)-YD(I))
**
2
*
H
END DO
SL2 = DSQRT(SUM)
C
IF ( IT.EQ.1 ) THEN
IND=0
IF ( SL2.LT.DELTA ) THEN
IND=1
Q = 1.D0/Q
END IF
ALPHA = ALPHA
*
Q
GO TO 100
ELSE
ALPHA = ALPHA
*
Q
IF ( IND.EQ.0 .AND. SL2.GT.DELTA ) GO TO 100
IF ( IND.EQ.1 .AND. SL2.LT.DELTA ) GO TO 100
END IF
C
C SOLUTION
C
WRITE ( 01,
*
) (FA(I), I=1,N)
WRITE ( 01,
*
) IT, ALPHA, SL2
CLOSE ( 01 )
C
STOP
END
DOUBLE PRECISION FUNCTION AK(X)
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C COEFFICIENT AT THE HIGHER DERIVATIVES
C
AK = 0.05D0
C
RETURN
END
DOUBLE PRECISION FUNCTION AQ(X)
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C COEFFICIENT AT THE LOWEST TERM
C
AQ = 0.D0
C
RETURN
172 Chapter 6 Right-hand side identification
END
DOUBLE PRECISION FUNCTION AF(X)
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C RIGHT-HAND SIDE
C
AF=X
IF ( X.GE.0.5D0 ) AF = 0.5D0
C
RETURN
END
The equation coefficients and the right-hand side are set in the subroutine functions
AK, AQ, AF.
6.1.6 Examples
As an example, the model problem (6.1), (6.2) with
k(x) = 0.05, q(x) = 0, f (x) =
x, 0 < x < 0.5,
0.5, 0.5 < x < 1,
was considered. The input-data inaccuracies were modeled by perturbing the differ-
ence solution of the problem with exact right-hand side at grid nodes:
y
δ
(x) = y(x) + 2δ(σ(x) − 1/2), x = x
i
, i = 0, 1,...,N − 1,
where σ(x) is a random quantity normally distributed over the interval from 0 to 1.
Note first some direct difference differentiation possibilities in calculating the right-
hand side of (6.1). Figure 6.1 shows data calculated on various difference grids. The
data were obtained for the perturbed solution of the direct problem on a grid with
N = 512; here, the noise level was defined by the parameter δ = 0.001. Acceptable
results could also be obtained using numerical differentiation on rather crude grids
with h ≥ 1/16.
The solution of the same identification problem obtained by the Tikhonov method is
shown in Figure 6.2. Here, a calculation grid with N = 200 was used. The effect due
to the inaccuracy level can be traced considering the data calculated with the input-
data inaccuracy δ = 0.01; these data are shown in Figure 6.3. A substantial loss of
accuracy near the right boundary x = l is observed. This loss is related with the fact
that the calculated operator D is defined on the set of functions satisfying the boundary
conditions (6.2). If the right-hand side to be found also satisfies these conditions, an
improved accuracy in finding the approximate solution can be expected. The latter
is illustrated by Figure 6.4 that shows the solution of the identification problem with
somewhat modified right-hand side.
Section 6.1 Right-hand side reconstruction from known solution: stationary problems 173
Figure 6.1 Numerical differentiation with δ = 0.001
Figure 6.2 Solution of the identification problem obtained with δ = 0.001
174 Chapter 6 Right-hand side identification
Figure 6.3 Solution of the identification problem obtained with δ = 0.01
Figure 6.4 Solution of the identification problem obtained with δ = 0.01 and for a
modified right-hand side
Section 6.2 Right-hand side identification in the case of parabolic equation 175
6.2 Right-hand side identification in the case of a parabolic
equation
Below, we consider an inverse problem in which it is required to reconstruct the right-
hand side of a one-dimensional parabolic equation from the known solution. Specific
features of the problem are discussed, related with the evolutionary character of the
problem and with the possibility to sequentially determine the right-hand side with
increasing time.
6.2.1 Model problem
As a model problem, consider the problem in which it is required to reconstruct the
right-hand side of a one-dimensional parabolic equation. Let us begin the considera-
tion with the statement of the direct problem.
The solution u(x, t) is defined in the rectangle
Q
T
= × [0, T ], ={x | 0 ≤ x ≤ l}, 0 ≤ t ≤ T.
The function u(x, t) satisfies the equation
∂u
∂t
=
∂
∂x
k(x)
∂u
∂x
+ f (x, t), 0 < x < l, 0 < t ≤ T (6.45)
under the standard constraint k(x) ≥ κ>0.
For simplicity, homogeneous boundary and initial conditions are assumed:
u(0, t) = 0, u(l, t) = 0, 0 < t ≤ T, (6.46)
u(x, 0) = 0, 0 ≤ x ≤ l. (6.47)
In the direct problem (6.45)–(6.47), the solution u(x, t) is to be found from the
known coefficient k(x) and from the known right-hand side f (x, t). In the inverse
problem, the unknown quantity is the right-hand side f (x, t) (source power), with
the solution u(x, t) assumed known. The right-hand side can be calculated by the
following explicit formula obtained from (6.45):
f (x, t) =
∂u
∂t
−
∂
∂x
k(x)
∂u
∂x
, 0 < x < l, 0 < t ≤ T. (6.48)
The input data are given with some inaccuracy; this circumstance makes the direct
use of formula (6.48) difficult. Most significant here is the effect resulting from the
solution inaccuracy. Let us know, instead of the exact solution u(x, t) of problem
(6.45)–(6.47), a perturbed solution u
δ
(x, t), and in some norm the parameter δ defines
the inaccuracy level in the solution, i. e.,
u
δ
(x, t) − u(x, t )≤δ. (6.49)