Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics
Подождите немного. Документ загружается.
186 Chapter 6 Right-hand side identification
A(I) = AK(X1) / (2
*
H
*
H)
B(I) = AK(X2) / (2
*
H
*
H)
C(I) = A(I) + B(I) + 1.D0 / TAU
FF(I) = A(I)
*
V(I-1,K-1)
+ + (1.D0 / TAU - A(I) - B(I) )
*
V(I,K-1)
+ + B(I)
*
V(I+1,K-1)
+ + (FA(I,K) + FA(I,K-1)) / 2
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(N) = 0.D0
C(N) = 1.D0
FF(N) = 0.D0
C
C SOLUTION OF THE PROBLEM AT THE NEXT TIME LAYER
C
ITASK = 1
CALL PROG3 ( N, A, C, B, FF, YY, ITASK )
DOI=1,N
V(I,K) = YY(I)
END DO
END DO
C
C SOLUTION OF THE CONJUGATE PROBLEM
C
C INITIAL CONDITION
C
T = TMAX
DOI=1,N
W(I,M) = 0.D0
END DO
C
C NEXT TIME LAYER
C
DO K = M-1, 1, -1
T=T-TAU
C
C DIFFERENCE-SCHEME COEFFICIENTS
C
DOI=2,N
X1 = (X(I) + X(I-1)) / 2
X2 = (X(I+1) + X(I)) / 2
A(I) = AK(X1) / (2
*
H
*
H)
B(I) = AK(X2) / (2
*
H
*
H)
C(I) = A(I) + B(I) + 1.D0 / TAU
FF(I) = A(I)
*
W(I-1,K+1)
+ + (1.D0 / TAU - A(I) - B(I) )
*
W(I,K+1)
+ + B(I)
*
W(I+1,K+1)
+ + (V(I,K) + V(I,K+1)) / 2
END DO
C
C BOUNDARY CONDITIONS AT THE LEFT AND RIGHT END POINTS
C
Section 6.2 Right-hand side identification in the case of parabolic equation 187
B(1) = 0.D0
C(1) = 1.D0
FF(1) = 0.D0
A(N) = 0.D0
C(N) = 1.D0
FF(N) = 0.D0
C
C SOLUTION OF THE PROBLEM AT THE NEXT TIME LAYER
C
ITASK = 1
CALL PROG3 ( N, A, C, B, FF, YY, ITASK )
DOI=1,N
W(I,K) = YY(I)
END DO
END DO
C
C DISCREPANCY
C
DOK=1,M
DOI=1,N
RK(I,K) = W(I,K) - GU(I,K)
END DO
END DO
C
C ITERATION PARAMETER
C
C INITIAL CONDITION
C
T = 0.D0
DOI=1,N
W(I,1) = 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
X1 = (X(I) + X(I-1)) / 2
X2 = (X(I+1) + X(I)) / 2
A(I) = AK(X1) / (2
*
H
*
H)
B(I) = AK(X2) / (2
*
H
*
H)
C(I) = A(I) + B(I) + 1.D0 / TAU
FF(I) = A(I)
*
W(I-1,K-1)
+ + (1.D0 / TAU - A(I) - B(I) )
*
W(I,K-1)
+ + B(I)
*
W(I+1,K-1)
+ + (RK(I,K) + RK(I,K-1)) / 2
END DO
C
C BOUNDARY CONDITIONS AT THE LEFT AND RIGHT END POINTS
C
188 Chapter 6 Right-hand side identification
B(1) = 0.D0
C(1) = 1.D0
FF(1) = 0.D0
A(N) = 0.D0
C(N) = 1.D0
FF(N) = 0.D0
C
C SOLUTION OF THE PROBLEM AT THE NEXT TIME LAYER
C
ITASK = 1
CALL PROG3 ( N, A, C, B, FF, YY, ITASK )
DOI=1,N
W(I,K) = YY(I)
END DO
END DO
C
C QUICKEST DESCEND METHOD
C
SUM1 = 0.D0
SUM2 = 0.D0
DOK=1,M
DOI=1,N
SUM1 = SUM1 + RK(I,K)
*
RK(I,K)
SUM2 = SUM2 + W(I,K)
*
W(I,K)
END DO
END DO
TAUK = SUM1/SUM2
C
C NEXT APPROXIMATION
C
DOK=1,M
DOI=1,N
FA(I,K) = FA(I,K) - TAUK
*
RK(I,K)
END DO
END DO
C
C EXIT FROM THE ITERATION CYCLE
C
SUM = 0.D0
DOK=1,M
DOI=1,N
SUM = SUM + (V(I,K) - YD(I,K))
**
2
END DO
END DO
SL2 = DSQRT(SUM
*
TAU
*
H)
WRITE ( 01,
*
) IT, TAUK, SL2
IF (SL2.GT.DELTA .AND. IT.LT.ITMAX) GO TO 100
C
C EXACT SOLUTION
C
SUM = 0.D0
DOK=1,M
T = (K-1)
*
TAU
DOI=1,N
F(I,K) = AF(X(I),T)
SUM = SUM + (FA(I,K) - F(I,K))
**
2
Section 6.2 Right-hand side identification in the case of parabolic equation 189
END DO
END DO
STL2 = DSQRT(SUM
*
TAU
*
H)
WRITE ( 01,
*
) ((FA(I,K), I = 1,N), K = 1,M)
WRITE ( 01,
*
) ((F(I,K), I = 1,N), K = 1,M)
WRITE ( 01,
*
) STL2
CLOSE ( 01 )
STOP
END
DOUBLE PRECISION FUNCTION AK(X)
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C COEFFICIENT AT THE HIGHER DERIVATIVES
C
AK = 1.D0
C
RETURN
END
DOUBLE PRECISION FUNCTION AF ( X, T )
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
C RIGHT-HAND SIDE
C
AF = 10.
*
T
*
(1.D0 - T)
*
X
*
(1.D0-X)
C
RETURN
END
The equation coefficients, dependent on the spatial variable, and the right-hand side
are to be set in the subroutine-functions AK and AF.
6.2.5 Computational experiments
The calculations were carried out on a uniform grid with the number of modes N =
100, N
0
= 100, the calculation domain being a unit square (l = 1, T = 1). The inverse
problem was solved within the framework of a quasi-real experiment in the case of
k(x) = 1, f (x, t) = 10t (1 − t)x(1 − x ).
The solution of the problem at the inaccuracy level δ = 0.001 is shown in Figure 6.5
(the total number of iterations is 3). Shown are the level curves for the exact (dashed
lines) and approximate solutions drawn with the step = 0.1. The effect of the
input-data inaccuracy on the right-hand side reconstruction accuracy is illustrated by
Figure 6.6, 6.7 that show data obtained with greater and lesser input-data inaccuracy.
For the problem to be solved with δ = 0.0001, a total of 23 iterations were required.
190 Chapter 6 Right-hand side identification
Figure 6.5 Solution of the identification problem obtained with δ = 0.001
Figure 6.6 Solution of the identification problem obtained with δ = 0.0001
Of course, the identification accuracy essentially depends on the exact solution.
In particular, noted above was the necessity to narrow the class of sought right-hand
sides due to homogeneous conditions given on some part of the calculation-domain
boundary. Figure 6.8 shows data obtained, with δ = 0.001, for the problem with the
right-hand side
f (x, t) = 2x(1 − x).
Section 6.3 Reconstruction of the time-dependent right-hand side 191
Figure 6.7 Solution of the identification problem obtained with δ = 0.01
Figure 6.8 Solution of the identification problem obtained with δ = 0.001
6.3 Reconstruction of the time-dependent right-hand side
In the theory of inverse heat-transfer problems, also problems are considered in which
it is required to reconstruct unknown heat sources from additional temperature mea-
surements made at some single points. In a similar manner, some important applied
problems are formulated that arise in hydrogeology. The solution can be made unique
by narrowing the class of admissible right-hand sides of the parabolic equations. In
many cases, one can assume the time-dependent right-hand side to be an unknown
quantity.
192 Chapter 6 Right-hand side identification
6.3.1 Inverse problem
In this section, we construct a computational algorithm for approximate solution of a
simplest one-dimensional (over space) inverse problem in which it is required to re-
construct the time-dependent right-hand side of a parabolic equation from the known
spatial distribution. Such a linear inverse problem belongs to the class of classically
well-posed mathematical physics problems under some special assumptions concern-
ing the points in space where additional measurements were performed, namely, under
the condition that the source acts at the observation points.
We assume the state of the system to be defined by the equation
∂u
∂t
=
∂
∂x
k(x)
∂u
∂x
+ f (x, t), 0 < x < l, 0 < t ≤ T (6.85)
with a sufficiently smooth positive coefficient k(x). We restrict ourselves to the prob-
lem with the simplest first-kind homogeneous boundary conditions:
u(0, t) = 0, u(l, t) = 0, 0 ≤ t ≤ T. (6.86)
Also given is the initial condition:
u(x, 0) = u
0
(x), 0 < x < l. (6.87)
The direct problem is formulated in the form (6.85)–(6.87).
We consider the inverse problem in which, apart from u(x, t ), also unknown is the
right-hand side f (x, t) of equation (6.85). We assume that the function f (x, t) can be
represented as
f (x, t) = η(t)ψ(x), (6.88)
where the function ψ(x) is given, and unknown is the time-dependent source, the
function η(t) in the representation (6.88). This dependence can be reconstructed from
additional observations of u(x, t) made at some internal point 0 < x
∗
< l:
u(x
∗
, t) = ϕ(t). (6.89)
We arrive at a simple right-hand side identification problem for the parabolic equation
(6.85)–(6.89).
6.3.2 Boundary value problem for the loaded equation
We consider the solution of the identification problem under the following constraints:
1) ψ(x
∗
) = 0;
2) ψ(x) is a sufficiently smooth function (ψ ∈ C
2
[0, 1]);
3) ψ(x) = 0 on the boundary of the calculation domain.
Section 6.3 Reconstruction of the time-dependent right-hand side 193
Particular attention should be given to the first assumption stating that the source to
be reconstructed acts at the observation point x
∗
. It is this assumption that makes the
identification problem under consideration a well-posed problem, providing for con-
tinuous dependence of the solution on initial data, right-hand side, and data measured
at the internal point.
The last assumption is of no fundamental significance; it is used just to simplify the
consideration.
We seek the solution of the inverse problem in the form
u(x, t) = θ(t)ψ(x) + w(x, t), (6.90)
where
θ(t) =
t
0
η(s) ds. (6.91)
Substitution of (6.90), (6.91) into (6.85), (6.88) yields the following equation for
w(x, t):
∂w
∂t
=
∂
∂x
k(x)
∂w
∂x
+ θ(t)
∂
∂x
k(x)
∂ψ
∂x
. (6.92)
With representation (6.90), the condition (6.89) yields the following representation for
the unknown function θ(t):
θ(t) =
1
ψ(x
∗
)
(ϕ(t) − w(x
∗
, t)). (6.93)
Substitution of (6.93) into (6.92) yields the sought loaded parabolic equation
∂w
∂t
=
∂
∂x
k(x)
∂w
∂x
+
1
ψ(x
∗
)
(ϕ(t) − w(x
∗
, t))
∂
∂x
k(x)
∂ψ
∂x
. (6.94)
Under the adopted assumptions about the right-hand side, the boundary condition at
the domain boundary is
w(0, t) = 0,w(l, t) = 0, 0 ≤ t ≤ T. (6.95)
With allowance for (6.91), for the auxiliary function θ(t) we have:
θ(0) = 0. (6.96)
This allows us to use the initial condition
w(x, 0) = u
0
(x), 0 < x < l. (6.97)
In this way, the inverse problem (6.85)–(6.89) can be formulated as a boundary
value problem for the loaded equation (6.94)–(6.97) with the representation (6.91),
(6.93) for the unknown time-dependent source.
194 Chapter 6 Right-hand side identification
6.3.3 Difference scheme
Note major points in the numerical solution of the identification problem. The compu-
tational algorithm is based on the approximate solution of the initial boundary problem
for the loaded equation. Today, the numerical solution methods for such non-classical
mathematical physics problems still remain poorly developed.
We assume a uniform grid ¯ω with a grid size h to be introduced over the coordinate
x. We denote the grid nodes as x
i
= ih, i = 0, 1,...,N, Nh = 1, and let v = v
i
=
v(x
i
). For simplicity, we assume that the observation point x = x
∗
coincides with an
internal node i = k.
We pass from one time layer t
n
= nτ , n = 0, 1,...,N
0
, N
0
τ = T , τ>0tothe
next time layer t
n+1
using a purely implicit difference scheme for equation (6.94). At
the internal grid nodes over space we have:
w
n+1
− w
n
τ
= (aw
n+1
¯x
)
x
+
1
ψ
k
(ϕ
n+1
− w
n+1
k
)(aψ
¯x
)
x
. (6.98)
In the case of problems with a sufficiently smooth coefficient k(x) we put, for instance,
a
i
= k(0.5(x
i
+ x
i−1
)). We approximate the conditions (6.95) and (6.97); then, we
have:
w
n+1
0
= 0,w
n+1
N
= 0, n = 0, 1,...,N
0
− 1, (6.99)
w
0
i
= u
0
(x
i
), i = 1, 2,...,N − 1. (6.100)
In accordance with (6.93), from the solution of the difference problem (6.98)–
(6.100) we determine
θ
n+1
=
1
ψ
k
(ϕ
n+1
− w
n+1
k
), n = 0, 1,...,N
0
− 1 (6.101)
and, then, supplement these relations with the condition (see (6.96)) θ
0
= 0. With
(6.91), for the sought time dependence of the right-hand side we use the simplest
numerical-differentiation procedure:
η
n+1
=
θ
n+1
− θ
n
τ
, n = 0, 1,...,N
0
− 1. (6.102)
To realize the implicit scheme, it is necessary to dwell here on the matter of solution
of the difference problem.
6.3.4 Non-local difference problem and program realization
In spite of non-locality of the difference problem on the next time layer, we encounter
no substantial problem in the computational realization of scheme (6.98)–(6.100). We
write equation (6.98) at internal nodes in the form
w
n+1
i
τ
− (aw
n+1
¯x
)
x,i
+
1
ψ
k
(aψ
¯x
)
x,i
w
n+1
k
= g
n
i
(6.103)
Section 6.3 Reconstruction of the time-dependent right-hand side 195
with given right-hand side g
n
i
and given boundary conditions (6.99). We seek the
solution of system (6.99), (6.103) in the form
w
n+1
i
= y
i
+ w
n+1
k
z
i
, i = 0, 1,...,N. (6.104)
Substitution of (6.104) into (6.103) allows us to formulate the following difference
problems for the auxiliary functions y
i
and z
i
:
y
i
τ
− (ay
¯x
)
x,i
= g
n
i
, i = 1, 2,...,N − 1, (6.105)
y
0
= 0, y
N
= 0, (6.106)
z
i
τ
− (az
¯x
)
x,i
+
1
ψ
k
(aψ
¯x
)
x,i
= 0, i = 1, 2,...,N − 1, (6.107)
z
0
= 0, z
N
= 0. (6.108)
Afterwards, with (6.104) we can find w
n+1
k
:
w
n+1
k
=
y
k
1 − z
k
. (6.109)
Correctness of the algorithm is guaranteed by the fact that the denominator in
(6.109) is never zero. For the difference problem (6.107), (6.108) the following a priori
estimate can be established based on the maximum principle for difference schemes:
max
0≤i≤N
|z
i
|≤τ max
0<i<N
1
ψ
k
(aψ
¯x
)
x,i
.
As a result, for a sufficiently small τ = O(1) we have: |z
i
| < 1, i. e., the time step size
must be sufficiently small.
The difference problems (6.105), (6.106) and (6.107), (6.108) are standard prob-
lems, their numerical solution presenting no difficulties. In the one-dimensional case
under consideration, the usual three-point sweep algorithm can be employed.
In fact, the computational difficulty of the used computational algorithm is equiva-
lent to double solution of the direct problem. For this reason, the considered method
can be classified to optimal methods. Below, we give the text of a program that nu-
merically solves the inverse problem of interest.
Program PROBLEM7
C
C PROBLEM7 - RIGHT-HAND SIDE IDENTIFICATION
C ONE-DIMENSIONAL NON-STATIONARY PROBLEM
C UNKNOWN TIME DEPENDENCE
C
IMPLICIT REAL
*
8 ( A-H, O-Z )
C
PARAMETER ( DELTA = 0.0005D0, N = 101, M = 101 )