Open access peer-reviewed chapter

# Stability of Algorithms in Statistical Modeling

Written By

Alexander A. Kronberg and Tatiana K. Kronberg

Reviewed: 28 June 2022 Published: 04 October 2022

DOI: 10.5772/intechopen.106136

From the Edited Volume

## Applied Probability Theory - New Perspectives, Recent Advances and Trends

Edited by Abdo Abou Jaoudé

Chapter metrics overview

View Full Metrics

## Abstract

In this paper, we investigate algorithms stability for calculation of multidimensional integrals using the statistical modeling methods. We considered issues of the algorithms optimization and we give sufficient conditions for the stability. We apply our approach to both calculation of integral from the regression function and the moments integral calculation. In all our numerical experiences, we used the mt19937 pseudorandom number generator.

### Keywords

• statistical modeling
• pseudorandom numbers
• optimal density
• integral estimation
• Monte Carlo methods

## 1. Introduction

One of the main problem of the statistical modeling method (the Monte Carlo method) is the problem of quality for pseudorandom numbers. In the paper, we consider a task of multidimensional integrals calculation by the statistical modeling method and give sufficient conditions for the stability of this task to quality of pseudorandom numbers. Included results of various numerical experiences with the mt19937 pseudorandom number generator. In our work, we discuss important issues of algorithms optimization in the statistical modeling. In particular, we apply the new approach to the following: a task of finding of integral functionals from solution of boundary-value problems for both the linear [1] or nonlinear [2] elliptic equations (the estimations are given near to a boundary).

The paper is organized as follows: In Section 2, we give the sufficient conditions of stability. Calculation of an integral of very large dimensions is discussed in Section 3. Rare events effect is the subject of Section 4. In Section 5, we describe calculation of integral moments. In Section 6, we apply our approach to calculate an integral of the regression function. In Section 7, we give the conclusion of our studies.

## 2. Sufficient conditions of stability

Let

I=DfxdxE1

be the Riemann integral. Here D is a domain of the s-dimensional Euclidean space Rs. If the dimension s is large enough then we must use a statistical modeling method. In this case, our integral has form of the mathematical expectation for a random value η=fξ/pξ:

Dfxdx=Dpxfxpxdx=Efxpx=Eη.E2

Here px is a density of random variable ξD. We put px0 for fx0, and we say that there exists integral

Dfxdx.

A variance of the random value η:

σ2=varη=Eη2Eη2=Dpxfxpx2dxI2=Df2xpxdxI2.E3

We estimate the mathematical expectation Eη by the sum i=1Nηi/N, where ηi are independent realizations of the random value η. Suppose σ2 is finite, and N is large enough; then from the classical central limit theorem, it follows that the random value i=1Nηi/N has distribution close to the normal distribution with a mathematical expectation I, and mean-square deviation σ/N. This property above is useful to estimate error, e.g., using the 3σ rule. So we have

I1Ni=1Nηi3σNE4

with probability 0,997, approximately.

Let us ηi be realizable sampling values; then the value σ is estimated as the following:

σ21Ni=1Nηi21Ni=1Nηi2.E5

Suppose Eη4 is finite and in Eq. (3) we replace the σ by its approximate value. Then it changes the estimation of error in calculation of the integral in order of O1N.

In practice, when we simulate random variables ξ, we receive simulation with some density qx instead of simulation with the origin density px. Now, we investigate the stability of the theoretical estimation Eq. (2). Let us consider the following expression:

DpxfxpxdxDpxqxpxdx=DfxpxpxqxdxεDfxpxdx,E6

where ε=supxDpxqx. By If/p denote the integral fxpxdx. The both values ε and If/p provide the guaranteed proximity of the real estimation to the theoretical one of the integral.

Example 1. The inequality Eq. (6) is reduced to the equality if D=D1D2,D1D2=. For xD1 we get pxqxε>0, and fx0 in D2. If the condition If/p=+ holds, then the error of the real estimation of the integral will be infinity for any ε>0.

Hence, a quality of pseudorandom variables (i.e., smallness of ε) does not yet guaranties the smallness of the error in general, as the integral If/p have to be both finite and not very great in magnitude.

Suppose we simultaneously make the estimations for both If/p and the origin integral Eq. (2) using the same density px. Then we need to ask boundedness of the integral Dfxp2xdx to get the guaranteed stability of the estimation for the integral If/p, and so on. The qualitative comparison of simulation with both densities p1x and p2x can be provided not only by a magnitude of the variance estimation (here we do not pay attention to the complexity of random values simulation) but also magnitudes of both the integrals If/p1 and If/p2. On the other hand we have the Schwarz inequality:

DfxpxpxqxdxDpxqx2dxDf2xp2xdx.E7

The sufficient condition for the estimation to be stability is that Df2xp2x to be finite and not very great. The Schwarz inequality is reduced to the equality if and only if

λfxpx=pxqx,E8

where λ is a real number. From the above we get

qx=pxλfxpx,Dqxdx=DpxdxλDfxpxdx.E9

Therefore, the equality in Eq. (7) is reached under the necessary condition Dfxpxdx=0, when Dpxdx=Dqxdx=1.

Example 2. The condition above is realized, e.g., if

D=11,D1=10,D1=01,fx=1inD1,fx=1inD2,px=px.

Let us ηk be fξ/pξk. Now we consider an estimation:

Eηk=Dpxfxpxkdx=Dfkxpk1xdx,k1.E10

This expectation is actually estimated by the integral: Dqxfx/pxkdx:

Dfkxpk1xdxDqxfkxpkxdx=DpxqxfkxpkxdxεDfkxpkxdx.E11

The last integral is assumed to be a finite, and not very large. These conditions are desirable. In Eq. (11) the equality is reached like to the Example 1.

For all cases above, the stability will be observed if fx/pxM<+ for not very great M. From the Schwarz inequality we have:

Dpxqxfkxpkxdx=DpxqxpβxfkxpkβxdxDpxqx2p2βxdxDf2kxp2kβxdx,E12

where β is a real number. We get a family of proximity measures for the distribution densities:

Dpxqx2p2βxdx.E13

For β=0,5 we obtain expression

χ2pq=Dpxqx2p(x)dx.E14

that is well known in the mathematical statistics.

For β=0,5 we get Dpxqx2pxdx and have the obvious inequalities

Dpxqx2pxdxsupxDpxDpxqx2dx,Dpxqx2pxdxsupxDpxqx2Dpxdx=supxDpxqx2.

In Eq. (6) the equality is satisfied if and only if

pqpβ=λfkpkβ,pq=λfkpβpkβ,E15

i.e., the necessary condition is Dfkxpk2βxdx=0. This is realized in the Example 2.

Let us remark that for k=1 and β=1 we have

Dpxqxpxfxdx=DpqpfpDpq2pdxDf2pdx=χ2pqEη2.E16

We assume that integrals Dfx+εixdx are known and the subintegral functions fx+εix>0 close to a function fx0. Suppose also

1δ1εi<fxfx+εix<1+δ2εi,Jfδ3εiJf+εiJf+δ4εi,δj0,asεi0,j=1,2,3,4.Dqxfxpxdx=DqxfxJf+εifx+εixdx=Eη̂,E17

where the random value η̂ has a form:

η̂=fξ̂Jf+εifξ̂+εiξ̂.E18

Here ξ̂ is distributed with the density qx. Keeping the above factors in mind we get the following:

Jfδ31δ1η1+δ2Jf+δ4;

1. Regardless of qx, i.e., regardless of quality of a pseudorandom number generator we have Eη̂Jf,varη̂0, asε0;

2. All moments Eη̂k of the random variable are finite.

Let we calculate Dfxpxdx using the density

p1x=fxIf/ppx,

then the estimation variance equals to zero.

Suppose we calculate the integral Dfxdx with the density p1x. In this case it would be interesting to know both the values Dfxp1dx and Df2p12dx.

Proposition 1.If/p1=If/p.

Proposition 2.If2/p12=I2f/pIp2.

Now we consider the density

p2x=f2xpxIf2/p.

Using the density above for the estimation of the integral Df2pdx we obtain the estimation variance equals to zero.

Further we estimate the integral Dfxdx with the density p2: η̂=fξ2/pξ2, ξ2 is distributed with p2x. Suppose η=fξ/pξ, ξ is distributed with px; then If2/p2=If2/p.

Proposition 3.varη=varη̂.

Proposition 4.Dfxpxdx2volDDf2xp2xdx, where volD is the volume of the domain D.

Proposition 5. If volD=1, p3x=f2xIf2,η3=fξ3pξ3, where ξ3 is a random variable distributed with the density of p3x; η4=fξ4pξ4, where ξ4 is a random variable distributed with the density of px1, then varη3=varη4.

In actual practice normalization constants are usually unknown for both p1x and p2x. But using densities close to them we can get the approximate equalities in the Prepositions 1, 2, 3.

Now we consider

If=0110x1x2x10dx1dx10,E19

where the integration domain D=0110 is the 10-dimensional unit cube, the subintegral function fx is equals to x1x2x10. To realize algorithms of the statistical modeling at a computer it is necessary to set a number N of realizations for random variable η=fξ/pξ, where ξ is distributed with the density pξ. In fact we realize the discrete set of numbers ξi,i=1,,N, which we can consider to be realizations of some distribution qNx.

In all our numerical computations we use the pseudorandom number generator: generator type mt19937 [3]. For s=10,px1 we have

If=1/2109,7656104,If/p=If,If2/p=If2/p2=1/3101,6935105,σ=1/3101/2203,998103.

Table 1 shows the empirical estimations Î and σ̂ for I and σ, respectively. Taking pxi=3xi2 over each coordinate we get

NÎσ̂
1,000,0009,7851043,991103
9,000,0009,7681044,001103
81,000,0009,7661043,997103
100,000,0009,7651043,995103

### Table 1.

The results of numerical calculations for the integral Eq. (19) with the uniform density.

If2/p=1/310,If/p=,If2/p2=,σ=1/3101/220.

The value σ is the same as one for px1.

For N=100000000 the computer code outputs an error because of machine zero divide. The reason of this event is η=110(1/3α3i, where αi are the pseudorandom numbers with the uniform density in 01 [3]. If the formulas of random numbers simulation generate division by very small numbers then such formulas are one more source of the algorithms instability in the statistical modeling. As seen in Table 2 the value of integral If is successfully estimated, but the empirical estimations of σ̂ are sufficiently different from the theoretical value σ. This result is explained by the following: Eη4=,If/p=,If2/p2=. Taking pxi=2,6xi1,6 we obtain σ1,223103, If/p0,67556, the finite value of Eη4, and If2/p2=. Although the last estimation is infinite, but Table 3 shows that both values If and σ̂ are successfully calculated. The value of σ̂ is very close to σ.

NÎσ̂
1,000,0009,7661043,276103
9,000,0009,7541043,619103
81,000,0009,7651043,759103
100,000,000Infnan

### Table 2.

The results of numerical calculations for the integral Eq. (19) with the density pxi=3xi2.

NÎσ̂
1,000,0009,7691041,212103
9,000,0009,7601041,220103
81,000,0009,7661041,223103

### Table 3.

The results of numerical calculations for the integral Eq. (19) with the density pxi=2,6xi1,6.

## 3. Integrals of very large dimensions

We are coming now to the question of calculation of an integral

If=00ex1+x2++xsdx1dx2dx10=1E20

with the distribution density px=λseλx1+x2++xs. For λ2, the estimation variance η is infinity. For 0<λ<2, the variance will be finite. For λ>1, we obtain If/p= and If2/p2=. However, as seen in Table 4, the results of calculations for N=10000 allow us to make the conclusion below. If we have the pseudorandom generator of the high quality and a good px then we can calculate the very high dimensional integrals.

λsÎσ̂σ
1,0110000,9990,3200,324
1,0110,0000,9951,391,31
100520,00010040,8070,805
100340,0000,9940,6490,658
100180,0000,9910,6050,661

### Table 4.

The results of numerical calculations for the integral Eq. (20).

## 4. Special integrals

Let us consider the following class of the integrals:

0+sfx1x2xsdx1dx2dxs,E21

We define the behavior of the subintegral function as follows: fx1x2xs to be [label = ()]

1. close to 1 in the cube 0as as 0<a<1;

2. much less than unit as a<xi<b,ba;

3. equil to zero as bxi<.

Note that very often the integration of functions can be reduced to the linear combination of the integrals similar to Eq. (21) using various replacements of variables.

Below, let us perform a theoretical and numerical analysis how to integrate a model function from our class. The model function is assumed to be f1 as 0xia, otherwise f=0. We take both distribution densities set p1xi=λeλxi,0xi< and p2xi=ω+11xiω,0xi<1 to be examined. Our goal is to determine what of two densities provides the best accuracy of the integral computation with given model function.

If we simulate a random point ξ=ξ1ξs with densities p1xi=λeλxi then the integral estimation is given by

ηs=i=1sλeλξi,varηs=Eηs2Eηs2=1λ0aeλxdxsas=eλa1λ2sas.E22

Testing the variance varηs for the extremum over λ we get the minimum condition

λaeλa2eλa+2=0.E23

Let A be λa; then the equation Eq. (23) is reduced to

AeA2eA+2=0.E24

The equation above has the unique root at A1,593620. It follows that λmin=A/a. For such λmin the relative error with the 3σ rule is given by

3σNas=3eA1sA2s10,5/N.E25

Suppose N=9106,s=10,λ=A/a; then the theoretical value of the relative error is approximately 8,72103. The numerical estimation of the relative error is approximately 8,69103 as a010,001. Thus, the numerical estimation of one gives a good fit to the predicted value over a wide range of a.

Now we discuss the use of the density p2xi=ω+11xiω. First, we estimate the second moment of a random value ηs:

Eηs2=11a1ω1ω2s.E26

The parameter ω is chosen to be A/a; then the expression above is rewritten as follows

Eηs2=11A2/a211a1A/as.E27

Let us consider Eηs2 as a0:

lima0Eηs2=lima0a2a2A211a1A/as==lima0a2A211a1aA/as=lima0a2A211a1a1/aAs==lima0a2A211aeAs=lima0a2A21eAseA1A2a2s.E28

Comparison between Eqs. (22) and (28) allows to make the following conclusion. If w is chosen to be A/a then the asymptotics of variances, as a0, are the same in the densities set of p1xi,p2xi. In the numerical simulation the relative accuracy of 8,68103 is reached as N=9106,s=10,ω=A/a,a0010,001.

Let us turn now to the integral

If=0110fx1x2x10dx1dx2dx10,E29

where

fx1x2x10=1,0xi1/4,0,otherwise.E30

We put px1; then If=If/p=If2/p=1/4109,5367107. For N=640000 all realizations are turned out to be equal to zero, i.e., q640000x=0 as 0xi1/4. In this case, we have

0110pxfxpxdx0110q640000fxpxdx=1/4100=1/410.E31

In accordance with N, the realizations numbers of ηi are turned out to be equal to 1 as N=810000; equal to 3 as N=4000000; equal to 15 as N=16000000.

We now take pxi=ω+11xiω in the unit cube 0110ω>1. Such choice provides the gross realizations of points in 01/410 and as consequence, we get benefit in quality of random values (simultaneously, we have decrease of the estimation variance, and as consequence decrease of the statistical error with the 3σ rule.) Table 5 shows the calculations results for N=9000000. Note that σ̂ reaches the minimum as ω=5. In this case, we have p̂x=61x5,If/p̂=If2/p̂3,491011. Making the more detailed research for both ω=5 and the theoretical value σ5,83106 we get the results represented in Table 6. If the function fx1x2x10 close to some constant in 01/410 and small out of this interval then we can advise to use p̂=61x5 to calculate the integral in 0110.

wÎrule “3σ
39,531078,61109
49,581076,36109
59,551075,82109
69,551076,49109
79,511078,06109

### Table 5.

The results of numerical calculations for the integral Eq. (29) at various ω values.

NÎσ̂
640,0009,611076,15106
810,0009,561076,00106
1,000,0009,511075,91106
4,000,0009,521075,74106
16,000,0009,551075,84106

### Table 6.

The results of numerical calculations for the integral Eq. (29) at ω=5.

## 5. Moments calculation

We are now concerned with the following issue: to find the kth moments of a random value τ with the distribution density px:

Eτk=abxkpxdx.E32

In fact we have realizations of the random value ξ with a distribution density qx. With px replaced by qx in Eq. (32) we get an error

abxkpxdxabxkqxdx=abxkpxqxdx.E33

Suppose b= and ξmax are the maximum value of the random variable over the all realizations for fixed N; then value of ξmax gives shift ξmaxxkpxdx that increases both monotonically and without limit. The condition qx=0 as x>ξmax determines the lower limit of the last integral.

Many solutions of the boundary-value problems for the elliptic and parabolic Equations [4, 5] have a form of the expectations for the random value moments. Meaning of these expectations is the first exit time of the Wiener process trajectories to the domain boundary.

Let a domain be the three-dimensional ball with the radius r=1 and the Wiener trajectories start from the ball center; then a function of distribution of the first exit time for the Wiener trajectory is, in particular, given by [5].

Ft=1+2k=11kexpk2π2t/2,t0+.E34

From the above, we obtain the distribution density:

pt=2k=11k+1μk2expμk2t,μ=π2/2.E35

Assuming τ is distributed with this density and calculating the expectation of the kth moment we get

Eτk=0tkptdt.E36

In Table 7, we put the calculations results for N=1000000. The kth moment expectation can be represented in a form.

MomentSimulationTheory
13,3041013,333101
21,5531011,556101
39,8481029,841102
48,0761028,063102
58,2911028,193102
69,8431029,969102
71,3191011,414101
82,0701012,293101
94,5181014,182101
107,2861018,474101
119,0211024,251103
121,1831031,637104
138,3891036,634104

### Table 7.

The results of numerical calculations for the moments by the first way.

Eτk=0qxtkptqtdt,E37

where qt is some density in 0. Taking qt=λexpλt, for λ=π2/2 we get that the number of realizations ξi>1 will be almost twice as small as in the case of the modeling with the original pt. In this situation we should obtain degradation of the estimation for the high moments. The calculations results with qt for N=1000000 are represented in Table 8. However, in realizations at a computer we get the obvious improvement in quality of the moments estimation for all k from 5 to 20. Consider the choice of modeling strategy with regards to the variance. Suppose ξ and η be estimations of the statistical modeling for a value J, i.e., Eξ=Eη=J with the variances of σ12ξ,σ22η and the realizations of ξ=ξ1++ξN/N, η=η1++ηN/N. It would seem that for σ1ξ<σ2η the real estimation of ξ will be occurred close to the origin value of J. But this statement does not need to be always true. Without loss of generality it can believed that J=0. Additionally, if N is large enough then ξ, η are chosen be normal random variables with N0σ1 and N0σ2, respectively. The following theorem holds.

MomentSimulation
58,188102
61,013101
71,468101
82,247101
93,950101
108,283101
182,833103
191,056104
205,118104

### Table 8.

The results of numerical calculations for the moments by the second way.

Proposition 6. Let ξ,η be normal random variables, and ξN0σ1,ηN0σ2 then Pξ>η=2πarctanσ1σ2.

Proof:

P|ξ>|η|=1σ12π0ey22σ12dy1σ22π0yex22σ22dx+1σ22π0yex22σ22dx++1σ12π0ey22σ12dy1σ22π0yex22σ22dx+1σ22π0yex22σ22dx==2σ12π0ey22σ12dy2σ22π0yex22σ22dx=42πσ1σ20ey22σ12dy0yex22σ22dx.

Using Taylor expansion

ey22σ22=n=01ny2n2nσ22nn!,

we get

0yn=01nx2n2nσ22nn!dx=n=01ny2n+12n+12nσ22nn!,0n=01ny2n+12n+12nσ22nn!ey22σ12dy=n=01n12n+12nσ22nn!0y2n+1ey22σ12dy==n=01n12n+12nσ22nn!2n+1σ12n+2n!2.

Note that the last equality is obtained with the help of the formula:

0x2n+1epx2dx=n!2pn+1,p>0.

In our case, p is 12σ12. We continue the equalities chain which is broken at (*):

=n=01n12n+1σ12n+2σ22n=n=01nσ122n+1σ1σ22n==2πn=01nσ1/σ22n+12n+1=2πarctanσ1σ2.

For the k-moment calculation we take

qkt=λk+1tkeλtk!,λ=π22

and get the results shown in Table 9.

MomentSimulation
13,297101
21,556101
39,842102
48,066102
58,194102
69,968102
71,414101
82,293101
94,181101
108,474101
184,251103
191,637104
206,634104

### Table 9.

The results of numerical calculations for the moments by the third way.

## 6. Integral from the regression function

Now, we consider the issue of calculation of an integral

Dfxdx,E38

where the function fx has no an analytical expression. Suppose there exists a random variable ξxw such that its expectation is equals to Eξxw=fx for some fixed x. The random variable ξxw may be realized neither as result of the physical measurements or some calculations (e.g., using the modeling statistical method). In this case the optimal density is given by [6].

px=fxdx+λ,E39

where dx is the variance of the random variable ξxw. Note that one should use the optimal density from [1] if complexity in calculations (experimental measurements) is much different from each other for any x. We determine the parameter λ from the condition Dfx/dx+λdx=1. Really in practice, we find a priori or a posteriori approaches to both fx and dx. By f¯x and d¯x denote these approaches. Then the approach to the optimal px will look like

p¯x=f¯xd¯x+λ¯andDf¯xd¯x+λ¯dx=1.E40

The parameter λ is often turned out to be find enough complicity [6]. If the domain D is the interval 0H for small H then it is suppose to use the quasioptimal density p¯x.

Example 3. We now consider the following issue: Suppose fx=x, dx=1/x,D=0H. The optimal density is given by

px=x1/x+λ=xxλx+1cx3/2.E41

We take the quasioptimal density in the form p¯x=5H5/2x3/2/2. In this case, for px1 the estimation variance of the random value η=fxw/px:

varη=0Hdxpxdx+01f2pxdxI2E42

is equals to . Taking px=2x/H2 we get varη=2/H. But if the function fx was precisely known for the same density px then varη=0. If we choose the quasioptimal density p¯x=5H5/2x3/2/2 then the estimation variance of η is equals to 17H4/12+4/15H. For H0 the variance behaves approximately as 4/15H. It is much the better than 2/H. For H=1 the estimation variance with the density px=2x/H2 is equals to 2, and the estimation variance with the quasioptimal density p¯x is equals to 101/60.

Suppose we practically realize calculation of the integral Eq. (38) with dx=1/x; then one should discard the interval 0δ and to calculate δHfxdx because of the values ξxw can be the intolerably large. Also one should replace fx by f̂x:

f̂x=0,0xδ,x,δ<xH.E43

The shift is 0δxdx=δ2/2 and choosing δ1/N4 we get the total error δ2/2+3σ/N of oder O1/N.

In applications the estimation variance for the integral functionals (e.g., field flow calculation neither across the arc or the surface) from the solutions of the boundary-value problems for both the linear [1] or nonlinear [2] elliptic equations is of interest. For the above variance is dxB/x2,fxa0+a1x+a2x2+, where x is the distance to the domain boundary. Suppose fxa1x+a2x2+; then the optimal density is given by

px=a1x+a2x2+B/x2+λ.E44

The quasioptimal density has the form p¯x=3x2/H3 for small H in 0H. In applications, this case is of our main interest. Taking δ=1/N4 like in the Example 3 we get the asymptotics of decrease for the total error as O1/N.

Suppose dxB/x2,fxa0+a1x+a2x2+, and a00 then there is no density kind of p¯x=w+1xw,x0H with the finite variance. The density px=fx/dx+λ will be give the estimation with the infinity variance. Instead of calculation of the integral 0Hfxdx we will be calculate the integral δHfxdx. For this integral we already can choice the quasioptimal density with the finite variance of the estimation: p¯x=2x/H2δ2. For δOlnN/N the total error will have the asymptotics OlnN/N.

Example 4. Suppose dx=1/x2,fx=1,H=1; then the asymptotics of the variance with the quasioptimal density has kind of 25lnδ.

In conditions of Example 4, choice of the optimal density in the form

px=x1δ21x2,xδ1E45

yields the following result: the estimation variance will have asymptotics 2lnδ for δ0.

Remark. If we know that a value of fx in the interval 0δ close to the number f0, then in Eq. (43) to use

f̂x=f0,0xδ,x,δ<xH,E46

more efficiently and also to take 0Hfxdxf0δ+δHfxdx.

## 7. Conclusion

In the paper we describe the sufficient conditions of the stable calculations for the multidimensional integrals by the Monte Carlo method. We get the results of numerous numerical computations using the mt19937 pseudorandom number generator. The article results can be also useful in the practical solution of the boundary value problem, for both the elliptic and parabolic equations. The earlier suggested approach to the optimal choice of the density [1, 6] often needs to solve a complicated secondary task. In the paper we suggest the approach to choice of the quasioptimal densities that is of considerable interest in applied problems solution.

## References

1. 1. Kronberg AA. On the numerical finding of some functionals on solving of Laplace’s equation. Izv. Vuz. Matematika. 1987;31:30-36
2. 2. Kronberg AA. Algorithms for solving elliptic boundary value problems. UUSR Computational Mathematics and Mathematical Physics. 1982;22:122-131. DOI: 10.1016/0041-5553(82)90103-3
3. 3. Matsumoto M, Nishimura T. Mersenne twister: A 623-dimensionally equidistributed uniform preudorandom number. ACM Transactions on Modeling and Computer Simulation (TOMACS). 1998;8:G1-G25. DOI: 10.1145/272991.272995
4. 4. Dynkin EB. Markov Processes. Vol. 2. Berlin: Springer-Verlag; 1965. p. 274
5. 5. Haji-Sheiki A, Sparrow EM. The solution of heat conduction problems by probability methods. Transmission ASME Series C: Journal of Heat Transfer. 1967;89:121-131. DOI: 10.1115/1.3614330
6. 6. Kronberg AA. On the numerical solution of Laplace’s equation in unbounded domain. Izv. Vuz. Matematika. 1984;28:45-48

Written By

Alexander A. Kronberg and Tatiana K. Kronberg

Reviewed: 28 June 2022 Published: 04 October 2022