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

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

36 Chapter 2 Boundary value problems for ordinary differential equations

Figure 2.2 Exact solution of the convection-diffusion problem

In the case of Pe  1, we have a process with dominating diffusion, whereas in the

case of Pe  1 convection prevails. In the former case, we arrive at regularly per-

turbed problems (small parameter Pe

−1

at low-order derivatives), whereas in the case

of strongly dominating convection, we obtain singularly perturbed problems (small

parameter Pe

−1

at hogh-order derivatives). Typical for singularly perturbed problems

is the occurrence of regions where the solution displays considerable variations, such

regions, for instance, being boundary layers and internal transition layers.

Computational algorithms for approximate solution of convection-diffusion prob-

lems can be understood considering the problem with constant coefﬁcients and homo-

geneous boundary conditions:

k(x) = κ, v(x) = 1, f (x) = 1,μ

= 0,μ

= 0. (2.67)

The exact solution of problem (2.64), (2.65), (2.66) is

u(x) = x −l

exp (x/κ) − 1

exp (l/κ) − 1

Speciﬁc features inherent to the problem with dominating convection (low diffusion

coefﬁcients) can be ﬁgured out considering Figure 2.2 that shows the solution of the

problem with l = 1 and κ = 1, 0.1, and 0.01. As the diffusion coefﬁcient decreases,

near the right boundary there forms a boundary layer (region with large solution gra-

dients).

Section 2.4 Program realization and computational examples 37

2.4.2 Difference schemes

In consideration of one-dimensional convection-diffusion problems (2.64), (2.65), we

choose using difference schemes written at the internal nodes in the form (2.51), (2.52).

Consider the difference schemes (2.51), (2.52) in which

> 0,β

> 0,γ

> 0, i = 1, 2,...,N − 1. (2.68)

Let us formulate a criterion for monotonicity of the difference scheme or, in other

words, formulate conditions under which the difference scheme (2.51), (2.52) satisﬁes

the difference principle of maximum.

Theorem 2.5 (Maximum principle) Let in the difference scheme (2.51), (2.52),

(2.68) μ

≥ 0, μ

≥ 0 and ϕ

≥ 0 for all i = 1, 2,...,N − 1 (or, alternatively,

≤ 0, μ

≤ 0 and ϕ

≤ 0 for i = 1, 2,...,N − 1). Then, provided that

≥ α

+ β

, i = 1, 2,...,N − 1, (2.69)

we have: y

≥ 0,i = 1, 2,...,N − 1 (y

≤ 0,i = 1, 2,...,N − 1).

Proof. Let us follow the line of reasoning assuming the opposite. Let conditions (2.69)

be fulﬁlled, but the difference solution of problem (2.51) with non-negative right-hand

side and non-negative boundary conditions be not non-negative at all nodes of the grid.

We designate as k the grid node at which the solution assumes the least negative value.

If such a value is attained at several nodes, then we choose the node where y

k−1

> y

We write the difference equation at this node:

−α

k−1

+ γ

− β

k+1

= ϕ

The right-hand side is non-negative, and for the left-hand side, in view of (2.68) and

(2.69), we have:

− α

k−1

+ γ

− β

k+1

= α

− y

k−1

) + (γ

− α

− β

+ β

− y

k+1

)>0.

The obtained contradiction shows that y

≥ 0 at all nodes i = 1, 2,...,N − 1.

To approximately solve the problem (2.64), (2.65), we use the simplest difference

scheme with central-difference approximation of the convective term. At the inter-

nal nodes of the computational grid, we approximate the differential equation (2.64),

accurate to the second order, with the difference equation

−(ay

¯x

)

+ by

◦

= ϕ, x ∈ ω, (2.70)

where, for instance, a(x) = k(x − 0.5h) and b(x) = v(x) in problems with smooth

coefﬁcients.

38 Chapter 2 Boundary value problems for ordinary differential equations

We write the difference scheme (2.52), (2.70) in the form (2.51), (2.52) with

i−1/2

,β

=−

i+1/2

i−1/2

+ k

i+1/2

), i = 1, 2,...,N − 1.

The sufﬁcient conditions for monotonicity (2.68), (2.69) (conditions for the max-

imum principle) for the scheme with the convection term approximated with cen-

tral differences will be fulﬁlled only in the case of sufﬁciently small grid steps (low

convective-transfer coefﬁcients). Disregarding the sign of velocity, we write the con-

straints as

≡

min{k

i−1/2

, k

i+1/2

}

< 2. (2.71)

Here Pe

is the mesh Peclet number. Hence, to obtain a monotone difference scheme,

we have to use a sufﬁciently ﬁne computation grid. With constraints (2.71), we may

speak only of conditional monotonicity of the difference scheme (2.52), (2.70).

Absolutely monotone difference schemes for convection-diffusion problems can be

constructed using ﬁrst-order approximations of the convection term with directed dif-

ferences. We use the settings

b(x) = b

(x) + b

−

(x),

(x) =

(b(x) +|b(x)|) ≥ 0, b

−

(x) =

(b(x) −|b(x)|) ≤ 0

for the non-positive and non-negative parts of the mesh function b(x), x ∈ ω. Instead

of (2.70), we use the difference scheme

−(ay

¯x

)

+ b

¯x

+ b

−

= ϕ, x ∈ ω. (2.72)

For scheme (2.52), (2.72), we have representation (2.5), (2.52) with coefﬁcients

i−1/2

,β

=−

−

i+1/2

i−1/2

+ k

i+1/2

), i = 1, 2,...,N − 1.

We see immediately that, here, sufﬁcient conditions for monotonicity (conditions

(2.68) and (2.69)) are fulﬁlled.

The scheme with directed differences is monotone with any grid steps. Its deﬁ-

ciencies are related with the fact that, unlike the scheme with central differences, this

scheme has ﬁrst approximation order.

Section 2.4 Program realization and computational examples 39

2.4.3 Program

Realization of the difference schemes (2.52), (2.70) and (2.52), (2.72) implies numer-

ical solution of a system of linear algebraic equations with a tridiagonal matrix. As it

was noted previously, such problems can be solved using the sweep algorithm. This

algorithm is correct under the condition of diagonal prevalence (see Lemma 2.4). For

non-monotone difference schemes of type (2.52), (2.70), one has to use more general

algorithms. Here we speak of non-monotone sweep, which can be considered as the

Gauss method with the choice of major element in solving the tridiagonal system of

linear algebraic equations.

We use the subroutine SGTSL from the well-known linear algebra package LIN-

PACK.

SUBROUTINE SGTSL (N, C, D, E, B, INFO)

***

BEGIN PROLOGUE SGTSL

***

PURPOSE Solve a tridiagonal linear system.

***

LIBRARY SLATEC (LINPACK)

***

CATEGORY D2A2A

***

TYPE SINGLE PRECISION (SGTSL-S, DGTSL-D, CGTSL-C)

***

KEYWORDS LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE, TRIDIAGONAL

***

AUTHOR Dongarra, J., (ANL)

***

DESCRIPTION

C SGTSL given a general tridiagonal matrix and a right hand

C side will find the solution.

C On Entry

C N INTEGER

C is the rank of the tridiagonal matrix.

C C REAL(N)

C is the subdiagonal of the tridiagonal matrix.

C C(2) through C(N) should contain the subdiagonal.

C On output, C is destroyed.

C D REAL(N)

C is the diagonal of the tridiagonal matrix.

C On output, D is destroyed.

C E REAL(N)

C is the superdiagonal of the tridiagonal matrix.

C E(1) through E(N-1) should contain the superdiagonal.

C On output, E is destroyed.

C B REAL(N)

C is the right hand side vector.

C On Return

C B is the solution vector.

C INFO INTEGER

C = 0 normal value.

C = K if the Kth element of the diagonal becomes

40 Chapter 2 Boundary value problems for ordinary differential equations

C exactly zero. The subroutine returns when

C this is detected.

***

REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.

C Stewart, LINPACK Users’ Guide, SIAM, 1979.

***

ROUTINES CALLED (NONE)

***

REVISION HISTORY (YYMMDD)

C 780814 DATE WRITTEN

C 890831 Modified array declarations. (WRB)

C 890831 REVISION DATE from Version 3.2

C 891214 Prologue converted to Version 4.0 format. (BAB)

C 900326 Removed duplicate information from DESCRIPTION section.

C (WRB)

C 920501 Reformatted the REFERENCES section. (WRB)

***

END PROLOGUE SGTSL

INTEGER N,INFO

REAL C(

),D(

),E(

),B(

)

INTEGER K,KB,KP1,NM1,NM2

REAL T

***

FIRST EXECUTABLE STATEMENT SGTSL

INFO = 0

C(1) = D(1)

NM1=N-1

IF (NM1 .LT. 1) GO TO 40

D(1) = E(1)

E(1) = 0.0E0

E(N) = 0.0E0

DO30K=1,NM1

KP1=K+1

C FIND THE LARGEST OF THE TWO ROWS

IF (ABS(C(KP1)) .LT. ABS(C(K))) GO TO 10

C INTERCHANGE ROW

T = C(KP1)

C(KP1) = C(K)

C(K) = T

T = D(KP1)

D(KP1) = D(K)

D(K) = T

T = E(KP1)

E(KP1) = E(K)

E(K) = T

T = B(KP1)

B(KP1) = B(K)

B(K) = T

10 CONTINUE

C ZERO ELEMENTS

IF (C(K) .NE. 0.0E0) GO TO 20

INFO = K

GO TO 100

20 CONTINUE

Section 2.4 Program realization and computational examples 41

T = -C(KP1)/C(K)

C(KP1) = D(KP1) + T

D(K)

D(KP1) = E(KP1) + T

E(K)

E(KP1) = 0.0E0

B(KP1) = B(KP1) + T

B(K)

30 CONTINUE

40 CONTINUE

IF (C(N) .NE. 0.0E0) GO TO 50

INFO = N

GO TO 90

50 CONTINUE

C BACK SOLVE

NM2=N-2

B(N) = B(N)/C(N)

IF (N .EQ. 1) GO TO 80

B(NM1) = (B(NM1) - D(NM1)

B(N))/C(NM1)

IF (NM2 .LT. 1) GO TO 70

DO 60 KB = 1, NM2

K=NM2-KB+1

B(K) = (B(K) - D(K)

B(K+1) - E(K)

B(K+2))/C(K)

60 CONTINUE

70 CONTINUE

80 CONTINUE

90 CONTINUE

100 CONTINUE

RETURN

END

Below we present a program solving the difference convection-diffusion problem in

which the convective term in (2.64) is approximated in two ways.

Program PROBLEM1

C PROBLEM1 - ONE-DIMENSIONAL STATIONARY

CONVECTION-DIFFUSION PROBLEM

REAL KAPPA

PARAMETER ( KAPPA = 0.01, ISCHEME = 0,N=20)

DIMENSION Y(N+1), X(N+1), PHI(N+1)

+ ,ALPHA(N+1), BETA(N+1), GAMMA(N+1)

C PARAMETERS:

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

C KAPPA - DIFFUSIVITY

C ISCHEME - CENTRAL-DIFFERENCE APPROXIMATIONS (ISCHEME = 0),

C SCHEME WITH DIRECTED DIFFERENCES;

C N + 1 - NUMBER OF NODES;

C Y(N+1) - EXACT SOLUTION

42 Chapter 2 Boundary value problems for ordinary differential equations

XL=0.

XR=1.

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

THE COMPUTED DATA

C MESH

H = (XR - XL)/N

DOI=1,N+1

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

END DO

C BOUNDARY CONDITION AT THE LEFT END

GAMMA(1) = 1.

BETA(1) = 0.

PHI(1) = 0.

C BOUNDARY CONDITION AT THE RIGHT END

ALPHA(N+1) = 0.

GAMMA(N+1) = 1.

PHI(N+1) = 0.

C ELEMENTS OF THE TRIDIAGONAL MATRIX

IF (ISCHEME.EQ.0) THEN

C SCHEME WITH CENTRAL-DIFFERENCE APPROXIMATIONS

DO I = 2,N

ALPHA(I) = V(X(I))/(2.

H) + KAPPA/(H

BETA(I) = - V(X(I))/(2.

H) + KAPPA/(H

GAMMA(I) = ALPHA(I) + BETA(I)

PHI(I) = F(X(I))

END DO

ELSE

C SCHEME WITH DIRECTED DIFFERENCE APPROXIMATIONS

DO I = 2,N

VP = 0.5

(V(X(I)) + ABS(V(X(I))))

VM = 0.5

(V(X(I)) - ABS(V(X(I))))

ALPHA(I) = VP/H + KAPPA/(H

BETA(I) = - VM/H + KAPPA/(H

GAMMA(I) = ALPHA(I) + BETA(I)

PHI(I) = F(X(I))

END DO

END IF

C SOLUTION OF THE DIFFERENCE PROBLEM

DO I = 1,N+1

ALPHA(I) = - ALPHA(I)

BETA(I) = - BETA(I)

END DO

Section 2.4 Program realization and computational examples 43

CALL SGTSL (N+1, ALPHA, GAMMA, BETA, PHI, INFO)

IF (INFO.NE.0) STOP ! POOR MATRIX

C EXACT SOLUTION

AL=XR-XL

DO I = 1,N+1

Y(I) = X(I) - AL

EXP((X(I)-AL)/(2.

KAPPA))

+ SINH(X(I)/(2.

KAPPA)) / SINH(AL/(2.

KAPPA))

END DO

C APPROXIMATE-SOLUTION INACCURACY

EC = 0.

EL2=0.

DO I = 1,N+1

AA = ABS(PHI(I) - Y(I))

IF (AA.GT.EC) EC = AA

EL2=EL2+AA

END DO

EL2 = SQRT(EL2

WRITE (01,

) KAPPA

WRITE (01,

) N+1, EC, EL2

WRITE (01,

) PHI

CLOSE (01)

STOP

END

FUNCTION F(X)

C RIGHT-HAND SIDE OF THE EQUATION

F=1.

RETURN

END

FUNCTION V(X)

C CONVECTIVE-TRANSFER COEFFICIENT

V=1.

RETURN

END

2.4.4 Computational experiments

Let us illustrate convergence of the above difference schemes with computational data

obtained on a sequence of consecutively reﬁned grids. The calculations were per-

formed for problem (2.64), (2.65) with parameters given as in (2.67). The approximate

solution is compared with the exact solution in the mesh norms in L

(ω) and C(ω).

44 Chapter 2 Boundary value problems for ordinary differential equations

Table 2.1 shows data obtained for the convection-controlled problem (κ = 0.01)

solved by scheme (2.72).

N 10 20 40 80 160

C 0.1599 0.2036 0.1579 0.09219 0.05070

0.03631 0.03451 0.02339 0.01340 0.007199

Table 2.1 Accuracy of the scheme with directed differences

Theoretical conclusions about the convergence behavior of the schemes with the

grid step tending to zero ﬁnd conﬁrmation here. The scheme with directed differences

converges with the ﬁrst order but, in the case of interest, this asymptotic convergence

is achieved with sufﬁciently ﬁne grids (starting from the grid with N = 80). Figure 2.3

shows plots of the approximate solution.

Figure 2.3 Approximate solution (scheme with directed differences)

Analogous data obtained for scheme (2.70) are given in Table 2.2.

N 10 20 40 80 160

C 0.4353 0.1932 0.05574 0.01212 0.003003

0.1074 0.03056 0.007142 0.001677 0.0004062

Table 2.2 Accuracy of the scheme with central-difference approximations

Section 2.5 Exercises 45

The inaccuracy decreases approximately fourfold on the twofold reﬁned grid

(second-order convergence). It should be noted that, in this example, the condition

for monotonicity (2.71) becomes violated on grids with N = 10, 20, 40. Violation

of this condition results in a substantially distorted solution behavior (see Figure 2.4).

We may even say that, here, non-monotone schemes cannot be used.

Figure 2.4 Scheme with central-difference approximations

2.5 Exercises

Exercise 2.1 On the solutions of equation (2.1), approximate the third-kind boundary

conditions (2.4) accurate to the second order.

Exercise 2.2 On the non-uniform grid

¯ω ={x | x = x

= x

i−1

+ h

, i = 1, 2,...,N , x

= 0, x

= l}

construct and examine the difference scheme

−

i+1

+ h



i+1

− y

i+1

− a

− y

i−1



+ c

= ϕ

, i = 1, 2,...,N − 1,

= μ

, y

= μ

making it possible to solve approximately problem (2.1), (2.2).