An economic method of computing LPτ-sequences

252

1. A. Antonov and V. M. Saleev Figure 3 shows the spectrum for a sphere with ER = 0.1. The continuous

line shows the

result of a calculation

of a basic sphere with 1;R = 0.1 for A = 1, the dotted line shows the

result of a calculation

of a basic sphere with XR = 4 for h = 0.025, and the dashed line the

result of a calculation

of a basic sphere with CR = 10 for h = 0.01. If the dotted values may be

regarded as partially satisfactory

(for hc < mc2), then the dashed values are obviously false.

The table shows, for some values of k, the sampling values of the reIative probable errors Ik in percentages (for ,2’= 1000). It can be seen that if condition the probable errors do not increase (in comparison

(3) is satisfied, then for X % 1

with X = 1). They increase when condition

(3)

is not satisfied. Translated by J. Berry. REFERENCES I. SOBOL,1. M. ~umer~ca~,~o~te Carlo mer~uds (Chislenn~e metody Monte-Karfo), “Nauka”, Moscon,1973. 2.

MORTON, K. 8’. ScaIiq neutron tracks in Monte Carlo shielding caIcuIations. J. Nucl. Energy, 5, 3/4, 320-324,1957.

3.

BUSLENPO, N. P. et al. Tl7e method of stafisticaf trials (the Monte Carlo method) (Metod statisticheskikh ispytanii (I\?etod Monte-Karlo)). SMB. Fizmatgiz, Moscow, 1962.

4.

SOBOL, I. M. The Monte Carlo method for calculating the criticality in the multi~roup appro~~ation. in: The Monte Carlo method in the problem of radio&n transfer (Metod Monte-Karfo v probleme perenosa izluchenii), 234-254, Atomizdat, Moscow, 1967.

5.

KHISA>fUTDINOV, A. I. Reduction of the variance of a probability estimate by the Monte Carlo method. Zh. @hisI. h;lar. mat. FL, 10, 6, 1547-1549, 1970.

6.

POZDNYAKOV. L. A., SOBOL, I. M. and SWNYAEV. R. A. The effect of multiple Compton scattering on the spectrum of X-ray radiation. Calculations by the Monte Carlo method. Astron. zh.,

54,6.1246-1258.1977

U.S.S.R. Comput. Maths Math. PJ7ys. Vol. 19, pp. 252-256 $ Pergamon Press Ltd. 1980. Printed in Great Britain.

AN ECONOMIC METHOD OF COMPUTING LP7-SEQUENCES* I. A. ANTONOV and V. M. SALEEV Moscow (Received

18 Ju(I~ 1977; revised version 24 Jwle 1978)

A NEW signi~cantly faster algorithm for computing sequence (LPr-sequence) is presented.

points of I. M. Sobol’s unifo~iy

dist~buted

In the number-theoretic sense Sobol’s LP+equences [I] possess more of the best qualities of uniform distributiveness than any other sequence of points in the multidimensional cube.

lZh. @bhisl. Mat. mat. Fi:., 19, 1, 243-245,

1979.

Short communications

The use of points of an LP,-sequence distributed

in this cube guarantees

algorithms and a more uniform problems

253

instead of independent

random points uniformly

a higher accuracy of the calculations

scanning of the parameter

in some Monte Carlo

space in the solution of optimization

[2,3] .

In this paper we describe a new algorithm for computing

LP,-sequences,

considerably

faster than that proposed in [ 11. The order of succession of the points is changed, but the uniformity

characteristics

remain as before.

1. The algorithm

btQo,Ql.. construct

. . be an arbitrary n-dimensional

LP,-sequence.

From the same points we

a new sequence Qo’, Q1’. . . . , putting

where F (i) is the Gray code [3] ?corresponding

The operation

9.

to the number i:

denotes digital mod 2 addition

in the binary system (the instruction

“excluding OR” for the ES computer), and E(z) is the integral part of the number z. It will be proved below that PO’, Q1 ‘, . . . is also an LP,-sequence.

F(i-I)

The number F(i) is interesting because the binary. representations and F(i) differ in only one digit. whose number 1 is 2=1-log_

and is easily calculated

on various computers

[rii-:‘l

S:r(il]

(31

by a small number of operations.

described was composed for the Dnepr-21 computer,

r(l-1)

was calculated in two operations: the first operation = const. and the constant is known beforehand.

*r(i)

By [I] , if in the binary system i = e,,,

This notation

The new table of constants

of the directing points fill,

of directing points (directing

From (I), (3) and (4) it follows that (\I:. ,_-l, *Q,,

(4)

‘=g

‘,

of the point Q is calculated f12), . . by formula (4).

numbers)

is given in [5 ] .

=I‘(‘].

From this we obtain for calculating the coordinates simple recurrence formula: for all I
determined

>‘I ,,I,-f”,,

means that each of the Cartesian coordinates coordinates

of normalization

‘1, then

0 =r,l”‘i,,

from the corresponding

A program using

by which I for any value of

the algorithm the quantity/

of two adjacent numbers

il.:(“,

(qil’,

, sin’) of the point Qi’ a

254

I. A. Antonov and V. M. Saleev Herei=1,2,...

;fori=Oalltheqo~‘=...=4~,,‘=0.

Therefore, the proposed algorithm for computing

4ij’ consists of the following:

1) by formula (2) we compute r(i); 2) by formula (3) we compute i and find the directing numbers

Vj(l) in the table;

3) by formula (5) we compute 4ij’.

2. Properties of the sequence Q’o, Q'I , We recall that formula (2) defines a one-to-one

transformation

... of the natural series into

itself, any segment of the natural series of the form 0
into itself:

0 G r (i) G 2r-I.

Lemma In the transformation (2) the set of all binary segments of the natural series of length ? is transformed one-to-one into itself. hoof:

An arbitrary

binary segment of the natural series of length Zs (see [ 1 ] ) consists of

numbers i which in the binary system are written in the csAie,.

I=?,,

form (see [4] )

(6)

cl,

where the ck are fixed, and the ek assume the values 0 and 1 Since E(i/2) =

c.

C,+1e,

e2, then (2) implies that r(f)=;,,

where

.

c--j?<.

G,.

(7)

fork= 1,2,. ,s-1 and P,=c,+c,_~, i,=c:,+~,._, fork=s+l,. ... It is easy to verify, that in (7) all the Fk are fiied and the Fk assume the values

P~=ei~~~,_!

-I,and$ =c~. 0 and 1. Therefore the numbers

of (7) also form a binary segment of length Zs.

Moreover, a binary segment consisting of numbers of the form (6) is part of the segment 0 Q i < Y-1,

which passes into itself. This “large” segment consists of 2pwS binary segments

of length 2s, which after the transformation

(2) again constitute

the same segment. It is obvious

that the images of these 2”-$ segments are distinct and that the mapping is one-to-one.

73reorem 1

If 009 Pl, . . . is xr LP,-sequence,

then Q,‘, Q1 ‘, . . . is also an LPr-sequence.

Proof: By [I], a :equence of points is said to be an LP,-sequence if each binary segment of it containing not les: Lhan 2r+l points is a &-mesh, and the binary segments of the sequence Qo’r Ql’, . . . are simultaneously binary segments of the sequence Qo, Q1.

255

Short communications

7lleorem 2 . t possesses the properties

If Qo. Ql,

A or A’ of [5], then Qu’, Ql’, . . . also possesses

these properties. To prove this it is sufficient

to note that in the formulation of these properties of [5] only the set of binary segments of length 2n or 2*” occurs, and these properties are identical for both sequences.

3. Estimation In computing

the coordinates

vectors Qi’ for i = 1,2, . . . , A’ in floating

of the n-dimensional

point form on the Dnepr-21 computer and for the calculation

of the speed of response

it is necessary to use

h ‘I,.\ ) =.\‘(.1~- 131

operations,

of the vectors Qi

,=I

where m is the minimum

number

ratio of the times for computing

of binary digits necessary for writing down the number A’. The -2’numbers by the two algorithms is

!

-=-

which equals 3 approximately

li (.Y)

(.Y)

t’ (.I i

(8)

’

h”(.Yi

for N = 64 and increases with A:

The value of t/t’ was calculated experimentally Therefore the new algorithm is considerably

for ,1; = 0.5.106 and was found to equal 5.6.

more efficient than the one already known.

As a comparison we also note that the calculation of Nn pseudorandom numbers in floating point form by Lehmer’s algorithm on the Dnepr-2 1 computer requires 6Nn operations. which include ,Wr slow multiplication operations. In the realization of the new algorithm only fast operations are used. An analysis of the system of instructions of various computers has shown that the ratio of (8) is valid to within 30% for the majority of them (for example, the Minsk-32, BESM-4, BESM-6, and Ural-2 computers). The difference on the solution

in the speed of computing

time of some problems.

cube of functions

LP,-sequences

may have a considerable

effect

For example, multiple integrals over a unit Jr-dimensional

of the form

.f(r:

. . s.i=(l,;.l-

III

(l-;.,Ti)-z,

(9)

,==I

were calculated

on the Dnepr-2 1 computer.

For )I = 8 and .I’= 2 16 integral (9) was computed

the sequence Qi’ in 25 min, and using the sequence Qi in 2 hours 20 min. In conclusion

the authors sincerely thank I. M. Sobol for his help. Translated by J. Berry.

using

A. 1. Tkachenko

256

REFERENCES 1.

SOBOL. I. M. Multidimensional formy i funktsii Khaara),

quadrature forms and Haar functions “Nauka”, Moscow. 1969.

2.

VOROh’TSOV, probability

3.

SOBOL, I. hl. and STATNIKOV, R. B. The LP-search and problems of optimal design. In: Random problems (Probl. sluchainopo poiska), No. 1, 117- 135, “Zinatne”, Riga, 1972.

4.

BLRLlXASlP,

5.

SOBOL. I. 51. Uniformly, distributed Fiz., 16. 5, 1332-1337, 1976.

Yu. V. and POLLYAK, simulation of systems.

E. R. Algebraic

(Mnogomernye

kvadraturnye

Yu. G. On the use of quasirandom sequences in the direct Avtomatika i vychisl. tekhn., No. 6, 23-27, 1971.

coding rheoy

(Algebraicheskaya

sequences

with an additional

L;.S.S.R. Comput. hfaths. Math. Ph)*s. \‘ol. 19, pp. 256-259 C Pergamon Press Ltd. 1980. Printed in Great Britain.

teoriya

kodirovaniya),

uniform

property.

“Mu”,

search

Moscow,

1971

Zh. v%hisl. Mat. mat.

0041-5553/79/0201/0256507.50/0

SOME POSSIBILITIES FOR THE REDUCTION OF ROUNDING ERRORS IN THE INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS*

A.I.TKACHENKO Kiev

(Received

20 June 1977)

THE POSSIBILITY of compensating the basic rounding errors affecting the accuracy of the numerical integration of ordinary. differential equations, without the use of double digit?, is discussed. In the integration of ordinary differential equations in digital computers with a digital grid of short length the rounding errors may be the main source of computational errors. Usually only some of the rounding errors arising in the execution of individual operations at the last integration step, have a substantial effect on the resultant error introduced at a given step. By compensating for these errors the accuracy of the solution can be significantly increased without having recourse to calculations

with doubled digitry.

1. Let the Cauchy problem ~1’- f(x, JS), J(XO) = ~‘0, where fix, ~1) is a sufficiently smooth function of its arguments, be solved with step h < 1 in the fried point mode with II binary digits of the mantissa by the Runge-Kutta

method:

(1)

Here pi, qi, sil are constant coefficients (possibly rounded). The asterisk denotes multiplication with rounding; the index zero indicates rounded quantities. The error caused by rounding introduced at a step is estimated by the expression *Zh. yVchis1. Mat. mat. Fiz., 19, 1, 245-248,

1979.