Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics
Подождите немного. Документ загружается.
386 Chapter 8 Other problems
Y(I,J) = 0.D0
END DO
END DO
DOK=2,M
C
C DIFFERENCE-SCHEME COEFFICIENTS IN THE DIRECT PROBLEM
C
DOI=1,N1
G(I) = UK(I,K)
END DO
CALL FDS (A(1), A(N+1), A(2
*
N+1), F, Y,
+ H1, H2, N1, N2, TAU, G)
C
C SOLUTION OF THE DIFFERENCE PROBLEM
C
IDEFAULT(1) = 0
IREPT = 0
CALL SBAND5 (N, N1, A(1), Y, F, EPSR, EPSA)
C
C INPUT DATA FOR THE CONJUGATE PROBLEM
C
DOI=1,N1
FIK(I,K) = Y(I,1)
END DO
END DO
C
C CONJUGATE STATE
C PURELY IMPLICIT DIFFERENCE SCHEME
C
C INITIAL CONDITION (ATt=T)
C
DOI=1,N1
DOJ=1,N2
Y(I,J) = 0.D0
END DO
END DO
DO K = M,2,-1
C
C DIFFERENCE-SCHEME COEFFICIENTS IN THE DIRECT PROBLEM
C
DOI=1,N1
G(I) = FIK(I,K)
END DO
CALL FDSS (A(1), A(N+1), A(2
*
N+1), F, Y,
+ H1, H2, N1, N2, TAU, G)
C
C SOLUTION OF THE DIFFERENCE PROBLEM
C
IDEFAULT(1) = 0
IREPT = 0
CALL SBAND5 (N, N1, A(1), Y, F, EPSR, EPSA)
C
C DISCREPANCY
C
DOI=1,N1
R(I,K) = Y(I,N2) - FS(I,K)
END DO
END DO
C
C QUICKEST DESCEND METHOD
Section 8.3 Identification of the boundary condition in two-dimensional problems 387
C
C AUXILIARY PROBLEM
C
C INITIAL CONDITION
C
DOI=1,N1
DOJ=1,N2
Y(I,J) = 0.D0
END DO
END DO
DOK=2,M
C
C DIFFERENCE-SCHEME COEFFICIENTS IN THE DIRECT PROBLEM
C
DOI=1,N1
G(I) = R(I,K)
END DO
CALL FDS (A(1), A(N+1), A(2
*
N+1), F, Y,
+ H1, H2, N1, N2, TAU, G)
C
C SOLUTION OF THE DIFFERENCE PROBLEM
C
IDEFAULT(1) = 0
IREPT = 0
CALL SBAND5 (N, N1, A(1), Y, F, EPSR, EPSA)
DOI=1,N1
AR(I,K) = Y(I,1)
END DO
END DO
C
C ITERATION PARAMETER
C
SUM1 = 0.D0
SUM2 = 0.D0
DOI=1,N1
DOK=2,M
SUM1 = SUM1 + R(I,K)
*
R(I,K)
SUM2 = SUM2 + AR(I,K)
*
AR(I,K)
END DO
END DO
SS = SUM1/SUM2
C
C NEXT APPROXIMATION
C
DOI=1,N1
DOK=2,M
UK(I,K) = UK(I,K) - SS
*
R(I,K)
END DO
END DO
C
C EXIT FROM THE ITERATIVE PROCESS BY THE DISCREPANCY CRITERION
C
SUM = 0.D0
DOI=1,N1
DOK=2,M
SUM = SUM + (FID(I,K) - FIK(I,K))
**
2
*
H1
*
TAU
END DO
END DO
SUM = SUM/((X1R-X1L)
*
TMAX)
SL2 = DSQRT(SUM)
388 Chapter 8 Other problems
IF ( SL2.GT.DELTA ) GO TO 100
C
C SOLUTION
C
WRITE ( 01,
*
) ((UK(I,K), I=1,N1), K=2,M)
WRITE ( 01,
*
) ((FID(I,K), I=1,N1), K=2,M)
CLOSE (01)
STOP
END
C
SUBROUTINE FLUX (U, X1, T, N1, M)
C
C BOUNDARY CONDITION
C
IMPLICIT REAL
*
8 ( A-H, O-Z )
DIMENSION U(N1,M), X1(N1), T(M)
DOI=1,N1
DOK=1,M
U(I,K) = X1(I)
*
T(K)/T(M)
C IF (T(K).GT.0.5D0
*
T(M)) U(I,K) = 2.D0
*
X1(I)
*
(T(M)-T(K))
END DO
END DO
C
RETURN
END
C
SUBROUTINE FDS (A0, A1, A2, F, U, H1, H2, N1, N2, TAU, G)
C
C GENERATION OF DIFFERENCE-SCHEME COEFFICIENTS
C FOR THE PARABOLIC EQUATION WITH CONSTANT COEFFICIENTS
C IN USING THE PURELY IMPLICIT SCHEME
C
IMPLICIT REAL
*
8 ( A-H, O-Z )
DIMENSION A0(N1,N2), A1(N1,N2), A2(N1,N2)
+ ,F(N1,N2), U(N1,N2), G(N1)
C
DOJ=2,N2-1
DOI=2,N1-1
A1(I-1,J) = 1.D0/(H1
*
H1)
A1(I,J) = 1.D0/(H1
*
H1)
A2(I,J-1) = 1.D0/(H2
*
H2)
A2(I,J) = 1.D0/(H2
*
H2)
A0(I,J) = A1(I,J) + A1(I-1,J) + A2(I,J) + A2(I,J-1)
+ + 1.D0/TAU
F(I,J) = U(I,J)/TAU
END DO
END DO
C
C BOUNDARY CONDITION OF THE SECOND KIND
C
DOJ=2,N2-1
A2(1,J) = 1.D0/(H2
*
H2)
A2(1,J-1) = 1.D0/(H2
*
H2)
A0(1,J) = A1(1,J) + A2(1,J) + A2(1,J-1)+ 1.D0/TAU
F(1,J) = U(1,J)/TAU
END DO
C
DOJ=2,N2-1
A2(N1,J) = 1.D0/(H2
*
H2)
A2(N1,J-1) = 1.D0/(H2
*
H2)
Section 8.3 Identification of the boundary condition in two-dimensional problems 389
A0(N1,J) = A1(N1-1,J) + A2(N1,J) + A2(N1,J-1)+ 1.D0/TAU
F(N1,J) = U(N1,J)/TAU
END DO
C
DOI=2,N1-1
A1(I,1) = 1.D0/(H1
*
H1)
A1(I-1,1) = 1.D0/(H1
*
H1)
A0(I,1) = A1(I,1) + A1(I-1,1) + A2(I,1) + 1.D0/TAU
F(I,1) = U(I,1)/TAU
END DO
C
DOI=2,N1-1
A1(I,N2) = 1.D0/(H1
*
H1)
A1(I-1,N2) = 1.D0/(H1
*
H1)
A0(I,N2) = A1(I,N2) + A1(I-1,N2) + A2(I,N2-1) + 1.D0/TAU
F(I,N2) = U(I,N2)/TAU + G(I)/H2
END DO
C
A0(1,1) = A1(1,1) + A2(1,1) + 1.D0/TAU
F(1,1) = U(1,1)/TAU
C
A0(N1,1) = A1(N1-1,1) + A2(N1,1) + 1.D0/TAU
F(N1,1) = U(N1,1)/TAU
C
A0(1,N2) = A1(1,N2) + A2(1,N2-1) + 1.D0/TAU
F(1,N2) = U(1,N2)/TAU + G(1)/H2
C
A0(N1,N2) = A1(N1-1,N2) + A2(N1,N2-1) + 1.D0/TAU
F(N1,N2) = U(N1,N2)/TAU + G(N1)/H2
C
RETURN
END
C
SUBROUTINE FDSS (A0, A1, A2, F, U, H1, H2, N1, N2, TAU, G)
C
C GENERATION OF DIFFERENCE-SCHEME COEFFICIENTS
C FOR THE CONJUGATE STATE
C
IMPLICIT REAL
*
8 ( A-H, O-Z )
DIMENSION A0(N1,N2), A1(N1,N2), A2(N1,N2)
+ ,F(N1,N2), U(N1,N2), G(N1)
C
DOJ=2,N2-1
DOI=2,N1-1
A1(I-1,J) = 1.D0/(H1
*
H1)
A1(I,J) = 1.D0/(H1
*
H1)
A2(I,J-1) = 1.D0/(H2
*
H2)
A2(I,J) = 1.D0/(H2
*
H2)
A0(I,J) = A1(I,J) + A1(I-1,J) + A2(I,J) + A2(I,J-1)
+ + 1.D0/TAU
F(I,J) = U(I,J)/TAU
END DO
END DO
C
C BOUNDARY CONDITION OF THE SECOND KIND
C
DOJ=2,N2-1
A2(1,J) = 1.D0/(H2
*
H2)
A2(1,J-1) = 1.D0/(H2
*
H2)
A0(1,J) = A1(1,J) + A2(1,J) + A2(1,J-1)+ 1.D0/TAU
390 Chapter 8 Other problems
F(1,J) = U(1,J)/TAU
END DO
C
DOJ=2,N2-1
A2(N1,J) = 1.D0/(H2
*
H2)
A2(N1,J-1) = 1.D0/(H2
*
H2)
A0(N1,J) = A1(N1-1,J) + A2(N1,J) + A2(N1,J-1)+ 1.D0/TAU
F(N1,J) = U(N1,J)/TAU
END DO
C
DOI=2,N1-1
A1(I,1) = 1.D0/(H1
*
H1)
A1(I-1,1) = 1.D0/(H1
*
H1)
A0(I,1) = A1(I,1) + A1(I-1,1) + A2(I,1) + 1.D0/TAU
F(I,1) = U(I,1)/TAU + G(I)/H2
END DO
C
DOI=2,N1-1
A1(I,N2) = 1.D0/(H1
*
H1)
A1(I-1,N2) = 1.D0/(H1
*
H1)
A0(I,N2) = A1(I,N2) + A1(I-1,N2) + A2(I,N2-1) + 1.D0/TAU
F(I,N2) = U(I,N2)/TAU
END DO
C
A0(1,1) = A1(1,1) + A2(1,1) + 1.D0/TAU
F(1,1) = U(1,1)/TAU + G(1)/H2
C
A0(N1,1) = A1(N1-1,1) + A2(N1,1) + 1.D0/TAU
F(N1,1) = U(N1,1)/TAU + G(N1)/H2
C
A0(1,N2) = A1(1,N2) + A2(1,N2-1) + 1.D0/TAU
F(1,N2) = U(1,N2)/TAU
C
A0(N1,N2) = A1(N1-1,N2) + A2(N1,N2-1) + 1.D0/TAU
F(N1,N2) = U(N1,N2)/TAU
C
RETURN
END
In the subroutine FLUX, boundary conditions of the second kind for the test direct
problem (exact inverse-problem solution) are introduced. In the subroutines FDS and
FDSS, difference-scheme coefficients for the ground and conjugate states are gener-
ated, respectively.
8.3.6 Computational experiments
The input data for the inverse problem are obtained as the solution of the direct problem
(8.91)–(8.95) with the boundary condition (8.94) given as
μ(x
1
, t) = x
1
t
T
.
The direct problem is being solved in the rectangle with l
1
= 1, l
2
= 0.5atT = 1
on the grid ω with N
1
= N
2
= 50 and τ = 0.01. The solution of the direct problem at
Section 8.3 Identification of the boundary condition in two-dimensional problems 391
the times t = 0.25, 0.5, 0.75, 1 are shown in Figures 8.12–8.15.
Figure 8.12 Direct-problem solution (t = 0.25)
Figure 8.13 Direct-problem solution (t = 0.5)
392 Chapter 8 Other problems
Figure 8.14 Direct-problem solution (t = 0.75)
Figure 8.15 Direct-problem solution (t = 1)
The found solution of the direct problem at the grid nodes closest to γ
∗
(the func-
tions ϕ
n
(x
1
) in (8.115)) at the same times are shown in Figure 8.16. These data were
randomly perturbed to be subsequently considered as input data for the inverse prob-
lem:
¯ϕ
n
(x
1
) = ϕ
n
(x
1
) + 2δ(σ(x
1
, t
n
) − 1/2), x
1
∈ ω, n = 1, 2,...,N
0
.
Section 8.3 Identification of the boundary condition in two-dimensional problems 393
Here, σ(x
1
, t
n
) is a random function distributed normally over the interval [0, 1]. The
approximate solution of the inverse problem obtained with the inaccuracy level defined
by the parameter δ = 0.005 is shown in Figure 8.17. The same figure shows the exact
and found boundary conditions at the times t = 0.25, 0.5, 0.75, and 1. The end-time
reconstruction accuracy (at t = 1) is seen to be very poor. In fact, here there is no other
thing to expect since the changes in the boundary condition at times around t = T do
not affect the solution at the observation points (on the other portion of the boundary):
the perturbations do not arrive at the observation points. The effect due to the accuracy
in setting the input data can be figured out considering Figures 8.18–8.19.
Figure 8.16 Unperturbed initial data for the inverse problem
Figure 8.17 Inverse-problem solution obtained with δ = 0.005
394 Chapter 8 Other problems
Figure 8.18 Inverse-problem solution obtained with δ = 0.01
Figure 8.19 Inverse-problem solution obtained with δ = 0.0025
8.4 Coefficient inverse problem for the nonlinear parabolic
equation
In this section, we consider the inverse problem in which it is required to determine
a coefficient that depends on the solution of the one-dimensional parabolic equation
from the solution observed at some internal point/points. Algorithms of functional and
Section 8.4 Coefficient inverse problem for the nonlinear parabolic equation 395
parametric optimization are discussed. Primary attention is paid to the construction
and realization of computational algorithms.
8.4.1 Statement of the problem
In many applied problems, there arises a problem in which it is required to identify
coefficients in a partial differential equation. Coefficient inverse problems for second-
order parabolic equations are typical of heat- and mass-transfer problems and problems
encountered in hydrogeology. The inverse problems in which it is required to identify
coefficients in linear equations are nonlinear problems. The latter circumstance sub-
stantially hampers the construction of computational algorithms for the approximate
solution of coefficient problems and makes complete and rigorous substantiation of
their convergence hardly possible. That is why the emphasis here is placed on maxi-
mum possible approbation of numerical methods aimed at obtaining most informative
solution examples for inverse problems.
Very often, of primary interest are problems in which it is required to find nonlinear
coefficients that depend on the solution. Let us formulate a simplest of such a problem.
Suppose that 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(u)
∂u
∂x
= 0, 0 < x < l, 0 < t ≤ T. (8.132)
We assume that k(u) ≥ κ>0. Consider the boundary value problem with the first-
kind boundary conditions
u(0, t) = 0, u(l, t) = g(t), 0 < t ≤ T (8.133)
and the homogeneous initial conditions
u(x, 0) = 0, 0 ≤ x ≤ l. (8.134)
The direct problem is formulated in the form (8.132)–(8.134).
In the coefficient inverse problem, the unknown function k(u) is to be found, for
instance, from additional observations performed at some internal points z
m
∈ ,
m = 1, 2,...,M. Considering the fact that these data are given inaccurate, we put
u(z
m
, t) ≈ ϕ
m
(t), 0 < t ≤ T, m = 1, 2,...,M. (8.135)
It is required to find the functions u(x, t) and k(u) from the conditions (8.132)–(8.135).
In the consideration of coefficient inverse problems similar to problem (8.132)–
(8.135), much attention is paid to problems of solution unicity in the case of exact