E Koorambas
email: elias.koor@gmail.com

The mass gap problem in the general gauge model framework

Keywords: number(s): 11.15.-q, 12.38.-t Gauge field theories, Quantum Chromodynamics, Yang-Mills problem

We demonstrate that Wu's version of the quantum chromodynamics (QCD) predicts mass gap Δ > 0 for the compact simple gauge group SU (3). This provides a solution to the second part of the Yang-Mills problem.

Introduction

The laws of quantum physics are for the world of elementary particles what Newton's laws of classical mechanics are for the macroscopic world. Almost half a century ago, Yang and Mills introduced a remarkable new framework to describe elementary particles by using geometrical structures. Since 1954, the Yang-Mills (YM) theory [START_REF] Young | Conservation of isotopic spin and isotopic gauge invariance[END_REF] has been the foundation of contemporary elementary particle theory. Although the predictions of the YM theory have been tested in many experiments [START_REF] Caso | [END_REF], its mathematical foundations remain unclear. The success of the YM theory in describing the strong interactions between elementary particles depends on a subtle quantum-mechanical property, termed the mass gap: quantum particles have positive masses even though classical waves travel at the speed of light. The mass gap has been discovered experimentaly and confirmed through computer simulations, but is yet to be understood theoretically [START_REF] Jaffe | Quantum Yang-Mills Theory[END_REF], [START_REF] Creutz | Monte Carlo Study of Quantized SU (2) Gauge Theory[END_REF], [START_REF]Quarks and strings on a lattice[END_REF]. The theoretical foundations of the YM theory and the mass gap may require the introduction of new fundamental ideas in physics and mathematics [START_REF] Young | Conservation of isotopic spin and isotopic gauge invariance[END_REF]. As explained by Arthur Jaffe and Edward Witten in the Caly Mathematical Institute (CMI) problem description [START_REF] Jaffe | Quantum Yang-Mills Theory[END_REF], the YM theory is a generalization of Maxwell's theory of electromagnetism, in which the basic dynamical variable is a connection on a G-bundle over four-dimensional space-time. The YM theory is the key ingredient in the Standard Model of elementary particles and their interactions. A solution to the YM problem, therefore, would both place the YM theory on a firm mathematical footing and demonstrate a key feature of the physics of strong interactions.

The foundation of the YM theory and the mass gap problem could be formulated as follows: Prove that, for any compact simple gauge group G, there exists a quantum YM theory of R 4 , and that this theory predicts a mass gap Δ > 0. Or, more explicitly, given a simple Lie group G (e.g. SU [START_REF] Jaffe | Quantum Yang-Mills Theory[END_REF] or SU (3)), show that A) there exists a full renormalized quantum version of YM theory on R 4 Assuming that the quantum chromodynamics (QCD) is the valid theory [1] , the first part of CMI problem description, the problem of the foundation of the YM theory, is a technical rather than a fundamental physical problem. The goal of this paper is to investigate the possibility of a solution to the second part of the official (CMI) problem description, namely the mass gap problem, by utilizing Wu's formulation of QCD ( [START_REF] Wu | Gauge Model with Massive Gravitons[END_REF], [START_REF] Wu | Quantum Gauge Theory of Gravity[END_REF], [START_REF] Wu | Gauge Theory of Gravity[END_REF], [START_REF] Wu | N Renormalizable Quantum Gauge Theory of Gravity[END_REF], [START_REF] Wu | [END_REF], [10], [11], [START_REF] Wu | Some discussions on strong interaction[END_REF]).

The physics of the mass gap problem

In the standard QCD model, processes that probe the short-distance structure of hadrons predict that quarks inside the hadrons interact weakly [START_REF] Caso | [END_REF], [START_REF] Dokshitzer | XXIXth International Conference on High Energy Physics, ICHEP98[END_REF], [START_REF] Ellis | [END_REF]. Since coupling, g, is small, the classical QCD analysis is a good first approximation for these processes of interaction [START_REF] Caso | [END_REF], [START_REF] Dokshitzer | XXIXth International Conference on High Energy Physics, ICHEP98[END_REF], [START_REF] Ellis | [END_REF]. However, for YM theories in general, it is a general requirement of the renormalization group equations of the quantum field that coupling, g s increases in reverse law to the hadrons' momentum transfer, until momentum transfer becomes equal to the vector boson. Since spontaneous symmetry-breaking via the Higgs mechanism to give the gluons mass is not present in QCD ( [16], [17], [18], [19], [20], [21], [22], [23]), QCD contains no mechanism to stop the increase of coupling, g s . In consequence, quantum effects become more and more dominant as distances increase. Analysis of the behavior of QCD at long distances, which includes deriving the hadrons spectrum, requires foundation of the full YM quantum theory. This analysis is proving to be very difficult ( [START_REF] Jaffe | Quantum Yang-Mills Theory[END_REF], [START_REF] Creutz | Monte Carlo Study of Quantized SU (2) Gauge Theory[END_REF], [START_REF]Quarks and strings on a lattice[END_REF]). [1] Note: Lattice QCD is a non-perturbative approach to solving the QCD theory [START_REF] Callaway | A 1982 Microcanonical Ensemble Formulation of Lattice Gauge Theory[END_REF], [START_REF] Callaway | Lattice gauge theory in the microcanonical ensemble[END_REF]. Analytic solutions in low-energy QCD are hard due to the nonlinear nature of the strong force. This formulation of QCD in discrete space-time introduces a momentum cut off at the order 1/s, where s is the lattice spacing [START_REF] Dürr | Ab Initio Determination of Light Hadron Masses[END_REF], [START_REF] Bazavov | Nonperturbative QCD simulations with 2+1 flavors of improved staggered quarks[END_REF].

The mass gap problem can be expressed as follows:

Why

The mass gap in the framework of Wu's QCD

As we have seen, the Higgs mechanism ( [16], [17], [18], [19], [20], [21], [22], [23]) is excluded in the QCD theory (gluons are massless because the QCD Lagrangian has no spinless fields and, therefore, no obvious possibility for spontaneous symmetry-breaking). It follows that the only possible mechanism that can introduce the mass term in the gluon field while remaining consistent with the occurrence of mass gap is the Wu mass generator mechanism ( [START_REF] Wu | Gauge Model with Massive Gravitons[END_REF][START_REF] Wu | Quantum Gauge Theory of Gravity[END_REF][START_REF] Wu | Gauge Theory of Gravity[END_REF][START_REF] Wu | N Renormalizable Quantum Gauge Theory of Gravity[END_REF][START_REF] Wu | [END_REF][10][11][START_REF] Wu | Some discussions on strong interaction[END_REF], [START_REF] Koorambas | [END_REF][42][43][START_REF] Koorambas | Recent Developments in Bosons Research[END_REF][START_REF] Koorambas | [END_REF][START_REF] Koorambas | Quantum[END_REF]). The Wu mass generator mechanism introduces mass terms without violating the SU (3) gauge symmetry.

The quark field is denoted as

( ), ( 1, 2, 3) a xa   (1)
where a is color index and the flavor index is omitted.

Ψ is defined as follows: A  is a purely quantum phenomenon: 1 a A  is generated through non-smoothness of the scalar phase of the fundamental spinor fields [42], [START_REF] Kuzmenko | Proceedings of Institute of Mathematics of NAS of Ukraine[END_REF].

From the viewpoint of the gauge field 1 a A  generation described here, the gauge principle is an "automatic" consequence of the non-smoothness of the field trajectory in the Feynman path integral [42], [START_REF] Kuzmenko | Proceedings of Institute of Mathematics of NAS of Ukraine[END_REF]. The Lagrangian of the model is:

5 11 1 2 2 1 2 2 1 2 1 [ ( ) ] ( ) 4 1 ( ) [(cos sin )(cos sin )] 42 igA m T r A A K m T r A A T r aA aA aA aA KK                        (3)
where

1 1 1 1 1 [ , ] A A A ig A A             , 2 2 2 2 2 [ , ] A A A igta A A             (4)
This Lagrangian can be proved to have strict SU (3) c gauge symmetry. Since A  and 2

A 
are not eigenvectors of mass matrix, we can apply the following transformations:

1 1 2 2 1 2 cos sin sin cos G aA aA G aA aA           (5)
After these, the Lagrangian given by equation (3) changes into

i i i i gI ig aG ig aG m G G m G G G G                         (6) Where 0 , ( 1, 2)
i i m m m G G G m           , /2 i m m i GG    (7)
The Lagrangian gI  only contains interaction terms of gauge fields. gI  includes a massive gluon field, 1 G , with mass m ,and a massless gluon field, 2 G . Since there exist two sets of gluons, there may exist three sets of glueballs in mass spectrum, (g 1 g 1 ), (g 1 g 2 ) and (g 2 g 2 ), with the same spin-parity but different masses (see [START_REF] Wu | Some discussions on strong interaction[END_REF]).

Transformation ( 5) is pure algebraic operations which do not affect the gauge symmetry of the Lagrangian [START_REF] Wu | Gauge Model with Massive Gravitons[END_REF][START_REF] Wu | Quantum Gauge Theory of Gravity[END_REF][START_REF] Wu | Gauge Theory of Gravity[END_REF][START_REF] Wu | N Renormalizable Quantum Gauge Theory of Gravity[END_REF][START_REF] Wu | [END_REF][10][11][START_REF] Wu | Some discussions on strong interaction[END_REF]. They can, therefore, be regarded as redefinitions of gauge fields. The local gauge symmetry of the Lagrangian is still strictly preserved after field transformations. In other words, the symmetry of the Lagrangian before transformations is absolutely the same with the symmetry of the Lagrangian after transformations. We do not introduce any kind of symmetry breaking in this paper. In Quantum Field Theory (QFT) the number of particles is not constant, due to the particles' creation and annihilation. Therefore, while we can define the mass-state of the system as the mass-state of all particles or equivalent, we cannot define the mass-state of a specific particle. In Wu's version of QCD, nevertheless, the mass-state of the system has non-zero mass. A possible interpretation of the existence of both massive and massless gluon states is that Wu's mass generator mechanism introduces the phenomenon of mass gap. This is achieved in the following steps: First, in the Wu QCD version the non-zero gluon mass corresponds to the massive vector propagator [24] and the zero gluon mass corresponds to the massless vector propagator of the short-distance QCD [START_REF] Caso | [END_REF], [START_REF] Dokshitzer | XXIXth International Conference on High Energy Physics, ICHEP98[END_REF], [START_REF] Ellis | [END_REF]. For the given wave vectors 1 k and 2 k the propagators of massless and massive gluons are given by:

2 2 2 2 2 () g ki ki                 (8) 11 1 1 1 2 2 2 1 () kk i k g k M i M                (9)
where [ 2 ] 1 2 ,1 2

     , ( 10a
) 1 1 1 ( , ) g e e     , 2 2 2 ( , ) g e e     , ( 10b
) 1 1 , e x      2 2 e x      (10c) 1 g  , 2
g  are metrics of the space-time regions and 2 , in the momentum space, in terms of its hypersurface of the present ( [25]) at the origin of the propagators paths as follows

e  2
12 1 2 1 2 1 1 2 2 1 2 ( , ) 0 ( ) ( ) lim inf { m in ( ( ), ( )} k k S i k k k k                        (11a)
where

1 1 1 1 1 2 2 1 1 ( ) ( ) , k k k m i              (11b) 2 2 2 2 2 2 2 1 ( ) ( ) k k ki              (11c)
are functions in the momentum space. The subset

1 \ S S r  
of the hypersurface of the present is disconnected, since it is the disjoint union of the two propagator's half-planes of the present that corresponds to the different 1 and 2  indices. [START_REF] Jaffe | Quantum Yang-Mills Theory[END_REF] Note: In Wu's version of QCD, the letters propagator in the momentum space

U S k k R ak bk c           (11f)
1 2 1 2 1 2 ( ) ( ) 0 k k                  (11g)
Equation (11g) shows that the meet of the propagators

1 1 () k    and 2 2 () k   
is zero, because the propagator's half-planes of the present are disjoint in the momentum space. After some calculations, the Fourier transformation to the position space of equation (11g) is given by ( 4)

1 ( 4 ) 2 2 2 1 1 1 1 2 2 1 2 ( ) ( ) ( ) ( ) 0 x y x y x y x y                         (12)
where

and

, ( 14
) 1 1 1 4 () 11 1 1 1 1 1 4 2 2 2 0 1 ( ) lim ( ) (2 ) ik x y kk dk e x y i g k M i M                             (15)
) ( 4 ) 1 4 4 ( 4 ) 2 4 4 2 2 1 1 2 1 1 1 2 2 1 2 1 2 ( ) ( ) ( ) ( ) 0 x y x y d x d x x y x y d x d x                           2 2 2 4 () 2 2 2 2 2 42 0 2 ( ) lim ( ) (2 ) ik x y dk e x y g i ki                     16
Here, we integrate equation ( 16) in the following steps. First, by fixed 12 , yy, the equation ( 16) becomes

1 4 2 4 1 1 1 2 2 2 1 2 ( ) ( ) 0 x b d x x b d x                      (17)
Next, by considering an insulating color source, in Feynman gauge, the propagators magnitudes are constant in the position space. Therefore the integral of the equation ( 12) over 4 1 dx,

4 2 dx, gives (1)
( 2 )

1 2 1 1 2 2 1 4 2 4 ( ) ( ) 0 x b V x b V                     (18)
where

(1) 3 4 1 1 V d x dt    , ( 2
) 3 4 2 1 V d x dt    (19)
(1) 4

V

, ( 2 ) 4 V are the space-time volumes.

Equation ( 18) is written in style of (11g) as follows

12 1 2 1 2 1 1 2 2 1 2 ( , ) ( ) ( ) lim inf {m in( ( ), ( )} 0 x x S x x x x                        (20a)
where

1 1 1 1 1 1 1 ( ) ( ) , mx e x x x             (20b) 2 2 2 2 2 2 1 ( ) ( ) x x x            (20c)
are functions in the position space. The subset 2 \ S S r    of the hypersurface of the present is disconnected, since it is the disjoint union of the two propagator's half-planes of the present that corresponds to the different 1 and 2  indices. lie outside each other's light cone [25]; these propagators, therefore, are not casually connected [START_REF] Weinberg | Gravitational and Cosmology[END_REF]. Propagators ). m is the mass of the particle that is constrained to the above space-like portion of space-time [START_REF] Bilaniuk | Meta' Relativity[END_REF], [START_REF] Feinberg | Possibility of Faster-Than-Light Particles[END_REF]. This particle of mass m either travels forward in time with imaginary mass or travels backward in time with real mass [START_REF] Bilaniuk | Meta' Relativity[END_REF], [START_REF] Feinberg | Possibility of Faster-Than-Light Particles[END_REF].

1 2 U U       and 2 1 2 \ ( , ) \ {0} S S r U U              (20d) 2 1 11 12 11 12 { {( , ) : 0} (0, ), U S x x R ax bx c             (20e) 2 1 21 22 21 22 { {( , ) : 0} ( , 0) U S x x R ax bx c             (
Next, a basic principle in nuclear physics which states that the combined operations of charge conjugation (C), time reversal (T), space inversion (P) in any order is an exact symmetry of the strong force [START_REF] Kostelecky | Data Tables for Lorentz and CPT Violation[END_REF], [START_REF] Sozzi | Discrete symmetries and CP violation[END_REF]. The CPT symmetry is conserved only if our theory respects the Lorentz invariant and the microcasuality principle [START_REF] Griffiths | Introduction to Elementary Particles[END_REF], [START_REF] Streater | PCT, spin and statistics, and all that[END_REF].

Here, the observations of particles that are constrained to the space-like portion of space-time, without violating the microcasuality principle are allowed by the following postulate:

The particle with space-like four-momentum that travels backward in time is equivalent to the mass gap (Δ) that is observed at the rest frame.

The gap in the particle's space-like four-momentum that travels backward in time is equivalent to the superluminal particle that travels forward in time and is observed at the laboratory frame.

The postulate predicts the observations of either superluminal particles that travel forward in time [START_REF] Giani | On velocities beyond the speed of light c[END_REF], [START_REF] Adam | Measurement of the neutrino velocity with the OPERA detector in the CNGS beam[END_REF] or a mass gap (Δ) at the rest frame [START_REF] Young | Conservation of isotopic spin and isotopic gauge invariance[END_REF], [START_REF] Jaffe | Quantum Yang-Mills Theory[END_REF].

From the above, equations (11g) and (20a) predict that the Wu Hamiltonian is not supported in an arbitrary space-like region 2 S  that corresponds to an arbitrary energy scale

( 2 4 2 2 / c S   
). Therefore, the Wu Hamiltonian in space-time is bounded by Δ.

Finally, in Wu's version of QCD, the vacuum state has zero energy (E=0), corresponding to G 2 (massless gluon). The excitation state above the vacuum energy has non-zero energy E m , corresponding to G 1 (massive gluon). The latter predicts that the Wu Hamiltonian is not supported in the region   0,  with 0  and is bounded below by  close to the gluon mass.

Renormalization of the Model

In order to quantize the suggested gauge field theory of nuclear forces in the path integral formulation, we first have to select gauge conditions [START_REF] Wu | [END_REF], [42]. To fix the degree of freedom of the gauge transformation, we must select two gauge conditions simultaneously: one for the massive gluons 1 G  and another for the massless gluons 2 G  [START_REF] Wu | [END_REF], [42] For instance, if we select temporal gauge condition for massless gluons 2 G  , 2 0, G  (21) there still exists a remainder gauge transformation degree of freedom, because the temporal gauge condition is unchanged under the following local gauge transformation:

  1 1 22 1 / sin , G U G U ig U U        (22)
where

0, t U  ( ). U U x   (23)
In order to make this remainder gauge transformation degree of freedom completely fixed, we have to select another gauge condition for gluons 1 G  . For instance, we can select the following gauge condition for gluons 1 G  , 1 0. G    (24) If we select two gauge conditions simultaneously, when we quantize the theory in path integral formulation there will be two gauge fixing terms in the effective Lagrangian. The effective Lagrangian can then be written as:

1 1 2 2 1 1 2 2 12 11 , 22 a a a a eff f f f f f f M M aa            (25)
where [START_REF] Callaway | A 1982 Microcanonical Ensemble Formulation of Lattice Gauge Theory[END_REF] then the propagator for massive gluons 1 G  is:

2 1 1 ( ), a a f f G   1 2 2 ( ). a a f f G   (26) If we select 21 , a fG   
2 2 2 2 12 ( ) / ( ) (1 ) / ( ) . g g F i k i k m i g a k k k a m                          (28)
If we let k approach infinity, then In this case, and according to the power-counting law, the gauge field theory of nuclear forces with massive gluons suggested in this paper is a kind of renormalizable theory [START_REF] Wu | [END_REF], [42]. A strict proof on the renormalizability of the gauge field theory can be found in [START_REF] Wu | [END_REF].

Discussion

Hadron colliders, such as the Tevatron, can place the strongest limits on the new colored particles (squarks and gluinos in supersymmetry or KK excitations of quarks and gluons in models with universal extra dimensions [35][36]. Typically, limits of∼ 200 GeV are obtained for new colored particles. Here, the predicted gluon mass is close to the mass gap state Δ, thus gluon mass is not subject to the Hadron collider's limits [35][36].

In the Lorentz gauge, the free Lagrangian of Wu's QCD field theory (equation ( 6), with 0 gI ) yields to the following field equations:

2 1 2 1 1 2 2 2 2 2 2 2 22 1 1 ( ) ( ), 1 ()
G G r G r c t L G Gr ct             (30)
L the length of the system given by

L mc   (31)
Where  the Planck constant, c is the speed of light, m the mass of the system. If the source is a point color at the origin, only the time-components 02 2 C  , 01 1

C  are nonvanishing. The short-distance solution of the eq (30) that corresponds to the massless gluons is given as follows:

20 () 4 | | s a Gr r   (32)
where s a is the quark-gluon coupling.

Gluons that respect the solution (32) are mediators of the long-range nuclear force at high energy scale. The long-distance solution of the eq (30) corresponds to the massive gluons and is given as follows:

| |/ 10 () 4 | | rL s a G r e r    (33)
Gluons that respect the solution [START_REF] Griffiths | Introduction to Elementary Particles[END_REF] are mediators of the short-range nuclear force at low energy scale. Here, the massive and massless gluons lie outside each other's light cone, because of eq (20a). The term [START_REF] Kostelecky | Data Tables for Lorentz and CPT Violation[END_REF] in the eq [START_REF] Bazavov | Nonperturbative QCD simulations with 2+1 flavors of improved staggered quarks[END_REF] and ( 33) reveals as a mass gap. Supposing that the range of the strong force is about the radius of the proton, equation [START_REF] Kostelecky | Data Tables for Lorentz and CPT Violation[END_REF] predicts mass gap of about 1GeV.

We find that, in Wu's version of QCD, free system's energies are indeed bounded in spacetime. Consequently, in space-time, the Wu has its minimum at 0 E  . Τhe first state m  is the vacuum vector. We now notice that, in space-time, the spectrum of the Wu-YM Hamiltonian is not supported in region   0,  , with 0  .

According to the Wu-YM theory, therefore, over compact gauge groups SU (3) there is always a vacuum vector m  and a mass gap Δ. Wu's mass generator mechanism can always be applied to any YM type of action. Thus, any YM theory accounts for a mass gap.

The proposed model is a massive theory without a symmetry breaking mechanism and is therefore a model with a quantum mass gap. Thus, the massless YM theory is altered with a quantum massive model and the original problem could be solved.

Conclusion

We have shown that the Wu gauge model can introduce mass terms while being consistent with the existing mass gap Δ. In this model the mass gap is introduced by Wu's mass generator mechanism.

Since the gluon field 1

G is massive, whereas the gluon field 2 G is massless, we conclude that system's free energies in Wu's version of QCD are bounded in space-time.

We also notice that, in space-time, the spectrum of the Wu-YM Hamiltonian is not supported in region   0,  with 0  . Wu's mass generator mechanism can always be applied to any YM type of action. Therefore any YM theory accounts for a mass gap.

The proposed model is a massive theory without a symmetry breaking mechanism and is therefore a model with a quantum mass gap. Thus, the massless YM theory is altered with a quantum massive model and the original problem could be solved.

1 aA  and 2 aA 2 aA

 122 ψ α form the fundamental representative space of SU (3) c. In this version of the Wu gauge model, the gauge fields  are introduced .Gauge field  is introduced to ensure the local gauge invariance of the theory. The generation of gauge field 1 a

ex  , 2 x

 2  are vectors of the basis induced by the selection of local coordinates 1 on the neighborhood of the given points. Furthermore, due to equation (10a), the wave vectors 1 k  and 2 k  lie in different regions in the momentum space in their own right. Here, we express the meet of the propagators



  are the propagator's half-planes of the present.Substituting equations (11d), (11c), (11d) to equation (11a), we calculate the meet of the

  20f) the 2 :0 r ax by c    is line such that the points 11 12 of S  lie on different sides of 2 r and 1 2 , UU   are the propagator's half-planes of the present. Equation (20a) shows that the propagator's half-planes of the present are disjoint in the position space, because the meet of

  of the disjoint space-time regions. The arbitrary spacelike region2 S corresponds to an arbitrary energy scale (

  -YM Hamiltonian wu H has spectrum bounded below. Let us denote this energy minimum of the Wu-YM Hamiltonian wu H by m E , and let us denote the lowest energy state by m  . By shifting the original Hamiltonian H by m

				E	, the
	new Hamiltonian w u H	H E 	m


Jean-Paul Morillon
email: paul.morillon@univ-reunion.fr

J.-P Morillon

Solving Wentzell-Dirichlet Boundary Value Problem with superabundant data using reflecting random walk simulation

Keywords: Monte Carlo method for linear BVP, Wentzell boundary condition, Reflected diffusion, Probabilistic representation, Stochastic numerical method PACS 90-08, 90C-15

In this paper, we are interested in numerical solution of some linear boundary value problems with Wentzell's boundary part and superabundant data on this part, by the means of simulation of reflected random walks. We use a probabilistic interpretation of solution, assuming that the diffusion coefficient and the boundary data are sufficiently smooth, and applying Itô's formula. From this stochastic representation of solution, we extend the algorithm obtained for mixed standard boundary conditions to the case of diffusion-reflection on the boundary, so called Wentzell's boundary condition. We then obtain numerical results by applying the stochastic method based upon this generalized algorithm.

1 Introduction A large number of numerical studies have been devoted to mixed boundary value problems (BVPs) using some deterministic methods: numerical solution can be performed by the means of usual or mixed finite element methods after a variational formulation of a class of standard BVPs.

Probabilistic methods can also be relevant: see, for instance, Monte Carlo methods for Neumann and Dirichlet BVPs in [START_REF] Kushner | Probabilistic methods for finite difference approximations to degenerate elliptic and parabolic equations with Neumann and Dirichlet boundaryconditions[END_REF] where the equations are discretized and specific treatments in the neighborhood of the boundary are applied.

The aim of this paper is to investigate probabilistic numerical approach to the mixed BVP with generalized Wentzell condition on a boundary part, by the means of simulation of reflected random walks. In particular, the BVPs with Dirichlet, Neumann or Robin reflection conditions are known to have integral representations [START_REF] Freidlin | Functional Integration and Partial Differential Equations[END_REF] which are used to obtain algorithms which compute the solution by implementing a reflected and/or absorbed random walk simulation [START_REF] Souza De Cursi | Numerical methods for linear boundary value problems based on Feyman-Kac representations[END_REF][START_REF] Morillon | Numerical solutions of linear mixed boundary value problems using stochastic representations[END_REF]. In this paper,we use a probabilistic interpretation of solution, assuming that the diffusion coefficient and the boundary data are sufficiently smooth, and applying Itô's formula. From this stochastic representation of solution, we extend the algorithm obtained for mixed standard boundary conditions to the case of diffusion-reflection on the boundary, so called Wentzell's boundary condition. We then obtain numerical results by applying the stochastic method based upon this generalized algorithm.

The limitations and the main advantages of these stochastic methods are already presented, for instance, in [START_REF] Souza De Cursi | Numerical methods for linear boundary value problems based on Feyman-Kac representations[END_REF][START_REF] Morillon | Numerical solutions of linear mixed boundary value problems using stochastic representations[END_REF]. In particular, the ensuing numerical methods do not require the storage of a grid discretization in computer memory, nor specific treatments related to domain local geometry. Programming is short, easy to check step by step. A contrario, the usual finite element method (FEM) leads to approximate Wentzell-Dirichlet BVP by means of a refinement of FEM on a boundary layer of thickness ε [START_REF] Raynaud | Approximation par éléments finis de problèmes de transmission raide[END_REF].

In this paper, we first introduce a Wentzell-Dirichlet BVP. Let G be a bounded open set of R 2 with a regular boundary Γ = ∂ G partitioned into two boundary parts Γ D and Γ V : Γ = Γ D ∪ Γ V and Γ D ∩ Γ V = / 0. A mixed Wentzell-Dirichlet BVP is the following one:

             - 1 2 a G ∆ u = f in G γ 0 u = g D on Γ D ∂ u ∂ ν + 1 2 a V ∂ 2 u ∂ τ 2 = g V on Γ V (1)
where the function u, defined from G ⊂ R 2 into R, is to be determined, and the data are the sufficiently smooth functions:

f : G -→ R a G : G -→ R + g D : Γ D -→ R a V : Γ V -→ R + g V : Γ V -→ R
with a G > a 0 > 0, a V > a 0 > 0 (a 0 is a strictly positive constant), and where ν is the inward unit normal vector and τ is the unit tangential vector defined on Γ V . If σ 1 is the Lebesgue measure on R, then we here assume: σ 1 (Γ D ) = 0 and σ 1 (Γ V ) = 0. When Γ is a polygonal boundary, we will impose the Wentzell condition at least on one straight line segment, and the Dirichlet condition on the other segments.

The boundary condition on Γ V = Γ \ Γ D so called Wentzell's boundary condition [START_REF] Ventcel | On lateral conditions for multidimensional diffusion processes[END_REF] corresponds to the heat exchange between the body G and a thin layer of thickness ε with a thermal conductivity a V /ε.

The body is either partially or completely covered with a highly heat-conductive thin shell, and the quantity u is then the temperature. The "heat particle" which reached the boundary at x ∈ Γ V , diffuses on Γ V just as on a boundary layer, and moves back into G along the line normal to Γ V just as the rebound effect for Neumann boundary condition; then the diffusion on G takes over again [START_REF] Goldstein | Derivation and physical interpretation of general boundary conditions[END_REF].

This type of non-standard BVP stems from various physical phenomenae. Let us note that the same type of boundary conditions appears in absorbing boundary conditions of wave-scattering problems and, more generally, in multilayer structure problems.

Our numerical approach of Wentzell BVP is based on the limit of the following problem:

                                       - 1 2 a G ∆ u = f in G \ Γ h V - 1 2 1 h a V ∆ u = g h V in Γ h V γ 0 u = g D on Γ D ∂ u ∂ ν = g ν V on Γ V + h/2 u + = u - on G ∩ Γ h V ∂ u ∂ ν + = 1 h ∂ u ∂ ν - on G ∩ Γ h V (2)
when h tends to 0, so that the domain G ∪ Γ h V of the BVP (2) converges to the initial domain G and the BVP (2) approximates to the Wentzell-Dirichlet BVP [START_REF] Bally | The Euler scheme for stochastic differential equations: error analysis with Malliavin calculus[END_REF], and where

Γ h V = x ∈ R 2 | d(x, Γ V ) < h/2 is a boundary layer so that Γ h V is a neighbor- hood of Γ V , and Γ V + h/2 = x ∈ R 2 \ G ∪ Γ h V | d(x, Γ V ) = h/2 is a boundary part of G ∪ Γ h V .
The two last equalities are the transmission conditions between domain G and boundary layer Γ h V . The extended domain G∪Γ h V is the union of G and boundary layer Γ h V , the boundary of G ∪ Γ h V is the union of Γ D and Γ V + h/2. Therefore, for the functional approach of the BVP (1), we distinguish the diffusion on a boundary layer Γ V and the Neumann condition on this boundary, so we consider a mixed Wentzell-Dirichlet BVP with superabundant data:

                     - 1 2 a G ∆ u = f in G γ 0 u = g D on Γ D ∂ u ∂ ν + 1 2 a V ∂ 2 u ∂ τ 2 = g V on Γ V ∂ u ∂ ν = g ν V on Γ V (3)
where the Neumann condition on Γ V is additional.

When Γ = N j=1
Γ j is a polygonal boundary where Γ j is an open straight line seg- ment, we consider the Wentzell condition at least on one straight line segment Γ j , and we obtain the stochastic representations even if the boundary ∂ G of domain G is not necessary of class C 3 , nor sufficiently smooth : ∂ G can show corners or edges (see Section 2.). Many situations exist in which solutions to BVP are represented as function space integrals. Such stochastic representations can be used to evaluate them numerically using Monte Carlo methods [START_REF] Graham | Stochastic Simulation and Monte Carlo Methods -Mathematical Foundations of Stochastic Simulation[END_REF]. For instance, a second-order elliptic operator subject to normal reflection or to oblique derivative boundary conditions generates a diffusion process. For general definition of diffusions with reflection in a domain, see [START_REF] Stroock | Diffusion processes with boundary conditions[END_REF], and for efficient approximation, see for instance [START_REF] Bossy | A Symmetrized Euler Scheme for an Efficient Approximation of Reflected Diffusions[END_REF] and references therein. We consider here Wentzell's boundary condition corresponding to such diffusion phenomenon that a Markovian particle moves both by jumps and continuously in the state space until it dies at the time when it reaches the set where the particle is definitely absorbed.

The schedule of this paper is the following one : in Section 2, we introduce the probabilistic formulation of BVPs of Wentzell's type by means of Itô's formula associated with stochastic differential equations (SDEs). For these SDEs, we need to consider the reflecting Brownian motion and a local time on a boundary. For theoretical results on reflecting Brownian motion and on local time, see [START_REF] Hsu | Reflecting Brownian motion, boundary local time and the Neumann problem[END_REF][START_REF] Lions | Stochastic differential equations with reflecting boundary conditions[END_REF].

Section 3 is devoted to the resolution algorithm. More precisely, the stochastic methods is introduced, and the process discretization is described. Section 4 is devoted to some numerical experiments and applications of algorithm obtained in the previous section. Finally, we present some comments in the last section.

Representation of solutions of Wentzell boundary value problem for the Laplacian

Wentzell's BVP is a model for the heat transfer between a solid G and its environment when the boundary ∂ G is covered with a thin layer of a material with higher conductibility.

When G is a bounded convex domain in R d , there exists a unique solution in H 2 (G) when the data are sufficiently smooth [START_REF] Lemrabet | Problème aux limites de Ventcel dans un domaine non régulier[END_REF]. From a probabilistic point of view, when G is an open bounded set of R d with C 3 boundary, there exists a unique solution in C 2 b G when the data and their derivatives are bounded and sufficiently smooth [START_REF] Bony | Semi-groupe de Feller sur une variété à bord compacte[END_REF].

This section is devoted to the stochastic representation of solution of Wentzell's BVP by means of stochastic integrals. Let us recall that the representations according to the Robin, Neumann and mixed boundary conditions are already presented in [START_REF] Morillon | Numerical solutions of linear mixed boundary value problems using stochastic representations[END_REF][START_REF] Morillon | Résolution de systèmes de réaction-diffusion; Méthodes déterministes, stochastiques et fonctionnelles à ε -près[END_REF], and that the theoretical background can be view, for instance, in [START_REF] Freidlin | Functional Integration and Partial Differential Equations[END_REF][START_REF] Cessenat | Méthodes probabilistes pour les équations de la physique[END_REF].

Let us consider an open bounded region G in R d (the space dimension d being fixed: d = 2 in this paper). Let us denote x ∈ R d the space variable and ν the unit inward normal defined on the boundary ∂ G.

We will show how, from stochastic differential equations, we can represent solution of Wentzell's BVP. We will use the expectation of integrals, i.e. functionals of trajectories and random variable which are solutions of SDE. The application of Itô's formula associated with SDE leads to a representation of solutions of Wentzell's BVP.

From a functional point of view, the stochastic representations of solutions of differential problems with general boundary conditions have only been obtained according to regularity hypotheses on the geometry of domain G, and on the boundary data (see [START_REF] Freidlin | Functional Integration and Partial Differential Equations[END_REF] for the probabilistic representation of solutions of first and second BVP). In particular, the solutions are assumed of class C 2 b G , i.e. with bounded first and second order continuous derivatives.

Nevertheless, we herein are using the stochastic representations obtained in applying Itô's formula even when the data are not necessary smooth. This procedure previously proved to be effective for the discretized equations [START_REF] Haji-Sheikh | Handbook of Numerical Heat Transfer[END_REF][START_REF] Cessenat | Méthodes probabilistes pour les équations de la physique[END_REF][START_REF] Kushner | Numerical Methods for Stochastic Control Problems in Continuous time[END_REF], for the Dirichlet problem [START_REF] Souza De Cursi | Numerical methods for linear boundary value problems based on Feyman-Kac representations[END_REF], and for the Robin, Neumann or mixed BVPs [START_REF] Morillon | Numerical solutions of linear mixed boundary value problems using stochastic representations[END_REF][START_REF] Morillon | Résolution de systèmes de réaction-diffusion; Méthodes déterministes, stochastiques et fonctionnelles à ε -près[END_REF] in bi-or tridimensional domains.

More precisely, Itô's formula is used even if solutions of BVP belong to Sobolev spaces of order 1 and are not necessary of class C 2 b G . In the sequel, we omit the arguments (X x t (ω) ) in the functions a G , a V , a j for brevity.

Wentzell-Dirichlet problem

Let us consider the following bidimensional BVP:

                     - 1 2 a G ∆ u = f in G γ 0 u = g D on Γ D ∂ u ∂ ν + 1 2 a V ∂ 2 u ∂ τ 2 = g ν V + g τ V on Γ V ∂ u ∂ ν = g ν V on Γ V (4)
where the superabundant data g N of problem (3) corresponds to g ν V and verifies g V = g ν V + g τ V , and where g ν V and g τ V : Γ V -→ R are also sufficiently smooth functions. Let us introduce the Markov process defined by the two-dimensional Laplacian operator within the domain G, by the one-dimensional Laplacian operator on the boundary part Γ V , and by reflection on the boundary in the direction of inward normal ν.

For sake of simplicity, the initial point x ∈ G being fixed, we denote: X t ≡ X x t , ξ t ≡ ξ x t , and θ ≡ θ x , where X t is the corresponding Markov process, ξ t the reflection process, and θ the hitting time of

R d \ (G ∪ Γ V ).
For the bidimensional case, we have: ∀t ≥ 0,

X t = x + √ a G W t + t 0 1 Γ V (X s ) √ a V τ(X s ) dβ ξ s + t 0 1 Γ V (X s )ν(X s ) dξ s , X 0 = x, ξ 0 = 0,
where: -W t is the standard Wiener process with values in R 2 , -β t is the standard Wiener process with values in R, independent of W t , -ξ t is a non decreasing process which increases only when t hits the set Λ defined by: Λ = s > 0 X x s ∈ Γ V , so that: ξ 0 = 0 and

t 0 1 G (X s ) dξ s = 0, ξ t is the so-called local time on Γ V , -1 A is the indicator of set A,
τ is the unit tangential vector, and ν is the unit inward normal vector. We apply Itô formula to the function u (X t ):

u (X t ) -u (X 0 ) = t 0 1 G (X s ) ( ∇u(X s ), √ a G dW s ) + t 0 1 2 a G ∆ u(X s )1 G (X s ) ds + t 0 1 Γ V (X s ) ∇u(X s ), √ a V τ(X s )dβ ξ τ s + t 0 1 2 a V ∂ 2 u ∂ τ 2 (X s )1 Γ V (X s ) dξ τ s + t 0 1 Γ V (X s ) (∇u(X s ), ν(X s ) dξ ν s )
where ξ τ s is the sojourn time when the process X s diffuses on Γ V and ξ ν s is the local time after diffusion when the process X s is reflected inward of G, then we use the data of problem (4):

u (X t ) -u (x) = t 0 1 G (X s ) (∇u(X s ), √ a G dW s ) - t 0 f (X s )1 G (X s ) ds + t 0 1 Γ V (X s ) ∇u(X s ), √ a V τ(X s )dβ ξ τ s + t 0 g τ V (X s )1 Γ V (X s )dξ τ s + t 0 g ν V (X s )1 Γ V (X s ) dξ ν s
On taking expectation of both sides of this equality, we get:

E [ u(X t ) ] -u(x) = -E t 0 f (X s )1 G (X s ) ds + E t 0 g τ V (X s )1 Γ V (X s ) dξ τ s + E t 0 g ν V (X s )1 Γ V (X s ) dξ ν s The hitting time θ of R d \ G ∪ Γ V is defined by: θ ≡ θ x = inf {t > 0 | X x t / ∈ G ∪ Γ V } .
and is denoted by θ to avoid confusion with the tangential vector τ. We obviously assume: Prob {θ x < ∞} = 1 for every x ∈ G since a G is strictly elliptic (see a similar result in [START_REF] Friedman | Stochastic differential equations and applications[END_REF], page 144).

The application of Itô's formula to u ( X t ) leads to the following representation of the solution of (4) in the form:

u(x) = E [R ] , x ∈ G, (5)
with the random variable:

R = g D (X θ ) + θ 0 f (X t )1 G (X t ) dt - θ 0 g τ V (X t )1 Γ V (X t ) dξ τ t - θ 0 g ν V (X t )1 Γ V (X t ) dξ ν t (6)
Let us mention that some functional results regarding SDE with reflection on Γ V , and associated BVP with general boundary conditions are given in [START_REF] Karoui | Processus de réflexion[END_REF][START_REF] Karoui | Un problème de réflexion et ses applications au temps local et aux équations différentielles stochastiques sur R, cas continu[END_REF].

If we adopt the modified sojourn time on Γ V : σ τ t = inf{ s | ξ τ s > t }, then we can write:

E θ 0 g τ V (X t )1 Γ V (X t ) dξ τ t = E ξ τ θ 0 g τ V (X σ τ t )1 Γ V (X σ τ t ) dt where X σ τ t is a Markov process with values in Γ V , t ∈ [0, θ ], ξ τ t ∈ [ξ τ 0 , ξ τ θ ] = [0, ξ τ θ ], and σ τ t ∈ σ τ ξ τ 0 , σ τ ξ τ θ = [0, θ ].
The numerical method given below in Section 3 will be based on the above ex-

pectation E ξ τ θ 0 g τ V (X σ τ t )1 Γ V (X σ τ t )
dt where the Markov process X σ τ t diffuses on Γ V , then purely and simply reflected by Γ V , so the Wentzell diffusion boundary condition is considered as a degenerated elliptic equation in a boundary layer.

Wentzell-Dirichlet BVP in a polygonal domain

Let us denote G a bounded open set of R 2 with a polygonal boundary

Γ = N k=1 Γ k , where Γ k is an open straight line segment ] S k-1 , S k [ or ] S N , S 1 [.
Let us consider the Wentzell-Dirichlet problem:

                                     - 1 2 a G ∆ u = f in G γ 0 u = g D on N 1 i=1 Γ k i ∂ u ∂ ν j + 1 2 a j ∂ 2 u ∂ τ 2 j = g ν j + g τ j on N j=N 1 +1 Γ k j ∂ u ∂ ν = g ν j on N j=N 1 +1 Γ k j a j ∂ u ∂ τ j -a j+1 ∂ u ∂ τ j+1 = 0 on S j (7)
with the superabundant data g N of problem (3) corresponds to g ν j , j = N 1 + 1, . . ., N and verifies g j = g ν j + g τ j , and where the unknown u is defined in G ∪ Γ with real values, N 1 verifies N 1 < N, ν j is the unit inward normal vector, τ j the unit tangential vector on Γ k j , j = N 1 + 1, . . ., N, the vertices S j exist and are defined by S j = Γ k j ∩ Γ k j+1 if j and j + 1 ∈ {N 1 + 1, . . ., N}, and the data a G , a j , j = N 1 + 1, . . . , N, are strictly positive functions so that: a G ≥ a 0 > 0 and a j ≥ a 0 > 0, j = N 1 + 1, . . ., N where a 0 is positive constant.

In order to avoid solutions defined up to an additive constant when the data are compatible, we assume that the interior of Γ D is not empty.

The elliptic part of the boundary condition on Γ k j is the tangential derivative of order 2. The Wentzell boundary condition is applied on the open set

Γ V ≡ N j=N 1 +1 Γ k j .
The point boundary conditions are the transmission conditions of the boundary operator.

Let us mention [START_REF] Lemrabet | Problème aux limites de Ventcel dans un domaine non régulier[END_REF] where functional results to Wentzell BVP in a polygonal domain G. In particular, the solution u belongs to a Sobolev space of order 1.

Let us introduce the Markov process X t with values in G ⊂ R 2 using analogous notation: ∀t ≥ 0,

X t = x + √ a G W t + N ∑ j=N 1 +1 t 0 1 Γ k j (X s ) √ a j τ j (X s ) dβ ξ s + N ∑ j=N 1 +1 t 0 1 Γ k j (X s )ν j (X s ) dξ s + N ∑ j=1 χ P(Γ V ) Γ k j ∪ Γ k j+1 t 0 1 S j (X s )a ρ j (ξ s ) (S j )τ ρ j (ξ s ) (S j ) dξ s , X 0 = x, ξ 0 = 0,
where χ P(Γ V ) (B ) equals to 1 if B is a part of Γ V , 0 otherwise, and ρ k is the uniform distribution on {k, k + 1} so that:

Prob {ρ k (ξ s ) = k} = Prob {ρ k (ξ s ) = k + 1} = 1 2 , k = 1, . . . , N -1, and Prob {ρ N (ξ s ) = N} = Prob {ρ N (ξ s ) = 1} = 1 2 .
As above, this process X t is absorbed by Γ D at the hitting time θ defined by:

θ = inf{ t > 0 | X t ∈ Γ D } .
The application of Itô's formula to u ( X t ) leads to the following representation of the solution of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] in the form:

u(x) = E [R ] , x ∈ G, (8)
with the random variable:

R = g D (X θ ) + θ 0 f (X t )1 G (X t ) dt - N ∑ j=N 1 +1 θ 0 g τ j (X t )1 Γ k j (X t ) dξ τ t - N ∑ j=N 1 +1 θ 0 g ν j (X t )1 Γ k j (X t ) dξ ν t ( 9
)
The problem ( 7) can be written in the following condensed form, using appropriate values of coefficients α, β , and γ:

       - 1 2 a G ∆ u = f in G β ∂ u ∂ ν + 1 2 α ∂ 2 u ∂ τ 2 -γu = g on Γ = ∂ G (10)
Let us note the analogy between problems ( 7) and [START_REF] Freidlin | Functional Integration and Partial Differential Equations[END_REF]. Nevertheless, from the functional point of view, some difficulties occur since the coefficients α, β and γ are discontinuous. Some partial results on this subject can be found in [START_REF] Bony | Semi-groupe de Feller sur une variété à bord compacte[END_REF][START_REF] Qingji | Solving boundary value problems by means of stochastic calculus[END_REF].

In the next section of this paper, we will focus on one stochastic method based upon the probabilistic representation of solution of BVP [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF].

Resolution algorithm for Wentzell-Dirichlet mixed boundary value problem

The main purpose of this section is to give an approach of the representations of Section 2 by means of realization of random processes, in order to compute the solutions of BVP with a Wentzell boundary part.

Let us recall that some efficient algorithms according to the Dirichlet, Neumann, Robin and mixed boundary conditions are already presented in [START_REF] Souza De Cursi | Numerical methods for linear boundary value problems based on Feyman-Kac representations[END_REF][START_REF] Morillon | Numerical solutions of linear mixed boundary value problems using stochastic representations[END_REF] for the Laplacian and in [START_REF] Morillon | Résolution de systèmes de réaction-diffusion; Méthodes déterministes, stochastiques et fonctionnelles à ε -près[END_REF] for the parabolic operator. Our present aim is then to extend the above techniques to our problem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF].

In order to estimate an expectation E[R], we simulate a large number (say NT ) of independent random variables with the same law as R. The mean of the NT realizations R i converges almost surely to E[R] (see [START_REF] Graham | Stochastic Simulation and Monte Carlo Methods -Mathematical Foundations of Stochastic Simulation[END_REF] where the almost sure convergence of the Monte Carlo methods is proved) :

1 NT NT ∑ i=1 R i a.s. -→ E[R] as NT → +∞.
Thanks to CLT, the rate of convergence of this approximation is known, it is of order 1/ √ NT .

In order to simulate the diffusion process X, we approach X by means of an Euler's scheme as in [START_REF] Bally | The Euler scheme for stochastic differential equations: error analysis with Malliavin calculus[END_REF][START_REF] Bally | The law of the Euler scheme for stochastic differential equations. I. Convergence rate of the distribution function[END_REF][START_REF] Morillon | Résolution de systèmes de réaction-diffusion; Méthodes déterministes, stochastiques et fonctionnelles à ε -près[END_REF] or the Box Muller algorihm. We need also the simulation of random walks by using some predefined generators of normally distributed pseudorandom numbers.

Let us consider the Wentzell-Dirichlet BVP [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] in a polygonal domain.

The solution u(x) to BVP ( 7) is represented as the mean of a random variable R, involving function of a non instantaneously reflected random walk, then absorbed random walk.

The representation (8)-( 9) is used and u(x) is computed as the mean of NT approximated values of R, i.e. R 1 , . . . , R NT , in order to obtain an estimation of u(x):

u(x) 1 NT NT ∑ i=1 R i .
Each value R i results from a simulation of a random walk on G, first reflected by the Wentzell boundary part Γ V , then absorbed by the Dirichlet boundary part Γ D . Each random walk is then finite. From now on, we shall assume that the hitting time θ is not too large and, in particular, the Dirichlet boundary part Γ D versus the Wentzell boundary part Γ V is significant.

The standard Wiener processes are simulated as follows. The mean of space step h > 0 being fixed, the process W t can be simulated by the sequence of random vectors W 0 , W 1 , . . . , W j defined by means of the recurrent formulae:

W 0 = x ∈ G, and W j+1 = W j + hD j , j ∈ N,
where D j is a random vector so that, if (e i , i = 1, 2 ) is the canonical basis of R 2 , D j verifies:

D j = r 1 e 1 + r 2 e 2 in G
where r 1 and r 2 are independent normally distributed pseudorandom numbers, and:

D j = r 3 τ + |r 4 | ν on Γ V
where r 3 and r 4 are independent normally distributed pseudorandom numbers. The absolute value |r 4 | leads to a reflection inward G.

Since the simulated process happens to leave the domain G, its trajectory resumes to the interior after a diffusion on Γ V or stops when the trajectory hits the Dirichlet boundary part Γ D .

The process is not instantaneously reflected by the boundary part Γ V when diffusion occurs.The elementary boundary reflections without diffusion are detailed in [START_REF] Morillon | Résolution de systèmes de réaction-diffusion; Méthodes déterministes, stochastiques et fonctionnelles à ε -près[END_REF], where the basic reflection algorithm is given, and general obliquely reflected diffusions (without diffusion on the boundary) are studied in [START_REF] Bossy | A Symmetrized Euler Scheme for an Efficient Approximation of Reflected Diffusions[END_REF]. The diffusion on the Wentzell boundary part is governed by the Wiener process β t which is simulated as the process W t but on the boundary layer:

Γ h V = x ∈ R 2 | d(x, Γ V ) < h/2 .
In the bidimensional case, the process diffuses on Γ V , along the tangential line given by τ, and is reflected inward the interior G along the normal direction given by ν.

Each simulation gives a discrete realization ( X j , 0 ≤ j ≤ NP ) of a random walk in G with NP steps, so that:

X 0 ∈ G \ Γ D , . . ., X NP-1 ∈ G \ Γ D , and X NP ∈ Γ h D = x ∈ R 2 | d(x, Γ D ) < h/2
, which trajectory is the polygonal line connecting successively the points X j and X j+1 , j = 0, . . . , NP -1, in order to form a random path.

The random processes are discretized in time according to the Euler scheme.

For each simulated random walk (X j , 0 ≤ j ≤ NP ), we compute the score:

R i = g D (X NP ) + ∆t • NP-1 ∑ j=0 f ( X j ) -∆ s • NP-1 ∑ j=0 g τ V ( X j ) 1 Γ h V ( X j ) -∆W • NP-1 ∑ j=0 g ν V (X j ) 1 Γ h V ( X j )
where ∆t is the time of a step in G, ∆ s is the local time of a step along τ on Γ h V , and ∆W is the step increment along ν.

Then this procedure is iterated NT times in order to obtain a sample of R, and an estimation of its arithmetic mean u(x).

The integer NT corresponds to the number of simulated random walks, and consequently to the number of absorptions by the Dirichlet boundary part Γ D , more precisely by Γ h D . This integer NT is the maximum index of the main loop of the associated program. The second (and latter) loop follows, step by step, each random walk by adding up successively the values of the source function f , and the values of the Wentzell condition g V by diffusion and the values of the Neumann condition by straight reflection, then the value of the absorption effect g D . Finally the score of each corresponding R i is added up.

For each realization of the random variable R, we generate:

1. a first random walk in the open set G, depending on the Wiener process W t (see, for instance, the simplest random walk in [START_REF] Souza De Cursi | Numerical methods for linear boundary value problems based on Feyman-Kac representations[END_REF]), 2. a second random walk on the boundary part Γ V , depending on the Wiener process β t which is simulated as the process W t extended on the boundary layer Γ h V , and 3. a reflection by Γ V along the inward normal ν.

If h is the fixed average of stepsize of a simulated bidimensional walk X j , the time and space increments are defined by:

∆t = h 2 a G in G, ∆ s = h 2 2a V on Γ h V , ∆W = h inward along ν.
If a V ≡ 0, then the random walk is instantaneously reflected by the boundary Γ V as in [START_REF] Morillon | Numerical solutions of linear mixed boundary value problems using stochastic representations[END_REF] for the Neumann BVP. Else a V ≡ 0, we have a diffusion on Γ V and a mirror effect: the random walk diffuses on Γ V and resumes to the interior; this procedure goes on until the trajectory reaches the Dirichlet boundary part Γ D .

Remark 1 When X t hits the boundary part Γ V by means of the Wiener process W t , then (W • τ)τ is replaced by β τ and (W • ν)ν by |W • ν|ν, i.e. if X t ∈ Γ V , then

X t+dt = X t + √ a V τ(X t ) dβ ξ t + √ a V | dW t • ν | ν
where • is the scalar product in the euclidean space R 2 and:

| dW t • ν | = sgn [ (X t + dW t ) • ν ] ( dW t • ν )
(see [START_REF] Freidlin | Functional Integration and Partial Differential Equations[END_REF], page 85 for the relation between the positive process | x +W t | in R + and the local time ξ t at zero).

Remark 2 In our algorithm (see the pseudo-code shown below, page 13), we can use antithetic random variables with the same distribution to estimate u and improve its variance.

Two processes W

(1) t and W

(2) t are simulated by the two sequences of random vectors W 0 , W

1 , . . . , W

j and W 0 , W (2) 1 , . . . , W (1)
j defined by means of the recurrent formulae:

W 0 = x ∈ G, W (1)
j+1 = W (1) j + hD
j = -D (1) j in G
where r 1 and r 2 are independent normally distributed pseudorandom numbers, and:

D (1) j = r 3 τ + |r 4 | ν and D (2) j = -r 3 τ + |r 4 | ν on Γ V
where r 3 and r 4 are independent normally distributed pseudorandom numbers, independent of r 1 and r 2 .

When a process W (i)

t is absorbed by Γ D , the second W ( j) t , j = i, continues until it also reaches Γ D .

For each simulated random walk (X j ; 0 ≤ j ≤ NP ), we compute the score:

R i = g D ( X NP ) + h 2 a G NP-1 ∑ j=0 f (X j ) 1 G ( X j ) - h 2 2a V NP-1 ∑ j=0 g τ V (X j ) 1 Γ h V ( X j ) -h NP-1 ∑ j=0 g ν V ( X j ) 1 Γ h V ( X j )
We then estimate u(x) by computing the arithmetic mean of these NT scores R i .

A pseudo-code description of the global algoritm can be written as follows: Declarations -Initializations, diffusion coefficient a V on Wentzell boundary part Γ V included Loop Points of domain G grid Set initial point x of the domain grid Main Loop i = 1, . . ., NT random walks from x to Dirichlet boundary part Γ D Let j = NT + i for the antithetic random walks Initialize for each random walk: Start points of antithetic random walks: W i = x; W j = x Source function f effect in domain G: ya i = 0.0 and ya j = 0.0 Wentzell diffusion effect on Γ V : maV i = 0.0 and maV j = 0.0 Neumann effect on Γ V : maN i = 0.0 and maN j = 0.0 While the random walks W i and W j within the domain G (Loop step by step)

Source function f effect at points W i and W j :

ya i = ya i + f (W i ) and ya j = ya j + f (W j )
Set D i a normally distributed pseudo-random vector Set r a normally distributed pseudo-random number

If W i on Wentzell boundary Γ V W i = W i + hrτ Diffusion effect: maV i = maV i + g V (W i )
Reflection of random walk inside G with Neumann effect:

maN i = maN i + g N (W i ) endIf If W j on Wentzell boundary Γ V W j = W j -hrτ Diffusion effect: maV j = maV j + g V (W j )
Reflection of random walk inside G with Neumann effect: maN j = maN j + g N (W j ) endIf W i = W i + hD i , i.e. a step of first random walk W j = W j -hD i , i.e. a step of second random walk endWhile the random walks W i and W j within the domain G (Loop step by step) Comment: Dirichlet condition: absorption of

W i or W j by Γ D If W i ∈ Γ D maU(x, i) = g D (W i ) -0.5h 2 /a V maV i -hmaN i endIf If W j ∈ Γ D maU(x, j) = g D (W j ) -0.5h 2 /a V maV j -hmaN j endIf Comment: Continue with the unabsorbed walk W i or W j If W j ∈ Γ D While W i / ∈ Γ D ya i = ya i + f (W i )
Set D i and r normally distributed pseudo-random vector and number

If W i on Wentzell boundary Γ V W i = W i + hrτ Diffusion effect: maV i = maV i + g V (W i ) Reflection of random walk inside G with Neumann effect: maN i = maN i + g N (W i ) ; endIf W i = W i + hD i A step of first random walk endWhile W i / ∈ Γ D If W i ∈ Γ D maU(x, i) = g D (W i ) -0.5h 2 /a V maV i -hmaN i endIf endIf W j ∈ Γ D If W i ∈ Γ D While W j / ∈ Γ D ya j = ya j + f (W j )
Set D i and r normally distributed pseudo-random vector and number If W j on Wentzell boundary Γ V W j = W j + hrτ Diffusion effect: maV j = maV j + g V (W j ) Reflection of random walk inside G with Neumann effect:

maN j = maN j + g N (W j ) ; endIf W j = W j + hD j A step of second random walk endWhile W j / ∈ Γ D If W j ∈ Γ D maU(x, j) = g D (W j ) -0.5h 2 /a V maV j -hmaN j endIf endIf W i ∈ Γ D
Add source effect and boundary conditions: 1) maU(x, i) = maU(x, i) + h 2 /a G ya i for the first random walk W i 2) maU(x, j) = maU(x, j) + h 2 /a G ya j for the antithetic random walk W j end Main Loop i = 1, . . ., NT random walks from x to Dirichlet boundary part Γ D Calculate estimated solution (arithmetic mean): The validity and suitability of the present approach, i.e. the stochastic representation in Section 2 and the approximate solution in Section 3, are investigated here by considering and evaluating a set of BVPs with Wentzell's boundary condition according to the above numerical scheme used sequential calculus on personal computers (Intel i7 CPU 2.67 GHz, and Athlon 2850e 1.80 GHz) by the means of MATLAB R applications.

u(x) = 1 2NT ∑ 2NT k=1 maU(x, k) U(x) = mean{maU(x, k) | 1 ≤ k ≤ 2NT }
In all the following examples, the solution u is approximated by

u c in NPT points Q 1 , Q 2 , . . . , Q NPT in G.
The quality of the stochastic approximation is controlled as follows:

ER max = max { ER i : 1 ≤ i ≤ NPT } ER i = 10 | u(Q i ) -u c (Q i ) | 1 + 10 | u(Q i ) |
so that ER i , i = 1, . . . , NPT , furnish relative errors or absolute errors, according to the magnitude of the theoretical value.

The following examples show BVP [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], one defined in a rectangular domain G = ]0, a[ × ]0, b[, and two in the unit square G = ]0, 1[ 2 . The fourth example shows BVP (4) defined in quadrant of ring G = D(0; 3) \ D(0; 1) ∩ (R * + × R * + ). For sake of simplicity, the diffusion coefficient a G in G equals to 1.

Wentzell-Dirichlet mixed BVP in a rectangular domain

Let us consider problem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] in G = ]0, a[ × ]0, b[ with the following data:

f ≡ 0 in G Γ D = ([0, a] × {0}) ∪ ([0, a] × {b}) g D (x, 0) = 0 on [0, a] × {0} g D (x, b) = cos 2 (πx/a) on [0, a] × {b} Γ V = ({0} × ]0, b[) ∪ ({a} × ]0, b[) g τ V (x, y) = a V π 2 a 2 sinh(2πy/a) sinh(2πb/a) on Γ V g ν V (x, y) = 0 on Γ V
such that this Wentzell-Dirichlet BVP admits a unique solution u written in the form:

u(x, y) = y 2b + cos(2πx/a) sinh(2πy/a) 2 sinh(2πb/a) .
Let us first consider the point (a/2, b/2) and the solution u(a/2, b/2) for different values of diffusion coefficient a V on Γ V : a V ∈ {1, . . ., 20}.

The results of the mean value of the approximate solution u(a/2, b/2) are obtained with the following data: The arithmetic mean of these 2 NT scores R i gives: u(a/2, b/2) 0.2455, the standard deviation σ = 0.1775 is obtained by means of antithetic variables, and 95% confidence intervals u c -1.96σ/ √ 2NT , u c + 1.96σ/ √ 2NT , around u c = 0.2455, is of amplitude 0.028.

The empirical error for the Wentzell diffusion coefficient a V :

ER a V (a/2, b/2) = 10 | u(a/2, b/2) -u c (a/2, b/2) | 1 + 10 | u(a/2, b/2) |
is less than 2 × 10 -3 and verifies:

max 1≤a V ≤20
ER a V (a/2, b/2) 1.9 × 10 -3 and 1 20 20

∑ a V =1 ER a V (a/2, b/2) 0.8 × 10 -3 .
See figure 1 where numerical solutions and empirical errors are represented versus a V . The effect of Dirichlet condition is isolated, so we can appreciate the effect of diffusion on Γ V . Let us secondly consider the points Q i j = Q(x i , y j ), where

x i = 0.1 × a × i, 1 ≤ i ≤ 9, y j = 0.1 × b × j, 1 ≤ j ≤ 9,
the approximated solution u c (Q i j ), the empirical errors ER i j at the point Q(x i , y j ), and ER max when the boundary diffusion coefficient a V belongs to {5, 10, 15, 20}.

The values of ER max given in Table 1 are obtained with the number of simulated random walks: NT = 2 × 10 5 and the data above, and the empirical errors ER i j are given in Figure 2.

Γ D = ([0, 1] × {0}) ∪ ([0, 1] × {1}) and Γ V = ({0} × ]0, 1[) ∪ ({1} × ]0, 1[) ,
and with the following data:

First problem P 1 in the unit square:

f ≡ -2 in G g D (x, 0) = x 2 + 1 ∀x ∈ [0, 1] g D (x, 1) = x 2 ∀x ∈ [0, 1] g τ V ≡ a V on Γ V g ν V (0, y) = 0 ∀y ∈ ]0, 1[ g ν V (1, y) = -2 ∀y ∈ ]0, 1[
such that this first Wentzell-Dirichlet BVP P 1 admits a unique solution u written in the form:

u(x, y) = x 2 + (y -1) 2 .
Second problem P 2 in the unit square:

f ≡ 0 in G g D (x, 0) = cos(πx) ∀x ∈ [0, 1] g D (x, 1) = cos(πx) exp(π) ∀x ∈ [0, 1] g τ V (0, y) = 1 2 a V π 2 exp(πy) ∀y ∈ ]0, 1[ g τ V (1, y) = -1 2 a V π 2 exp(πy) ∀y ∈ ]0, 1[ g ν V ≡ 0 on Γ V
such that this second Wentzell-Dirichlet BVP P 2 admits a unique solution u: u(x, y) = cos(πx) exp(πy).

Let us consider the points Q i j = Q(x i , y j ), where x i = 0.1 × i, 1 ≤ i ≤ 9, y j = 0.1 × j, 1 ≤ j ≤ 9, the approximated solution u c (Q i j ), the empirical errors ER i j at the point Q(x i , y j ), and ER max when the boundary diffusion coefficient a V belongs to {5, 10, 15, 20}.

The value of ER max and the mean of ER given in Table 2 for problem P 1 and in Table 3 for problem P 2 are obtained with the following data:

1. Number of simulated random walks: NT = 6 ×10 4 , i.e. 2NT = 12 ×10 4 by means of antithetic variables 2. Stepsize of each random walk: h = 0.01 and the empirical errors ER i j are given in Figure 3 for problem P 1 and in Figure 4 for problem P 2 .

For the BVP P 2 , the maximum error ER max is obviously where the exact solution equals to zero (the results in Figure 4 where the zero points are clearly identified, confirm this effect due to the choice of the expression of ER max ); for these points, the maximum error is equal to 10 times the absolute error and the 95% confidence intervals given in Table 4 always contain the exact solution u = 0. Let us finally consider the BVP P 1 with a coefficient a V different on each part of Γ V ; the results are analogous as for instance:

a V = 4 on Γ 1 V = ]0, 1[ × {0} 8 on Γ 2 V = ]0, 1[ × {1}
with the following data: The empirical errors ER i j for a V non constant are given in Figure 5.

Γ V = (]1, 3[ × {0}) ∪ ({0} × ]1, 3[) g τ V ≡ a V on Γ V g ν V ≡ 0 on Γ V
such that this Wentzell-Dirichlet BVP admits a unique solution u written in the form: u(x, y) = x 2 + y 2 . Let us consider the points Q i j = Q(x i , y j ) ∈ G, where x i = 0.2 × i, 0 ≤ i ≤ 15, y j = 0.2 × j, 0 ≤ j ≤ 15, and 1 < x 2 i + y 2 j < 9, i.e. a total number of 163 points in G. The estimated solution u c (x i , y j ) and empirical error ER i j when the boundary diffusion coefficient a V equals to 5, are obtained with the following data:

1. Number of simulated random walks: NT = 10 4 , i.e. 2NT = 2 × 10 4 by means of antithetic variables 2. Stepsize of each random walk: h = 0.01 and are given in Figure 6.

Remark 3 When the number of realizations is NT = 10 4 , the average CPU time used by the MATLAB R application implemented with a processor Intel i7 2.67 GHz (6.9 operations per second) for the numerical resolution of the BVP P 2 in the unit square, equals to 160 seconds when the results are displayed at each point, and 97 second per point without displaying.

The results can be improved by using more realizations and smaller stepsize h, but this increases CPU time involved. Let us note that the reflections increase the CPU time involved, in particular when the initial point of the random walk is near one part of the boundary Γ V .

Concluding remarks

In the present paper, we have presented stochastic method for the numerical approximation of some boundary value problems with Wentzell's boundary part and superabundant data on this part. The representations of solution are deduced from the application of Itô's formula even when the boundary data are non smooth. These probabilistic interpretations have given some easily implementable algorithms: they consist in the simulation of diffused-reflected random walks in order to generate samples from a distribution having a mean equal to the solution of the Wentzell BVP considered. This algorithm do not require to store in memory any mesh of the region G, nor deal with the corresponding numbering-lists. Programming is short, easy to check step by step.

As for the classical Monte Carlo schemes, this stochastic method admits an expected rate of convergence of about 1/ √ NT . Relative error and empirical variance can be computed. The variance reduction techniques as antithetic variates can be applied. The essential properties of Monte Carlo methods are maintained; on the one hand, this stochastic method of the solution approximation does not depend on the choice of the point, on the other hand, their implementations on a parallel computer are intrinsically possible.

The development of the stochastic methods by the means of an integral representation makes possible not only the solution of these systems but extensions to more general models with general linear elliptic or parabolic operators and with wide classes of boundary conditions.

jj

  are random vectors so that, if (e i , i = 1, 2 ) is the canonical basis of R 2 , D = r 1 e 1 + r 2 e 2 and D[START_REF] Bally | The law of the Euler scheme for stochastic differential equations. I. Convergence rate of the distribution function[END_REF]

  Display per point approximated solution, empirical errors, standard deviation and 95% confidence intervals endLoop Points of domain G grid 4 Numerical experiments

1 .Fig. 1

 11 Fig. 1 Numerical solution u(a/2, b/2) and empirical error versus a V

Fig. 2

 2 Fig. 2 Empirical errors in rectangular domain for a V = 5, 10, 15, 20

Fig. 3 BVP P 1 :

 31 Fig. 3 BVP P 1 : Empirical errors in the unit square when a V = 5, 10, 15, 20

Fig. 4 BVP P 2 :

 42 Fig. 4 BVP P 2 : Empirical errors in the unit square when a V = 5, 10, 15, 20

Fig. 5 2 .

 52 Fig.5Empirical error when a V is different on each part of Γ V

4. 3

 3 Wentzell-Dirichlet mixed BVP in quadrant of a ringLet us consider problem (4) in G = (x, y) ∈ R 2 | 1 < x 2 + y 2 < 9, x > 0, y > 0 with the following data:f ≡ -2 in G Γ D = (x, y) ∈ R 2 | x 2 + y 2 = 1, x >= 0, y >= 0 ∪ (x, y) ∈ R 2 | x 2 + y 2 = 9, x >= 0, y >= 0 g D ≡ 4 on (x, y) ∈ R 2 | x 2 + y 2 = 1, x >= 0, y >= 0 g D ≡ 9 on (x, y) ∈ R 2 | x 2 + y 2 = 9,x >= 0, y >= 0

Fig. 6

 6 Fig. 6 Exact and estimated solutions -Empirical error on quadrant of a ring

Table 1

 1 Maximum empirical errors for different a V a V ER max at point (x i , y i )

	5	0.0251 at (0.6; 2.4)
	10	0.0180 at (1.8; 2.7)
	15	0.0194 at (1.6; 2.4)
	20	0.0208 at (0.6; 2.7)
											a V =5																			a V =10
																				0.025																	0.018
																																					0.016
			0.03																	0.02		0.02															0.014
			0.025
			0.02																			0.015															0.012
																				0.015
		ER ij	0.015																		ij ER	0.01															0.01
			0.01
			0.005																	0.01		0.005															0.008
			0																			0															0.006
			2.7																			2.7
			2.4	2.1	1.8	1.5 y j	1.2	0.9	0.6	0.3	0.2	0.4	0.6	0.8	1.0 x i	1.2	1.4	1.6	1.8	0 0.005		2.4	2.1	1.8	y j 1.5	1.2	0.9	0.6	0.3	0.2	0.4	0.6	0.8	x i 1.0	1.2	1.4	1.6	1.8	0 0.002 0.004
											a V =15									x 10 -3										a V =20
																																					0.02
																				18
																																					0.018
																				16
			0.02																			0.025															0.016
			0.015																	14		0.02															0.014
																				12		0.015
		ER ij	0.01																	10	ER ij	0.01															0.012
			0.005																	8		0.005															0.01
																																					0.008
			0																			0
			2.7																	6		2.7															0.006
			2.4																			2.4
				2.1	1.8	1.5 y j	1.2	0.9	0.6	0.3	0.2	0.4	0.6	0.8	1.0 x i	1.2	1.4	1.6	1.8	2 4			2.1	1.8	y j 1.5	1.2	0.9	0.6	0.3	0.2	0.4	0.6	0.8	x i 1.0	1.2	1.4	1.6	1.8	0 0.002 0.004

Table 2

 2 BVP P 1 : Maximum empirical error and error mean for different a V

	a V	ER max							at point (x i , y i ) mean(ER)
	5	0.0248 at (0.2; 0.2)								0.0086
	10	0.0182 at (0.1; 0.1)								0.0055
	15	0.0336 at (0.1; 0.2)								0.0074
	20	0.0214 at (0.1; 0.3)								0.0076
																				a V =5																			a V =10
																																							0.018
																																							0.016
				0.025																													0.02			0.02			0.014
				0.02																																0.015			0.012
		ER(x i ,y j )	0.01 0.015																													0.015	) j ,y i ER(x	0.01			0.01
				0.005																													0.01			0.005			0.008
				0																																0			0.006
				0.9	0.8	0.7	0.6	0.5 y j	0.4	0.3	0.2	0.1	0.1	0.2	0.3	0.4	0.5 x i	0.6	0.7	0.8	0.9	0.005			0.9	0.8	0.7	0.6	y j 0.5	0.4	0.3	0.2	0.1	0.1	0.2	0.3	0.4	x i 0.5	0.6	0.7	0.8	0.9	0 0.002 0.004
																				a V =15																			a V =20
																																							0.02
																																	0.03
																																							0.018
																																			0.025
			0.03																													0.025		0.02			0.016
																																							0.014
		ER(x i ,y j )	0.02																													0.02	ER(x ) j ,y i	0.01 0.015			0.012
			0.01																													0.015		0.005			0.01
				0																																0			0.008
				0.9																													0.01			0.9			0.006
				0.8	0.7	0.6	0.5 y j	0.4	0.3	0.2	0.1	0.1	0.2	0.3	0.4	0.5 x i	0.6	0.7	0.8	0.9	0 0.005			0.8	0.7	0.6	y j 0.5	0.4	0.3	0.2	0.1	0.1	0.2	0.3	0.4	x i 0.5	0.6	0.7	0.8	0.9	0 0.002 0.004

Table 3

 3 BVP P 2 : Maximum empirical error and error mean for different a V

	a V	ER max	at point (x i , y i ) mean(ER)
	5	0.3250 at (0.5; 0.1)	0.0497
	10	0.6922 at (0.5; 0.4)	0.0531
	15	0.6208 at (0.5; 0.6)	0.0569
	20	0.7145 at (0.5; 0.2)	0.0612

Table 4

 4 BVP P 2 : 95% confidence intervals when empirical error is maximum

	a V	95% confidence interval at point (x i , y i )
	5	[-0.025899 ; 0.09089]	at (0.5; 0.1)
	10	[-0.027408 ; 0.16584]	at (0.5; 0.4)
	15	[-0.011347 ; 0.13551]	at (0.5; 0.6)
	20	[-0.001648 ; 0.14454]	at (0.5; 0.2)
				a v =5
			0.7
			0.6
			0.5
		ER(xi,yj)	0.3 0.4
			0.2
			0.1
			0
			0.9
			0.8
			0.7
			0.6	0.5		0.7	0.8	0.9
				0.4		0.6
				0.3		0.5
				0.2	0.3	0.4
				0.1	0.2
				0.1
			y j			x i

  1. Number of simulated random walks: NT = 1 × 10 4 , i.e. 2NT = 2 × 10 4 by means of antithetic variables

								a V = 4 on Γ V 1 and a V =8 on Γ V 2
																		0.022
																		0.02
		0.025																0.018
		0.02																0.016
	ER(x i ,y j )	0.01 0.015																0.012 0.014
		0.005																0.01
		0																0.008
		0.9	0.8															0.006
			0.7	0.6	y j 0.5	0.4	0.3	0.2	0.1	0.1	0.2	0.3	0.4	x i 0.5	0.6	0.7	0.8	0.9	0.002 0.004

Acknowledgements I warmly thank the reviewers for their suggestions and comments.