# A New Algorithm for System of Integral Equations

Authors: Abdujabar Rasulov; Adem Kilicman; Zainidin Eshkuvatov; Gulnora Raimova
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We develop a new algorithm to solve the system of integral equations. In this new method no need to use matrix weights. Beacause of it, we reduce computational complexity considerable. Using the new algorithm it is also possible to solve an initial boundary value problem for system of parabolic equations. To verify the efficiency, the results of computational experiments are given.


---

## Body

## 1. Introduction

The theory and application of integral equations are an important subject within pure and applied mathematics and they appear in various types in many fields of science and engineering. The integral equations can also be represented as convolution integral equations; see Srivastava and Buschman [1]. In the applications, the number of computational problems can be reduced to the solution of a system of integral equations (system of IEs) of the second kind; see [2–4]. However, solving systems of integrodifferential equations are very important and such systems might be difficult analytically, so many researchers have attempted to propose different numerical methods which are accurate and efficient. For example, numerical expansion methods for solving a system of linear IDEs by interpolation and Clenshaw Curtis quadrature rules were presented in [5], where the integral system was transferred into a matrix equation by the interpolation points. Pour in [6] studied an extension of the Tau method to obtain the numerical solution of Fredholm integrodifferential equations systems ad applied Chebyshev basis to solve IDEs. Similarly, Arikoglu and Ozkol [7] obtained solutions of integral and integrodifferential equation systems by using differential transform method where the approch provides very good approximation to the exact solution.Recently, the solution of the system has been estimated by many different basic functions, such as orthonormal bases and wavelets; see, for example [8, 9], and the hybrid Legendre Block-Pulse functions, that is, a combination of the Block-Pulse functions on [
0,1
] and Legendre polynomials was proposed. In addition, the Bessel matrix method was introduced in [10] for solving a system of high order linear Fredholm differential equations with variable coefficients. In the literature there are several methods to solve the different type of integral equations; see [11–16]. One of the novel methods is known as the vector Monte Carlo algorithms to solve the system of IEs. Among the vector Monte Carlo algorithms the following are well known:(i)
an algorithm for solving the system of transfer equations with polarization;
(ii)
a vector algorithm for solving multigroup transfer equations;
(iii)
a Monte Carlo technique combined with the finite sum method and vector Monte Carlo method for solving metaharmonic equations.In the use of this method one can easily see that the variance of the vector estimate largely depends on the form of transitional density. Thus appropriate choice of the density leads to the reduction of the complexity calculations, which is defined as the product of the variance and the computational time. To determine the density is difficult as to solve the problem itself, although in some cases it is possible to obtain a minimal criterian of uniform optimality of the method. The transitional density that corresponds to minimum complexity of algoritm is said to be optimal for a given problem.In Mikhailov [17], vector Monte Carlo algorithms are used to solve system of IEs. The distinguished feature of that vector algorithm is that its “weight” appears in the form of a matrix weight. This matrix weight is multiplied by the kernel matrix of the system of IEs dividing by a transition density function in the Markov chain simulation, so that a number of computational problems can be reduced to the solution of a system of IEs of second kind. By introducing a suitable discrete-continuous measure of the integration, we can write the system of IEs in the form of a single integral equation, and this allows us to use standard algorithms of the Monte Carlo method. However, it is more expedient to make use of the matrix structure of the system and solve the problem by the Monte Carlo method with vector weights. The following vector Monte Carlo algorithms are well known: an algorithm for solving the system of transfer equations with polarization taken into account, a vector algorithm for solving multigroup transfer equations, a Monte Carlo technique combined with the finite sum method, and vector Monte Carlo method for solving metaharmonic equation.In this study, a new algorithm is proposed for the numerical solution of system of IEs but in this algorithm we do not use matrix weights. The proposed algorithm has usual advantages of ordinary Monte Carlo method. The new algorithm is considerably reduced to computational complexity. Using this new algorithm we have solved an initial boundary value problem for system of parabolic equations. The paper is organized as follows. In Section2, we present the description of the problem and proposed a new Monte Carlo algorithm for the solution of system of IEs. In Section 3, we discuss the application of the method to the solution of system of parabolic equations. In Section 4, we will construct biased and ɛ-biased estimators for the solution. In Section 5, the results of computational experiments are given, followed by the conclusion in Section 6.


## 2. Description of the Problem and a New Approach for the Solution of System of IEs

Let us consider second kind nonhomogeneous system of IEs of the form(1)
φ
i
(
x
)
=
∑
j
=
1
n
∫
X
k
i
,
j
(
x
,
y
)
φ
j
(
y
)
d
y
+
h
i
(
x
)
,
i
=
1
,
…
,
n
,
where x
∈
X
⊆
R
m
,
m
≥
1 or in vector form
(2)
Φ
=
K
Φ
+
H
,
here operator K
:
L
∞
→
L
∞ where L
∞ is the space of bounded function almost everywhere and
(3)
H
=
(
h
1
,
…
,
h
n
)
∈
L
∞
,
K
=
(
K
i
j
)
∈
L
∞
,
Φ
=
(
φ
1
,
φ
2
,
.
.
.
,
φ
n
)
∈
L
∞
,
where the norm of H is
(4)
∥
H
∥
L
∞
=
vrai
sup
⁡
1
≤
i
≤
n
,
x
∈
X
|
h
i
(
x
)
|
.
Suppose the spectral radius ρ
(
K
) satisfy the inequalities
(5)
ρ
(
K
)
=
lim
⁡
n
→
∞
∥
K
n
∥
1
/
n
<
1
,
ρ
(
K
)
≤
∥
K
n
∥
1
/
2
,
where K
n
φ
=
K
n
-
1
K
φ.Let Markov chain{
x
n
,
n
=
0,1
,
…
,
N
} with transition density p
(
x
,
y
) be
(6)
g
(
x
)
=
1
-
∫
X
p
(
x
,
y
)
d
y
≥
0
,
where g
(
x
) is the probability of absorption at the point x
N, where N is the random number of the last moment and in initial moment x
0
=
x.A standard vector algorithm of Monte Carlo forΦ
(
x
) is
(7)
Φ
(
x
)
=
(
M
ξ
x
,
ξ
x
)
=
H
(
x
)
+
∑
n
=
1
N
Q
n
H
(
x
n
)
,
Q
0
=
I
,
Q
n
+
1
=
Q
n
K
(
x
n
,
x
n
+
1
)
P
(
x
n
,
x
n
+
1
)
,
n
=
0,1
,
2
,
…
,
where I is a unit matrix, K
(
x
,
y
) is a kernel matrix {
k
i
j
(
x
,
y
)
}, and p
(
x
n
,
x
n
+
1
) is the transition density function at the points (
x
n
,
x
n
+
1
). The condition for unbiasedness is
(8)
p
(
K
1
)
<
1
or
K
1
=
∥
K
∥
<
1
.
We will assume also that the spectral radius of the operator K
1 obtained from K by the substitution k
i
,
j
→
|
k
i
,
j
| is less than one. Then, by using standard methods of Monte Carlo theory we can show that
(9)
Φ
(
x
)
=
E
ξ
x
,
ξ
x
=
∑
n
=
0
N
Q
n
H
(
x
n
)
,
Q
0
=
{
δ
i
,
j
}
i
,
j
=
(
1
,
…
,
n
)
,
Q
n
=
Q
n
-
1
K
(
x
n
-
1
,
x
n
)
p
(
x
n
-
1
,
x
n
)
,
where Q
n can be considered as matrix weight and
(10)
K
(
x
n
-
1
,
x
n
)
p
(
x
n
-
1
,
x
n
)
=
{
k
i
j
(
x
n
-
1
,
x
n
)
p
(
x
n
-
1
,
x
n
)
}
,
i
,
j
=
{
1,2
,
…
,
n
}
.
The Monte Carlo method is used to estimate linear functionals of the form
(11)
(
F
,
Φ
)
=
∫
X
F
′
(
x
)
Φ
(
x
)
d
x
,
where F
′
(
x
)
=
(
f
1
(
x
)
,
f
2
(
x
)
,
…
f
n
(
x
)
) with
(12)
∥
F
∥
L
1
=
∑
j
=
1
n
∫
X
|
f
j
(
x
)
|
d
x
<
∞
.
Let the point x
0 be distributed with initial probability density π
(
x
) such that
(13)
π
(
x
)
≠
0
if
F
′
(
x
)
Φ
(
x
)
≠
0
.
Then, obviously (see Mikhailov [17]),
(14)
(
F
,
Φ
)
=
E
[
F
′
(
x
0
)
π
(
x
0
)
ξ
x
0
]
=
E
[
∑
n
=
0
N
F
′
(
x
0
)
π
(
x
0
)
Q
n
H
(
x
n
)
]
=
E
[
∑
n
=
0
N
H
′
(
x
n
)
Q
n
′
F
′
(
x
0
)
π
(
x
0
)
]
.
The random vector with weight Q
n
(
1
)
=
Q
n
′
(
F
′
(
x
0
)
/
π
(
x
0
)
) is computed by the formula
(15)
Q
n
(
1
)
=
K
′
(
x
n
-
1
,
x
n
)
p
(
x
n
-
1
,
x
n
)
Q
n
-
1
(
1
)
.
Precisely such a vector algorithm, corresponding to the representation I
=
(
Φ
*
,
H
), has been formulated in the work of Mikhailov [17]. Below on the contrary to vector algorithms we will propose a new algorithm for the solution of system of integral equations. Our method does not use matrix weights.Suppose we have to find the solution of the inhomogeneous system of IEs (1) of the second kind at the pointx
∈
X. We will define two types of Markov chain {
i
K
} and {
x
K
} by the following way.(a) Definition of the First Homogeneous Markov Chain.Now we simulate the Markov chain i
0
,
i
1
⋯
∈
N with n
+
1 state. Initial state i
0 will simulate according to initial distribution π
=
(
π
1
,
…
,
π
n
,
0
) and the next i
1 with the transition matrix
(16)
A
=
A
(
x
)
=
∥
α
i
,
j
(
x
)
∥
i
,
j
=
1
n
+
1
,
∑
j
=
1
n
α
i
,
j
(
x
)
=
1
-
g
i
(
x
)
,
g
i
(
x
)
=
α
i
,
n
+
1
(
x
)
,
i
=
1
,
…
,
n
.
Here α
n
+
1
,
n
+
1
(
x
)
=
1 and
(17)
A
(
x
)
=
(
α
11
(
x
)
,
…
,
α
1
n
(
x
)
,
g
1
(
x
)
α
21
(
x
)
,
…
,
α
2
n
(
x
)
,
g
2
(
x
)
α
n
1
(
x
)
,
…
,
α
n
n
(
x
)
,
g
n
(
x
)
0
,
…
,
0
,
1
)
.
Let N be a random absorption moment with N
=
{
max
⁡
k
,
i
k
≠
n
+
1
}, a life time of chain.(b) A second homogeneous Markov chain
{
x
k
} with space phase X is defined by the following way.Firstly, we define the transition density matrix as(18)
P
(
x
1
⟶
x
)
=
∥
P
i
,
j
(
x
1
⟶
x
)
∥
i
,
j
=
1
n
+
1
,
P
(
x
1
⟶
x
)
=
(
P
11
(
x
1
⟶
x
)
,
…
,
P
1
n
(
x
1
⟶
x
)
,
0
P
21
(
x
1
⟶
x
)
,
…
,
P
2
n
(
x
1
⟶
x
)
,
0
P
n
1
(
x
1
⟶
x
)
,
…
,
P
n
n
(
x
1
⟶
x
)
,
0
0
,
…
,
0
,
1
)
.
Let an initial point x
0
=
x; using π
(
x
) we will simulate initial moment i
0, then according to the transition matrix A
(
x
0
) we are able to simulate again the next state of chain i
1. It means with the probability α
i
0
,
i
1
(
x
0
)
,
P
(
i
1
=
k
)
=
α
0
k
(
x
0
).The next phase coordinates of the chainx
1 simulated according to p
i
0
,
i
1
(
x
0
,
x
1
). The probability of absorption of the trajectory is g
i
0
(
x
0
). Let (
i
k
,
x
k
) be known then the next value of i
k
+
1 will be defined according to the matrix A
(
x
k
) and next random point x
k
+
1 simulated according to the probability density function P
i
k
,
i
k
+
1
(
x
k
1
→
x
) and so on.Letξ
x
0
=
ξ
N
(
i
0
,
i
1
,
…
i
N
;
x
0
,
x
1
,
…
x
N
) be some random variable which is defined by the set of trajectory Markov chains. The mathematical expectations of random variable will be
(19)
E
ξ
x
0
=
∑
k
=
0
∞
‍
∑
i
0
,
.
.
i
k
=
1
n
‍
∫
⋯
∫
‍
︸
k
π
i
0
α
i
0
,
i
1
(
x
0
)
p
i
0
,
i
1
(
x
0
x
1
)
·
α
i
1
,
i
2
(
x
1
)
p
i
1
,
i
2
(
x
1
,
x
2
)
⋯
α
i
k
-
1
,
i
k
(
x
k
-
1
)
×
p
i
k
-
1
-
1
,
i
k
(
x
k
-
1
,
x
k
)
g
i
k
(
x
k
)
ξ
k
(
i
0
,
…
i
k
,
x
0
,
…
x
k
)
d
x
1
⋯
d
x
k
.
Let us consider calculation of the functional (
Φ
,
F
), where F
T
=
(
f
1
(
x
)
,
.
.
.
,
f
n
(
x
)
) column vector. Let us compute the functional (
Φ
,
F
)
=
∑
i
=
1
n
φ
i
(
x
)
f
i
(
x
). For doing this task we introduce two well-known estimators according to the Monte Carlo theory. First of them is analog of absorption estimator
(20)
ξ
1
(
x
0
)
=
f
i
0
(
x
0
)
π
i
0
(
x
0
)
k
i
0
i
1
(
x
0
,
x
1
)
α
i
0
,
i
1
(
x
0
)
p
i
0
,
i
1
(
x
0
,
x
1
)
·
·
·
k
i
n
-
1
,
i
n
(
x
n
-
1
,
x
n
)
α
i
n
-
1
,
i
n
(
x
n
-
1
)
p
i
n
-
1
,
i
n
(
x
n
-
1
,
x
n
)
·
h
i
n
(
x
n
)
g
i
n
(
x
n
)
and the second one is analog of collision estimator
(21)
ξ
2
(
x
0
)
=
∑
j
=
1
N
‍
f
i
0
(
x
0
)
π
i
0
(
x
0
)
K
i
0
i
1
(
x
0
,
x
1
)
α
i
0
,
i
1
(
x
0
)
p
i
0
,
i
1
(
x
0
,
x
1
)
·
·
·
K
i
j
-
1
,
i
j
(
x
j
-
1
,
x
j
)
α
i
j
-
1
,
i
j
(
x
j
-
1
)
p
i
j
-
1
,
i
j
(
x
j
-
1
,
x
j
)
·
h
i
j
(
x
j
)
.Theorem 1.
Iff
i
0
(
x
i
0
)
≠
0 then π
i
0
(
x
i
0
)
≠
0 and if k
i
i
i
j
(
x
i
,
x
j
)
≠
0 then
(22)
α
i
,
j
(
x
i
)
p
i
i
i
j
(
x
i
,
x
j
)
≠
0
,
f
o
r
a
n
y
1
≤
i
,
j
≤
n
.
In this case E
ξ
1,2
=
(
Φ
,
F
).The proof of the theorem is similar to the theorem Ermakov [23], and therefore proof is omitted. Now we will apply the obtained results to the solution system of parabolic equations.


## 3. Application to System of Parabolic Equations

In this section we consider initial boundary problem for system of parabolic equations. LetD be bounded domain in R
m with enough smooth boundary ∂
D of Ω
=
D
×
[
0
,
T
] and Ω is the cylinder in R
m
+
1 with parallel spin axis t. The basement is the domain D on the surface t
=
0 and T is the fixed constant. The functions
(23)
y
0
i
(
x
)
∈
C
(
D
¯
)
,
y
i
(
x
,
t
)
∈
C
(
∂
D
×
[
0
,
T
]
)
,
f
i
(
x
,
t
)
∈
C
(
Ω
¯
)
,
where C
(
D
¯
) stands for a continuous function on the closed domain D.Now consider the following initial boundary value problem (BVP) for system of parabolic equations:(24)
∂
u
1
(
x
,
t
)
∂
t
-
a
1
Δ
u
1
(
x
,
t
)
+
c
11
u
1
(
x
,
t
)
-
c
12
u
2
(
x
,
t
)
-
⋯
-
c
1
n
u
n
(
x
,
t
)
=
f
1
(
x
,
t
)
∂
u
2
(
x
,
t
)
∂
t
-
a
2
Δ
u
2
(
x
,
t
)
+
c
22
u
2
(
x
,
t
)
-
c
21
u
1
(
x
,
t
)
-
⋯
-
c
2
n
u
n
(
x
,
t
)
=
f
2
(
x
,
t
)
⋮
∂
u
n
(
x
,
t
)
∂
t
-
a
n
Δ
u
n
(
x
,
t
)
+
c
n
n
u
n
(
x
,
t
)
-
c
n
1
u
1
(
x
,
t
)
-
⋯
-
c
(
n
-
1
)
n
u
n
-
1
(
x
,
t
)
=
f
n
(
x
,
t
)
,
where the coefficients a
i
>
0
,
c
i
j
>
0
,
(
i
=
1
,
…
,
n
,
j
=
1
,
…
,
n
), and (
x
,
t
)
∈
Ω with initial and boundary conditions
(25)
u
i
(
x
,
t
)
=
y
i
(
x
,
t
)
,
x
∈
∂
D
,
t
∈
[
0
,
T
]
,
i
=
1
,
n
¯
u
i
(
x
,
0
)
=
y
o
i
(
x
)
,
x
∈
D
,
i
=
1
,
n
.
¯
Further suppose f
i
(
x
,
t
)
,
y
o
i
(
x
)
,
y
i
(
x
,
t
), and coefficients  a
i
,
c
i
j
(
i
,
j
=
1
,
n
¯
) are given such that there exists unique solution Ladyzhenskaya et al. [18] and Lions [19] of the initial BVR (24)-(25) and
(26)
u
i
(
x
,
t
)
∈
C
(
D
-
×
[
0
,
T
]
)
∩
C
2,1
(
D
¯
×
[
0
,
T
]
)
(
i
=
1
,
n
¯
)
,
where C
2,1 is the set of continuous functions in the given region with continuous derivatives u
x
,
u
x
x
, and u
t.Now we construct unbiased estimator for the problem (24)-(25) in the arbitrary point (
x
,
t
)
∈
Ω on the trajectory some random process. For that we use mean value formula and construct some special system of integral equations for u
i
(
x
,
t
) in special constructed domains (spheroid or balloid with the center (
x
,
t
)).According to Section2 below we will propose a new nonstationary Markov chain on which trajectory will construct unbiased estimators for the obtained system of integral equations.In our algorithm we do not used matrix weight; it means the computational complexity of new algorithm is much better. The basis for the constructing of algorithms will be the formula of parabolic mean for the heat conductivity equations. As we know the fundamental solutionZ
(
x
,
t
,
y
,
τ
) for heat equation u
t
-
a
Δ
u
=
0 is given by
(27)
Z
(
x
,
t
;
y
,
τ
)
=
(
4
π
a
(
t
-
τ
)
)
-
m
/
2
exp
⁡
(
-
|
x
-
y
|
2
4
a
(
t
-
τ
)
)
.
Firstly, we define a special domain using a fundamental solution of the heat equation Q
r
(
x
,
t
) which depends on r
>
0 and points (
x
,
t
)
∈
R
m
+
1as
(28)
Q
r
(
x
,
t
)
=
{
(
y
,
τ
)
:
Z
(
x
,
t
;
y
,
τ
)
>
(
4
π
a
r
)
-
m
/
2
,
τ
<
t
}
.
The domain Q
r
(
x
,
t
), we call balloid and ∂
Q
r
(
x
,
t
), spheroid with the center in the point (
x
,
t
). From the definition balloid Q
r
(
x
,
t
), described by following inequality (Kupcov [20]):
(29)
Q
r
(
x
,
t
)
=
{
(
y
,
τ
)
:
|
x
-
y
|
2
<
2
m
a
(
t
-
τ
)
ln
⁡
r
t
-
r
,
τ
<
t
}
.
Each section with the sectional plain of balloid when τ
=
constant will be m-dimensional ball B
(
x
,
R
(
t
-
τ
)
) with the center x and with the radius
(30)
R
(
t
-
τ
)
=
2
m
a
(
t
-
τ
)
ln
⁡
r
t
-
r
.
Let (
x
,
t
)
∈
Ω and
(31)
r
=
r
(
x
,
t
)
=
min
⁡
{
R
1
2
(
x
)
e
2
a
m
,
t
}
,
where R
1
(
x
) is the minimum distance from point x until the boundary; that is,
(32)
R
1
(
x
)
=
inf
⁡
{
|
x
-
x
′
|
,
x
∈
∂
D
,
x
′
∈
D
¯
}
.
In this case Q
r
(
x
,
t
)
⊂
Ω
.
¯ By further using Greens function and fundamental solution we will transfer from the system of differential equations into the system of integral equations. In the book [21] special balance equation analogies were constructed as in [22], which connected the value of function u
(
x
,
t
) with its integral from the spheroid and balloid with the center in the point (
x
,
t
).Lemma 2 (Kurbanmuradov [22]).
Let the functionu
(
x
,
t
) satisfy the following equation:
(33)
∂
u
(
x
,
t
)
∂
t
-
a
Δ
u
(
x
,
t
)
=
f
(
x
,
t
)
,
(
x
,
t
)
∈
Ω
.
Then the following formula of mean  is true (mean value formula):
(34)
u
(
x
,
t
)
=
a
∬
∂
Q
r
(
x
,
t
)
(
1
-
t
-
τ
r
)
×
(
-
∂
Z
(
x
,
t
;
y
,
τ
)
∂
n
y
)
u
(
y
,
τ
)
d
s
d
τ
+
1
r
∬
Q
r
(
x
,
t
)
Z
r
(
x
,
t
;
y
,
τ
)
u
(
y
,
τ
)
d
y
d
τ
+
F
r
(
x
,
t
)
,
where
(35)
F
r
(
x
,
t
)
=
1
r
∬
Q
r
(
x
,
t
)
(
r
-
(
t
-
τ
)
)
Z
r
(
x
,
t
;
y
,
τ
)
f
(
y
,
τ
)
d
y
d
τ
,
Z
r
(
x
,
t
;
y
,
τ
)
=
Z
(
x
,
t
;
y
,
τ
)
-
(
4
π
a
r
)
-
m
/
2
,
here ds is the element of small area of sphere ∂
B
(
x
,
R
(
t
-
τ
)
). In further using these results we will get special integral representation.


### 3.1. Transforming a System and Obtaining Integral Representation

Let us define the family of domainsQ
r
(
i
)
(
x
,
t
), which depends on positive parameters r
>
0 and point (
x
,
t
)
∈
R
m
+
1,where
(36)
Q
r
(
i
)
(
x
,
t
)
=
{
(
y
,
τ
)
:
Z
r
i
(
x
,
t
;
y
,
τ
)
>
0
,
τ
<
t
}
,
where Z
r
i
(
x
,
t
;
y
,
τ
)defined analogous Z
r
(
x
,
t
;
y
,
τ
) changing a for a
i (see above Lemma 2). The domain Q
r
i
(
x
,
t
) we will call a balloid with radius r which a center in a point (
x
,
t
) and a boundary ∂
Q
r
(
i
)
(
x
,
t
)
=
{
(
y
,
τ
)
:
Z
r
i
(
x
,
t
;
y
,
τ
)
=
0
,
τ
≤
t
} is spheroid. Here
(37)
r
=
r
(
x
,
t
)
=
min
⁡
{
R
1
2
(
x
)
e
2
a
1
m
,
…
,
R
1
2
(
x
)
e
2
a
n
m
,
t
}
.
Let (
x
,
t
)
∈
Ω and D
:
R
1
(
x
)
=
inf
⁡
{
|
x
-
x
′
|
,
x
∈
∂
D
,
x
′
∈
D
¯
} where R
1
(
x
) is the distance from the point (
x
,
t
) to the boundary of domain. In this case Q
r
i
(
x
,
t
)
⊂
Ω
.
¯ Appling the expression (34) to each of the equations we will get the following system of integral equations (
i
=
1,2
,
…
,
n
):
(38)
u
i
(
x
,
t
)
=

a
i
∬
∂
Q
r
i
(
x
,
t
)
(
1

-

t
-
τ
r
)
(

-

∂
Z
(
i
)
(
x
,
t
;
y
,
τ
)
∂
n
y
)
u
i
(
y
,
τ
)
d
s
d
τ
+
1
r
∬
∂
Q
r
i
(
x
,
t
)
(
1
-
(
r
-
(
t
-
τ
)
c
i
i
)
)
×
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
u
i
(
y
,
τ
)
d
y
d
τ
+
1
r
∬
∂
Q
r
i
(
x
,
t
)
(
r
-
(
t
-
τ
)
)
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
×
∑
j
=
1
,
n
;
i
≠
j
‍
c
i
j
u
j
(
y
,
τ
)
d
y
d
τ
+
1
r
∬
Q
r
i
(
x
,
t
)
(
r
-
(
t
-
τ
)
)
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
f
i
(
y
,
τ
)
d
y
d
τ
,
i
=
1
,
…
,
n
,
where
(39)
Z
(
i
)
(
x
,
t
;
y
,
τ
)
=
(
4
π
a
i
(
t
-
τ
)
)
-
m
/
2
exp
⁡
(
-
|
x
-
y
|
2
4
a
i
(
t
-
τ
)
)
,
Z
r
i
(
x
,
t
;
y
,
τ
)
=
Z
(
i
)
(
x
,
t
;
y
,
τ
)
-
(
4
π
a
i
(
t
-
τ
)
)
-
m
/
2
.
The derived system (38) is similar to system IEs which was considered in Section 2. That is way we can use the method which was given in Section 2.


### 3.2. The Probabilistic Representation of the Solution

After some transformation we will get for separate terms of the system (38) as follows:
(40)
I
1
(
i
)
(
x
,
t
)

=

a
i
∬
∂
Q
r
i
(
x
,
t
)
(
1
-
t
-
τ
r
)
(

-

∂
Z
(
i
)
(
x
,
t
;
y
,
τ
)
∂
n
y
)
u
i
(
y
,
τ
)
d
s
d
τ
=
(
1
-
q
m
)
∫
0
∞
q
1
(
ρ
)
d
ρ
∫
S
1
(
0
)
q
2
(
ω
)
u
i
(
y
(
i
)
(
ρ
,
ω
)
,
τ
(
ρ
)
)
d
s
=
(
1
-
q
m
)
E
u
i
(
y
(
i
)
(
ξ
,
ω
)
,
τ
(
ξ
)
)
,
where
(41)
q
1
(
ρ
)
=
ρ
m
/
2
exp
⁡
⁡
(
-
ρ
)
(
1
-
exp
⁡
(
-
2
ρ
m
)
)
×
(
(
1
-
q
m
)
Γ
(
1
+
m
2
)
)
-
1
,
q
2
(
ω
)
=
1
σ
m
=
Γ
(
m
2
)
(
2
π
m
/
2
)
-
1
,
q
m
=
(
1
+
2
m
)
-
(
1
+
m
/
2
)
,
y
(
i
)
(
ξ
,
ω
)
=
x
+
4
r
ξ
a
i
exp
⁡
(
-
2
ξ
m
)
ω
,
τ
(
ξ
)
=
t
-
r
exp
⁡
(
-
2
ξ
m
)
,
ξ is a random variable with density functions q
1
(
ρ
)
,
ω random point on the surface S
1
(
0
), which has a density function q
2
(
w
)
,
S
1
(
0
) unit sphere, ds element of surface, σ
m square of the surface unit sphere, and Γ
(
·
) Gamma function.Let us consider the second terms of (38)
(42)
I
2
(
i
)
(
x
,
t
)
=
1
r
∬
Q
r
i
(
x
,
t
)
(
r
-
(
t

-

τ
)
)
Z
r
i
(
x
,
t
;
y
,
τ
)
f
i
(
y
,
τ
)
d
y
d
τ
=
r
q
m
∫
0
1
q
3
(
ν
)
d
ν
∫
0
∞
q
4
(
z
)
d
z
·
∫
S
1
(
0
)
(
1
-
ν
2
/
m
exp
⁡
(
-
2
z
m
+
2
)
)
×
f
i
(
y
1
(
i
)
(
z
,
ν
,
w
)
,
τ
1
(
z
,
ν
)
)
d
S
w
=
r
q
m
E
{
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
×
f
i
(
y
1
(
i
)
(
ξ
1
,
v
,
ω
)
,
τ
1
(
ξ
1
,
v
)
)
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
}
,
where q
3
(
ν
)
=
(
1
-
ν
)
ν
2
/
m
-
1
(
B
(
2
,
2
/
m
)
)
-
1. Then the density of Beta distribution with parameters (
2
,
2
/
m
)
(43)
q
4
(
z
)
=
exp
⁡
(
-
z
)
z
m
/
2
-
1
(
Γ
(
m
2
)
)
-
1
,
and the density of Gamma distribution with parameters m
/
2,ω—unit random vector,
(44)
y
1
i
(
ξ
1
,
v
,
ω
)

=

x
+
[
4
m
m

+

2
r
a
i
ξ
1
v
2
/
m
exp
⁡
(
-
2
ξ
1
m

+
2
)
]
1
/
2
ω
,
τ
1
(
ξ
1
,
v
)
=
t
-
r
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
,
where ξ
1 is the random variable with density function q
4
(
z
) and v is another random variable with the density function q
3
(
ν
).Let(45)
r
=
r
(
x
,
t
)
=
min
⁡
{
e
R
1
2
(
x
)
2
m
a
1
,
…
,
e
R
1
2
(
x
)
2
m
a
n
;
1
c
11
,
…
,
1
c
n
n
;
t
}
,
then Q
r
(
i
)
¯
(
x
,
t
)
∈
Ω and the function
(46)
p
1
(
i
)
(
x
,
t
;
y
,
τ
)
=
[
1
-
(
r
-
(
t
-
τ
)
)
c
i
]
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
r
q
m
(
1
-
r
q
1
m
c
i
i
)
×
I
Q
r
(
i
)
(
x
,
t
)
(
y
,
τ
)
is the transition density in Q
r
i
(
x
,
t
) with fixed point (
x
,
t
), where
(47)
q
1
m
=
1
-
1
2
(
m
+
2
m
+
4
)
1
+
m
/
2
.
Let (
y
2
(
i
)
,
τ
2
(
i
)
) be a random point of balloid Q
r
i
(
x
,
t
) which has the following density function p
1
(
i
)
(
x
,
t
;
y
,
τ
)
(
i
=
1
,
…
,
n
) in the fixed point (
x
,
t
).In this case(48)
1
r
∬
Q
r
(
i
)
(
x
,
t
)
(
1
-
(
r
-
(
t
-
τ
)
)
c
i
i
)
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
u
i
(
y
,
τ
)
d
y
d
τ
=
q
m
(
1
-
r
c
i
i
q
1
m
)
E
u
i
(
y
2
(
i
)
,
τ
2
(
i
)
)
,
i
=
1
,
.
.
.
,
n
.
The obtained results we will put to (38) and we will get the probabilistic representation of problem (24)-(25). It follows there from that we could to following proposition.Theorem 3.
For the solution of initial BVP (24)-(25) the following probabilistic representation is valid:
(49)
u
i
(
x
,
t
)
=
(
1
-
q
m
)
E
u
i
(
y
(
i
)
(
ξ
,
ω
)
,
τ
(
ξ
)
)
+
q
m
(
1
-
r
c
i
i
q
1
m
)
E
u
i
(
y
2
(
i
)
,
τ
2
(
i
)
)
+
q
m
r
E
{
∑
j
=
1
,
j
≠
i
n
c
i
j
u
j
(
y
1
(
i
)
(
ξ
1
,
v
,
ω
,
)
,
τ
1
(
ξ
1
,
v
)
)
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
×
∑
j
=
1
,
j
≠
i
n
c
i
j
u
j
(
y
1
(
i
)
(
ξ
1
,
v
,
ω
,
)
,
τ
1
(
ξ
1
,
v
)
)
}
+
q
m
r
E
(
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
aaaaaaaaaa
×
f
i
(
y
1
(
i
)
(
ξ
1
,
v
,
ω
)
,
τ
1
(
ξ
1
,
v
)
)
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
)
,
(
i
=
1
,
n
¯
)
,
where (
y
(
i
)
(
ξ
,
ω
)
,
τ
(
ξ
)
) is defined by (40) and (
y
(
i
)
(
ξ
,
v
,
ω
)
,
τ
(
ξ
,
v
)
) is determined by (41).The proof of Theorem3 is the consequence of the above mentioned reasoning.By further using the presentation (44) we will construct a random process in Ω and propose the Monte Carlo algorithm for the solution of system IEs.


### 3.3. Description Random Process and the Algorithm Simulation

Let(50)
r
=
r
(
x
,
t
)
=
min
⁡
{
e
R
1
2
(
x
)
2
m
a
1
,
…
,
e
R
1
2
(
x
)
2
m
a
n
;
1
c
11
,
…
,
1
c
n
n
;
t
}
,
The functions
(51)
p
0
(
i
)
(
x
,
t
;
y
,
τ
)
=
1
1
-
q
m
(
1
-
t
-
τ
r
)
(
-
∂
Z
(
i
)
(
x
,
t
;
y
,
τ
)
∂
n
y
)
×
I
∂
Q
r
(
i
)
(
x
,
t
)
(
y
,
τ
)
,
p
1
(
i
)
(
x
,
t
;
y
,
τ
)
=
(
1
-
(
r
-
(
t
-
τ
)
)
c
i
i
)
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
r
q
m
(
1
-
r
q
1
m
c
i
i
)
×
I
∂
Q
r
(
i
)
(
x
,
t
)
(
y
,
τ
)
,
p
2
(
i
)
(
x
,
t
;
y
,
τ
)
=
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
r
q
m
I
∂
Q
r
(
i
)
(
x
,
t
)
(
y
,
τ
)
,
are the transition density functions in Q
r
(
i
)
¯
(
x
,
t
) at a fixed point (
x
,
t
)
(
i
=
1
,
n
¯
). We will define in Ω a random process as was proposed in Section 2.Let us define a transition matrix as(52)
A
(
x
,
t
)
=
[
α
11
α
12
…
α
1
(
n
+
1
)
α
21
α
22
…
α
2
(
n
+
1
)
…
…
…
…
α
n
1
α
n
2
…
α
n
(
n
+
1
)
0
0
…
1
]
,
where α
i
i
=
1
-
q
m
c
i
i
q
1
m
r
(
x
,
t
), and let
(53)
β
i
=
q
m
c
i
i
q
1
m
r
(
x
,
t
)
(
n
-
1
)
n
,
α
i
j
=
β
i
c
i
j
M
i
,
(
i
,
j
=
1
,
n
¯
;
i
≠
j
)
,
M
i
=
∑
j
=
1
,
n
;
j
≠
i
c
i
j
,
(
i
=
1
,
n
¯
)
,
α
i
(
n
+
1
)
=
q
m
q
1
m
r
(
x
,
t
)
c
i
i
n
,
(
i
=
1
,
n
¯
)
.
Now we will define the density function of transition matrix P
(
x
,
t
;
y
,
τ
):
(54)
P
(
x
,
t
;
y
,
τ
)
=
[
p
11
(
x
,
t
;
y
,
τ
)
p
12
(
x
,
t
;
y
,
τ
)
.
.
.
p
1
(
n
+
1
)
(
x
,
t
;
y
,
τ
)
p
21
(
x
,
t
;
y
,
τ
)
p
22
(
x
,
t
;
y
,
τ
)
.
.
.
p
2
(
n
+
1
)
(
x
,
t
;
y
,
τ
)
.
.
.
.
.
.
.
.
.
.
.
.
p
n
1
(
x
,
t
;
y
,
τ
)
p
n
2
(
x
,
t
;
y
,
τ
)
.
.
.
p
n
(
n
+
1
)
(
x
,
t
;
y
,
τ
)
0
0
.
.
.
1
]
,
where
(55)
p
i
j
(
x
,
t
;
y
,
τ
)
=
p
2
(
i
)
(
x
,
t
;
y
,
τ
)
,
(
i
=
1
,
n
¯
;
j
=
1
,
n
+
1
¯
;
i
≠
j
)
,
p
i
i
(
x
,
t
;
y
,
τ
)
=
(
1
-
q
m
)
p
0
(
i
)
(
x
,
t
;
y
,
τ
)
+
q
m
(
1
-
r
q
1
m
c
i
i
)
p
1
(
i
)
(
x
,
t
;
y
,
τ
)
1
-
r
(
x
,
t
)
q
1
m
c
i
i
q
m
.
Then we will fix the initial point (
x
0
,
t
0
)
=
(
x
,
t
) and the number of equations i
0
∈
{
1
,
…
,
n
}. Let an initial moment at the point (
x
0
,
t
0
)
=
(
x
,
t
); we will have one particle. For one step a particlei
k
→
i
k
+
1 moves from its position according to the transition matrix A
(
x
k
,
t
k
) and moves with probability α
i
k
,
i
k
+
1
(
x
k
,
t
k
) from the point (
x
k
,
t
k
) to the point (
x
k
+
1
,
t
k
+
1
). The next (
x
k
+
1
,
t
k
+
1
) point will be simulated using the density function p
i
k
i
k
+
1
(
x
k
,
t
k
;
y
,
τ
).The probability of breaking of trajectory in the point(
x
n
,
t
n
) is
(56)
g
(
x
n
,
t
n
)
=
{
1
,
(
x
n
,
t
n
)
∈
Ω
;
α
i
n
,
n
+
1
(
x
n
-
1
,
t
n
-
1
)
,
(
x
n
,
t
n
)
∈
Ω
.
The next coordinate of the particle will be defined in the following way.(1) If the density function of the point(
x
n
+
1
,
t
n
+
1
) equals p
0
(
i
)
(
x
n
,
t
n
;
y
,
τ
) in the fixed point (
x
n
,
t
n
) then
(57)
x
n
+
1
=
x
n
+
2
(
r
(
x
n
,
t
n
)
ξ
n
a
i
)
1
/
2
exp
⁡
(
-
ξ
n
m
)
ω
n
,
t
n
+
1
=
t
n
-
r
(
x
n
,
t
n
)
exp
⁡
(
-
2
ξ
n
m
)
,
where {
ξ
n
}
n
=
0
∞
,
{
ω
n
}
n
=
0
∞ the sequence of independent random variables with the density function q
1
(
ρ
) and independent isotropic vectors. The value r
(
x
n
,
t
n
) will be defined as (50).(2) If the density function of the point(
x
n
+
1
,
t
n
+
1
) is equal to p
1
(
i
)
(
x
n
,
t
n
;
y
,
τ
) at a fixed (
x
n
,
t
n
) then
(58)
x
n
+
1
=
x
n
+
2
(
(
-
2
ξ
n
′
m
+
2
)
m
m
+
2
r
(
x
n
,
t
n
)
×
ξ
n
′
(
v
n
′
)
2
/
m
a
i
exp
⁡
⁡
(
-
2
ξ
n
′
m
+
2
)
)
1
/
2
ω
n
,
t
n
+
1
=
t
n
-
r
(
x
n
,
t
n
)
(
v
n
′
)
2
/
m
exp
⁡
(
-
2
ξ
n
′
m
+
2
)
,
where {
ξ
n
′
}
n
=
0
∞
,
{
v
n
′
}
n
=
0
∞ is a sequence of independent random variables, which will be obtained from the algorithm below (Algorithm 4) (Neumann acceptance rejection method).Algorithm 4.
(a) We firstly simulateξ, Gamma distributed random variable with the parameters (
m
/
2
), secondly simulate γ, uniformly distributed random variable on (0, 1), and thirdly simulate ν, Beta distributed random variable with the parameters (2
,
2
/
m).
(b) Ifγ
>
1
-
c
i
i
r
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
/
(
m
+
2
)
)
) then we will go to (a) and so on; otherwise v
′
=
v
,
ξ
′
=
ξ.
(3) If the density function at the point(
x
n
+
1
,
t
n
+
1
) equals p
2
(
i
)
(
x
n
,
t
n
;
y
,
τ
) under fixed point (
x
n
,
t
n
), then
(59)
x
n
+
1
=
x
n
+
2
(
(
-
2
ξ
n
m
+
2
)
m
m
+
2
r
(
x
n
,
t
n
)
ξ
n
(
v
n
)
2
/
m
×
a
i
exp
⁡
(
-
2
ξ
n
m
+
2
)
)
1
/
2
ω
n
,
t
n
+
1
=
t
n
-
r
(
x
n
,
t
n
)
(
v
n
)
2
/
m
exp
⁡
(
-
2
ξ
n
m
+
2
)
,
where {
ξ
n
}
n
=
0
∞
,
{
v
n
}
n
=
0
∞
,
{
ω
n
}
n
=
0
∞ is sequence of independent Gamma distributed random variables with parameters (
2
/
m
), Beta distributed random variables with parameters (
2
,
2
/
m
), and independent isotropic vectors, respectively.If at the momentn was held break, then we will put (
x
n
+
k
,
t
n
+
k
)
=
(
x
n
,
t
n
)
,
k
=
0,1
,
2
,
.
.. obviously the sequence of coordinate of the particle forms Markov chains. The random process which was described above was considered in Ermakov et al. [23] for the solution of initial BVP for the heat equation and adapted in Kurbanmuradov [22] for the heat equation with variable coefficients.Now we prove the auxiliary Lemma5.Lemma 5.
With the probability one Markov chain{
x
n
,
t
n
}
n
=
0
∞ converges when n
→
∞ to the random point of boundary (
x
∞
,
t
∞
)
∈
∂
Ω, or it is absorbed inside of the domain.Proof.
Since{
t
n
} is decreasing sequence and t
n
≥
0, it has a limit t
∞
=
lim
⁡
n
→
∞
t
n. Let R
n
-
σ be algebra, which was generated by random variables
(60)
{
ω
k
}
k
=
0
n
-
1
,
{
v
k
}
k
=
0
n
-
1
,
{
ξ
k
}
k
=
0
n
-
1
,
{
v
k
′
}
k
=
0
n
-
1
,
{
ξ
k
′
}
k
=
0
n
-
1
.
From the definition R
n and (57)–(59) it follows that x
n is measurable relatively R
n. It is obvious that the coordinates of vector process formed limited martingale relatively {
R
m
}
m
=
1
∞:
(61)
E
(
x
n
+
1
R
n
)
=
E
{
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
x
n
R
n
(
1
-
q
m
)
×
[
x
n
+
2
(
r
(
x
n
,
t
n
)
ξ
n
a
i
n
)
1
/
2
exp
⁡
(
-
ξ
n
m
)
ω
n
]
+
q
m
(
1
-
q
1
m
r
(
x
n
,
t
n
)
c
i
n
i
n
)
×
[
exp
⁡
(
-
ξ
n
′
m
+
2
)
ω
n
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
ξ
n
′
v
n
′
2
/
m
)
1
/
2
×
exp
⁡
(
-
ξ
n
′
m
+
2
)
ω
n
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
-
1
n
×
[
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
ς
n
v
n
2
/
m
)
1
/
2
×
exp
⁡
(
-
ς
m
+
2
)
ω
n
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
ς
n
v
n
2
/
m
)
1
/
2
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
x
n
R
n
}
=
(
1
-
q
m
)
[
(
ξ
n
1
/
2
exp
⁡
(
-
ξ
n
m
)
ω
n
)
x
n
+
2
(
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
×
E
(
ξ
n
1
/
2
exp
⁡
(
-
ξ
n
m
)
ω
n
)
]
+
q
m
(
1
-
q
1
m
c
i
n
i
n
r
(
x
n
,
t
n
)
)
×
[
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
×
E
(
(
ξ
n
′
v
n
′
2
/
m
)
exp
⁡
(
-
ξ
n
m
+
2
)
ω
n
)
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
-
1
n
×
[
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
×
E
(
(
ξ
n
′
v
n
′
2
/
m
)
exp
⁡
(
-
ξ
n
m
+
2
)
ω
n
)
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
-
1
n
×
[
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
×
E
(
(
ς
n
v
n
2
/
m
)
1
/
2
exp
⁡
(
-
ς
n
m
+
2
)
ω
n
)
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
x
n
=
(
1
-
q
m
)
x
n
+
q
m
(
1
-
q
1
m
c
i
n
i
n
r
(
x
n
,
t
n
)
)
x
n
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
x
n
=
x
n
,
where {
x
n
} are limited martingale; it converges with the probability one Shiryaev [24].
Let(
x
∞
,
t
∞
)
=
lim
⁡
n
→
∞
(
x
n
,
t
n
) be the limit vector. We show that (
x
∞
,
t
∞
)
∈
∂
Ω. If t
∞
=
0 then the process is broken inside the domain. Let t
∞
>
0. As far as the process converges, according to the formulas (50)–(58) we have
(62)
E
(
x
0
,
t
0
)
|
x
n
+
1
-
x
n
|
⟶
0
,
E
(
x
0
,
t
0
)
|
x
n
+
1
-
x
n
|
=
E
(
x
0
,
t
0
)
{
r
(
x
n
,
t
n
)
h
(
r
(
x
n
,
t
n
)
)
}
,
where h
(
r
) is strictly positive. Applying Lebesque Theorem (about the limited convergence) we get
(63)
E
(
x
0
,
t
0
)
(
r
(
x
∞
,
t
∞
)
)
1
/
2
=
0
.
It means r
(
x
∞
,
t
∞
)
=
0. Then from the definition of r
(
x
,
t
) and using the formulae (49) we obtain
(64)
R
1
(
x
∞
)
=
0
,
(
x
∞
,
t
∞
)
∈
∂
Ω
.
Lemma is proven.


## 3.1. Transforming a System and Obtaining Integral Representation

Let us define the family of domainsQ
r
(
i
)
(
x
,
t
), which depends on positive parameters r
>
0 and point (
x
,
t
)
∈
R
m
+
1,where
(36)
Q
r
(
i
)
(
x
,
t
)
=
{
(
y
,
τ
)
:
Z
r
i
(
x
,
t
;
y
,
τ
)
>
0
,
τ
<
t
}
,
where Z
r
i
(
x
,
t
;
y
,
τ
)defined analogous Z
r
(
x
,
t
;
y
,
τ
) changing a for a
i (see above Lemma 2). The domain Q
r
i
(
x
,
t
) we will call a balloid with radius r which a center in a point (
x
,
t
) and a boundary ∂
Q
r
(
i
)
(
x
,
t
)
=
{
(
y
,
τ
)
:
Z
r
i
(
x
,
t
;
y
,
τ
)
=
0
,
τ
≤
t
} is spheroid. Here
(37)
r
=
r
(
x
,
t
)
=
min
⁡
{
R
1
2
(
x
)
e
2
a
1
m
,
…
,
R
1
2
(
x
)
e
2
a
n
m
,
t
}
.
Let (
x
,
t
)
∈
Ω and D
:
R
1
(
x
)
=
inf
⁡
{
|
x
-
x
′
|
,
x
∈
∂
D
,
x
′
∈
D
¯
} where R
1
(
x
) is the distance from the point (
x
,
t
) to the boundary of domain. In this case Q
r
i
(
x
,
t
)
⊂
Ω
.
¯ Appling the expression (34) to each of the equations we will get the following system of integral equations (
i
=
1,2
,
…
,
n
):
(38)
u
i
(
x
,
t
)
=

a
i
∬
∂
Q
r
i
(
x
,
t
)
(
1

-

t
-
τ
r
)
(

-

∂
Z
(
i
)
(
x
,
t
;
y
,
τ
)
∂
n
y
)
u
i
(
y
,
τ
)
d
s
d
τ
+
1
r
∬
∂
Q
r
i
(
x
,
t
)
(
1
-
(
r
-
(
t
-
τ
)
c
i
i
)
)
×
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
u
i
(
y
,
τ
)
d
y
d
τ
+
1
r
∬
∂
Q
r
i
(
x
,
t
)
(
r
-
(
t
-
τ
)
)
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
×
∑
j
=
1
,
n
;
i
≠
j
‍
c
i
j
u
j
(
y
,
τ
)
d
y
d
τ
+
1
r
∬
Q
r
i
(
x
,
t
)
(
r
-
(
t
-
τ
)
)
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
f
i
(
y
,
τ
)
d
y
d
τ
,
i
=
1
,
…
,
n
,
where
(39)
Z
(
i
)
(
x
,
t
;
y
,
τ
)
=
(
4
π
a
i
(
t
-
τ
)
)
-
m
/
2
exp
⁡
(
-
|
x
-
y
|
2
4
a
i
(
t
-
τ
)
)
,
Z
r
i
(
x
,
t
;
y
,
τ
)
=
Z
(
i
)
(
x
,
t
;
y
,
τ
)
-
(
4
π
a
i
(
t
-
τ
)
)
-
m
/
2
.
The derived system (38) is similar to system IEs which was considered in Section 2. That is way we can use the method which was given in Section 2.


## 3.2. The Probabilistic Representation of the Solution

After some transformation we will get for separate terms of the system (38) as follows:
(40)
I
1
(
i
)
(
x
,
t
)

=

a
i
∬
∂
Q
r
i
(
x
,
t
)
(
1
-
t
-
τ
r
)
(

-

∂
Z
(
i
)
(
x
,
t
;
y
,
τ
)
∂
n
y
)
u
i
(
y
,
τ
)
d
s
d
τ
=
(
1
-
q
m
)
∫
0
∞
q
1
(
ρ
)
d
ρ
∫
S
1
(
0
)
q
2
(
ω
)
u
i
(
y
(
i
)
(
ρ
,
ω
)
,
τ
(
ρ
)
)
d
s
=
(
1
-
q
m
)
E
u
i
(
y
(
i
)
(
ξ
,
ω
)
,
τ
(
ξ
)
)
,
where
(41)
q
1
(
ρ
)
=
ρ
m
/
2
exp
⁡
⁡
(
-
ρ
)
(
1
-
exp
⁡
(
-
2
ρ
m
)
)
×
(
(
1
-
q
m
)
Γ
(
1
+
m
2
)
)
-
1
,
q
2
(
ω
)
=
1
σ
m
=
Γ
(
m
2
)
(
2
π
m
/
2
)
-
1
,
q
m
=
(
1
+
2
m
)
-
(
1
+
m
/
2
)
,
y
(
i
)
(
ξ
,
ω
)
=
x
+
4
r
ξ
a
i
exp
⁡
(
-
2
ξ
m
)
ω
,
τ
(
ξ
)
=
t
-
r
exp
⁡
(
-
2
ξ
m
)
,
ξ is a random variable with density functions q
1
(
ρ
)
,
ω random point on the surface S
1
(
0
), which has a density function q
2
(
w
)
,
S
1
(
0
) unit sphere, ds element of surface, σ
m square of the surface unit sphere, and Γ
(
·
) Gamma function.Let us consider the second terms of (38)
(42)
I
2
(
i
)
(
x
,
t
)
=
1
r
∬
Q
r
i
(
x
,
t
)
(
r
-
(
t

-

τ
)
)
Z
r
i
(
x
,
t
;
y
,
τ
)
f
i
(
y
,
τ
)
d
y
d
τ
=
r
q
m
∫
0
1
q
3
(
ν
)
d
ν
∫
0
∞
q
4
(
z
)
d
z
·
∫
S
1
(
0
)
(
1
-
ν
2
/
m
exp
⁡
(
-
2
z
m
+
2
)
)
×
f
i
(
y
1
(
i
)
(
z
,
ν
,
w
)
,
τ
1
(
z
,
ν
)
)
d
S
w
=
r
q
m
E
{
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
×
f
i
(
y
1
(
i
)
(
ξ
1
,
v
,
ω
)
,
τ
1
(
ξ
1
,
v
)
)
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
}
,
where q
3
(
ν
)
=
(
1
-
ν
)
ν
2
/
m
-
1
(
B
(
2
,
2
/
m
)
)
-
1. Then the density of Beta distribution with parameters (
2
,
2
/
m
)
(43)
q
4
(
z
)
=
exp
⁡
(
-
z
)
z
m
/
2
-
1
(
Γ
(
m
2
)
)
-
1
,
and the density of Gamma distribution with parameters m
/
2,ω—unit random vector,
(44)
y
1
i
(
ξ
1
,
v
,
ω
)

=

x
+
[
4
m
m

+

2
r
a
i
ξ
1
v
2
/
m
exp
⁡
(
-
2
ξ
1
m

+
2
)
]
1
/
2
ω
,
τ
1
(
ξ
1
,
v
)
=
t
-
r
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
,
where ξ
1 is the random variable with density function q
4
(
z
) and v is another random variable with the density function q
3
(
ν
).Let(45)
r
=
r
(
x
,
t
)
=
min
⁡
{
e
R
1
2
(
x
)
2
m
a
1
,
…
,
e
R
1
2
(
x
)
2
m
a
n
;
1
c
11
,
…
,
1
c
n
n
;
t
}
,
then Q
r
(
i
)
¯
(
x
,
t
)
∈
Ω and the function
(46)
p
1
(
i
)
(
x
,
t
;
y
,
τ
)
=
[
1
-
(
r
-
(
t
-
τ
)
)
c
i
]
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
r
q
m
(
1
-
r
q
1
m
c
i
i
)
×
I
Q
r
(
i
)
(
x
,
t
)
(
y
,
τ
)
is the transition density in Q
r
i
(
x
,
t
) with fixed point (
x
,
t
), where
(47)
q
1
m
=
1
-
1
2
(
m
+
2
m
+
4
)
1
+
m
/
2
.
Let (
y
2
(
i
)
,
τ
2
(
i
)
) be a random point of balloid Q
r
i
(
x
,
t
) which has the following density function p
1
(
i
)
(
x
,
t
;
y
,
τ
)
(
i
=
1
,
…
,
n
) in the fixed point (
x
,
t
).In this case(48)
1
r
∬
Q
r
(
i
)
(
x
,
t
)
(
1
-
(
r
-
(
t
-
τ
)
)
c
i
i
)
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
u
i
(
y
,
τ
)
d
y
d
τ
=
q
m
(
1
-
r
c
i
i
q
1
m
)
E
u
i
(
y
2
(
i
)
,
τ
2
(
i
)
)
,
i
=
1
,
.
.
.
,
n
.
The obtained results we will put to (38) and we will get the probabilistic representation of problem (24)-(25). It follows there from that we could to following proposition.Theorem 3.
For the solution of initial BVP (24)-(25) the following probabilistic representation is valid:
(49)
u
i
(
x
,
t
)
=
(
1
-
q
m
)
E
u
i
(
y
(
i
)
(
ξ
,
ω
)
,
τ
(
ξ
)
)
+
q
m
(
1
-
r
c
i
i
q
1
m
)
E
u
i
(
y
2
(
i
)
,
τ
2
(
i
)
)
+
q
m
r
E
{
∑
j
=
1
,
j
≠
i
n
c
i
j
u
j
(
y
1
(
i
)
(
ξ
1
,
v
,
ω
,
)
,
τ
1
(
ξ
1
,
v
)
)
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
×
∑
j
=
1
,
j
≠
i
n
c
i
j
u
j
(
y
1
(
i
)
(
ξ
1
,
v
,
ω
,
)
,
τ
1
(
ξ
1
,
v
)
)
}
+
q
m
r
E
(
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
aaaaaaaaaa
×
f
i
(
y
1
(
i
)
(
ξ
1
,
v
,
ω
)
,
τ
1
(
ξ
1
,
v
)
)
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
1
m
+
2
)
)
)
,
(
i
=
1
,
n
¯
)
,
where (
y
(
i
)
(
ξ
,
ω
)
,
τ
(
ξ
)
) is defined by (40) and (
y
(
i
)
(
ξ
,
v
,
ω
)
,
τ
(
ξ
,
v
)
) is determined by (41).The proof of Theorem3 is the consequence of the above mentioned reasoning.By further using the presentation (44) we will construct a random process in Ω and propose the Monte Carlo algorithm for the solution of system IEs.


## 3.3. Description Random Process and the Algorithm Simulation

Let(50)
r
=
r
(
x
,
t
)
=
min
⁡
{
e
R
1
2
(
x
)
2
m
a
1
,
…
,
e
R
1
2
(
x
)
2
m
a
n
;
1
c
11
,
…
,
1
c
n
n
;
t
}
,
The functions
(51)
p
0
(
i
)
(
x
,
t
;
y
,
τ
)
=
1
1
-
q
m
(
1
-
t
-
τ
r
)
(
-
∂
Z
(
i
)
(
x
,
t
;
y
,
τ
)
∂
n
y
)
×
I
∂
Q
r
(
i
)
(
x
,
t
)
(
y
,
τ
)
,
p
1
(
i
)
(
x
,
t
;
y
,
τ
)
=
(
1
-
(
r
-
(
t
-
τ
)
)
c
i
i
)
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
r
q
m
(
1
-
r
q
1
m
c
i
i
)
×
I
∂
Q
r
(
i
)
(
x
,
t
)
(
y
,
τ
)
,
p
2
(
i
)
(
x
,
t
;
y
,
τ
)
=
Z
r
(
i
)
(
x
,
t
;
y
,
τ
)
r
q
m
I
∂
Q
r
(
i
)
(
x
,
t
)
(
y
,
τ
)
,
are the transition density functions in Q
r
(
i
)
¯
(
x
,
t
) at a fixed point (
x
,
t
)
(
i
=
1
,
n
¯
). We will define in Ω a random process as was proposed in Section 2.Let us define a transition matrix as(52)
A
(
x
,
t
)
=
[
α
11
α
12
…
α
1
(
n
+
1
)
α
21
α
22
…
α
2
(
n
+
1
)
…
…
…
…
α
n
1
α
n
2
…
α
n
(
n
+
1
)
0
0
…
1
]
,
where α
i
i
=
1
-
q
m
c
i
i
q
1
m
r
(
x
,
t
), and let
(53)
β
i
=
q
m
c
i
i
q
1
m
r
(
x
,
t
)
(
n
-
1
)
n
,
α
i
j
=
β
i
c
i
j
M
i
,
(
i
,
j
=
1
,
n
¯
;
i
≠
j
)
,
M
i
=
∑
j
=
1
,
n
;
j
≠
i
c
i
j
,
(
i
=
1
,
n
¯
)
,
α
i
(
n
+
1
)
=
q
m
q
1
m
r
(
x
,
t
)
c
i
i
n
,
(
i
=
1
,
n
¯
)
.
Now we will define the density function of transition matrix P
(
x
,
t
;
y
,
τ
):
(54)
P
(
x
,
t
;
y
,
τ
)
=
[
p
11
(
x
,
t
;
y
,
τ
)
p
12
(
x
,
t
;
y
,
τ
)
.
.
.
p
1
(
n
+
1
)
(
x
,
t
;
y
,
τ
)
p
21
(
x
,
t
;
y
,
τ
)
p
22
(
x
,
t
;
y
,
τ
)
.
.
.
p
2
(
n
+
1
)
(
x
,
t
;
y
,
τ
)
.
.
.
.
.
.
.
.
.
.
.
.
p
n
1
(
x
,
t
;
y
,
τ
)
p
n
2
(
x
,
t
;
y
,
τ
)
.
.
.
p
n
(
n
+
1
)
(
x
,
t
;
y
,
τ
)
0
0
.
.
.
1
]
,
where
(55)
p
i
j
(
x
,
t
;
y
,
τ
)
=
p
2
(
i
)
(
x
,
t
;
y
,
τ
)
,
(
i
=
1
,
n
¯
;
j
=
1
,
n
+
1
¯
;
i
≠
j
)
,
p
i
i
(
x
,
t
;
y
,
τ
)
=
(
1
-
q
m
)
p
0
(
i
)
(
x
,
t
;
y
,
τ
)
+
q
m
(
1
-
r
q
1
m
c
i
i
)
p
1
(
i
)
(
x
,
t
;
y
,
τ
)
1
-
r
(
x
,
t
)
q
1
m
c
i
i
q
m
.
Then we will fix the initial point (
x
0
,
t
0
)
=
(
x
,
t
) and the number of equations i
0
∈
{
1
,
…
,
n
}. Let an initial moment at the point (
x
0
,
t
0
)
=
(
x
,
t
); we will have one particle. For one step a particlei
k
→
i
k
+
1 moves from its position according to the transition matrix A
(
x
k
,
t
k
) and moves with probability α
i
k
,
i
k
+
1
(
x
k
,
t
k
) from the point (
x
k
,
t
k
) to the point (
x
k
+
1
,
t
k
+
1
). The next (
x
k
+
1
,
t
k
+
1
) point will be simulated using the density function p
i
k
i
k
+
1
(
x
k
,
t
k
;
y
,
τ
).The probability of breaking of trajectory in the point(
x
n
,
t
n
) is
(56)
g
(
x
n
,
t
n
)
=
{
1
,
(
x
n
,
t
n
)
∈
Ω
;
α
i
n
,
n
+
1
(
x
n
-
1
,
t
n
-
1
)
,
(
x
n
,
t
n
)
∈
Ω
.
The next coordinate of the particle will be defined in the following way.(1) If the density function of the point(
x
n
+
1
,
t
n
+
1
) equals p
0
(
i
)
(
x
n
,
t
n
;
y
,
τ
) in the fixed point (
x
n
,
t
n
) then
(57)
x
n
+
1
=
x
n
+
2
(
r
(
x
n
,
t
n
)
ξ
n
a
i
)
1
/
2
exp
⁡
(
-
ξ
n
m
)
ω
n
,
t
n
+
1
=
t
n
-
r
(
x
n
,
t
n
)
exp
⁡
(
-
2
ξ
n
m
)
,
where {
ξ
n
}
n
=
0
∞
,
{
ω
n
}
n
=
0
∞ the sequence of independent random variables with the density function q
1
(
ρ
) and independent isotropic vectors. The value r
(
x
n
,
t
n
) will be defined as (50).(2) If the density function of the point(
x
n
+
1
,
t
n
+
1
) is equal to p
1
(
i
)
(
x
n
,
t
n
;
y
,
τ
) at a fixed (
x
n
,
t
n
) then
(58)
x
n
+
1
=
x
n
+
2
(
(
-
2
ξ
n
′
m
+
2
)
m
m
+
2
r
(
x
n
,
t
n
)
×
ξ
n
′
(
v
n
′
)
2
/
m
a
i
exp
⁡
⁡
(
-
2
ξ
n
′
m
+
2
)
)
1
/
2
ω
n
,
t
n
+
1
=
t
n
-
r
(
x
n
,
t
n
)
(
v
n
′
)
2
/
m
exp
⁡
(
-
2
ξ
n
′
m
+
2
)
,
where {
ξ
n
′
}
n
=
0
∞
,
{
v
n
′
}
n
=
0
∞ is a sequence of independent random variables, which will be obtained from the algorithm below (Algorithm 4) (Neumann acceptance rejection method).Algorithm 4.
(a) We firstly simulateξ, Gamma distributed random variable with the parameters (
m
/
2
), secondly simulate γ, uniformly distributed random variable on (0, 1), and thirdly simulate ν, Beta distributed random variable with the parameters (2
,
2
/
m).
(b) Ifγ
>
1
-
c
i
i
r
(
1
-
v
2
/
m
exp
⁡
(
-
2
ξ
/
(
m
+
2
)
)
) then we will go to (a) and so on; otherwise v
′
=
v
,
ξ
′
=
ξ.
(3) If the density function at the point(
x
n
+
1
,
t
n
+
1
) equals p
2
(
i
)
(
x
n
,
t
n
;
y
,
τ
) under fixed point (
x
n
,
t
n
), then
(59)
x
n
+
1
=
x
n
+
2
(
(
-
2
ξ
n
m
+
2
)
m
m
+
2
r
(
x
n
,
t
n
)
ξ
n
(
v
n
)
2
/
m
×
a
i
exp
⁡
(
-
2
ξ
n
m
+
2
)
)
1
/
2
ω
n
,
t
n
+
1
=
t
n
-
r
(
x
n
,
t
n
)
(
v
n
)
2
/
m
exp
⁡
(
-
2
ξ
n
m
+
2
)
,
where {
ξ
n
}
n
=
0
∞
,
{
v
n
}
n
=
0
∞
,
{
ω
n
}
n
=
0
∞ is sequence of independent Gamma distributed random variables with parameters (
2
/
m
), Beta distributed random variables with parameters (
2
,
2
/
m
), and independent isotropic vectors, respectively.If at the momentn was held break, then we will put (
x
n
+
k
,
t
n
+
k
)
=
(
x
n
,
t
n
)
,
k
=
0,1
,
2
,
.
.. obviously the sequence of coordinate of the particle forms Markov chains. The random process which was described above was considered in Ermakov et al. [23] for the solution of initial BVP for the heat equation and adapted in Kurbanmuradov [22] for the heat equation with variable coefficients.Now we prove the auxiliary Lemma5.Lemma 5.
With the probability one Markov chain{
x
n
,
t
n
}
n
=
0
∞ converges when n
→
∞ to the random point of boundary (
x
∞
,
t
∞
)
∈
∂
Ω, or it is absorbed inside of the domain.Proof.
Since{
t
n
} is decreasing sequence and t
n
≥
0, it has a limit t
∞
=
lim
⁡
n
→
∞
t
n. Let R
n
-
σ be algebra, which was generated by random variables
(60)
{
ω
k
}
k
=
0
n
-
1
,
{
v
k
}
k
=
0
n
-
1
,
{
ξ
k
}
k
=
0
n
-
1
,
{
v
k
′
}
k
=
0
n
-
1
,
{
ξ
k
′
}
k
=
0
n
-
1
.
From the definition R
n and (57)–(59) it follows that x
n is measurable relatively R
n. It is obvious that the coordinates of vector process formed limited martingale relatively {
R
m
}
m
=
1
∞:
(61)
E
(
x
n
+
1
R
n
)
=
E
{
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
x
n
R
n
(
1
-
q
m
)
×
[
x
n
+
2
(
r
(
x
n
,
t
n
)
ξ
n
a
i
n
)
1
/
2
exp
⁡
(
-
ξ
n
m
)
ω
n
]
+
q
m
(
1
-
q
1
m
r
(
x
n
,
t
n
)
c
i
n
i
n
)
×
[
exp
⁡
(
-
ξ
n
′
m
+
2
)
ω
n
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
ξ
n
′
v
n
′
2
/
m
)
1
/
2
×
exp
⁡
(
-
ξ
n
′
m
+
2
)
ω
n
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
-
1
n
×
[
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
ς
n
v
n
2
/
m
)
1
/
2
×
exp
⁡
(
-
ς
m
+
2
)
ω
n
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
ς
n
v
n
2
/
m
)
1
/
2
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
x
n
R
n
}
=
(
1
-
q
m
)
[
(
ξ
n
1
/
2
exp
⁡
(
-
ξ
n
m
)
ω
n
)
x
n
+
2
(
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
×
E
(
ξ
n
1
/
2
exp
⁡
(
-
ξ
n
m
)
ω
n
)
]
+
q
m
(
1
-
q
1
m
c
i
n
i
n
r
(
x
n
,
t
n
)
)
×
[
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
×
E
(
(
ξ
n
′
v
n
′
2
/
m
)
exp
⁡
(
-
ξ
n
m
+
2
)
ω
n
)
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
-
1
n
×
[
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
×
E
(
(
ξ
n
′
v
n
′
2
/
m
)
exp
⁡
(
-
ξ
n
m
+
2
)
ω
n
)
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
-
1
n
×
[
x
n
+
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
×
E
(
(
ς
n
v
n
2
/
m
)
1
/
2
exp
⁡
(
-
ς
n
m
+
2
)
ω
n
)
2
(
m
m
+
2
r
(
x
n
,
t
n
)
a
i
n
)
1
/
2
]
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
n
x
n
=
(
1
-
q
m
)
x
n
+
q
m
(
1
-
q
1
m
c
i
n
i
n
r
(
x
n
,
t
n
)
)
x
n
+
c
i
n
i
n
r
(
x
n
,
t
n
)
q
m
q
1
m
x
n
=
x
n
,
where {
x
n
} are limited martingale; it converges with the probability one Shiryaev [24].
Let(
x
∞
,
t
∞
)
=
lim
⁡
n
→
∞
(
x
n
,
t
n
) be the limit vector. We show that (
x
∞
,
t
∞
)
∈
∂
Ω. If t
∞
=
0 then the process is broken inside the domain. Let t
∞
>
0. As far as the process converges, according to the formulas (50)–(58) we have
(62)
E
(
x
0
,
t
0
)
|
x
n
+
1
-
x
n
|
⟶
0
,
E
(
x
0
,
t
0
)
|
x
n
+
1
-
x
n
|
=
E
(
x
0
,
t
0
)
{
r
(
x
n
,
t
n
)
h
(
r
(
x
n
,
t
n
)
)
}
,
where h
(
r
) is strictly positive. Applying Lebesque Theorem (about the limited convergence) we get
(63)
E
(
x
0
,
t
0
)
(
r
(
x
∞
,
t
∞
)
)
1
/
2
=
0
.
It means r
(
x
∞
,
t
∞
)
=
0. Then from the definition of r
(
x
,
t
) and using the formulae (49) we obtain
(64)
R
1
(
x
∞
)
=
0
,
(
x
∞
,
t
∞
)
∈
∂
Ω
.
Lemma is proven.


## 4. Construction Unbiased andɛ-Biased Estimators

Let(
x
k
,
t
k
)
k
=
0
∞ be the trajectories of random process which was described above. We will define on it the sequence of the random variables {
η
n
(
i
0
)
}
n
=
0
∞. Let
(65)
Θ
0
=
1
,
Θ
n
=
Θ
n
-
1
×
V
i
n
-
1
i
n
(
x
n
-
1
,
t
n
-
1
;
x
n
,
t
n
)
,
where V
i
j
(
x
n
-
1
,
t
n
-
1
;
x
n
,
t
n
) is defined as follows:
(66)
V
i
j
(
x
n
-
1
,
t
n
-
1
;
x
n
,
t
n
)
=
n
M
i
(
n
-
1
)
c
i
i
q
1
m
(
1
-
v
n
2
/
m
exp
⁡
(
-
2
ς
n
m
+
2
)
)
,
(
i
,
j
=
1
,
n
¯
;
i
≠
j
)
,
V
i
i
(
x
n
-
1
,
t
n
-
1
;
x
n
,
t
n
)
=
1
;
(
i
=
1
,
n
¯
)
;
V
i
(
n
+
1
)
(
x
n
-
1
,
t
n
-
1
;
x
n
,
t
n
)
=
n
c
i
i
q
1
m
(
1
-
v
n
2
/
m
exp
⁡
(
-
2
ς
n
m
+
2
)
)
;
(
i
=
1
,
n
¯
)
.
Here {
ς
n
}
n
=
0
∞
,
{
v
n
}
n
=
0
∞ are the sequences of independent Gamma function with the parameters (m
/
2) and Beta function with parameters (2,2
/
m) distributed random variables, respectively. We will define the sequence
(67)
η
n
(
i
0
)
=
Θ
n
×
F
(
x
n
,
t
n
)
=
Θ
n
{
u
j
(
x
n
,
t
n
)
,
i
n
=
j
,
j
≠
n
+
1
,
f
i
n
-
1
(
x
n
,
t
n
)
,
i
n
=
n
+
1
.
If at the moment n happen break n, we will put
(68)
η
n
+
k
(
i
0
)
=
η
n
(
i
0
)
,
(
x
n
+
k
,
t
n
+
k
)
=
(
x
n
,
t
n
)
,
k
=
1,2
,
…
,
where algebra R
n
-
σ generated until the moment n.Theorem 6.
Let the sequence be form martingale{
η
n
(
i
0
)
}
n
=
1
∞ with R
n, respectively. If
(69)
∑
j
=
1
,
…
,
n
;
j
≠
i
c
i
j
<
(
n
-
1
)
c
i
i
q
1
m
n
,
(
i
=
1
,
n
¯
)
,
max
⁡
(
x
,
t
)
∈
Ω
|
f
i
(
x
,
t
)
|
≤
c
0
,
(
c
0
=
c
o
n
s
t
,
i
=
1
,
n
¯
)
.
Then the sequence will be {
η
n
(
i
0
)
} uniformly integrable martingale.Proof.
From the definitionη
n
(
i
0
) the R
n is measurable. In this case
(70)
E
(
η
n
+
1
(
i
0
)
R
n
)
=
E
(
Θ
n
+
1
×
F
(
x
n
+
1
,
t
n
+
1
)
R
n
)
=
E
(
Θ
n
×
V
i
n
i
n
+
1
(
x
n
,
t
n
;
x
n
+
1
,
t
n
+
1
)
×
F
(
x
n
+
1
,
t
n
+
1
)
R
n
)
=
Θ
n
E
(
V
i
n
i
n
+
1
(
x
n
,
t
n
;
x
n
+
1
,
t
n
+
1
)
×
F
(
x
n
+
1
,
t
n
+
1
)
)
=
Θ
n
[
∑
j
=
1
,
n
α
i
n
j
∬
Q
r
(
i
n
)
¯
(
x
n
,
t
n
)
p
i
n
j
(
x
n
,
t
n
;
y
,
τ
)
u
j
(
y
,
τ
)
d
y
d
τ
+
α
i
n
n
+
1
∬
Q
r
(
i
n
)
¯
(
x
n
,
t
n
)
p
i
n
n
+

1
(
x
n
,
t
n
;
y
,
τ
)
f
i
n
(
y
,
τ
)
d
y
d
τ
∑
j
=
1
,
n
α
i
n
j
∬
Q
r
(
i
n
)
¯
(
x
n
,
t
n
)
p
i
n
j
(
x
n
,
t
n
;
y
,
τ
)
u
j
(
y
,
τ
)
d
y
d
τ
]
=
η
n
(
i
0
)
.
As far as the sequence is martingale {
η
n
(
i
0
)
}. We can show the uniformly integrability of η
n
(
i
0
). To do that it is enough to show |
η
n
(
i
0
)
|
<
∞.
Since(71)
u
i
(
x
,
t
)
∈
C
(
D
¯
×
[
0
,
T
]
)
∩
C
2,1
(
D
-
×
[
0
,
T
]
)
,
and Ω is bounded domain |
u
i
(
x
,
t
)
|
≤
const, for any (
x
,
t
)
∈
Ω. From the condition of theorem |
Θ
n
|
≤
1, it is followed |
η
n
(
i
0
)
|
≤
const. It means {
η
n
(
i
0
)
} is uniformly integrable.
The theorem is proved.Now we will construct computable (realizable) estimatorη
n
(
i
0
). We will take ε-neighborhoods of the domain (
∂
Ω
)
ε
=
{
D
×
[
0
,
ε
]
}
∪
{
(
∂
D
)
ε
×
[
0
,
T
]
}.LetN
1 be a breaking moment of process inside of domain. N
ε is the first passage moment (
∂
Ω
)
ε. N
=
min
⁡
{
N
1
,
N
ε
}, stopping moment of process {
(
x
n
,
t
n
)
}. In this case the probability of absorbing will be
(72)
g
(
x
n
,
t
n
)
=
{
1
,
(
x
n
,
t
n
)
∈
(
∂
Ω
)
ε
α
i
n
,
(
n
+
1
)
(
x
n
-
1
,
t
n
-
1
)
,
(
x
n
,
t
n
)
∈
Ω
¯
∖
(
∂
Ω
)
ε
.
From Lemma 5 it follows that N
<
∞. It could be proved that the mathematical expectation of the stopping time {
(
x
n
,
t
n
)
} of Markov process is finite.Theorem 7.
Let the condition of Theorem6 be satisfied; then the estimator η
n
(
i
0
) will be unbiased estimator with finite variance, where u
i
0
(
x
,
t
) is i
0th component of solution vector u
(
x
,
t
).Proof.
Sinceη
n
(
i
0
) is uniformly integrable martingale and N is Markov moment, according to the theorem in ([24, 25]) for the martingale {
η
n
(
i
0
)
} we obtain
(73)
E
(
x
0
,
t
0
)
η
N
(
i
0
)
=
E
(
x
0
,
t
0
)
η
1
(
i
0
)
.
From the definition η
n
(
i
0
) holds E
(
x
0
,
t
0
)
η
1
(
i
0
)
=
u
i
0
(
x
,
t
). From condition of Theorem 6  E
(
η
N
)
2
(
i
0
)
<
∞ is valid accordingly the variance is finite. The theorem is proved.Further, fromη
N
(
i
0
) we could construct biased but computable (realizable) estimator η
N
*
(
i
0
). Let, for x
∈
∂
D, t
∈
[
0
,
T
] and ψ
i
(
x
,
t
)
=
y
i
(
x
,
t
) for x
∈
D, ψ
i
(
x
,
0
)
=
y
0
i
(
x
)
,
(
x
*
,
t
*
) closed to the point (
x
,
t
) of boundary ∂
Ω. η
N
* will be obtained with changing
(74)
u
i
(
x
N
,
t
N
)
in
η
N
(
i
0
)
to
ψ
i
(
x
N
*
,
t
N
*
)
.
Let us evaluate bias η
N
*
(
i
0
). It is clear that
(75)
|
E
(
x
,
t
)
η
N
(
i
0
)
-
u
i
0
(
x
,
t
)
|
≤
E
(
x
,
t
)
|
η
N
*
(
i
0
)
-
η
N
(
i
0
)
|
.
If N
=
N
1, in this case the process is broken when do not reach the boundary (
∂
Ω
)
ε and η
N
*
(
i
0
)
=
η
N
(
i
0
). If N
=
N
ε then (
x
N
,
t
N
)
∈
(
∂
Ω
)
ε.LetA
i
(
ε
) be a module of continuity function u
i
(
x
,
t
). In this case it is true:
(76)
|
η
N
*
(
i
0
)
-
η
N
(
i
0
)
|
≤
|
Θ
N
|
A
(
ε
)
,
where A
(
ε
)
=
max
⁡
i
A
i
(
ε
), since |
Θ
N
|
≤
1 then E
(
x
,
t
)
|
η
N
*
(
i
0
)
-
u
i
0
(
x
,
t
)
|
≤
A
(
ε
). Finiteness of variance followed from E
η
N
2
(
i
0
)
<
∞. The proposed algorithm we could generalize for the case with variable coefficients c
i
j
=
c
i
j
(
x
,
t
) and one could get the same results.


## 5. Computational Example

LetD
∈
R
3 be bounded domains, Ω
=
D
×
[
0
,
T
]. We will consider for some mode linitial boundary value problem
(77)
∂
u
i
(
x
,
t
)
∂
t
-
a
i
Δ
u
i
(
x
,
t
)
+
c
i
i
u
i
(
x
,
t
)
-
∑
j
=
1,4
;
j
≠
i
c
i
j
u
j
(
x
,
t
)
=
f
i
(
x
,
t
)
,
(
i
=
1,4
¯
)
for (
x
,
t
)
∈
Ω with the initial boundary conditions
(78)
u
i
(
x
,
t
)
=
y
i
(
x
,
t
)
,
x
∈
∂
D
,
t
∈
[
0
,
T
]
,
i
=
1,4
¯
,
u
i
(
x
,
0
)
=
y
0
i
(
x
)
,
x
∈
D
,
i
=
1,4
¯
.
As domain is chosen as the simple ball, D
=
{
(
x
1
,
x
2
,
x
3
)
:
x
1
2
+
x
2
2
+
x
3
2
≤
R
2
}.The coefficients(79)
(
a
1
a
2
a
3
a
4
)
=
(
0.5
0.7
0.1
1.0
)
,
{
c
i
j
}
i
,
j
=
1
,
.
.
.
,
4
=
[
2
0.4
0.5
0.2
0.7
3
0.4
0.6
0.3
0.1
1
0.1
0.2
0.3
0.3
1.5
]
The initial and boundary conditions
(80)
y
01
(
x
)
=
x
1
2
+
x
2
2
+
x
3
2
;
y
1
(
x
,
t
)
=
R
2
exp
⁡
⁡
(
t
)
;
y
02
(
x
)
=
(
x
1
x
2
x
3
)
2
;
y
2
(
x
,
t
)
=
exp
⁡
⁡
(
t
)
(
x
2
x
3
)
2
(
R
2
-
x
2
2
-
x
3
2
)
;
y
03
(
x
)
=
exp
⁡
⁡
(
x
1
+
x
2
+
x
3
)
;
y
3
(
x
,
t
)
=
exp
⁡
⁡
(
t
+
x
1
+
x
2
+
x
3
)
;
y
04
(
x
)
=
1
;
y
4
(
x
,
t
)
=
exp
⁡
(
t
x
1
x
2
x
3
)
.
Left hand sides
(81)
f
1
(
x
,
t
)
=
exp
⁡
⁡
(
t
)
[
4
(
x
1
2
+
x
2
2
+
x
3
2
)
-
0.3
(
x
1
x
2
x
3
)
2
-
3
]
-
0.2
(
x
1
+
x
2
+
x
3
+
t
)
-
0.5
exp
⁡
(
x
1
x
2
x
3
t
)
;
f
2
(
x
,
t
)
=
exp
⁡
⁡
(
t
)
[
3.5
(
x
1
x
2
x
3
)
2
-
1.4
(
x
1
x
2
)
2
+
(
x
2
x
3
)
2
+
(
x
1
x
3
)
2
-
0.4
(
x
1
2
+
x
2
2
+
x
3
2
)
]
-
0.3
(
x
1
+
x
2
+
x
3
+
t
)
-
0.2
exp
⁡
(
x
1
x
2
x
3
t
)
;
f
3
(
x
,
t
)
=
4.7
exp
⁡
(
x
1
+
x
2
+
x
3
+
t
)
-
exp
⁡
(
t
)
[
0.5
(
x
1
2
+
x
2
2
+
x
3
2
)
+
0.4
(
x
1
x
2
x
3
)
2
]
-
0.7
exp
⁡
(
x
1
x
2
x
3
t
)
;
f
4
(
x
,
t
)
=
exp
⁡
(
x
1
x
2
x
3
t
)
×
[
(
x
1
x
2
x
3
)
+
3.5
-
t
2
(
(
x
1
x
2
)
2
+
(
x
2
x
3
)
2
+
(
x
1
x
3
)
2
)
]
-
exp
⁡
⁡
(
t
)
[
0.3
(
x
1
2
+
x
2
2
+
x
3
2
)
-
0.4
(
x
1
x
2
x
3
)
2
]
-
0.6
(
x
1
+
x
2
+
x
3
+
t
)
.
The exact solutions are known:
(82)
u
1
(
x
,
t
)
=
exp
⁡
⁡
(
t
)
(
x
1
2
+
x
2
2
+
x
3
2
)
,
u
2
(
x
,
t
)
=
exp
⁡
⁡
(
t
)
(
x
1
x
2
x
3
)
2
,
u
3
(
x
,
t
)
=
exp
⁡
⁡
(
x
1
+
x
2
+
x
3
+
t
)
,
u
4
(
x
,
t
)
=
exp
⁡
(
x
1
x
2
x
3
t
)
.


## 6. Conclusions

It is known that the distinguishing feature of the vector algorithm is that its “weight” appears to be a matrix weight. This matrix weight is multiplied by the kernel matrix of the system of integral equations divided by a transition density function after each transition in the Markov chain simulation. In this case the computational complexity is higher enough than simple Monte Carlo method. On the contrary to the vector algorithms we proposed a new Monte Carlo algorithm for the solution of system of integral equations. This method has the simple structure of the computation algorithm and the errors do not depend on the dimension of domain and smoothness of boundary. One can solve the problem at one point and we do not use matrix weights. Proposed algorithm applied to the solution of system of the parabolic equations. To do so we derived corresponding system of integral equations and construct a special probabilistic representation. This probabilistic representation uses for simulation the random process and construction the unbiased andε-biased estimator for the solution of systems IEs.Numerical experiments show that the computational complexity of our algorithm is reduced. In the future the proposed algorithm might be generalized for the case with variable coefficientsc
i
j
=
c
i
j
(
x
,
t
). The results of numerical experiments are shown with the probability almost one; the approximate solution tends to the exact solution of the problem. In the given example the exact solution is known; therefore we can make sure that all the estimators really are in the confidence intervals (see Table 1).Table 1
The results of computational experiments.
(
x
0
,
t
0
)
i
0
N
t
ε
U
e
MC
3
σ
err
(
0.1
;
-
0.8
;
0.4
;
1.2
)
1
5000
0.005
2.689294
2.881862
0.2821
0.1925
(
0.1
;
-
0.8
;
0.4
;
1.2
)
2
5000
0.005
0.003399
0.018739
0.224699
0.015339
(
0.1
;
-
0.8
;
0.4
;
1.2
)
3
5000
0.005
2.459603
2.483835
0.21815
0.02423
(
0.1
;
-
0.8
;
0.4
;
1.2
)
4
5000
0.005
0.962327
0.95895
0.03377
0.1598
(
-
0.5
;
0.25
;
-
0.4
;
1.3
)
1
5000
0.005
1.733942
1.554352
0.32319
0.17939
(
-
0.5
;
0.25
;
-
0.4
;
1.3
)
2
5000
0.005
0.009173
0.009841
0.215628
0.000668
(
-
0.5
;
0.25
;
-
0.4
;
1.3
)
3
5000
0.005
1.91554
1.980669
0.18057
0.064528
(
-
0.5
;
0.25
;
-
0.4
;
1.3
)
4
5000
0.005
1.067159
1.084443
0.162095
0.01728
(
x
0
;
t
0
) is the point which solved BVP; N
t is the quantity of samples (trajectories); ε is neighborhood area; i
0 is the number of equation; U
e is exact solution at the point (
x
0
;
t
0
). MC is Monte Carlo solutions; 3
σ is confidence interval; err is the difference between exact ant MC solutions err
⁡
=
|
u
i
0
(
x
0
;
t
0
)
-
MS
|.



---

*Source: 236065-2014-08-14.xml*