Monte Carlo method for constructing an unbelised assessment of diffusion problems Текст научной статьи по специальности «Математика»

Section 2. Mathematics

https://doi.org/10.29013/ESR-20-1.2-7-12

Tozhiev Tokhirjon Halimovich, PhD, Associate Professor E-mail: r_rasulov51@mail.ru Ibragimov Shavkat Mamirovich, Senior Lecture,

Kayumov Akhror Muminovich,

teacher

Abdullaev Shokhboz Solizhonovich,

teacher

Solizhonov Zhasur Valerievich,

Bachelor Ferghana State University

MONTE CARLO METHOD FOR CONSTRUCTING AN UNBELISED ASSESSMENT OF DIFFUSION PROBLEMS

Abstract. In this article is used the Monte Carlo method to solve Cauchy problems in a diffusion equation. A Markov chain is constructed with absorption on the trajectory of this chain, an unbiased estimate is constructed for solving problems. Examples of estimates of solution functionals are given and the behavior of their variances is discussed. The resulting integral equation is solved using branched Markov chains with absorption.

Keywords: diffusion, Monte Carlo method, approximate solutions, unbiased estimates, algorithm, algorithm efficiency.

Introduction. In the research of dynamic systems in computer modeling, Monte Carlo theories and methods are often used to solve linear boundary value problems of mathematical physics. This has been sufficiently researched in [1-9]. Tasks usually come down to solving equivalent integral equations.

Algorithms of the Monte Carlo method for problems with initial conditions have been little researched even in the case of ordinary differen-

tial equations. The only research with examples of Monte Carlo solutions for systems of ordinary differential equations was considered in [5]. In these cases, Monte Carlo methods may be more efficient than classical computational methods. If we keep in mind that the application of the Monte Carlo method, which has been shown to be effective in multidimensional problems, the development of numerical integration methods in a system with many noise is very relevant.

Statement of the problem. Consider the Cau-chy problem in the classical formulation in (n + 1) -dimensional space Rn+1 in a layer Q = Rn x [O,T]:

, w N du(x,t) * ti d2u(x,t) (Lu)(x,t) = —*—'--} a11-}—± +

V 1 >\ > ¿—t a.Ja.j

(1) - (2), we construct an algorithm for its numerical implementation.

To simplify the solutions of the problem, we introduce (in block notation) the following matrices

dt

j=1

dx dxj

+Z Z*

j=1 m=1

mx!

du (x ,t )

dxm

= f (x ,t ),

(1)

n x n as a = a(5) =

' L 0

u(x,0) = f(x), x e Rn, (2)

where aij u ¡3™ (i = 1,2,..., k, j = 1, 2,..., k, m = = 1, 2,..., l) are constant coefficients, n = k +1. Let a = (aij ) is a matrix with dimension k x k, f = (fm) is a matrix with dimension l x k .

Further suppose that: 1. a = (aij ) - is a symmetric matrix and is positive definite. 2. ¡3 = (¡m) a matrix such that the Gram matrix (© = (fia(3T )-1) with dimension l x l is positive definite, i.e. she is non-degenerate.

The functions f (x,t) and <(x) will be considered continuous in Rn. Under the assumption of the existence and uniqueness of the solution of problems

C (5 ) =

V

k

-5 A

im1m2 A

V mlm3 J f

q = q (5 ) = a 1 (5 ),

= exp

1i J d (p) = 1

0 0

V

-5$ 0,

( „1/2 T p Ik

0

3/2 j

P 1

i y

Here, m1 = — [a'1 + 3fiTa>pA is k x k dimension-

45L J

3

al block, m2 =—-/3T® a block with a size k x l, m^

25

3

a block with a size l x k, m3 = —a a block with a size

3 s

l x l, matrices Ik ,o, /, Il also have similar dimensions. Here in after, Ir is the identity matrix with a size r x r, d(p) is the diagonal matrix, s = t-t .

The matrices a(s) and q(s) for s > 0 are positive definite. Then the fundamental solution of equation (1) with a singularity at a point (y,t ) has the form

_____j_ I T

Z(x,t;y,r) = n 2 (HI2 (t_t) 2 exp\_(y _Cx) d

t _T

ad

t _T

(y _ Cx )!

(3)

where ||a|| = det (a), y = k + 3l. Using the fundamental solution Z (x ,t; y ,t ) in space Rn+1, we define areas (sharoid) as follows. Let the number (parameter) be positive, i.e. r > 0. Further area

Br (x,t) = j (y,t): (y - Cx)Td

Br (x,t) = \ (y,t ) : Z(x,t;y,r) >

n

a 2 r Y,t >t

will be called a sharoid with radius r whose center is located at a point (x,t). Then Br,(x,t) it can be rewritten as

( 1 ^

ad

V t-TJ

1

t-T

(y - Cx )<Y ln-^- ,t >t

(4)

t-t

From (4) it is seen that t satisfies the following toward the center (x,t). Therefore, there is r > 0 conditions: t < t,t > t - r2. Each section of a spheroid with a horizontal plane t = con5t, t - r2 < con5t < t,

such that (x,t)eQ Br(x,t)<=Q.

is an n-dimensional ellipsoid centered at a point ,t . If r ^ 0, then Bp(x ,t )<= Br (x ,t ). When

Let r > 0 such that Br (x,t) <= Q . Then, to solve problem (l)-(2), the following relation holds:

u (x ,t ) = (Er ,u )(x ,t) + f (x ,t), (5)

where

Br ( x, t) and dBr ( x ,t ) monotonously contracted

n ^ [ n ^

(Er,u)(x,t) = (-]2 jIn[- 12 JJ u(y(A,0),t(A))ht (9)4babTH(d)d5dX

\n J 0 y^J 5, (0)

0

2

biased estimate of the solution u (x ,t) to problem (1) - (2).

For u (x ,t) = 1, applying formula (5) we get

jx'-'

■nil

\\

V vx/y

dxj | HT (0)4babTH(9)ds = 1

si (0)

It follows that the core of integral equation (5) can be considered as the distribution density. Consider the integral

f (x,t)= J J Z(x,t;y,z)-n 2 ||a||2 rY f (y,z)dydT.

Br (x ,t) L _

Here S1 (0) is the ( n -1) - dimensional unit sphere with the usual orthogonal coordinates:

e = (e2,e3,...,en), 0<ej <n, for 2< j<n-1,

0 <9n < 2n . H(0) e S1 (0) is unit n - dimensional vector

1

y(1,0) = e^^x + (y lnfiYf d(r2!2)b-1H(0)

V V1)) ,(6)

t(1) = t - r212

1

It follows from (6) that if r (x ,t) < t2, then when

(x ,t )eQ we have Br (x ,t )<=Q. j = .

We proceed to the construction of a Markov ^ 2) 0

chain {(xj,tj)} on which we will construct an un- 1 w \ • a. c xt

U ' !\j=0 where T(n) is the gamma function. Now we can

imagine (Eru )(x ,t )as follows

J = K'

( ( 1 ^ 2 ■n -

V VXJJ

dX •

Making a change of variables X = e Y, we get

f

w n

J e ~pp2 d p

r

' + -

V 2

Y

(Eru)(x,t) = JP'(p)dp J P2 (H)u

0 S'(0)

f r

y

V v

e r ,6

r

,T

dS

V JJ

where P1 (p) is the density of the gamma-distributed can model a random vector with a distribution den-random variable with the parameter 1 + n, ... sity (7) by the Neumann method.

P2 (H) =

HT 4babTH Y°n '

Below is shown the algorithm of modeling. Algorithm: 1. An isotropic vector

In what follows, we will model a random vector a = (a>1,a2,...,an) and y2 a random variable distrib-

with a distribution density

P2 (H ) = P2 ((, H 2,...., Hn ) =

tKHH, xi (o)(H )

(7)

where XSl (0)(H) is the indicator of the set S1 (0), an is the surface of the unit sphere, K j are the elements of n x n dimensional matrix K = 4babT . Since H1,H2,....,Hn are the coordinates of the unit vector, H2 + H2 +... + H\ = 1 we get that HTKH < ¡u1, where /u1 is the largest eigenvalue of the matrix K. Then we

uted uniformly in (0,1) are modeled. 2. E = i±v.. 3. If

f n }

-- > E, then a it is accepted, otherwise

paragraph (l) is repeated.

Let {''} is a sequence of independent gamma distributed random variables with a parameter {coj} is a sequence of independent vectors with a distribution density P2 (H ). Now we define the Markov chain in Q the following recurrence relations:

V 2 y

x0 = x, t0 = t, tj = tj-1 - exp

y J

f n \

xj = xj +t, -'exp

j I '

Y )

j Yb-

t>; / v*^in

CO■

j / j~im i m='

Xk+p = Xk+p tj-1 exp

f m A

+ ij-exp

Y J c=1 p j

3A

Y J

% &+p ,

w

,v i

(8)

where i = 1,2,...,k, p = 1,2,...,l, j = 1,2,...,

i

r (x1 -1,tj-1 ) = (-1 )2 and relation (8) is obtained from (6). Now (5) we can write in the form

u(x1 -1,tj-1) = E(xj-1,tj-1)(x1 ,tj) + f (x1 -1,tj-1). (9)

Theorem 1. a) the sequence {r¡¡ forms a martingale with respect to the sequence of a - algebras (3¡}j=0; b) if uf,0 (x,t)<+oo and u^,0 (x,t)<+oo , then n it is quadratically integrable.

Proof: First, we prove that {n }j=0 it forms a mar-

We define a sequence ofrandom variables {n }°_0 tingale. From the definition 3 it is clear that n is

, , . ,, measurable 3,, in then using the conditional expec-

by the following equality, l

n =£ h (x1 ,t1 )f (y1 ) + u (xl ,tl), (10)

1=1 """(x,t)L'IM' ^lJ_ ■ 0 where (y1 ,t1 )isa random point of a sharoid Br (x,t)

^ ' It follows that {n} _ it forms a martingale rela-

at fixed (J ),having a distribution densityin it as tively ^. Let us p°rove that E(x,t)n2 . It is

tation property and formulas (9, 11) we get

E{xt)[n] = lf h(xj,tj)f(yjTj) + u(xlt) = n

1

Z(xj,tj;y,r)-n 2 ||a||2 r~

h (xj ,tj ) Z(x,t;y,r)-n 2 ||a||2 r~

n 1

dydx.

enough to show that

l-1

I = E(x,t)(Z^(xj,tj)f (yj,tj))2 <

j=0

Dividing I into two terms r2 (x,t) < t, we obtain from the final condition that h (x ,t ) <

where h (x ,t )= JJ

Br (x ,t)

Assuming that u(x,t) = t and applying formula this implies the proof of the theorem. (9) for 1 = 1 from (8) we get

f

Y

KY + 2 y

t,

h(x,t) = r2 (x,t)( I(y + 2))) 2J. (11)

Let is a sequence of a algebras generated

by random variables |1,|2, ..., ^ and a sequence of vectors o0,©1,...,©' and random points (yV),(y 1,T1),...,(yl_1,t!_1) .We denote the solu- f(x,t) =

tions to problem (1) - (2) corresponding to those where given f,y as uf„(x,t).

Now we show one of the methods for estimating a random node of a quantity

_ r -n 1

f (x,t)= JJ z(x,t;y,r)-n 2IHI2 r— f (y,r)dydr.

Br (x ,t )L _

Lemma 1. For the function f (x ,t), the relation

' y ^

vY + 2 y

r 2Bf (y (Ç,Z,œ),r(Ç,Z)),

y ($>C>®)=e

-r exp I--- Z

Y+2

px + $

Y + 2

a

r exp

V V

r 2$ ^ 2/

T(Ç>Z) =t - exP

Y + 2

Z21 Y

Y + 2 J

Proof: We introduce the region Br = {(y,t) : yTa(1 fr)ad(1 /t )y < y /2lnr2/T, t > 0

(12)

Here | is gamma when a distributed random The resulting mirror image of the spheroids variable with a parameter r((2), Z is beta distrib- Br(0,0) relative to the plane t = 0. These areas will

is a random unit vector.

\/r>

uted random variable with a parameter ( (Y 2 ), © also be calle2d sharoides (radius r ). Then we have

f (x ,t ) = JÍ7jJ[r rT7'2 exp (- yTd (1/r)ad (1/r)y )-

V=1

B

x f (e~*x + y ,t -r^dydr We make a change ofvariables and some integral transformation and obtain the proof of Lemma 1.

Consider the question of the computational re-alizability of estimate (10). Take E small enough and consider (SQ)£ = (iT *[(U]) . Let

Ns = min {l: (xl, tl) e (5Q)£} bethe moment of the

first hit of the process (xl,tl) in (5Q)S, i.e. Ns is the moment the process stops (Markov moment). Lemma 2. There is an inequality:

n

/ \1+-( Y = 2 1 2 t

{ Y J s

Proof. Having taken u(x,t) = t and applying formulas (10) and (11) we get

N„-1

t = U1>0(X,t) > E(xt) X h (xj ,tj ) =

í \ Y

Y + 2

j=i

Ns -1

Ext) S '2 (xj ,tj )

j=i

It follows r(x,t) from the definition that r2(xj,tt) = {tj}>£ . Hence we get that

t >

í \ Y

1+n/2

Y + 2

£E(x ,t )N.

Consequently

Ex,t)N <

Y + 2

,1+n/2

—. The lemma is proved.

s

\ Y J

Theorem 2. Let the con ditions of Theorem 1 be

Proof. It follows from Theorem 1 that n it is qua-dratically integrable and hence it is uniformly integrable and Ns < +QO. Further, the moment the process is stopped is a Markov moment. Therefore, according to the Doob theorem "On the transformation of free choice" [2] and formula (9) EnN = En, = u(x,t) i.e. nN¡¡ is an unbiased estimate for u(x,t). From nN and n the definition of random variables it can be seen that DnN < Dn. From nNe with the standard method, a mixed, but practically realizable estimate . Let <p1(x,0) = $(x) x e Rn and (x,t ) be the point 5Q closest to (x,t).

In

the

N,-1

r¡Nc = ^ h(xj ,tj )f (yj Tj ) + u(xN« ,tN«

assessment

replace

j=0

w(x,tNs ) and u(x,t* s ) with get

N* -1

= X h(x' )f (yj T ) + ç>i(x,t* )

* j=0

Theorem 3. Let u(x,t) the Lipschitz condition and A(s) the modulus of continuity be satisfied. Then the random variable is a biased estimate for u(x,t). is limited parameter s function.

Proof: Since E(x t )rqN = u(x ,t ) then

1 u(x,t) - E(x,tl=| E(x,t)1n. - E(xtvN. H

1=1 E(x t)u(x,tN ) - E(x>t)^i(x,t) I< £(x>t) I u(x,tN* ) -- u(x,t*N" )I= A(s)

The theorem is proved.

satisfied. Then nNs is an unbiased estimate for u(x,t). Its dispersion is finite.

References:

1. Kolmogorov A. N. Is ber die analitichen Metho der in der Wahrscheinlichkeizsrechnung. Mathemat. Annalen. 1931.- 104 s.- S. 415-458.

2. Doob J. L. Classical Potential Theory and its. Probabilistic Counterpart. Springer - Varlag. 1984.- 846 p.

3. Ermakov S. M., Rasulov A. S., Bakoev M. T., Vaselovskaya A. Z. Selected Monte Carlo Method Algorithms / - Tashkent: University, 1992.- 132 p.

4. Mikhailov G. A. and Voytishek A. V Numerical Statistical Modeling. Monte Carlo Methods (Akademiya,-M., 2006). [in Russian].

5. Akhtar M. N., Durad M. H. and Ahmed A. "Solving initial value ordinary differential equations by Monte Carlo method," Proc. IAM 4, 2015.- P. 149-174.

6. Tozhiev T. Kh. Ibragimov Sh. M. Stochastic approximation methods for solving diffusion problems. "Fundamental and Applied Scientific Research: Actual Issues, Achievements, and Innovations" collection of articles of the XVI International Scientific and Practical Conference.- Penza: ICSN "Science and Enlightenment". 2018.- P. 13-15.

7. Ermakova S. M., Tovstika T. M. Monte Carlo Method for Solving ODE Systems. Vestnik St. Petersburg University, Mathematics, 2019.- Vol. 52.- No. 3.- P. 272-280.

8. Bakoev M. The solution of the mixed problem for the Kolmogorov equation // Problems of computational and applied mathematics.- No. 2. 2019.- P. 60-70. (01.00.00., No. 9).

9. Tozhiev T. H., Ibrahimov Sh. M. Numerical solutions of the Cauchy problem for the generalized equation of nonisotropic diffusion. Scientific Bulletin of NamSU 2019.- No. 10.- P. 33-42.