Goujon Person
email: marjorie.person@univ-reunion.fr

SOLAR ENERGY RESEARCH AND DEVELOPMENT PROGRAM ON THE EXPLOITATION OF THE SOLAR RESOURCE ON THE REUNION ISLAND AND ITS INTEGRATION INTO AN ELECTRICAL POWER GRID

Reunion Island stands out by its specific energy context due to its growing population and a strong economic development. The main objective of this French region is to ensure that its exceptional potential for renewable power generation could meet an increasing energy demand. Regarding the energy development, the Regional strategy is to reach the energy independence by 2030 (SPL "Energies Réunion") without using fossil fuels. In this context, the LE 2 P is leading a solar resource research programme with the aim to use solar resource as a stable source of energy and ensure its management in a reliable and efficient way for its integration into an electrical power grid. This paper aims to describe the methodology used for the development of this ambitious and challenging R&D program, projects articulation within the LE 2 P roadmap's and the collaborations and partnerships initiated in UE and OI zone, permitting Eco city development especially in the IO area.

INTRODUCTION

Using solar resource as an electrical and stable source of energy address three scientific and technical points not still completely solved. The first point concerns the solar irradiation variability characterization and its prediction. Next point is to manage the solar PV production through storage systems with the objective to reduce the intermittent nature of the solar resource impact and thus to transform the solar resource into a stable source of energy. The final point is to dynamically manage this stable source of energy into a power grid so that it could able to meet the energy needs and to cope with the load fluctuations.

The particularities of the Reunion Island energy context of lies in the fact that as a part of Island areas not interconnected to the metropolitan power grid, this French region is dependant from fossil fuels while having a High potential in renewable energies. In this context, the French government and the region have set up the SPL "Energies Réunion" project, with objectives to reach the energetic autonomy by 2030 and to make Reunion Island an exceptional territory of all the environmental innovations interesting mobility, constructions, tourism, urban development and mainly energy and its practices. Another specific aspect is related to the Reunion Island socio economic conditions with high unemployment, a low income per capita, a growing population and increasing energy needs. Taking into account the whole Reunion context, the LE 2 P has defined a solar energy development roadmap structured around the three points mentioned above. The paper is composed of three main parts. In the first one, we will describe the projects with their specific objectives, their connections with the scientific and technical points mentioned and how they hang together within the LE 2 P Solar energy roadmap which comes along two lines. The first line is composed of two steps. The first step is the solar irradiation prediction and its management through storage systems in order to get a stable source of energy. The second step deals with the energy delivery management through an electrical power grid so that the tandem production consumption could be balanced, as does a smart grid. The second line aims to bring together several scientific communities around the analysis and the characterization of the solar radiations through the extension of our existing ground sensors network to other territories of the southwest Indian Ocean area. Specific objectives associated to this line are first to provide solar radiation measurements with a high degree of quality through a BSRN (Baseline Surface Radiation Network) application and secondly to develop specific data storage system like Data Warehouse able to dynamically manage an increasing volume of data. In the second one, we will presented the specific methodology along several R&D lines, requires reaching the LE 2 P roadmaps objectives.

The third part will focus on the micro-grid concept, seen as an alternative approach for integrating small scale distributed energy resources into low-voltage electricity system. Finally, in the last part we will discuss about R&D lines identified in the solar energy roadmap along two points; the scientific and technical outlooks that could be open up then the scientific collaborations involved and those will be initiated.

THE LE 2 P SOLAR ENERGY DEVELOPMENT ROADMAP

The solar energy development scientific research program undertaken by the LE 2 P since 2011 is based on the RCI_GS project (Intelligent Sensor Network for a better knowledge of the solar resource all over the Reunion Island) dedicated to the solar irradiation measurement and the characterization of its spatial and temporal variability. The approach used to measure the solar radiations (Global, Diffuse and Direct) was to develop a ground sensors network distributed all over the Reunion Island. The RCI_GS project has been completed successfully, and from the experience gained in solar metrology, sensor network running and statistical classification methods, the LE 2 P has defined a roadmap along two main steps (Figure 1). The initial step consists in turning an intermittent source of energy into a stable one through the solar resource characterization and prediction and the management of the resulted PV electrical production through storage units. The second step deals with the dynamic management of the energy delivery through a power grid. The concept that will be used for this part is the micro grid concept. A Micro grid can be defined as an integration platform for supply-side (micro-generators) and demand-side (storage units and controllable loads) located in a local distribution grid. This emerging concept appears to meet growing customer needs for electric power with an emphasis on reliability, power quality and contribution to different economic, environmental ant technical benefits. In order to complete all the objectives (yellow on figure 1), the roadmap has been mapped out in several scientific and technical stages (green on figure 1) each of these associated to different projects (in blue on the figure 1). The results of the RCI_GS project have led to characterize the solar irradiation variability on the Reunion Island through the development of a ground sensors network and classification methods. Class 2 Intermittent bad day: a sunny beginning until mid-morning around 9:00-9:30 AM and a cloudy afternoon.

Class 3 Disturbed days: corresponds to a variable weather day with a high variability of the solar direct fraction.

Class 4 intermittent good day: similar to the class 2 but with a stronger sunny regime during all morning till early afternoon.

Class 5 clear sky day: regime of good weather throughout the day. The main objective of the AXISoi project (Correlation analysis of the solar radiation in the Indian Ocean) is to test the robustness of the classification method developed for the RCI_GS project by applying it to South Africa. The results achieved through the RCI_GS and AXISoi projects enable the extension of the RCI_GS ground sensors network to other territories of the Indian Ocean southwest area, which is the aim of the GEO_SUN project with South Africa as initial partner. One objective besides the project is to bring together a large number of scientific laboratories around the solar radiation characterization and prediction in order to meet a solar radiation mapping across the Indian Ocean southwest area. To complete successfully the project and meet the objectives requires leading two R&D actions: -To ensure a high degree of quality to the solar radiations measurements by -To develop specific means able to store and dynamically manage an increasing volume of data as a Data Warehouse (DWH). Besides the DWH would be used as a decision and learning process support. -To add an R&D line dedicated to the satellite imagery with the objectives:

o To assist in the choice of the ground station locations to maximize the representativeness in term of local climate. o To fuse the satellite data with the ground stations data in order to meet a spatiotemporal cartography with a sufficiently fine resolution. The GEO_SUN project will result in the end to predict accurately the solar resource and its management through storage units. This project should enable the LE 2 P to complete a crucial step and will allow preparing the implementation of the future projects linked to the LE 2 P roadmap final step. One of them is the RIVAGES project (Reunion Vanilla Islands Solar Resource), which will consist to expand the GEO_SUN ground sensors network through the southwest Indian Ocean area. Another future project that prefigures the solar resource management through a power grid is named GYSOMATE, which means "solar resource dynamic management for the energy consumption control through a smart city network composed of production and storage units".

THE R&D PROGRAM AND METHODOLOGY

The methodology used to complete the R&D program is structured around three main areas of study: 1) the solar irradiation prediction, 2) the electrical solar energy production management through storage units 3) Energy delivery management into an urban micro grid.

1) The Solar irradiation prediction

The solar irradiation prediction requires first to measure solar radiations with enough accuracy and to characterize its variability both spatially and temporally [START_REF] Badosa | Scales of spatial and temporal variation of solar irradiance on Reunion tropical island[END_REF][START_REF] Bessafi | Classification des journées en fonction des radiations solaires sur l'île de la Réunion. 45 e journées de la Statistique[END_REF]. Two approaches are studied: the first has been used for the RCI_GS project and consist in developing classification methods to identify several classes of solar radiation timelines, another approach consists in investigating multifractality of Solar Irradiance time series through wavelet based multifractal spectral analysis. For the long-term solar irradiance prediction, the first approach base on classification methods is coupled with an artificial neural network technic. The architecture the most commonly used is the multilayer perceptron (MLP). One MLP consist of three or more layers (one input layer, one or more hidden layer and one output layer) in which each layer is fully connected to the next one. Except for the input nodes, each node is a neuron or a processing element with a non-linear activation function.

Figure 2 shows a one hidden layer MLP. Short-term solar irradiation prediction requires a different approach based on two points. The first point lies in the fact that solar irradiation fluctuations are within several time scales. The second point that has to be taking into account is that solar irradiation fluctuations do not exhibit a simple monofractal scaling behavior. There might exist some crossover time scales separating regimes with different scaling exponents. The approach used in this case is based on investigating the multifractality of the solar irradiance data through wavelet based multifractal spectral analysis. The qualitative measure of the degree of multifractality of the solar irradiance signal will help us to decide which signal processing tools that can be used to extract the features of the signal. ..) after conversion of basic electricity produced. This intermediate energy is subsequently converted into electrical energy according to the needs of the user. Among the storage means most commonly used, the storage of energy, hydrogen is distinguished by its ability to store a long time and in large quantities [START_REF] Benchrifa | Rôle de l'hydrogène dans le stockage de l'électricité à base des énergies renouvelables[END_REF]. Moreover it's transportable as natural gas.

Towards exploring the use of Hydrogen as an energy vector, two ways for clean production at a large scale have reached a technological maturity: hydrogen production by water electrolysis and by thermo chemical biomass conversion [START_REF] Ni | An overview of hydrogen production from biomass[END_REF].

Electrolysis, especially PEM electrolysis way is the most promising [START_REF] Barbir | PEM electrolysis for production of hydrogen from renewable energy sources[END_REF].

The excess electricity from solar energy over the amount used will be used to produce and store hydrogen in order to use it later as needed. Contrary to fuel cells, electrolysers produce hydrogen and oxygen by water electrolysis according to the chemical reaction:

H 2 O → H 2 +1/2 O 2
The electrolyse enthalpy is about 285 KJ/Mole at 298 K and 1 Bar. The theoretical potential threshold to ignite the decomposition process at 298K is 1, 48 V, using cells with potential in the range of 1, 7 to 2, 1 V.

Concerning water thermolysis, the thermal waters dissociation requires high temperatures (2500K), effective techniques to separate hydrogen and oxygen and prevent recombination. Another technical solution to store energy is the battery, which is an electrochemical storage system. It can be used to get back the chemical energy created by electrochemical reaction as electrical energy. These reactions are activated in a unit cell with two electrodes immersed in an electrolyte. The battery is based on a reversible electrochemical system, which makes it rechargeable.

3) Energy delivery management into an urban micro-grid.

The main problems that smart grids are meant to solve are those of smaller operating margins accompanied by a more volatile supply. Operating margin can be defined as the excess electrical power over the amount used at any moment. Volatile supply can be seen as the result of using sources whose outputs cannot be precisely predicted or controlled. Together these two issues lead to an oversupply or an under supply at any given moment.

Actual specifications tend toward an energy consumption reduction or to meet a particular load curve while balancing two factors, energy surety and economics.

With these objectives, micro grid is an emergent concept providing numerous advantages relatively to insular contexts with areas little or not connected to an electrical grid distribution. A micro grid is defined as a small grid that can operate as a part of a larger grid or that can operate independently of the larger grid.

Rolling out micro grids offer several advantages:

-On a technical side, micro grid allows optimized management of the electrical production from renewable energies at a local scale. It helps to maintain the voltage stability of the public supply system and reduces its load when disconnected from the micro grid.

-

On an economical side, the energy production and consumption closeness allows to optimize the energy transport. In addition, micro grids result in an increase of the flexibility and the quality of energy loads. -On a societal side, the micro grid increases the reliability and the stability of the power grid and provides responses to the energy needs changing of a territory.

-On an environmental side, the micro grid concept is the solution to integrate renewable sources of energy in an electrical distribution network.

MICRO-GRIDS:

A NEW PARADIGM FOR ELECTRICITY DELIVERY.

CHALLENGES AND APPLICATIONS.

Micro-grids are much more than backup generation for residential customers. A micro grid is "a group of interconnected loads and Distributed Energy Resources (DERs) with clearly defined electrical boundaries that acts as a single controllable entity with respect to the grid. DERs consist of distributed generation and Distributed Energy Storage (DES) installed at utility facilities, e.g., distribution substations, distributed generation sites or consumer premises. Micro grids may be operated in two modes: 1) Interconnected to the grid: The micro grid can import, export or have zero power exchange with the grid. This type of operation is designed as normal conditions. The objectives are to improve the grid performance and efficiency by using local DERs in order to defer capacity investments, reduce system losses and improve local reliability.

2) Disconnected from the grid: The micro grid is allowed to operate islanded mode. This operating mode requires the DERs within the micro-grid to be dispatched in a coordinated fashion to provide voltage and frequency regulation. Successful islanded operations may entail the implementation of energy demand management to achieve production load balance. The micro grid is expected to return to interconnected mode after the contingency has been addressed. A conceptual illustration of a micro rid within the context of a smart grid is shown in Figure 3.

Figure 2: A conceptual illustration of a micro grid Electric power can be generated at a distribution level in a micro grid. It usually includes a variety of small power generating sources, as well as energy storage systems such as batteries, flywheels, and super-capacitors. The power generating sources may include renewable sources such as solar panels and wind turbines, which are typically located close to the consumer sites. A micro grid can be coupled with the utility power grid through a single connection, known as point of common coupling (PCC). The electrical energy can flow in either direction through this coupling, based on the available energy generated within the micro grid and the demands of the consumers within the micro grid. A micro-grid, when disconnected from the main grid, is known as an "islanded micro grid." In an islanded micro grid operation, distributed generations (DGs) continue to power the users of the micro grid without requiring obtaining electric power from the utility grid. The micro grid connection and disconnection operations are decided by the PCC.

CONCLUSION AND PERSPECTIVES

This paper aims at highlighting the Reunion regional strategy regarding the solar energy sector through the LE 2 P R&D program including several inter-linked projects drawing up a roadmap regarding the solar energy exploitation and its use as a stable electrical source of energy. The Laboratory expertise is clearly presented by describing all the future and current projects leaded in the Laboratory. We emphasize about the key role played by the GEO_SUN project conducted in partnership with UKZN (University of Kwazulu Natal, South Africa) both on a scientific and technical side and in terms of influence on the regional, national and international scientific community. With this in mind, the LE 2 P has decided to initiate an important action of valorisation and communication showing our willingness to implement a long-term collaboration with South Africa as well as the southwest Indian Ocean Islands through the planed future projects RIVAGES and GYSOMATE.

Figure 1 LE 2 P

 12 Figure 1 LE 2 P Solar Energy Development Roadmap These classifications methods have led to identify several classes of solar radiation timelines, covering all the sunshine profiles existing on the Island [1, 2]; Class 1 Cloudy day: It corresponds to a very low sunshine level all day.Class 2 Intermittent bad day: a sunny beginning until mid-morning around 9:00-9:30 AM and a cloudy afternoon.Class 3 Disturbed days: corresponds to a variable weather day with a high variability of the solar direct fraction.Class 4 intermittent good day: similar to the class 2 but with a stronger sunny regime during all morning till early afternoon.Class 5 clear sky day: regime of good weather throughout the day. The main objective of the AXISoi project (Correlation analysis of the solar radiation in the Indian Ocean) is to test the robustness of the classification method developed for the RCI_GS project by applying it to South Africa. The results achieved through the RCI_GS and AXISoi projects enable the extension of the RCI_GS ground sensors network to other territories of the Indian Ocean southwest area, which is the aim of the GEO_SUN project with South Africa as initial partner. One objective besides the project is to bring together a large number of scientific laboratories around the solar radiation characterization and prediction in order to meet a solar radiation mapping across the Indian Ocean southwest area. To complete successfully the project and meet the objectives requires leading two R&D actions: -To ensure a high degree of quality to the solar radiations measurements by -To develop specific means able to store and dynamically manage an increasing volume of data as a Data Warehouse (DWH). Besides the DWH would be used as a decision and learning process support. -To add an R&D line dedicated to the satellite imagery with the objectives:

Figure 2 :

 2 Figure 2: Multilayer Perceptron architecture

ACKNOWLEDGEMENT

This work is supported by Regional Reunion, Island Council and the French government through the ERDF (European Regional Development Fund).


C Chauvière
email: cedric.chauviere@uca.fr

H Djellout
email: hacene.djellout@uca.fr

An efficient spectral method for the numerical solution to stochastic differential equations

Keywords: Stochastic differential equations, Numerical approximation, Spectral expansion, Stochastic collocation

We consider a new approach for the numerical approximation of stochastic differential equations driven by white noise. The proposed method shares some features with the stochastic collocation techniques and, in particular, it takes advantage of the assumption of smoothness of the functional to be approximated, to achieve fast convergence. The solution to the stochastic differential equation is represented by means of Lagrange polynomials. The coefficients of the polynomial basis are functions of time and they can be computed by solving a system of deterministic ordinary differential equations. Numerical examples are presented to illustrate the accuracy and the efficiency of the proposed method.

Introduction

Stochastic differential equations (SDE) driven by a white noise are present in a wide variety of areas. Historically, they were first used to model the motion of Brownian particles at a microscopic level and a link between microscopic and macroscopic quantities could be made by Einstein and Langevin at the beginning of the 20th century. Later, the same ideas were used to model various physical phenomena. For example, non-Newtonian fluid models have been designed by considering the action of white noise to account for the effect of the thermal agitation on polymer chains [START_REF] Öttinger | Stochastic processes in polymeric fluids[END_REF]. Then, starting in the mid seventies, there has been an upsurge of interest in the study of SDE by the mathematical finance community after the work of Black and Scholes [START_REF] Black | The pricing of options and corporate liabilities[END_REF]. Another important aspect of stochastic problems lies in the study of the propagation of white noise into a physical system described by partial differential equations [START_REF] Hou | Wiener chaos expansions and numerical solutions of randomly forced equations of fluid mechanics[END_REF][START_REF] Zhang | Numerical solution of the Stratonovich-and Ito-Euler equations: application to the stochastic piston problem[END_REF]. All these topics have in common a need for fast and accurate numerical methods to solve their underlying SDE.

One of the most common scheme of approximation is the Euler-Maruyama method and its closest cousin, the Milstein scheme [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. Both are limited to a weak order of convergence (in time) equal to 1. More important is the strong order of convergence, which is also equal to 1 for Milsten but drops to 1/2 for Euler-Maruyama. Unfortunately, designing general numerical schemes with higher order of convergence is not quite as easy as it is for its deterministic Runge-Kutta counterpart [START_REF] Burrage | Comment on "numerical methods for stochastic differential equations[END_REF] and different strategies have to be thought of. One of them, called "exact sampling" consists in generating realizations of the path for the SDE [START_REF] Beskos | Exact simulation of diffusions[END_REF][START_REF] Beskos | Retrospective exact simulation of diffusion sample paths with applications[END_REF]. This method is based on an acceptance/rejection sampling algorithm. It is "exact" in the sense that the marginal distribution of the simulated values coincides with the marginal distribution of the continuoustime process at grid points. Nevertheless one should bear in mind that all the methods quoted above are plagued with a slow rate of convergence typical to Monte-Carlo simulations (proportional to 1/ √ K , where K is the number of realizations of the sample paths). Although antithetic sampling or control variable techniques may be used to improve the convergence, such control variable approaches usually remain an ad hoc art [START_REF] Iacus | Simulation and inference for stochastic differential equations[END_REF].

In order to remedy this shortcoming, there has recently been increasing interest in considering the stochastic collocation method, or its modal counterpart the Polynomial Chaos expansion, to compute approximate solutions of SDEs. Since the work of Wiener in 1938 [START_REF] Wiener | The Homogeneous Chaos[END_REF] and the seminal book of Ghanem and Spanos in the 90's [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF], Polynomial Chaos (PC) based on orthonormal polynomials have been very successful at solving a wide variety of random problems involving random variables. Later, its analogous based on Lagrange polynomials was popularized by the paper of Hesthaven and Xiu [START_REF] Xiu | High-order collocation methods for differential equations with random inputs[END_REF]. Quickly these two techniques became the methods of choices for stochastic problems involving a limited number of random variables. The interested reader might refer to the book of Xiu [START_REF] Xiu | Numerical methods for stochastic computations[END_REF] for an overview of the possibilities offered by the PC and the stochastic collocation method. Building on the success of these polynomialbased approximations, numerical experiments were carried out in an attempt to apply the same ideas in the context of SDE driven by white noise. The least we can say is that the transition from random variables to stochastic process is not just a simple adaptation to existing techniques. The two main reasons are the infinite dimensional character of random process, leading to the so-called "curse of dimensionality" and the imperative need that the marginal distribution of the continuous-time process be properly taken into account. When this approach is considered, the Brownian motion is represented as a truncated Karhunen-Loeve series expansion and the stochastic differential equation reduces to a random differential equation [START_REF] Hou | Wiener chaos expansions and numerical solutions of randomly forced equations of fluid mechanics[END_REF][START_REF] Zhang | Numerical solution of the Stratonovich-and Ito-Euler equations: application to the stochastic piston problem[END_REF] involving a finite number of independent random variables. It has been shown by Wong and Zakai [22] that such approximation converges to the Stratanovich SDE and therefore, Itô SDE should be first converted into its corresponding Stratanovich equation. However, numerical experiments have shown that such a procedure failed to handle long-time simulations and treating non-linearities may be a problem. Discussion of these issues and possible remedies may be found in the recent book of Zhang and Karniadakis [START_REF] Zhang | Numerical methods for stochastic partial differential equations with white noise[END_REF].

Here, we propose a radically new approach that shares some features with the stochastic collocation method (SCM). In its simplest form, the SCM was successful in representing a random variable (r.v.) Y of unknown distribution (not necessary Gaussian) as a function of a r.v. X of known distributions (Gaussian, uniform, exponential...) under the form

Y = n i=1 y i L i (X), (1)
where {L i (x)} 1≤i≤n form a Lagrange polynomial basis and {y i } 1≤i≤n are n real numbers. If a functional relationship of the form Y = f (X) exists, then the coefficients y i are simply f (x i ) and equation ( 1) is the projection of f onto a Lagrange polynomial basis. If f is sufficiently regular and the collocations points {x i } 1≤i≤n are carefully chosen, one can observe spectral convergence to the exact solution [START_REF] Gottlieb | Numerical analysis of spectral methods: theory and applications[END_REF][START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF]. In this paper, following the same idea, we seek to express an unknown process (Y t ) t∈[0,T ] (not necessary Gaussian) as a function of a known process (G t ) t∈[0,T ] (Wiener, Ornstein-Uhlenbeck...) under the form

Y t = n i=1 y i (t)L i (G t ), for t ∈ [0, T ],
where {L i (x)} 1≤i≤n still form a Lagrange polynomial basis but now, {y i (t)} 1≤i≤n are n deterministic functions of time. Similarly to the SCM, macroscopic quantities (mean, variance, autocorrelation function...) can easily be computed from the coefficients of such representation as long as the probability density function of G t is known. In doing so, there is no need to discretize the process under the form of a finite number of independent random variables and the curse of dimensionality is no longer a problem. Another benefit of this approach is that it does not require a Monte-Carlo simulation to compute macroscopic quantities. We will use such representation for the numerical solution to SDE and the coefficients y i (t) will be shown to be the solution of a system of ordinary differential equations. In this paper, we assume that (G t ) t∈[0,T ] is a Gaussian process of the form

G t = t 0 k(s)dW s , (2)
where (W s ) s∈[0,T ] is a Wiener process and k : [0, T ] → R is a real function computed in tandem with y i (t). This paper is organized as follows. In section 2, we describe the basic ideas of the stochastic collocation method when applied to random variables problems. In an effort to make the paper self-contained, we also restate some classical results of polynomial approximation theory and quadrature rules. Then, we explain how the stochastic collocation method can be adapted in the context of random processes. Section 3 describes the type of SDEs addressed in this paper.

Based on Itô's lemma, in section 4 we derive an ordinary differential equation for the functions y i (t). Section 5 offers simple worked out examples when the coefficients y i (t) and the function k(t) can be computed explicitly. We then move on to more complex numerical examples, followed by some conclusions.

The stochastic collocation method (SCM)

Polynomial approximation and quadrature rules

The two main ingredients of the SCM (interpolating polynomials and quadrature rule) are borrowed from the deterministic community and we start to recall some of their main features. Let P n denote the linear space of polynomials of degree at most n. Given (n + 1) distinct points x 0 , x 1 , ..., x n and (n + 1) corresponding values y 0 , y 1 , ..., y n , there exists a unique polynomial Q ∈ P n such that Q(x i ) = y i . Lagrange provides an elegant way of constructing such polynomials under the form

Q(x) = n i=0 y i L i (x),
where

L i (x) = n j=0 j =i x -x j x i -x j , i = 0, ..., n (3)
form a polynomial basis of degree n.

Let f ∈ C n+1 (I) such that f (x i ) = y i and I is the smallest interval containing all the points x i . Following standard results in numerical analysis, the interpolating error satisfies

∀x ∈ I, f (x) -Q(x) = f (n+1) (ξ x ) (n + 1)! n+1 i=0 (x -x i ),
where ξ x ∈ I. From this equality, we see that the quality of the approximation not only depends on the regularity of f but also on the choice of the points x i . This can be made more precise using the following theorem, see [START_REF] Gottlieb | Numerical analysis of spectral methods: theory and applications[END_REF][START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF].

Theorem 2.1. Assume that f is a continuous function on the interval I = [-1, 1] and Q is a n th-order interpolating polynomial of f based on grid points

x i , then f -Q ∞ ≤ (1 + Λ n ) f -Q * ∞
, where Q * is the best n th-order interpolating polynomial of f and

Λ n = max x∈[-1,1] n i=0 |L i (x)| is called the Lebesgue constant.
This theorem shows that for a continuous function, the Lebesgue constant, and thus the quality of the approximation, only depends on the choice of the points x i . Unfortunately, designing optimal set of points is not a trivial task and for numerical simulations, one is often led to choose points that correspond to quadrature rules since they usually exhibit good approximating properties. We recall that a quadrature rule approximates an integral as a weighted sum as follows

I f (x)p(x)dx n i=0 ω i f (x i ), (4)
where ω i are the quadrature weight and x i are the quadrature points. The positive function p typically corresponds to a probability density of support

I. For example, if p(x) = 1 √ 2π
e -x 2 /2 and I = R, then the corresponding cubature rule is called Gauss-Hermite quadrature [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF]. It can be shown that Gauss quadrature rules (4) are exact when f is a polynomial of degree less or equal to (2n -1). For the construction of general cubature rules, one can consult [START_REF] Fejér | Mechanische Quadraturen mit positiven Cotesschen Zahlen[END_REF][START_REF] Dahlquist | Numerical methods in scientific computing[END_REF][START_REF] Davis | Methods of numerical integration[END_REF][START_REF] Rahman | Extended Polynomial Dimensional Decomposition for Arbitrary Probability Distributions[END_REF], for example. Equipped with those two ingredients, we can now move on to the stochastic collocation method.

The SCM applied to random variables

Denote by (Ω, F, P) the probability space, where as usual Ω is the set of all possible outcomes, F is a σ-algebra over Ω, and P is a function F → [0, 1] that gives a probability measure on F. Consider an R-valued random variable X of support I ⊂ R that describes input uncertainties. We assume that the probability law of X is known and that it is defined by a probability density function p(x), x ∈ I. Let Y = f (X) be a random model described by a function f : I → R. The beauty of the SCM resides in the interpolation points that are chosen to coincide with the quadrature points of (4). Given the set of collocation points {x i } 0≤i≤n , the function f can be evaluated at those points and approximated as

f (x) n i=0 f (x i )L i (x).
The quality of the approximation solely depends on the regularity of f and the choice of the collocation points. It is independent of the nature of the input argument x, should it be a real variable, a random variable or a stochastic process. In the case of random variables, we can write the approximation

Y = f (X) n i=0 f (x i )L i (X). ( 5
)
Macroscopic quantities such as the mean or the variance can be computed directly from (5) without having to resort to Monte-Carlo simulations. In addition, using the property of Lagrange polynomials L i (x j ) = δ ij (δ ij is the Kronecker delta) together with ( 4) and ( 5), it is easy to show that the formulas for the mean and the variance simplify to

E(Y ) n i=0 ω i f (x i ) (6)
and

var(Y ) n i=0 ω i f (x i ) 2 -E(Y ) 2 (7)
In the next subsection, we will see how these ideas can be applied in a similar way in the context of stochastic processes.

The SCM applied to stochastic processes

Let (G t ) t∈[0,T ] be a Gaussian stochastic process defined on a filtered probability space (Ω, F, P, (F t ) t∈[0,T ] ). We further assume that for all fixed t ∈ [0, T ], G t is a zero mean r.v. with variance σ 2 (t). Let f (t, x) : [0, T ] × R → R be a real function and (Y t ) t∈[0,T ] be a stochastic process (not necessarily Gaussian) defined by the functional relationship

Y t = f (t, G t ), t ∈ [0, T ].
Applying the same reasoning as before, for any fixed instant of time t ∈ [0, T ], Y t can be approximated as

Y t = f (t, G t ) n i=0 f i (t)L i (G t ), (8)
where the Lagrange polynomials are defined by (3) together with (n+1) distinct points {x i } 0≤i≤n . It can be easily seen that the coefficients f i (t) are obtained by the evaluation of f (t, x i ). The probability density function of G t writes

p Gt (t, x) = 1 √ 2πσ(t) exp - x 2 2σ(t) 2
and therefore the computation of the mean and the variance of Y t will require quadrature rules for integrals of the form R g(t, x)p Gt (t, x)dx, for some real functions g : [0, T ] × R → R. Denoting by (ω i , x i ) 0≤i≤n the Gauss-Hermite quadrature points and weights associated to a standard normal distribution, then a simple change of variable reveals that (ω i , σ(t)x i ) 0≤i≤n are the points and weights associated to the density p Gt (t, x). The most natural way to transcribe the SCM to random process would then be to consider timedependent Lagrange polynomials of the form

L i (t, x) = n j=0 j =i x -σ(t)x j σ(t)x i -σ(t)x j , i = 0, ..., n.
By so doing, we would obtain formulas similar to ( 6)-( 7) for the mean and the variance of

Y t i.e. E(Y t ) n i=0 ω i f (t, σ(t)x i ),
and

var(Y t ) n i=0 ω i f (t, σ(t)x i ) 2 -E(Y t ) 2 .
Although this choice of collocation point seems the most natural one, numerical results for the resolution of SDE using such representation were disappointing in terms of convergence and stability. This is because we need a set of interpolating points that are also capable of accurately representing derivatives of a function and not only the function itself. Accordingly, in this paper, we have used the Chebyshev nodes defined on the interval [a, b] by

x i = a + b 2 + b -a 2 cos 2i + 1 2(n + 1) π , i = 0, .., n (9)
to construct the Lagrange basis. Such basis offers good convergence properties since this set of points gives a near optimal Lebesgue constant [START_REF] Trefethen | Approximation theory and approximation practice[END_REF]. However, the formulas for the mean and the variance of Y t do not simplify and we give their new expressions in the sequel. Let ( ω i , x i ) 0≤i≤m denote cubature points and nodes for some general quadrature rule, for any fixed t ∈

[0, T ] R g(t, x)dx m i=0 ω i g(t, x i ). (10)
Then, assuming an approximation of Y t under the form (8), we have

E(Y t ) = R n k=0 f k (t)L k (x) p Gt (t, x)dx, E(Y t ) m i=0 n k=0 ω i f k (t)L k ( x i )p Gt (t, x i ), (11)
and similarly

var(Y t ) = R n k=0 f k (t)L k (x) 2 p Gt (t, x)dx -E(Y t ) 2 , var(Y t ) m i=0 ω i n k=0 f k (t)L k ( x i ) 2 p Gt (t, x i ) -E(Y t ) 2 . ( 12
)
For stochastic processes, the computation of the autocorrelation function E(Y t Y s ), (t, s) ∈ [0, T ] 2 is also of particular importance. This requires the knowledge of the joint density of the vector (G t , G s ) for t = s. It is given in the following lemma : Lemma 2.2. Let (G t ) t∈[0,T ] be a centered Gaussian process under the form (2). For fixed t, s into [0, T ] with s > t, the vector

(G t , G s ) is a centered Gaussian vec- tor of covariance matrix Γ(t, s) = σ(t) 2 σ(t) 2 σ(t) 2 σ(s) 2 , where σ(t) 2 = t 0 k 2 (u)du.
Furthermore, the joint density of the vector (G t , G s ) takes the form

p (Gt,Gs) (t, s, x, y) = 1 2π det(Γ(t, s)) e -1 2 det(Γ(t, s)) (σ(s) 2 x 2 +σ(t) 2 (y 2 -2xy)) , with det(Γ(t, s)) = σ(t) 2 (σ(s) 2 -σ(t) 2 ).
Proof 2.3. Let (G t ) t∈[0,T ] be the process given by (2). For s > t, we have

G t G s = 1 0 1 1 G t G s -G t = 1 0 1 1     t 0 k(u)dW u s t k(u)dW u    
The vector

    t 0 k(u)dW u s t k(u)dW u   
 is a centered Gaussian vector, since its component are independent centered Gaussian random variables. Therefore G t G s is a centered Gaussian vector as a linear transformation of a centered Gaussian vector.

Using Itô isometry we have

var(G t ) = E t 0 k(u)dW u 2 = t 0 k(u) 2 du := σ(t) 2 .
For s > t, using the independence of G t and G s -G t and the fact that E(G t ) = 0, we obtain

E(G t G s ) = E(G t (G s -G t )) + E(G 2 t ) = E(G 2 t ) = σ(t) 2 ,
which gives the components of the covariance matrix. The expression of the joint density arises directly from this result.

Using the approximation (8) and the quadrature rule [START_REF] Gottlieb | Numerical analysis of spectral methods: theory and applications[END_REF], for any fixed s and t into [0, T ], and assuming s > t, we have

E(Y t Y s ) = R 2 n k=0 f k (t)L k (x) n l=0 f l (s)L l (y) p Gt (t, s, x, y)dxdy E(Y t Y s ) m i=0 m j=0 n k=0 n l=0 ω i ω j f k (t)f l (s)L k ( x i )L l ( x j )p Gt,Gs (t, s, x i , x j ) (13)
In the numerical section, means, variances and autocorrelation functions of considered processes will be evaluated with formulas [START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF], [START_REF] Higham | An algorithmic introduction to numerical simulation of stochastic differential equations[END_REF], and (13).

Itô stochastic differential equation

In this paper, we are interested in the numerical solution to SDEs involving a random white noise term dW t that write under the standard form

dX t = a(t, X t )dt + b(t, X t )dW t , t > 0 X t=0 = X 0 , (14)
where, as usual, (W t ) t∈[0,T ] denotes Brownian motion on a filtered probability space (Ω, F, P, (F t ) t∈[0,T ] ), X 0 is F 0 -measurable and a(t, x) : [0, T ] × R → R and b(t, x) : [0, T ] × R → R are two deterministic functions. Equation ( 14) should be interpreted as its equivalent integral form

X t = X 0 + t 0 a(s, X s )ds + t 0 b(s, X s )dW s ,
where the first integral on the right-hand side is a Riemann integral and the second integral is considered as an Itô integral, see [START_REF] Oksendal | Stochastic differential equations[END_REF]. We assume that standard conditions are satisfied to ensure a unique strong solution to equation [START_REF] Iacus | Simulation and inference for stochastic differential equations[END_REF] i.e.

• The initial condition satisfies E(X 2 0 ) < ∞ and is independent of (W t ) t∈[0,T ]

• For all t ∈ [0, T ] and x, y ∈ R, a(t, x) and b(t, x) are continuous

• For all t ∈ [0, T ] and x, y ∈ R, a(t, x) and b(t, x) satisfy Lipschitz condition with respect to the second variable

|a(t, x) -a(t, y)| + |b(t, x) -b(t, y)| ≤ K |x -y|
In the numerical section, we will consider SDE of the form (14) that satisfy the above conditions.

Description of the proposed method

We now seek an approximate solution to [START_REF] Iacus | Simulation and inference for stochastic differential equations[END_REF] under the form

X t = f (t, G t ) = n i=0 f i (t)L i (G t ), (15)
where {L i (x)} 0≤i≤n are the Lagrange polynomials defined by (3) and f is some (unknown) smooth function. (G t ) t∈[0,T ] is a Gaussian process that we assume to be of the form

G t = t 0 k(s)dW s , (16)
where k : [0, T ] → R is a continuous real function such that k(0) = 0. The value of k(0) is unimportant and we take arbitrary k(0) = 1. In case k(0) = 1, we use the ratio k(s)/k(0) to define the stochastic process G t in ( 16) instead of k(s). Such representation covers a large variety of Gaussian processes: for example, the Wiener process for k(s) = 1, the Ornstein-Uhlenbeck process for k(s) = e θs with θ a real parameter, the Brownian bridge process for k(s) = α α -s (0 ≤ s < α) with α a real positive parameter... Using Itô isometry, the standard deviation σ(t) of G t can be shown to be

σ(t) = t 0 k(s) 2 ds. (17)
The unknown of the problem are the (n+1) coefficients f i (t) of equation ( 15) together with the function k(s) of equation [START_REF] Kloeden | Numerical solution of SDE through computer experiments[END_REF]. Once those quantities have been computed, the zero mean Gaussian process (G t ) t∈[0,T ] is fully specified, and the process (X t ) t∈[0,T ] is also completely represented by its approximation [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. It is then possible to compute macroscopic quantities with equations ( 11), [START_REF] Higham | An algorithmic introduction to numerical simulation of stochastic differential equations[END_REF] and [START_REF] Hou | Wiener chaos expansions and numerical solutions of randomly forced equations of fluid mechanics[END_REF]. The following theorem specifies the expression of f i (t) and k(t). Theorem 4.1. Let consider the Itô stochastic differential equation ( 14) and let assume that its exact solution can be expressed as an explicit function f (t, G t ) with (G t ) t∈[0,T ] a Gaussian process defined by ( 16) and f : [0, T ]×R → R a function whose second oder derivatives are continuous. Then f satisfies following system of partial differential equations for all (t, x) in [0, T ] × R :

         ∂f ∂t (t, x) = a(t, f (t, x)) - 1 2 ∂b(t, f (t, x)) ∂x × b(t, f (t, x)) ∂f ∂x (t, x)k(t) = b (t, f (t, x)) (18)
Proof 4.2. f being assumed to have second oder derivatives continuous, Itô lemma can be applied to obtain

dX t = df (t, G t ) = ∂f ∂t (t, G t )dt + ∂f ∂x (t, G t )dG t + 1 2 ∂ 2 f ∂x 2 (t, G t )dG t dG t . ( 19
)
In addition, from equation ( 16), we have

dG t = k(t)dW t and dG t dG t = k(t) 2 dW t dW t = k(t) 2 dt
Replacing those expressions into [START_REF] Rahman | Extended Polynomial Dimensional Decomposition for Arbitrary Probability Distributions[END_REF], we get

dX t = ∂f ∂t (t, G t ) + k(t) 2 2 ∂ 2 f ∂x 2 (t, G t ) dt + k(t) ∂f ∂x (t, G t )dW t .
Identifying this equation with ( 14) after substituting X t by f (t, G t ), the following functional relationship should hold for all (t, x) ∈ [0, T ] × R :

         a(t, f (t, x)) = ∂f ∂t (t, x) + k(t) 2 2 ∂ 2 f ∂x 2 (t, x) b(t, f (t, x)) = k(t) ∂f ∂x (t, x) (20)
Differentiating the second equation in [START_REF] Trefethen | Approximation theory and approximation practice[END_REF] with respect to x and replacing k(t) ∂f ∂x (t, x) by b(t, f (t, x)), we have

∂b(t, f (t, x)) ∂x × b(t, f (t, x)) = k(t) 2 ∂ 2 f ∂x 2 (t, x).
And replacing anew the above second order derivative into the first equation of the system (20), we obtain the equivalent system [START_REF] Öttinger | Stochastic processes in polymeric fluids[END_REF].

We now seek an approximate solution to the system (18) together with given initial conditions that will be specified a bit further in the sequel. The approximate solution is assumed to take the form

f (t, x) = n i=0 f i (t)L i (x) (21)
and the expression of f i (t) = f (t, x i ) is obtained as the solution to an ordinary differential equation, as stated in the following theorem.

Theorem 4.3. Assume that an approximation of f writes under the form [START_REF] Wiener | The Homogeneous Chaos[END_REF].

Then the coefficients f i (t) with i = 0, .., n satisfy the following system of decoupled ordinary differential equations

     f i (t) = - 1 2 ∂b ∂x (t, f i (t)) × b(t, f i (t)) + a(t, f i (t)) f i (0) = f (0, x i ) (22)
We obtain the initial conditions f (0, x i ) by solving the following differential equation :

     df dx (0, x) = b(0, f (0, x)) f (0, 0) = X 0 (23)
Furthermore, assuming that for all (t, x) in [0, T ] × R : ∂f ∂x (t, x) = 0, then the function k(t) of equation ( 16) is given by

k(t) = b (t, f (t, x)) ∂f ∂x (t, x) b t, n i=0 f i (t)L i (x) n i=0 f i (t)L i (x) (24)
Proof 4.4. Following the standard collocation method, we replace the approximation of f (21) into equation ( 18)

           n i=0 f i (t)L i (x) = a t, n i=0 f i (t)L i (x) - 1 2 ∂b ∂x t, n i=0 f i (t)L i (x) × b t, n i=0 f i (t)L i (x) k(t) × n i=0 f i (t)L i (x) = b t, n i=0 f i (t)L i (x)
We then impose the equations to be satisfied at collocations points x j and use the property of Lagrange polynomial L i (x j ) = δ ij to end up with :

                 f i (t) = a(t, f i (t)) - 1 2 ∂b ∂x (t, f i (t)) × b(t, f i (t)) k(t) = b t, n i=0 f i (t)L i (x) n i=0 f i (t)L i (x) (25)
The differential equation for f i (t) should be provided with initial conditions f i (0) = f (0, x i ). The initial conditions to the SDE ( 14) cannot be directly used since it only provides X 0 = f (0, 0). However, the value of f (0, x i ) can be obtained by considering the second equation of the system (20) and setting t = 0. Since k(0) = 1, this leads to the differential equation ( 23) that can be solved to get the value of f (0, x i ) for i = 0, .., n.

Remark 4.5. The second equation of ( 18) gives an a-posteriori criterion to check whether the hypothesis ( 16) is valid. Indeed, if

k(t) = b (t, f (t, x)) ∂f ∂x (t, x)
is a function of x, then the considered SDE has a solution that cannot be expressed as on explicit function of G t defined by ( 16).

Worked out examples

Before moving on to more complex numerical examples, we apply the proposed method to some classical SDE for which explicit calculations can be performed by hand.

The geometric Brownian motion

As a first example, we consider the geometric Brownian motion process which satisfies the following SDE :

dX t = µX t dt + ρX t dW t , t > 0 X t=0 = X 0 , ( 26
)
where µ and ρ are some real parameters. An analytic solutions can be found for this classical SDE and it writes, for t ≥ 0 :

X t = X 0 exp µ - ρ 2 2 t + ρW t . ( 27
)
The real functions a(t, x) and b(t, x) are define on R + × R by

a(t, x) = µx b(t, x) = ρx .
The differential equation ( 23) to be solved for the initial conditions of ( 22) is

     df dx (0, x) = ρf (0, x) f (0, 0) = X 0 .
Its solution is f (0, x) = X 0 exp (ρx) and therefore, the initial conditions of the differential equation for the coefficients f i (t) will be

f i (0) = f (0, x i ) = X 0 exp(ρx i ). Since ∂b(t, x) ∂x = ρ, the differential equation for f i (t) with i = 0, .., n simplifies to    f i (t) = µ - ρ 2 2 f i (t) f i (0) = X 0 exp (ρx i )
and its solution is

f i (t) = X 0 exp(ρx i ) exp µ - ρ 2 2 t .
We now have all the elements available to compute the function k(t) that will specify the process (G t ) t∈[0,T ] . Replacing f i (t) by its expression in the second equation of ( 25), we get after simplifications

k(t) = ρ n i=0 f i (t)L i (x) n i=0 f i (t)L i (x) = ρ n i=0 exp(ρx i )L i (x) n i=0 exp(ρx i )L i (x)
. At this point, this fraction cannot be further simplified. However, we can recognize the projection of the real function x → exp(ρx) onto a Lagrange basis at the numerator and the projection of the derivative of the same function at the denominator. Therefore, we have

k(t) = ρ n i=0 exp(ρx i )L i (x) n i=0 exp(ρx i )L i (x) ρ exp(ρx) (exp(ρx)) = ρ exp(ρx) ρ exp(ρx) = 1.
As a conclusion, the expression of the process (G t ) simplifies to

G t = t 0 k(s)dW s = t 0 dW s = W t ,
which is nothing other than a Wiener process. Therefore, an approximate solution to ( 26) can write

X t = X 0 exp µ - ρ 2 2 t n i=0 exp(ρx i )L i (W t ). (28)
Comparing the exact solution ( 27) with the approximate one (28), we notice that the term

n i=0 exp(ρx i )L i (W t ) is the projection of the real function x → exp(ρx)
onto a Lagrange basis evaluted at x = W t .

The Ornstein-Uhlenbeck process

As a second example, we consider the Ornstein-Uhlenbeck process which satisfies the following SDE:

dX t = θ(µ -X t )dt + ρdW t , t > 0 X t=0 = X 0 ,
where θ, µ and ρ are some real parameters. The analytic solutions for this SDE writes, for t ≥ 0

X t = X 0 e -θt + µ(1 -e -θt ) + ρe -θt t 0 e θs dW s (29)
The real functions a(t, x) and b(t, x) are define on R + × R by

a(t, x) = θ(µ -x) b(t, x) = ρ ,
and the differential equation ( 23) to be solved for the initial conditions of ( 22)

reduces to df dx (0, x) = ρ f (0, 0) = X 0 .
Its solution writes f (0, x) = ρx + X 0 , which gives initial conditions

f i (0) = f (0, x i ) = ρx i + X 0 .
Since for all (t, x) in R + × R we have ∂b(t, x) ∂x = 0, the differential equation for f i (t) with i = 0, .., n comes down to

f i (t) = θ (µ -f i (t)) f i (0) = ρx i + X 0 and its solution is f i (t) = µ + (X 0 + ρx i -µ) exp(-θt).
As before, we compute the function k(t) by replacing f i (t) by its expression in the second equation of ( 25)

k(t) = ρ n i=0 (µ + (X 0 + ρx i -µ)e -θt ) L i (x) = ρ (µ + (X 0 -µ)e -θt ) n i=0 L i (x) + ρe -θt n i=0 x i L i (x)
.

The denominator of the above expression can be simplified by using two properties of Lagrange's polynomials, i.e. for all x in R,

n i=0 L i (x) = 0 and n i=0 x i L i (x) = 1; leading to k(t) = ρ ρe -θt = e θt .
The same properties can also be used to simplify the approximate solution of the SDE as follows :

X t = n i=0 µ + (X 0 + ρx i -µ)e -θt L i (G t ) = X 0 e -θt + µ(1 -e -θt ) + ρe -θt n i=0 x i L i (G t ), = X 0 e -θt + µ(1 -e -θt ) + ρe -θt G t ,
where in the last equation, we have used the additional property: for all x in R, n i=0

x i L i (x) = x. Then, the proposed numerical method applied to the Ornstein-Uhlenbeck SDE gives the solution

X t = X 0 e -θt + µ(1 -e -θt ) + ρe -θt G t with G t = t 0 k(s)dW s = t 0 e θs dW s .
This is precisely the exact solution of the SDE (see equation ( 29)). This comes as no surprise, since the exact solution can be expressed as a degree one polynomial function of G t when G t takes the form ( 16) with k(s) = e θs .

Numerical examples

In this section, we tackle two SDEs and we explain how each step of the proposed method can be numerically addressed.

SDE with known analytical solution

As a first numerical example, we consider the following SDE:

   dX t = 1 2 (tanh(X t )) 3 dt + 1 cosh(X t ) dW t , t > 0 X t=0 = arcsinh(1)
of known analytical solution (see [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF][START_REF] Kloeden | Numerical solution of SDE through computer experiments[END_REF])

X t = arcsinh e t/2 1 + t 0 e -s/2 dW s . (30)
We start by computing (n+1) Chebyshev collocation points (see [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF]) {x i } 0≤i≤n in the interval [a, b] = [-6, 6] (the choice of this interval will be justified in the sequel). The next step is to solve the ordinary differential equation [START_REF] Xiu | Efficient collocational approach for parametric uncertainty analysis[END_REF] twice: once in the interval [0, 6] to get the values of f (0, x i ) when x i ∈ [0, 6] and once in the interval [-6, 0] to get the values of f (0, x i ) when x i ∈ [-6, 0]. We have to proceed that way since the initial condition to this differential equation is given at x = 0. For the equations considered in this paper, the choice of the numerical solver was found to be of little importance (4th order Kunge-Kutta for example) since high accuracies could be obtained with mild computing efforts. The (n + 1) ordinary differentials equations [START_REF] Wong | On the relation between ordinary and stochastic differential equations[END_REF] are then solved numerically to get an approximate solution of f i (t) at points t l ∈ [0, T ]. Notice that they all take the same form but have different initial conditions. Here again, a standard ODE solver may be used. The function of time k(t) that serves to define the process (G t ) t∈[0,T ] can be evaluated with equation ( 24) at all discrete points t l in time and at all collocation points x j as follows:

k(t l ) = b (t l , f i (t l )) n i=0 f i (t l )L i (x j )
If the exact solution to the SDE can be expressed as an explicit function of the Gaussian process defined by [START_REF] Kloeden | Numerical solution of SDE through computer experiments[END_REF] this expression should be independent of x j . This gives an a-posteriori criterion for establishing whether the considered SDE can be processed with our numerical method.

The computation of the mean and the variance of X t requires an evaluation of the standard deviation σ(t) of G t . According to [START_REF] Oksendal | Stochastic differential equations[END_REF] and after a simple change of variable, this is done by evaluating the function (for all points t l ∈ [0, T ])

σ(t l ) = t l 1 0 k(t l s) 2 ds.
The above integral is computed with Fejér quadrature rule [START_REF] Dahlquist | Numerical methods in scientific computing[END_REF] based on the roots of Chebychev polynomials of first kind. The probability density function of G t at instants of time t l ∈ [0, T ] is then given by

p Gt (t l , x) = 1 √ 2πσ(t l ) e - x 2 2σ(t l ) 2 ,
and it is possible to compute the mean, the variance and the autocorrelation function of the process (X t ) t∈[0,T ] with formulas ( 11), ( 12) and ( 13), respectively. The integrals on R are truncated to the interval [-6, 6] and they are also evaluated with Fejér quadrature rules. The truncation of the integrals do not affect the accuracy of the method due to the fast decay of p Gt (t l , x) with respect to |x|. In this example we can make this assertion more precise, since σ(t) can be computed explicitly:

σ(t) = t 0 (e -s/2 ) 2 ds = √ 1 -e -t . (31)
For all t ≥ 0, the maximum value of σ(t) is one, so for all t ≥ 0, if |x| ≥ 6,

0 ≤ 1 √ 2π e - x 2 2σ(t) 2 ≤ 1 √ 2π e - x 2 2 ≤ 1 √ 2π e - 6 2 2 6 × 10 -9 .
The interval of time for the simulations is chosen to be [0, T ] = [0, 4] and the time step for representing the results are ∆t = t l+1 -t l = 10 -2 . Notice that for this method, it is irrelevant to study the weak or strong convergence in time. Indeed, the differential equations ( 22) and ( 23) can be solved up to machine's precision with standard ODE solvers and the major part of the error comes from the interpolating error of Lagrange's polynomials. What we consider as an "exact" mean or an "exact" variance of X t can be numerically computed from the solution (30) using high precision quadrature rules. We start by rewriting (30) as

X t = arcsinh e t/2 (1 + G t ) ,
with (G t ) t∈[0,T ] a Gaussian process defined by

G t = t 0
e -s/2 dW s . For that, a sample path for the process (G t ) t∈[0,T ] needs to be generated from one trajectory of a Wiener process (W t ) t∈[0,T ] [START_REF] Higham | An algorithmic introduction to numerical simulation of stochastic differential equations[END_REF] as follows:

G t l+1 = G t l + e -t l /2 W t l+1 -W t l ,
for the path of the exact solution; and

G t l+1 = G t l + k(t l ) W t l+1 -W t l ,
for the path of the approximated solution. We see that the two trajectories are almost indistinguishable (here n + 1 = 80). For the same level of discretization, we also represent the mean and the variance (exact and approximate) on Figures 2 and3, respectively. Here again, the two curves match perfectly well. Finally, Figure 4 shows the autocorrelation function of (X t ) t∈[0,T ] , as computed by [START_REF] Hou | Wiener chaos expansions and numerical solutions of randomly forced equations of fluid mechanics[END_REF] for (s, t) ∈ [0, 4] × [0, 4].

SDE without known analytical solution

As a second numerical example, we consider the following SDE

     dX t = 1 4 X t + 1 + X 2 t 3 + 2X t 1 + X 2 t dt + X t + 1 + X 2 t dW t , t > 0 X t=0 = 0
for which no solution can be derived analytically. This SDE is particularly challenging for Euler-like numerical methods, due to the variance of the solution that grows extremely fast with respect to t. As a consequence, the Monte-Carlo simulations to compute the mean and the variance of X t requires a very large number of trajectories to get a converged solution. At first, we limit the interval of study to [0, T ] = [0, 1] to be able to compare our numerical method with the Milstein scheme with a time step ∆t = 10 -2 . The mean and the variance are then obtained with a 10 6 samples Monte-Carlo simulation. The Milstein scheme has a weak and a strong order of convergence in time equal to 1 and it writes :

X t l+1 = X t l + a(t l , X t l )∆t + b(t l , X t l )∆W + 1 2 b(t l , X t l ) ∂b ∂x (t l , X t l ) ∆W 2 -∆t X t0 = 0
, with ∆W = W t l+1 -W t l . Figure 5 and Figure 6 show the mean and the variance of X t for our method with n + 1 = 60 collocation points and for the Milstein/Monte-Carlo simulation. We can see a good match of the two curves for times t ≤ 0.6 but it quickly deteriorates for the Milstein scheme as we advance in time.

We now study the convergence of our scheme for long time simulations by 20 solution is a function of G t = W t . However, this function cannot be analytically computed. Finally, Figure 9 shows the autocorrelation function of X t for (s, t) ∈ [0, 1] × [0, 1]. We have limited the interval for (s, t) to [0, 1] × [0, 1] otherwise, the plot would be a flat land with an upsurge for the highest values of s and t.

Conclusion

A new class of high order numerical method for the solution to SDEs is proposed in this paper. It mainly relies on two standard mathematical and numerical ingredients: polynomial interpolation and cubature rules. However, it is based on a new paradigm: the white noise is no longer discretized under the form of a truncated Karhunen-Loeve series, as in Polynomial chaos expansions, and the stochastic process (X t ) t∈[0,T ] is no longer expressed as a finite number of independent random variables. Instead, (X t ) t∈[0,T ] (not necessarily Gaussian) is directly sought as a polynomial function of a Gaussian process (G t ) t∈[0,T ] and the coefficients of the expansion are obtained by solving ordinary differential equations (ODE). The macroscopic quantities (mean, variance,...) can be computed directly from the coefficients of the decomposition. Thus, two major drawbacks of competing methods can be avoided: the curse of dimensionality for Polynomial chaos expansions and the resort to Monte-Carlo simulations, with slow convergence, for Euler-like or exact sampling techniques. Assuming that the ODEs can be solved with high degree of accuracy, then the error of our method mainly comes from the interpolation error.

Several examples were studied to illustrate the efficiency and the accuracy of the proposed method. First considering the Ornstein-Uhlenbeck and the geometric Brownian motion SDEs, hands calculations have shown that our method leads to an approximate solution that is the projection of the exact solution onto a Lagrange polynomial basis. Then, more complex SDEs were numerically solved, showing the high rate of convergence and the efficiency of our approach.

Extension of the method to stochastic processes (G t ) t∈[0,T ] that are not necessarily Gaussian would allow to tackle a wider class of SDEs. This is the subject of our ongoing research.

Figure 1 :

 1 Figure 1: Sample trajectory for Xt computed from equation (35) (exact) and from its approximation equation (36)

Figure 2 :

 2 Figure 2: "Exact" and approximate mean of Xt computed with equation (32) and equation (11) respectively

Figure 3 :

 3 Figure 3: "Exact" and approximate variance of Xt computed with equation (33) and equation (12) respectively

Figure 4 :

 4 Figure 4: Autocorrelation function of Xt computed with equation (13)

Figure 5 :

 5 Figure 5: Comparison of the mean of Xt computed with our stochastic collocation method and with a Milstein scheme

Figure 6 :

 6 Figure 6: Comparison of the variance of Xt computed with our stochastic collocation method and with a Milstein scheme

Figure 7 :Figure 8 :

 78 Figure 7: Mean of Xt for different degrees of Lagrange's polynomials

Figure 9 :

 9 Figure 9: Autocorrelation function of Xt computed with equation (13)

Then for any t ∈ [0, T ] the probability density function of G t is

with σ(t) given by (31). Denoting by µ Xt (t) and σ 2 Xt (t) the exact mean and variance of X t , we have

Such integrals are evaluated after truncating their support to [-10, 10] and using Fejér quadrature rule with 10 3 quadrature points. The values obtained will be considered as "exact" values for µ Xt (t) and σ Xt (t) in the sequel. Let mean exact and mean approx be two vectors containing the values of the exact mean and the approximated mean at times {t l } 1≤l≤400 . The relative error on the mean is defined as

where the norms above are the Euclidean norms. The relative error on the variance is defined in a similar way. Table 1 shows such errors for different levels of discretization when the number of collocations points vary from n + 1 = 40 up to n + 1 = 180.

n + 1 40 60 80 100 ε mean 6.82 × 10 -3 3.24 × 10 -3 6.96 × 10 -4 3.17 We can see the high rate of convergence of the proposed method. We do not provide CPU times since it only takes a few second to run for the coarse level of discretization and up to a few minutes for the finest level of discretisation. A sample trajectory is represented on Figure 1 for the exact process X t given by X t = arcsinh e t/2 (1 + G t ) ,

and its approximation given by