A New Algorithm for System of Integral Equations

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.

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 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 Section 2, 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 where or in vector form here operator where is the space of bounded function almost everywhere and where the norm of is Suppose the spectral radius satisfy the inequalities where .

Let Markov chain with transition density be where is the probability of absorption at the point , where is the random number of the last moment and in initial moment .

A standard vector algorithm of Monte Carlo for is where is a unit matrix, is a kernel matrix , and is the transition density function at the points . The condition for unbiasedness is We will assume also that the spectral radius of the operator obtained from by the substitution is less than one. Then, by using standard methods of Monte Carlo theory we can show that where can be considered as matrix weight and The Monte Carlo method is used to estimate linear functionals of the form where with Let the point be distributed with initial probability density such that Then, obviously (see Mikhailov [17]), The random vector with weight is computed by the formula Precisely such a vector algorithm, corresponding to the representation , 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 point . We will define two types of Markov chain and by the following way.

(a) Definition of the First Homogeneous Markov Chain. Now we simulate the Markov chain with state. Initial state will simulate according to initial distribution and the next with the transition matrix Here and Let be a random absorption moment with , a life time of chain.

(b) A second homogeneous Markov chain with space phase is defined by the following way.

Firstly, we define the transition density matrix as Let an initial point ; using we will simulate initial moment , then according to the transition matrix we are able to simulate again the next state of chain . It means with the probability .

The next phase coordinates of the chain simulated according to . The probability of absorption of the trajectory is . Let be known then the next value of will be defined according to the matrix and next random point simulated according to the probability density function and so on.

Let be some random variable which is defined by the set of trajectory Markov chains. The mathematical expectations of random variable will be Let us consider calculation of the functional , where column vector. Let us compute the functional . For doing this task we introduce two well-known estimators according to the Monte Carlo theory. First of them is analog of absorption estimator and the second one is analog of collision estimator

Theorem 1. If then and if then In this case .

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. Let be bounded domain in with enough smooth boundary of and is the cylinder in with parallel spin axis . The basement is the domain on the surface and is the fixed constant. The functions where stands for a continuous function on the closed domain .

Now consider the following initial boundary value problem (BVP) for system of parabolic equations: where the coefficients , and with initial and boundary conditions Further suppose , and coefficients   are given such that there exists unique solution Ladyzhenskaya et al. [18] and Lions [19] of the initial BVR (24)-(25) and where is the set of continuous functions in the given region with continuous derivatives and .

Now we construct unbiased estimator for the problem (24)-(25) in the arbitrary point on the trajectory some random process. For that we use mean value formula and construct some special system of integral equations for in special constructed domains (spheroid or balloid with the center ).

According to Section 2 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 solution for heat equation is given by Firstly, we define a special domain using a fundamental solution of the heat equation which depends on and points as The domain , we call balloid and , spheroid with the center in the point . From the definition balloid , described by following inequality (Kupcov [20]): Each section with the sectional plain of balloid when will be -dimensional ball with the center and with the radius Let and where is the minimum distance from point until the boundary; that is, In this case 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 with its integral from the spheroid and balloid with the center in the point .

Lemma 2 (Kurbanmuradov [22]). Let the function satisfy the following equation: Then the following formula of mean  is true (mean value formula): where here ds is the element of small area of sphere . 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 domains , which depends on positive parameters and point ,where where defined analogous changing a for (see above Lemma 2). The domain we will call a balloid with radius which a center in a point and a boundary is spheroid. Here Let and where is the distance from the point to the boundary of domain. In this case Appling the expression (34) to each of the equations we will get the following system of integral equations : where 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: where is a random variable with density functions random point on the surface , which has a density function unit sphere, ds element of surface, square of the surface unit sphere, and Gamma function.

Let us consider the second terms of (38) where . Then the density of Beta distribution with parameters and the density of Gamma distribution with parameters ,—unit random vector, where is the random variable with density function and is another random variable with the density function .

Let then and the function is the transition density in with fixed point , where Let be a random point of balloid which has the following density function in the fixed point .

In this case 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: where is defined by (40) and is determined by (41).

The proof of Theorem 3 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 The functions are the transition density functions in at a fixed point . We will define in Ω a random process as was proposed in Section 2.

Let us define a transition matrix as where , and let Now we will define the density function of transition matrix : where Then we will fix the initial point and the number of equations . Let an initial moment at the point ; we will have one particle. For one step a particle moves from its position according to the transition matrix and moves with probability from the point to the point . The next point will be simulated using the density function .

The probability of breaking of trajectory in the point is The next coordinate of the particle will be defined in the following way.

(1) If the density function of the point equals in the fixed point then where the sequence of independent random variables with the density function and independent isotropic vectors. The value will be defined as (50).

(2) If the density function of the point is equal to at a fixed then where 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 , secondly simulate , uniformly distributed random variable on (0, 1), and thirdly simulate , Beta distributed random variable with the parameters ().
(b) If then we will go to (a) and so on; otherwise .
(3) If the density function at the point equals under fixed point , then where is sequence of independent Gamma distributed random variables with parameters , Beta distributed random variables with parameters , and independent isotropic vectors, respectively.

If at the moment was held break, then we will put . 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 Lemma 5.

Lemma 5. With the probability one Markov chain converges when to the random point of boundary , or it is absorbed inside of the domain.

Proof. Since is decreasing sequence and , it has a limit . Let be algebra, which was generated by random variables From the definition and (57)–(59) it follows that is measurable relatively . It is obvious that the coordinates of vector process formed limited martingale relatively : where are limited martingale; it converges with the probability one Shiryaev [24].
Let be the limit vector. We show that . If then the process is broken inside the domain. Let . As far as the process converges, according to the formulas (50)–(58) we have where is strictly positive. Applying Lebesque Theorem (about the limited convergence) we get It means . Then from the definition of and using the formulae (49) we obtain Lemma is proven.

4. Construction Unbiased and -Biased Estimators

Let be the trajectories of random process which was described above. We will define on it the sequence of the random variables . Let where is defined as follows: Here are the sequences of independent Gamma function with the parameters () and Beta function with parameters () distributed random variables, respectively. We will define the sequence If at the moment happen break , we will put where algebra generated until the moment .

Theorem 6. Let the sequence be form martingale with , respectively. If Then the sequence will be uniformly integrable martingale.

Proof. From the definition the is measurable. In this case As far as the sequence is martingale . We can show the uniformly integrability of . To do that it is enough to show .
Since and is bounded domain , for any . From the condition of theorem , it is followed . It means is uniformly integrable.
The theorem is proved.

Now we will construct computable (realizable) estimator . We will take -neighborhoods of the domain .

Let be a breaking moment of process inside of domain. is the first passage moment . , stopping moment of process . In this case the probability of absorbing will be From Lemma 5 it follows that . It could be proved that the mathematical expectation of the stopping time of Markov process is finite.

Theorem 7. Let the condition of Theorem 6 be satisfied; then the estimator will be unbiased estimator with finite variance, where is th component of solution vector .

Proof. Since is uniformly integrable martingale and is Markov moment, according to the theorem in ([24, 25]) for the martingale we obtain From the definition holds . From condition of Theorem 6   is valid accordingly the variance is finite. The theorem is proved.

Further, from we could construct biased but computable (realizable) estimator . Let, for , and for , closed to the point of boundary . will be obtained with changing Let us evaluate bias . It is clear that If , in this case the process is broken when do not reach the boundary and . If then .

Let be a module of continuity function . In this case it is true: where , since then . Finiteness of variance followed from . The proposed algorithm we could generalize for the case with variable coefficients and one could get the same results.

5. Computational Example

Let be bounded domains, . We will consider for some mode linitial boundary value problem for with the initial boundary conditions As domain is chosen as the simple ball, .

The coefficients The initial and boundary conditions Left hand sides The exact solutions are known:

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 coefficients . 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.

Conflict of Interests

The authors declare that they have no conflict of interests regarding publication of this paper and they do not have direct financial relation that might lead to conflict of interests for any of the authors.

Acknowledgments

The authors would like to thank the anonymous referees, especially the second referee, for their constructive comments and suggestions.