The Stationary Problem of Two-Phase Filtration by the Monte-Carlo
Method: Solutions and Applications
M. TASTANOV, A. UTEMISSOVA, F. MAIYER, R. YSMAGUL
Kostanay Regional University, Baytursynov Street 47, Kostanay 110000,
KAZAKHSTAN
Abstract: - This article is devoted to solving the problems of applying Monte-Carlo algorithms to filtration
problems. The "random walk by spheres" and "random walk by boundary" algorithms of Monte-Carlo methods
are used to solve the stationary problem of filtration of two immiscible inhomogeneous incompressible liquids
in a porous medium. Estimates constructed using the "random walk by spheres" and "random walk by
boundary" algorithms of Monte-Carlo methods will be mostly 𝜀 biased. Unbiased estimates are in most cases
unrealizable on a computer, since with a probability of 1 they do not go to the boundary of the region, and
therefore are of little use. In practice, they are usually limited to only the first two points of evaluation.
Key-Words: Monte-Carlo method, continuity equation, Dirichlet problem, Markov chains, algorithms of random
walk by spheres, random walk by boundary, estimation of the solution and derivatives of the solution.
Received: October 17, 2021. Revised: September 14, 2022. Accepted: October 16, 2022. Published: November 4, 2022.
1 Introduction
Numerical implementation of mathematical models
challenges the construction of effective algorithms
for the solution of engineering problems related to
filtration of liquids. To date, such challenges are
approached by state-of-the-art computational
methods based on the Monte Carlo algorithm. The
algorithms of the Monte-Carlo methods well
implement, firstly, multidimensional problems, and
secondly, with the help of the algorithms of the
Monte-Carlo methods, it is possible to find a
solution at a single point in a complex area, which is
a very relevant problem in underground hydraulic
mechanics (for example, determining the point of
greatest pressure in the area). In this paper, solutions
and derivatives of solutions to the above two-phase
filtration problems are evaluated by Monte- Carlo
methods. Estimates constructed using the "random walk
by spheres" and "random walk by boundary" algorithms
of Monte-Carlo methods will be mostly 𝜀 −biased.
Unbiased estimates in most cases are unrealizable on a
computer, since with a probability of 1 they do not go to
the boundary of the region, and therefore are of little
use. In practice, they are usually limited to only the first
two points of evaluation [1-5].
We consider the mathematical model of filtration
process of two immiscible liquids (water and oil)
through the porous medium.
2 Formulation of the Problem in
Saturations and Pressure (s, p)
From the equation of continuity follows that
(1) ,0 21 VVVVdiv
Here
V
velocity vector of mixture filtration
and
V
- total flow intensity of mixture.Let
introduce the new function of “reduced”
pressure
(2) ,
1
02
1hgd
k
k
s
p
pp c
where
,,
0201 ghgkkk
we present vector
V
through
,and Sp
and it is independent of s:
)3(,),,()( 0
KkKpsVfpKV
where
Here
)(
00 skKK ii
symmetric tensor of
phase permeability,
)(
0xK
filtration tensor for
homogeneous liquid. In the same way, using
(2), we have
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
770
Volume 21, 2022
.
).()( 02
1
02
1111
d
k
k
s
p
s
k
k
s
p
pKgppKV c
sx
c
By
0201
0201 kk
kk
s
p
ac
and
d
k
k
s
p
Kf c
sx02
1
0
we receive
(4) ).,(
10101 psVfpKsaKV
Using (3) we find
)(
1
11 fVKKpK
and
considering
and 00011KkKKkK
we receive
).(
1
01
1
1sbKkKK
Then (4) will be as follows
(5) . -
001 fbfF,F VbsaKV
In the equation of continuity (1), substituting for
the first phase of expression (4), we come to the
system relative to {s, p}:
(7) .0),()(
(6) ),,()( 1010
psVdivfpKdiv
psVdivfpKsaKdiv
t
s
m
For the initial boundary value problem (6)-(7),
we consider the filtration flow in the fixed finite
area with piecewise-smooth border ∂Ω [6-7].
Let
nTSTQ ii
],,0[],,0[
external
normal to ∂Ω. We determine boundary data in
relation to functions s, p. The impermeability
conditions on ∂Ω0 for both phases are as follows
(8) ].,0[),(,0 001 TStxnVnV
The boundary conditions are respectively
rearranged to the form
(9) ,],0[),(),,(),,( 2200 TStxtxsstxpp
(10) ],,0[),(),,()( 110 TStxtxRnVnfpK
(11) ),(),,()( 11010 StxtxRbnVnfpKsaK
At R(x,t)=0 equality (10) and (11) are
equivalent (8), and ∂Ω0 add in ∂Ω1 and assume that
∂Ω1 consists of several components, on parts as
R=0.
Therefore,
.
21
The initial
condition is set only for saturation s(x,t):
(12) ,),(
~
)0,( 0 xxsxs
Components ∂Ω1or ∂Ω2 on ∂Ω0 can absent, i.e.
.or 21
In case of
1
the
law of conservation of mass on the area Ω need for
following necessary condition:
(13) ].,0[,0),(),( TtdxtxRdxtxp
Coefficients of the equation and boundary
conditions (8) – (12) have the form:
,
)(
)(
),()(
,
)(
)()(),(
),(
1
1
01
0201
0201
KK
sk
sk
b(s)bskskk(s)
sk
sksk
s
sxp
sxaa c
)),()()((
),()(
020100
00
skskxKkKK
skxKK ii
,2,1),()(),(
],)([
,),()(),(),(
),,(),(),(
,))(,(),(),(
),(),(),(),(
000
120
1
02010
21
1222
0
1
1
ixKsksxK
gpK
kkksxfsbsxfsxFF
sxKsxKsxKK
gsxKsxpsxK
sxfsxKsxKsxff
i
cx
cx
(14) .
)(
)(
),(
)( 02
1
000
d
k
k
s
sxp
Ksff s
s
s
Further assume that total velocity of filtration
),()(or
21 psVfpKVVVV
is independent of
s
, in other words. if
coefficients
),( and )()(
0sxfsksKK
are
independent of
,s
and system of equation (6), (7) is
separated and can be consistently determined
).,( and txsV i
From
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
771
Volume 21, 2022
.g)-(
)(
)(),(
),( 1222
02
1
0
KpKd
k
k
s
sxp
Kksxf cs
c
s
x
follows that this assumption is right:
1) as k=const for miscible liquid. For
immiscible liquids the significant deviation k from
a constant is observed only near limit values s=0,1
given saturation;
2) as
),(spp cc
in other words
;00or 0
s
p
p
x
pc
xcx
i
c
3) if gravity is not considered, i.e. g=0 or liquids
have identical density
.
21
The assumptions 2)
and 3) provide the condition fulfillment
.0
s
f
3 Formulation of the Stationary
Problem
Stationary boundary value problem for (s(x),p(x))
will be as follows:
(20) ,on V
(19) ,on
(18) ,in 0)(
(17) ,on V
(16) ,on )(ss
(15) ,in 0)(
1
20
11
20
10
Rn
pp
VdivfpKdiv
Rbn
x
VdivFbVsaKdiv
If coefficients
),( and )()(
0sxfskxKK
are independent of s, it follows from (18) that
(21)
,0)())()((
))()()(()(
01
01
xhxpxKdivC
xfxpxKCdivfpKdiv
where
).()(,)(
1xfdivxhconstskC
Further we suppose that we have homogeneous
and isotropic earth. Then from (21) we receive
orxhxзCconstCxK )()( and )( 320
),()( 1xhxp
where
.
)(
)(and
3
1213 C
xh
xhCCC
(22))()()()()( 301 xfxpCxfxpxKCV
if the medium is homogeneous and isotropic.
Therefore, we obtain the problem:
(25) ,on V
(24) ,on ,)(
(23) ,in ),()(
1
20
1
Rn
pxp
xhxp
(23) - (25) is Dirichlet problem for Poisson
equation [8-11]. Further we consider the case,
when
.
1
Then from the conditions (13) we
receive
.0)()(
dxxRdxxp
As
21
then
,
2
we
further suppose that
.
2
Then from (23)
(25) we obtain Dirichlet problem for Poisson
equation:
(27) .on ),0()(
(26) ,in ),()(
0
1
pxp
xhxp
Solution, the first and second derivatives from
solutions of this problem, in other words
)( and )(),(xpxpxp
are estimated by the
Monte-Carlo method.
4 Evaluation of the Solution and
Derivatives of the Solution by Methods
Monte Carlo
We formulate algorithms of Monte-Carlo
methods for evaluating the solution and
derivatives of the solution of the problem (26)
(27).
Using a special Fredholm integral equation
of the 2nd kind with a degenerate kernel, the form
of which is determined by the "ball" Green
function of problem (26), (27) we construct an
algorithm for solving the Dirichlet problem for the
Poisson equation. It is known that this Green
function will take the form
.)(,
11
4
1
),( drdrr
drr
drG
Then, for a ball of radius d0 centered at a
point x0, the solution of the problem (26), (27) at a
point x0can be represented as
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
772
Volume 21, 2022
00
0
0
(28) ,)(),()(
),(
)( 10
0
)(
0
dxr
dr
xS
drrhxrGdyyp
n
xrG
xp
where
).(00 xdd
The first integral in (28) is the
integral over the surface of the sphere.
),(0
xS
the
second one is all over the sphere
.
00 dxr
Relation (28) can be considered as a
conjugate (according to the terminology accepted
in the theory of Monte-Carlo methods) Fredholm
integral equation of the 2nd kind with a
generalized kernel representing a uniform
probability distribution on the sphere
);(0
xS
after
the introduction of this kernel, the first integral in
(28) becomes three-dimensional. Let's return to
equation (28). It is not difficult to understand that
the standard Monte-Carlo algorithms apply to
such integral equations if the kernel features are
included in the transition density of the modeled
Markov chain. In this case, from this point
0
x
we
should move to the surface of the sphere.
);(0
xS
we call such a chain "random walk by
spheres". The relation (28) must be supplemented
with the following equality:
(29) , );()( 00 xxpxp
which means that the kernel of the integral
equation vanishes if the first argument
.x
Thus, after reaching the border, the chain should
be cut off, adding the value to the estimate
)(
0xp
with the appropriate weight.
These considerations lead to an unbiased
probabilistic evaluation of the solution at the
point.
,
0
x
which is unrealizable, since with
probability 1, "random walk by spheres" do not
reach the boundary in a finite number of steps.
This is due to the fact that the norm of the
integral operator in the space we are considering
L1 equal to 1.
Next, we will construct a "biased", but
realizable estimate of the solution of problem
(28) and estimate the amount of bias.
Suppose that the solution of the Dirichlet
problem is known at every point of the set
surrounding area of the border
.
Then for the function p(r) we can write the
following integral equation:
(30) ),()(),()( rrdrprrkrp
where
r
rr
rrk r
if,0
if),(
),(
(31)
if ),(
, if,)(
11
4
1
)(
0
rrp
rrdrh
drr
rdrr
Here
)(),(rrdd r
generalized density
corresponds to a uniform probability distribution
on the sphere S(r). It is obvious that the
Neumann series for equation (30) converges,
since the norm of the integral operator. K in a
natural space for this task L1 will be less than 1.
Really,

).(-1 )(
)()(
),(),(
rdr
rdrdrr
rdrdrrkrrk
r
rr
Here
max
2
max
2
where,
4
)( d
d
- the exact
upper bound of the radii of spheres lying entirely
in Ω. From here we obtain
.1)(1
)(
2
1
L
K
Consequently, this
relation ensures the convergence of the Neumann
series and, thereby, the possibility of using
Monte-Carlo methods for equation (30). Equation
(30) has the form of a conjugate integral
equation, so to estimate p(x0) we can apply the
ratio
(32) ,)()(,)(
1
00
N
ii
xxMxp
where
)( i
x
determined from (31). Here {xi} -
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
773
Volume 21, 2022
Markov chain, which we determine as follows:
)()( 0
xrrq
the density of the initial
distribution;
)(),( rrrq r
transition density
from
)(;in rgrr
the probability of a chain
termination, determined by the expression:
. r ,1
, r ,0
)(
rg
Such a chain is called "random walk by
spheres".
We assumed that the solution p(r) known in
.
This allows us to get ε - biased evaluation
of the solution. Indeed, instead of exact values
p(r) in
we can use approximate values, for
example, taking them from the nearest border
points, i.e., assuming:
, ),()( *
0
rrprp
.)(, **
rdrrr
Then
,))()(( *
crrMpp N
where
c is some constant that is finite due to the
limitation in the area
derivatives of the
solution.
It turns out that the integral expressing
, by )(
rr
it can be estimated by the
Monte-Carlo method by one random "node"
distributed with density
.0,
)/1(6
2
sin
2
1
)( 2
,, dx
d
dxx
xg
Appropriate random evaluation
)(r
has the
form
),(
6
),,( 1
2
1
rh
d
r
with
).(),,(
1rrM
We construct an algorithm of Monte-Carlo
methods for estimating the solution of problem
(26), (27) at a given point
0
x
:
1) from the point
0
x
modeling the chain
n
x
of Markov before the first hit in the
is the
neighborhood of the boundary
,
we
determine the point
.
*
N
x
Here
*
N
x
is the
point of the boundary closest to the last state
N
x
point of the boundary, N is the random number of
the last state of the Markov chain;
2) respectively, the density
dx
d
dxx
xg
0,
)/1(6
2
sin
2
1
)( 2
,,
in each random sphere S(x), entirely lying in the
area Ω, the value of the function is calculated
).(
6
),,( 1
2
rh
d
rФ
Here
r
radius
vector of a point
rrx
,
random distance
between the centers of the spheres,
isotropic
random unit vector,
,)( drdrr
the
density
)(
,, xg
is given in the polar coordinate
system.
3) We calculate a random variable
by the
formula
1
1
*
0000 ).(),,(),,(
N
nNnnn xpxФxФ
The desired estimate of the solution
)( 0
xp
at
a given point
0
x
is obtained by averaging over all
trajectories of the quantity
,
that is
(33) ,)(
1
)(
1
0
L
ii
L
xp
where L is the number of trajectories outgoing
from the point
.
0
x
The theorem. Variance of a random variable
uniformly bounded by ε, therefore,
. constD
As proof, due to the limited initial pressure
and initial saturation, i.e. due to the limited
function
)( and )( 00 xsxp
it is enough to assume
that
.0)(
0xp
We obtain (it is sufficient to
consider the case of the Dirichlet problem for the
Poisson equation, i.e. in the Helmholtz equation
we use с=0).
if ),(
, if,)(
11
4
1
)(
0
rrp
rrdrh
d
rr
rdrr
Because of
,in 1 00 cc
n
Q
it is enough
to consider the case c=0, in which
.1
n
Q
In this
case, the algorithm under study is a direct
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
774
Volume 21, 2022
modeling for the integral equation (30) [9].
The corresponding variance is expressed by
the following Ermakov - Zolotukhin formula
[10]:
(34) ,]2[, *
ffD
where
f
density of the centers of spheres, and
*
f
solution of the problem for a given
value
.
It can be shown that
.)( 1
d
c
df
(see, for
example, §3.4, [11]). At the same time
,)(, 2
32
*dcdcf
by virtue of the
determination φ(x), i.e.(31).
From here we obtain the statement of the theorem
by integrating the variance (34) along a fairly
narrow "border layer".
Now we estimate the derivative with respect
to the
),3,2,1( ixi
solution of problem (26) -
(27) at that point
0
x
. For this purpose, using
,
,0 x
x
we denote a point
3
0
2
0
1
0,, xxxx
by
(temporarily, the upper index will correspond to
the number of the independent variable, and the
lower one - to a fixed point). Then the solution to
problem (26), (27) at the point
x
x,0
for a sphere of
radius
)( 0
xd
has the form
(35)
.)(),(
)(),()(
0
0
10
)(
,0
drrhdrG
dyypxwxp
dxr
x
xS
x
The spherical Green's function
),( 0
drGx
and
its normal derivative
),( xw
are expressed by
the formulas:
,
2)(
))((
4
1
),(
,
2)(
1
4
1
),(
2
3
1
0
212
0
212
00
1
0
2
0
221
0
,0
0
axdxd
xdd
xw
xxddrx
d
xr
drG
p
x
x
where ω the unit vector of the direction from
),( and , to 10
xaaxx
the cosine of the angle
between ω and the axis
p
xx,
is the magnitude
of the projection of the vector
0
xr
onto the
axis
.)(, 00 dxdxx
We differentiate relation
(35) with respect to x, taking into account the
expressions for the definition
),( and ),( 0xwdrGx
, by setting x=0, we obtain:
)36(
)()(
4
1
3
)(
)(
4
1
)(
0
00
*
01
1
0
1
3
0
3
0
3
0
3
00
j
dxr
N
j
jj
N
jj
i
dxr
i
i
xprdrh
xr
xrd
dd
l
Mrdrh
xrd
xrdx
x
xp
Here
),(
ii xll
is the cosine of the angles
between
and coordinate axes
Mixi,3,2,1,
is the expectation operator. The first integrals in
expressions (36) can be estimated by the Monte-
Carlo method from one random "node" with
density
(37) ,
)(
)( 00
3
0
3
0
3
0
3
01
1dxr
xrd
xrdx
rfx
To evaluate all three derivatives by
i
x
it is
advisable to use density
(38) ,)(
~
)( 2
0
2
0
3
0
3
0
2
1
3
1
2
xrd
xrd
rfrf i
x
ix
where
., 0
2
1
32100 xrxxxdxr
In this case, the quality of the algorithm is
estimated by the sum of the variance of the
estimates of the corresponding integrals.
Integrals standing under the sign of
mathematical expectation in (36) can also be
evaluated by one random "node" distributed with
density
(39) ,)( 2
2
3
3
d.xr
xrd
xrd
rf
In work [8] algorithms for modeling random
variables with given distribution laws, namely
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
775
Volume 21, 2022
with densities, are given:
1)
.0,
)/1(6
2
sin
2
1
)( 2
,, dx
d
dxx
xg
2)
,)( 2
3
3
3
d.xr
xrd
xrd
rf
Random variables with such densities are
encountered when evaluating the solution p(x)
and derivatives of the solution
).3,2,1(,
)(
i
x
xp
i
Dirichlet problems for the Poisson equation.
First, let's consider the algorithm for
constructing a random variable with density
.0,
)/1(6
)( 3dx
d
dxx
xg
Replacement
d
x
y
leads to
(40) 10),1(6)( yyyyf
This expression corresponds to the distribution
of the second ordinal statistic of three sample
values of a random variable evenly distributed in
the interval (0, 1). And the guiding cosines
),3,2,1(, ili
i.e. the cosines of the angles
between the unit vector
),(
and coordinate
axes
321 ,, xxx
are modeled using the algorithm:
;21;21;21)1 321
that ,1if,)2 22 qq
;
1
;
1
;)3
2
3
2
21 q
l
q
ll
otherwise, paragraph 2 is done. Here
3,2,1for (0,1) i
i
sample values of a
random variable α, evenly distributed in (0,1).
Now let's consider the algorithm for modeling a
random variable η, distributed with density
.,)( 2
3
3
3
dxr
xrd
xrd
rf
Moving on to the polar coordinate system
with the center
x
and, making a replacement
,
d
xr
y
we find the density
.10,
3
)1(4
)(
3
y
y
yfn
Hence, the modeling
formula is obtained
(41)
421
Note that formula (36) makes it possible to
evaluate and
.(x)p
5 Assessment of
).(xV
Using the estimates
)3,2,1(,
i
x
(x)p
(x)p
i
from
formula (36), we find
.(x)p
Further, formula
(22) for a homogeneous isotropic medium gives
us
(42) ).()( 3 xf(x)pCxV
To determine the saturation s(x), we consider
equations (15) - (17).
Let us denote by the vector
W
(43) .FVbW
We remind that
,)( 1
1
KKsbb
.),( 0
1
0201 KkkksxFF
).,()(),()( 012 sxfsbsxfgpcx
from (14) and
,
01 KkK
we assumed that
constk
and
.
20 constCK
Then
.
4constCb
We also assumed that
),( sxf
does not depend on s, i.e.
).(),( xfsxf
If the
condition is satisfied 2) with respect to the
filtration tensor for a homogeneous liquid
,)(det
)(
1
if is, that ),(00 constxK
xm
xK
).3,2,1( ,0then
i
x
p
i
c
From here and from (4) it follows
).(
0
02
1
0xfd
k
k
s
p
Kf c
sx
Therefore, under some of the above
assumptions, we obtain the estimate
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
776
Volume 21, 2022
(44) .)( FVbxWW
Now, substituting (44) into (15) and taking into
account
,
20 constCK
we obtain
(46) on ),(
(45) in ,0)(
0
2
xss
WsaCdif
6 Conclusion
In this paper, the issues of applying Monte-Carlo
algorithms to filtration problems are
investigated. We were able to apply the
algorithms of "random walk by spheres" and
"random walk by boundary" Monte-Carlo
methods to solve the stationary problem of
filtration of two immiscible inhomogeneous
incompressible liquids in a porous medium. We
have constructed a 𝜀 biased estimate of the
solution and derivatives of the solution (the first
derivative of saturation, the first and second
derivatives of pressure) of the stationary problem
of two-phase filtration of incompressible
immiscible liquids in a porous medium. Using
the same Monte-Carlo algorithms, it is possible
to solve the filtration problem taking into
account temperature (i.e. the energy equation is
added to the filtration equations). The scientific
novelty of the research consists in the fact that,
theoretically and practically, the method of
statistical testing (Monte-Carlo) for solving
boundary value problems of stationary and non-
stationary filtration in areas of arbitrary
configuration, which has positive possibilities of
implementation on modern computers, is being
developed.
References:
[1] Temam R, Navier-Stokes equations. Theory
and numerical analysis. (Mir, Мoscow, 1981),
pp. 408.
[2] Ladyzhenskaya О.А.,Mathematical problems
in the dynamics of a viscous
incompressible fluid. (Nauka, Мoscow, 1970),
pp. 288.
[3] Belonosov S.M. and Chernous K.A., Boundary
value problems for the Navier- Stokes
equations. (Nauka, Мoscow, 1985), pp. 311.
[4] Ermakov S.M., Mikhailov G.A. Statistical
modeling. 3rd ed., supplemented. (M.:
Fizmatlit, 2000), pp. 298.
[5] Mikhailov G.A. Optimization of Monte Carlo
weight methods (Nauka, Мoscow, Gl. ed.
phys.-mat. lit., 2007), pp. 280.
[6] Ermakov S.M. Monte Carlo method and
related issues (Nauka, Мoscow, 2015), pp.
492.
[7] Antontsev S.N., Kazhihov А.V., Monahov
V.N., Boundary value problems in the
mechanics of inhomogeneous fluids. (Nauka,
Novosibirsk, 1983), pp. 319.
[8] Mishustin B.А., On the solution of the
Dirichlet problem for the Laplace equation by
the method of statistical tests. ZhVMiMF,
1967, Vol. 7, 5, pp. 1178-1187.
[9] Smagulov Sh.С. and Shakenov К.К.,On the
solution of the stationary problem of two-phase
filtration by the Monte Carlo method” in
Materials of the International conference
Actual problems of mathematics and
mathematical modeling of ecological systems
dedicated to the 60th anniversary of
academician Sultangazin U.M. (“Gylym-
Almaty”, 1996.), pp. 55.
[10] Shakenov К.К., Sultanova М.S., Tastanov
М.G. “Numerical solution the one Musket
Leverett’s model by MonteCarlo methods and
probability difference method” in Materials of
International conference “Informational
technology and mathematical modeling in
science, technology and education” dedicated
to 75th anniversary of academician A.
Djaynakov (Bishkek, Kyrgyzstan, 2016), pp.
220-226.
[11] Tastanov M.G., Utemissova A.A., Maier F.F.,
Ysmagul R.S. Solution of problems for the
relaxation model proceeding under linear
Darsy’s law using Monte- Carlo and
probabilistic-difference method. Global and
Stochastic Analysis, 2020, Vol. 7 No. 1,pp. 87-
98.
[12] Ermakov S.M., Monte Carlo method and
related issues.Second edition,expanded. (Nauka,
Moskow, 1975), pp. 346-363.
[13] Ermakov S.M. and Zolotukhin V.G.,
Application of the Monte Carlo method for
calculating protection against nuclear
radiation. Questions of physics of reactor
protection. (Atomizdat, Moscow, 1963),
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
777
Volume 21, 2022
pp.171-182.
[14] Yelepov B.S., Kronberg A.A., Mikhailov G.A.,
Sabelfeld K.K., Solving boundary value
problems by the Monte Carlo method. (Nauka,
Novosibirsk, 1980), pp.176
[15] N. Shimcin. Monte-Carlo Methods Lecture
Notes, Spring, 2015.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.e
n_US
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.88
M. Tastanov, A. Utemissova, F. Maiyer, R. Ysmagul
E-ISSN: 2224-2880
778
Volume 21, 2022