In engineering, social and other areas, a lot of problems can be converted to Volterra integral equations to solve, such as elastic system in aviation, viscoelastic and electromagnetic material system and biological system, and some differential equations are often transformed into integral equations to solve in order to simplify the calculation. For example, the drying process in airflow, pipe heating, gas absorption and some other physical processes can be reduced to the Goursat problem. Then, some of the Goursat problem can be described by Volterra integral equations . Another example, when one-dimensional situations are concerned and the coolant flow is incompressible, the definite solution problem of the transpiration cooling control with surface ablation appears as Volterra integral equations of second kind . In practice, the analytical solutions for this kind of integral equations are difficult to obtain. Therefore, it is more practical to research the numerical method for solving this kind of integral equations.
The main aim of this paper is to propose a numerical algorithm based on Monte Carlo method for approximating solutions of the following system of Volterra integral equations
where, are known kernel functions, the functions,
are given and defined in, and are the unknown functions to be determined. One of the earliest methods for solving integral equations using Monte Carlo method was proposed by Albert , and was later developed . Literatures - employed Monte Carlo method to solve numerical solutions of Fredholm integral equations of the second kind. But very few studies are devoted to employing Monte Carlo method to solve Volterra integral equations and the system of Volterra integral equations. In this paper, we present and discuss a relaxed Monte Carlo approach with importance sampling to solve numerically systems of Volterra integral equations. Due to less accuracy and lower efficiency of Monte Carlo method, in this paper, combination of Monte Carlo and quadrature formula will be used to deal with Equation (1) and importance sampling is applied to accelerate the convergence and improve the accuracy of Monte Carlo method. Some numerical examples are given to show the efficiency and the feasibility of proposed Monte Carlo method.
2. Discretizing System of Integral Equations
Here, Newton-Cotes quadrature formula is used to discretize Equation (1). Dividing the interval into subintervals with step length, defining, where . For convenience, denoting the notation, , , where,. As,. Thus the following linear algebra system can be obtain
, , , is the weight of Newton-
Cotes quadrature formula. The matrix of coefficients of Equation (2) is. If we assume that there exists a unique solution of (2), the solution would be a numerical approximation of (1). This process will produce an error which is determined by numerical quadrature formula and can be reduced by increasing the number of nodes for a given quadrature formula. For a large number of nodes, Equation (2) is too large to solve directly. It is well known that Monte Carlo technique has a unique advantage for large systems or high-dimensional problems. At the same time, this method can obtain function values at some specified points or their linear combination that is just what researchers need. But for determined numerical methods, in order to obtain function value at a certain point, it is often necessary that find function values for all nodes. Here relaxed Monte Carlo method is used to Equation (2) based on a random sample from Markov chain with discrete state. According to theory of importance sampling, probability transition kernel is selected to suggest a possible move. To obtain solution of the linear algebraic system (2), the following iterative formula is considered
where is the iterative matrix, is a diagonal matrix with elements, , relaxation parameter of the iterative formula is chosen such that it minimizes the norm of for accelerating the convergence, and. The iterative formula (3) can define a Neumann series, as following
Set iterative initial value, is the exact solution of Equation (2), the truncation error and convergency of the iterative formula (3) can be obtained by the following expression
This conclusion can be proved by using theories in numerical analysis. Here, the iterative matrix satisfies
To achieve a desirable norm in each row of, a set of relaxation parameters, will be used in
place of a single value. According to the arguments of Faddeev and Faddeeva  , the relaxed Monte Carlo method will converge if
here denotes the row norm of the given matrix.
3. Relaxed Monte Carlo Method with Importance Sampling
For Neumann series (4), we have
In order to obtain the approximation solution of linear system (2) and system of integral Equation (1), the kth iteration of will be evaluated by means of computing the following series
Construct the Markov chain
on the state space. Let the initial probability and the transition probability of Markov chain respectively
and they must satisfy and for any. According to non-after-effect
property of Markov chain, one can get
For and the mth series (5), let, estimators are established in the following form
The weight function of Markov chain is defined as follows
By expressions (7) and (8), the following conclusion can be gotten.
Theorem 3.1 For the given, we have
This theorem is easy to prove.
In the light of the expression (7), the following estimator is defined
Due to Theorem 3.1, the conclusion (11) is easy to prove.
To estimate, by the transition probability, random paths of Markov chain are simulated
the length of Markov chain is defined by, is the precision of truncation error and given in advance. Then one can evaluate the sample mean
If the standard deviation is bounded, according to the Central Limit Theorem, we would obtain
So the precision of the estimator in the sense of probability can be measured by its variance.
Based upon the minimum variance of estimator, by the variance expression
the transition probability of Markov chain should be chosen in the following form
This form of leads to that more samples are taken in regions which have higher function values. This is importance sampling.
According to the obtained approximation of global approximation functions of solutions and of Equation (1) would be achieved
4. Numerical Examples
In this section, we employ the proposed relaxed Monte Carlo method with importance sampling (say RMCIS) to compute the numerical solution of some examples and compare it with their exact solutions. The numerical results are presented in Table 1 and Table 2, where AE means absolute error for. We plot the
Table 1. Numerical results of Example 1 with.
Example 1 Consider the equations  
Table 2.Numerical results of Example 2 with.
Figure 1. The figure of average absolute errors (MAE) for Example 1 at eleven points, (a) for and (b) for.
Figure 2.The figure of average absolute errors (MAE) for Example 2 at eleven points, (a) for and (b) for.
where. The exact solutions are and
. The numerical results are listed in Table 1.
Example 2 Consider the equations 
where. The exact solutions are and, The nu-
merical results are listed in Table 2.
In this paper, a relaxed Monte Carlo numerical method is provided to solve a system of linear Volterra integral equations. The most important advantage of this method is simplicity and easy-to-apply in programming, in comparison with other methods. The implementation of current approach RMCIS is effective. The numerical examples that have been presented in the paper and the compared results support our claims.
This research was supported by National Natural Science Foundation of China under Grant No. 11361036, Specialized Research Fund for the Doctoral Program of Higher Education under Grant No.20131514110005, Natural Science Foundation of Inner Mongolia under Grant No. 2015MS0104.
 Yang, Z.J., Xu, Y.H. and Yang, X.S. (1993) The Second-Type Volterra Integral Equations of Transpiration Cooling System and Its Numerical Solution. Journal of University of Science and Technology of China, 23, 310-317.
 Albert, G.E. and Meyer, M.A. (1956) Symposium of Monte Carlo Methods. A general Theory of Stochastic Estimates of the Neumann Series for Solution of Certain Fredholm Integral Equations and Related Series. Wiley, New York.
 Farnoosh, R. and Ebrahimi, M. (2008) Monte Carlo Method for Solving Fredholm Integral Equations of the Second Kind. Applied Mathematics and Computation, 195, 309-315. http://dx.doi.org/10.1016/j.amc.2007.04.097
 Hong, Z.M., Yan, Z.Z. and Chen, J.R. (2012) Monte Carlo Method for Solving Fredholm Integral Equations of the Second Kind. Transport Theory and Statistical Physics, 41, 513-528. http://dx.doi.org/10.1080/00411450.2012.695317
 Wei, J. and Zhong, C. (2013) Solving a System of Linear Volterra Integral Equations Using the New Reproducing Kernel Method. Applied Mathematics and Computation, 219, 10225-10230. http://dx.doi.org/10.1016/j.amc.2013.03.123
 Rabbani, M., Maleknejad, K. and Aghazadeh, N. (2007) Numerical Computational Solution of the Volterra Integral Equations System of the Second Kind by Using an Expansion Method. Applied Mathematics and Computation, 187, 1143-1146. http://dx.doi.org/10.1016/j.amc.2006.09.012
 Babolian, E. and Mordad, M. (2011) A Numerical Method for Solving Systems of Linear and Nonlinear Integral Equations of the Second Kind by Hat Basis Functions. Computers and Mathematics with Applications, 62, 187-198. http://dx.doi.org/10.1016/j.camwa.2011.04.066