Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics
Подождите немного. Документ загружается.
46 Chapter 2 Boundary value problems for ordinary differential equations
Exercise 2.3 Approximate the first-kind boundary conditions (2.4) accurate to the sec-
ond order so that to solve equation (2.1) on the half-step extended grid
x
i
= x
0
+ ih, i = 0, 1,...,N, x
0
=−h/2, x
N
= l + h/2.
Exercise 2.4 Show that the difference scheme
−(ay
¯x
)
x
= ϕ(x), x ∈ ω,
where
a
i
=
1
h
x
i
x
i−1
dx
k(x)
−1
,
ϕ
i
=
1
h
2
a
i+1
x
i+1
x
i
dx
k(x)
x
x
i
f (s) ds + a
i
x
i
x
i−1
dx
k(x)
x
i
x
f (s) ds
for the approximate solution of the equation
−
d
dx
k(x)
du
dx
= f (x), 0 < x < l
with boundary conditions (2.2) is a perfectly accurate scheme.
Exercise 2.5 Construct an absolutely monotone difference scheme of second approx-
imation order for the boundary value problem (2.64), (2.65) by passing to the equation
−
d
dx
˜
k(x)
du
dx
+˜q(x)u =
˜
f (x), 0 < x < l.
Exercise 2.6 For the boundary value problem (2.1) with homogeneous boundary con-
ditions
u(0) = 0, u(l) = 0,
construct a finite element scheme with solution representation in the form (2.22) based
on minimization of the functional (the Ritz method)
J (v) =
1
2
l
0
k(x)
dv
dx
2
+ q(x)v
2
(x)
dx −
l
0
f (x)v(x) dx.
Exercise 2.7 Prove the difference Friedrichs inequality (Lemma 2.1) in the case of
even N .
Exercise 2.8 Prove the difference Friedrichs inequality for mesh functions given on
the non-uniform grid
¯ω ={x | x = x
i
= x
i−1
+ h
i
, i = 1, 2,...,N, x
0
= 0, x
N
= l},
where h
i
> 0, i = 1, 2,...,N .
Section 2.5 Exercises 47
Exercise 2.9 Prove that for any mesh function y(x) given on the uniform grid ¯ω and
vanishing at x = 0 and x = l there holds the following inequality (embedding theorem
for mesh functions):
y(x)
∞
≤
√
l
2
y
x
+
.
Exercise 2.10 Find eigenfunctions and eigenvalues of the difference problem
y
xx
+ λy = 0, x ∈ ω,
y
0
= 0, y
N
= 0.
Exercise 2.11 Suppose that in solving the boundary value problem
−
d
dx
k(x)
du
dx
= 0,
u(0) = 1, u(1) = 0
one uses the difference scheme
−k(x)y
xx
− k
◦
x
y
◦
x
= 0.
Show that this scheme diverges in the class of discontinuous coefficients. Perform
numerical experiments.
Exercise 2.12 Let the following conditions be fulfilled:
|α
i
| > 0, |β
i
| > 0,δ
i
=|γ
i
|−|α
i
|−|β
i
| > 0,
i = 1, 2,...,N − 1.
Prove that under such conditions for the solution of the problem
−α
i
y
i−1
+ γ
i
y
i
− β
i
y
i+1
= ϕ
i
, i = 1, 2,...,N − 1,
y
0
= 0, y
N
= 0
there holds the estimate
y
∞
≤ϕ/δ
∞
.
Exercise 2.13 Let in the difference scheme (2.51), (2.52), (2.68) μ
1
≥ 0, μ
2
≥ 0
and ϕ
i
≥ 0 for all i = 1, 2,...,N − 1 (or μ
1
≤ 0, μ
2
≤ 0 and ϕ
i
≤ 0 for i =
1, 2,...,N − 1). Then, under the conditions
γ
i
≥ α
i+1
+ β
i−1
, i = 2, 3,...,N − 2,
γ
1
>α
2
,γ
N −1
>β
N −2
(see Theorem 2.5) we have y
i
≥ 0, i = 1, 2,...,N −1(y
i
≤ 0, i = 1, 2,...,N −1).
48 Chapter 2 Boundary value problems for ordinary differential equations
Exercise 2.14 Construct a sweep method for solving the system of linear equations
−α
0
y
N
+ γ
0
y
0
− β
0
y
1
= ϕ
0
,
−α
i
y
i−1
+ γ
i
y
i
− β
i
y
i+1
= ϕ
i
, i = 1, 2,...,N − 1.
−α
N
y
N −1
+ γ
N
y
N
− β
N
y
0
= ϕ
N
that arises, for instance, when one solves the boundary value problem for equation
(2.1) with periodic boundary conditions.
Exercise 2.15 To numerically solve the problem (2.64), (2.65), (2.67), construct a
difference scheme on a piecewise-constant grid using uniform grids on the segments
[0, x
∗
] and [x
∗
, l] with grid densening in the vicinity of x = l. Write a program
and perform computational experiments to examine the rate of convergence of the
difference scheme.
3 Boundary value problems for
elliptic equations
Among stationary mathematical physics problems, problems most important for appli-
cations are boundary value problems for second-order elliptic equations. For a model
two-dimensional problem, we consider the matter of construction of its discrete ana-
logues with the use of regular rectangular grids on the basis of finite-difference ap-
proximations. The convergence of approximate solution to the exact solution is exam-
ined in mesh Hilbert spaces using the properties of positive definiteness of the elliptic
difference operator. The convergence in the homogeneous norm is considered based
on the maximum principle for elliptic difference equations. In solving the difference
equations that arise upon discretization of boundary value problems for elliptic equa-
tions, most frequently used are iteration methods. Here, alternate-triangular iteration
methods deserve a detailed consideration. We present a program that solves the Dirich-
let problem for the second-order elliptic equation with variable coefficients and give
examples of performed calculations.
3.1 The difference elliptic problem
We consider the matter of approximation for model boundary value problems involv-
ing second-order elliptic equations. The present approach is primarily based on using
the integro-interpolation method (balance method). It is shown that it is possible to
construct an appropriate difference problem on the basis of finite element approxima-
tion. We briefly discuss the issue of finite-difference approximation of boundary value
problems in irregular domains.
3.1.1 Boundary value problems
In the present consideration, primary attention will be paid to two-dimensional bound-
ary value problems in which the calculation domain has a simplest form:
=
x | x = (x
1
, x
2
), 0 < x
α
< l
α
,α= 1, 2
,
i.e., is shaped as a rectangle. The main object in the present consideration is the
second-order elliptic equation
−
2
α=1
∂
∂x
α
k(x)
∂u
∂x
α
+ q(x)u = f (x), x ∈ . (3.1)
The following constraints are imposed on the equation coefficients:
k(x) ≥ κ>0, q(x) ≥ 0, x ∈ .
50 Chapter 3 Boundary value problems for elliptic equations
A typical simplest example of the second-order elliptic equation is given by the
Poisson equation:
−u ≡−
2
α=1
∂
2
u
∂x
2
α
= f (x), x ∈ . (3.2)
In the latter case, in (3.1) we have k(x) = 1 and q(x) = 0.
Equation (3.1) is to be supplemented with some boundary conditions. In the case of
the Dirichlet problem, the boundary conditions are
u(x) = μ(x), x ∈ ∂. (3.3)
In more complex cases, second- or third-order boundary conditions can be given on
the domain boundary or on a part of the boundary, for instance,
k(x)
∂u
∂n
+ σ(x)u = μ(x), x ∈ ∂, (3.4)
where, recall, n is the external normal to .
3.1.2 Difference problem
We use a grid uniform in both directions. For grids over particular directions x
α
,
α = 1, 2 we use the notation
ω
α
=
x
α
| x
α
= i
α
h
α
, i
α
= 0, 1,...,N
α
, N
α
h
α
= l
α
.
Here,
ω
α
=
x
α
| x
α
= i
α
h
α
, i
α
= 1, 2,...,N
α
− 1, N
α
h
α
= l
α
,
ω
+
α
=
x
α
| x
α
= i
α
h
α
, i
α
= 1, 2,...,N
α
, N
α
h
α
= l
α
.
For a grid in the rectangle ,weput
ω = ω
1
× ω
2
=
x | x = (x
1
, x
2
), x
α
∈ ω
α
,α= 1, 2
,
ω = ω
1
× ω
2
.
If the coefficients in (2.1) are smooth, then the difference scheme can be constructed
based on immediate passage from differential to difference operators. Similarly to
the one-dimensional case, for the boundary value problem (3.1), (3.3), we put into
correspondence to the differential equation the difference equation
Ay = ϕ(x), x ∈ ω, (3.5)
with
A =
2
α=1
A
(α)
, x ∈ ω, (3.6)
Section 3.1 The difference elliptic problem 51
Figure 3.1 Control volume
A
(α)
y =−(a
(α)
y
¯x
α
)
x
α
+ θ
α
c(x)y,α= 1, 2, x ∈ ω,
where θ
1
+ θ
2
= 1.
For the coefficients at higher-order derivatives we put
a
(1)
(x) = k(x
1
− 0.5h
1
, x
2
), x
1
∈ ω
+
1
, x
2
∈ ω
2
,
a
(2)
(x) = k(x
1
, x
2
− 0.5h
2
), x
1
∈ ω
1
, x
2
∈ ω
+
2
.
For the lowest coefficient and for the right-hand side in (3.5), (3.6) we have:
c(x) = q(x), ϕ(x) = f (x), x ∈ ω.
To approximate the boundary conditions (3.3) at the boundary nodes of ∂ω (
ω =
ω ∪ ∂ω), we use
y(x) = μ(x), x ∈ ∂ω. (3.7)
The most efficient approach to the construction of difference schemes implies using
the integro-interpolation method. First of all, the latter is the case with irregular grids
(triangular ones, for instance), but also with simplest rectangular grids. In the case of
the uniform rectangular grid under consideration, equation (3.1) is integrated over a
control volume around each inner node x of ω (see Figure 3.1):
x
=
s | s = (s
1
, s
2
), x
1
− 0.5h
1
≤ s
1
≤ x
1
+ 0.5h
1
,
x
2
− 0.5h
2
≤ s
2
≤ x
2
+ 0.5h
2
.
Like in the one-dimensional case, we arrive at the difference scheme (3.5), (3.6) with
a
(1)
(x) =
1
h
2
x
2
+0.5h
2
x
2
−0.5h
2
1
h
1
x
1
x
1
−h
1
ds
1
k(s)
−1
ds
2
,
52 Chapter 3 Boundary value problems for elliptic equations
a
(2)
(x) =
1
h
1
x
1
+0.5h
1
x
1
−0.5h
1
1
h
2
x
2
x
2
−h
2
ds
2
k(s)
−1
ds
1
.
The right-hand side and the lower terms are approximated with the expressions
ϕ(x) =
1
h
1
h
2
x
1
+0.5h
1
x
1
−0.5h
1
x
2
+0.5h
2
x
2
−0.5h
2
f (x) dx,
c(x) =
1
h
1
h
2
x
1
+0.5h
1
x
1
−0.5h
1
x
2
+0.5h
2
x
2
−0.5h
2
q(x) dx.
In a similar manner, difference schemes for irregular grids can be constructed.
3.1.3 Problems in irregular domains
Certain difficulties arise in numerical solution of boundary value problems for elliptic
equations in complex irregular calculation domains. So far we have dealt with prob-
lems in a rectangular region of a (regular calculation domain). We will not discuss
here the latter scope of rather complex problems at much length; instead, we only give
a short summary of major lines in this field.
Traditionally, in the approximate solution of stationary mathematical physics prob-
lems irregular calculation grids are widely used. Among irregular grids, two main
classes of grids can be distinguished.
Structured grids. A most important example of such grids are irregular quadrangular
grids, which in many respects inherit the properties of standard rectangular grids
or, in other words, present grids topologically equivalent to rectangular grids.
Unstructured grids. Here, the mesh pattern has a variable structure. It is impossible
to relate the calculation grid with some regular rectangular grid. In particular,
the scheme can be written at each point with different numbers of neighbors.
Approximation on structured grids can be constructed based on the noted closeness
of such grids to standard rectangular grids. The latter can be most easily done by
introducing new independent variables.
A second possibility bears no relation with formal introduction of new coordinates
and can be realized through approximation of the initial problem on such an irregular
grid. Of course, in the case of irregular grids the use of simplest approaches for the
construction of difference schemes on the basis of undetermined coefficients, although
possible, seems not to be a structurally reasonable strategy. In the latter case, one can
use the balance method. In both cases, i.e., in the case of general unstructured grids
and in the case of structured grids, one can construct difference schemes using finite
element approximations.
In a number of cases, for irregular grids it is possible to construct a matched grid
formed by the nodes of an ordinary non-uniform rectangular grid and by the boundary
Section 3.1 The difference elliptic problem 53
Figure 3.2 Method of fictitious domains
nodes. For the boundary to be formed by nodes, one has to use strongly non-uniform
grids. Problems related with construction of difference schemes on a matched grid can
be solved in the ordinary way.
A matched difference scheme can be constructed only for rather a narrow class of
calculation domains. That is why, normally, other approaches are used to solve the
problem related with calculation-domain irregularity. A simplest method here implies
using an ordinary rectangular grid in the calculation domain with a boundary condition
transferred to the node adjacent to the boundary. In fact, we deal here with the passage
from the initial problem (3.1), (3.3) to the problem in another domain whose boundary
is adapted to the grid (approximation of boundary).
The most natural and rather universal method here is an approach using a grid
formed by the nodes of a regular (uniform) grid (inner nodes) and additional irreg-
ular boundary nodes at the domain boundary. The latter nodes are formed by intersec-
tions of the lines drawn through the nodes of the regular grid and the boundary of the
calculation domain.
More convenient in the approximate solution of boundary value problems in irregu-
lar domains is the method of fictitious domains. This method is based on extension of
the initial calculation domain to some regular domain
0
( ⊂
0
), for instance, to
a rectangle in two-dimensional problems (Figure 3.2). Afterwards, the problem in
0
can be solved by ordinary difference methods. It is also necessary to so extend the
solution of the initial problem into the fictitious domain
1
=
0
\ that the differ-
ence solution of the problem in the extended domain
0
would give an approximate
solution in the initial domain .
The problem in the extended domain involves small (high) coefficients of the differ-
ential equation. As a result, it becomes necessary to examine the matter of accuracy
and computational realization of various iteration methods for solving resulting prob-
lems.
54 Chapter 3 Boundary value problems for elliptic equations
A method alternative to the method of fictitious domains in approximate solution
of boundary value problems in irregular calculation domains is the decomposition
method, in which partition of the calculation domain into simple subdomains is used.
This approach has been widely discussed in literature as applied to the development of
solution algorithms for boundary value problems using advanced parallel computers.
In each of the subdomains, individual boundary value problems are to be solved, and
the solutions are to be matched via boundary conditions. In the subdomains, special
grids can be used, matched or not with the grids in other subdomains. That is why
difference domain-decomposition schemes can be regarded on the difference level as
difference methods on composite grids.
3.2 Approximate-solution inaccuracy
Let us outline possible approaches for the examination of the rate of convergence of
approximate solution to the exact solution. The present consideration is based on the
examination of the properties of self-adjointness, positive definiteness, and monotonic-
ity of elliptic difference operators.
3.2.1 Elliptic difference operators
Note some basic properties of difference operators that arise in approximate solution of
the model problem (3.1), (3.3). Like in the previous consideration for one-dimensional
problems, let us reformulate (by modifying the right-hand side at boundary nodes) the
difference problem (3.5), (3.7) so that to obtain a problem with homogeneous boundary
conditions.
For mesh functions vanishing on the set of boundary nodes ∂ω, we define a Hilbert
space H = L
2
(ω) with the scalar product and the norm defined as
(y,w) ≡
x∈ω
y(x)w(x)h
1
h
2
, y≡(y, y)
1/2
.
As a most important property, we distinguish the property of self-adjointness of the
elliptic difference operator (3.5). The latter property of A stems, with regard to (3.6),
from self-adjointness of the one-dimensional operators A
(α)
, α = 1, 2. Taking the
latter into account, we obtain:
(Ay ,w)=
x
2
∈ω
2
h
2
x
1
∈ω
1
A
(1)
y(x)w(x)h
1
+
x
1
∈ω
1
h
1
x
2
∈ω
2
A
(2)
y(x)w(x)h
2
=
x
2
∈ω
2
h
2
x
1
∈ω
1
w(x)A
(1)
y(x)h
1
+
x
1
∈ω
1
h
1
x
2
∈ω
2
w(x)A
(2)
y(x)h
2
= (y, Aw).
Section 3.2 Approximate-solution inaccuracy 55
For two-dimensional mesh functions vanishing on ∂ω, we define the following dif-
ference analogue to the norm in W
1
2
(ω):
∇y
2
≡
x
1
∈ω
+
1
x
2
∈ω
2
(y
¯x
1
)
2
h
1
h
2
+
x
1
∈ω
1
x
2
∈ω
+
2
(y
¯x
2
)
2
h
1
h
2
.
Since
(Ay , y) = (A
(1)
y, y) + ( A
(2)
y, y)
=
x
2
∈ω
2
h
2
x
1
∈ω
+
1
a
(1)
(y
¯x
1
)
2
h
1
+
x
1
∈ω
1
h
1
x
2
∈ω
+
2
a
(2)
(y
¯x
2
)
2
h
2
and a
α
(x) ≥ κ, α = 1, 2, then for the two-dimensional difference operator (3.6) the
following inequality holds:
(Ay , y) ≥ κ∇y
2
. (3.8)
To estimate the two-dimensional difference operator of diffusion transport, we use
the Friedrichs inequality for two-dimensional mesh functions.
Lemma 3.1 For mesh functions y(x) vanishing on ∂ω the following inequality holds:
y
2
≤ M
0
∇y
2
, M
−1
0
=
8
l
2
1
+
8
l
2
2
. (3.9)
Proof. We take into account the Friedrichs inequality for one-dimensional mesh func-
tions (Lemma 2.1); then, we obtain the inequality
x
1
∈ω
+
1
x
2
∈ω
2
(y
¯x
1
)
2
h
1
h
2
+
x
1
∈ω
1
x
2
∈ω
+
2
(y
¯x
2
)
2
h
1
h
2
≥
8
l
2
1
+
8
l
2
2
x
1
∈ω
1
x
2
∈ω
2
y
2
h
1
h
2
.
The latter inequality yields (3.9).
From (3.8) and (3.9), the following lower estimate for A follows:
A ≥ κλ
0
E,λ
0
= M
−1
0
. (3.10)
Let us give now an upper estimate of the elliptic difference operator of interest.
Lemma 3.2 For the difference operator A there holds the inequality
A ≤ M
1
E (3.11)
with the constant
M
1
=
4
h
2
1
max
x∈ω
a
(1)
(x) + a
(1)
(x
1
+ h
1
, x
2
)
2
+
4
h
2
2
max
x∈ω
a
(2)
(x) + a
(2)
(x
1
, x
2
+ h
2
)
2
+ max
x∈ω
c(x).
Proof. This statement can be proved analogously to the one-dimensional case (see
Lemma 2.2).