Georgian Mathematical Journal 1(1994), No. 6, 575-586

COMPUTATION OF POISSON TYPE INTEGRALS

T. BUCHUKURI

Abstract. We consider problems occurring in computing the Poisson integral when the point at which the integral is evaluated approaches the ball surface. Techniques are proposed enabling one to improve computation effectiveness.

0. Let B^{+} be the ball from ℝ^{3} with center at the origin and radius ρ:

B^{+} ≡ {x ∈ ℝ^{3} | |x| < ρ}

and B^{−} be the unbounded domain:

B^{−} ≡ {x ∈ ℝ^{3} | |x| > ρ}, S ≡ ∂B^{+} = ∂B^{−} = {x ∈ ℝ^{3} | |x| = ρ}.

The solution of the Dirichlet problem for the Laplace equation in the domain B^{+} is expressed by the Poisson integral. The solutions of the Dirichlet problem in B^{−} and of the Neumann problem in B^{+} and B^{−} are expressed by integrals of the same kind. We shall call these expressions Poisson type integrals.

It is well known (see [2-5]) that in case of ball solutions of the basic boundary value problems of elasticity, thermoelasticity, elastic mixtures, fluid flow can also be expressed by simple combinations of Poisson type integrals. Such representations prove convenient in constructing numerical solutions. The latter solutions possess the advantages of the method of boundary integral equations, namely: they decrease the problem dimension by one and allow us to evaluate the solution at any point without using solution values at other points of the domain under consideration. The method can equally be used for the domains B^{+} and B^{−}; in both cases we must compute integrals extended over the surface S. The solutions are represented as combinations of Poisson type integrals whose kernels are expressed by means of elementary functions and densities are given boundary conditions.

575

1991 Mathematics Subject Classification. 65D30, 65D32.

1072-947X/1994/1100-0575$07.00/0 \begin{matrix}c \\ ◯ \end{matrix} 1994 Plenum Publishing Corporation

Methods of computing such integrals do not actually differ from those commonly used in computing two-dimensional integrals. Nevertheless they need a certain modification so that they could lead to effective algorithms for computing Poisson type integrals. In particular, difficulties arise when integrals are computed at points close to the surface S because at these points the kernels of integrals have strong singularities.

1. In computing integrals which are solutions of the considered boundary value problems we come across the same difficulty as in the case of the Poisson integral regarded as the simplest one among integrals of this kind.

Let u be the solution of the Dirichlet problem for the ball B^{+}:

∀x ∈ B^{+} : Δu(x) = 0, ∀y ∈ S : u^{+}(y) = f(y),

where f is a given function on S. Then u can be represented by the Poisson formula

u(ρ_{0}, ϑ_{0}, φ_{0}) = \frac{1 − τ^{2}}{4π} ∫_{0}^{π} ∫_{0}^{2π} \frac{sin ϑ}{(1 − 2τ cos γ + τ^{2})^{3/2}} f̃(ϑ, φ)dφdϑ. (1)

Here (ρ_{0}, ϑ_{0}, φ_{0}) and (ρ, ϑ, φ) are the spherical coordinates of the points x^{0} ∈ B^{+} and y ∈ S;

x_{1}^{0} = ρ_{0} cos φ_{0} sin ϑ_{0}, x_{2}^{0} = ρ_{0} sin φ_{0} sin ϑ_{0}, x_{3}^{0} = ρ_{0} cos ϑ_{0};
y_{1} = ρ cos φ sin ϑ, y_{2} = ρ sin φ sin ϑ, y_{3} = ρ cos ϑ,

γ is the angle between the vectors x^{0} and y,

τ = ρ_{0}/ρ , f̃(ϑ, φ) ≡ f(ρ cos φ sin ϑ, ρ sin φ sin φ, ρ cos ϑ).

To compute the integral (1) it is convenient to represent it as an iterated integral and to use any of the quadrature formulae for one-dimensional integrals (in our case this will be Simpson’s formula; see, for example, [6]).

2. In computing the integral

∫_{a}^{b} f(t) dt  (2)

by Simpson’s method, its value is approximately replaced by the sum

S(f, a, b, m) ≡ \frac{h}{3}(f(a) + f(b) + 2 ∑_{k=1}^{m−1} f(a + 2kh) +

+ 4 ∑_{k=0}^{m−1} f(a + (2k + 1)h),  (3)

where h ≡ \frac{b−a}{2m}.

Note that (3) contains the value of f at the 2m + 1 points: t = a, t = a + h, …, t = a + (2m + 1)h, t = b.

Denote the error of Simpson’s formula by R(h):

R(h) ≡ ∫_{a}^{b} f(t)dt − S(f, a, b, m).  (4)

If f ∈ C^{4}([a, b]), then [6]

R(h) = −\frac{(b − a)h^{4}}{180} f^{(4)}(ξ)  (5)

for some ξ ∈]a, b[. Thus to estimate the error we obtain

|R(h)| ≤ \frac{(b − a)h^{4}}{180} max_{a≤ξ≤b} |f^{(4)}(ξ)|.  (6)

The estimate (6) is rather crude and its application may lead to a substantial increase of m in the sum (3). This happens particularly when f^{(4)} sharply increases in the neighborhood of some point of [a, b]. It will be shown below that the same situation occurs in computing Poisson type integrals.

In practice, the error is frequently estimated using Runge’s principle [6] which is as follows: If the condition

σ ≡ |S(f, a, b, 2m) − S(f, a, b, m)| ≤ ε  (7)

is fulfilled for some m, then S(f, a, b, 2m) is taken as an approximate value of integral (2) and the number ε for the error. As has been established, for Simpson’s formula the error can be estimated by R(h) ≈ \frac{σ}{15}.

3. To compute the integral (1) the above-mentioned procedure of replacing the one-dimensional integral by the sum (3) should be applied twice. Let us estimate the error for either of the cases so that the computational error for the integral (3) be not greater that ε.

Denote by S̃(f, a, b, δ) sum (3) for m such that

| ∫_{a}^{b} f(t) dt − S(f, a, b, m)| ≤ δ.  (8)

From (6) it follows that for the fulfilment of (8) it is sufficient to take

m ≥ \frac{(b − a)^{5/4}}{2 · 180^{1/4}} ( max_{a≤ξ≤b} |f^{(4)}(ξ)|)^{1/4}δ^{−1/4}.  (9)

578 T. BUCHUKURI

Denote by F the function

F (ϑ, φ) ≡ K(τ, ϑ, φ)f̃(ϑ, φ),

K(τ, ϑ, φ) ≡ \frac{sin ϑ}{(1 − 2τ cos γ + τ^{2})^{3/2}}

and by I the integral

I(ϑ) ≡ ∫_{0}^{2π}F (ϑ, φ) dφ.

Now (1) can be rewritten as

u(ρ_{0}, ϑ_{0}, φ_{0}) = \frac{1 − τ^{2}}{4π} ∫_{0}^{π} I(ϑ) dϑ.  (1′)

Due to (8) and (9) we have

|I(ϑ) − Ĩ(ϑ)| ≤ δ_{1},  (10)

where

Ĩ(ϑ) ≡ S̃(F (ϑ, ·), 0, 2π, δ_{1}) = S(F (ϑ, ·), 0, 2π, m_{1}),
m_{1} ≥ \frac{(2π)^{5/4}}{2 · 180^{1/4}} ( max_{ϑ,φ} |\frac{∂^{4}F (ϑ, φ)}{∂φ^{4}} |)^{1/4}δ_{1}^{−1/4}.  (11)

Let ũ(ρ_{0}, ϑ_{0}, φ_{0}) be an approximate value of integral (1′):

ũ(ρ_{0}, ϑ_{0}, φ_{0}) ≡ \frac{1 − τ^{2}}{4π} S̃(Ĩ, 0, π, δ_{2}).

Then

|u(ρ_{0}, ϑ_{0}, φ_{0}) − ũ(ρ_{0}, ϑ_{0}, φ_{0})| ≤ \frac{1 − τ^{2}}{4π} | ∫_{0}^{π} I(ϑ) dϑ −

S̃(I, 0, π, δ_{2})|dϑ + \frac{1 − τ^{2}}{4π} |S̃(I, 0, π, δ_{2}) − S̃(Ĩ, 0, π, δ_{2})|,  (12)
S̃(I, 0, π, δ_{2}) ≡ S(I, 0, π, m_{2}),
m_{2} ≥ \frac{π^{3/4}}{2 · 180^{1/4}} ( max_{ϑ,φ} |\frac{∂^{4}F (ϑ, φ)}{∂ϑ^{4}} |)^{1/4}δ_{2}^{−1/4}.  (13)

The first term on the right-hand side of (12) is less than (1 − τ^{2})δ_{2}/(4π). Let us estimate the second term. Note that

|S(f, a, b, m)| ≤ (b − a) max_{a≤x≤b} |f(x)|.

Therefore

|S̃(I, 0, π, δ_{2}) − S̃(Ĩ, 0, π, δ_{2})| ≤ π max_{0≤ϑ≤π} |I(ϑ) − Ĩ(ϑ)| ≤ πδ_{1}.

Now (12) yields

|u(ρ_{0}, ϑ_{0}, φ_{0}) − ũ(ρ_{0}, ϑ_{0}, φ_{0})| ≤ \frac{1 − τ^{2}}{4π} (δ_{2} + πδ_{1}).

If

δ_{1} = δ_{2} = \frac{4πε}{(1 − τ^{2})(π + 1)},

we finally obtain

|u(ρ_{0}, ϑ_{0}, φ_{0}) − ũ(ρ_{0}, ϑ_{0}, φ_{0})| ≤ ε.

From (11) and (13) it follows that the number of nodes N required that a given accuracy ε be achieved is estimated as follows:

N = m_{1}m_{2} = c_{0}c(F)(1 − τ)^{1/2}ε^{−1/2}.  (14)

where c_{0} does not depend on F and

c(F) = ( max_{ϑ,φ} |\frac{∂^{4}F (ϑ, φ)}{∂φ^{4}} | · max_{ϑ,φ} |\frac{∂^{4}F (ϑ, φ)}{∂ϑ^{4}} |)^{1/4}.  (15)

Below c will denote an arbitrary positive constant.

4. In practice it is important to know the values of the Poisson integral at points close to the boundary. But when the point x^{0} = (ρ_{0}, ϑ_{0}, φ_{0}) approaches the boundary S = ∂B(0, ρ), in computing the values of the Poisson integral at the point x^{0} by the above-described method, we observe a sharp increase of the number of nodes at which the function f is evaluated, and, accordingly, nearly the same increase of the computational time. We shall show why this happens.

Let A ≡ (1 − 2τ cos γ + τ^{2})^{1/2}, then

A^{2} = (1 − τ)^{2} + 4τ sin^{2} \frac{γ}{2} > (1 − τ)^{2}

and 1/A < 1/d, where d ≡ 1 − τ. Moreover, if

sin \frac{γ}{2} < \frac{d}{2\sqrt{1 − d}},

then A^{2} < 2(1 − τ)^{2} and 1/A > 1/\sqrt{2} d. Hence we conclude that 1/A has order 1/d for small γ.

Let us now estimate fourth order derivatives of F.

Assume that

∀(ϑ, φ) ∈ [0, π] × [0, 2π] : |\frac{∂^{α+β}f̃(ϑ, φ)}{∂ϑ^{α}∂φ^{β}} | ≤ c, α + β ≤ 4.  (16)

Since γ is the angle between the vectors x^{0} and y, we have either |ϑ−ϑ_{0}| ≤ γ and |φ − φ_{0}| ≤ γ or 2π − γ < |φ − φ_{0}| < 2π. Therefore

| sin(φ − φ_{0})| ≤ \frac{d}{2\sqrt{1 − d}}, | sin(ϑ − ϑ_{0})| ≤ \frac{d}{2\sqrt{1 − d}}.

Hence we conclude that if 0 ≤ α ≤ 4, then

\frac{∂^{α}K(τ, ϑ, φ)}{∂φ^{α}} ∼ \frac{1}{d^{α+3}} , \frac{∂^{α}K(τ, ϑ, φ)}{∂ϑ^{α}} ∼ \frac{1}{d^{α+3}} ,

and therefore

|\frac{∂^{4}F (ϑ, φ)}{∂φ^{4}} | ∼ \frac{1}{d^{7}} , |\frac{∂^{4}F (ϑ, φ)}{∂ϑ^{4}} | ∼ \frac{1}{d^{7}}  (17)

for

sin \frac{γ}{2} ≤ \frac{d}{2\sqrt{1 − d}}.

From (14), (15) and (17) we obtain

Theorem 1. The number of nodes required that a given accuracy be achieved in formula (1) admits the estimate

N(d) = O(d^{−3}).  (18)

This theorem explains the phenomenon of an increasing number of nodes as the point x^{0} approaches the boundary S.

5. Now let us find how for a given accuracy we can decrease the number of nodes and, accordingly, the computational time when the point x^{0} is near the boundary. This can be accomplished in different ways. One of the ways of reducing nodes is the so-called method of separation of singularities. To realize this method note that the identity

\frac{1 − τ^{2}}{4π} ∫_{0}^{π}∫_{0}^{2π} \frac{sin ϑ}{(1 − 2τ cos γ + τ^{2})^{3/2}} dφ dϑ = 1  (19)

is obviously valid and therefore (1) can be rewritten as

u(ρ_{0}, ϑ_{0}, φ_{0}) = f̃(ϑ_{0}, φ_{0}) + \frac{1 − τ^{2}}{4π} ∫_{0}^{π}∫_{0}^{2π}F_{1}(ϑ, φ) dφ dϑ,  (20)

with

F_{1}(ϑ, φ) = \frac{sin ϑ}{(1 − 2τ cos γ + τ^{2})^{3/2}} (f̃(ϑ, φ) − f̃(ϑ_{0}, φ_{0})).

Lagrange’s theorem implies

|f̃(ϑ, φ) − f̃(ϑ_{0}, φ_{0})| ≤ cγ.

This makes it possible to improve the error estimate, since

|\frac{∂^{4}F_{1}(ϑ, φ)}{∂φ^{4}} | ∼ \frac{1}{d^{6}} , |\frac{∂^{4}F_{1}(ϑ, φ)}{∂ϑ^{4}} | ∼ \frac{1}{d^{6}}

for

sin \frac{γ}{2} < \frac{d}{2\sqrt{1 − d}},

and therefore

N(d) = O(d^{−5/2}).  (21)

A further improvement of the method can be effectively achieved by the rotation of the coordinate system.

6. Let x^{0} ∈ B^{+}\{0}. Denote by A_{x^{0}} the following transformation of the Cartesian coordinates y = (y_{1}, y_{2}, y_{3}).

(A_{x^{0}}(y))_{1} = (\frac{(x_{1}^{0})^{2}}{|x^{0}|(|x^{0}| + x_{3}^{0})} − 1)y_{1} + \frac{x_{1}^{0}x_{2}^{0}}{|x^{0}|(|x^{0}| + x_{3}^{0})} y_{2} + \frac{x_{1}^{0}}{|x|} y_{3},

(A_{x^{0}}(y))_{2} = \frac{x_{1}^{0}x_{2}^{0}}{|x^{0}|(|x^{0}| + x_{3}^{0})} y_{1} + (\frac{(x_{2}^{0})^{2}}{|x^{0}|(|x^{0}| + x_{3}^{0})} − 1)y_{2} + \frac{x_{2}^{0}}{|x|} y_{3},

(A_{x^{0}}(y))_{3} = \frac{x_{1}^{0}}{|x|} y_{1} + \frac{x_{2}^{0}}{|x|} y_{2} + \frac{x_{3}^{0}}{|x|} y_{3}

for x_{3}^{0} ≥ 0 and

(A_{x^{0}}(y))_{1} = (1 − \frac{(x_{1}^{0})^{2}}{|x^{0}|(|x^{0}| − x_{3}^{0})})y_{1} − \frac{x_{1}^{0}x_{2}^{0}}{|x^{0}|(|x^{0}| − x_{3}^{0})} y_{2} + \frac{x_{1}^{0}}{|x|} y_{3},

(A_{x^{0}}(y))_{2} = −\frac{x_{1}^{0}x_{2}^{0}}{|x^{0}|(|x^{0}| − x_{3}^{0})} y_{1} + (1 − \frac{(x_{2}^{0})^{2}}{|x^{0}|(|x^{0}| − x_{3}^{0})})y_{2} + \frac{x_{2}^{0}}{|x|} y_{3},

(A_{x^{0}}(y))_{3} = \frac{x_{1}^{0}}{|x|} y_{1} + \frac{x_{2}^{0}}{|x|} y_{2} + \frac{x_{3}^{0}}{|x|} y_{3}.

Now we have

Theorem 2. Let x^{0} ∈ B^{+}\{0}. Then u(x^{0}) = u(ρ_{0}, ϑ_{0}, φ_{0}) from (1) can be written as

u(ρ_{0}, ϑ_{0}, φ_{0}) = f̃(ϑ_{0}, φ_{0}) + \frac{1 − τ^{2}}{4π} ∫_{0}^{π} \frac{sinϑ}{(1 − 2τ cos γ + τ^{2})^{3/2}} ×

× ( ∫_{0}^{2π}g̃(ϑ, φ) dφ − 2πf̃(ϑ_{0}, φ_{0})) dϑ,  (22)

where

g(y) ≡ f(A_{x^{0}}(y)), g̃(ϑ, φ) = g(ρ cos φ sin ϑ, ρ sin φ sin ϑ, ρ cos ϑ).

Proof. As one can easy verify, for any fixed x^{0} ∈ B^{+}\{0} the transformation A_{x^{0}} is the rotation of the coordinate system transforming the point x^{0} = (0, 0, |x^{0}|) into the point x^{0}. Therefore v(x) = u(A_{x^{0}}(x)) is a harmonic function taking on S the boundary value g(y). Moreover, in terms of spherical coordinates we shall have

v(ρ_{0}, 0, 0) = u(ρ_{0}, ϑ_{0}, φ_{0}), g̃(0, 0) = f̃(ϑ_{0}, φ_{0}).

Now apply (20) to the function v at the point (ρ_{0}, 0, 0). Note that in the case under consideration γ = ϑ. Therefore the kernel K does not depend on φ and can be put outside the internal integral. Then (20) gives (22). ■

Now let us estimate the number of nodes required for the computation of u(ρ_{0}, ϑ_{0}, φ_{0}) by formula (22).

Theorem 3. The number of nodes required that a given accuracy be achieved in computing u(ρ_{0}, ϑ_{0}, φ_{0}) by (22) admits the estimate

N(d) = O(d^{−1}).  (23)

Proof. We introduce the notation

I_{1}(ϑ) ≡ ∫_{0}^{2π}g̃(ϑ, φ) dφ − 2πf̃(ϑ_{0}, φ_{0}),

Ĩ_{1}(ϑ) ≡ S̃(g̃(ϑ, ·), 0, 2π, δ_{1}) − 2πf̃(ϑ_{0}, φ_{0}),

ũ(ρ_{0}, ϑ_{0}, φ_{0}) = f̃(ϑ_{0}, φ_{0}) + \frac{1 − τ^{2}}{4π} S̃(K(τ, ·)Ĩ_{1}, 0, π, δ_{2}).

COMPUTATION OF POISSON TYPE INTEGRALS 583

Then

|u(ρ_{0}, ϑ_{0}, φ_{0})−ũ(ρ_{0}, ϑ_{0}, φ_{0})|≤ \frac{1 − τ^{2}}{4π} ∫_{0}^{π} K(τ, ϑ)|I_{1}(ϑ)−Ĩ_{1}(ϑ)|dϑ+

+\frac{1 − τ^{2}}{4π} |∫_{0}^{π} K(τ, ϑ)Ĩ_{1}(ϑ)dϑ − S̃(K(τ, ·)Ĩ_{1}, 0, π, δ_{2})| .

Since

\frac{1 − τ^{2}}{4π} ∫_{0}^{π} K(τ, ϑ) dϑ = \frac{1}{2π} ,

we obtain

|u(ρ_{0}, ϑ_{0}, φ_{0}) − ũ(ρ_{0}, ϑ_{0}, φ_{0})| ≤ \frac{δ_{1}}{2π} + \frac{1 − τ^{2}}{4π} δ_{2}.

Choosing δ_{1} = πε, δ_{2} = 2πd^{−1}ε, we have

|u(ρ_{0}, ϑ_{0}, φ_{0}) − ũ(ρ_{0}, ϑ_{0}, φ_{0})| ≤ ε.

Let

S̃(g̃(ϑ, ·), 0, 2π, δ_{1}) = S(g̃(ϑ, ·), 0, 2π, m_{1}),
S̃(K(τ, ·), Ĩ_{1}, 0, π, δ_{2}) = S(K(τ, ·), Ĩ_{1}, 0, π, m_{2}),

Then by (11) and (16) m ≥ cε^{−1/4}.

Now we shall estimate m_{2}. Taking into account (9), we obtain

m_{2} ≥ c ( max_{0≤ϑ≤π} |\frac{∂^{4}}{∂ϑ^{4}} (K(τ, ϑ)Ĩ_{1}(ϑ))|)^{1/4} d^{1/4}ε^{−1/4},

Let us show that |Ĩ_{1}(ϑ)| ≤ cϑ. Indeed, f̃(ϑ_{0}, φ_{0}) = g̃(0, φ_{0}) and therefore

|Ĩ_{1}(ϑ)| = |S̃(g̃(ϑ, ·) − g̃(0, ·), 0, 2π, δ_{1})| ≤ |S̃(cϑ, 0, 2π, δ_{1})| ≤ c_{1}ϑ. With regard to this estimate we have

max_{0≤ϑ≤π} |\frac{∂^{4}}{∂ϑ^{4}} (K(τ, ϑ)Ĩ_{1}(ϑ))| ≤ cd^{−5}.

Therefore

N(d) = m_{1}m_{2} = O(d^{−1}). ■

7. A further improvement of the technique of computing the Poisson integral may be accomplished by giving up the uniform distribution of nodes on the integration interval. From the analysis of (22) it is obvious that when the condition (16) is fulfilled for not too large c, the approach of the point x^{0} to the boundary S does not affect in any essential way the computation of the internal integral. The computation of the external integral becomes, however, more difficult because the kernel

K(τ, ϑ) = \frac{sin ϑ}{(1 − 2τ cos ϑ + τ^{2})^{3/2}}

and its derivatives sharply increase as x^{0} approaches the boundary.

In case the computation is performed with a constant step (this implies the uniform distribution of nodes), there occurs a loss of accuracy on the part of the integration interval on which the derivatives of K sharply increase (see the estimates (6) and (7)), i.e. near the boundary. This means that we can improve the computational effectiveness by taking the lesser division step, the greater K and its derivatives are. We shall describe the technique realizing this idea.

Suppose we must compute the integral

∫_{a}^{b} f(t) dt  (24)

to within ε (ε > 0) when f sharply increases near the limit a.

A positive nonincreasing function δ determined on the segment [a, b] so that

∫_{a}^{b} δ(t) dt = 1

will be called a node distribution function.

Let, further,

Δ(x) ≡ ∫_{a}^{x} δ(t)dt.

Divide [a, b] by the points a = a_{0} < a_{1} < a_{2} < ⋯ < a_{n} = b into n parts such that the condition

| ∫_{a_{k}}^{a_{k+1}}f(t) dt − S(f, a_{k}, a_{k+1}, 2)| ≤ ε(Δ(a_{k+1}) − Δ(a_{k})),

where S is the sum determined by (3), be fulfilled on each part [a_{k}, a_{k+1}]. For this, by Runge’s principle it is sufficient that

|S(f, a_{k}, a_{k+1}, 1) − S(f, a_{k}, a_{k+1}, 2)| ≤ ε(Δ(a_{k+1}) − Δ(a_{k})). (25)

Then, denoting by S^{∗} (f, a, b, ε) the sum

S^{∗} (f, a, b, ε) ≡ ∑_{k=0}^{n−1} S(f, a_{k}, a_{k} + 1, 2),  (26)

we will have

| ∫_{a}^{b} f(t) dt− S^{∗} (f, a, b, ε)| ≤ ∑_{k=0}^{n−1} | ∫_{a_{k}}^{b_{k}}f(t) dt − S(f, a_{k}, a_{k+1}, 2)| ≤

≤ ε ∑_{k=0}^{n−1} (Δ(a_{k+1}) − Δ(a_{k})) = ε ∫_{a}^{b} δ(t) dt = ε.

Therefore S^{∗} (f, a, b, ε) is the desired approximation of the integral (24). Let the points a_{0}, ⋯, a_{n} be chosen such that for any k = 0, ⋯, n:

|S(f, a_{k}, a_{k+1}, 1) − S(f, a_{k}, a_{k+1}, 2)| ∼ ε(Δ(a_{k+1}) − Δ(a_{k})).

Then the estimate (5) and the equality

Δ(a_{k+1}) − Δ(a_{k}) = (a_{k+1} − a_{k})δ(ξ_{k}), a_{k} < ξ_{k} < a_{k+1},

imply

\frac{1}{15} \frac{(a_{k+1} − a_{k})^{5}}{180 · 4^{4}} |f^{(4)}(η_{k})| ∼ ε(a_{k+1} − a_{k})δ(ξ_{k}),

Hence

(a_{k+1} − a_{k})^{4} ∼ \frac{180 · 4^{4}εδ(ξ_{k})}{15|f^{(4)}(η_{k})|}.  (27)

The relation (27) provides us with the criteria for choosing points a_{k}. The integration step h_{k} = (a_{k+1} − a_{k})/4 depends on the values of δ and f^{(4)} on a given segment [a_{k}, a_{k+1}] of the integration interval. Therefore we can choose points a_{k} by induction so that (25) be fulfilled.

By an appropriate choice of δ we can obtain various degrees of dependence of the step on the function f. In particular, if

δ(x) ≡ \frac{f^{(4)}(x)}{∫_{a}^{b} f^{(4)}(t) dt},

then the actual step h_{k} will be nearly independent of f, i.e. we shall have the integration ”with an almost constant step”.

Consider the simplest case when δ is a constant function on [a, b]. Then

δ(x) = \frac{1}{b − a}, Δ(x) = \frac{x − a}{b − a}  (28)

and condition (25) becomes

|S(f, a_{k}, a_{k+1}, 1) − S(f, a_{k}, a_{k+1}, 2)| ≤ \frac{ε(a_{k+1} − a_{k})}{b − a}.

In that case

h_{k} = O ( \frac{1}{\sqrt[4]{|f^{(4)}(a_{k})|}}) .

Thus by an appropriate choice of the node distribution function we can improve the computational algorithm. In practice, a noticeable effect can be achieved even in the simplest case (28).

References

1. S.G. Mihlin, Linear partial differential equations. (Russian) Vysshaya Shkola, Moscow, 1977.

2. R.K. Chichinadze, Boundary value problems for the Stokes equations in the case of a domain with a spherical boundary. (Russian) Trudy Tbiliss. Mat. Inst. Razmadze 87(1987), 136-149.

3. ”—-”, On boundary value problems of the elasticity theory for the ball. (Russian) Trudy Tbiliss. Mat. Inst. Razmadze 96(1991), 97-122.

4. ”—-”, The solution of the third and fourth boundary value problems of the elasticity theory for the ball. (Russian) Bull. Acad. Sci. Georgia 145(1992), No. 2, 493-498.

5. R.K. Chichinadze and T.G. Gegelia, Solution in quadratures of the problems of thermoelasticity for the sphere and spherical cavity. Z. Anal. Anwendugen 8(1989), No. 6, 515-553.

6. L. Collatz, Funktionalanalysis und numerische Mathematik. Springer-Verlag, Berlin, 1964.

(Received 23.07.1993)

Author’s address:

A. Razmadze Mathematical Institute

Georgian Academy of Sciences

1, Z. Rukhadze St., Tbilisi, 380093

Republic of Georgia