A Dynamical Sparse Grid Collocation Method for Differential Equations Driven by White Noise

TL;DR

This paper introduces a sparse grid stochastic collocation method for efficiently simulating long-time behavior of stochastic differential equations driven by white noise, maintaining accuracy with limited computational resources.

Contribution

It develops a novel evolving sparse quadrature rule approach with a restarting scheme to handle long-time SDE simulations efficiently.

Findings

Accurately captures long-time solutions with small degrees of freedom.

Demonstrates effectiveness on low-dimensional nonlinear SDEs.

Maintains fixed computational complexity over time.

Abstract

We propose a sparse grid stochastic collocation method for long-time simulations of stochastic differential equations (SDEs) driven by white noise. The method uses pre-determined sparse quadrature rules for the forcing term and constructs evolving set of sparse quadrature rules for the solution variables in time. We carry out a restarting scheme to keep the dimension of random variables for the forcing term, therefore also the number of quadrature points, independent of time. At each restart, a sparse quadrature rule for the current solution variables is constructed based on the knowledge of moments and the previous quadrature rules via a minimization procedure. In this way, the method allows us to capture the long-time solutions accurately using small degrees of freedom. We apply the algorithm to low-dimensional nonlinear SDEs and demonstrate its capability in long-time simulations numerically.

Tables

Table 1. Table 1: The accuracy of quadrature rules for two independent Gaussian variables using different transformations.

	Without scaling	With scaling On	Hermite poly.
cond(A)	7.01E+9	1.91E+3	2.98E+2
${\Vert b-\tilde{b}\Vert }_{\infty }/{\Vert b\Vert }_{\infty }$	3.25E-8	3.95E-13	$3.19$E-$14$

Table 2. Table 2: Cumulants obtained by Algorithm 1 and Fokker-Planck equation at T = 8 𝑇 8 T=8 for the Ornstein-Uhlenbeck process with random damping.

	${\kappa }_1$	${\kappa }_2$	${\kappa }_3$	${\kappa }_4$	${\kappa }_5$	${\kappa }_6$
DSGC: $N=4$	2.09E-5	4.39	1.75E-4	6.06	-7.38	76.00
DSGC: $N=6$	2.09E-5	4.39	1.75E-4	6.06	1.96E-3	33.85
Fokker-Planck	0	4.39	0	6.06	0	33.85

Table 3. Table 3: Cumulants obtained by Algorithm 1 and Fokker-Planck equation at T = 4 𝑇 4 T=4 for the cubic nonlinearity.

	${\kappa }_1$	${\kappa }_2$	${\kappa }_3$	${\kappa }_4$	${\kappa }_5$	${\kappa }_6$
DSGC: $N=4$	1.26E-2	7.35E-1	4.12E-2	-3.22E-1	2.08E-1	4.75E-1
DSGC: $N=6$	-2.01E-3	7.34E-1	8.20E-3	-3.39E-1	-7.27E-2	9.35E-1
DSGC: $N=8$	3.48E-4	7.33E-1	-2.91E-3	-3.39E-1	2.40E-3	9.52E-1
Fokker-Planck	0	7.33E-1	0	-3.39E-1	0	9.64E-1

Table 4. Table 4: Errors of the variance at T = 1 𝑇 1 T=1 , and relative timings of DSGC and MC methods using different degrees of freedom.

Algorithm	${\sigma }_u=0.5$	${\sigma }_u=1$	${\sigma }_u=2$	Time ratio
DSGC $N=3$	1.2E-4	2.2E-3	5.7E-3	0.5
DSGC $N=4$	2.9E-5	1.5E-3	1.9E-3	0.8
DSGC $N=5$	1.5E-5	1.2E-3	1.3E-3	1.3
DSGC $N=6$	1.1E-5	6.3E-4	6.4E-4	2.0
MC $M_{samp}=$1E+4	1.21E-2	2.01E-2	2.46E-2	0.125
MC $M_{samp}=$2E+4	7.7E-3	1.30E-2	1.74E-2	0.25
MC $M_{samp}=$4E+4	5.7E-3	1.00E-2	1.20E-2	0.5
MC $M_{samp}=$8E+4	3.5E-3	6.03E-3	7.8E-3	1.0

Table 5. Table 5: Errors of the mean and the variance at T = 1 𝑇 1 T=1 , and relative timings of DSGC, DgPC, and MC methods using different degrees of freedom.

Algorithm	${ϵ}_{\text{mean}}$	${ϵ}_{\text{var}}$	Time ratio
DSGC $N=3$	4.6E-4	3.01E-2	0.095
DSGC $N=5$	9.2E-6	4.5E-3	0.2
DgPC $N=2$	3.8E-4	3.7E-3	2.28
DgPC $N=3$	3.2E-4	1.8E-3	6.71
MC $M_{samp}=$4E+4	4.3E-3	6.2E-3	0.5
MC $M_{samp}=$8E+4	2.6E-3	4.8E-3	1.0

Equations

$$I_Q(u)(\xi ):=\sum _{q=1}^Q{w}^{q}u({\xi }^q),$$

$$\mathbb{E}[u(\mathbf{\xi })]\approx I_Q^{\otimes K}:=\sum _{{\alpha }_1=1}^Q\dots \sum _{{\alpha }_K=1}^{Q}u({\xi }_1^{{\alpha }_1},\dots ,{\xi }_K^{{\alpha }_K})w_1^{{\alpha }_1}\dots w_K^{{\alpha }_K}.$$

$$\mathbb{E}[u]\approx \sum _{\lambda \le |\mathbf{\alpha }|\le \lambda +K-1}(-1{)}^{\lambda +K-|\mathbf{\alpha }|-1}(\genfrac{}{}{0ex}{}{K-1}{|\mathbf{\alpha }|-\lambda })(I_{{\alpha }_1}\otimes \dots \otimes I_{{\alpha }_K})(u),$$

$$du(t)=\mathcal{L}(u(t))dt+\sigma (u(t))dW(t),u(0)=u_0,t\in [0,T],$$

$${\tau }_i=i\Delta \tau ,i=0,\dots ,M_T,$$

$$u({\tau }_{i+1})$$

$$\mathbb{E}{\left[W(t)-\sum _{k=1}^K{\xi }_k{\int }_0^t{m}_k(s)ds\right]}^2\to 0,K\to \infty ,$$

$$u({\tau }_{i+1};u_0,\mathbf{\xi })$$

$$u({\tau }_{i+1};u_0^p,{\mathbf{\xi }}^q)$$

$$\mathbb{E}[u(t;u_0,\mathbf{\xi }{)}^{\mathbf{\alpha }}]\approx \sum _{p=1}^{Q_{u_0}}\sum _{q=1}^{Q_{\mathbf{\xi }}}{w}_0^p{\mathbf{w}}^{q}u(t;u_0^p,{\mathbf{\xi }}^q{)}^{\mathbf{\alpha }},$$

$$u_{j+1}({\tau }_{j,i+1};u_j^p,{\mathbf{\xi }}_j^q)$$

$$+\sigma (u_{j+1}({\tau }_{j,i};u_j^p,{\mathbf{\xi }}_j^q))\sum _{k=1}^K{\xi }_{j,k}^q{\int }_{{\tau }_{j,i}}^{{\tau }_{j,i+1}}{m}_{j,k}(s)ds,$$

$$J_{d,N}:=\{\mathbf{\alpha }=({\alpha }_1,\dots ,{\alpha }_d)\in {\mathbb{N}}_0^d,|\mathbf{\alpha }|\le N\}.$$

$$\underset{w_{j+1}\in {\mathbb{R}}^{{\tilde{Q}}_{u_{j+1}}}}{\min }||w|{|}_1,\text{mbox}subjecttoAw=b,$$

$$A:=\left[\begin{array}{ccc}{T}_0(u_{j+1}^1)& \dots & T_0(u_{j+1}^{{\tilde{Q}}_{u_{j+1}}})\\ ⋮& & ⋮\\ T_M(u_{j+1}^1)& \dots & T_M(u_{j+1}^{{\tilde{Q}}_{u_{j+1}}})\end{array}\right],$$

$$b:={\left[\begin{array}{ccc}\mathbb{E}[T_0(u_{j+1})]& \dots & \mathbb{E}[T_M(u_{j+1})]\end{array}\right]}^T.$$

$$\mathbb{E}[u_{j+1}^{\mathbf{\alpha }}]={\mathbb{E}}_{(u_j,{\mathbf{\xi }}_j)}[u_{j+1}^{\mathbf{\alpha }}(u_j,{\mathbf{\xi }}_j)]\approx \sum _{p=1}^{Q_{u_j}}\sum _{q=1}^{Q_{{\mathbf{\xi }}_j}}{w}_j^p{\mathbf{w}}_j^q\prod _{i=1}^d(u_{j+1,i}(u_j^p,{\mathbf{\xi }}_j^q){)}^{{\alpha }_i}.$$

$$m_{j,1}(t)=\frac{1}{\sqrt{t_{j+1}-t_j}},m_{j,k}(t)=\sqrt{\frac{2}{t_{j+1}-t_j}}\cos \left(\frac{(k-1)\pi (t-t_j)}{t_{j+1}-t_j}\right),k\ge 2,$$

$$\displaystyle\epsilon_{\mbox{mean}}(t):=\bigg{|}\frac{\mu_{\mbox{pce}}(t)-\mu_{\mbox{exact}}(t)}{\mu_{\mbox{exact}}(t)}\bigg{|},\quad\epsilon_{\mbox{var}}(t):=\bigg{|}\frac{\sigma^{2}_{\mbox{pce}}(t)-\sigma^{2}_{\mbox{exact}}(t)}{\sigma^{2}_{\mbox{exact}}(t)}\bigg{|},$$

$$du(t)=b_u({\mu }_u-u(t))dt+{\sigma }_{u}dW(t),t\in [0,T],u(0)=u_0,$$

$$du=-(u^2+1)udt+{\sigma }_{u}dW,u(0)=1.$$

$$du(t)=b_u({\mu }_u-u(t))dt+{\sigma }_u\sqrt{u(t)}dW(t),$$

$$u({\tau }_{i+1})$$

$$+{\sigma }_u\sqrt{u({\tau }_i)}(W({\tau }_{i+1})-W({\tau }_i))+\frac{{\sigma }_u^2}{4}(W({\tau }_{i+1})-W({\tau }_i){)}^2,$$

$$du(t)$$

$$dv(t)$$

Full text

Keywords: Stochastic collocation, stochastic differential equations, long-time integration, white noise

1 Introduction

Stochastic collocation methods have been popular tools for computing expectations of functionals of solutions to stochastic differential equations (SDEs) driven by white noise [38, 23, 24, 20]. Undoubtedly, the most commonly used approach to compute expectations has been Monte Carlo (MC) method, which is a sampling technique in the random parameter space. However, random sampling-based methods are notorious for their slow convergence, which hinders their computational efficiency in simulations; especially in long-time simulations. Although several extensions have been proposed to speed up the convergence, their speed does not compete with fast convergence of spectral methods; see [16, 21, 24, 33, 25] and references therein.

In this work, we focus on stochastic collocation methods which use deterministic quadrature nodes in the random space to approximate expectations of functionals of solutions. These methods aim to achieve the ease of implementation of MC methods and fast convergence behavior of stochastic Galerkin methods at the same time. Similar to Galerkin methods, their convergence depends on the smoothness of the solution in the random input variables. Their effectiveness relies upon the dimensionality of the random parameter space and they work reasonably well if the stochastic system involves moderately large random dimensions. It has been shown, especially in uncertainty quantification literature, that they provide a strong alternative to MC methods for differential equations with time-independent random forcing; see e.g., [50, 52, 4, 35, 34, 53, 25].

For equations driven by white noise, collocation methods with pre-determined quadrature rules have appeared in the recent literature [15, 14, 55, 56, 30, 26]. Manuscripts [55, 56] combine a sparse grid collocation (SGC) method for white noise with weak-sense time integration methods to compute second order moments of the solutions of SDEs and stochastic partial differential equations (SPDEs). It has been proved and observed numerically in [55] that straightforward application of collocation methods leads a failure in long-time integration. The underlying reason of this failure is that in the presence of random forcing in time, the number of stochastic variables needed to maintain a prescribed accuracy increases with time. This means, for collocation-based methods, that the number of collocation points should grow with time. The manuscript [55] then introduces a recursive approach for long times on top of a straightforward implementation of SGC to compute moments of linear SPDEs; see also [27, 32]. Similar long-time integration issues have been noted before in polynomial chaos expansion-based methods in [48, 12, 6, 39, 40]. In contrast to collocation methods with pre-determined quadrature rules, optimal quantization methods [41, 29, 42] aim to find optimal discrete approximations, e.g. Voronoi quantizations, to the solutions of SDEs, which are tailored for the underlying dynamics. These quantizations are usually obtained by approximating Brownian motion by a finite-dimensional random process, e.g. truncated Karhunen-Loeve expansion (KLE) of Brownian motion, and deriving systems of ordinary differential equations for the quantizers. Although resulting quadrature rules are adapted to the measures under consideration, slow convergence rates limit use of optimal quantization methods in numerical applications without further acceleration techniques [41].

A possible remedy to the long-time integration problem in the setting of SDEs is to exploit Markov property of the dynamics and allow the algorithm to forget about the past. This observation is crucially leveraged in our previous works [39, 40] to construct dynamical polynomial chaos expansions, called Dynamical generalized Polynomial Chaos (DgPC), in terms of the solution variables and the random forcing; see also [54] for a different approach. DgPC captures the solutions of SDEs and SPDEs accurately in the long-time with less computational time compared to MC methods; especially for complicated SPDE dynamics. In this work, using similar principles of [39, 40], we propose a dynamical collocation-based method in time to alleviate long time integration of SDEs. We propagate optimal quadrature rules for the solution in time and use pre-determined quadrature rules for the random forcing. In this sense, the method can be considered as an extended combination of the proposed methods in [55, 41]. The method uses a restarting scheme similar that of [39, 40] to construct sparse quadrature rules for the solution variables on-the-fly. It then estimates expectations of functionals of future solution variables by using sparse quadrature rules of the random forcing variables and the solution variables at the current time. By constructing quadrature rules with small number of nodes and employing frequent restarts, the algorithm can utilize small and time-independent degrees of freedom while maintaining accuracy in the long-time. We demonstrate the efficiency of the proposed method numerically in several cases using low-dimensional nonlinear SDEs forced by white noise.

The outline of the paper is as follows. Section 2 introduces the necessary background material on sparse quadrature methods. A simple stochastic collocation method is introduced in Section 3 and its drawbacks are discussed. Our methodology, called Dynamical Sparse Grid Collocation (DSGC), is described in detail in Section 4. Numerical experiments for nonlinear SDEs are presented in Section 5. Some conclusions are offered in the last section.

2 Sparse Quadrature Methods

Approximation of an output $u(\boldsymbol{\xi})\in\mathbb{R}^{d}$ of a model depending on a stochastic input variables $\boldsymbol{\xi}\in\mathbb{R}^{K}$ requires the solution of the corresponding forward propagation operator. This approximation procedure usually involves computation of expectations of the functional $u(\boldsymbol{\xi})$. For instance, for polynomial chaos expansions (PCEs), the coefficients of the representation can be obtained by a spectral projection onto associated orthogonal polynomial basis $T_{\boldsymbol{\alpha}}(\boldsymbol{\xi})$ and are given in the form of expectations $\mathbb{E}[u(\boldsymbol{\xi})\,T_{\boldsymbol{\alpha}}(\boldsymbol{\xi})]$ with respect to the probability measure of $\boldsymbol{\xi}$. The dimension $K$ of $\boldsymbol{\xi}$ is usually large in most applications, which in turn requires efficient estimations of high-dimensional integrals. In the following, we focus on quadrature-based collocation methods to compute expectations.

Assuming the probability distribution of the input variables $\boldsymbol{\xi}$ is known, e.g. it is given in the Askey family [51], and the components are independent, one can construct multi-dimensional quadrature rules by employing tensorization of one-dimensional quadrature rules. For the sake of brevity, we consider evaluation of $K$-dimensional integrals of the output $u(\boldsymbol{\xi})$ with respect to a Gaussian measure $p_{\boldsymbol{\xi}}(\boldsymbol{\xi})\propto\exp(-||\boldsymbol{\xi}||_{2}/2)$.

In one dimension, we define the Gauss-Hermite quadrature rules $I_{Q}$ consisting of the weights and nodes $\{w^{q},\xi^{q}\}_{q=1}^{Q}$, $Q\in\mathbb{N}$ and $\xi\in\mathbb{R}$:

[TABLE]

where $\xi^{q}$’s are the roots of the Hermite polynomials $H_{Q}$ and the weights are given by $w^{q}=1/(Q^{2}(H_{Q-1}(\xi^{q}))^{2})$. It is known that $I_{Q}$ is exact if $u$ is a polynomial of degree less than or equal to $2Q-1$; [9].

For the multi-dimensional case, the integral $\mathbb{E}[u(\boldsymbol{\xi})]$ can be approximated by the tensor product formula

[TABLE]

Here we use the multi-index notation: $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{K})\in\mathbb{N}^{K}$ with $|\boldsymbol{\alpha}|=\sum_{k=1}^{K}\alpha_{k}$. We denote by $Q_{\boldsymbol{\xi}}$ the total number of resulting quadrature nodes. Approximations based on this tensor product suffer from curse of dimensionality and computational costs scale exponentially with dimension $K$, i.e. $Q_{\boldsymbol{\xi}}=Q^{K}$.

If the dimension of the random variables is moderately high, a sparse quadrature rule first proposed by Smolyak, can be used to reduce the number of quadrature nodes while maintaining the accuracy [45]. Following [49, 15], we write the sparse grid approximation to the multi-dimensional integral with the level $\lambda$

[TABLE]

where $\boldsymbol{\alpha}\geq\mathbf{1}$. This quadrature rule is exact for multivariate polynomials of total degree up to $2\lambda-1$ and greatly reduces the number of evaluations compared to the tensor product rule above [15, 37]. In this work, we employ isotropic Smolyak sparse grid quadrature rules for Gaussian measures, meaning that the level $\lambda$ is the same for each dimension. We also note that the weights of this sparse quadrature rule can be negative.

For stochastic equations with random input data, sparse grid collocation methods have been shown to be effective, especially equations with time-independent forcing and parameters, and in moderately high dimensions; see [50, 52, 53, 4, 35, 34] and references therein. For stochastic equations with time-dependent noise, there are mainly two methods which use pre-determined set of quadrature nodes and weights: cubatures on Wiener space [30, 26] and sparse grid Smolyak quadrature rules [15, 14, 55, 56]. In [55], applications of sparse grid collocation methods to equations driven by white noise have been studied in conjunction with Monte-Carlo methods and it is noted that sparse grid collocation works only for short time integration. Motivated by this observation, in the next section, we introduce a sparse, deterministic collocation method based on quadrature rules for SDEs driven by white noise using a spectral expansion of the white noise and confirm the findings in [55].

3 A stochastic collocation method for SDEs

We consider the following $d$-dimensional SDE:

[TABLE]

where $\mathcal{L}(u)$ is a general function of $u\in\mathbb{R}^{d}$ (and possibly other deterministic or stochastic parameters), $W(t)$ is a Brownian motion, $\sigma(u)$ is the amplitude which depends on the solution, and $u_{0}$ is an initial condition with a prescribed probability density. The stochastic integral is considered in the Ito sense. We assume that the solution exists and is a second order stochastic process. We also use the notation $W(t,\omega)$ for Brownian motion $W(t)$.

We introduce a sparse quadrature-based collocation method for (3.1) as follows. First, we consider a time discretization of the interval $[0,T]$

[TABLE]

where $\Delta\tau=T/M_{T}$. Then we approximate the solution $u(t;u_{0},\{W({\tau});\tau\leq t\})$ of (3.1) by the Euler method [24]

[TABLE]

Given a complete countable orthonormal system $m_{k}(t)\in L^{2}[0,T]$, $k\in\mathbb{N}$, we project the Brownian motion $W(t,\omega)$ onto $L^{2}[0,T]$ by defining $\xi_{k}(\omega):=\int_{0}^{T}m_{k}(t)\,dW(t,\omega)$. Then, the random vector $(\xi_{1},\xi_{2},\ldots)$ consists of a countable number of independent and identically distributed standard Gaussian random variables, and the convergence

[TABLE]

holds for all $t\leq T$; see [22, 28, 31]. This convergence property allows us to approximate the increments $(W(\tau_{i+1})-W(\tau_{i}))$ of the random forcing by the finite dimensional random variable $\boldsymbol{\xi}=(\xi_{1},\xi_{2},\ldots,\xi_{K})$:

[TABLE]

where ${u}(\tau_{0};u_{0},\boldsymbol{\xi})=u_{0}$.

Now, let $\{w_{0}^{p},u_{0}^{p}\}_{p=1}^{Q_{u_{0}}}$ and $\{\mathbf{w}^{q},\boldsymbol{\xi}^{q}\}_{q=1}^{Q_{\boldsymbol{\xi}}}$ be pre-determined quadrature rules for the random variables $u_{0}$ and $\boldsymbol{\xi}$. For instance, we can use Smolyak sparse grid quadrature (2.1) for Gaussian $\boldsymbol{\xi}$ and for $u_{0}$, any accurate quadrature rule can be considered. Here these quadrature rules denote any enumerations of their multi-dimensional versions. Then equation (3.3) naturally gives rise to a non-intrusive approximation to (3.1) by the following equation:

[TABLE]

for $p=1,\ldots,Q_{u_{0}}$ and $q=1,\ldots,Q_{\boldsymbol{\xi}}$. Equation (3.4) dictates the evolution of the initial particles $u_{0}^{p}$ under the trajectories of the forcing particles $\boldsymbol{\xi}^{q}$. In contrast to Monte Carlo methods, the random forcing is sampled deterministically and approximated by its finite-dimensional approximation via the spectral projection (3.2), and the samples of $u_{0}$ are taken as quadrature points. Thus, the method is sample-error free.

Remark 3.1**.**

The method introduces three level of approximations. First, SDE is discretized in time. Then, Brownian motion increments are approximated by their finite-dimensional approximations. Finally, the semi-discrete equation (3.3) is approximated by its fully discrete version (3.4) using quadrature rules. Hence, there are three degrees of freedom that are of interest: $\Delta\tau$, $K$, and $\lambda$. Note also that although we used the Euler method in the formulation, this is not required. Any higher order method can be used to discretize the SDE in time.

Remark 3.2**.**

For any fixed $K$, the finite dimensional approximation (3.2) of Brownian motion entails a continuous finite-variation process on $[0,T]$. Thus, the integral with respect to this finite-variation process can be understood in the Stieltjes sense. Then the main questions are when and in what sense the approximate solution (3.4) converges to the true solution. Here we are only intereseted in the numerical efficiency of the method and we refer to [41, 29, 42] for theoretical discussions on the convergence of approximations for SDEs with smooth coefficients in a similar setting. Nevertheless, our numerical experiments show clear convergence in moments; see Section 5.

The approximate solution $u(t;u_{0}^{p},\boldsymbol{\xi}^{q})$ of (3.4) readily entails approximations to statistical moments by computing

[TABLE]

where $\boldsymbol{\alpha}$ is a multi-index, see (4.2).

A similar collocation strategy in a Monte-Carlo setting using weak-integration is employed in [55, 56] to compute first two moments of the solution. It is noted that the efficiency of the collocation strategy depends on the strength of the noise and the length of the time interval. Indeed, in order to maintain a prescribed accuracy the expansion (3.2) requires the number of stochastic variables to be increasing with time. Thus, the number of quadrature points needed to maintain an accuracy becomes quickly overwhelming for long times.

Here is a simple motivating demonstration of the long-time integration problem in case of the Ornstein-Uhlenbeck process. We set $\mathcal{L}(u)=10(0.1-u)$ and $\sigma=4$, and take a deterministic initial condition $u_{0}=1$. We also take $K=8,16,32,64$, and consider different final times $T=1,2,4,8,16$. Since the solution stays Gaussian, we use a sparse Gauss-Hermite rule for $\boldsymbol{\xi}$ with level $\lambda_{\boldsymbol{\xi}}=1$. To make a fair comparison, we use the same time discretization method (second order Runge-Kutta method) with the time step $\Delta\tau=1$E-$3$ in each scenario. Figure 1 shows that as the time increases from $T=1$ to $T=16$, the convergence behavior in the variance of this simple method decreases from $O(K^{-3})$ to $O(K^{-1})$. This clearly indicates that the degrees of freedom required for this simple collocation method to maintain a desired accuracy should increase with time.

4 Proposed methodology

4.1 Motivation

As our previous works [39, 40] and other manuscripts [32, 22, 54, 55, 56] noted before, the main difficulty for differential equations driven by white-noise forcing that should be addressed is that the number of variables need to be accounted for grows in time. In the settings that dynamics satisfy a Markov property, any adapted algorithm can and should utilize this fact to “forget” the past variables and as a consequence keep the dimension of the random variables independent of time. Basically, this formulation leverages the intrinsic sparsity of the dynamics in the sense that a sparse representation of future times can be obtained using the solution at the current time and the corresponding random forcing variables. In [39, 40], utilizing these observations, we introduced a PCE-based method which utilizes a restarting procedure to compute orthogonal polynomials of the solution dynamically on-the-fly and constructs evolving polynomial expansions to represent future solution in a reasonably sparse representation. The method performed on par with standard MC methods in most cases and, in case of complex dynamics, it can be executed faster than standard MC [40]. Although, PCE-based restart methods provide a viable means to keep the number of random variables in reasonable levels and to compute statistical properties of long-time evolutions for fairly complicated equations, the methods in [39, 40] are intrusive. In this connection, we present a non-intrusive method which exploits the same forgetting strategies to construct dynamical quadrature-based collocation rules for the solution in time. For other forgetting strategies, we refer the reader to [32, 26, 54, 55, 56].

4.2 Formulation

The solution $u(t,\omega)$ of the $d$-dimensional SDE (3.1) is a nonlinear functional of the initial condition and countable number of variables $\boldsymbol{\xi}$. The restart methods in [39, 40] crucially depend on the following two observations: (i) for sufficiently small time lags $\varepsilon>0$, the solution $u(t+\varepsilon)$, can be efficiently captured by low-order polynomials in $u(t)$; and (ii) the solution $u(t+\varepsilon)$ depends on Brownian motion only on the interval $[t,t+\varepsilon]$ by Markov property. Based on these observations, the algorithms in[39, 40] introduced an increasing sequence of restart times $0<t_{j}<t_{j+1}<T$, with $\Delta t=t_{j+1}-t_{j}$, and constructed a new PCE basis at each $t_{j}$ based on the solution $u(t_{j})$ and all additional random variables $\boldsymbol{\xi}$.

In the following, we propose to utilize quadrature rules of the solution $u(t)$ and the variables $\boldsymbol{\xi}$ to represent the solution $u(t+\epsilon)$ at future times with a sufficiently small number of quadrature points. By the fact (i), $u(t+\epsilon)$ can be captured by low order polynomials in $u(t)$, therefore, the quadrature level required to integrate polynomials in $u(t)$ can be selected small. Moreover, the number of quadrature points for $\boldsymbol{\xi}$ on each time interval $[t_{j},t_{j+1}]$ can also be made small by selecting a short time horizon. Notice also that by the properties of Brownian motion, the quadrature rule for $\boldsymbol{\xi}$ can be read from tables or computed only once in the offline stage. The challenge is then to compute efficient, sparse quadrature rules for the solution $u(t)$ in time. These rules are not straightforward to compute and have to be computed online since the probability distribution of $u(t_{j})$ is arbitrary and evolving.

Let a sequence of increasing restart times $t_{j}$ be given. We denote by $u_{j}$ the approximation of $u(t_{j})$ given by the algorithm at $t_{j}$ and by $\boldsymbol{\xi}_{j}=(\xi_{j,1},\ldots,\xi_{j,K})$ the variables for the random forcing on the interval $[t_{j},t_{j+1}]$. Suppose for now that the sparse quadrature rules $\{w^{p}_{j},u_{j}^{p}\}_{p=1}^{Q_{u_{j}}}$ and $\{\mathbf{w}^{q}_{j},\boldsymbol{\xi}_{j}^{q}\}_{q=1}^{Q_{\boldsymbol{\xi}_{j}}}$ have already been constructed for $u_{j}$ and $\boldsymbol{\xi}_{j}$ at time $t_{j}$, respectively. An analog of the equation (3.4) can be written for the approximate solution $u_{j+1}(t;u_{j},\boldsymbol{\xi}_{j})$ as

[TABLE]

where $m_{j,\cdot}(t)$ is a complete orthonormal system for $L^{2}[t_{j},t_{j+1}]$ and ${\tau_{j,\cdot}}$’s denote a time discretization for the interval $[t_{j},t_{j+1}]$.

The question here is how to construct an efficient quadrature rule for the next solution variable $u_{j+1}$ using (4.2). The evolution of the particles $\{u_{j}^{p},\boldsymbol{\xi}_{j}^{q}\}$ via the equation (4.2) entails a set $\{u^{p\times q}_{j+1}\}$ of particles of the approximate solution $u_{j+1}$, i.e quadrature nodes at $t_{j}$ follow the trajectories of the dynamics and give rise to an initial set of nodes at $t_{j+1}$. The challenge is to find a small subset of these nodes and corresponding weights which accurately integrates polynomials in $u_{j+1}$. Following [2, 3] and [44], we construct such a sparse quadrature for $u_{j+1}$ by using the following $L^{1}$-minimization procedure.

We define multi-indices $\boldsymbol{\alpha}$ up to degree $N$

[TABLE]

and let $M:=|J_{d,N}|={d+N\choose{d}}$ be the total number of elements in the set. Let the particles $u^{p\times q}_{j+1}:=u_{j+1}(t_{j+1};u_{j}^{p},\boldsymbol{\xi}_{j}^{q})$, where $p=1,\ldots,Q_{u_{j}}$ and $q=1,\ldots,Q_{\boldsymbol{\xi}_{j}}$, serve as an initial set of quadrature nodes for $u_{j+1}$. Let also the set $\{T_{\boldsymbol{\alpha}}(u)\,:\boldsymbol{\alpha}\in\mathcal{J}_{d,N}\}$ be any orthonormal basis of polynomials up to degree $N$ in dimension $d$. Then to extract a sparse quadrature rule, following [2, 3], we solve the convex optimization problem:

[TABLE]

where $w\in\mathbb{R}^{\tilde{Q}_{u_{j+1}}}$ with $\tilde{Q}_{u_{j+1}}=Q_{u_{j}}\times Q_{\boldsymbol{\xi}_{j}}$, and the constraints $Aw=b$ necessitate the exactness of the quadrature rule up to degree $N$. We enumerate the basis $T_{\boldsymbol{\alpha}}(u)$ and denote it by $\{T_{k}(u)\,:\,k=0,\ldots,M\}$, then define

[TABLE]

and the right-hand side vector consisting of the exact moments

[TABLE]

Here we assume that the moments $\mathbb{E}[u_{j+1}^{\boldsymbol{\alpha}}]$, for each $t_{j+1}$ and $|\boldsymbol{\alpha}|\leq N$, are finite. Also, we typically have that $M$ is much smaller than $\tilde{Q}_{u_{j+1}}$. Then a sparse subset denoted by $\{w^{p}_{j+1},u^{p}_{j+1}\}_{p=1}^{Q_{u_{j+1}}}$ with $Q_{u_{j+1}}\leq M\ll\tilde{Q}_{j+1}$, can be extracted having at most $M$ nodes from the solution of the optimization procedure (4.3); see section 4.3. Once a quadrature rule for $u_{j+1}$ is constructed, the algorithm restarts and evolves the nodes on the next time interval according to (4.2).

Remark 4.1**.**

For probability measures on $\mathbb{R}^{d}$, the existence of an exact quadrature rule with positive weights is guaranteed; [46, 43, 8, 11]. In general, we do not enforce positivity condition for the weights $w$ since a sparse optimal solution with positive weights may not exist; see Section 4.3 for further details. We note that (4.3) is not the only construction to find optimal quadrature rules; see also [17, 8, 2, 3] and references therein. Furthermore, the convergence of exact quadrature rules for compactly supported probability measures has been studied extensively and the results can be found in classical literature [11, 9].

Remark 4.2**.**

We do not tensorize quadrature rules for each component of $d$-dimensional random vector $u_{j}$. It is quite possible that components of $u_{j}$ exhibit correlation; therefore the tensorization is not optimal. However, since we impose constraints on multivariate moments $\mathbb{E}[u_{j}^{\boldsymbol{\alpha}}]$, the algorithm automatically establishes a quadrature rule for the full vector $u_{j}$. Moreover, if the dimension $d$ is high, the number of constraints $M$ in (4.3) can be reduced by considering a sparse version of the multi-index set $\mathcal{J}_{d,N}$; see [39, 40, 22, 5].

Remark 4.3**.**

This remark concerns the differences and similarities of our approach to the methods in [55, 56] and optimal quantization methods in [41, 29].

Our method uses pre-determined quadrature rules for the random forcing variables and does not a have sampling error similar to the method in [55]. The main difference is that the paper first discretizes the stochastic equations in time and then considers quadrature rules for the random forcing variables in each time-step, i.e. the dimension of randomness depends on the resolution of the time-integration and might grow rapidly with fine discretization. In contrast, our method discretizes in the random space by considering the $L^{2}$-approximation (3.2) of Brownian motion with a fixed degree of freedom $K$. Although a different restarting mechanism is used in [55], their method can only compute moments up to second order, whereas approximations to higher order moments are available in our method by quadrature rules provided approximations to higher moments converge.

Our method finds optimal quadrature rules adapted to the evolving solution in a similar sense to optimal quantization methods. The quantization methods in [41, 29] aim to discretize the paths of an infinite dimensional randomness by a random vector in finite dimension. Finite dimensional approximations are obtained by the solution of a minimization procedure and are typically given by Voronoi cell collocation. For Brownian motion, quantizers based on its KLE are considered. Then the evolution of these quantizers for SDEs is obtained by solving a corresponding integral equation, which is similar to (4.2). The convergence order is logarithmic in the number of quantizers, which is a poor rate of convergence for practical applications. In contrast, our method utilizes Gaussian-type quadratures tailored for the solution and the random forcing, which leverage the smoothness of the response to provide fast convergence. Moreover, frequent restarts allow us to mitigate dimensionality and use low-dimensional approximations to Brownian motion forcing.

Our methodology is given in Algorithm 1 below.

4.3 Implementation

4.3.1 Offline stage

In the offline stage, certain quadratures need to be computed. First, we compute sparse quadrature rule for Gaussian $\boldsymbol{\xi}_{0}$ by using the Smolyak sparse grid with the level $\lambda_{\boldsymbol{\xi}_{0}}$, which builds upon the standard 1D Gauss-Hermite rule; see Section 2. By the independent increment property of Brownian motion, all $\boldsymbol{\xi}_{j}$, $j\geq 0$, may have the same quadrature rule assuming $K$ is fixed. Note that although it is not necessary, we keep the number of variables in $\boldsymbol{\xi}$ the same throughout the evolution. If the distribution of the initial condition $u_{0}$ is known, a sparse Gauss quadrature is constructed with the level $\lambda_{u_{0}}$. If its distribution is arbitrary, then the optimization procedure (4.3) can be used with Monte Carlo initialization. We make use of the C++ library “UQ Toolkit” to compute Gauss rules [10].

4.3.2 Moments and orthogonal polynomials

At the restart time $t_{j+1}$, estimation of the right-hand side vector in the constraints in (4.3) requires the calculation of multivariate moments $\mathbb{E}[u_{j+1}^{\boldsymbol{\alpha}}]$, where $|\boldsymbol{\alpha}|\leq N$ and $u_{j+1}=(u_{j+1,1},\ldots,u_{j+1,d})$. As noted before in Section 3, these moments are computed using the available quadrature rules $\{w_{j}^{p},u_{j}^{p}\}_{p=1}^{Q_{u_{j}}}$ and $\{\mathbf{w}^{q}_{j},\boldsymbol{\xi}_{j}^{q}\}_{q=1}^{Q_{\boldsymbol{\xi}_{j}}}$ from time $t_{j}$:

[TABLE]

The optimization procedure does not depend on the particular choice of the polynomials $\{T_{\boldsymbol{\alpha}}(u)\,:\,\boldsymbol{\alpha}\in\mathcal{J}_{d,N}\}$, e.g. even monomials can be used. However, the choice of $T_{\boldsymbol{\alpha}}$ certainly affects the condition number of the constraint matrix, which in turn affects the stability of the numerical minimization algorithm. To better condition the constraints and improve convergence of the optimization algorithms, we make use of couple of linear transformations as preconditioning steps. Similar transformation techniques are applied before in moment-constrained maximum entropy problem [1].

A linear transformation is applied to the solution $u_{j}$ so that it becomes mean zero. Then a further transformation makes its components uncorrelated, i.e. its covariance matrix becomes identity. Even with these transformations, a scaling issue related to the moments arises in the constraint equations. For instance, for a standard Gaussian random variable $\xi$, we have $\mathbb{E}[\xi^{12}]/\mathbb{E}[\xi^{2}]=10395$, i.e. the twelfth moment is larger than the second moment by $5$ orders magnitude. To alleviate this scaling issue, we further scale $u$ by its largest moment so that the maximum moment becomes $1$.

A more direct preconditioning can be applied by a judicious selection of the orthonormal basis $T_{\boldsymbol{\alpha}}$. A basis can be selected using an educated guess in the offline stage, which does not require any online computation. However, since the measure of $u_{j}$ is evolving in time, this may not be optimal in the long-time in Algorithm 1. An optimal choice for $T_{\boldsymbol{\alpha}}$ is the set of polynomials which are orthogonal with respect to the measure of $u_{j}$. Unfortunately, a corresponding orthogonal system for $u_{j}$ is not available a priori in our algorithm, but it can be computed online if a further preconditioning is required. Although the computation of orthogonal polynomials is an ill-posed problem, a Gram-Schmidt procedure based on the knowledge of multivariate moments can be used in the computation. In numerical simulations, the choice of the orthonormal system will be explicitly stated. We refer to [11, 39, 40] and references therein for detailed discussions on how to compute orthogonal polynomials of $u_{j}$ dynamically.

Here is a simple demonstration of the effects of these transformations. Let $\xi_{1}$ and $\xi_{2}$ be two independent $N(3,1)$ variables and the maximum degree be $N=8$. Then, the number of constraints becomes $M=45$. We use $500$ samples from normal distribution to initialize the optimization procedure and keep the samples same for each scenario to make a fair comparison. A sparse quadrature rule with $M$ nodes is extracted according to the algorithm discussed in the next section, and, afterwards, the right-hand side vector $b$ is computed numerically using this quadrature rule to check the accuracy. The numerical approximation of $b$ is denoted by $\tilde{b}$ in the following.

Table 1 shows condition numbers of the linear system in (4.3) and the accuracy $||b-\tilde{b}||_{\infty}/||b||_{\infty}$ of the associated quadrature rule. First two scenarios use monomials as the polynomial basis $T_{\boldsymbol{\alpha}}$ and the last one uses Hermite polynomials, which are the associated orthogonal system in this example. Note finally that condition numbers are independent of the sparse extraction procedure. Clearly, scaling transformations or a careful selection of the polynomials basis leads to at least $5$-digit gain of accuracy in this example.

4.3.3 Sparse quadrature rule

Algorithm 1 constructs dynamical quadrature rules in time for the solution $u_{j}$. However, the implementation of the optimization algorithm (4.3) to construct an efficient quadrature rule is not straightforward. From (4.3), we observe that at each restart $t_{j}$, the optimization procedure is initialized with $\tilde{Q}_{u_{j}}=Q_{u_{j-1}}\times Q_{\boldsymbol{\xi}_{j-1}}$ number of nodes for $u_{j}$. Therefore, the number of quadrature nodes may grow with the number of restarts provided the solution of the optimization procedure is not already sparse. The challenge is then to compute a sparse quadrature rule containing a smaller set of nodes and weights while keeping the exactness of the original quadrature rule. To this end, after finding the optimal solution of (4.3), we further employ an extraction routine.

A straightforward sparsification of the optimal solution of (4.3) would be cutting-off the weights that are greater than a certain threshold. Depending on the numerical minimization algorithm, this may not be possible. As discussed in [2, 3], an application of simplex algorithm yields a sparse optimal solution whereas interior-point methods give fully populated solution [36]. This behavior is due to the fact that the simplex algorithm acts on the vertices of a polytope whereas interior point methods result in sparse vectors only in the limit, which is not achieved in numerics. We choose to use CVX, a package for specifying and solving convex programs [19, 18]. Under the hood, CVX uses SDPT3 which employs interior-point methods to compute the optimal solution [47]. The following procedure is used to extract a sparse quadrature rule; see [2, 3, 36].

At time $t_{j}$, the constraints matrix $A$ is of dimension $M\times\tilde{Q}_{u_{j}}$, where $\tilde{Q}_{u_{j}}$ is much bigger than $M$. Then, the dimension of the nullspace of $A$ is at least $\tilde{Q}_{u_{j}}-M$. The key observation is that any vector $z\in\mbox{null}(A)$ can be added to the weights vector without changing the equality constraints, i.e. $A(w+z)=b$. Thus, by selecting vectors carefully, we can construct an iterative routine to make most of the weights zero. Based on these observations, we follow the approach given in [2, 3, 36] and employ the following routine, Algorithm 2, at each restart $t_{j}$.

Algorithm 2 allows us to find a sparse quadrature rule $\{w_{j}^{p},u_{j}^{p}\}_{p=1}^{Q_{u_{j}}}$ with the number of nodes $Q_{u_{j}}$ satisfying $Q_{u_{j}}\leq M\ll\tilde{Q}_{u_{j}}$. Thus, the number of nodes can be made independent of time and frequent restarts can be used in Algorithm 1. A one way to find a vector $z$ in the nullspace of $A$ is to compute a basis for the nullspace. In the implementation, we make use of the QR method to quickly select a vector at each iteration. Finally, we note that this routine does not necessarily yield nonnegative weights.

An application of this procedure to the two dimensional Gaussian random variable discussed before in Table 1 reduces the size of the quadrature rule from $500$ to $45$ while the accuracy remains almost the same as $3.19$E-$14$.

5 Numerical Experiments

In this section, we present several numerical simulations of our algorithm using low-dimensional nonlinear SDE models. We refer to [6, 39, 40] for details on the following numerical examples.

For the rest of the section, $T\in\mathbb{R}$ stands for the endpoint of the time interval while $\Delta t=T/n$ denotes the time-step after which restarts occur at $t_{j}=j\Delta t$. Furthermore, following [28, 22, 39, 40], we choose the following orthonormal basis for $L^{2}[t_{j},t_{j+1}]$:

[TABLE]

where $t\in[t_{j},t_{j+1}]$. We utilize either a first- or a second-order time integration method.

In our numerical examples, the dynamics converge to an invariant measure. To demonstrate the convergence behavior of our algorithm, we compare our results to the second order exact statistics or Monte Carlo simulations with sufficiently high sampling rate $M_{samp}$ where the exact solution is not available. Comparisons involve the following relative pointwise errors:

[TABLE]

where $\mu$ and $\sigma^{2}$ represent the mean and the variance, respectively. In some cases, we give estimations of the first six cumulants of the invariant measure.

Example 5.1**.**

As a first example we consider an OU process

[TABLE]

where the damping parameter is random and uniformly distributed in $[1,3]$, i.e. $b_{u}\sim U(1,3)$. This is an example of 2-dimensional non-Gaussian dynamics that may be seen as a coupled system for $(u,b_{u})$ with $db_{u}=0$.

For the first simulation, we consider the time domain $[0,4]$. The mean-return parameter $\mu_{u}$ is set to be $0.2$. The initial condition is normally distributed u_{0}\sim N(1,0.04)\rotatebox[origin={c}]{90.0}{\models}W(t) and $\sigma_{u}=4$. We use the Gauss-Hermite rule for the initial condition with the level $\lambda_{u_{0}}=3$, whereas for the damping parameter $b_{u}$, we use the Gauss-Legendre quadrature rule with a varying $\lambda_{b_{u}}$. For Brownian motion, we use $2$-dimensional approximation with the product Gauss-Hermite rule of the level $\lambda_{\boldsymbol{\xi}}=2$. We also take the set $T_{\boldsymbol{\alpha}}(u)$ as the normalized Hermite polynomials. For time integration, we utilize a second-order weak Runge-Kutta method with $\Delta\tau=5$E-$4$.

In Figure 2, we compare second order statistics obtained by our method to the exact solutions using $N=1,2$ and $3$, $\lambda_{b_{u}}=2,4,$ and $8$, and $\Delta t=0.4,0.2,0.1,$ and $0.05$. The results are obtained by calculating the moments in $u$ variable and then taking averaging with respect to the measure of $b_{u}$. The rows of the figure correspond to $N$-, $\lambda_{b_{u}}$-, and $\Delta t$-convergence of the method while keeping the other two degrees of freedom constant, respectively. For each different restart step $\Delta t$, we keep the time-integration step $\Delta\tau$ the same.

$N$-convergence of the method, Figure 3a and Figure 3b, shows that first two moments can be captured accurately with $N=2$. The level of accuracies are $O(10^{-9})$ and $O(10^{-5})$ for the mean and the variance, respectively. We notice that using a larger number quadrature level $\lambda_{b_{u}}$ and more frequent restarts also help to reduce the relative errors. We also observe that the convergence behavior of the method in terms of the size of time interval is at least quadratic; i.e. $O(\Delta t^{2})$.

Table 2 exhibits the first six cumulants of the equation (5.1) with $\mu_{u}=0.0$ in the long-time $T=8$. The limiting stationary measure can be obtained by solving the corresponding standard Fokker-Planck equation; see [6, 39, 24, 38]. The first $6$ cumulants, $\kappa_{i}$, $i=1,\ldots 6$, are obtained by computing moments in the solution variable and then averaging with respect to the damping parameter. The table shows that increasing the degree $N$ of polynomials in the constraints in (4.3) clearly helps to accurately capture the higher cumulants in the long-time. We also note that the approximations for the cumulants $\kappa_{5}$ and $\kappa_{6}$ become accurate when $N=6$ is used while they are inaccurate for $N=4$. This type of convergence behavior is related to the fact that the equation is linear in this case.

Example 5.2**.**

We now consider a nonlinearity in the equation so that the damping term includes a cubic term:

[TABLE]

We take $T=4$, $\Delta t=0.04$, and $\sigma_{u}=2$. For Brownian motion we use $K=2$ dimensional $\boldsymbol{\xi}$ vector with a sparse Gauss-Hermite rule with the level $\lambda_{\boldsymbol{\xi}}=3$. We also take the set $T_{\boldsymbol{\alpha}}(u)$ as Hermite polynomials.

In this case, we observe from Table 3 that the Algorithm 1 can be executed accurately in the long-time by increasing the number of degrees of freedom $N$. Comparing Table 3 and Table 2, we see that the accuracy in the higher cumulants are slightly decreased and it is harder to capture higher moments as this is a more complicated dynamics having a nonlinearity.

Example 5.3**.**

In this example, we consider a multiplicative noise term and the Cox-IngerSoll-Ross (CIR) model

[TABLE]

which is used to describe the evolution of interest rates [7]. The same model is also used in the Heston model to model random volatility. We impose the condition $2b_{u}\mu_{u}\geq\sigma_{u}^{2}$ so that the process stays positive. The process has a stationary distribution in the long time and the second order statistics can be computed analytically. Note also that the noise amplitude is non-Lipschitz.

We apply the first-order Milstein method [24]

[TABLE]

and then approximate Brownian motion increments by their finite-dimensional approximations. Since the noise term involves a function of $u$, it is necessary to use a time-integration method which takes this dependence into account. We set $T=3,u_{0}=1$, $\mu_{u}=0.6$, $\sigma_{u}=0.5$, $b_{u}=2$, $K=2$, $\lambda_{\xi}=4$, and $\Delta t=0.1$. We compute the associated orthogonal polynomials when we solve the optimization procedure.

The first row of Figure 3 shows the evolution of the mean and the variance of the model (5.2) obtained by the analytical solution and Algorithm 1 with $N=4$. We observe that in the long time, the statistics become stationary and are captured well by the algorithm. The second row shows that the algorithm can be performed accurately in the long-time with a relatively small noise-level $\sigma_{u}=0.5$.

In Table 4, we compare the DSGC algorithm and MC method in terms of accuracy and timing. For a fair comparison and to minimize the error of time integration, both methods use Milstein’s method with the time step $\Delta\tau=1$E-$4$. We compare the convergence behaviors of the variance using $N$ and $M_{samp}$ for DSGC and MC, respectively. To assess the robustness of the methods, we take $b_{u}=4$ and consider different magnitudes of the noise: $\sigma_{u}=0.5,1,$ and $2$. Since MC estimates are noisy, we repeat each MC simulation $20$ times and average the errors. We first observe that the accuracy of both methods drops when $\sigma_{u}$ is increased. Also, if a high accuracy is needed, DSGC algorithm can be executed faster than MC in this scenario, e.g., for $\sigma_{u}=1$ and $\sigma_{u}=2$, MC should use at least $M_{samp}=128$E+$4$ to get to the same level of accuracy of DSGC with $N=4$, which implies that DSGC is $20$ times faster. We also see that DSGC is very efficient for small magnitude of noise $\sigma_{u}=0.5$. The elapsed times for DSGC seem to scale quadratically with $N$. Note also that the convergence of MC method is guaranteed, however, although we observe numerical convergence in DSGC, we do not have a convergence rigorous result for the approximate solution of (5.3). Similar settings in [41, 42, 29] establish convergence results for SDEs with sufficiently smooth coefficients, which do not apply in this case.

Example 5.4**.**

This example concerns the following $2$-dimensional nonlinear system

[TABLE]

where $a_{u}\geq 0$, $b_{u},b_{v}>0$ are damping parameters, $\sigma_{u},\sigma_{v}>0$ are constants, and $W_{u}$ and $W_{v}$ are two real independent Brownian motions. The second equation is an Ornstein-Uhlenbeck process and it acts as a multiplicative noise in the first equation. Depending on the regime of the parameters, the dynamics of the solution $u$ exhibit intermittent non-Gaussian behavior; see detailed discussions in [13, 6].

Following [6], we consider the system parameters as $a_{u}=1$ , $b_{u}=1.2$, $b_{v}=0.5$,$\sigma_{u}=0.5$, $\sigma_{v}=0.5$ and take the initials u_{0}=N(1,\sigma_{u}^{2}/8b_{u})\,\rotatebox[origin={c}]{90.0}{\models}\,v_{0}=N(0,{\sigma_{v}^{2}/8b_{v}}). In this regime, the dynamics of $u$ are characterized by unstable bursts of large amplitude [6].

We consider a long-time $T=8$ and use $K=4$ variables for the random forcing terms with the quadrature level $\lambda_{\boldsymbol{\xi}}=2$. A sparse Gauss-Hermite quadrature rule is used for the random initial conditions. Since we use a small quadrature level $\lambda_{\boldsymbol{\xi}}$, we restart frequently and take $\Delta t=0.02$ in the following simulations. The ability to use small degrees of freedom with frequent restarts is one of the main advantages of the method.

The first simulation concerns the choice of the polynomials $T_{\boldsymbol{\alpha}}(u,v)$. Figure 4 shows the condition numbers of the constraint matrix $A$ in (4.3) for different choices of polynomials and varying $N$. We observe that in each case, the lowest condition numbers correspond to ones which are obtained by orthogonal polynomials with respect to the joint distribution of $(u,v)$. Although the computation of orthogonal polynomials for large degree of polynomials is unstable, once the computation is carried out, it yields a well-conditioned constraint matrix for fixed $N$. Moreover, increasing the degree of polynomials $N$ leads to overall larger condition numbers, which might affect the stability of the numerical minimization procedures for large $N$.

Using polynomials $T_{\boldsymbol{\alpha}}(u,v)$ which are orthogonal with respect to the joint distribution of the solution, we next demonstrate $N$-convergence of the method in terms of the mean and the variance in Figure 5. We observe from the figure that the mean and the variance can be captured in the long-time up to $O(10^{-3})$ of accuracy. Throughout the time evolution, the number of quadrature points of the joint distribution of $(u,v)$ is ${N+2\choose 2}$. For instance, for $N=5,6,7$, and $8$, this number becomes $21,28,36$, and $45$, respectively. It is useful to note that although increasing $N$ gives better errors, it also makes the numerical optimization procedure relatively unstable.

Finally, we take a relatively short time $T=1$, and compare the computational times and the accuracy of DSGC, DgPC ([40]), and MC methods. All methods use a second order time-integration method with the same time step $\Delta\tau=2$E-$4$. MC simulations are repeated 20 times to get a stable estimate of the error. Note that it is getting computationally harder to obtain a stable estimate in MC for longer times since the estimates are noisy. Both DSGC and DgPC use the restart step $\Delta t=0.05$.

From Table 5, we first notice that for the same level of accuracy, DSGC performs better than both methods in terms of computational time. This behavior is due to the fact that the method can utilize a small number of samples and frequent restarts to achieve a high accuracy in a fast manner. We also observe that DgPC with $N=2$ offers a sufficient accuracy for many applications with a longer computational time in this scenario. Thus, DSGC offers a relatively fast way to compute moments of the solution in this low-dimensional case. Finally, we note that although DSGC is faster than DgPC in this case, it has been shown in [40] that DgPC offers a faster way to compute the solutions of SPDEs compared to MC methods because sampling methods suffer from high number of samples and complicated dynamics. In that case, DSGC might have an advantage over MC methods since it can use small number of samples by leveraging the regularity of the solution; see also concluding remarks about extensions to high dimensional cases.

6 Conclusion

We presented a collocation method, called Dynamical Sparse Grid Collocation method, which exploits the Markovian structure of SDEs to keep the degrees of freedom under control and obtain accurate long-time computations. The method uses a restart mechanism to construct a sparse quadrature rule for the solution on-the-fly and incorporates that rule with quadrature rules for the random forcing to capture the solution at later times. This fundamental idea is very similar to those of algorithms presented in our earlier works [39, 40] with one major difference: the current algorithm is non-intrusive. Being a non-intrusive method, one can leverage optimized legacy differential equation solvers and carry out evolution of the solution particles in parallel.

The main computational difficulties are how to keep the number of quadrature nodes independent of time, and to compute a sparse quadrature rule for the solution variable in a stable and efficient manner on-the-fly. We made use of a $L^{1}$-minimization procedure with the constraints that the quadrature rule integrates polynomials in the solution variable up to certain degrees of freedom exactly. We also discussed how to extract a quadrature rule with a few nodes and presented several different polynomials bases used in the optimization. Numerical results for different nonlinear low-dimensional dynamics confirmed the ability of the algorithm reaching long times.

For high-dimensional SDEs or SPDEs, the algorithm needs further modifications. Dynamical sparse quadrature rules for the solution will have higher number of nodes, which, in turn, may hinder the efficiency of the algorithm; especially in case of complex dynamics. Thus, further dimensionality reduction techniques and simple parallelization can be considered. Moreover, the stability of the constrained optimization in high-dimensional cases should also be investigated. We demonstrated such an extension from SDEs to SPDEs in our previous works [39, 40]. The DgPC method was coupled with Karhunen–Loeve expansion and a few dominating random modes at each restart were selected and incorporated into a PCE to represent future solutions. A similar extension can be done for DSGC by using the KLE of random solutions in time and constructing sparse quadrature rules for the selected finite number of KLE random modes.