A New Monte Carlo Approach for the Solution of Semi-Linear Neumann Boundary Value Problem

Abstract

In this paper we will consider the Neumann boundary-value problem for the Laplace operator with a polynomial nonlinearity on the right-hand side. We utilize Greens formula and an elliptic mean-value theorem. This allows us to derive a probabilistic representation as special integral equation which equates the value of the function u(x) at the centre-point x, with its integral over the domain D, and on boundary of the domain G. According to this probabilistic representation, a Markov branching process is constructed and an unbiased estimator of the solution of the nonlinear problem is formed by taking the mathematical expectation over these branching processes. The unbiased estimator which we derive has a finite variance. The complexity of the proposed algorithms is investigated, and the ways decreasing of the computational work are proposed. Some particular cases are considered in detail. The effectiveness of two estimators “by absorption” and “by collision” is studied in the nonlinear case.

1. Introduction

In this work, we will propose new probabilostic algorithms for the solution of the Laplace operator with a polynomial nonlinearity on the right-hand side. The complexity of the proposed algorithms is investigated, and the ways decreasing of the computational work are proposed. For the solution of the semi-linear Neumann problem, a new and effective Monte Carlo algorithm is derived. Some particular cases are considered in detail. The effectiveness of two estimators “by absorption” and “by collision” is studied in the nonlinear case.

2. Solutions of semi-linear Neumann problems

Assume $D$ is a bounded convex domain in $R^n,n\ge 3$ with a regular smooth boundary $\mathrm{\Gamma }$, and $x=\left(x_1,\dots ,x_n\right)\in D$. Consider the following boundary value problem:

$$\begin{array}{l}\mathit{Lu}\left(x\right)\equiv -\mathrm{\Delta }u\left(x\right)+g\left(x\right)u\left(x\right)=f\left(x,u\left(x\right)\right),\\ {\left.\frac{\partial u\left(x\right)}{\partial n_x}\right|}_{\mathrm{\Gamma }}=0.\end{array}$$E1

where $n_x$ is the external normal on the boundary $\mathrm{\Gamma }$ at the point $x$, and the function $g\left(x\right)>0$ and $g\left(x\right)\in C\left(\overline{D}\right)$. We shall determine the class of functions

$$v=\left\{\omega \left(x\right):B^{\ast }M(x,0)+\psi \left(\mathit{xt}\right)\ge 0,\text{in}D\right\};$$E2

Here $M\left(x,y\right)$ is the Levy function for operator L (see: [1]), $\psi \left(x\right)$ is a continuous function, which has continuous derivatives of second order in $D\cup \mathrm{\Gamma }$, and $B$ is a positive constant. It is known, that this class of functions is not empty.

Suppose that the function $f\left(x,\omega \left(x\right)\right)$ satisfies the following three conditions:

1. $\partial f\left(x,\omega \right)/\partial \omega \ge 0$ is an increasing function at $\omega \ge 0$;

2. There exists such $k$, that $f\left(x,\omega \right)=o\left(r^{k-2}\right)$ and $\frac{\partial f\left(x,\omega \right)}{\partial \omega }=o\left(r^{k-1}\right)$, here $r$ is the distance from the origin to the point $x$;

3. The limit $\underset{x\to 0}{lim}\left[\omega \left(x\right)h^{-1}(x,y)\right]=B$ holds, where $y=\left(y_1,\dots ,y_n\right)$$\in D$ and

$$h\left(x,y\right)=\frac{1}{{\left(-1\right)}^{n/2}\left(n-2\right){\omega }_n}{\left[{\left(-1\right)}^n\sum _{i=1}^n{x}_i-y_i){}^2\right]}^{\left(2-n\right)/n},$$E3

and ${\omega }_n$ is a surface element of the unit sphere.

The next assertion is then true.

Lemma 1 The unique positive solution of(1)exists within the class of functions,$v$, if the above conditions a), b), and c) are satisfied.

This assertion follows from the results of [2], therefore the proof is omitted. Thus, there exists a unique function $\omega \left(x\right)=u\left(x\right)\ge 0$ (from lemma 1) that satisfies conditions (a), (b), and (c) and which is the solution of problem (1).

Let us construct a computational scheme for the solution of problem (1) using proposed algorithms from the work [3]. Suppose $\mathrm{\Gamma }$ is a sufficiently smooth boundary and the function $v\left(x\right)\in C^2\left(D\right)\cap C^1\left(D\right)$. Then, we can use Green’s identity to obtain

$$\underset{D}{\int }\left(v\left(y\right)\mathit{Lu}\left(y\right)-u\left(y\right)\mathit{Lv}\left(y\right)\right)\mathit{dy}=\underset{\mathrm{\Gamma }}{\int }u\left(y\right)\frac{\partial v\left(y\right)}{\partial n_y}-v\left(y\right)\frac{\partial u\left(y\right)}{\partial n_y}\mathit{dy}\mathrm{\Gamma }$$E4

Let us $c^2>\underset{x\in \overline{D}}{max}g\left(x\right)$. then using formula (4) we can derive the integral equation for $u\left(x\right)$

$$\begin{array}{ll}u\left(x\right)=& \frac{1}{{\omega }_n}\left(\underset{D}{\int }\frac{exp\left(-\mathit{c\rho }\right)}{{\rho }^{n-2}}f\left(y\right),u\left(y\right)\right)\mathit{dy}+\\ & +\underset{D}{\int }\frac{A_2\left(\rho \right)-\frac{g\left(y\right)}{n-2}\rho exp\left(-\mathit{c\rho }\right)}{{\rho }^{n-1}}u\left(y\right)\mathit{dy}+\underset{\mathrm{\Gamma }}{\int }{A}_1\left(\rho \right)\frac{cos{\phi }_{\mathit{xy}}}{{\rho }^{n-1}}{d}_y\mathrm{\Gamma }.\end{array}$$E5

Here ${\phi }_{\mathit{xy}}$ is the angle between the normal $n_y$ with the vector $y-x$, $\rho =\mid x-y\mid$,

$$A_1\left(\rho \right)=\left(1+\frac{\mathit{c\rho }}{n-2}\right)exp\left(-\mathit{c\rho }\right),A_2\left(\rho \right)=\frac{c^2\rho -c\left(n-3\right)}{n-2}exp\left(-\mathit{c\rho }\right).$$E6

It is easy to show that ${\left(A_1\left(\rho \right)\right)}^{\prime }=-A_2\left(\rho \right)$, $A_1\left(0\right)=1$. Thus,

$$\underset{0}{\overset{\rho }{\int }}{A}_2\left(\rho \right)\mathit{d\rho }=1-A_1\left(\rho \right)$$E7

Furthermore, if $x\in \mathrm{\Gamma }$, then ${\omega }_n$ should be changed into ${\omega }_n/2$.

The linear case, when $f\left(x,u\left(x\right)\right)=f\left(x\right)$ was considered by A.Sipin. In his work [4] sequences of unbiased estimators were developed and the question of simulating the required special densities was discussed. Furthermore, using the results of [4], we will present an algorithm to construct of unbiased estimators for the solution of problem (1).

For $u\left(x\right)$, after transforming Eq. (1) as discussed below, we obtain the following integral equation

$$\begin{array}{ll}u\left(x\right)=& \underset{D}{\int }\left(1-\pi \left(x,y\right)\right)\mathrm{\Pi }\left(x,y\right)u\left(y\right)\mathit{d\mu }\left(y\right)+\\ & ++\frac{1+{\chi }_{\mathrm{\Gamma }}\left(x\right)}{{\omega }_n}\underset{D}{\int }f\left(y,u\left(y\right)\right)\frac{exp\left(-\mathit{c\rho }\right)}{{\rho }^{n-2}}\mathit{d\lambda }\left(y\right).\end{array}$$E8

Here

$$\begin{array}{c}{\chi }_{\mathrm{\Gamma }}\left(x\right)=\left\{\begin{array}{cc}0,& x\notin \mathrm{\Gamma },\\ 1,& x\in \mathrm{\Gamma }.\end{array}\right.\pi \left(x,y\right)=\left\{\begin{array}{ll}0,& y\in \mathrm{\Gamma }\\ \frac{g\left(y\right)\rho exp\left(-\mathit{c\rho }\right)}{\left(n-2\right)A_2\left(\rho \right)},& y\in D\end{array}\right.\\ \mathrm{\Pi }\left(x,y\right)=\frac{1+{\chi }_{\mathrm{\Gamma }}\left(x\right)}{{\omega }_n}\left\{\begin{array}{ll}\left(A_1\left(\rho \right)cos{\phi }_{\mathit{xy}}\right)/{\rho }^{n-1},& y\in \mathrm{\Gamma },\\ A_2\left(\rho \right)/{\rho }^{n-2},& y\in D\end{array}\right.\end{array}$$E9

The measure, $\mu$, is determined by the $\sigma$-algebra of Borel subsets $\overline{D}$ with the equality $\mu \left(E\right)=\lambda \left(E\right)+S\left(E\cap \mathrm{\Gamma }\right)$, where Lebesque $\lambda$-measure is in $R^n$, $S$ is a surface element. It is easy to see that $\mathrm{\Pi }\left(x,y\right)\ge 0$. Because the domain $D$ is convex, equality (7) implies that $\mathrm{\Pi }\left(x,y\right)$ is a transition probability density. Assume $q\left(x\right)\not\equiv 0$,$c^2>\underset{x\in \overline{D}}{max}g\left(x\right)$. Then, $0\le \pi \left(x,y\right)\le 1$. Let now consider

$$f\left(x,u\left(x\right)\right)=\sum _{i=0}^m{b}_i\left(x\right)u^i\left(x\right),x\in \overline{D},$$E10

where $b_i\left(x\right)\ge 0$. We choose ${\alpha }_i,i=0,\dots ,m+1$, to satisfy the condition:

$$\sum _{i=0}^{m+1}{\alpha }_i=1,\sum _{i=1}^{m}i{\alpha }_i+{\alpha }_{m+1}\le 1$$E11

The iterative method for the solution of Eq. (8) will converge only with small coefficients, $b_i\left(x\right)$. Therefore, we assume that parameters ${\alpha }_i$ satisfy the following condition for each $i$, $0\le i\le m+1$:

$$\parallel b_i\left(x\right){\parallel }_{C\left(\overline{D}\right)}\le 2{\alpha }_i\left(n-1\right)/d^2,$$E12

here $d$ is diameter of the domain $D$, and ${\alpha }_i\not\equiv 0$, $i=0,\dots ,m$. Using (10) and (11) we can transform Eq. (8), as follows

$$\begin{array}{ll}u\left(x\right)=& {\alpha }_{m+1}\underset{\overline{D}}{\int }{k}_{m+1}\left(x,y\right)\mathrm{\Pi }\left(x,y\right)u\left(y\right)\mathit{d\mu }\left(y\right)+\\ & +\frac{1+{\chi }_{\mathrm{\Gamma }}\left(x\right)}{{\omega }_n}\sum _{i=1}^m{\alpha }_i\underset{D}{\int }\frac{K_i\left(x,y\right)}{{\rho }^{n-2}}{u}^i\left(y\right)\mathit{d\lambda }\left(y\right)+\\ & +\frac{1+{\chi }_{\mathrm{\Gamma }}\left(x\right)}{{\omega }_n}{\alpha }_0\underset{D}{\int }\frac{K_0\left(x,y\right)}{{\rho }^{n-2}}\mathit{d\lambda }\left(x\right)\end{array}$$E13

where

$$\begin{array}{cc}{K}_{m+1}\left(x,y\right)=& \left(1-\pi \left(x,y\right)\right)/{\alpha }_{m+1},k_0\left(x,y\right)=b_0\left(y\right)exp\left(-\mathit{c\rho }\right)/{\alpha }_0,\\ & K_i\left(x,y\right)=b_i\left(y\right)exp\left(-\mathit{c\rho }\right)/{\alpha }_i,i=\overline{1,\dots ,m}\end{array}$$E14

We can use a method based on branching random process to Eq. (13), by applying the method for integral equations with polynomial nonlinearities.

We define a branching Markov process in $\overline{D}$, and a sequence of unbiased estimators on the trajectory of the process as follows. Assume at the initial time and at $x$ there is only one particle. With probability ${\alpha }_0$ the particle moves to the point $y\in D$, where $y=x+\rho \overline{\omega }$ where $\overline{\omega }$ is an isotropic vector in $R^n$ (if $x\in D$) and $\rho$ is distributed uniformly in $\left[0;d\right]$ with density $2\rho /d^2$. If $x\in \mathrm{\Gamma }$, then $\overline{\omega }$ is an isotropic vector in the half-space $G$, in which $D$ lies. At the point $y$, the particle is absorbed. In this case, the estimator is

$${\xi }_1\left(x\right)=\frac{d^2}{2\left(n-2\right)}exp\left(-\mathit{c\rho }\right){\tilde{b}}_0\left(x+\rho \overline{\omega }\right)/{\alpha }_0,$$E15

where

$${\tilde{b}}_0\left(x\right)=\left\{\begin{array}{ll}{b}_0\left(x\right),& x\in \overline{D}\\ 0,& x\in \overline{D}\end{array}\right.$$E16

With probability ${\alpha }_i$, $i=1,\dots ,m$, the particle at point $x$ moves to $y=\left(x+\rho \overline{\omega }\right)$ (here $\rho$ is as above) and at $y$ we generate $i$ particles.

In this case, the estimator will be

$${\xi }_1\left(y\right)=(\frac{d^2}{2\left(n-2\right)}exp\left(-\mathit{c\rho }\tilde{b}\left(x+\rho \overline{\omega }\right)/{\mathit{\text{alpha}}}_i\right)u^i\left(x+\rho \overline{\omega }\right),$$E17

here

$${\tilde{b}}_i\left(x\right)=\left\{\begin{array}{ll}{b}_i\left(x\right),& x\in \overline{D}\\ 0,& x\overline{\in }\overline{D}\end{array}\right.$$E18

With probability ${\alpha }_{m+1}$, the particle moves to the point $y$, which is distributed with density $\mathrm{\Pi }\left(x,y\right)$, and generates one particle. The estimator is thus

$${\xi }_1\left(x\right)=\left(\left(1-\pi \left(x,y\right)\right)/{\alpha }_{m+1}\right)u\left(y_1\right).$$E19

The algorithm for simulating a random point with density functions $\mathrm{\Pi }\left(x,y\right)$ is the following [3]:

1. Let $n\ge 0$ and $x_n$ be known.

2. If $x_n\in G$, then ${\mathrm{\Pi }}_{x_n}$ denote a half-space in which $D\in {\mathrm{\Pi }}_{x_n}$, and is defined by the tangent plane to $D$ at the point $x_n.$

3. Let $\omega$ is an isotropic vector in ${\mathit{Pi}}_{x_n}$, if $x_n\in G$ and an isotropic vector in $R^m$, if $x_n\in D$

4. $\rho =\rho \left(x_n,\omega \right)=\mathit{sup}\left\{t>0:x_n+\mathit{t\omega }\right\}$ is the length of a ray, that emitting from the point $x_n$ and lies in $D$ pointing in the $\omega$ direction. $x_n+\mathit{\rho \omega }$ is always a point on the boundary.

5. The Markov chain jumps from $x_n$ to $x_{n+1}$, where $x_{n+1}=\xi =x_n+\mathit{\rho \omega }$ with probability $A_1\left(\rho \right)$, and with probability $1-A_1\left(\rho \right)\xi$ is distributed on the segment $\left[x_n,x_n+\mathit{\rho \varpi }\right]$ with the probability density $p_1\left(r\right)=A_2\left(\rho \right)/\left(1-A_1\left(\rho \right)\right)$.

The simulation method, which was defined above gives a vector $\xi$, which is distributed via$\mathrm{\Pi }\left(x,y\right)$.

It is not difficult (see Ref. [3]), that ${\xi }_t\left(x\right)$ is a martingale relative to the $\sigma$-algebra ${\left\{F_t\right\}}_{t=1}^{\infty }$, where $F_t$ was generated by the process up to time $t$. Because ${\xi }_t\left(x\right)\ge 0$ (from lemma 1), ${\xi }_z\left(x\right)$ is a supermartingale, then (see Ref. [5], p. 112) with probability one if tends to the integrable random variable ${\xi }_{\infty }\left(x\right)$. To show that $\mathit{Ex}\mid {\xi }_{\infty }\left(x\right)\mid <+\infty$ and ${\mathit{xi}}_n\left(x\right)=\mathit{Ex}\left\{{\xi }_{\infty }\left(x\right)/F_n\right\}$, we must show that ${\xi }_t\left(x\right)$ is uniformly integrable. We will show that our estimator is bounded and that ${\xi }_{\infty }\left(x\right)$ is an unbiased estimator for $u\left(x\right)$. It is known [5] that if $\mathit{Ex}{\mathit{xi}}_t^2\left(x\right)<+\infty$, then ${\xi }_t\left(x\right)$ is uniformly integrable.

Theorem 1 If condition(12)is satisfied and

$${\alpha }_{m+1}\ge \underset{\left(x,y\right)\in D}{max}\left\{\left(c^2-g\left(y\right)\rho +c\left(n-3\right)\right)/(\left(c^2{p}^2+c\left(n-3\right)\right)\right\}$$E20

then ${\xi }_t\left(x\right)$ is a supermartingale and tends to integrable random variable ${\xi }_{\infty }\left(x\right)$ almost surely.

Proof: From the definition of our estimator, it is easy to see that if we choose ${\alpha }_{m+1}$ to satisfy the inequality below

$${\alpha }_{m+1}\ge \underset{x,y\in D}{max}\left\{\frac{\left(c^2-g\left(y\right)\right)\rho +c\left(n-3\right)}{c^2\rho +c\left(n-3\right)}\right\}$$E21

Then, $E_x{\xi }_t^2\left(x\right)<+\infty$ from condition (12). From (21) one can see that the inequality $0<{\alpha }_{m+1}\le 1$ holds. The parameter $c$ we could choose, so that the inequality (21) is less than one because $g\left(y\right)>0$. The rest of the ${\alpha }_i$ can be chosen differently from zero.$\square$.

Thus, ${\xi }_{\infty }\left(x\right)$ is an unbiased estimator for $u\left(x\right)$ at $x\in D$. Here, the estimator ${\xi }_{\infty }\left(x\right)$ is more amenable to practical computation. It differs from the estimator used for the nonlinear Dirichlet problem, and it is not necessary to create an $\epsilon$ -neighborhood of the boundary.

3. The case when $n=2$ and $g\left(x\right)=c^2$

Let $D$ be a bounded domain in $R^3$, with a regular smooth boundary $\mathrm{\Gamma }$. Consider the following nonlinear problem

$$\begin{array}{l}-\mathrm{\Delta }u\left(x\right)+c^{2}u\left(x\right)=b\left(x\right)u^2+f\left(x\right)\\ {\left.\frac{\partial u\left(x\right)}{\partial n_x}\right|}_{\mathrm{\Gamma }}=0\end{array}$$E22

Since $A_2\left(\rho \right)=c^2\rho exp\left(-\mathit{c\rho }\right)$ the volume integral in (5) vanishes. In this particular case, we obtain the following integral equation

$$\begin{array}{ll}u\left(x\right)=& \frac{1+{\chi }_{\mathrm{\Gamma }}\left(x\right)}{4\pi }\underset{D}{\int }\left(exp\left(-\mathit{c\rho }\right)/\rho \right)b\left(y\right)u^2\left(y\right)\mathit{dy}+\\ & +\frac{1+{\chi }_{\mathrm{\Gamma }}\left(x\right)}{4\pi }\underset{D}{\int }\left(exp\left(-\mathit{c\rho }\right)/\rho \right)f\left(y\right)\mathit{dy}+\\ & +\frac{1+{\chi }_{\mathrm{\Gamma }}\left(x\right)}{4\pi }\underset{\mathrm{\Gamma }}{\int }\left(1+\mathit{c\rho }\right)exp\left(-\mathit{c\rho }\right)\left(cos{\phi }_{\mathit{xy}}\right)/{\rho }^2)u\left(y\right)\mathit{dy}\mathrm{\Gamma }\end{array}$$E23

We denote.

$\pi \left(x,y\right)=1-exp\left(-\mathit{c\rho }\right)\left(1+\mathit{c\rho }\right)$, and $p_1\left(x,y\right)=\frac{1+X_{\mathrm{\Gamma }}\left(x\right)}{4{\mathit{\pi \rho }}^2}cos{\phi }_{\mathit{xy}}$ and transform the above integral equation to get

$$\begin{array}{ll}u\left(x\right)=& \frac{1+X_{\mathrm{\Gamma }}\left(x\right)}{2}{d}^2\underset{D}{\int }\left(2\rho exp\left(-\mathit{c\rho }\right)/4\pi d^2{\rho }^2\right)b\left(y\right)u^2\left(y\right)\mathit{dy}\\ & +\frac{1+{\chi }_{\mathrm{\Gamma }}\left(x\right)}{2}{d}^2\underset{D}{\int }\left(2\rho exp\left(-\mathit{c\rho }\right)/4\pi d^2{\rho }^2\right)f\left(y\right)\mathit{dy}+\\ & +\underset{\mathrm{\Gamma }}{\int }\left(1-\pi \left(x,y\right)\right)p_1\left(x,y\right)u\left(y\right)d_y\mathrm{\Gamma }\end{array}$$E24

From Eq. (24) it is obvious how to construct the branching Markov process and the sequence of unbiased estimators analogous to those determined above [6, 7]. Here the difference is that we do not introduce artificial parameters ${\alpha }_i$. Usually, when the parameters ${\alpha }_i$ are introduced, the estimation will change, and sometimes, if we are not careful, the variance of the estimators may be unbounded.

Assume $b\left(x\right)\ge 0$, $f\left(x\right)\ge 0$ are such that

$$max\left\{\parallel b\left(x\right){\parallel }_{C\left(\overline{D}\right)},\parallel f\left(x\right){\parallel }_{C\left(\overline{D}\right)}\right\}\le \underset{x,y\in D}{min}\left\{\sqrt{\pi \left(x,y\right)}/d^2\right\}$$E25

Here $d$ is the diameter of the domain $D$.

Now, we will determine the branching Markov process in $\overline{D}$ and the corresponding sequence of unbiased estimators. Assume that at $x$ there is only one particle. Furthermore, the point $y$ is distributed with density $p_1\left(x,y\right)$$y\in \mathrm{\Gamma }$. With probability $1-\pi \left(x,y\right)$, we generate a particle at $y$. In this case, the estimator is ${\xi }_1\left(x\right)\equiv 1\cdot u\left(y\right)$. With probability $\pi \left(x,y\right)/2$, the particle moves from $x$ to the point $y_1=x+\rho \overline{\omega }$ and is absorbed. Here, $\overline{\omega }$ an isotropic vector in $R^3$, and $\rho$ is distributed in $\left[0,d\right]$ with density $2\rho /d^2$. In this case, the estimator will be

$${\xi }_1\left(x\right)=\left(d^{2}exp\left(-\mathit{c\rho }\right)\pi \left(x,y\right)f\left(y_1\right)\right),$$E26

where

$$\tilde{f}\left(y_1\right)=\left\{\begin{array}{ll}f\left(y_1\right),& y\in \overline{D},\\ 0,& y\notin \overline{D}.\end{array}\right..$$E27

Two particles are created at $y_1$ with probability $\pi \left(x,y\right)/2$. In this case, the estimator is

$${\xi }_1\left(x\right)=\left(d^{2}exp\left(-\mathit{c\rho }\right)\overline{b}\left(y_1\right)/\pi \left(x,y\right)\right)u^2\left(y_1\right),$$E28

where

$$\tilde{b}\left(y_1\right)=\left\{\begin{array}{ll}b\left(y_1\right),& y_1\in \overline{D},\\ 0,& y_1\notin \overline{D}.\end{array}\right.$$E29

New particles behave in a manner similar to their parents. ${\xi }_t\left(x\right)={\xi }_{t-1}\left(x\right)$ if the process terminates up to time $t$, and ${\xi }_t\left(x\right)$ is obtained from ${\xi }_{t-1}\left(x\right)$ by the replacing $u\left(y\right)$ in the estimators by ${\xi }_1\left(x\right)$. Obviously, ${\xi }_t\left(x\right)$ is a martingale relative to the family of $\sigma -$ algebras, ${\left\{F_t\right\}}_{t=1}^{\infty }$, where $F_t$ is generated by the process up to time $t$.

From the conditions (25), it is easy to see that $E_x{\xi }_t^2\left(x\right)<+\infty$, so that ${\xi }_t\left(x\right)$ is a uniformly integrable martingale. Therefore, ${\xi }_t\left(x\right)$ tends to ${\xi }_{\infty }\left(x\right)$ almost surely, and $E_x\left|{\xi }_{\infty }\left(x\right)\right|<+\infty$, and ${\xi }_n\left(x\right)=E_x\left\{{\xi }_{\infty }\left(x\right)/F_n\right\}$. Thus, ${\xi }_{\infty }\left(x\right)$ is the unbiased estimator for solution to problem (22) at the point $x\in D$.

Note that $m_1$ is the average number of particles generated over one time stop. In our case, $m_1=1$ and so the constructed branching Markov process is critical in the sense of neutronics.

The next assertion is true.

Theorem 2 With probability one, the constructed branching Markov process terminates.

The proof is analogous to that of theorem [8] and is thus omitted.

We should note that the above method for constructing unbiased estimators may be used for solving mixed problems when parts of the boundary are given a Dirichlet condition and other parts of Neumann condition. It is only important, that pieces of the boundary, where the normal derivative exists, should be convex surfaces. In this case, we use the estimators with a combination of the above proposed algorithms with the process “walk on spheres with branching” (see Ref. [9]).

It is not difficult to see, that these estimators use namely absorption. Depending on the complexity of algorithms, one may use the other types of estimators (e.g., collisions estimators) (see Ref. [10]).

4. Comparison of the complexity of the estimators absorption and collision in the quasi-linear case

The comparison of the complexity absorption and collision estimators in the quasi-linear case is a difficult problem. However, in some particular cases, one can analyze the complexity in this case.

To optimize the calculation procedure we can sometimes choose different Markov chains and estimators. In this case, the choice of the estimator is most important and could increase the efficiency of the algorithm. It is often useful to compare two estimators and chose one of them base on the particular problem being solved. We note that the computational cost may depend on the type of equation and the type of Markov process. It may turn out that one estimator is best in one situation while the other maybe best for a difficult equation.

For quasi-linear equations we will give some recommendations in choosing estimators that minimize the complexity of a calculation (see Ref. [11]).

Let consider the following nonlinear integral equation

$$\begin{array}{ll}\phi \left(x\right)=& f\left(x\right)+\\ & +\sum _{m=1}^n\underset{D}{\int }\dots \underset{D}{\int }{K}_m\left(x,y_1,\dots ,y_m\right)\prod _{i=1}^m\phi \left(y_1\right){\mathit{dy}}_1\dots {\mathit{dy}}_i\end{array}$$E30

where $D\subset R^n$, $K_i(x,y_1,\dots ,y_m)$, and $f\left(x\right)$ are given functions in $n$ times $D\times \cdots \times D$, and $\phi \left(x\right)$ is the unknown function. These problems, which are connected with the construction of unbiased estimators for the functional $\left(h,\phi \right)$, where $h\left(x\right)$ is a given function and was discussed in Refs [10, 12]. We assume that $K_i\left(x,y_1,\dots ,y_i\right)$, and $f\left(x\right)$ satisfy the condition in [10, 12].

It’s known that the solution of nonlinear integral equations by Monte Carlo is connected with branching Markov processes, determined by the initial distribution density $P_0\left(x\right)\ge 0$, and the transition density $P_i^1\left(x,y_1,\dots ,y_i\right)$. Here

$$\underset{0}{\int }{P}_0\left(x\right)\mathit{dx}=1,$$E31

$$\sum _{i=1}^n\underset{D}{\int }\dots \underset{D}{\int }{P}_i^1\left(x,y_1,\dots ,y_i\right){\mathit{dy}}_i\dots {\mathit{dy}}_i=1-g^1\left(x\right),$$E32

where $g^1\left(x\right)\in \left[0,1\right)$ is the probability of terminating the branching Markov process trajectory. Realizations of this process are connected with “trees,” on which were constructed different unbiased estimators.

We consider only estimators using absorption or collision and will give recommendations for their use based on complexity.

As is known [13], the integral equation

$$a\left(x\right)=\sum _{m=1}^n\underset{D}{\int }\dots \underset{D}{\int }{P}_m^1\left(x,y_1,\dots ,y_m\right)\prod _{i=1}^{m}a\left(y_i\right){\mathit{dy}}_1,\dots {\mathit{dy}}_i+g^1\left(x\right)$$E33

has a unique positive solution for all $x\in D$ which is determined through a Neumann series. With probability one, the number of generations in the branching Markov process is bounded.

If the branching Markov process in the branches splits into $m$- new particles with given probabilities ${\alpha }_m$, and the probability of absorption is $g^1\left(x\right)={\alpha }_0$, then probabilities ${\alpha }_m$ satisfy the conditions.

$$\sum _{i=1}^n\underset{D}{\int }\dots \underset{D}{\int }{\alpha }_i{P}_i\left(x,y_1,\dots ,y_i\right){\mathit{dy}}_1\dots {\mathit{dy}}_i=1-{\alpha }_0,$$E34

where the functions $P_1\left(x,y_1,\dots ,y_i\right)$ satisfy

$$\underset{D}{\int }\dots \underset{D}{\int }{P}_i\left(x,y_1,\dots ,y_i\right){\mathit{dy}}_1\dots {\mathit{dy}}_i=1.$$E35

It is not difficult to see that (see sec.1) that the function $M\left(x\right)$, which is the mean number of particles in the branching process, satisfies the linear integral equation

$$M\left(x\right)={\alpha }_0+\sum _{i=1}^n{\alpha }_i\underset{D}{\int }\dots \underset{D}{\int }{P}_i\left(x,y_1,\dots ,y_i\right)\sum _{j=1}^{i}M\left(y_j\right){\mathit{dy}}_1{\mathit{dy}}_i.$$E36

Here $M\left(x\right)$ the iterative solution to this equation with the initial value ${\alpha }_0$. It’s known (see Section 1) that the constructed branching Markov process will be terminated in condition $b={\sum }_{m=1}^{n}m{\alpha }_m\le 1$, this is why we choose ${\alpha }_m$ so that $b\le 1$, when $b$ is the average number of particles generated in one step.

We solve the integral Eq. (30) at the point $x\in D$. When the function $h$ is a Dirac $\delta$-function, that is $\int h\left(x-y\right)\phi \left(x\right)\mathit{dx}=\phi \left(y\right)$. In this case, at the initial time, we have one particle at the point $x$ almost, surely. The branching process is constructed on a tree of branches then we will construct unbiased estimators using absorption, ${\xi }_n\left(x\right)$, and collision, ${\xi }_c\left(x\right)$. To compare the complexity of these we compute on this tree.

Obviously, the number of arithmetic operations depends on the average number of particles in the tree. Furthermore, we will denote $D_n$ as the variance of the absorption estimator, $D_c$ and the variance of the collision estimator. The complexity of the absorption estimator is $T_n=t^{\prime }\cdot D_n$ and the complexity collision estimator $D_c={t^{\prime }}^{\prime }\cdot D_c$. Here $t^{\prime }$ and ${t^{\prime }}^{\prime }$ as the average CPU time for calculation estimators, which will be needed for calculation $T_n$ and $T_c$.

First, we will investigate the variance of both estimators.

Let $F^{+}$ and $F^{-}$ be a set of positive and negative bounded functions. We shall consider the following integral equation

$$\begin{array}{ll}\psi \left(x\right)=& \frac{f^2\left(x\right)}{{\alpha }_0}+\\ & +\sum _{m=1}^n\underset{D}{\int }\dots \underset{D}{\int }\frac{K_m^2\left(x,y_1,\dots ,y_m\right)}{{\alpha }_m{p}_m\left(x,y_1,\dots ,y_m\right)}\prod _{i=1}^m\psi \left(y_i\right){\mathit{dy}}_1,\dots ,{\mathit{dy}}_i\end{array}$$E37

$$\begin{array}{ll}{\psi }_1\left(x\right)=& f\left(x\right)\left(2\phi \left(x\right)-f\left(x\right)\right)+\\ & +\sum _{m=1}^n\underset{D}{\int }\dots \underset{D}{\int }\frac{K_m^2\left(x,y_1,\dots ,y_m\right)}{{\alpha }_m{p}_m\left(x,y_1,\dots ,y_m\right)}\prod _{i=1}^m{\psi }_1\left(y_i\right){\mathit{dy}}_1,\dots ,{\mathit{dy}}_i\end{array}$$E38

We suppose that the minimal solutions of Eqs. (37), (38) exist and that Neumann series converges for (37) with the initial function ${\psi }^{\left(0\right)}\left(x\right)=\frac{f^2\left(x\right)}{{\alpha }_0}$ and for (38) with the initial function ${\psi }^{\left(0\right)}\left(x\right)=f\left(x\right)\left(2\phi \left(x\right)-f\left(x\right)\right)$. Then the expression for the variances will be

$$D_n=\left(h^2/p_0,\psi \right)-{\left(h,\phi \right)}^2,D_c=\left(h^2/p_0,{\psi }_1\right)-{\left(h,\phi \right)}^2.$$E39

Let the function $G\left(x\right)\in F^{+}$, where $G\left(x\right)=f\left(x\right)\left(2\phi \left(x\right)-f\left(x\right)\right)$. Moreover, if $f\left(x\right)\ge 0$ then $G\left(x\right)\in F^{+}$, in this case, the following hold:

Theorem 3 Let$f\left(x\right)\in F^{+}$. If the probability of terminating the trajectory of our branching processes is

$${\alpha }_0\in \left[0,minG\left(x\right),x\in D\right],$$E40

then $D_c\le D_n$.

Proof: Since the values of $D_n$ and $D_c$ depend on solutions of Eqs. (37), (38), we will investigate the minimal solutions of these equations. The minimal solutions of these equations continuously depend on right-hand sides of the equations. Consequently, if ${\alpha }_0<min\left(f\left(x\right)\left(2\phi \left(x\right)-f\left(x\right)\right)\right)$, then ${\psi }_1\left(x\right)\le \psi$. From here, $D_c\le D_n$$\square$.

Consequence: If it follows that ${\alpha }_0\notin \left[0,minG\left(x\right),x\in D\right]$, then $D_n\le D_c$.

Theorem 4 Suppose that$f\left(x\right)\in F^{-}$. If the probability of terminating a trajectory of our branching processes is

$${\alpha }_0\in \left[0,minG\left(x\right),x\in D\right],$$E41

then $D_n\le \mathit{Dc}$.

The proof is analogous to the above.

Thus, using a priori information about the variance of estimators (solutions of the Eqs. (37), (38)), one can give recommendations on choosing optimal estimators.