Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics
Подождите немного. Документ загружается.
206 Chapter 6 Right-hand side identification
w
0
= 0,w
N
0
= 0, x ∈ ω, (6.131)
in which
ψ =
1
T
Aϕ.
For two-dimensional mesh functions, we define a Hilbert space H = L
2
(Q
T
),in
which the scalar product and the norm are
(v, w)
∗
=
N
0
−1
n=1
(v
n
,w
n
)τ, v
∗
=
(v, v)
∗
.
We scalarwise multiply equation (6.130) in H = L
2
(Q
T
) by w; this yields:
w
¯
tt
,w
∗
+
Aw
◦
t
,w
∗
=−(ψ, w)
∗
. (6.132)
Here, the following standard settings adopted in the theory of difference schemes are
used:
w
¯
t
=
w
n
− w
n−1
τ
,w
t
=
w
n+1
− w
n
τ
,
w
¯
tt
=
w
n+1
− 2w
n
+ w
n−1
τ
2
,w
◦
t
=
w
n+1
− w
n−1
2τ
.
Further consideration is based on the directly verifiable property of skew symmetry
of the operator of central difference derivative, so that
Aw
◦
t
,w
∗
= 0.
Analogously to (6.121), from (6.132) we obtain
N
0
n=1
(w
¯
t
)
2
τ,1
= (ψ, w)
∗
. (6.133)
With the difference Friedrichs inequality taken into account, we obtain
N
0
−1
n=1
w
2
τ ≤ M
0
N
0
n=1
(w
¯
t
)
2
τ,
where
M
0
= T
2
/8.
In view of
(ψ, w)
∗
≤ψ
∗
w
∗
,
from (6.133) we obtain the estimate
w
∗
≤ M
0
ψ
∗
. (6.134)
Section 6.4 Identification of time-independent right-hand side: parabolic equations 207
Figure 6.12 Mesh pattern of the difference problem
Then, representation (6.129) and inequality (6.134) yield the desired estimate
u
∗
≤ϕ
∗
+
M
0
T
Aϕ
∗
. (6.135)
Estimate (6.135) is fully consistent with the a priori estimate (6.123) for the solution
of the differential problem.
The right-hand side can be found from the solution of problem (6.127), (6.128) at
half-integer nodes using the expression
f
n+1/2
=
u
n+1
− u
n
τ
+ A
u
n+1
+ u
n
2
, (6.136)
matched with the difference approximation (6.127) of equation (6.114).
6.4.4 Solution of the difference problem
The a priori estimate (6.135) obtained previously in the consideration of the problem
for the inaccuracy can be used to establish convergence of the difference solution found
from (6.127), (6.128) to the solution of the boundary value problem (6.112)–(6.114)
with second order over time and space. Some problems arise in the computational
realization and in finding the solution of the difference problem (6.127), (6.128). For
problems with the constant coefficient k(x), fast algorithms can be constructed around
the variable separation method. For more general problems (in particular, for problems
with k(x) = const), iteration methods are preferable. Figure 6.12 shows the mesh
pattern for the difference problem (6.127), (6.128).
Consider a possible approach to solving the difference problem in question. With
(6.129), the problem reduces to finding the mesh function w as the solution of (6.130),
(6.131). We write this difference problem as
¯
Aw = ψ,
¯
A = A
0
+ A
1
. (6.137)
208 Chapter 6 Right-hand side identification
The operators A
0
and A
1
are defined by
A
0
w
n
=−
w
n+1
− 2w
n
+ w
n−1
τ
2
, (6.138)
A
1
w
n
=−A
w
n+1
− w
n−1
2τ
, n = 1, 2,...,N
0
− 1 (6.139)
on the set of mesh functions satisfying the conditions (6.131).
A specific feature of the difference problem (6.137) consists in that the operator
¯
A
in this problem is not a self-adjoint operator. This circumstance substantially hampers
the construction of efficient iteration methods. As it was shown above, the operator
A
1
defined by (6.139) is a skew-symmetric operator (A
1
w, w) = 0). Among the main
properties of the operator A
0
(see (6.138), the following property deserves mention:
A
0
= A
∗
0
≥ M
−1
0
E.
Thus, the operators A
0
and A
1
are respectively the self-adjoint and skew-symmetric
parts of
¯
A:
A
0
=
1
2
(
¯
A +
¯
A
∗
), A
1
=
1
2
(
¯
A −
¯
A
∗
).
Among the general approaches to the solution of difference problems with non-
self-adjoint operators, methods using preliminary symmetrization are worth noting. In
this way, one can pass from an initial problem with a non-self-adjoint operator to the
problem with a self-adjoint operator. An example here is the Gauss symmetrization,
in which case instead of (6.137) to be solved is the equation
¯
A
∗
¯
Aw =
¯
A
∗
ψ. (6.140)
A second idea widely used in computational practice, in the approximate solution
of problems with non-self-adjoint operator (bad problem) is related with the choice,
as the iteration-method reconditioner, the self-conjugate part of the problem operator
(B = A
0
) (good problem at each iteration step). In the solution of problem (6.137),
we will use an iteration method developed around the mentioned passage to a problem
with self-adjoint operator, with the operator A
0
used as the reconditioner.
At the first stage, the initial problem is to be reconditioned as
A
−1
0
(A
0
+ A
1
)w = A
−1
0
ψ.
At the second stage, the operator A
0
+ A
∗
1
, A
∗
1
=−A
1
(compare (6.140)) is used to
perform symmetrization:
˜
Aw =
˜
f , (6.141)
˜
A = ( A
0
+ A
∗
1
)A
−1
0
(A
0
+ A
1
),
˜
f = (A
0
+ A
∗
1
)A
−1
0
ψ.
Section 6.4 Identification of time-independent right-hand side: parabolic equations 209
The problem (6.141) can be solved using the iterative conjugate method that takes,
with the so-chosen reconditioner operator (B = A
0
), the form
A
0
w
k+1
= α
k+1
(A
0
− τ
k+1
˜
A)w
k
+ (1 − α
k+1
)A
0
w
k−1
+ α
k+1
τ
k+1
˜
f ,
k = 1, 2,..., (6.142)
A
0
w
1
= (A
0
− τ
1
˜
A)w
0
+ τ
1
˜
f .
The iteration parameters α
k+1
and τ
k+1
can be calculated by the formulas
τ
k+1
=
( ˜w
k
, r
k
)
(
˜
A ˜w
k
, ˜w
k
)
, k = 0, 1,...,
α
k+1
=
1 −
τ
k+1
τ
k
( ˜w
k
, r
k
)
( ˜w
k−1
, r
k−1
)
1
α
k
−1
, k = 1, 2 ..., α
1
= 1,
(6.143)
where r
k
=
˜
Aw
k
−
˜
f is the discrepancy and ˜w = A
−1
0
r
k
is the correction.
The iteration method (6.142), (6.143) is embodied in the program PROBLEM8.
Program PROBLEM8
C
C PROBLEM8 - IDENTIFICATION OF TIME-INDEPENDENT RIGHT-HAND SIDE
C ONE-DIMENSIONAL NON-STATIONARY PROBLEM
C (UNKNOWN SPATIAL DISTRIBUTION)
C
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
PARAMETER ( DELTA = 0.D0,N=51,M=51)
DIMENSION X(N), Y(N), FT(N), FY(N), PHI(N), PHID(N)
+ ,A(M), B(M), C(M), F(M) ! M .GE. N
+ ,V(N,M), PSI(N), FP(N,M), FR(N,M), WORK(N,M)
+ ,W(N,M), WOLD(N,M), RK(N,M), ARK(N,M), FTILDE(N,M)
C
C PARAMETERS:
C
C XL, XR - LEFT AND RIGHT END POINTS OF THE SEGMENT;
C N - NUMBER OF GRID NODES OVER THE SPATIAL VARIABLE;
C TMAX - MAXIMAL TIME;
C M - NUMBER OF GRID NODES OVER TIME;
C DELTA - INPUT-DATA INACCURACY;
C PHI(N) - EXACT SOLUTION AT THE END TIME;
C PHID(N) - DISTURBED END-TIME SOLUTION;
C FT(N) - EXACT SPATIAL DEPENDENCE OF THE RIGHT-HAND SIDE;
C FY(N)) - CALCULATED SPATIAL DEPENDENCE OF THE RIGHT-HAND SIDE;
C EPS - REQUIRED RELATIVE ACCURACY IN THE ITERATIVE APPROACH
C TO THE SOLUTION.
C
XL = 0.D0
XR = 1.D0
TMAX = 1.D0
C
OPEN (01, FILE = ’RESULT.DAT’) ! FILE TO STORE THE CALCULATED DATA
C
210 Chapter 6 Right-hand side identification
C GRID
C
H =(XR-XL)/(N-1)
TAU=TMAX/(M-1)
C
C DIRECT PROBLEM
C
C EXACT RIGHT-HAND SIDE
C
DOI=1,N
X(I) = (I-1)
*
H
FT(I) = AF(X(I))
END DO
C
C SOLUTION OF THE DIRECT PROBLEM
C
C INITIAL CONDITION
C
T = 0.D0
DOI=1,N
Y(I) = 0.D0
END DO
C
C NEXT TIME LAYER
C
DOK=2,M
T=T+TAU
C
C DIFFERENCE-SCHEME COEFFICIENTS
C
DOI=2,N-1
X1 = (X(I) + X(I-1)) / 2
X2 = (X(I+1) + X(I)) / 2
A(I) = AK(X1) / (2.D0
*
H
*
H)
B(I) = AK(X2) / (2.D0
*
H
*
H)
C(I) = A(I) + B(I) + 1.D0 / TAU
F(I) = A(I)
*
Y(I-1)
+ + (1.D0 / TAU - A(I) - B(I))
*
Y(I)
+ + B(I)
*
Y(I+1)
+ + AF(X(I))
END DO
C
C BOUNDARY CONDITIONS AT THE LEFT AND RIGHT END POINTS
C
B(1) = 0.D0
C(1) = 1.D0
F(1) = 0.D0
A(N) = 0.D0
C(N) = 1.D0
F(N) = 0.D0
C
C SOLUTION OF THE PROBLEM AT THE NEXT TIME LAYER
C
ITASK = 1
CALL PROG3 ( N, A, C, B, F, Y, ITASK )
END DO
C
C SOLUTION AT THE END TIME
Section 6.4 Identification of time-independent right-hand side: parabolic equations 211
C
DOI=1,N
PHI(I) = Y(I)
END DO
C
C DISTURBING OF MEASURED QUANTITIES
C
DOI=2,N-1
PHID(I) = PHI(I) + 2.D0
*
DELTA
*
(RAND(0)-0.5D0)
END DO
PHID(1) = 0.D0
PHID(N) = 0.D0
C
C INVERSE PROBLEM
C
EPS = 1.D-8
C
C AUXILIARY FUNCTION AND INITIALIZATION
C
DOI=1,N
DOK=1,M
V(I,K) = (K-1)
*
PHID(I) / (M-1)
FP(I,K) = 0.D0
FR(I,K) = 0.D0
RK(I,K) = 0.D0
ARK(I,K) = 0.D0
WORK(I,K) = 0.D0
FTILDE(I,K) = 0.D0
END DO
END DO
C
C RIGHT-HAND SIDE OF THE SYMMETRIZED EQUATION
C
DOI=2,N-1
X1 = (X(I) + X(I-1)) / 2
X2 = (X(I+1) + X(I)) / 2
A(I) = AK(X1) / (H
*
H)
B(I) = AK(X2) / (H
*
H)
PSI(I) = - A(I)
*
PHID(I-1) + (A(I) + B(I))
*
PHID(I)
+ - B(I)
*
PHID(I+1)
PSI(I) = PSI(I) / TMAX
END DO
DOI=2,N-1
DOK=2,M-1
FP(I,K) = PSI(I)
END DO
END DO
CALL A0I ( N, M, TAU, FR, FP, A, B, C, F, Y )
DOI=2,N-1
DOK=2,M-1
X1 = (X(I) + X(I-1)) / 2
X2 = (X(I+1) + X(I)) / 2
AA = AK(X1) / (H
*
H)
BB = AK(X2) / (H
*
H)
WORK(I,K)=-AA
*
(FR(I-1,K+1) - FR(I-1,K-1))
+ + (AA + BB)
*
(FR(I,K+1) - FR(I,K-1))
+-BB
*
(FR(I+1,K+1) - FR(I+1,K-1))
WORK(I,K) = WORK(I,K) / (2.D0
*
TAU)
212 Chapter 6 Right-hand side identification
END DO
END DO
DOI=2,N-1
DOK=2,M-1
FTILDE(I,K) = FP(I,K) + WORK(I,K)
END DO
END DO
C
C ITERATIVE CONJUGATE-GRADIENT METHOD
C
NIT = 0
NITMAX = 5000
AL = 1.D0
C
C INITIAL APPROXIMATION
C
DOI=1,N
DOK=1,M
W(I,K) = 0.D0
WOLD(I,K) = 0.D0
END DO
END DO
C
C ITERATION CYCLE
C
100 NIT = NIT + 1
C
C DISCREPANCY
C
CALL ATILDE ( N, M, H, TAU, W, RK, WORK, A, B, C, F, Y )
DOI=2,N-1
DOK=2,M-1
RK(I,K) = RK(I,K) - FTILDE(I,K)
END DO
END DO
C
C CORRECTION
C
CALL A0I ( N, M, TAU, ARK, RK, A, B, C, F, Y )
C
C ITERATION PARAMETERS
C
CALL ATILDE ( N, M, H, TAU, ARK, FR, WORK, A, B, C, F, Y )
RR = 0.D0
RA = 0.D0
DOI=2,N-1
DOK=2,M-1
RR = RR + RK(I,K)
*
ARK(I,K)
RA = RA + FR(I,K)
*
ARK(I,K)
END DO
END DO
IF (NIT.EQ.1) RR0 = RR
Section 6.4 Identification of time-independent right-hand side: parabolic equations 213
TAUK=RR/RA
IF (NIT.GT.1)
+ AL = 1.D0 / (1.D0 - TAUK
*
RR / (TAUKOLD
*
RROLD
*
ALOLD))
C
C NEXT APPROXIMATION
C
DOI=2,N-1
DOK=2,M-1
AA=AL
*
W(I,K) + (1.D0 - AL)
*
WOLD(I,K)
+ - TAUK
*
AL
*
ARK(I,K)
WOLD(I,K) = W(I,K)
W(I,K) = AA
END DO
END DO
RROLD = RR
TAUKOLD = TAUK
ALOLD = AL
C
C END OF ITERATIONS
C
IF (RR .GE. EPS
*
RR0 .AND. NIT .LT. NITMAX) GO TO 100
C
C APPROXIMATE SOLUTION OF THE INVERSE PROBLEM
C
DOI=1,N
DOK=1,M
W(I,K) = W(I,K) + V(I,K)
END DO
END DO
DOI=2,N-1
X1 = (X(I) + X(I-1)) / 2
X2 = (X(I+1) + X(I)) / 2
A(I) = AK(X1) / (2.D0
*
H
*
H)
B(I) = AK(X2) / (2.D0
*
H
*
H)
FY(I) = (W(I,2) - W(I,1)) / TAU
+ - A(I)
*
(W(I-1,2) + W(I-1,1))
+ + (A(I)+B(I))
*
(W(I,2) + W(I,1))
+ - B(I)
*
(W(I+1,2) + W(I+1,1))
END DO
C
WRITE ( 01,
*
) NIT
WRITE ( 01,
*
) (FT(I), I = 2,N-1)
WRITE ( 01,
*
) (FY(I), I = 2,N-1)
WRITE ( 01,
*
) (PHI(I), I = 1,N)
WRITE ( 01,
*
) (PHID(I), I = 1,N)
CLOSE (01)
STOP
END
C
DOUBLE PRECISION FUNCTION AK(X)
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C COEFFICIENT AT THE HIGHER DERIVATIVES
C
AK = 1.D-1
C
RETURN
END
C
214 Chapter 6 Right-hand side identification
DOUBLE PRECISION FUNCTION AF(X)
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C RIGHT-HAND SIDE
C
AF = 1.D0
IF (X.GT.0.5D0) AF = 0.D0
C
RETURN
END
C
SUBROUTINE A0I ( N, M, TAU, V, F, A, B, C, FF, YY )
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C SOLUTION OF THE DIFFERENCE PROBLEM A0V=F
C
DIMENSION V(N,M), F(N,M)
+ ,A(M), C(M), B(M), FF(M), YY(M)
C
DOI=2,N-1
C
C DIFFERENCE-SCHEME COEFFICIENTS
C
DOK=2,M-1
A(K) = 1.D0 / (TAU
*
TAU)
B(K) = A(K)
C(K) = A(K) + B(K)
FF(K) = F(I,K)
END DO
C
C BOUNDARY CONDITIONS AT THE LEFT AND RIGHT END POINTS
C
B(1) = 0.D0
C(1) = 1.D0
FF(1) = 0.D0
A(M) = 0.D0
C(M) = 1.D0
FF(M) = 0.D0
C
C SOLUTION OF THE PROBLEM
C
ITASK = 1
CALL PROG3 ( M, A, C, B, FF, YY, ITASK )
DOK=2,M-1
V(I,K) = YY(K)
END DO
END DO
RETURN
END
C
SUBROUTINE ATILDE ( N, M, H, TAU, V, F, WORK, A, B, C, FF, YY )
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C CALCULATION OF F = ATILDE V
C
Section 6.4 Identification of time-independent right-hand side: parabolic equations 215
DIMENSION V(N,M), F(N,M), WORK(N,M)
+ ,A(M), C(M), B(M), FF(M), YY(M)
C
DOI=2,N-1
DOK=2,M-1
X1 = (I-0.5D0)
*
H
X2 = (I+0.5D0)
*
H
AA = AK(X1) / (H
*
H)
BB = AK(X2) / (H
*
H)
F(I,K)=-AA
*
(V(I-1,K+1) - V(I-1,K-1))
+ + (AA + BB)
*
(V(I,K+1) - V(I,K-1))
+-BB
*
(V(I+1,K+1) - V(I+1,K-1))
F(I,K) = F(I,K) / (2.D0
*
TAU)
END DO
END DO
CALL A0I ( N, M, TAU, WORK, F, A, B, C, FF, YY )
C
DOI=2,N-1
DOK=2,M-1
X1 = (I-0.5D0)
*
H
X2 = (I+0.5D0)
*
H
AA = AK(X1) / (H
*
H)
BB = AK(X2) / (H
*
H)
F(I,K)=-AA
*
(WORK(I-1,K+1) - WORK(I-1,K-1))
+ + (AA + BB)
*
(WORK(I,K+1) - WORK(I,K-1))
+-BB
*
(WORK(I+1,K+1) - WORK(I+1,K-1))
F(I,K) = F(I,K) / (2.D0
*
TAU)
END DO
END DO
C
AA = 1.D0 / (TAU
*
TAU)
DOI=2,N-1
DOK=2,M-1
F(I,K)=-AA
*
(V(I,K-1) - 2.D0
*
V(I,K) + V(I,K+1))
+ - F(I,K)
END DO
END DO
RETURN
END
The subroutine A0I solves the equation A
0
v = f , and the subroutine ATILDE
calculates the value of f =
˜
Av from the given v.
6.4.5 Computational experiments
The presented program PROBLEM8 solves the inverse problem (6.110)–(6.113) in the
case of
k(x) = 0.1, f (x) =
1, 0 < x < 0.5,
0, 0.5 ≤ x < 1,
l = 1, T = 1,