Branching diffusion representation of semi-linear elliptic PDEs and estimation using Monte Carlo method

TL;DR

This paper introduces a probabilistic branching diffusion approach to solve semi-linear elliptic PDEs with polynomial non-linearity, enabling solution estimation via Monte Carlo methods without assuming prior existence of solutions.

Contribution

It extends branching diffusion representations to PDEs with non-linear gradient terms and provides explicit conditions for their validity, advancing probabilistic solution techniques.

Findings

Probabilistic representation for PDE solutions using branching processes.

Explicit conditions under which the representations are valid.

Successful Monte Carlo estimation of solutions in multi-dimensional cases.

Abstract

We study semi-linear elliptic PDEs with polynomial non-linearity and provide a probabilistic representation of their solution using branching diffusion processes. When the non-linearity involves the unknown function but not its derivatives, we extend previous results in the literature by showing that our probabilistic representation provides a solution to the PDE without assuming its existence. In the general case, we derive a new representation of the solution by using marked branching diffusion processes and automatic differentiation formulas to account for the non-linear gradient term. In both cases, we develop new theoretical tools to provide explicit sufficient conditions under which our probabilistic representations hold. As an application, we consider several examples including multi-dimensional semi-linear elliptic PDEs and estimate their solution by using the Monte Carlo method.

Tables

Table 1. Table 5.1: Numerical results for the unbiased estimator of probabilistic representation of one-dimensional semi-linear PDE ( 33 ).

$x$	Estimate	99% conf. interval	Std. Dev./Mean	Relative error	Run time (secs)
0	1.4144	[1.4134, 1.4153]	0.2644	0.0112%	13
-0.2	1.3859	[1.3852, 1.3866]	0.1872	0.0358%	26

Table 2. Table 5.2: Numerical results for the unbiased estimator for different values of r 𝑟 r and x = 0 . 𝑥 0 x=0.

$r$	99% conf.	Std. Dev./Mean	Run time
	interval		(secs)
0.4	[1.413654, 1.417886]	0.5793	17
0.5	[1.410199, 1.417666]	0.9134	23

Table 3. Table 5.3 : Numerical results for the unbiased estimator of probabilistic representation for domain ( − 0.9 , 0.9 ) . 0.9 0.9 (-0.9,0.9).

$x$	Estimate	99% conf. interval	Std. Dev./Mean	Run time (secs)
0	0.9597	[0.9595, 0.9600]	0.0746	123
-0.2	0.9612	[0.9611, 0.9614]	0.0722	121

Table 4. Table 5.4: Numerical results for the unbiased estimator of probabilistic representation of one-dimensional semi-linear PDE ( 35 ).

$x$	Estimate	99% conf. interval	Std. Dev./Mean	Relative error	Run time (secs)
0	0.9999	[0.9989, 1.0001]	0.3931	0.0097%	10
-0.1	1.0198	[1.0190, 1.0205]	0.2745	0.0351%	12

Table 5. Table 5.5: Numerical results for the unbiased estimator of probabilistic representation of one-dimensional semi-linear PDE ( 35 ).

$x$	Estimate	99% conf. interval	Std. Dev./Mean	Relative error	Run time (secs)
0.0	1.0006	[0.9971, 1.0040]	1.3378	0.0575%	12
-0.1	1.0203	[1.0168, 1.0234]	1.3641	0.0209%	17

Table 6. Table 5.6: Numerical results for the unbiased estimator of probabilistic representation of semi-linear PDE ( 36 ) in d = 2 𝑑 2 d=2 with β = 2 . 𝛽 2 \beta=2.

$x$	Estimate	99% conf. interval	Std. Dev./Mean	Relative error	Run time (secs)
(0.0,0.0)	-0.0001	[-0.0021, 0.0019]	–	–	30
(0.1,0.0)	0.0998	[0.0978, 0.1018]	5.4883	0.5726%	34
(0.2,0.1)	0.3095	[0.3075, 0.3114]	1.7342	0.0426%	30
(0.2,0.2)	0.4245	[0.4226, 0.4265]	1.2646	0.4144%	25

Table 7. Table 5.7: Numerical results for the unbiased estimator of probabilistic representation of semi-linear PDE ( 36 ) in d = 4 𝑑 4 d=4 with β = 4 . 𝛽 4 \beta=4.

$x$	Estimate	99% conf. interval	Std. Dev./Mean	Relative error	Run time (secs)
(0,0,0,0)	0.0000	[-0.0010, 0.0012]	–	–	190
(0.1,0,0,0)	0.1009	[0.0999, 0.1019]	2.6185	0.5584%	309
(0.1,0.1,0,0)	0.2021	[0.2012, 0.2031]	1.2986	0.2924%	413
(0.1,0.1,0.1,0)	0.3097	[0.3087, 0.31072]	0.8750	0.1305%	510

Equations

$$\displaystyle{\mathcal{L}}u+f\bigl{(}u,Du\bigr{)}=0~{}~{}\text{in }{\mathcal{O}},\quad~{}u=h~{}~{}\text{on }\partial{\mathcal{O}},$$

$${\partial }_{t}u+\mathcal{L}u+\beta \left(\sum _{l\in \mathbb{N}}{p}_l{u}^l-u\right)=0\text{in }[0,T)\times {\mathbb{R}}^d,u(T,\cdot )=g\text{in }{\mathbb{R}}^d,$$

$${\partial }_{t}u+\mathcal{L}u+f(u,Du)=0\text{in }[0,T)\times {\mathbb{R}}^d,u(T,\cdot )=g\text{in }{\mathbb{R}}^d,$$

$$\mathcal{L}u+\beta \left(\sum _{l\in \mathbb{N}}{c}_l{u}^l-u\right)=0\text{in }\mathcal{O},u=h\text{on }\partial \mathcal{O},$$

$$f(x,y,z):=\sum _{l=(l_0,l_1,\dots ,l_m)\in L}{c}_l(x)y^{l_0}\prod _{i=1}^m(b_i(x)\cdot z{)}^{l_i},$$

$$\displaystyle{\mathcal{L}}u+\beta\left(f\bigl{(}u,Du\bigr{)}-u\right)=0~{}~{}\text{in }{\mathcal{O}},\quad~{}u$$

$$X_t^x=x+{\int }_0^t\mu (X_s^x)ds+{\int }_0^t\sigma (X_s^x)d{W}_s,\mathbb{P}-\text{a.s.},$$

$$\displaystyle T^{\emptyset}:=\tau\wedge\eta^{x}\quad\text{where}~{}~{}\eta^{x}:=\inf\Big{\{}t\geq 0;\,X^{x}_{t}\notin{\mathcal{O}}\Big{\}}.$$

$$X_t^k=X_{T_{k^{-}}}^{k^{-}}+{\int }_{T_{k^{-}}}^t\mu (X_s^k)ds+{\int }_{T_{k^{-}}}^t\sigma (X_s^k)d{W}_{s-T_{k-}}^k,\mathbb{P}-\text{a.s.}$$

$$\displaystyle T_{k}:=\left(T_{k^{-}}+\tau^{k}\right)\wedge\inf\Big{\{}t\geq T_{k^{-}};\,X^{k}_{t}\notin{\mathcal{O}}\Big{\}}.$$

$${\mathcal{F}}_n:=\sigma \left({\tau }^k,I^k,W^k,k\in \{\varnothing \}\cup \bigcup _{i=1}^n{\mathbb{N}}^i\right),n\in \mathbb{N}.$$

$$\beta \left(\sum _{l\in L}|l|p_l-1\right)-\frac{{\lambda }_1}{2}\le 0,$$

$$\mathcal{L}u+\beta \left(f(u)-u\right)=0\text{in }\mathcal{O},u=h\text{on }\partial \mathcal{O},$$

$${\psi }^x:=\prod _{\begin{array}{c}k\in {\mathcal{K}}^x\\ X_{T_k}^k\notin \mathcal{O}\end{array}}h(X_{T_k}^k)\prod _{\begin{array}{c}k\in {\mathcal{K}}^x\\ X_{T_k}^k\in \mathcal{O}\end{array}}\frac{c_{I^k}(X_{T_k}^k)}{p_{I^k}}.$$

$${\psi }_n^x:=\prod _{\begin{array}{c}k\in {\cup }_{i=0}^n{\mathcal{K}}_i^x\\ X_{T_k}^k\notin \mathcal{O}\end{array}}h(X_{T_k}^k)\prod _{\begin{array}{c}k\in {\cup }_{i=0}^n{\mathcal{K}}_i^x\\ X_{T_k}^k\in \mathcal{O}\end{array}}\frac{c_{I^k}(X_{T_k}^k)}{p_{I^k}}\prod _{k\in {\mathcal{K}}_{n+1}^x}u(X_{T_{k^{-}}}^k),$$

$$\displaystyle u(x)={\mathbb{E}}\Bigl{[}e^{-\beta{\eta^{x}}}h(X^{x}_{\eta^{x}})+\int^{\eta^{x}}_{0}\beta e^{-\beta s}f(\cdot,u)(X^{x}_{s})\,ds\Bigr{]}.$$

$$u(x)$$

$$=\mathbb{E}\left[h(X_{{\eta }^x}^x)1_{\tau \ge {\eta }^x}+\frac{c_I(X_{\tau }^x)}{p_I}{u}^I(X_{\tau }^x)1_{\tau <{\eta }^x}\right].$$

$$\displaystyle u(X^{k}_{T_{k^{-}}})={\mathbb{E}}\left[h(X^{k}_{T_{k}})\mathbf{1}_{X^{k}_{T_{k}}\notin{\mathcal{O}}}+\frac{c_{I^{k}}(X^{k}_{T_{k}})}{p_{I^{k}}}u^{I^{k}}(X^{k}_{T_{k}})\mathbf{1}_{X^{k}_{T_{k}}\in{\mathcal{O}}}\,\bigg{|}\,{\mathcal{F}}_{0}\right].$$

$$u(x)=\mathbb{E}\left[h(X_{{\eta }^x}^x)1_{\tau \ge {\eta }^x}+\frac{c_I(X_{\tau }^x)}{p_I}\prod _{i=0}^{I-1}{\psi }_i^{X_{\tau }^x}{1}_{\tau <{\eta }^x}\right].$$

$${\psi }_i^{X_{\tau }^x}:=\prod _{\begin{array}{c}k=(i,\dots )\in {\mathcal{K}}^x\\ X_{T_k}^k\notin \mathcal{O}\end{array}}h(X_{T_k}^k)\prod _{\begin{array}{c}k=(i,\dots )\in {\mathcal{K}}^x\\ X_{T_k}^k\in \mathcal{O}\end{array}}\frac{c_{I^k}(X_{T_k}^k)}{p_{I^k}}.$$

$$\displaystyle{\mathbb{E}}\Big{[}\prod_{i=0}^{I-1}\psi^{X^{x}_{\tau}}_{i}\,\Big{|}\,X^{x}_{\tau},I\Big{]}\mathbf{1}_{\tau<\eta^{x}}=u^{I}\left(X^{x}_{\tau}\right)\mathbf{1}_{\tau<\eta^{x}}.$$

$$\displaystyle u(x)={\mathbb{E}}\Bigl{[}e^{-\beta{\eta^{x}}}h(X^{x}_{\eta^{x}})+\int^{\eta^{x}}_{0}\beta e^{-\beta s}f(\cdot,u)(X^{x}_{s})\,ds\Bigr{]}.$$

$$\displaystyle u(x)={\mathbb{E}}\Bigl{[}e^{-\beta({\eta^{x}}\wedge\delta)}u(X^{x}_{{\eta^{x}}\wedge\delta})+\int^{{\eta^{x}}\wedge\delta}_{0}\beta e^{-\beta s}f(\cdot,u)(X^{x}_{s})\,ds\Bigr{]}.$$

$$\mathbb{E}\left[e^{-\beta ({\eta }^x\wedge \delta )}\phi (X_{{\eta }^x\wedge \delta }^x)\right]=\phi (x)+\mathbb{E}\left[{\int }_0^{{\eta }^x\wedge \delta }{e}^{-\beta s}\left(\mathcal{L}\phi -\beta \phi \right)(X_s^x)ds\right].$$

$$\mathbb{E}\left[\frac{1}{\delta }{\int }_0^{{\eta }^x\wedge \delta }{e}^{-\beta s}\left(\mathcal{L}\phi -\beta \phi +\beta f(\cdot ,u)\right)(X_s^x)ds\right]\ge 0.$$

$$\mathcal{L}\phi (x)+\beta \left(f(x,\phi (x))-\phi (x)\right)\ge 0.$$

$$\mathcal{L}v+\beta \left(\sum _{l\in L}\left|c_l\right|v^l-v\right)=0\text{in }\mathcal{O},v=\left|h\right|\text{on }\partial \mathcal{O}.$$

$$u(x)=\mathbb{E}\left[\prod _{\begin{array}{c}k\in {\mathcal{K}}^x\\ X_{T_k}^k\notin \mathcal{O}\end{array}}\frac{e^{-\beta \Delta T_k}h(X_{T_k}^k)}{\bar{F}(\Delta T_k)}\prod _{\begin{array}{c}k\in {\mathcal{K}}^x\\ X_{T_k}^k\in \mathcal{O}\end{array}}\frac{\beta e^{-\beta \Delta T_k}{c}_{I^k}(X_{T_k}^k)}{p_{I^k}\rho (\Delta T_k)}\right],$$

$$\mathcal{L}\phi -\beta \phi =0\text{in }\mathcal{O},\phi =h\text{on }\partial \mathcal{O}.$$

$${\partial }_t\chi -\mathcal{L}\chi$$

Full text

Key words: automatic differentiation formula, branching diffusion processes, elliptic, exit time, Monte Carlo method, partial differential equation, semi-linear

AMS subject classifications (2010): 35J61, 60H30, 60J85, 65C05

1   Introduction

In this paper, we are interested in the following class of semi-linear elliptic partial differential equations (PDEs): given a bounded domain ${\mathcal{O}}\subset\mathbb{R}^{d}$, $f:{\mathcal{O}}\times\mathbb{R}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$, $h:\partial{\mathcal{O}}\rightarrow\mathbb{R}$,

[TABLE]

where ${\mathcal{L}}$ is the infinitesimal generator of a diffusion. This class of PDEs arises naturally in many fields of science (see, e.g., [2, 25] and references therein). In most cases, the related PDEs only involve non-linearity in the unknown function but not in its derivatives. For instance, the Poisson–Boltzmann equation or the steady states in nonlinear equations of the Klein–Gordon or Schrödinger type. Besides, the stationary Burgers equation or Hamilton-Jacobi-Bellman equations when the diffusion coefficient is uncontrolled are examples of PDEs where the non-linearity involves both the unknown function and its first order derivatives. In this paper, we provide existence results for this class of PDEs when the generator $f$ is polynomial and derive a new probabilistic representation of the solution which is well suited for numerical application, especially in high dimensions, using the Monte Carlo method.

The classical probabilistic approach for semi-linear PDEs relies on the theory of backward stochastic differential equations (BSDEs) initiated by Pardoux and Peng [30]. The case of elliptic PDE with Dirichlet condition was first studied by Darling and Pardoux [10] for Lipschitz generator and later more general results were provided by Pardoux [29] and Briand et al. [8]. These results rely on the so-called monotonicity assumption, which requires $y\mapsto f(x,y,z)-\mu y$ to be decreasing for some constant $\mu\in\mathbb{R}$, and the assumption that $f$ has a linear or quadratic growth in $z$. We do not put such restrictions in our setting and propose a method to investigate a different class of elliptic PDEs with a generator which is a polynomial in $y$ and $z$.

Our probabilistic representation is based on branching diffusion processes. These processes describe the evolution of a population of independent and identical particles moving according to a diffusion process. They were first introduced by Skorokhod [34] and later studied, more thoroughly and systematically, in a series of papers by Ikeda et al. [20, 21, 22]. In particular, these authors, using branching diffusion processes, established a probabilistic representation of semi-linear parabolic PDEs of the form

[TABLE]

where $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$, $\beta>0$ and $(p_{l})_{l\in\mathbb{N}}$ is a probability mass function. More precisely, they consider a process where each particle moves according to a diffusion with generator ${\mathcal{L}}$ and dies at an exponentially distributed random time with parameter $\beta$ to give birth to $l$ offsprings with probability $p_{l}$. Then, the solution of PDE (2) is expressed as ${\mathbb{E}}[\prod_{i=1}^{N_{T}}g(X_{T}^{i})]$ where $N_{T}$ and $(X_{T}^{i})_{i=1,\ldots,N_{T}}$ denote the number of particles and their positions at time $T$ respectively. The special case $p_{2}=1$, which corresponds to the celebrated Fisher–KPP equation, has been particularly well-studied (see, e.g., McKean [26]).

In Henry-Labordère et al. [18], the authors extended the probabilistic representation of (path-dependent) PDE (2) where the probability mass function $(p_{l})_{l\in\mathbb{N}}$ is replaced by an arbitrary real-valued sequence of functions. See also Rasulov et al. [32]. This result was further extended in Henry-Labordère et al. [19] to parabolic semi-linear PDE of the form

[TABLE]

where $f$ is a $(d+1)$–variate polynomial. The main idea consists of introducing marked particles which carry a weight function given by the Bismuth-Elworthy-Li formula to account for non-linearity in gradient of the solution. Under conditions of “small non-linearity” or small maturity $T$, the authors showed that their probabilistic representation provides a continuously differentiable viscosity solution to PDE (3). They also performed numerical simulations to illustrate the accuracy of their method to solve PDEs using the Monte Carlo method.

The aim of this paper is to extend the results in Henry-Labordère et al. [18, 19] to the case of elliptic PDEs with Dirichlet condition. The first result on the link between semi-linear elliptic PDE and branching diffusion processes was obtained by Watanabe [36] who derived a criterion for extinction of branching Brownian motion absorbed at the boundary of a domain. Recently, Bossy et al. [4] extended this result to derive a probabilistic representation for PDE of the form

[TABLE]

where $(c_{l})_{l\in\mathbb{N}}$ is a sequence of real-valued functions. They used it to compute a Monte Carlo approximation of the solution to Poisson-Boltzmann equation in dimension three. One of the critical assumptions in their work is the existence of a smooth solution to PDE (4). In this paper, we do not require this assumption as we show directly that the probabilistic representation provides a continuous viscosity solution to PDE (4). The main difficulty compared to [18] is to ensure the continuity and integrability of the probabilistic representation. In particular, the arguments of small maturity used in [18] cannot be exploited here and instead, we resolve this problem by developing new tools which ensure that our probabilistic representation holds in small domain.

In the second part of our paper, we perform a rigorous analysis for the case of semi-linear elliptic PDE involving non-linearity in gradient of the unknown function. To the best of our knowledge, this is the first paper in the literature to provide a representation for this class of PDE using branching diffusion processes. In contrast with Henry-Labordère et al. [19], Malliavin calculus cannot be used in this setting since the exit time of a diffusion from a domain is not differentiable in the Malliavin sense. However, Delarue [11] and Gobet [16] have established a suitable automatic differentiation formula, based on the work of Thalmaier [35], which allows us to derive a probabilistic representation analogous to [19]. Under the assumption that the non-linear gradient term vanishes at the boundary, it provides a continuously differentiable solution to PDE (1) when $f$ is a multivariate polynomial.

Importantly, our probabilistic representation provides a means to evaluate solutions of semi-linear elliptic PDEs by using the Monte Carlo method whose accuracy relies on the dimensionless central limit theorem. For this reason, it is particularly well suited for numerical applications in high dimension where deterministic methods usually fail to provide accurate estimates. Regarding the probabilistic approach, although different numerical methods for BSDEs have been introduced (see, e.g., [6, 37]), to the best of our knowledge, they have never been used in the literature for elliptic PDEs. 111A related problem has been studied by Bouchard and Menozzi [5] for the case of parabolic PDEs in a cylindrical domain with finite time horizon. In addition, the classical approach to solve BSDEs induces a discretization bias and involves a sequence of conditional expectation estimators. The available methods to estimate conditional expectations, such as least-squares regression (see Gobet et al. [17]), induce an additional approximation bias and are computationally expensive in high dimension.

We consider several numerical examples, including multi-dimensional elliptic PDEs, to illustrate the accuracy of our method. One of the critical steps in the algorithm consists of simulating the first exit time and position of a diffusion from a domain. There are three classical approaches to perform this task – walk on spheres scheme, walk on parallelepipeds scheme and Euler discretization method (see Bouchard et al. [7] and references therein). For instance, the walk on spheres scheme was used by Bossy et al. [4] to estimate the solution of Poisson-Boltzmann equation in dimension three. In this paper, we consider numerical examples which involve Brownian motion over a rectangular domain and therefore, we generate unbiased samples for exit time and position using the walk on squares scheme as introduced in Faure [12] and in Milstein and Tretyakov [27], and implemented in the library developed by Lejay [24].

The rest of the paper is organized as follows. In the next section, we provide a precise formulation of the problem and introduce branching diffusion processes used to derive our probabilistic representation. Then, we consider the case of semi-linear PDE with linear gradient term in Section 3 and give explicit sufficient conditions to ensure that the probabilistic representation provides a continuous viscosity solution to the PDE. In Section 4, we perform the analysis for semi-linear PDE with non-linear gradient term. In particular, we establish appropriate automatic differentiation formulas. Finally, in Section 5, we present several numerical examples to illustrate the applicability of our results in different settings. We provide the proof of the automatic differentiation formula in Appendix A.

1.1   Notations

Any element $x\in\mathbb{R}^{d},d\geq 1,$ is a column vector with $i$th component $x_{i}$ and Euclidean norm $|x|.$ $x\cdot y$ denotes the usual dot product and $d(x,y)$ denotes the Euclidean distance for any $x,y\in\mathbb{R}^{d}.$ Given an open set ${\mathcal{O}}\subset\mathbb{R}^{d},$ $\bar{{\mathcal{O}}}$ denotes its closure, $\partial{\mathcal{O}}$ its boundary and $\mathrm{diam}({\mathcal{O}})$ its diameter. For any $g:{\mathcal{O}}\rightarrow\mathbb{R}$, the supremum norm is defined as $\|g\|_{\infty}=\sup\{|g(x)|:x\in{\mathcal{O}}\}.$ ${\mathcal{C}}({\mathcal{O}})$ (resp. ${\mathcal{C}}(\bar{{\mathcal{O}}})$) denotes the set of continuous functions on ${\mathcal{O}}$ (resp. $\bar{{\mathcal{O}}}$). ${\mathcal{C}}^{k}({\mathcal{O}})$ (resp. ${\mathcal{C}}^{k}(\bar{{\mathcal{O}}})$) denotes the set of functions with continuous derivatives of all orders less than or equal to $k$ (resp. which further have continuous extensions to $\bar{{\mathcal{O}}}$). ${\mathcal{C}}^{k,\alpha}({\mathcal{O}})$ (resp. ${\mathcal{C}}^{k,\alpha}(\bar{{\mathcal{O}}})$) are subspaces of ${\mathcal{C}}^{k}({\mathcal{O}})$ (resp. ${\mathcal{C}}^{k}(\bar{{\mathcal{O}}})$) consisting of functions whose partial derivatives of order less than or equal to $k$ are locally (resp. globally) Hölder continuous with exponent $0<\alpha<1.$ Given $T>0$, ${\mathcal{C}}^{1,k}((0,T]\times{\mathcal{O}})$ denotes the set of functions with continuous first derivative with respect to the first component and continuous derivatives of all orders less than or equal to $k$ with respect to the second component. For a smooth function $u(t,x),$ $Du$ and $D^{2}u$ stand for, respectively, its gradient (as a row vector) and Hessian matrix with respect to its second component. Furthermore, if $g:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ is a differentiable function, its gradient $Dg=(\partial_{x_{1}}g(x),\ldots,\partial_{x_{m}}g(x))$ takes values in $\mathbb{R}^{n}\otimes\mathbb{R}^{m}.$ For $d\geq 1,\mathbb{M}^{d}$ denotes the set of all $d\times d$ matrices and $I_{d}\in\mathbb{M}^{d}$ the identity matrix.

2   Problem Definition

2.1   A Class of Semi-Linear Elliptic PDE

Let $(\mu,\sigma):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\times\mathbb{M}^{d}$ denote the drift and diffusion coefficient. Then for a non-negative integer $m,$ we consider a subset $L\subset\mathbb{N}^{m+1}$ and a generator function $f:\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ defined as

[TABLE]

where $(c_{l})_{l\in L},c_{l}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ and $(b_{i})_{i=1,\ldots,m},b_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ are sequences of functions. Furthermore, for every $l=(l_{0},l_{1},\ldots,l_{m}),$ we denote $|l|:=\sum^{m}_{i=0}l_{i}.$ Given a bounded domain ${\mathcal{O}}\subset\mathbb{R}^{d}$, we consider the following semi-linear elliptic PDE:

[TABLE]

where $\beta$ is a positive constant, $h:\partial{\mathcal{O}}\rightarrow\mathbb{R}$ is the Dirichlet boundary condition and ${\mathcal{L}}$ is the infinitesimal generator associated to a diffusion process with parameters $(\mu,\sigma)$. Under a set of general assumptions, we provide a probabilistic representation of solution $u$ using the theory of branching diffusion processes. We list the assumptions on parameters of PDE (5) which are needed for our results at the outset.

Assumption 2.1**.**

(i)* The functions $b_{i}$, $c_{l}$ are continuous on $\bar{{\mathcal{O}}}$.

(ii) The function $h$ is continuous on $\partial{\mathcal{O}}$.

(iii) The function $(x,y,\zeta)\mapsto\sum_{l\in L}c_{l}(x)y^{l_{0}}\prod^{m}_{i=1}\zeta_{i}^{l_{i}}$ is continuous on $\bar{{\mathcal{O}}}\times\mathbb{R}\times\mathbb{R}^{m}$.*

Part (iii) above implies that the multivariate power series $(y,\zeta)\mapsto\sum_{l\in L}c_{l}(x)y^{l_{0}}\prod^{m}_{i=1}\zeta_{i}^{l_{i}}$ has an infinite radius of convergence for all $x\in\bar{{\mathcal{O}}}$. It is a simplifying assumption which avoids the use of localization arguments in the rest of the paper.

2.2   Marked Branching Diffusion Processes with Absorption

A (age-dependent) marked branching diffusion process is characterized by diffusion parameters $(\mu,\sigma)$, a probability density function $\rho:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$, and a probability mass function $(p_{l})_{l\in L}.$ In this process, we start with one particle of mark [math] at position $x\in\mathbb{R}^{d}$ which undergoes a diffusion with parameters $(\mu,\sigma)$ during its lifetime distributed according to $\rho.$ At the end of its lifetime (arrival time), the particle dies and gives rise to $|l|=\sum_{i=0}^{m}l_{i}$ offsprings with probability $p_{l}$, among which $l_{0}$ have mark 0, $l_{1}$ have mark 1, and so on. After their birth, each offspring performs the same but an independent branching diffusion process as the parent particle. Additionally, we consider that particles are absorbed at the boundary, i.e., they die without giving rise to any offspring when they leave the domain ${\mathcal{O}}$. In order to construct the above process, we denote by ${\mathbb{K}}:=\{\emptyset\}\cup\bigcup_{n\geq 1}{{\mathbb{N}}^{n}}$ the set of labels and we consider the probability space $(\Omega,{\mathcal{F}},{\mathbb{P}})$ equipped with

• •

a sequence of i.i.d. positive random variables $(\tau^{k})_{k\in{\mathbb{K}}}$ distributed with density function $\rho$

• •

a sequence of i.i.d. random elements $(I^{k})_{k\in{\mathbb{K}}}$ with ${\mathbb{P}}(I^{k}=l)=p_{l},l\in L$

• •

a sequence of independent $d$-dimensional Brownian motions $(W^{k})_{k\in{\mathbb{K}}}$

Furthermore, we consider the sequences $(\tau^{k})_{k\in{\mathbb{K}}},(I^{k})_{k\in{\mathbb{K}}}$ and $(W^{k})_{k\in{\mathbb{K}}}$ to be mutually independent. The age-dependent branching process is constructed as follows:

Start from one particle at position $x\in\mathbb{R}^{d}$ and index it by label $\emptyset.$ This is the common ancestor of all particles, the only particle of generation 0. For notational convenience, we simply write $\tau$, $I$ and $W$ instead of $\tau^{\emptyset}$, $I^{\emptyset}$ and $W^{\emptyset}$. The dynamic of particle $\emptyset$ is given as follows:

• •

The position $X^{\emptyset}=X^{x}$ of the particle during its lifetime is given by

[TABLE]

• •

The arrival time of the particle is given by

[TABLE]

• •

At the arrival time, if $\eta^{x}\leq\tau$, then the particle $\emptyset$ dies without giving rise to any offspring, else the particle $\emptyset$ dies and gives rise to $|I|$ offsprings which belong to the first generation, and are indexed by label $i$ for $i=0,\ldots,|I|-1.$

• •

Given $I=(I_{0},I_{1},\ldots,I_{m}),$ we have $|I|=\sum^{m}_{i=0}I_{i}$ offspring particles, among which the first $I_{0}$ have mark 0, $I_{1}$ have mark 1, and so on, so that each particle has mark $i$ for $i=0,\ldots,m.$ 2. 2.

For generation $n\geq 1,$ let the label for a particle be given as $k=(k_{1},\ldots,k_{n-1},k_{n})\in{\mathbb{N}}^{n}.$ Furthermore, denote by $k^{-}:=(k_{1},\ldots,k_{n-2},k_{n-1})$ the parent particle of $k$. The particle $k$ starts from $X^{k^{-}}_{T_{k^{-}}}$ at time $T_{k^{-}}$:

• •

The position $X^{k}$ of the particle during its lifetime is given by

[TABLE]

• •

The arrival time of the particle is given by

[TABLE]

• •

At the arrival time, if $X^{k}_{T_{k}}\notin{\mathcal{O}}$, then the particle $k$ dies without giving rise to any offspring, else the particle $k$ dies and gives rise to $|I^{k}|$ offsprings which belong to the $(n+1)$th generation, and are indexed by label $(k_{1},\ldots,k_{n-1},k_{n},i)$ for $i=0,\ldots,|I^{k}|-1.$

• •

Given $I^{k}=(I^{k}_{0},I^{k}_{1},\ldots,I^{k}_{m}),$ we have $|I^{k}|=\sum^{m}_{i=0}I^{k}_{i}$ offspring particles, among which the first $I^{k}_{0}$ have mark 0, $I^{k}_{1}$ have mark 1, and so on, so that each particle has mark $i$ for $i=0,\ldots,m.$

In addition, we denote by ${\mathcal{K}}^{x}_{n}$ the collection of particles of the $n$th generation and by ${\mathcal{K}}^{x}=\bigcup_{n\in\mathbb{N}}{{\mathcal{K}}^{x}_{n}}$ the collection of all particles. Finally, we introduce the filtration

[TABLE]

We make the following assumptions to ensure that the branching diffusion process above is well-defined and for further developments.

Assumption 2.2**.**

(i)* The probability distribution $(p_{l})_{l\in L}$ satisfies $p_{l}>0$ and $\sum_{l\in L}{|l|p_{l}}<\infty$.

(ii) The probability density $\rho$ is strictly positive.

(iii) The coefficients $(\mu,\sigma)$ are Lipschitz on $\bar{{\mathcal{O}}}$.*

Part (i) of the assumption above ensures that the number of particles remains finite in finite time in the underlying branching process, see, e.g., Athreya and Ney [1, Theorem 4.1.1]. In addition, under assertion (iii) , there exists a unique solution (up to the boundary) to the stochastic differential equation (SDE) corresponding to $(\mu,\sigma)$ and so the branching diffusion process is well-defined.

Assumption 2.3**.**

The branching diffusion process goes extinct almost surely.

It is clearly sufficient to assume that $\sum_{l\in L}{|l|p_{l}}\leq 1$ for Assumption 2.3 to hold. However, since particles are absorbed at the boundary, we can derive much weaker conditions, especially if the domain is small. For instance, in the case of branching Brownian motion with exponential lifetime of parameter $\beta$, Assumption 2.3 is equivalent to

[TABLE]

where $\lambda_{1}$ is the first positive eigenvalue of the Laplacian operator in the domain ${\mathcal{O}}$, see Sevast′yanov [33] or Watanabe [36]. Also see Remark 3.3 below for further developments.

3   Semi-Linear PDEs with Linear Gradient Term

In this section, we are concerned with semi-linear PDEs with polynomial non-linearity involving the unknown function but not its derivatives, i.e., we assume that $m=0$ in Section 2 so that all the particles have the same mark [math] and PDE (5) reads as

[TABLE]

where $f(x,y)=\sum_{l\in L}{c_{l}(x)y^{l}}$, $L\subset\mathbb{N}$. Throughout this section, we suppose that Assumption 2.1–2.3 remain valid. Additionally, we set $\rho$ as the density of exponential distribution with parameter $\beta$.

3.1   Probabilistic Representation

We consider a branching diffusion process starting from $x\in{\mathcal{O}}$ as in Section 2.2 and we introduce the following random variable:

[TABLE]

Proposition 3.1**.**

Suppose that PDE (7) has a solution $u\in{\mathcal{C}}^{2}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ such that the sequence $(\psi^{x}_{n})_{n\in\mathbb{N}}$ given by

[TABLE]

is uniformly integrable, then we have $u(x)={\mathbb{E}}[\psi^{x}]$.

Proof.

By applying Itô’s formula to $(e^{-\beta t}u(X^{x}_{t}))_{t\geq 0}$, we obtain the following Feynman-Kac representation:

[TABLE]

See, e.g., Freidlin [13, Theorem 2.2.1] for a detailed proof. Then, we can write

[TABLE]

Since an empty product is equal to $1$ by convention, it follows that $u(x)={\mathbb{E}}[\psi^{x}_{0}]$. Next, since each offspring has the same dynamic as the parent particle, we can repeat the above calculations to write for $k\in{\mathcal{K}}^{x}_{1},$

[TABLE]

We use the result in (10) and plug it back in (8) to obtain $u(x)={\mathbb{E}}[\psi^{x}_{1}]$ by conditional independence of particles in ${\mathcal{K}}^{x}_{1}$ given ${\mathcal{F}}_{0}$. Similarly, we can show by iteration that for any $n\in\mathbb{N},$ we have $u(x)={\mathbb{E}}[\psi^{x}_{n}]$. To conclude, it remains to observe that $\psi_{n}^{x}$ converges to $\psi^{x}$ almost surely in view of Assumption 2.3. Thus, if we suppose that $(\psi^{x}_{n})_{n\in\mathbb{N}}$ is uniformly integrable, as $n\rightarrow\infty,$ we get $u(x)={\mathbb{E}}[\psi^{x}]$. ∎

Proposition 3.1 shows that there is at most one solution to PDE (7) satisfying an appropriate integrability condition and such a solution admits a probabilistic representation using branching diffusion processes. It is not completely satisfactory since one needs to prove first the existence of a classical solution and even then the uniform integrability condition is hard to derive as illustrated by the example of Section 5.1. In order to overcome these limitations, we next establish a result which shows directly that the probabilistic representation provides a (viscosity) solution of PDE (7). This is the main result of this section. It is stated under abstract assumptions for which we will provide explicit sufficient conditions in the subsequent sections.

Theorem 3.1**.**

Suppose $u:x\mapsto{\mathbb{E}}[\psi^{x}]$ is well-defined and continuous on $\bar{{\mathcal{O}}}$. Then $u$ is a viscosity solution of PDE (7).

Proof.

We first observe that by definition of $u$, it holds

[TABLE]

where

[TABLE]

Furthermore, the branching property, which says that each offspring starts the same but an independent branching diffusion as the parent particle, yields that conditioned on $X^{x}_{\tau}$ and $I$, $(\psi^{X^{x}_{\tau}}_{0},\ldots,\psi^{X^{x}_{\tau}}_{I-1})$ are independent random variables identical in law to $\psi^{X_{0}}$ where $X_{0}$ is distributed as $X^{x}_{\tau}$ and independent of ${\mathcal{F}}_{n}$ for all $n\in\mathbb{N}$. Through this argument, we deduce that

[TABLE]

Working backward along the lines of the proof of Proposition 3.1, we deduce that

[TABLE]

Next, for any $\delta>0,$ it follows from Markov property that

[TABLE]

The fact that $u$ is a viscosity solution of PDE (7) now follows from classical arguments. For the sake of completeness, let us show that $u$ is a viscosity subsolution. Consider $\varphi\in{\mathcal{C}}^{2}_{b}({\mathcal{O}})$ such that $x\in{\mathcal{O}}$ is a maximum point of $u-\varphi$ and $u(x)=\varphi(x)$. First, we observe by applying Itô’s formula that

[TABLE]

Since $u(x)=\varphi(x)$ and $u\leq\varphi$ otherwise, we deduce by using (12) that

[TABLE]

Since $y\mapsto f(y,u(y))$ is continuous and bounded, it follows from the mean value theorem and the dominated convergence theorem that

[TABLE]

Thus $u$ is a viscosity subsolution of PDE (7). The fact that $u$ is a viscosity supersolution results from similar arguments.

∎

Remark 3.1**.**

It will turn out from the study of the next sections that the explicit conditions we provide to ensure that $x\mapsto{\mathbb{E}}[\psi^{x}]$ is a solution to PDE (7) also entail that $x\mapsto{\mathbb{E}}[|\psi^{x}|]$ is a solution to

[TABLE]

In particular, we observe that the monotonicity assumption is not relevant in our setting. The limiting assumption in our approach is the integrability condition on $(\psi^{x})_{x\in{\mathcal{O}}}$ which ensures that $u:x\mapsto{\mathbb{E}}[\psi^{x}]$ is well-defined. See Section 3.3 for more details.

Remark 3.2**.**

We can also work with a general lifetime distribution $\rho$ and extend the arguments above, as done in Section 4, to derive another probabilistic representation which is given as

[TABLE]

where $\Delta T_{k}:=T_{k}-T_{k^{-}}$ is the lifetime of particle $k$ and $\bar{F}(t):=\int_{t}^{\infty}{\rho(s)\,ds}$, $t\geq 0$. However, this leads to more stringent assumptions when studying the integrability condition as in Section 3.3.

3.2   Continuity Assumption

Let us give explicit sufficient conditions for the continuity assumption in Theorem 3.1 to hold. We use two different approaches leading to slightly different conditions. Both approaches rely on appropriate integrability conditions on $(\psi^{x})_{x\in{\mathcal{O}}}$ which we discuss in the next section.

3.2.1   PDE Approach

Assumption 3.1**.**

(i)* The diffusion coefficient $\sigma$ is uniformly elliptic.222There exists $\lambda>0$ such that $\sigma\sigma^{*}(x)\geq\lambda I_{d}$ for all $x\in\bar{{\mathcal{O}}}$.

*(ii)The boundary $\partial{\mathcal{O}}$ is of class ${\mathcal{C}}^{1,\alpha}$.

Proposition 3.2**.**

Suppose Assumption 3.1 holds and $(\psi^{x})_{x\in{\mathcal{O}}}$ is uniformly bounded in $L^{1}$, then the map $u:x\mapsto{\mathbb{E}}[\psi^{x}]$ is continuous.

The proof follows immediately from (11) and the following lemma.

Lemma 3.1**.**

*Suppose Assumption 3.1 holds, then the following statements are satisfied:

(i) The map $\bar{{\mathcal{O}}}\ni x\mapsto{\mathbb{E}}[e^{-\beta{\eta^{x}}}h(X^{x}_{\eta^{x}})]$ is continuous.

(ii) For any $g:{\mathcal{O}}\rightarrow\mathbb{R}$ bounded measurable, the map $\bar{{\mathcal{O}}}\ni x\mapsto{\mathbb{E}}[\int_{0}^{\eta^{x}}{e^{-\beta s}g(X^{x}_{s})\,ds}]$ is continuous.*

Proof.

(i) Let us first study the continuity of $x\mapsto{\mathbb{E}}[e^{-\beta{\eta^{x}}}h(X^{x}_{\eta^{x}})]$. Under the assumptions above, it is well-known that there exists a smooth solution $\varphi\in{\mathcal{C}}^{2}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ to the following PDE:

[TABLE]

See, e.g., Gilbarg and Trudinger [15, Theorem 6.13]. By Itô’s formula, we deduce that $\varphi(x)={\mathbb{E}}[e^{-\beta{\eta^{x}}}$ $h(X^{x}_{\eta^{x}})]$ and the conclusion follows.

(ii)Let us now turn to the continuity of $x\mapsto{\mathbb{E}}[\int_{0}^{\eta^{x}}{e^{-\beta s}g(X^{x}_{s})\,ds}]$. We start by showing that $x\mapsto{\mathbb{E}}[g(X^{x}_{s})\mathbf{1}_{s<{\eta^{x}}}]$ is continuous for $s>0$. If $g$ is continuous and $g=0$ on $\partial{\mathcal{O}}$, it is known that the unique smooth solution $\chi\in{\mathcal{C}}^{1,2}((0,T]\times{\mathcal{O}})\cap{\mathcal{C}}([0,T]\times\bar{{\mathcal{O}}})$ of

[TABLE]

is of the form

[TABLE]

where $G$ is the so-called Green function of PDE (3.2.1) (see, e.g., Ladyženskaja et al. [23, Theorem 4.16.2]). Then, it follows from Itô’s formula that

[TABLE]

In particular, $y\mapsto G(s,x;0,y)$ appears as the density of $X^{x}_{s}$ on the event $\{s<{\eta^{x}}\}$, and so the identity above remains valid for any $g$ bounded measurable. Furthermore, $x\mapsto G(s,x;0,y)$ is continuous and satisfies

[TABLE]

for some constant $C>0$ depending on $T$, see [23, Equation 4.16.16]. Hence the desired result follows from the dominated convergence theorem. To conclude, it remains to observe that

[TABLE]

and to apply once again the dominated convergence theorem.

∎

3.2.2   Probabilistic Approach

Assumption 3.2**.**

*(i)**The diffusion coefficient $\sigma$ is uniformly elliptic on $\partial{\mathcal{O}}$.333There exists $\lambda>0$ such that $\sigma\sigma^{*}(x)\geq\lambda I_{d}$ for all $x\in\partial{\mathcal{O}}$.

*(ii)*The boundary $\partial{\mathcal{O}}$ satisfies an exterior cone condition.444See, e.g., Gilbarg and Trudinger [15, Problem 2.12] for a definition.

*(iii)The stopping time $\eta^{x}$ is finite almost surely.

Sufficient conditions for Part (iii) to hold are provided in Freidlin [13, Lemma 3.3.1]. For instance, it suffices to assume that there exists $1\leq i\leq d$ such that $\sum_{j=1}^{d}\sigma_{ij}^{2}(x)>0$ or $|\mu_{i}(x)|>0$ for all $x\in{\mathcal{O}}$.

Proposition 3.3**.**

Suppose Assumption 3.2 holds and $(\psi^{x})_{x\in{\mathcal{O}}}$ is uniformly integrable, then $u:x\mapsto{\mathbb{E}}[\psi^{x}]$ is continuous.

Proof.

Clearly it suffices to prove that $\bar{{\mathcal{O}}}\ni x\mapsto\psi^{x}$ is almost surely continuous, in the sense that for all sequence $(x_{n})_{n\in\mathbb{N}}$ converging to $x$, it holds

[TABLE]

This essentially follows from the almost sure continuity of $x\mapsto\eta^{x}$ stated in Lemma 3.2 below.

(i) Let us first show that the contribution of the first particle to $\psi^{x}$ is almost surely continuous, i.e.,

[TABLE]

To achieve this, let us consider the set

[TABLE]

In view of Lemma 3.2 below, it is clear that ${\mathbb{P}}(\Omega^{\emptyset})=1$. Furthermore, one easily checks that, for every $\omega\in\Omega^{\emptyset}$, both (15)–(16) hold.

(ii) Let us show next that the contribution of the particles of the first generation to $\psi^{x}$ is almost surely continuous. For every $i\in\mathbb{N}$, we denote by $(X^{i,x}_{s})_{s\geq 0}$ the unique solution of

[TABLE]

where the Brownian motion $B^{i}$ is defined by

[TABLE]

Clearly, $(X^{i,x}_{s})_{s\geq 0}$ coincides with the trajectory $X^{i}$ of particle $i$ during its lifetime. Let $\eta^{i,x}$ be the first exit time of $X^{i,x}$ from ${\mathcal{O}}$ and consider

[TABLE]

where $\bar{\tau}^{i}:=\tau+\tau^{i}$. Once again, ${\mathbb{P}}(\Omega^{i})=1$ and

[TABLE]

Thus, we conclude that on the set $\Omega^{\emptyset}\cap\bigcap_{i\in\mathbb{N}}{\Omega^{i}}$ of probability $1$, the contributions of the particles of generations [math] and $1$ to $\psi^{x}$ is almost surely continuous.

(iii) The desired result follows by iteration. We denote for all $n\geq 2$, $k\in\mathbb{N}^{n}$,

[TABLE]

where $\bar{\tau}^{k}:=\bar{\tau}^{k^{-}}+\tau$ and $(X^{k,x}_{s})_{s\geq 0}$ is the unique solution of

[TABLE]

with

[TABLE]

Then on the set $\bigcap_{k\in{\mathbb{K}}}\Omega^{k}$ of probability $1$, it holds that $\lim_{n\rightarrow\infty}\psi^{x_{n}}=\psi^{x}$. ∎

Lemma 3.2**.**

Suppose Assumption 3.2 holds, then the map $\bar{{\mathcal{O}}}\ni x\mapsto\eta^{x}$ is almost surely continuous, in the sense that for all sequence $(x_{n})_{n\in\mathbb{N}}$ converging to $x$, it holds

[TABLE]

Proof.

The proof is due to Darling and Pardoux [10], see also [28, Proposition 4.4]. More precisely, the authors show that if the stopping time

[TABLE]

is finite almost surely and ${\mathbb{P}}\left(\tau^{x}>0\right)=0$ for all $x\in\partial{\mathcal{O}}$, then $x\mapsto\tau^{x}$ is almost surely continuous. It remains to observe that, under Assumption 3.2 (i)–(ii), ${\mathbb{P}}(\tau^{x}=0)=1$ for all $x\in\partial{\mathcal{O}}$ (see Bass [3, Corollary 3.3.2] or Pinsky [31, Theorem 2.3.3]) and thus $\tau^{x}=\eta^{x}$ for all $x\in\bar{{\mathcal{O}}}$. ∎

3.3   Integrability Condition

In this section, we provide explicit sufficient conditions to verify the integrability conditions on $(\psi_{x})_{x\in{\mathcal{O}}}$ required for Propositions 3.2 and 3.3. Actually we study boundedness of $(\psi_{x})_{x\in{\mathcal{O}}}$ in $L^{q}$ for $q\geq 1$. In particular, the case $q=2$ ensures that the corresponding Monte Carlo estimator has finite variance. Let us introduce the constant

[TABLE]

Clearly, it holds $|\psi^{x}|\leq C_{0}^{|{\mathcal{K}}^{x}|}$ where $|{\mathcal{K}}^{x}|$ denotes the cardinality of the set ${\mathcal{K}}^{x}$, i.e., the total number of particles. In particular, if $C_{0}\leq 1$, then $|\psi^{x}|\leq 1$. To the best of our knowledge, this is the only condition that has been used so far in the literature to ensure the integrability of $\psi^{x}$.

In the rest of this section, we provide more technical but weaker conditions. First, we establish the desired result under minimal assumptions.

Proposition 3.4**.**

Suppose that there exists a non-negative function $v\in{\mathcal{C}}^{2}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ satisfying

[TABLE]

then we have ${\mathbb{E}}[|\psi^{x}|^{q}]\leq v(x)$. In particular, $(\psi^{x})_{x\in{\mathcal{O}}}$ is uniformly bounded in $L^{q}$.

Proof.

First, we observe that Itô’s formula yields

[TABLE]

Next, by repeating the arguments of Proposition 3.1, we get

[TABLE]

Then the conclusion follows easily from Fatou’s lemma. ∎

In practice, one can look for a constant supersolution in order to apply Proposition 3.4. For instance, for the case $q=1$, such a solution exists if and only if $\sum_{l\in L}{\left|c_{l}\right|}\|h\|_{\infty}^{l}\leq\|h\|_{\infty}.$ Otherwise, finding a supersolution has to be done on a case by case basis and might turn out to be a difficult task. For this reason, we provide an alternative result for which conditions are easier to verify than Proposition 3.4 and less stringent than assuming $C_{0}\leq 1$. Essentially, it signifies that the integrability condition is satisfied if the domain is sufficiently small.

The idea is to dispose of the spatial parameter $x$ by introducing a branching process which stochastically dominates the branching diffusion process. Denote

[TABLE]

Let us introduce a new probability mass function $(\tilde{p}_{l})_{l\in\tilde{L}}$ with $\tilde{L}:=L\cup\{0\}$ as follows:

[TABLE]

and the corresponding transition matrix $\tilde{P}=(\tilde{P}_{i,j})_{i,j\geq 0}$ given by

[TABLE]

Proposition 3.5**.**

Denote by $R$ the common radius of convergence of $f(s):=\sum_{l\in L}{p_{l}s^{l}}$ and $\tilde{f}(s):=\sum_{l\in\tilde{L}}{\tilde{p}_{l}s^{l}}$. If $R>1$ and $\sum_{l\in\tilde{L}}{l\tilde{p}_{l}}<1$, then $(\psi^{x})_{x\in{\mathcal{O}}}$ is uniformly bounded in $L^{q}$ provided that $C^{q}_{0}\leq\gamma$ where

[TABLE]

where $s^{*}$ is the solution (if any) to $s\tilde{f}^{\prime}(s)=\tilde{f}(s)$ and $s^{*}=R$ otherwise. In addition, if $R=\infty$, then $\gamma$ goes to infinity as $\mathrm{diam}({\mathcal{O}})$ goes to zero.

Proof.

(i) Denote by $(N^{x}_{n})_{n\in\mathbb{N}}$ and $(\tilde{N}_{n})_{n\in\mathbb{N}}$ the number of particles at the $n$th jump time (arrival time) in the branching diffusion process and in a branching process with offspring distribution $(\tilde{p}_{l})_{l\in\tilde{L}}$ respectively. In particular, $(\tilde{N}_{n})_{n\in\mathbb{N}}$ is a Markov chain with transition matrix $\tilde{P}$. Let us first show that $(N^{x}_{n})_{n\in\mathbb{N}}$ is stochastically dominated by $(\tilde{N}_{n})$, i.e., ${\mathbb{P}}(N^{x}_{n}\geq l)\leq{\mathbb{P}}(\tilde{N}_{n}\geq l)$ for all $n,l\geq 1$. We observe that

[TABLE]

It follows that ${\mathbb{P}}(N^{x}_{1}\geq l)\leq{\mathbb{P}}(\tilde{N}_{1}\geq l)$ for all $l\geq 1$. For the incremental step, we first observe that for all $l\geq 1$,

[TABLE]

where the second inequality follows from Markov property. In addition, it holds

[TABLE]

where the second inequality follows once again from Markov property. Combining both computations above, we deduce that

[TABLE]

Hence the desired result follows easily by induction and Markov property.

(ii) Let us show next that ${\mathbb{E}}[|\psi^{x}|]\leq{\mathbb{E}}[C_{0}^{\zeta}]$ where $\zeta$ the extinction time of $(\tilde{N}_{n})_{n\in\mathbb{N}}$. First we recall that $|\psi^{x}|\leq C_{0}^{|{\mathcal{K}}^{x}|}$ where $|{\mathcal{K}}^{x}|$ denotes the cardinality of the set ${\mathcal{K}}^{x}$. We observe further that $|{\mathcal{K}}^{x}|$ coincides with the extinction time of $(N^{x}_{n})_{n\in\mathbb{N}}$. Clearly, Step (i) above yields that ${\mathbb{P}}(|{\mathcal{K}}^{x}|\geq n)\leq{\mathbb{P}}(\zeta\geq n)$ for all $n\geq 1$, and thus

[TABLE]

(iii) We are now in a position to conclude the proof. It follows from Daley [9, Theorem 2] that the power series ${\mathbb{E}}[s^{\zeta}]$ converges on its radius of convergence $\gamma=\frac{s^{*}}{\tilde{f}(s^{*})}$. In addition, if $R=\infty$, then there exists a solution $s^{*}$ to $s\tilde{f}^{\prime}(s)=\tilde{f}(s)$. Using further $\tilde{f}(s)=1-\delta+\delta f(s)$, we deduce that

[TABLE]

It follows that $s^{*}$ goes to infinity as $\delta$ goes to zero, or equivalently, $\mathrm{diam}({\mathcal{O}})$ goes to zero. Then we have

[TABLE]

To conclude, it remains to observe that $\frac{f(s^{*})}{s^{*}f^{\prime}(s^{*})}$ is bounded away from $1$ since $s\mapsto\frac{f(s)}{sf^{\prime}(s)}$ is decreasing.

∎

Remark 3.3**.**

It follows immediately from Step (i) of the proof of Proposition 3.5 that if $\sum_{l\in\tilde{L}}{l\tilde{p}_{l}}\leq 1$, then the branching diffusion process goes extinct almost surely, i.e., Assumption 2.3 holds.

4   Semi-Linear PDEs with Non-Linear Gradient Term

In this section we study the case of semi-linear PDE with non-linearity in gradient of the solution, i.e., we assume that $m\geq 1$ in Section 2 so that the particles in the branching diffusion process carry different marks to account for it. Throughout this section, we suppose that Assumption 2.1–2.3 remain valid.

4.1   Probabilistic Representation

Our next assumption is the key automatic differentiation condition on the underlying diffusion process $X^{x}$. We will provide explicit conditions and formulas for it in the next sections.

Assumption 4.1**.**

(i)* The map $x\mapsto{\mathbb{E}}[e^{-\beta{\eta^{x}}}h(X^{x}_{\eta^{x}})]$ belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ and there exists a measurable function ${\mathcal{W}}_{\partial{\mathcal{O}}}(x,W)={\mathcal{W}}_{\partial{\mathcal{O}}}(x,(W_{r})_{r\in[0,{\eta^{x}}]})$ such that*

[TABLE]

(ii)* For any $g:{\mathcal{O}}\rightarrow\mathbb{R}$ bounded measurable, the map $x\mapsto{\mathbb{E}}[\int_{0}^{{\eta^{x}}}{e^{-\beta s}g(X^{x}_{s})\,ds}]$ belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ and there exists a measurable function ${\mathcal{W}}_{{\mathcal{O}}}(s,x,W)={\mathcal{W}}_{{\mathcal{O}}}(s,x,(W_{r})_{r\in[0,s]})$ such that*

[TABLE]

Let us define ${\mathcal{W}}(s,x,W)={\mathcal{W}}(s,x,(W_{r})_{r\in[0,s]})$ as follows:

[TABLE]

We consider a marked branching diffusion process starting from $x\in{\mathcal{O}}$ as in Section 2.2 and denote

[TABLE]

where $m_{k}$ and $\Delta T_{k}:=T_{k}-T_{k^{-}}$ stand for the mark and the lifetime of particle $k$ respectively. We next introduce the following random variable:

[TABLE]

where $\bar{F}(t):=\int_{t}^{\infty}{\rho(s)\,ds}$, $t\geq 0$.

Proposition 4.1**.**

Suppose Assumption 4.1 holds. Assume further that PDE (5) has a solution $u\in{\mathcal{C}}^{2}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ such that the functions $(b_{i}\cdot Du)_{i=1,\cdots,m}$ are bounded and the sequence $(\psi^{x}_{n})_{n\in\mathbb{N}}$ defined by

[TABLE]

is uniformly integrable. Then, we have $u(x)={\mathbb{E}}[\psi^{x}]$.

Proof.

The proof follows similar arguments as the proof of Proposition 3.1. Using Itô’s formula, we have the following Feynman-Kac representation:

[TABLE]

Then, we can write

[TABLE]

In other words, we have $u(x)={\mathbb{E}}\left[\psi_{0}^{x}\right]$. Furthermore, in the original Feynman-Kac formula, we obtain by differentiating and using Assumption 4.1,

[TABLE]

Since each offspring has the same dynamic as the parent particle, we can repeat the above calculations for $k\in{\mathcal{K}}^{x}_{1}$ and plug the results back in (20) to obtain $u(x)={\mathbb{E}}[\psi_{1}^{x}]$ by conditional independence of particles in ${\mathcal{K}}^{x}_{1}$ given ${\mathcal{F}}_{0}$. We conclude by iteration that, for any $n\in\mathbb{N}$, $u(x)={\mathbb{E}}[\psi_{n}^{x}]$ and, as $n\rightarrow\infty,$ $u(x)={\mathbb{E}}[\psi^{x}]$. ∎

Similar to Section 3.1, Proposition 4.1 provides a result of uniqueness for a class of semi-linear PDEs given an appropriate uniform integrability condition is satisfied. Next we establish a result of existence by showing that the probabilistic representation is a viscosity solution of PDE (5). This is the main result of this section.

Theorem 4.1**.**

Suppose Assumption 4.1 holds. If we further assume that $(\psi^{x})_{x\in{\mathcal{O}}}$ and, for $i=1,\ldots,m$, $(\psi^{x}b_{i}(x)\cdot{\mathcal{W}}(T_{\emptyset},x,W))_{x\in{\mathcal{O}}}$ are uniformly bounded in $L^{1}$, then $u:x\mapsto{\mathbb{E}}[\psi^{x}]$ belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ and solves PDE (5) in the viscosity sense.

Proof.

The proof follows similar arguments as the proof of Theorem 3.1. We observe first that, by definition of $u$, it holds

[TABLE]

where

[TABLE]

It further follows from the branching property that, conditioned on ${\mathcal{F}}_{0}$, $(\psi^{X^{x}_{\tau}}_{i})_{i=0,\ldots,\lvert I\rvert-1}$ are independent random variables, among which the first $I_{0}$ are identical in law to $\psi^{X_{0}}$, the next $I_{1}$ are identical in law to $\psi^{X_{0}}b_{1}(x)\cdot{\mathcal{W}}(T_{\emptyset},x,W)$ and so on, where $X_{0}$ is distributed as $X^{x}_{\tau}$ and independent of ${\mathcal{F}}_{n}$ for all $n\in\mathbb{N}$. Through this argument, we deduce that

[TABLE]

where $(v_{i})_{i=1,\ldots,m}$, $v_{i}:{\mathcal{O}}\mapsto\mathbb{R}$ are defined as

[TABLE]

Working backward along the lines of the proof of Proposition 4.1, we deduce that

[TABLE]

In particular, since $u$ and $(v_{i})_{i=1,\ldots,m}$ are bounded by assumption, it follows from Assumption 4.1 that $u$ belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ and

[TABLE]

Thus, for all $i=1,\ldots,m$, $v_{i}$ coincides with $b_{i}\cdot Du$ and (24) reads as

[TABLE]

The fact that $u$ is a viscosity solution of PDE (5) now follows by classical arguments. ∎

4.2   Automatic Differentiation Formula: the General Case

The aim of this section is to provide sufficient conditions to ensure that Assumption 4.1 holds and to derive explicit formula for ${\mathcal{W}}$. The automatic differentiation formula discussed in the following originates from Thalmaier [35] and was subsequently developed by Delarue [11] and Gobet [16].

Assumption 4.2**.**

(i)* The coefficients $(\mu,\sigma)$ belong to ${\mathcal{C}}^{1,\alpha}(\bar{{\mathcal{O}}})$.

(ii) The diffusion coefficient $\sigma$ is uniformly elliptic.

(iii) The boundary $\partial{\mathcal{O}}$ is of class ${\mathcal{C}}^{2}$.

(iv) The function $h$ can be extended to a function of class ${\mathcal{C}}^{1,\alpha}$ on $\bar{{\mathcal{O}}}$.*

We start by establishing a technical lemma. Fix a finite horizon $T>0$ and denote for any $s>0$,

[TABLE]

Lemma 4.1**.**

Under Assumption 4.2, it holds for all $x\in{\mathcal{O}}$ and $s>0$,

[TABLE]

In addition, if we denote

[TABLE]

then there exists $t<s$ such that $\zeta_{s}\leq{\eta^{x}}\wedge t$ and for all $q\geq 1$,

[TABLE]

where $C>0$ depends on $q$ and $T$ but not on $x$ or $s$.

Proof.

The proof essentially follows from Delarue [11]. Indeed, if ${\mathcal{O}}$ is a ball, both identities (25) and (26) are easily obtained by repeating the arguments of Propositions 2.3 and 2.4 in [11] while working with $(r,y)\mapsto d(y,\partial{\mathcal{O}})\sqrt{s-r}$ instead of $(r,y)\mapsto d(y,\partial{\mathcal{O}})(s-r)$. For an arbitrary domain, it suffices to work with a ${\mathcal{C}}^{2}$–extension of the distance to the boundary (see, e.g., Gilbarg and Trudinger [15, Lemma 14.16]). Additionally, it follows from (25) that $\zeta_{s}\leq{\eta^{x}}\wedge s$. Furthermore, it holds for all $t<s$,

[TABLE]

where $C:=\mathrm{diam}({\mathcal{O}})^{2}/4$. Thus for $t=(1-e^{-C})s$, we have $\zeta_{s}\leq{\eta^{x}}\wedge t$. ∎

Proposition 4.2**.**

Under Assumption 4.2, the assertions of Assumption 4.1 are satisfied with

[TABLE]

where $Y^{x}$ is the tangent process given by

[TABLE]

and $\sigma_{i}$ denotes the $i$th column of $\sigma$. In addition, it holds for all $q\geq 1$,

[TABLE]

This result is a slight extension of the automatic differentiation formula obtained by Delarue [11] and Gobet [16]. We postpone the proof to Appendix A. Note that estimate (27) is an easy consequence of Lemma 4.1. Indeed, successively using the Burkhölder-Davis-Gundy inequality, the boundedness of $\sigma^{-1}$ and the fact that $\sup_{0\leq r\leq T}\{|Y^{x}_{r}|\}$ has finite moments, we derive for all $s\leq T$,

[TABLE]

Remark 4.1**.**

The choice of weight function in the automatic differentiation formula is not unique. Actually we can replace $(\theta_{s}(r,X^{x}_{r}))_{r\geq 0}$ by any predictable process $(\tilde{\theta}_{s}(r))_{r\geq 0}$ satisfying appropriate integrability conditions and such that

[TABLE]

We refer to Delarue [11] for more examples.

We now provide sufficient conditions to ensure that $(\psi^{x})_{x\in{\mathcal{O}}}$ and $(\psi^{x}b_{i}(x)\cdot{\mathcal{W}}(T_{\emptyset},x,W))_{x\in{\mathcal{O}}}$, $i=1,\ldots,m$, are uniformly bounded in $L^{q}$ for $q\geq 1$. Let us define

[TABLE]

In view of (27), $C_{1,q}$ and $C_{2,q}$ are clearly finite if the following condition is satisfied:

[TABLE]

Unlike in the case of linear gradient term, we need to chose a lifetime distribution $\rho(t)\neq\beta e^{-\beta t}$ to ensure that the first two terms in the above expression are finite. For instance, we can take a Gamma distribution with shape parameter $0.5$ and rate parameter $\beta^{\prime}<\beta$, i.e., $\rho(t)=\sqrt{\beta^{\prime}/(\pi t)}e^{-\beta^{\prime}t}$. Another choice consists in taking a generalized gamma distribution of the form $\rho(t)=\beta^{\prime}/(2\sqrt{t})e^{-\beta^{\prime}\sqrt{t}}$ for any $\beta^{\prime}>0$.

Let us denote

[TABLE]

Similar to Section 3.3, we then introduce a new probability mass function $(\tilde{p}_{l})_{l\in\tilde{L}}$ with $\tilde{L}:=L\cup\{0_{m+1}\}$ where $0_{m+1}$ is the zero element of $\mathbb{N}^{m+1}$ as follows:

[TABLE]

and the corresponding transition matrix $\tilde{P}=(\tilde{P}_{i,j})_{i,j\geq 0}$ given by

[TABLE]

Proposition 4.3**.**

*Denote $C_{0,q}:=C_{1,q}\vee C_{2,q}$.

(i) If $C_{0,q}\leq 1$, then $(\psi^{x})_{x\in{\mathcal{O}}}$ and $(\psi^{x}b_{i}(x)\cdot{\mathcal{W}}(T_{\emptyset},x,W))_{x\in{\mathcal{O}}}$, $i=1,\ldots,m$, are bounded by $1$ in $L^{q}$.

(ii) Denote by $R$ the common radius of convergence of $f(s):=\sum_{l\in L}{p_{l}s^{|l|}}$ and $\tilde{f}(s):=\sum_{l\in\tilde{L}}{\tilde{p}_{l}s^{|l|}}$. If $R>1$ and $\sum_{l\in\tilde{L}}{|l|\tilde{p}_{l}}<1$, then $(\psi^{x})_{x\in{\mathcal{O}}}$ and $(\psi^{x}b_{i}(x)\cdot{\mathcal{W}}(T_{\emptyset},x,W))_{x\in{\mathcal{O}}}$, $i=1,\ldots m$, are uniformly bounded in $L^{q}$ provided that $C^{q}_{0,q}\leq\gamma$ where*

[TABLE]

where $s^{*}$ is the solution (if any) to $s\tilde{f}^{\prime}(s)=\tilde{f}(s)$ and $s^{*}=R$ otherwise. In addition, if $R=+\infty$, then $\gamma$ goes to infinity as $\mathrm{diam}({\mathcal{O}})$ goes to zero.

Proof.

Let us denote by $\mathcal{G}$ the following $\sigma$-algebra:

[TABLE]

Conditioning by $\mathcal{G}$, we obtain for all $x\in{\mathcal{O}}$,

[TABLE]

and the same inequality holds if we consider $\psi^{x}b_{i}(x)\cdot{\mathcal{W}}(T_{\emptyset},x,W)$ instead of $\psi^{x}$. Part (i) of Proposition 4.3 follows immediately. As for Part (ii), the desired result follows by the same arguments as in Proposition 3.5. ∎

4.3   Automatic Differentiation Formula: the One-Dimensional Case

The automatic differentiation formula given by Proposition 4.2, though quite general, has several drawbacks for numerical applications. Indeed one needs to compute a stochastic integral for each particle that holds a non-zero mark. In addition, the weight function ${\mathcal{W}}(s,x,W)$ explodes as $s$ goes to [math] or $x$ approaches the boundary. In this section, we provide a simpler automatic differentiation formula that is satisfied for a one-dimensional Brownian motion exiting from an interval.

Let us assume throughout this section that $d=1$, $b=0$, $\sigma=1$ and ${\mathcal{O}}=(-r,r)$ for some $r>0$. We denote

[TABLE]

We also introduce the function $\mathfrak{W}:{\mathcal{O}}\times\bar{{\mathcal{O}}}\rightarrow\mathbb{R}$ as

[TABLE]

Both lemmas below show that Assumption 4.1 is satisfied in this setting.

Lemma 4.2**.**

For any $g:{\mathcal{O}}\rightarrow\mathbb{R}$ bounded measurable, the map $x\mapsto{\mathbb{E}}[\int_{0}^{{\eta^{x}}}{e^{-\beta s}g(W^{x}_{s})\,ds}]$ belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ and

[TABLE]

Proof.

Let us denote

[TABLE]

We first assume that $g$ is continuous. Then $u$ satisfies the following ODE:

[TABLE]

By standard arguments, including variation of parameters, we deduce that

[TABLE]

where $G:\bar{{\mathcal{O}}}\times{\mathcal{O}}\rightarrow\mathbb{R}$ is given by

[TABLE]

By a monotone class argument, (29) remains valid if $g$ is only assumed bounded measurable. Then the continuity of $\chi$ follows immediately from the dominated convergence theorem. Further, a direct calculation yields that

[TABLE]

The desired result then follows by differentiation under the integral sign in (29). ∎

Lemma 4.3**.**

The map $x\mapsto{\mathbb{E}}[e^{-\beta{\eta^{x}}}h(W^{x}_{{\eta^{x}}})]$ belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ and

[TABLE]

Proof.

Let us denote

[TABLE]

It satisfies the following ODE:

[TABLE]

It follows that

[TABLE]

where $H:{\mathcal{O}}\times\partial{\mathcal{O}}\rightarrow\mathbb{R}$ is given by

[TABLE]

In particular, $\varphi$ clearly belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$. To conclude, it remains to observe that

[TABLE]

∎

Remark 4.2**.**

We observe that the weight function still explodes when $x$ approaches the boundary as in Section 4.2. In view of the proof of Lemma 4.2, this feature cannot be avoided since $\chi^{\prime}$ satisfies an ODE of type (28) with non-zero boundary conditions. As a consequence, we need to assume once again that the non-linear gradient term vanishes at the boundary of the domain in order to control the explosion of the weight function.

Since the weight function does not explode in time anymore, we can work with $\rho(t)=\beta e^{-\beta t}$ as in Section 3. Thus we have

[TABLE]

where ${\mathcal{W}}_{k}$ are constructed from the formulas of Lemmas 4.2–4.3. We now provide sufficient conditions to ensure that $(\psi^{x})_{x\in{\mathcal{O}}}$ and $(\psi^{x}b_{i}(x)\cdot{\mathcal{W}}(T_{\emptyset},x,W))_{x\in{\mathcal{O}}}$, $i=1,\ldots,m$, are uniformly bounded in $L^{q}$. Let us define

[TABLE]

Clearly, the constant $C_{1}$ is finite if the following condition is satisfied:

[TABLE]

Now we consider the family $(\bar{\psi}^{x})_{x\in{\mathcal{O}}}$ given by

[TABLE]

The tools developed in Section 3.3 allow to derive sufficient conditions to ensure that the family $(\bar{\psi}^{x})_{x\in{\mathcal{O}}}$ is bounded in $L^{q}$ . Since $|\psi^{x}|\leq|\bar{\psi}^{x}|$, the boundedness on $(\psi^{x})_{x\in{\mathcal{O}}}$ in $L^{q}$ follows immediately. The same holds for $(\psi^{x}b_{i}(x)\cdot{\mathcal{W}}(T_{\emptyset},x,W))_{x\in{\mathcal{O}}}$ for $i=1,\ldots,m$.

5   Numerical Examples with Related Technical Discussions

In this section we demonstrate the application of our theoretical results with the help of several examples. We restrict our attention to semi-linear elliptic PDEs over a rectangular domain where the underlying diffusion is a Brownian motion. This allows us to generate unbiased samples for exit time and position using the walk on squares scheme as implemented in the numerical library exitbm developed by Lejay [24]. We note that our theoretical results are applicable more widely but the numerical estimates for semi-linear PDEs driven by Brownian motion over rectangular domain do not involve any bias. We performed all the numerical computations on a machine with 2,5 GHz Intel Core i5 processor and 4GB RAM.

5.1   Example 1

We begin with the following semi-linear PDE in one-dimension (ODE)

[TABLE]

with explicit solution $u(x)=\frac{\sqrt{2}}{\cosh(x)}.$ To test the validity of our probabilistic representation in Theorem 3.1, we reformulate the above PDE, in the form (7), as following

[TABLE]

For domain ${\mathcal{O}}=(-r,r),$ we know from standard literature that the first eigenvalue of Laplacian is given as $\lambda_{1}=\tfrac{\pi^{2}}{4r^{2}}.$ Then, by choosing $\beta=1,$ $p_{1}=\tfrac{1}{2}$ and $p_{3}=\tfrac{1}{2},$ the criteria for almost sure extinction of the branching Brownian motion (6) holds for any $r\leq\tfrac{\pi}{\sqrt{8}}$. In addition, the integrability condition as presented in Proposition 3.4 is satisfied for $q=1$ by $v(x)=\frac{\sqrt{2}}{\cosh(x)}.$ Therefore, our probabilistic representation of the solution given as

[TABLE]

is valid for any $r\leq\tfrac{\pi}{\sqrt{8}}.$ Once more, we make use of Proposition 3.4 and see that for $q=2,$ it is satisfied by v(x)=\sqrt{1+\lambda}\cos\bigl{(}\sqrt{\lambda}x\bigr{)} with $\lambda=6$ and $r\leq 0.338.$ This implies that $(\psi^{x})_{x\in{\mathcal{O}}}$ is uniformly bounded in $L^{2}$ and we are certain that it has finite variance. In Table 5.1, for $r=0.3,$ we illustrate numerical results for different values of starting position $x$ with 1$\times 10^{6}$ Monte Carlo sample paths which exhibits the accuracy of our estimator.

We observe that even for larger values of $r,$ our probabilistic representation provides an accurate estimator but with a standard deviation estimate which converges slowly. In Table 5.2, we vary the size of the interval to test the accuracy and notice that as more branching particles are generated for larger intervals, the estimator has higher variance and run time.

Remark 5.1**.**

If $r\geq\operatorname{arcosh}(\sqrt{2})$, the solution provided by the probabilistic representation takes values in $[0,1]$ and thus it does not coincide with $u(x)=\tfrac{\sqrt{2}}{\cosh(x)}$. In particular, in view of Proposition 3.1, it turns out that the sequence $(\psi_{n}^{x})_{n\in\mathbb{N}}$ is not uniformly integrable for $\operatorname{arcosh}(\sqrt{2})\leq r\leq\tfrac{\pi}{\sqrt{8}}$. In Table 5.3, we present the results which support our observation.

5.2   Example 2

We consider another semi-linear PDE in one dimension. It is given as

[TABLE]

with solution $u(x)=1+2\tan^{2}(x).$ Then, for domain ${\mathcal{O}}=(-r,r)$, $\beta=1$ and probability mass function as $p_{0}=0.25$ and $p_{2}=0.75,$ the branching Brownian motion goes extinct almost surely for $r\leq\tfrac{\pi}{2}.$ Next, to verify that our probabilistic representation holds, we work with the conditions in Proposition 3.5. We first compute

[TABLE]

From standard calculations, we get

[TABLE]

This gives us $\tilde{p}_{0}=0.25+\tfrac{0.75}{\cosh(\sqrt{2}r)}$ and $\tilde{p}_{2}=0.75\left(1-\tfrac{1}{\cosh(\sqrt{2}r)}\right).$ In addition, the integrability constant $C_{0}=\max\{2,1+2\tan^{2}(r)\}.$ Next, in view of (18), the choice of threshold is given as $\gamma=\tfrac{1}{\sqrt{4\tilde{p}_{2}(1-\tilde{p}_{2})}}.$ Then, for our probabilistic representation of the solution to be valid, we need that $C_{0}<\gamma$, which gives us that $r\leq 0.31.$ As in the previous example, we would like to compute a confidence interval for the estimate of the solution. Once again, by checking the conditions in Proposition 3.5, we get that for our estimator $(\psi_{x})_{x\in{\mathcal{O}}}$ to be uniformly bounded in $L^{2},$ we must choose $r\leq 0.146.$ We set $r=0.14,$ and obtain the following results for $1\times 10^{6}$ Monte Carlo sample paths in Table 5.4 which exhibits the accuracy of our estimator.

In the course of our experiments, we chose different values of $r$ and observed that even for higher values, we are able to obtain accurate estimates of the solution. For $r=0.3,$ we present the results in Table 5.5.

5.3   Example 3

Next, we consider a multidimensional semi-linear PDE

[TABLE]

with an explicit solution $u(x)=\tan(\sum^{d}_{i=1}x_{i})$ where $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}.$ We reformulate the above PDE as

[TABLE]

with ${\mathcal{O}}=(-r,r)^{d}$ for $r>0$ and set $\beta=d$ and $p_{3}=1.$ With $\lambda_{1}=\tfrac{d\pi^{2}}{4r^{2}},$ the branching Brownian motion goes extinct for $r\leq\tfrac{\pi}{4}.$

We first consider $d=2.$ For $r<\tfrac{\pi}{8},$ the integrability constant $C_{0}<1$ and thus our estimator will have finite moments. In order to check if our probabilistic representation is valid for even a larger domain, we verify the condition in Proposition 3.5 for $q=1.$ We need to compute

[TABLE]

In the absence of an explicit formula for the Laplace transform of exit time of Brownian motion from a multidimensional rectangular domain, we estimate it by sampling the exit time using the numerical library developed by Lejay [24]. The integrability constant needs to be less than the threshold $\gamma$ in Proposition 3.5 for our probabilistic representation to be valid. In view of (18), the choice of threshold is given as $\gamma=\sqrt[3]{\tfrac{4}{27\tilde{p}_{0}^{2}(1-\tilde{p}_{0})}}.$ From our estimation procedures, we deduce that $C_{0}<\gamma$ for $r<0.484.$ Here, we exhibit our numerical results for $5\times 10^{5}$ Monte Carlo sample paths with $r=0.48$ and different values of $x$ in Table 5.6. We observe that our estimator remains accurate but with an increased variance which is an expected effect due to an increased variance in the estimation of exit time and position of multi-dimensional Brownian motion using the library exitbm.

Typically, finite difference methods are quite commonly used to numerically solve PDEs. However, such methods are practically implementable and provide stable estimates only when $d\leq 3.$ Thus, we study PDE (36) in the case of $d=4$ to illustrate the broader applicability of our method. We extend the modules in library exitbm for our purpose using the existing functions. For $r<\tfrac{\pi}{16},$ the integrability constant $C_{0}<1$ and thus the estimator will have finite moments. Furthermore, through numerical estimation procedures, we have that $C_{0}<\gamma$ for $r<0.242.$ In Table 5.7, we exhibit our results with $5\times 10^{5}$ Monte Carlo sample paths for $r=0.24$ which shows that our probabilistic representation provides an accurate estimator of the true solution.

6   Conclusions

In this work, we introduced a probabilistic representation for the solution of semi-linear elliptic PDEs with polynomial non-linearity by using the theory of branching diffusion processes. We performed a detailed analysis to derive explicit conditions under which our representations remain valid. In the linear gradient case, we essentially established that it holds provided that the domain is small enough. In the general case, we need to assume further that the nonlinear gradient term vanishes at the boundary of the domain to balance the explosion of the weight functions. As illustrated by our analysis of the one-dimensional case, this feature is inherent to any automatic differentiation formula. This restricts the choice of semi-linear elliptic PDEs for which we could use our probabilistic representation to obtain numerical estimates. For semi-linear elliptic PDEs of form (7), we illustrated the applicability of our theoretical results to obtain numerical estimates with the help of several examples including multi-dimensional. Finally, in the numerical implementation of semi-linear elliptic PDEs driven by general diffusion processes, the error analysis due to the discretization of the process has been left for future research.

Acknowledgement

The authors would like to thank Pierre Henry-Labordère, Zhenjie Ren, Xiaolu Tan and Nizar Touzi for their valuable comments and suggestions. The first author research was conducted while at CMAP, École Polytechnique and is part of the Chair Financial Risks of the Risk Foundation. The second author acknowledges the financial support of ERC 321111 Rofirm.

Appendix A Proof of Proposition 4.2

For the purpose of clarity, we split the proof in two lemmas.

Lemma A.1**.**

For any $g$ bounded measurable, the map $x\mapsto{\mathbb{E}}[\int_{0}^{{\eta^{x}}}{e^{-\beta s}g(X^{x}_{s})\,ds}]$ belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ and

[TABLE]

Proof.

The proof consists of four steps. The first three parts essentially follow by repeating the arguments in Gobet [16]. Denote

[TABLE]

First step. Let us collect first preliminary results on the regularity of the function $\chi$. Recall that, as established in the proof of Lemma 3.1, we have the following representation

[TABLE]

where $G$ is the Green function of PDE (3.2.1) (see, e.g., Ladyženskaja et al. [23, Theorem 4.16.2] or Friedman [14, Theorem 3.16]). Furthermore, in view of Theorem 4.16.3 in [23], for any $0<t<T$, we can differentiate under the integral sign and thus $\chi$ belongs to ${\mathcal{C}}^{1,2}_{b}([t,T]\times{\mathcal{O}})$ and satisfies

[TABLE]

In addition, under Assumption 4.2, $\chi$ belongs to $C^{1,3}([t,T]\times{\mathcal{O}})$ (see, e.g., Friedman [14, Theorem 3.10]).

Second step. Next we show that $N_{\theta}:=D\chi(s-\theta,X^{x}_{\theta})Y^{x}_{\theta}$ is a martingale up to time $\zeta_{s}$, i.e., $(N_{\zeta_{s}\wedge\theta})_{\theta\geq 0}$ is a martingale. Using Itô’s formula, we obtain

[TABLE]

In addition, due to Lemma 4.1, there exists $t<s$ such that $\zeta_{s}\leq t$. As the maps $D\chi$ and $D^{2}\chi$ are bounded in $[s-t,s]\times{\mathcal{O}}$, we conclude that $(N_{\zeta_{s}\wedge\theta})_{\theta\geq 0}$ is a martingale.

Third step. Then let us show that

[TABLE]

Indeed, Itô’s formula yields that

[TABLE]

Successively using Markov property and Itô’s isometry, we obtain

[TABLE]

Once more, successively using the results that $N_{\zeta_{s}\wedge r}={\mathbb{E}}[N_{\zeta_{s}}\,|\,{\mathcal{F}}_{r}]$ and $\int_{0}^{\zeta_{s}}{\theta_{s}(r,X_{r}^{x})\,dr}=1$, we deduce that

[TABLE]

Fourth step. We are now in a position to complete the proof. From Fubini’s theorem we obtain

[TABLE]

Finally, in view of the estimates in Proposition 4.2, we can apply differentiation under the integral sign in (38). This completes the proof. ∎

Lemma A.2**.**

The map $x\mapsto{\mathbb{E}}[e^{-\beta{\eta^{x}}}h(X^{x}_{{\eta^{x}}})]$ belongs to ${\mathcal{C}}^{1}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ and

[TABLE]

Proof.

This proof is similar to the proof of Lemma A.1 and we split it in three parts. Denote

[TABLE]

First step. Let us collect first preliminary results on the regularity of the function $\varphi$. Recall that $\varphi\in{\mathcal{C}}^{2}({\mathcal{O}})\cap{\mathcal{C}}(\bar{{\mathcal{O}}})$ satisfies the following PDE:

[TABLE]

In addition, since $h\in{\mathcal{C}}^{1,\alpha}(\bar{{\mathcal{O}}})$ by assumption, it is known that $\varphi\in{\mathcal{C}}^{1,\alpha}(\bar{{\mathcal{O}}})$ (see, e.g., Gilbarg and Trudinger [15, Theorem 8.34]). Furthermore, under Assumption 4.2, the function $\varphi$ belongs to $C^{3}({\mathcal{O}})$ (see, e.g., Gilbarg and Trudinger [15, Theorem 6.17]).

Second step. Next we show that $N_{s}:=e^{-\beta s}D\varphi(X^{x}_{s})Y^{x}_{s}$ is a martingale up to time $\eta^{x}$, i.e., $(N_{{\eta^{x}}\wedge s})_{s\geq 0}$ is a martingale. Using Itô’s formula, we obtain

[TABLE]

Thus $(N_{{\eta^{x}}\wedge s})_{s\geq 0}$ is a local martingale. To conclude, it remains to observe that for any $s\geq 0$,

[TABLE]

since $D\varphi$ is bounded and $\sup_{0\leq r\leq s}\{|Y^{x}_{r}|\}$ admits finite moments.

Third step. We are now in a position to complete the proof. Indeed, Itô’s formula yields for all $s\geq 0$,

[TABLE]

Using Itô’s isometry, we obtain

[TABLE]

Finally, successively using the results that $N_{\zeta_{T}\wedge r}={\mathbb{E}}[N_{\zeta_{T}}\,|\,{\mathcal{F}}_{r}]$ and $\int_{0}^{\zeta_{T}}{\theta_{T}(r,X_{r}^{x})\,dr}=1$, we deduce that

[TABLE]

∎