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

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

256 Chapter 7 Evolutionary inverse problems

Theorem 7.7 Let self-adjoint, time-independent and positively deﬁned operators S

and A be mutually permutable. Then, the solution of the difference optimal control

problem (7.60)–(7.63) satisﬁes the equation

n+1

− y

− Ay

= 0, n = 0, 1,...,2N

− 1, (7.76)

supplemented with the non-local condition (7.75).

Thus, for the difference optimal control problem we have established their relation

with problems with time-non-local conditions over the double time interval. Above,

the relation between non-local problems and optimal control problems was established

under the assumption that the initial operator A is a self-adjoint, time-independent

operator positively deﬁned in H. These restrictions are not related with the essence of

the problems, and for more general optimal control problems, in a similar manner, an

Euler equation in the form of a non-local problem can be constructed.

7.1.7 Program realization

For methods with non-locally perturbed initial conditions, there exists a problem with

computational realization of such regularizing algorithms. It should be noted that,

presently, simple and convenient computational schemes for numerical solution of

general non-local boundary value problems for mathematical physics equations are

lacking. In order to avoid computational difﬁculties that obscure the essence of

the problem, below we consider a simple inverse problem for the one-dimensional

parabolic equation with constant coefﬁcients. The solution of the difference problem

will be constructed using the variable separation method in the fast Fourier transform

technique.

Within the framework of a quasi-real experiment, we ﬁrst solve the direct, initially

boundary value problem:

∂u

∂t

−

∂

∂x

= 0, 0 < x < l, 0 < t ≤ T,

u(0, t) = 0, u(l, t) = 0, 0 < t ≤ T,

u(x, 0) = u

(x), 0 ≤ x ≤ l.

The solution of the problem at the end time (u(x, T )) is perturbed to subsequently use

this solution as input data for the inverse problem of reconstructing the initial solution

at t = 0.

The realization is based in the fast Fourier transform (subroutines SINT and

SINTI). A more detailed description was given above (program PROBLEM9).

The value of the regularization parameter is chosen considering the discrepancy.

The computational realization is based on using the sequence

= α

, q > 0

Section 7.1 Non-local perturbation of initial conditions 257

with given α

and q.

Program PROBLEM10

C PROBLEM10 - PROBLEM WITH INVERTED TIME

C ONE-DIMENSIONAL PROBLEM

C NON-LOCALLY DISTURBED INITIAL CONDITION

PARAMETER ( DELTA = 0.005,N=65,M=65)

DIMENSION U0(N), UT(N), UTD(N), UT1(N), U(N), U1(N)

+ ,X(N), ALAM(N-2), W(N-2), W1(N-2), WSAVE(4

C PARAMETERS:

C XL, XR - LEFT AND RIGHT END POINTS OF THE SEGMENT;

C N - NUMBER OF GRID NODES OVER SPACE;

C TMAX - MAXIMAL TIME;

C M - NUMBER OF GRID NODES OVER TIME;

C DELTA - INPUT-DATA INACCURACY LEVEL;

C Q - MULTIPLIER IN THE REGULARIZATION PARAMETER;

C U0(N) - INITIAL CONDITION TO BE RECONSTRUCTED;

C UT(N) - END-TIME SOLUTION OF THE DIRECT PROBLEM;

C UTD(N) - DISTURBED SOLUTION OF THE DIRECT PROBLEM;

C U(N) - APPROXIMATE SOLUTION OF THE INVERSE PROBLEM;

XL = 0.

XR = 1.

TMAX = 0.1

PI = 3.1415926

OPEN (01, FILE = ’RESULT.DAT’) !FILE TO STORE THE CALCULATED DATA

C GRID

H=(XR-XL)/(N-1)

TAU = TMAX / (M-1)

DOI=1,N

X(I) = XL + (I-1)

END DO

C DIRECT PROBLEM

C INITIAL CONDITION

DOI=1,N

U0(I) = AU(X(I))

END DO

C EIGENVALUES OF THE DIFFERENCE OPERATOR

DOI=1,N-2

ALAM(I) = 4./H

(SIN(PI

I/(2.

(N-1))))

END DO

C FORWARD FOURIER TRANSFORM

CALL SINTI(N-2,WSAVE)

DOI=2,N-1

258 Chapter 7 Evolutionary inverse problems

W(I-1) = U0(I)

END DO

CALL SINT(N-2,W,WSAVE)

C FOURIER INVERSION

DOI=1,N-2

QQ=1./(1.+TAU

ALAM(I))

W1(I) = QQ

(M-1)

W(I)

END DO

CALL SINT(N-2,W1,WSAVE)

DOI=2,N-1

UT(I) = 1./(2.

(N-1))

W1(I-1)

END DO

UT(1) = 0.

UT(N) = 0.

DOI=1,N-2

QQ=1./(1.+TAU

ALAM(I))

W1(I) = QQ

((M-1)/2)

W(I)

END DO

CALL SINT(N-2,W1,WSAVE)

DOI=2,N-1

UT1(I) = 1./(2.

(N-1))

W1(I-1)

END DO

UT1(1) = 0.

UT1(N) = 0.

C DISTURBING OF MEASURED QUANTITIES

DOI=2,N-1

UTD(I) = UT(I) + 2.

DELTA

(RAND(0)-0.5)

END DO

UTD(1) = 0.

UTD(N) = 0.

C INVERSE PROBLEM

C INPUT-DATA FOURIER TRANSFORM

DOI=2,N-1

W(I-1) = UTD(I)

END DO

CALL SINT(N-2,W,WSAVE)

C ITERATIVE PROCESS TO ADJUST THE REGULARIZATION PARAMETER

IT=0

ITMAX = 1000

ALPHA = 0.001

Q = 0.75

100IT=IT+1

C FOURIER INVERSION

DOI=1,N-2

QQ=1.+TAU

ALAM(I)

W1(I) = W(I) / (1. + ALPHA

(M-1))

END DO

CALL SINT(N-2,W1,WSAVE)

Section 7.1 Non-local perturbation of initial conditions 259

DOI=2,N-1

U(I) = 1./(2.

(N-1))

W1(I-1)

END DO

C DISCREPANCY

SUM = 0.D0

DO I = 2,N-1

SUM = SUM + (U(I)-UTD(I))

END DO

SL2 = SQRT(SUM)

IF ( IT.EQ.1 ) THEN

IND=0

IF ( SL2.LT.DELTA ) THEN

IND=1

Q = 1.D0/Q

END IF

ALPHA = ALPHA

GO TO 100

ELSE

ALPHA = ALPHA

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 SOLUTION

DOI=1,N-2

QQ=1.+TAU

ALAM(I)

W1(I) = QQ

((M-1)/2)

W(I) / (1. + ALPHA

(M-1))

((M-1)/2)

END DO

CALL SINT(N-2,W1,WSAVE)

DOI=2,N-1

U(I) = 1./(2.

(N-1))

W1(I-1)

END DO

U(1) = 0.

U(N) = 0.

DOI=1,N-2

QQ=1.+TAU

ALAM(I)

W1(I) = QQ

((M-1)/2)

W(I) / (1. + ALPHA

(M-1))

END DO

CALL SINT(N-2,W1,WSAVE)

DOI=2,N-1

U1(I) = 1./(2.

(N-1))

W1(I-1)

END DO

U1(1) = 0.

U1(N) = 0.

WRITE ( 01,

) (U0(I),I=1,N)

WRITE ( 01,

) (UT(I),I=1,N)

WRITE ( 01,

) (UT1(I),I=1,N)

WRITE ( 01,

) (UTD(I),I=1,N)

WRITE ( 01,

) (U(I),I=1,N)

WRITE ( 01,

) (U1(I),I=1,N)

260 Chapter 7 Evolutionary inverse problems

WRITE ( 01,

) (X(I),I=1,N)

WRITE ( 01,

) IT, ALPHA

CLOSE (01)

STOP

END

FUNCTION AU(X)

C INITIAL CONDITION

AU=X/0.3

IF (X.GT.0.3) AU = (1. - X) / 0.7

RETURN

END

The initial condition is set in the subroutine AU. With regard to the previous considera-

tion of the convergence, the approximate solution is calculated at the end time (t = T )

and at the time t = T /2.

7.1.8 Computational experiments

We solved a model problem whose input data were taken equal to the difference solu-

tion of the direct problem (l = 1, t = 0.1) with the initial condition

(x) =



x/0.3, 0 < x < 0.3,

(1 − x)/0.7, 0.3 < x < 1.

The computations were performed on a uniform grid with h = 1/64, τ = 1/640. The

results obtained by solving the direct problem are shown in Figure 7.1.

The effect due to inaccuracies was modeled by perturbing the input data (solution

of the direct problem at the time t = T ):

˜u(x, T ) = u(x, T ) + 2δ



σ(x) −



, x ∈ ω.

Here, σ(x) is a random function normally distributed over the interval [0, 1]. Fig-

ure 7.2 shows the exact and approximate solutions of the inverse problem at two char-

acteristic times, t = T and t = T/2, for δ = 0.005. Similar data obtained with a

greater (δ = 0.015) and lower (δ = 0.0025) inaccuracy are shown respectively in

Figures 7.3 and 7.4.

The data presented show that the solution reconstruction accuracy essentially de-

pends on the time. The end-time solution is reconstructed inaccurately, and the con-

vergence related with decreased input-data inaccuracy is indistinctly observed. A more

favorable situation is observed at t = T /2. It is by this time, as it was shown previ-

ously in the theoretical consideration of the problem, that the rate of convergence of

the approximations to the exact solution attains its highest.

Section 7.1 Non-local perturbation of initial conditions 261

Figure 7.1 Mesh solution of the direct problem

Figure 7.2 Solution of the inverse problem obtained with δ = 0.005

262 Chapter 7 Evolutionary inverse problems

Figure 7.3 Solution of the inverse problem obtained with δ = 0.01

Figure 7.4 Solution of the inverse problem obtained with δ = 0.0025

Section 7.2 Regularized difference schemes 263

7.2 Regularized difference schemes

An important class of solution methods for ill-posed evolutionary problems is related

with some perturbation applied to the initial equation. In the generalized inverse

method, such a perturbation can be applied immediately to the initial differential prob-

lem.

In a more natural approach, no auxiliary differential problem is considered; instead,

this approach uses a perturbation applied to the difference problem. In this way, regu-

larized difference schemes can be constructed.

7.2.1 Regularization principle for difference schemes

The regularization principle for difference schemes is presently recognized as a guid-

ing principle in improvement of difference schemes. For general two- and three-layer

schemes, recipes aimed at improving the quality of the difference schemes (their sta-

bility, accuracy and efﬁciency) can be formulated. Using this principle, one can ex-

amine stability and convergence of a broad class of difference schemes for boundary

value mathematical physics problems and develop solution algorithms for difference

problems.

Traditionally, the regularization principle is widely used to develop stable difference

schemes for well-posed problems involving partial differential equations. The same

uniform methodological base is used to construct difference schemes for conditionally

well-posed non-stationary mathematical physics problems. Weakly perturbing prob-

lem operators, one can exert control over the growth of the solution norm on passage

from one to the next time layer.

Absolutely stable difference schemes can be constructed, with the help of the regu-

larization principle, in the following manner:

1. For a given initial problem we construct some simplest difference scheme

(generating difference scheme) that does not possesses the required properties;

i. e., a scheme conditionally stable or even absolutely unstable.

2. We write the difference scheme in the standard (canonical) form for which

stability conditions are known.

3. The properties of the difference scheme (its stability) can be improved by per-

turbing the operators in the difference scheme.

Thus, the regularization principle for a difference scheme is based on using already

known results concerning conditional stability. Such criteria are given by the general

stability theory for difference schemes. From the latter standpoint, we can consider the

regularization principle as a means enabling an efﬁcient use of general results yielded

by the stability theory for difference schemes. The latter can be achieved by writing

difference schemes in rather a general canonical form and by formulation of easy-to-

check stability criteria.

264 Chapter 7 Evolutionary inverse problems

To illustrate the considerable potential the regularization principle offers in the real-

ization of difference schemes, we will use this general approach to construct absolutely

stable schemes for direct mathematical physics problems. As a model problem, below

we consider the ﬁrst boundary value problem for the parabolic equation. As a generat-

ing function, consider a conditionally stable explicit scheme. For an absolutely stable

scheme to be obtained, we perturb the operators in the explicit difference scheme. We

separately consider the variants with additively (standardly) and multiplicative (non-

standardly) perturbed operators of the generating difference scheme.

Consider a two-dimensional problem in the rectangle

 ={x | x = (x

, x

), 0 < x

< l

,α= 1, 2}.

In , we seek the solution of the parabolic equation

∂u

∂t

−



α=1

∂

∂x



k(x)

∂u

∂x



= 0, x ∈ , 0 < t < T , (7.77)

supplemented with the simplest homogeneous boundary conditions of the ﬁrst kind:

u(x, t) = 0, x ∈ ∂, 0 < t < T . (7.78)

Also, an initial condition is given:

u(x, 0) = u

(x), x ∈ . (7.79)

We assume that the coefﬁcient in (7.77) is sufﬁciently smooth and k(x) ≥ κ, κ>0,

x ∈ .

We put into correspondence to the differential problem (7.77)–(7.79) a differential-

difference problem by performing discretization over space. We assume for simplicity

that a grid with step sizes h

, α = 1, 2, uniform along each of the directions, is

introduced in the domain . As usually, ω is the set of internal nodes.

On the set of mesh functions y(x) such that y(x) = 0, x = ω, we deﬁne, by the

relation

y =−



α=1

¯x

)

, (7.80)

a difference operator , putting, for instance,

(x) = k(x

− 0.5h

, x

), a

(x) = k(x

, x

− 0.5h

In the Hilbert space H = L

(ω), we introduce the scalar product and the norm by the

relations

(y,w) =



x∈ω

y(x)w(x)h

, y=



(y, y).

Section 7.2 Regularized difference schemes 265

In H ,wehave = 

∗

≥ mE, m > 0. We pass from (7.77)–(7.79) to the differential-

operator equation

+ y = 0, 0 < t < T (7.81)

atagiven

y(0) = u

, x ∈ ω. (7.82)

Now, we are going to construct absolutely stable two-layer difference schemes for

problem (7.81), (7.82) around the regularization principle.

In compliance with the regularization principle, we ﬁrst choose some difference

scheme to start from. As such generating a scheme, we can consider the simplest

explicit scheme

n+1

− y

+ y

= 0, n = 0, 1,...,N

− 1, (7.83)

= u

, x ∈ ω, (7.84)

where N

τ = T .

We write the difference scheme (7.83), (7.84) in the canonical form for two-layer

operator-difference schemes

n+1

− y

+ Ay

= 0, t

∈ ω

(7.85)

with the operators

B = E, A = .

By Theorem 4.2, the condition

B ≥

A, t ∈ ω

(7.86)

is a condition necessary and sufﬁcient for the scheme (7.84), (7.85) to be stable in H

or, in other words, for the estimate

y

n+1



≤u



, t ∈ ω

to hold.

With the inequality  ≤E taken into account, we would like to obtain, from

the necessary and sufﬁcient stability conditions (7.86), the following constraints on the

time step size in the explicit scheme (7.83), (7.84):

τ ≤



In the case under consideration, =O(|h|

), where |h|

= h

+ h

, and for the

maximum admissible step size we have τ

= O(|h|