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

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

26 Chapter 2 Boundary value problems for ordinary differential equations

integrate the relation

=−

k(x)

Q(x) over the segment x

i−1

≤ x ≤ x

i−1

− u



i−1

Q(x)

k(x)

dx ≈ Q

i−1/2



i−1

k(x)

We denote





i−1

k(x)



−1

. (2.27)

Then, we obtain

i−1/2

≈−a

¯x,i

, Q

i+1/2

≈−a

i+1

x,i

For the right-hand side, we have:



i+1/2

i−1/2

f (x) dx.

In view of the chosen completion, we put



i+1/2

i−1/2

q(x)u(x) dx ≈ c

where



i+1/2

i−1/2

q(x) dx.

Finally, in calculating the mesh function a

by formulas (2.27) (compare (2.17)) we

arrive at the difference scheme (2.16).

The above conservative schemes belong to the class of homogeneous difference

schemes (here, the coefﬁcients of the difference equation are calculated at all grid

nodes by one and the same formulas). In much the same manner to the balance

method, difference schemes for more general problems can be constructed, including

problems with boundary conditions (2.3), (2.4), problems speciﬁed on non-uniform

grids, multi-dimensional problems, etc. That is why the integro-interpolation method

is distinguished as a generic one in construction of discrete analogues of mathematical

physics problems.

2.2 Convergence of difference schemes

The key point in consideration of discrete analogues of boundary value problems is

examination of closeness of approximate solution to the exact solution. Based on

a priori estimates of solution stability for the difference problem, here we give an

estimate for the rate of convergence in mesh Hilbert spaces in numerical solution of

the model boundary value problem (2.1), (2.2).

Section 2.2 Convergence of difference schemes 27

2.2.1 Difference identities

Recall some notions in use in the theory of difference schemes. We assume that a

uniform grid

¯ω ={x | x = x

= ih, i = 0, 1,...,N, Nh = l},

is introduced over the segment [0, l]. Here, ω is the set of internal nodes:

ω ={x | x = x

= ih, i = 1, 2,...,N − 1, Nh = l}.

For other parts of ¯ω we use the following settings:

={x | x = x

= ih, i = 1, 2,...,N, Nh = l},

−

={x | x = x

= ih, i = 0, 1,...,N − 1, Nh = l}.

On the set of nodes ω and ω

we deﬁne the scalar products

(y,w)≡



x∈ω

y(x)w(x)h,

(y,w)

≡



x∈ω

y(x)w(x)h.

Let us also give difference analogues to the differentiation formula for the product

of two functions and to the integration-by-parts formula. With the previously intro-

duced deﬁnitions of the operators of the right and left difference derivatives, one can

immediately check the validity of the following equalities:

(yw)

¯x,i

= y

i−1

¯x,i

+ w

¯x,i

= y

¯x,i

+ w

i−1

¯x,i

(yw)

x,i

= y

i+1

x,i

+ w

x,i

= y

x,i

+ w

i+1

x,i

(2.28)

These equalities give the difference analogue of the differentiation formula

(yw) = y

+ w

The analogues for the integration-by-parts formula



w dx = y(b)w(b) − y(a)w(a) −



are the ﬁnite-difference identities

,w) =−(y,w

¯x

)

+ y

− y

¯x

,w) =−(y,w

)

−

+ y

N −1

− y

(2.29)

28 Chapter 2 Boundary value problems for ordinary differential equations

Replacing y

with a

¯x,i

in (2.29), we obtain the ﬁrst Green difference formula:

((ay

¯x

)

,w) =−(ay

¯x

)

+ a

¯x,N

− a

x,0

. (2.30)

The second Green difference formula is

((ay

¯x

)

,w)− (y,(aw

¯x

)

= a

¯x,N

− y

¯x,N

) − a

x,0

− y

x,0

). (2.31)

Formulas (2.30) and (2.31) take a simpler form,

((ay

¯x

)

,w) =−(ay

¯x

)

, (2.32)

((ay

¯x

)

,w) = (y,(aw

¯x

)

), (2.33)

for mesh functions y(x) and w(x) vanishing at x = 0 and x = l (on ∂ω).

2.2.2 Properties of the operator A

For the model problem (2.1), (2.2), we have constructed the difference scheme (2.20),

in which the difference operator A is deﬁned on the set of mesh functions y(x), x ∈¯ω

vanishing on ∂ω by (see (2.16))

Ay =−(ay

¯x

)

+ cy, x ∈ ω, (2.34)

where, for instance, a(x) = k(x − 0.5h), c(x) = q(x), x ∈ ω

. We introduce the

norm in the mesh Hilbert space H by the relation y=(y, y)

1/2

Like in the differential case, the difference operator A in H is a self-adjoint operator:

A = A

∗

. (2.35)

The equality (Ay,w) = (y, Aw) readily follows from (2.33).

With the usual constraints k(x) ≥ κ>0, q(x) ≥ 0, the following lower estimate

for A holds:

A ≥ κλ

E. (2.36)

Here λ

is the lowest eigenvalue of the difference operator of the second derivative. In

the case of a uniform grid, we have:

sin

πh

≥

Let us obtain such a lower estimate for A based on the following difference Friedrichs

inequality.

Lemma 2.1 For all mesh functions vanishing on ∂ω, the following inequality holds:

y

≤ M

(y

¯x



)

, w

≡ ((w, w)

)

1/2

, M

= l

/8. (2.37)

Section 2.2 Convergence of difference schemes 29

Proof. For such a mesh function y

= y(x

) we have:



k=1

¯x,i

h, (2.38)

=−



k=i+1

¯x,i

h. (2.39)

To estimate the right-hand sides of (2.38) and (2.39), we use the inequality







≤



We put a

= y

¯x,i

1/2

and b

= h

1/2

; then, from (2.38) and (2.39) we obtain:

≤ x



k=1

¯x,k

)

h, (2.40)

≤ (l − x

)



k=i+1

¯x,k

)

h. (2.41)

Let us consider n = (N −1)/2 if, for instance, N is an odd number. The case in which

N is even will be considered separately. From (2.40) and (2.41) we have:

≤ x



k=1

¯x,k

)

h, 1 ≤ i ≤ n,

≤ (l − x

)



k=n+1

¯x,k

)

h, n + 1 < i ≤ N .

We multiply each of the above inequalities by h and add the inequalities together:



i=1

h ≤



i=1



k=1

¯x,k

)

h +



i=n+1

(l − x

)



k=n+1

¯x,k

)

h. (2.42)

Taking into account that



i=1

h =

n+1



i=n+1

(l − x

)h <

from (2.42) we obtain the estimate (2.37).

30 Chapter 2 Boundary value problems for ordinary differential equations

The above proof can be extended, with no substantial changes, to the case of a

non-uniform grid. The same note applies to multi-dimensional problems provided that

rectangular or general irregular grids are used.

For the difference operator (2.34) with bounded coefﬁcients k(x) and q(x) as for

an operator in the ﬁnite-dimensional space H , the following upper estimate is useful.

Note that it is not valid for the differential operator A (unbounded operator).

Lemma 2.2 The estimate

(Ay , y) ≤ M

(y, y)(A ≤ M

E), (2.43)

with the constant

max

1≤i≤N −1

+ a

i+1

+ max

1≤i≤N −1

holds.

Proof. Indeed, we have:

(Ay , y) = (a(y

¯x

)

, 1) + (cy, y) =



i=1

− y

i−1

)



i=1

We use the inequality (a + b)

≤ 2a

+ 2b

, the condition for positiveness of the

mesh functions a(x), c(x), x ∈ ω, and the conditions y(x) = 0, x /∈ ω; then, we

obtain:

(Ay , y) ≤



i=1

+ a

i+1



i=1

This yield the desired estimate (2.43).

2.2.3 Accuracy of difference schemes

To approximately solve the problem (2.1), (2.2), we use the difference scheme

Ay = ϕ(x), x ∈ ω, (2.44)

= μ

, y

= μ

. (2.45)

To examine the accuracy of (2.44), consider the problem for the inaccuracy

z(x) = y(x) − u(x), x ∈¯ω.

For this inaccuracy, from (2.44) and (2.45) we obtain the difference problem

Az = ψ(x), x ∈ ω, (2.46)

Section 2.3 Solution of the difference problem 31

which, in view of the exact approximation of boundary conditions (2.2), is considered

on the set of mesh functions z(x) vanishing at the boundary. In (2.46) ψ(x) is the

approximation inaccuracy:

ψ(x) = ϕ(x) − Au , x ∈ ω. (2.47)

Considering the case of sufﬁciently smooth coefﬁcients and a sufﬁciently smooth

right-hand side of (2.1), in using (2.34) for the approximation inaccuracy we obtain:

ψ(x) = O(h

), x ∈ ω. (2.48)

Let us formulate now a simplest statement concerning the accuracy of the difference

scheme (2.44), (2.45) used to solve the model one-dimensional problem (2.1), (2.2).

Theorem 2.3 For the inaccuracy of the difference solution found by formulas (2.44),

(2.45), the following a priori estimate holds:

z

¯x



≤

1/2

ψ. (2.49)

Proof. We scalarwise multiply the equation (2.46) for the inaccuracy by z(x), x ∈ ω

to ﬁnd that

(Az , z) = (ψ, z). (2.50)

For the left-hand side of this equality, we have

(Az , z) ≥ κ(z

¯x



)

whereas the right-hand side can be estimated from below using the Friedrichs inequal-

ity (2.37):

(ψ, z) ≤ψz≤M

1/2

ψz

¯x



Substitution into (2.50) yields the desired estimate (2.49).

From the expression (2.48) for the local approximation inaccuracy and from the a

priori estimate (2.49), convergence of the difference solution to the exact solution with

the second-order accuracy follows.

2.3 Solution of the difference problem

On discretization of one-dimensional convection-diffusion problems, we arrive at

three-point difference problems. The difference solutions can be found using tradi-

tional direct linear algebra methods. Below, the sweep method (Thomas algorithm) is

outlined, presenting, as everybody knows, the classical Gauss algorithm for matrices

with special banded structure.

32 Chapter 2 Boundary value problems for ordinary differential equations

2.3.1 The sweep method

Consider, as a model problem, the problem (2.1), (2.2). On a uniform grid ¯ω with

a step size h, we put into correspondence to this differential problem for equation

(2.1) the difference problem (2.34), (2.44), (2.45). We write the difference problem as

follows:

−α

i−1

+ γ

− β

i+1

= ϕ

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

= μ

, y

= μ

. (2.52)

From (2.34) and (2.44) we obtain the following expressions for the coefﬁcients:

i−1/2

,β

i+1/2

i−1/2

+ k

i+1/2

) + c

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

To ﬁnd the solution of (2.51), (2.52) by the sweep method (Thomas algorithm),we

use the representation

= ξ

i+1

+ ϑ

i+1

, i = 0, 1,...,N − 1 (2.53)

with still undetermined coefﬁcients ξ

, ϑ

. Substitution of y

i−1

= ξ

+ϑ

into (2.51)

yields:

(γ

− α

− β

i+1

= ϕ

+ α

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

Now, we use representation (2.53) and obtain

((γ

− α

)ξ

i+1

− β

i+1

= ϕ

+ α

+ (γ

− α

)ϑ

i+1

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

This equality is fulﬁlled with arbitrary y

i+1

(γ

− α

)ξ

i+1

− β

= 0,

+ α

+ (γ

− α

)ϑ

i+1

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

From here, we obtain the following recurrence formulas for the sweep coefﬁcients

,ϑ

i+1

− α

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

i+1

+ α

− α

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

To start the calculations, we write the boundary condition (2.52) at the left end in the

form of (2.53), i.e. as y

= ξ

+ ϑ

, so that

= 0,ϑ

= μ

. (2.56)

After the sweep coefﬁcients are calculated by the recurrence formulas (2.54)–(2.56),

we can found, using formula (2.53) and the second boundary condition (2.52), the

solution itself.

Section 2.3 Solution of the difference problem 33

2.3.2 Correctness of the sweep algorithm

Let us formulate conditions sufﬁcient for the use of the above sweep-method formulas.

We do not consider here the whole scope of problems that arise in substantiation of

the sweep method. Here, we restrict ourselves just to the matter of correctness of

the method, which is equivalent in the case of interest to the requirement of nonzero

denominator in (2.54), (2.55).

Lemma 2.4 Let the following conditions be fulﬁlled for system (2.51), (2.52):

|α

| > 0, |β

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

|γ

|≥|α

|+|β

|, i = 1, 2,...,N − 1. (2.58)

Then, algorithm (2.53)–(2.56) is correct, i.e. in the formulas (2.54), (2.55) we have

− α

= 0.

Proof. We are going to show that

|ξ

|≤1, i = 1, 2,...,N − 1. (2.59)

In view of (2.57), (2.58), under such constraints on the sweep coefﬁcients we have:

|γ

− α

|≥||γ

|−|α

||ξ

|| ≥ ||γ

|−|α

|| ≥ |β

| > 0.

We will prove (2.59) by induction. For i = 1, inequality (2.59) holds. Suppose that

inequality (2.59) holds for i; then, by (2.54) we obtain for i + 1:

|ξ

i+1

|β

|γ

− α

≤

|β

≤ 1.

Provided that inequality (2.59) takes place, we also have: γ

− α

= 0.

For our model problem (2.1), (2.2), in using the difference scheme (2.51), (2.52)

conditions (2.57), (2.58) for sweep correctness will be fulﬁlled if

> 0,β

> 0,γ

≥ α

+ β

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

These conditions fulﬁlled, the maximum principle holds for the difference solution of

problem (2.51), (2.52), i.e., the difference scheme is monotone. This matter will be

discussed in more detail below.

34 Chapter 2 Boundary value problems for ordinary differential equations

2.3.3 The Gauss method

There exist a variety of versions of the sweep method more accurately taking into

account particular speciﬁc features of particular problems. For instance, one can ad-

here the classical direct solution method for systems of linear algebraic equations, the

Gauss method (LU-decomposition, compact scheme of the Gauss method).

In matrix form, the difference problem (2.51), (2.52) looks as

Ay = ϕ, x ∈ ω (2.60)

with a properly deﬁned right-hand side (see (2.20)). In view of the established proper-

ties of the operator A, the matrix A is a positive-deﬁnite matrix. In the latter case, the

major (corner) minors of the matrix are all positive and, hence, an LU-decomposition

A = LU takes place.

In the notation of (2.60), for the matrix A we have the following representation:

A =

⎡

⎢

⎣

−β

0 ... 0

−α

−β

... 0

0 −α

... 0

000... γ

N −1

⎤

⎥

⎦

In view of the banded structure of the initial matrix A, for the elements of the LU-

decomposition we can conveniently put

L =

⎡

⎢

⎣

−1

00... 0

−α

−1

0 ... 0

0 −α

−1

... 0

000... ξ

−1

⎤

⎥

⎦

, U =

⎡

⎢

⎣

1 −ξ

0 ... 0

01−ξ

... 0

00 1... 0

00 0... 1

⎤

⎥

⎦

The diagonal elements of L are given by the following recurrence relations:

i+1

− β

i−1

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

= 0. (2.61)

The solution of the equation Lϑ = ϕ is given by

i+1

+ α

− β

i−1

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

= 0. (2.62)

Now, from the system Uy = ϑ we obtain the solution of problem (2.60):

= ξ

i+1

+ ϑ

i+1

, i = 0, 1,...,N − 1, y

= 0. (2.63)

Section 2.4 Program realization and computational examples 35

Formulas (2.61)–(2.63) for the LU-decomposition stem from the above-described

sweep algorithm (2.53)–(2.56) on the substitution

−→ β

i−1

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

And vice versa (which is more correct), sweep formulas (2.53)–(2.56) follow from

the formulas of the Gauss-method compact scheme applied to the difference problem

(2.51), (2.52).

2.4 Program realization and computational examples

The matter of numerical solution of a boundary value problem for the convection-

diffusion equation is considered, the problem being a model one in continuum me-

chanics. Two types of difference schemes are constructed, with central differences

of second approximation order and directed differences of ﬁrst order used to approx-

imate the convective-transfer operator. A program that solves the model problem and

computational results obtained are presented.

2.4.1 Problem statement

In the theoretical consideration of approximate solution methods for boundary value

problems for ordinary differential equations, we have restricted ourselves to the case of

the simplest boundary value problem (2.1), (2.2). The matter of practical implementa-

tion of the numerical methods will be illustrated below for an example of a somewhat

more general problem for an equation of type (2.5).

Consider the boundary value Dirichlet problem for the convection-diffusion equa-

tion:

−



k(x)



+ v(x)

= f (x), 0 < x < l, (2.64)

u(0) = μ

, u(l) = μ

(2.65)

with k(x) ≥ κ>0.

In many practical problems, the prevailing contribution can be either due to the dif-

fusion term (the term with diffusivity k(x) in (2.64)) or due to convective transfer (the

term with velocity v(x)). The importance of convective transfer can be evaluated by

the Peclet number, that arises when one nondimensionalizes the convection-diffusion

equation:

Pe =

. (2.66)

Here v

is the characteristic velocity, and k

is the characteristic diffusivity. In the

equation of motion for a continuum medium, an analogous parameter is the Reynolds

number.