I T Dimov
email: ivdimov@bas.bg

S Maire
email: maire@univ-tln.fr

R Georgieva

A Study of Biased and Unbiased Stochastic Algorithms for Solving Integral Equations

In this paper we propose and analyse a new unbiased stochastic method for solving a class of integral equations, namely the second kind Fredholm integral equations. We study and compare three possible approaches to compute linear functionals of the integral under consideration: i) biased Monte Carlo method based on evaluation of truncated Liouville-Neumann series, ii) transformation of this problem into the problem of computing a finite number of integrals, and iii) unbiased stochastic approach. Five Monte Carlo algorithms for numerical integration have been applied for approach (ii). Error balancing of both stochastic and systematic errors has been discussed and applied during the numerical implementation of the biased algorithms. Extensive numerical experiments have been performed to support the theoretical studies regarding the convergence rate of Monte Carlo methods for numerical integration done in our previous studies. We compare the results obtained by some of the best biased stochastic approaches with the results obtained by the proposed unbiased approach. Conclusions about the applicability and efficiency of the algorithms have been drawn.

Introduction

Integral equations are of high applicability in different areas of applied mathematics, physics, and engineering. In particular, they are widely used in mechanics, geophysics, electricity and magnetism, kinetic theory of gases, quantum mechanics, mathematical economics, and queuing theory. That is why it is reasonable to develop and study efficient and reliable approaches to solve integral equations.

An important advantage of Monte Carlo (MC) methods is that they allow to find directly an unknown linear functional of the solution of integral equations with a number of operations necessary to calculate the solution of an integral equation only at one point of the domain [START_REF] Sobol | Monte Carlo Numerical Methods[END_REF]. The existing MC methods for integral equations (MCM-int eq.) are based on probabilistic representations of the Liouville-Neumann series for the second kind Fredholm integral equation [START_REF] Curtiss | Monte Carlo methods for the iteration of linear operators[END_REF][START_REF] Dimov | Monte Carlo methods for applied scientists[END_REF].

The possible unbiased approaches deal with infinite series, while the biased MC approaches use probabilistic representations of truncated Liouville-Neumann series. A well known and widely used biased method is the Markov chain MC (see, for example [START_REF] Kalos | Monte Carlo Methods[END_REF]). Usually, the Markov chain stops after a fixed number of steps.

The unbiased stochastic algorithm proposed in this work is actually based on a probabilistic representation of the infinite series introducing some probability for absorbtion. The idea of the construction of the algorithm is similar to the one used in [START_REF] Dimov | A New Walk on Equations Monte Carlo Method for Linear Algebraic Problems[END_REF] for solving linear systems.

A possible approach to deal with the problem of approximation of linear functionals of the solution of an integral equation is to transform it into approximate evaluation of finite number of integrals (FNI) (linear functionals of iterative functions) [START_REF] Dimov | Efficient and overconvergent Monte Carlo Methods. Parallel algorithms[END_REF][START_REF] Dimov | Monte Carlo methods for applied scientists[END_REF]. We have already applied a variance reduction approach called importance separation method for each of integrals and extended the study of its properties for solving integral equations (see [START_REF] Georgieva | Importance Separation for Solving Integral Equations[END_REF]). Here we extend the study of the properties to five different MC algorithms for multidimensional numerical integration for solving integral equations. These algorithms are: Crude MC algorithm [START_REF] Sobol | Monte Carlo Numerical Methods[END_REF], based on SIMD-oriented fast Mersenne Twister pseudo-random number generator [START_REF] Saito | SIMD-oriented fast Mersenne Twister: a 128-bit pseudorandom number generator[END_REF]18], quasi-Monte Carlo (qMC) algorithm based on ΛΠ τ Sobol quasirandom sequences, MC algorithms (MCA-MSS-1 and MCA-MSS-2) based on modified Sobol quasirandom sequences [START_REF] Dimov | Multidimensional Sensitivity Analysis of Large-scale Mathematical Models[END_REF][START_REF] Dimov | Monte Carlo method for numerical integration based on Sobol' sequences[END_REF][START_REF] Dimov | Advanced algorithms for multidimensional sensitivity studies of large-scale air pollution models based on Sobol sequences, Computers and Mathematics with Applications 65[END_REF] and a stratified symmetrised MC algorithm (MCA-MSS-2-S) [START_REF] Dimov | Multidimensional Sensitivity Analysis of Large-scale Mathematical Models[END_REF].

The paper is organized as follows. The problem setting is given in Section 2. A description of biased stochastic approaches applied to solve the problem under consideration is presented in Section 3. In this section some error analysis based on the concept of balancing both stochastic and systematic errors is done. The new unbiased stochastic approach is presented in Section 4. Numerical results are described and analysed in Section 5. Some final conclusions are drawn in Section 6.

Formulation of the Problem

Consider the second kind Fredholm integral equation:

u(x) = ∫ Ω k(x, x ′ )u(x ′ )dx ′ + f (x) (1)
or u = Ku + f (K is an integral operator), where k(x, x ′ ) ∈ L 2 (Ω × Ω), f (x) ∈ L 2 (Ω) are given functions and u(x) ∈ L 2 (Ω) is an unknown function, x, x ′ ∈ Ω ⊂ R n (Ω is a bounded domain). We are interested in Monte Carlo methods for evaluation of linear functionals of the solution of the following type:

J(u) = ∫ φ(x)u(x)dx = (φ, u), (2)
where φ(x) ∈ L 2 (Ω) is a given function. We can apply successive approximation method for solving integral equations:

u (i) = i ∑ j=0 K (j) f = f + Kf + . . . + K (i-1) f + K (i) f, i = 1, 2, . . . ( 3
)
where u (0) (x) ≡ f (x). It is known that the condition ∥K∥ L2 < 1 is a sufficient condition for convergence of the Liouville-Neumann series [START_REF] Curtiss | Monte Carlo methods for the iteration of linear operators[END_REF][START_REF] Dimov | Monte Carlo methods for applied scientists[END_REF]. Thus, when this condition is satisfied, the following statement holds:

u (i) -→ u as i → ∞.
Therefore,

J(u) = (φ, u) = lim i→∞ (φ, u (i) ) = lim i→∞   φ, i ∑ j=0 K (j) f   = lim i→∞ i ∑ j=0 ( φ, K (j) f
) .

Usually it is assumed that φ

(x) ∈ L 2 (Ω), because k(x, x ′ ) ∈ L 2 (Ω×Ω), f (x) ∈ L 2 (Ω). In a more general setting k(x, x ′ ) ∈ X(Ω × Ω), f (x) ∈ X(Ω),
where X is a Banach space. Then the given function φ(x) should belong to the adjoint Banach space X * , and the problem under consideration may by formulated in an alternative way:

v = K * v + φ, (4)
where v, φ ∈ X * (Ω), and

K * (Ω × Ω) ∈ [X * → X * ].
In such a way one may compute the value 2). An important case for practical computations is the case when X ≡ L 1 , where L 1 is defined in a standard way:

J(v) = ∫ f (x)v(x)dx = (f, v) instead of (
∥ f ∥ L1 = ∫ Ω |f (x)|dx; ∥ K ∥ L1 ≤ sup x ∫ Ω |k(x, x ′ )|dx ′ .
In this case the function φ(x) from the functional (2) belongs to

L ∞ , i.e. φ(x) ∈ L ∞ since L * 1 ≡ L ∞ .
It is also easy to see that (K * v, u) = (v, Ku), and also (φ, u) = (f, v). This fact is important for practical computations since often the computational complexity for solving the adjoint problem is smaller than the complexity for solving the original one. The above consideration shows that if φ(x) is a Dirac δ-function δ(x -x 0 ), then the functional

J(u) = (φ(x), u(x)) = (δ(x -x 0 ), u(x)) = u(x 0 ) is the solution of the integral equation at the point x 0 , if u ∈ L ∞ .
An approximation of the unknown value (φ, u) can be obtained using a truncated Liouville-Neumann series (3) for a sufficiently large i:

(φ, u (i) ) = (φ, f ) + (φ, Kf ) + . . . + (φ, K (i-1) f ) + (φ, K (i) f ).
So, we transform the problem for solving integral equations into a problem for approximate evaluation of a finite number of multidimensional integrals. We will use the following denotation (φ, K (j) f ) = I(j), where I(j) is a value, obtained after integration over Ω j+1 = Ω × . . . × Ω, j = 0, . . . , i. It is obvious that the calculation of the estimate (φ, u (i) ) can be replaced by evaluation of a sum of linear functionals of iterative functions of the following type (φ, K (j) f ), j = 0, . . . , i, which can be presented as:

(φ, K (j) f ) = ∫ Ω φ(t 0 )K (j) f (t 0 )dt 0 = = ∫ G φ(t 0 )k(t 0 , t 1 ) . . . k(t j-1 , t j )f (t j )dt 0 . . . dt j , ( 5
)
where t = (t 0 , . . . ,

t j ) ∈ G ≡ Ω j+1 ⊂ R n(j+1) . If we denote by F j (t) the integrand function F (t) = φ(t 0 )k(t 0 , t 1 ) . . . k(t j-1 , t j )f (t j ), t ∈ Ω j+1 ,
then we will obtain the following expression for (5):

I(j) = (φ, K (j) f ) = ∫ G F j (t)dt, t ∈ G ⊂ R n(j+1) . ( 6
)
Thus, we will consider the problem for approximate calculation of multiple integrals of the type [START_REF] Dimov | Monte Carlo method for numerical integration based on Sobol' sequences[END_REF].

The above consideration shows that there are two classes of possible stochastic approaches. The first one is the so-called biased approach when one is looking for a random variable which mathematical expectation is equal to the approximation of the solution problem by a truncated Liouville-Neumann series (3) for a sufficiently large i. An unbiased approach assumes that the formulated random variable is such that its mean value approaches the true solution of the problem. Obviously, in the first class of approaches there are two errors: a systematic one (a truncation error) R sys and a stochastic (a probabilistic) one, namely R N , which depends on the number N of values of the random variable,or the number of chains used in the estimate. In the case of unbiased stochastic methods one should analyse the only probabilistic error. In the case of biased stochastic methods more careful error analysis is needed: balancing of both systematic and stochastic error should be done in order to minimize the computational complexity of the methods (for more details, see [START_REF] Dimov | Monte Carlo methods for applied scientists[END_REF]).

Let us now present the simplest stochastic approach, known also as Crude Monte Carlo method (CMCM) for evaluating integrals. It is well-known that CMCM reduces problem to the approximate calculation of mathematical expectation which coincides with the unknown functional defined by (2) (see, for example, [START_REF] Dimov | Monte Carlo methods for applied scientists[END_REF]). The Monte Carlo quadrature formula is based on the probabilistic interpretation of an integral. If {x l }, l = 1, . . . , N is a uniformly distributed sequence in G, then the Crude Monte Carlo approximation to the integral I(j) is

I(j) ≈ I N (j) = V(G) N N ∑ l=1 F j (x l )
with integration error ε N (j) = |I(j) -

I N (j)| ≈ √ V ar(F ) N
, where V(G) is the volume of G and V ar(F j ) is the variance of the corresponding random variable whose mathematical expectation coincides with the unknown functional. This random variable depends on the integrand F (t).

Biased Stochastic Approaches

In this section, a brief description of the biased stochastic approaches for solving the problem (1)-( 2) presented in Section 2 is given.

Monte Carlo Method for Integral Equations

It is known (see [START_REF] Sobol | Monte Carlo Numerical Methods[END_REF][START_REF] Dimov | Monte Carlo methods for applied scientists[END_REF]) that the approximate calculation of linear functionals of the solution of an integral equation (φ, u) brings to the calculation of a finite sum of linear functionals of iterative functions (φ, K j f ), j = 0, . . . , i. First, we construct a random trajectory (Markov chain) T i of length i starting from state x 0 in the domain Ω:

T i : x 0 -→ x 1 -→ . . . -→ x i
according to the initial π(x) and transition p(x, x ′ ) probabilities. The functions π(x) and p(x, x ′ ) satisfy the requirements for non-negativeness, to be acceptable to function φ(x) and the kernel k(x, x ′ ) respectively and

∫ Ω π(x)dx = 1, ∫ Ω p(x, x ′ )dx ′ = 1 for any x ∈ Ω ⊂ R n .
From the above supposition on the kernel k(x, x ′ ) and the well-known facts that

Eθ i [φ] = (φ, u (i) ), where θ i [φ] = φ(x 0 ) π(x 0 ) i ∑ j=0 W j f (x j ),
and

W 0 = 1, W j = W j-1 k(x j-1 , x j ) p(x j-1 , x j ) , j = 1, . . . , i
it follows that the corresponding Monte Carlo estimation of (φ, u (i) ) can be presented in the following form:

(φ, u (i) ) ≈ 1 N N ∑ l=1 θ i [φ] l .
Therefore, the random variable θ i [φ] can be considered as an biased estimate of the desired value (φ, u) for i sufficiently large with a statistical error of order O(N -1/2 ), where N is the number of chains and θ i [φ] l is the value of θ i [φ] taken over the l-th chain.

It is very important that the same trajectories of the type T i can be used for a biased approximate evaluation of (φ, u (i) ) for various functions φ(x). Furthermore, they can be used for various integral equations with the same kernel k(x, x ′ ), but with different right-hand sides f (x).

Monte Carlo Algorithms based on Modified Sobol ΛΠ τ Sequences

ΛΠ τ sequences are uniformly distributed sequences (u.d.s.). The term u.d.s. was introduced by Hermann Weyl [START_REF] Weyl | Ueber die Gleichverteilung von Zahlen mod Eins[END_REF]. For practical purposes an u.d.s. should satisfy the following requirements [START_REF] Sobol | Monte Carlo Numerical Methods[END_REF][START_REF] Sobol | On quadratic formulas for functions of several variables satisfying a general Lipschitz condition[END_REF]: (i) the best asymptote as N → ∞, (ii) well distributed points for small N , and (iii) a computationally inexpensive algorithm. Suitable distributions such as ΛΠ τ sequences are also called (t, m, s)-nets and (t, s)-sequences in base b ≥ 2 [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF]. I. M. Sobol [START_REF] Sobol | Monte Carlo Numerical Methods[END_REF] defines his Π τ -meshes and ΛΠ τ sequences, which are (t, m, s)-nets and (t, s)-sequences in base 2, respectively. Subroutines to compute these points can be found in [START_REF] Bradley | Algorithm 659: Implementing Sobol's quasi random sequence generator[END_REF], and with more details in [START_REF] Sobol | Construction and comparison of high-dimensional Sobol' generators[END_REF]. Here we consider two randomized quasi-Monte Carlo algorithms based on Sobol ΛΠ τ sequences, namely MCA-MSS-1 and MCA-MSS-2.

Summary of the MCM-MSS-1.

One of the algorithms (based on a procedure of shaking) was proposed recently in [START_REF] Dimov | Monte Carlo method for numerical integration based on Sobol' sequences[END_REF]. The idea is that we take a Sobol ΛΠ τ d-dimensional point (vector) x. Then x is considered as a centrum of a sphere with a radius ρ << 1. A random point ξ ∈ U d uniformly distributed on the sphere is taken. Consider a random variable θ defined as a value of the integrand at that random point, i.e., θ j = F j (ξ). Consider random points

ξ (l) (ρ) ∈ U d , l = 1, . . . , N . Assume ξ (l) (ρ) = x (l) + ρω (l)
, where ω (l) is a unique

uniformly distributed vector in U d . The radius ρ is relatively small ρ ≪ 1 2 dj , such that ξ (l) (ρ) is still in the same elementary l th interval E d l = d ∏ ν=1 [ a (l) j 2 dν , a (l) ν + 1 2 dν ] ,
where the pattern ΛΠ τ point x (l) is. We use a subscript i in

E d i to indicate that the l-th ΛΠ τ point x (i) is in it. So, we assume that if x (l) ∈ E d l , then ξ (l) (ρ) ∈ E d l too.
It was proven in [START_REF] Dimov | Monte Carlo method for numerical integration based on Sobol' sequences[END_REF] that the mathematical expectation of the random variable θ j = F j (ξ) is equal to the value of the integral [START_REF] Dimov | Monte Carlo method for numerical integration based on Sobol' sequences[END_REF], that is

Eθ j = ∫ U d F j (t)
dt. This result allows to define a randomized algorithm called Monte Carlo Algorithm based on Modified Sobol Sequences (MCA-MSS-1).

In [START_REF] Dimov | Advanced algorithms for multidimensional sensitivity studies of large-scale air pollution models based on Sobol sequences, Computers and Mathematics with Applications 65[END_REF] the probability error of the algorithm MCA-MSS-1 is analysed. It is proved that for integrands with continuous and bounded first derivatives, i.e.

F j ∈ W 1 (L; U d ), where L = ∥F j ∥, it holds err(F j , d) ≤ c ′ d ∥F j ∥ N -1 2 -1 d and r(F j , d) ≤ c ′′ d ∥F j ∥ N -1 2 -1 d ,

Summary of the MCM-MSS-2.

The second algorithm called MCA-MSS-2 is a modification of MCA-MSS-1 algorithm. It is proposed and analysed in [START_REF] Dimov | Multidimensional Sensitivity Analysis of Large-scale Mathematical Models[END_REF].

It is assumed that n = m d , m ≥ 1. The unit cube U d is divided into m d disjoint sub
-domains, such that they coincide with the elementary d-dimensional

subintervals U d = ∪ m d l=1 K l , where K l = ∏ d ν=1 [a (l) ν , b (l) ν ), with b (l) ν -a (l) ν =
1 m for all l = 1, . . . , d. In such a way in each d-dimensional sub-domain K l there is exactly one ΛΠ τ point x (l) . Assuming that after shaking, the random point stays inside K l , i.e., ξ (l) (ρ) = x (l) + ρω (l) ∈ K l one may try to exploit the smoothness of the integrand in case if the integrand F j belongs to W 2 (L; U d ). In each subdomain K j the central point is denoted by s (j) , where s (j) = (s

(j) 1 , s (j) 2 , . . . , s (j) d ).
Suppose two random points ξ (l) and ξ (l) ′ are chosen, such that ξ (l) is selected during our procedure used in MCA-MSS-1. The second point ξ (l) ′ is chosen to be symmetric to ξ (l) according to the central point s (l) in each cube K l . In such a way the number of random points is 2m d . One may calculate all function values f (ξ (l) ) and f (ξ (l) ′ ), for l = 1, . . . , m d and approximate the value of the integral in the following way:

I(f ) ≈ 1 2m d 2N ∑ l=1 [ f (ξ (l) ) + f (ξ (l) ′ ) ] . ( 7
)
This estimate corresponds to MCA-MSS-2.

In [START_REF] Dimov | Multidimensional Sensitivity Analysis of Large-scale Mathematical Models[END_REF] it is proved that the Monte Carlo algorithm MCA-MSS-2 has an optimal rate of convergence for functions with continuous and bounded second derivative [START_REF] Dimov | Monte Carlo methods for applied scientists[END_REF]. This means that if

F j ∈ W 2 (L; U d ), where L = ∥F j ∥, it holds err(F j , d) ≤ c ′ d ∥F j ∥ N -1 2 -2 d and r(F j , d) ≤ c ′′ d ∥F j ∥ N -1 2 -2 d ,

Remark 31 Note that both MCA-MSS-1 and MCA-MSS-2 have one control parameter, that is the radius ρ of the sphere of shaking. At the same time, to be able to efficiently use this control parameter one should increase the computational complexity. It happens because after shaking the random point may leave the multidimensional sub-domain, and one needs to check if the random point is still in the same sub-domain. It is clear that the procedure of checking if a random point is inside the given domain is a computationally expensive procedure

when one has to deal with a large number of points.

Stratified Symmetrised Monte Carlo Approach

The main idea of the stratified symmetrised Monte Carlo approach is to generate a random point ξ (l) ∈ K l uniformly distributed inside K l and after that take the symmetric point ξ (l) ′ according to the central point s (l) . Such a completely randomized approach simulates the concept of MCA-MSS-2 algorithm, but the shaking is with different radiuses ρ in each sub-domain. This algorithm is called MCA-MSS-2-S, because this approach looks like the stratified symmetrised Monte Carlo. Obviously, MCA-MSS-2-S is less expensive than MCA-MSS-2, but the drawback is that there is not such a control parameter like the radius ρ, which can be considered as a parameter randomly chosen in each sub-domain K j .

It is important to notice that all three above mentioned algorithms MCA-MSS-1, MCA-MSS-2 and MCA-MSS-2-S have optimal (unimprovable) rate of convergence for the corresponding functional classes, that is MCA-MSS-1 is optimal in W 1 (L; U d ) and both MCA-MSS-2 and MCA-MSS-2-S are optimal in W 2 (L; U d ).

Specific Modification of the Algorithms for Solving Integral Equations

In the current versions of MCA-MSS-1 and MCA-MSS-2 algorithms the following modifications of the algorithms have been made:

-When the function from the desired linear functional is the δ-function, the dimension of the corresponding quasirandom and pseudorandom sequences coincides with the number of transitions in Liouville-Neumann series. For any other function, the dimension is larger than the number of transitions by one. -Minimal distance between quasirandom points is calculated to compute each of the integrals with a different dimension. This procedure during the preprocessing approximation of the integrals increases the total computational time for MCA-MSS-1 and MCA-MSS-2 algorithms. -We do one additional check: it is verified not only if the generated pseudorandom point lies inside the corresponding integration domain, but also if it lies inside the corresponding elementary multidimensional interval. Usually some points are rejected and that is why the total number of samples corresponding to MCA-MSS-2 is smaller than the number of samples for MCA-MSS-1 for some prescribed fixed values of coefficient ratio. -One important feature of MCA-MSS-2 is that it is efficient when the number of divisions is the same for all directions, i.e. the integration domain is divided into multidimensional subcubes.

Error Analysis of Biased Stochastic Approaches

So far we have considered only biased stochastic approaches. We analyse the problem for approximate calculation of linear functional of the solution of integral equation (φ, u). Let us denote the i-th iterative approximation of u (i ≥ 0) by u (i) and the Monte Carlo approximation by ũ. It is known that two errors systematic one (a truncation error) R sys and stochastic (a probabilistic) one, R N appear. Actually, we approximate the truncated Liouville-Neumann series.

If ε is a given sufficiently small positive parameter then

|(φ, u) -(φ, ũ)| ≤ |(φ, u) -(φ, u (i) )| + |(φ, u (i) ) -(φ, ũ)| = ε i + ε N < ε,
where ε i is the truncation error, ε N is the probability error. We can obtain a lower bound for the number i of iterations using [START_REF] Curtiss | Monte Carlo methods for the iteration of linear operators[END_REF]:

i > ln ε 2|(φ, u (0) ) -(φ, u)| ln ∥K∥ L2 ,
where u (0) is the initial approximation of u. Using the Cauchy-Bunyakowski inequality it is easy to show that

|(φ, u (0) ) -(φ, u)| ≤ ∥φ∥ L2(Ω) ∥u (0) -u∥ L2(Ω) .
The second multiplier can be estimated:

∥u (0) -u∥ L2(Ω) ≤ (∥I∥ L2 + ∥K∥ L2 + . . . + ∥K∥ i L2 + . . .)∥r (0) ∥ L2(Ω) ,
where r (j) = f -u (j) -Ku (j) , j = 0, 1, . . . .

Taking the limit for i → ∞ one can obtain

∥u (0) -u∥ L2 ≤ ∥r (0) ∥ L2(Ω) 1 -∥K∥ L2 .
It holds if the condition ∥K∥ L2 < 1 is fulfilled.

Unbiased Stochastic Approach

In this section an Unbiased Stochastic Approach (USA) will be discussed and analysed. The unbiased approach is based on the presentation of the functional under consideration as an infinite series. We are following a similar strategy like in our previous work [START_REF] Dimov | A New Walk on Equations Monte Carlo Method for Linear Algebraic Problems[END_REF] accept that now the state space is continuous and the transition density function, as well as the probability of absorbtion is no more discrete. In this case the random trajectory (Markov chain) T i introduced in Section 3, subsection 3.1 is infinite. The USA method is a generalization to integral equations of the method developed in [START_REF] Dimov | A New Walk on Equations Monte Carlo Method for Linear Algebraic Problems[END_REF] for solving systems of linear algebraic equations. The discrete Markov chain used in [START_REF] Dimov | A New Walk on Equations Monte Carlo Method for Linear Algebraic Problems[END_REF] is replaced by a continuous one with transition probabilities depending on the kernel k(x, x ′ ). In subsection 4.1 we describe the case 0 ≤ k(x, x ′ ) ≤ 1, where 1 -k(x, x ′ ) is an absorbtion rate.

In subsection 4.2 we develop the USA algorithm for more general kernels, where the above condition may not be fulfilled.

Unbiased approach for a simplified problem

We start with a simple case assuming that 0

≤ k(x, x ′ ) ≤ 1 for any x, x ′ ∈ Ω ≡ [0, 1] d , x = (x 1 , . . . , x d ), x ′ = (x ′ 1 , . . . , x ′ d ). Assume that Y ≡ U [0, 1] d is a uniformly distributed random variable in [0, 1] d .
Then we use the following stochastic presentation:

u(x) = E{k(x, Y )u(Y ) + (1 -k(x, Y ))0} + f (x),
which suggests a representation for the function u(x) based on the same ideas used in [START_REF] Dimov | A New Walk on Equations Monte Carlo Method for Linear Algebraic Problems[END_REF] for solving linear systems of equations. The probability for absorbtion is 1 -k(x, Y ), the score is the function f (x), and if the trajectory is not absorbed then it continues at point Y . We consider a random trajectory (Markov chain) T i of length i starting from state X = x = x 0 in the domain Ω assuming that

∂ ≡ R d \ Ω, i.e.: T i : x 0 -→ x 1 -→ . . . -→ x i , where P (x i+1 ∈ [a, b]/x i ̸ = ∂) = k(x, U i )P (U i ∈ [a, b]), P (x i+1 = ∂/x i ̸ = ∂) = 1 -k(x, U i ) and P (x i+1 = ∂/x i = ∂) = 1, and f(∂) = 0.
Then we have:

u(x) = E ( ∞ ∑ i=0 f (x i )/x 0 = x ) = E ( τ ∑ i=0 f (x i )/x 0 = x )
, where τ = inf i≥0

(x i ) = ∂.
Obviously, one can write

u(x) = E ( ∞ ∑ i=0 f (x i+1 )/x 0 = x ) + f (x) = E ( E ( ∞ ∑ i=0 f (x i+1 )/x 1
))

+ f (x) = ∫ Ω k(x, y)u(y)dy + f (x) (8)
The algorithm in this simplified case may be described as follows.

Algorithm 41 :

Unbiased stochastic algorithm for 0 ≤ k(x, y) ≤ 1, x, y ∈ [0, 1] d 1.
Initialization Input initial data: the kernel k(x, y), the function f (x), and the number of random trajectories N . 2. Calculations:

2.1. Set score=0 2.2. Set x 0 = x, test = 1 Do j = 1, N Do While (test ̸ = 0) score = score + f (x 0 ) U = rand[0, 1] d , V = rand[0, 1] If k(x 0 , U ) < V then test = 0 else x 0 = U Endif Endwhile Enddo 3. Compute the solution: u(x) = score N

Unbiased approach for the general problem

Now we may describe the method for the case when the kernel k(x, y) may be negative and also we allow that |k(x,

y)| ≥ 1, but still k(x, y) ∈ L 2 (Ω × Ω). If -1 ≤ k(x, Y ) ≤ 0, then k(x, Y
) is not any more the probability to stop. Instead, |k(x, Y )| is used as a probability to stop. But now one needs to compute -u(Y ) instead of u(Y ). In this case, the score is -f (x) instead of f (x). If |k(x, Y )| ≥ 1, then the random walk continues and the score is k(x, Y )f (x). The above described procedure leads to the following method, which is a modification of the previous approach. Assume we deal with the Markov chain T i of length i starting from state X = x = x 0 in the domain Ω assuming that δ ≡ R d -Ω:

T i : x 0 -→ x 1 -→ . . . -→ x i ,
where in this case

P (x i+1 ∈ b]/x i ̸ = ∂) = min (|k(x i , x i+1 )|, 1) , P (x i+1 = ∂/x i ̸ = ∂) = 1 -min (|k(x i , x i+1 )|, 1) ,
and

P (x i+1 = ∂/x i = ∂) = 1.
Then we have:

u(x) = f (x) + E ( ∞ ∑ k=1 k ∏ i=1 k(x i-1 , x i )f (x i-1 )/x 0 = x ) ,
where k(x, x ′ ) is defined in the following way:

k(x, x ′ ) =    k(x, x ′ ) if |k(x, x ′ )| ≥ 1 1 if 0 ≤ k(x, x ′ ) ≤ 1 -1 if -1 ≤ k(x, x ′ ) ≤ 0
Now, we may present the algorithm in this more general case:

Algorithm 42 :

Unbiased stochastic algorithm for the general case:

1. Initialization Input initial data: the kernel k(x, y), the function f (x), and the number of random trajectories N . 2. Calculations:

2.1. Set score=0 2.2. Do j = 1, N 2.3. Set x 0 = x, test = 1, prod=1 Do While (test ̸ = 0) score = score + f (x 0 ) * prod U = rand(0, 1) d If |k(x 0 , U )| > 1 then prod = prod * k(x 0 , U ) else x 0 = U V = rand[0, 1] d If |k(x 0 , U )| < V then test = 0 else prod = prod * sign (k(x 0 , U )), x 0 = U Endif Endwhile Enddo

Compute the solution:

u(x) = score N 4.3 Some extensions of the unbiased approach 4.3.1 Weak approximation In many situations, like for instance in sensitivity analysis [START_REF] Dimov | Multidimensional Sensitivity Analysis of Large-scale Mathematical Models[END_REF], one is not necessary interested in pointwise value of the solution but in mean values ∫ Ω φ(x)u(x)dx. In some applications the given function φ(x) may be considered as a weight function. Often the support of φ(x) is a relatively small subdomain of Ω.

We can still use a double randomization technique for computing such a quantity. If φ(x) is a probability density function defined on Ω, then the integral under consideration can be presented as a mean value of u(x), i.e.,

∫ Ω u(x)φ(x)dx = E (u(Y ))
, where Y is a random variable which admits φ(x) as a density function. This will only lead to a slight modification of the algorithm. More precisely, the starting point x 0 is the random variable Y . It is not always easy to sample from φ, more importantly, φ is not always a density function. There is another possibility that can be used if φ is bounded. One can pick the starting point x 0 uniformly distributed at random and multiply the sum of the values of the random variable (the score) by φ(x 0 ).

Global approximation

Another important extension is the possibility to compute the solution of the Fredholm integral equation of second kind in all points of the domain Ω using the results of the same unbiased USA algorithm. One should take into account that during the random walk many points of the domain are visited. The information obtained during this process can be used (after a proper statistical treatment) to obtain an approximation of the solution at all point of the domain. The following procedure may be used to do that:

(i) The starting point x 0 is picked up at random uniformly.

(ii) Because of that, all the visited points among all the random trajectories are also uniformly distributed. One may state that for each visited point we obtain a score, which mean value is an unbiased estimate of the solution at this point.

(iii) We may express the approximation of the solution at any point x by using the Taylor expansion. We assume that we have values of the solution in a neighborhood of x of distance h at points x i , i = 1, . . . , N h and h ∈ R d is a d-dimensional vector, such that h = x i -x:

u(x i ) = u(x + h i ) + ∇u(x) h i + 1 2 (h) T [D 2 u(x)]h + O(h 3 ), where ∇u(x) = [ ∂u(x) ∂x 1 , . . . , ∂f (x) ∂x d ] is the gradient, [D 2 u(x)] = [ ∂ 2 u(x) ∂x i ∂x k ] d i,k=1
is the Hessian matrix, and (h) T means transposed vector h. Obviously, the mean value of the solution at points x i may be expressed in the following way:

1 N h N h ∑ i=1 u(x i ) = u(x) + 1 N h N h ∑ i=1 h i ∇u(x) + 1 N h N h ∑ i=1 h 2 i 2 [D 2 u(x)] + O(h 3 ), ( 9
)
where the Hessian matrix may be estimated from above by the maximum of all second derivatives at point x.

Obviously, it is important to chose the value of h, the size of the d-dimensional domain A N h,x around the point x in which we want to compute the solution. The value h should be small enough, such that the second order term in the Taylor expansion is relatively small. At the same time, the number of points lying in A N h,x should not be too small, since than the variance of the first order term is too big. The mathematical expectation of the first order term in ( 9) is zero, and its variance depends on the number of point as O( 1N h ). We can now estimate the value u(x) by the mean value 1

N h ∑ N h i=1 u(x i ).

Numerical Results and Discussion

Numerical results are presented and discussed within two subsections. In the first one we are focussing on the comparison of the three classes of methods described above, namely Markov chain MC for integral equations, computing finite number of integrals, and the newly proposed unbiased stochastic algorithm. In the second subsection we do further investigation of the USA algorithm.

Comparison of the different algorithms

The numerical algorithms are tested on the following example [START_REF] Farnoosh | Monte Carlo Method for Solving Fredholm Integral Equations of the Second Kind[END_REF]:

u(x) = ∫ Ω k(x, x ′ )u(x ′ )dx ′ + f (x), Ω ≡ [0; 1] ( 10
)
where k(x, x ′ ) = x 2 e x ′ (x-1) , ( 11)

f (x) = x + (1 -x)e x , ( 12
)
φ(x) = δ(x -x 0 ). ( 13
)
The solution of this test problem is u(x) = e x . We are interested in an approximate calculation of (φ, u), where φ(x) = δ(x -x 0 ), x 0 = 0.5.

Numerical experiments

We have performed the following biased stochastic algorithms: MC algorithm and qMC based on Sobol sequences for integral equations and several algorithms for integrals into which the problem of solving integral equations is transformed. For the above set of algorithms we use Crude MCA, quasi-Monte Carlo algorithm based on ΛΠ τ Sobol quasirandom sequences, MCA-MSS-1, MCA-MSS-2, and MCA-MSS-2-S.

The algorithms have been studied after 10 runs to average the final approximation and for various sample sizes chosen according to proper sample size for Sobol quasirandom sequences. The samples number is defined by the requirement for error balancing following the approach given in Section 3, subsection 3.5. The number of iterations i/d is fixed, but it is chosen according to the L 2 -norm of the kernel k(x, x ′ ). For approximate computation of any integral I(j), j = 0, . . . , i different number of samples are chosen to satisfy the error balancing requirements. Because of the bias of the first two classes of algorithms there is a reason to present the relative errors, as well as, the approximation of the corresponding Liouville-Neumann series with a fixed length in Table 2, Table 3, and Table 4, as well. We have used symbolic computations to determine the values of the systematic error R sys for number of iterations up to 3. In Table 6 we present the computed values for the systematic error R sys at point x 0 = 0.1 and x 0 = 0.5 just to have an idea about the magnitude of this kind of error. One can clearly see that the systematic error decreases when the number of iterations increase. At the same time, for such small number of iteration i the systematic error dominates in almost all cases. The reason for that is that the analysis of balancing of stochastic and systematic errors are done following the assumption that Crude MC method is used. Actually, the applied algorithms, especially randomized quasi-Monte Carlo algorithms MCA-MSS-1, MCA-MSS-2 based on shaking of Sobol points are much higher quality and their stochastic error R N is much smaller than the systematic error R sys . Some typical values of R N are, say 10 -6 for relatively small number of random trajectories. These algorithms are very efficient when the norm of the kernel is relatively small. Then one may increase (but not too much) the number of iterations i, which will allow for a small number of N (and low computational complexity) to have a relatively high accuracy.

One should take into account that the choice of sample size to compute integrals defined on subdomains (IDSs) for the algorithm MCA-MSS-2-S is a crucial question. The IDSs have different dimensionality and this fact affects the corresponding sample size. The following strategy is applied here: first, the number N i is chosen to be approximately equal to 2 k (to give a possibility for a fair comparison with results obtained using Sobol quasirandom sequences) and to be close to the i-th power of an integer. We need this because the integration domain is the i-th dimensional unit cube. Then, N 1 := N i , N l < N i , l = 2, 3, . . . , i -1, and N l is the l-th power of an integer that is the number of sunintervals on each direction. Finally, we have chosen to compute each IDS with the same number of samples. It defines serious restrictions about the sample size applied during the numerical experiments. Nevertheless, we succeeded to achieve the goal of studying properties of MCA-MSS-2-S for solving integral equations based on the computation of a finite number of integrals.

Discussion and applicability

The tendencies concerning the behavior of biased MC algorithms here are similar to the tendencies observed in [START_REF] Dimov | Multidimensional Sensitivity Analysis of Large-scale Mathematical Models[END_REF]. The following observations can be drawn based on the results of our computations:

-MCA-MSS-2 algorithm gives a slightly larger relative stochastic error R N in comparison with MCA-MSS-1 related to the linear functional under consideration, and smaller relative stochastic error (almost 10 times smaller for the most cases) related to the finite number of integrals. -MCA-MSS-2-S algorithm gives the smallest relative stochastic error and the shortest computational time for a fixed number of steps in Liouville-Neumann series in comparison with MCA-MSS-1 and MCA-MSS-2, but serious restrictions about the choice of sample size exist for its implementation. -The main disadvantage of MCA-MSS-1 and MCA-MSS-2 algorithms is the high computational complexity due to computing the minimal distance between the generated original Sobol sequences. -MCM-int. eq. algorithm gives the smallest relative errors and corresponding to the problem of approximation of a linear functional of the integral equation [START_REF] Georgieva | Importance Separation for Solving Integral Equations[END_REF] in comparison with other biased MC algorithms and similar computing time to the Crude MC and MCA-MSS-2-S algorithms. -MCM-int. eq. has smaller relative errors in comparison with unbiased USA algorithm especially for smaller sample sizes (512, 1024) or similar relative errors for larger ones. -USA algorithm is the fastest and the most reliable for larger sample sizes than all other algorithms under consideration. The reason for that could be the fact that the algorithm is unbiased and also the probability for absorbtion of the random trajectory is chosen to depend on the problem data, namely, k(x, x ′ ) and f (x).

-The average number of the steps of the unbiased USA algorithm is smaller than 2 due to the fast convergence because i) the initial point at each chain is 0.5, and ii) the most of the kernel values are ranged in (0, 1; 0, 3) in this case. Just to make a comparison: for MC algorithms 2, 3, 4 jumps are necessary to reach the prescribed accuracy (the L 2 norm of the kernel is approximately equal to 0.3, i.e. one can expect a fast convergent iterative process).

-USA algorithm behavior has been studied for other initial points, especially for x 0 = 0.99; the average number of iteration steps for USA algorithm increased and was approximately equal to 2 (more precisely, it is around 2.4) for all sample sizes considered.

Further investigation of the USA algorithm

There are two parts in this subsection. In the first one we study a new example, to show how to handle the different difficulties described in Subsection 4.3. In the second part we compute a global approximation of the solution of the first example [START_REF] Georgieva | Importance Separation for Solving Integral Equations[END_REF].

General USA algorithm

In order to validate the general version of the unbiased algorithm we consider a second numerical example, where the kernel can take negative values and also its absolute value can be greater than 1:

u(x) = ∫ 1 0 k(x, x ′ )u(x ′ )dx ′ + f (x), (14)
where

k(x, x ′ ) = β(2x ′ -1)x 2 e x ′ (x-1) , ( 15
) f (x) = e x + β(2e x -2 -xe x -x). ( 16
) φ(x) = δ(x -x 0 ), (17)
where x 0 take values 0.1, 0.5, and 0.99 in different numerical experiments. The solution of the above equation is still u(x) = e x . We use special parameter β to be able to deal with various difficulties, namely, with negative kernels and values of the kernel greater then one in absolute value. If the values of the parameter β increases, then we approach instability barrier when the variance is getting infinite. Following the idea in subsection 4.3, Section 4 we also compute two quantities:

I 1 = ∫ 1 0 u(x)dx, (18)
as well as

I 2 = ∫ 1 0 2xu(x)dx (19)
to illustrate our new unbiased approach.

For the first quantity the initial value is just picked up uniformly, while for the second one we have two possibilities:

(i) pick the starting point x 0 according to the density 2x; this leads to a random variable distributed as √ U (U is a uniformly distributed random variable) using the well-known inversion method;

(ii) one can also start uniformly at random from the point x 0 and multiply the sum of values of the random variable (score) by 2x 0 . Some results for the numerical solution of the second example ( 14)-( 17) are presented in Table 7. All simulations are done for a fixed number of random trajectories N = 10 6 . One can see, that for each point x 0 = 0.1, 0.5, 0.99 the standard deviation and the error increases with β. This result is in full agreement with the theory, since when the parameter β increases the norm of the kernel k(x, x ′ ) is getting closer and closer to 1. That is why for high values of β we can observe some numerical instability.

In the second part of Table 7 the results obtained for another functional of the solution, namely I 2 (see, ( 19)) by applying different techniques of sampling: the first example is when we sample uniformly Ū , the next example is when we pick the starting point x 0 according to the density 2x (see, approach (i) described above), and the last example is when we start uniformly at random from the point x 0 and multiply the sum of values of the random variable (score) by 2x 0 (see, approach (ii) described above). One can see from the results presented in the table that approach (ii) has outperformed approach (i). These results give a sort of hint how to chose the sampling technique. The results also show that with the increasing of β when we approach to the limit of the convergency the relative error and standard deviation are getting bigger and bigger.

Global approximation

We are now dealing with the algorithm described in Subsection 4.3 of Section 4 to compute an approximation in all the points of the domain. It is also important to see how the accuracy improves when the number of random trajectories N increases if we apply our approach for computing the solution in all the points of the domain. Results for the relative error, standard deviation σ and mean value of number of steps i in each random trajectory are presented in Table 8.

The numerical results show that the relative error decreases with the increasing the number of random trajectories N according to the expectations. When the unbiased USA algorithm is used for computing the values of the solution in all the points of the domain, the following observations can be drown:

when the solution is computed at any point x by using the Taylor expansion the accuracy is satisfactory; the relative error increases as expected, but it is still small enough; if the parameter h is too big (h = 0.2), then the procedure smoothes the information obtained from the points where we have an unbiased estimation for the solution; that is why the increase of the number of random trajectories does not improve the accuracy; it means that the information which is used is too integrated ; for relatively small values of h, say h = 0.05, when N is small we have very small number of points in which we have unbiased estimates of the solution; the information is fragmentized, and that is why the relative error is big; with increasing the number of random trajectories, the information used is getting more integrated, and the error decreases significantly.

Conclusions and Future Plans

A new unbiased stochastic method for solving Fredholm integral equations of second kind was proposed and analysed. We have studied three possible approaches to compute linear functionals of the solution of integral equation under consideration:

biased Monte Carlo method based on evaluation of truncated Liouville-Neumann series, transformation of this problem into the problem of computing a finite number of integrals, and unbiased stochastic approach.

Five Monte Carlo algorithms for numerical integration have been applied in the second case, namely ) .

We have shown that the procedure of balancing of both stochastic R N and systematic R sys errors is very important for the quality of the biased algorithms.

In almost all numerical experiments performed, the unbiased USA algorithm performs with relatively high accuracy and low computational complexity in comparison with the best available biased algorithms. The numerical examples studied in this work are one-dimensional since we are comparing Monte Carlo algorithms only. It is obvious that the same algorithms can be used in any fixed dimension. However, the algorithms based on evaluation of finite number of integrals will suffer more from the effect of high dimensionality, because they are based on quadrature points.

As a further step, we would like to focus on the analysis of the variance of the unbiased USA algorithm, especially to study how the unbalancing of the kernel may reflect to the accuracy of the results. It is also important to study the computational complexity of USA in terms of number of floating point operations, and how the complexity depends on the problem's parameters. Some comparison of the stochastic approaches under consideration with known deterministic algorithms could be also considered.

  where the constants c ′ d and c ′′ d do not depend on N .

  where the constants c ′ d and c ′′ d do not depend on N .

-

  Crude Monte Carlo method based on a high quality SIMD-oriented Mersenne Twister pseudorandom number generator, -Quasi-Monte Carlo based on ΛΠ τ Sobol quasirandom sequences, a stratified symmetrised MC algorithm MCA-MSS-2-S, a randomized Quasi-Monte Carlo based on a special procedure of shaking ΛΠ τ Sobol quasirandom points with a convergence rate O

Table 1 .

 1 Relative errors and computational time for computing u(x0) at point x0 = 0.5 using MCM for integral equations.

	i	ε	N	Crude MC	Sobol seq.
				Rel. err.	Time (s)	Rel. err.	Time (s)
	1	0.4	128	0.05960	< 0.0001	0.05916	< 0.0001
	2	0.14	1024	0.00566	< 0.0001	0.00670	< 0.0001
	2	0.08	4 096	0.00294	0.01	0.00712	< 0.0001
	3	0.02	32 768	0.00092	0.17	0.00204	0.01
	4	0.018	65 536	0.00076	0.49	0.00096	0.03

Table 2 .

 2 Relative errors and computational time for computing u(x0) as a finite number (i) of integrals from the Liouville-Neumann series at x0 = 0.5 using Crude MCM.

	i	N		Integral equation		FNI
			Sobol seq.	Crude MC	Sobol seq. Crude MC
			Rel. err. Time (s) Rel. err. Time (s)	Rel. RN	Rel. RN
	1	128	0.05158 < 0.0001	0.05197 < 0.0001	0.00024	0.00016
	2	1024	0.01439 < 0.0001	0.01405 < 0.0001	2e-05	0.00033
	2 4 096	0.01438 < 0.0001	0.01422 < 0.0001	5e-06	0.00015
	3 32 768	0.00399	0.01	0.00404	0.11	1e-06	4e-05
	4 65 536	0.00111	0.03	0.00108	0.63	1e-06	3e-05

Table 3 .

 3 Relative errors and computational time for computing u(x0) as a finite number of integrals from the Liouville-Neumann series at x0 = 0.5 using MCA-MSS-1 and MCA-MSS-2 algorithms for integrals.

	i N	ρ	MCA-MSS-1	MCA-MSS-2
			Int. eq.	FNI	Int. eq.	FNI
			Sobol CMCM Time Sobol CMCM	CMCM
			Rel. err. Rel. err. (s) Rel. RN Rel. RN Rel. RN Time (s) Rel. RN
	1 128 2e-03 0.0516 0.0515 <1e-04 0.0003 0.0004 0.0514 <1e-04 0.0005
	2 1017 2e-04 0.0144 0.0144 0.02	2e-06 2e-06 0.0144	0.02	2e-05
	2 4 081 6e-05 0.0144 0.0144 0.14	2e-06 3e-06 0.0144	0.16	7e-06
	3 32 749 8e-06 0.0040 0.0040 10.4	8e-07 8e-07 0.0040 10.73	1e-06
	4 65 521 4e-06 0.0011 0.0011 54.9	1e-07 1e-07 0.0011	55.4	2e-07

Table 4 .

 4 Relative errors and computational time for computing u(x0) as a finite number (i) of integrals from the Liouville-Neumann series at x0 = 0.5 using MCA-MSS-2-S algorithm for integrals.

	i	N	Integral equation	Finite number of int, FNI
			Crude MC	Crude MC
			Rel. err.	Time (s)	Rel. RN
	1	128	0.05182	< 0.0001	4e-09
	2	1 024	0.01437	< 0.0001	6e-08
	2	4 096	0.01437	0.01	2e-08
	3	46 656	0.00399	0.31	2e-08
	4	4 096	0.00111	0.04	1e-06
	4	531 441	0.00111	5.97	7e-09

Table 5 .

 5 Relative errors and computational time for computing u(x0) at different initial points using unbiased USA for integral equations.

	N		x0 = 0.5			x0 = 0.99
		i	Rel. err.	Time (s)	i	Rel. err.	Time (s)
	128	1.3	0.02618	0.001	2.38	0.00237	0.001
	512	1.28	0.01160	0.002	2.39	0.00714	0.002
	1 024	1.28	0.01020	0.001	2.38	0.00199	0.005
	2 048	1.27	0.00422	0.006	2.38	0.00213	0.005
	16 384	1.27	0.00142	0.024	2.38	9e -06	0.044
	65 536	1.27	0.00075	0.096	2.37	0.00041	0.142
	131 072	1.27	0.00037	0.154	2.37	0.00024	0.281
	262 144	1.27	0.00014	0.305	2.37	0.00026	0.558
	4 194 304	1.27	0.00003	4.888	2.37	0.00011	8.802

Table 6 .

 6 Relative systematic error Rsys computed for three number of iterations i =

	1, 2, 3 at two points x0
	x0	i = 1	i = 2	i = 3
	x = 0.1	2.28 e-03	0.635 e-03	0.176 e-03
	x = 0.5	51.8 e-03	14.3 e-03	3.99 e-03

Table 7 .

 7 Relative error (Err.), standard deviation (σ) and mean value of steps of absorbtion i for the USA approach for various values of β; N = 10 6 .

			β = 1			β = 2			β = 4
		Err.	σ	i	Err.	σ	i	Err.	σ	i
	x = 0.1	7e-05 0.105 1.003 1e-05 0.17 1.009 9e-04	0.62	1.02
	x = 0.5	2e-04 0.600 1.116 1e-04 0.97	1.27	0.004	8.43	1.63
	x = 0.99	5e-04 1.388 1.588 3e-04 2.76	2.05	0.008	32.9	2.48
	Ū meth.1	6e-04 0.900 1.178 8e-04 1.39	1.38	0.012 11.65 1.70
	I2 meth.(i)	1e-03	1.95	1.176 1e-03 2.63	1.38	0.015	20.5	1.70
	I2 meth.(ii) 2e-04	1.03	1.27	8e-04 1.70	1.58	0.013	14	2

Table 8 .

 8 Relative error, standard deviation and mean value of steps of absorbtion i for the USA approach for various values of number of random trajectories N .

		σ	i	N = 10 4	N = 10 5	N = 10 6	N = 10 7
	x = 0.1	0.14	1.009	3e-04	6e-04	5e-05	4e-05
	x = 0.5	0.76	1.27	3e-03	1e-03	7e-04	1e-05
	x = 0.99	0.92	2.37	2e-03	1e-03	8e-04	8e-05
	h = 0.2	0.91	1.41	5e-03	6e-03	7e-03	6e-03
	h = 0.1			8e-03	2e-03	2e-03	1.5e-03
	h = 0.05			1e-02	4e-04	8e-04	7e-05

Acknowledgment

This work has been supported by the EC FP7 Project AComIn (FP7-REGPOT-2012-2013-1), by ISITV, Université de Sud Toulon-Var, as well as by the Bulgarian NSF Grants DMU 03/61/2011 and DCVP 02/1/2010. We would like to thank Dr. Mehmet Ersoy who helped us with the symbolic computations.