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

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

116 Chapter 4 Boundary value problems for parabolic equations

of implicit approximations, at the next time layer we arrive at nonlinear difference

equations. For the approximate solution to be found, we have to use these or those

iteration methods for solving systems of nonlinear equations. To avoid this situation,

in computational practice they widely use linearized difference schemes in which the

solution at the next layer is found from a system of linear equations. We will illustrate

some possibilities that arise along this line with the example of difference schemes for

the nonlinear problem (4.117)–(4.119).

By analogy with (4.11), on the set of mesh functions given on ¯ω and vanishing on

∂ω, we deﬁne the operator

A(v)y =−(a(x,v)y

¯x

)

, x ∈ ω.

Here, the coefﬁcient a(x,v)is to be deﬁned, for instance, as

a(x,v) = k(x − 0.5h, 0.5(v(x) + v(x − h))),

a(x,v) =

(k(x − h,v(x − h)) + k(x, v)).

To the initial differential problem (4.117)–(4.119), we put into correspondence the

differential-difference problem

+ A(y)y = f (x, t, y), x ∈ ω, t > 0, (4.124)

y(x, 0) = u

(x), x ∈ ω. (4.125)

Let us discuss difference schemes for problem (4.124), (4.125). To begin with, con-

sider nonlinear difference schemes in which the solution at the next time layer is found

as the solution of a nonlinear difference problem. Such schemes can be constructed

similarly to difference schemes for the linear parabolic equation. For instance, a purely

implicit difference scheme for (4.124), (4.125) is

n+1

− y

+ A(y

n+1

= f (x, t

n+1

, y

n+1

n = 0, 1,...,N

− 1,

(4.126)

= u

(x), x ∈ ω. (4.127)

In the nonlinear scheme (4.126), (4.127), for the difference solution at the next time

layer to be found, a nonlinear difference problem must be solved. For the values of

n+1

to be determined, these or those iterative processes are used. Some important

speciﬁc features of the corresponding nonlinear difference problems deserve mention.

The ﬁrst speciﬁc feature is related with the fact that in the iterative embodiment

of implicit difference schemes a good initial approximation is always available. This

initial approximation can be chosen as the solution at the previous layer. The second,

Section 4.5 Program realization and computation examples 117

also beneﬁcial speciﬁc feature is the fact that the difference problem for y

n+1

involves

a small parameter, the time step size τ . This parameter essentially affects the rate of

convergence of the iterative process (the smaller τ , the higher is the rate of convergence

in the process).

We can also identify the class of linearized difference schemes whose speciﬁc fea-

ture consists in that, here, the solution at the next time layer is to be found from the

solution of a linear difference problem. A simplest scheme is characterized by that,

in it, the coefﬁcients are taken from the previous time layer. An example here is the

difference scheme

n+1

− y

+ A(y

, y

n+1

) = f (x, t

, y

). (4.128)

This difference scheme has, apparently, an approximation inaccuracy of order

O(τ +|h|

The main drawback of scheme (4.128) is often related with the fact that the right-

hand side (source) is taken at the previous time layer. Further development of the

difference scheme (4.128) is the scheme with quasi-linearized right-hand side:

n+1

− y

+ A(y

, y

n+1

) = f (x, t

n+1

, y

) +

∂ f

∂y

(x, t

n+1

, y

)(y

n+1

− y

). (4.129)

Scheme (4.129), which remains linear, has wider stability margins with respect to the

nonlinear right-hand side.

Linearized schemes can be constructed around the difference schemes “predictor-

corrector”, which conceptually border on the additive difference schemes (split

schemes). Let us restrict ourselves here to a simplest variant of the predictor-corrector

scheme:

˜y

n+1

− y

+ A(y

, y

) = f (x, t

, y

). (4.130)

Scheme (4.130) is used to calculate the coefﬁcients and the right-hand side; that is why

at the correction stage we can use the scheme

n+1

− y

+ A( ˜y

n+1

, y

n+1

) = f (x, t

n+1

, ˜y

n+1

). (4.131)

Alternatively, the linearization scheme (4.129) can be used at the correction stage.

The above schemes show considerable potential in constructing linearized differ-

ence schemes. In practical simulations, special methodical studies need to be per-

formed aimed at making a proper choice of difference schemes for a particular class

of nonlinear boundary value problems. In the case of nonlinear problems, theoretical

considerations give only obscure guiding lines.

Of course, linearized difference schemes for problem (4.124), (4.125) can also be

constructed around three-layer difference schemes. Without dwelling on a detailed

118 Chapter 4 Boundary value problems for parabolic equations

description, let us give here only two simplest examples. The problem (4.124), (4.125)

can be approximately solved using the three-layer symmetric difference schemes

n+1

− y

n−1

2τ

+ A(y

,σy

n+1

+ (1 − 2σ)y

+ σ y

n−1

) = f (x, t

, y

) (4.132)

with appropriate initial and boundary conditions. These linearized schemes are of the

second approximation order both over space and time.

4.5.3 Program

Below, we give the text of a program in which two linearized schemes, schemes (4.128)

and (4.129) were implemented. To ﬁnd the solution at the next time layer, the standard

sweep algorithm was employed.

Program PROBLEM3

C PROBLEM3 - ONE-DIMENSIONAL NON-STATIONARY PROBLEM

C QUASI-LINEAR 1D PARABOLIC EQUATION

IMPLICIT REAL

8 ( A-H, O-Z )

PARAMETER ( ISCHEME = 0, N = 1000, M = 1000 )

DIMENSION X(N+1), Y(N+1), A(N+1), B(N+1), C(N+1), F(N+1)

+ ,ALPHA(N+2), BETA(N+2)

C PROBLEM PARAMETERS:

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

C ISCHEME - PARAMETER FOR THE CHOICE OF THE SCHEME,

C N + 1 - NUMBER OF GRID NODES OVER SPACE;

C M + 1 - NUMBER OF GRID NODES OVER TIME;

XL = 0.D0

XR = 10.D0

TMAX = 2.5D0

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

C GRID

H =(XR-XL)/N

TAU = TMAX / M

DOI=1,N+1

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

END DO

C INITIAL CONDITION

T = 0.D0

XD = 2.5D0

DOI=1,N+1

Section 4.5 Program realization and computation examples 119

Y(I) = 0.D0

IF ( X(I).LT.XD ) Y(I) = 1.D0

END DO

WRITE ( 01,

) (Y(I),I=1,N+1)

C NEXT TIME LAYER

DOK=1,M

T=K

TAU

C DIFFERENCE-SCHEME COEFFICIENTS

IF ( ISCHEME.EQ.0 ) THEN

C LINEARIZED SCHEME WITH EXPLICIT RIGHT-HAND SIDE

DOI=2,N

X1 = (X(I) + X(I-1)) / 2

X2 = (X(I+1) + X(I)) / 2

U1 = (Y(I) + Y(I-1)) / 2

U2 = (Y(I+1) + Y(I)) / 2

A(I) = AK(X1,U1) / (H

B(I) = AK(X2,U2) / (H

C(I) = A(I) + B(I) + 1.D0 / TAU

F(I) = Y(I) / TAU + AF(X(I),T,Y(I))

END DO

END IF

IF ( ISCHEME.EQ.1 ) THEN

C LINEARIZED SCHEME WITH QUASI-LINEARIZED RIGHT-HAND SIDE

DOI=2,N

X1 = (X(I) + X(I-1)) / 2

X2 = (X(I+1) + X(I)) / 2

U1 = (Y(I) + Y(I-1)) / 2

U2 = (Y(I+1) + Y(I)) / 2

A(I) = AK(X1,U1) / (H

B(I) = AK(X2,U2) / (H

C(I) = A(I) + B(I) + 1.D0 / TAU - ADF(X(I),T,Y(I))

F(I) = Y(I) / TAU + AF(X(I),T,Y(I)) - ADF(X(I),T,Y(I))

Y(I)

END DO

END IF

C LEFT BOUNDARY CONDITION

B(1) = 0.D0

C(1) = 1.D0

F(1) = 1.D0

C RIGHT BOUNDARY CONDITION

A(N+1) = 0.D0

C(N+1) = 1.D0

F(N+1) = 0.D0

120 Chapter 4 Boundary value problems for parabolic equations

C SOLUTION OF THE PROBLEM AT THE NEXT TIME LAYER

CALL PROG ( N+1, A, C, B, F, ALPHA, BETA, Y )

C SOLUTION RECORDING

IF ( K/200

200.EQ.K ) THEN

WRITE ( 01,

) (Y(I),I=1,N+1)

END IF

END DO

CLOSE ( 01 )

STOP

END

SUBROUTINE PROG ( N, A, C, B, F, AL, BET, Y )

IMPLICIT REAL

8 ( A-H, O-Z )

C SWEEP METHOD

DIMENSION A(N), C(N), B(N), F(N), Y(N), AL(N+1), BET(N+1)

AL(1) = B(1) / C(1)

BET(1) = F(1) / C(1)

DOI=2,N

SS = C(I) - AL(I-1)

A(I)

AL(I) = B(I) / SS

BET(I) = ( F(I) + BET(I-1)

A(I))/SS

END DO

Y(N) = BET(N)

DO I = N-1, 1, -1

Y(I) = AL(I)

Y(I+1) + BET(I)

END DO

RETURN

END

DOUBLE PRECISION FUNCTION AK ( X, U )

IMPLICIT REAL

8 ( A-H, O-Z )

C COEFFICIENT AT THE HIGHER DERIVATIVE

AK = 0.2D0

RETURN

END

DOUBLE PRECISION FUNCTION AF ( X, T, U )

IMPLICIT REAL

8 ( A-H, O-Z )

C RIGHT-HAND SIDE OF THE EQUATION

AF = 5.D0

(1.D0 - U)

Section 4.5 Program realization and computation examples 121

RETURN

END

DOUBLE PRECISION FUNCTION ADF ( X, T, U )

IMPLICIT REAL

8 ( A-H, O-Z )

C DIRIVATIVE OF THE RIGHT-HAND SIDE OF THE EQUATION

ADF = 5.D0 - 10.D0

RETURN

END

The above program text refers to the case in which equation (4.117) is solved with

the coefﬁcients

k(x, u) = κ, f (x, t, u) = χu(1 − u), (4.133)

where κ = 0.2 and χ = 5. The latter nonlinear diffusion-reaction equation is known

as the Fisher equation (Kolmogorov–Petrovskii–Piskunov equation).

4.5.4 Computational experiments

Very often, nonlinearity brings about new effects related with an unusual behavior of

the solution. The latter is also the case with the quasi-linear parabolic equation (4.117).

A typical example here is the Fisher equation (4.117), (4.133).

For the linear parabolic-type equation (4.1), inﬁnite propagation velocity of per-

turbations is typical. In consideration of equations with nonlinear coefﬁcients and

nonlinear right-hand side, one can identify solutions with a ﬁnite velocity of perturba-

tions. In the case of the Fisher equation (4.117), (4.133), there exist solutions of the

traveling-wave type:

u(x, t) = ψ(ξ), ξ = x − ct,

where c = const is the wave velocity.

A. N. Kolmogorov, I. G. Petrovskii and N. S. Piskunov showed that, for instance,

with

u(x, 0) =



1, x < x

∗

0, x > x

∗

the solution of the Cauchy problem for equation (4.117), (4.133) is unique, presenting

a wave traveling at the velocity c = 2

√

κχ. Figure 4.1 shows the solution of the

122 Chapter 4 Boundary value problems for parabolic equations

boundary value problem for equation (4.117), (4.133) with κ = 0.2, χ = 5 and with

the following boundary and initial conditions:

u(0, t) = 1, u(10, t) = 1,

u(x, 0) =

1, x < x

∗

0, x > x

∗

0 < x < 10

where x

∗

= 2.5.

Figure 4.1 The solution at various times

In the calculations, a ﬁne calculation grid with h = 0.01,τ= 0.0025 and scheme

(4.129) were used. Transformation of the initially rectangular waveform into a solitary

wave traveling to the right is observed.

Some speciﬁc features of schemes in which linearization of the right-hand side of

the equation is used are illustrated by Figures 4.2 and 4.3 that show calculation data

obtained with a larger time step size (τ = 0.025) using the linearized schemes (4.128)

(Figure 4.2) and (4.129) (Figure 4.3). The scheme with explicit right-hand side yields

more accurate results, whereas in the case of the scheme with right-hand side partially

transferred to the next time layer the front of the approximate solution propagates

somewhat faster. It should be noted that the conclusion about certain advantages of

non-linearized schemes concerns only the considered problem; in other cases, a lin-

earized scheme may yield better results.

Section 4.5 Program realization and computation examples 123

Figure 4.2 Difference scheme with non-linearized right-hand side

Figure 4.3 Difference scheme with linearized right-hand side

124 Chapter 4 Boundary value problems for parabolic equations

4.6 Exercises

Exercise 4.1 Consider the Cauchy problem for the second-order evolutionary equa-

tion:

+ Au = f (t), t > 0, (4.134)

u(0) = u

(0) = v

. (4.135)

Show that in the case of A = A

∗

> 0 the following a priori estimate holds for the

solution of problem (4.134), (4.135):

u(t)

∗

≤ exp(t)



u



+v



exp(−θ)f (θ)

dθ



where

u

∗



+u

Exercise 4.2 Construct a two-layer weighted scheme for the problem (4.1), (4.3) with

the third-kind boundary conditions (4.4).

Exercise 4.3 Suppose that we seek an approximate solution of the equation

∂u

∂t

∂

∂x

+ f (x, t), 0 < x < l, 0 < t ≤ T,

supplemented with conditions (4.2), (4.3). Construct a two-layer scheme with an ap-

proximation order O(τ

+ h

Exercise 4.4 In the class of weighted schemes



n+1

− y

+ (1 − θ)

− y

n−1



+ A(σ

n+1

+ (1 − σ

− σ

+ σ

n−1

) = ϕ

construct a scheme with an approximation order O(τ

Exercise 4.5 Prove Theorem 4.3.

Exercise 4.6 Establish stability conditions for the scheme

n+1

− 4y

+ y

n−1

2τ

+ Ay

n+1

= 0

in the case of A = A

∗

> 0.

Section 4.6 Exercises 125

Exercise 4.7 Suppose that in an examination of the weighted scheme

n+1

− y

+ A(σ y

n+1

+ (1 − σ))y

= ϕ

(4.136)

the pure implicit scheme (with σ = 1) was found to be stable with respect to initial

data and right-hand side. Show that, in this case, all schemes with σ>1 are also

stable.

Exercise 4.8 Suppose that for the scheme

n+1

− y

+ Ay

= ϕ

, t

∈ ω

= u

with a constant operator A = A

∗

> 0 there holds the inequality

B ≥

Show that for the solution there holds the following stability estimate with respect

to initial data and right-hand side:

y

n+1



≤u



+ϕ



−1

+ϕ



−1



k=1



− ϕ

k−1



Exercise 4.9 Based on the maximum principle for difference schemes, establish con-

ditions of stability in C(ω) for the weighted scheme (4.136) with

Ay =−(ay

¯x

)

, x ∈ ω

used to solve the boundary value problem (4.1)–(4.3).

Exercise 4.10 Examine stability of the three-layer symmetric scheme

n+1

− y

n−1

2τ

+ A(σ y

n+1

+ (1 − 2σ)y

+ σ y

n−1

) = ϕ

used to solve the boundary value problem (4.1)–(4.3).

Exercise 4.11 In the solution of problem (4.134), (4.135), one can naturally use the

weighted scheme

n+1

− 2y

+ y

n−1

+ A(σ

n+1

+ (1 − σ

− σ

+ σ

n−1

) = ϕ

n = 1, 2,...,N

− 1.

Obtain stability conditions for this scheme.