We present a method to learn mean residence time and escape probability from data modeled by stochastic differential equations. This method is a combination of machine learning from data (to extract stochastic differential equations as models) and stochastic dynamics (to quantify dynamical behaviors with deterministic tools). The goal is to learn and understand stochastic dynamics based on data. This method is applicable to sample path data collected from complex systems, as long as these systems can be modeled as stochastic differential equations.

Stochastic dynamical systems are appropriate models for randomly influenced systems. Understanding the complex dynamical behaviors of these systems is a challenge in diverse areas of science and engineering. In deterministic dynamical systems, invariant manifolds and other invariant structures provide global information for dynamical evolution. For stochastic dynamical systems, better quantitative analysis and understanding is needed because of limitations in our current analytical skills or computation capability. Fortunately, researchers are increasingly using data-driven methods for system identification and the discovery of dynamics. In this work, we propose a new approach to determine some computable dynamical quantities from data, such as the mean residence time and escape probability, which offer insights into global dynamics under uncertainty. We demonstrate the algorithm to be effective and robust by reproducing known dynamics and evaluating errors for several prototypical stochastic dynamical systems with Brownian motions.

Stochastic dynamical systems arise in modeling molecular dynamics, mechanical and electrical engineering, climate dynamics, geophysical and environmental systems, among others. Advances in machine learning and data science are leading to new progresses in the analysis and understanding of complex dynamics for systems with massive observation data sets. Despite the rapid development of tools to extract governing equations from data, there has been slow progress in distilling quantities that may be used to explore stochastic dynamics. It is desirable to extract deterministic quantities that carry dynamical information of stochastic differential equations (SDEs). These deterministic quantities include moments for solution paths, probability density functions for solution paths, mean residence time, and escape probability.7,9 These concepts help us to understand various phenomena in complex systems under uncertainty.1 

In this present paper, we will present a new approach to extract the underlying deterministic quantities, mean residence time, and escape probability that describe certain aspects of stochastic dynamics. In fact, mean residence time for a stochastic dynamical system quantifies how long the system stays in a region, and escape probability describes the likelihood of a system transition from one regime to another. Fortunately, these deterministic quantities can be determined by solving an elliptic partial differential equation as in Duan,9 once the underlying stochastic differential equation model for the system evolution is discovered from data.

We assume that a data set is composed of sample paths governed by a stochastic differential equation (SDE) in Rn,

dXt=b(Xt)dt+σ(Xt)dBt,
(1)

where b is an n-dimensional vector function, σ is an n×m matrix function, and Bt is an m-dimensional Brownian motion. This is the customary, probabilistic way of writing the equation

dXtdt=b(Xt)+σ(Xt)dBtdt.

Often, b is called “drift” and σ is called “diffusion.” Assume that b and σ satisfy an appropriate local Lipschitz condition as follows:

b(x)b(y)+σ(x)σ(y)KNxy,

for xN, yN and N>0. Here, the Lipschitz constant KN depends on the positive number N. The generator for this SDE system is a linear second-order differential operator:

Ag=b(g)+12Tr[σσTH(g)],gH02(Rn),
(2)

where H denotes the Hessian matrix of a multivariate function and Tr denotes the trace of a matrix. We view A as a linear differential operator in Hilbert space L2(Rn) with domain of definition D(A)=H02(Rn). Furthermore, we assume that the generator A is uniformly elliptic. That is for xD and all ξRn, there exists a positive constant C such that

i,j=1n(σ(x)σT(x))i,jξiξjC|ξ|2.
(3)

We discuss mean residence time and escape probability as deterministic quantities that carry dynamical information for solution orbits of (1). For a bounded domain DRn (with boundary D), the first exit time for a solution orbit starting at xD is a stopping time defined as

τD(ω)inf{t>0:X0=x,XtD}.

The mean residence time is defined as u(x)EτD(ω). It is the mean residence time of a particle initially at x inside D until the particle first hits the boundary D or escapes from D. It turns out that the mean residence time u of stochastic system (1) can be determined by solving a deterministic partial differential equation as follows:

Au=1,
(4)
u|D=0,
(5)

where A is the generator defined as (2). To prove the existence and uniqueness of u, we recall the Dynkin’s formula (Theorem 7.4.1 in Oksendal17),

Ex[f(XτD)]=f(x)+Ex[0τDAf(Xs)ds],

for f in the domain of definition of the generator A, Ex is the expectation with respect to the probability law Qx induced by a solution process Xt starting at xRn. We take the continuous bounded boundary value ϕ=0 and the continuous inhomogeneous term g=1. By Theorem 9.3.3 in Oksendal,17 which is a consequence of Dynkin’s formula, we know that the linear expectation Ex[ϕ(XτD)]+Ex[0τDg(Xs)ds], which is just EτD=u(x), solves Au=1 with the boundary condition u|D=0. Thus, we can numerically compute mean residence time u by solving elliptic partial differential equations, so we know how long the stochastic system (1) stays in the region D. For more details, see Duan.9 

The escape probability is the likelihood that an orbit starting inside a domain D exits from this domain first through a specific part Γ of the boundary D. Let Γ be a subset of the boundary D. We define the escape probability p(x) from D through Γ as the likelihood that Xt starting at x exits from D first through p(x)=P{XτDΓ}. We will show that the escape probability p(x) solves a linear elliptic partial differential equation, with a specifically chosen Dirichlet condition as follows:

Ap=0,
(6)
p|Γ=1,
(7)
p|DΓ=0,
(8)

where A is the generator defined as (2). Taking

ϕ(x)={1,xΓ,0,xDΓ,

we have

E[ϕ(XτD(x))]={ω:XτD(x)Γ}ϕ(XτD(x))dP(ω)+{ω:XτD(x)DΓ}ϕ(XτD(x))dP(ω)=P{ω:XτD(x)Γ}=p(x).

This means that E[ϕ(XτD(x))] is the escape probability p(x) that we are looking for. We know Ap=0 together with p|Γ=1 and p|DΓ=0 by Theorem 9.2.14 in Okendal.17 For more details, see Duan.9 Moreover, we suppose that solution orbits, i.e., “particles” are initially uniformly distributed in D, then the average escape probability P that a trajectory will leave D through Γ is P=1DDp(x)dx.

To obtain the mean residence time u and escape probability p from sample path data XTrue [with the underlying model (1)], we propose the following machine learning algorithm: First, we collect sample path data XTrue by sample-wisely simulating (1) via the Euler method (and treat these data as our observation data), then we try to learn the stochastic dynamical system model from these data with the following model ansatz:

X˙=ΘΞ,
(9)

where the basis Θ consists of polynomial functions {1,x,y,,xNyN} together with the (generalized) time derivative of Brownian motion dBt/dt, and Ξ is the coefficient (or coordinate, or weight) under this basis. Here, denotes scalar product. Each column ξk=[ξ1k,ξ2k,,ξNk] of Ξ is a vector of coefficients, determining which terms are active in the right-hand side for one of the row equations in (1). Brunton et al.20 considered such a learned model, with a basis not containing noise (i.e., dBt/dt).

Furthermore, we set up a regression problem to determine coefficient Ξ=[ξ1,,ξn] by minimizing the mean-square discrepancy (other metrics are possible) between the data XTrue (many samples) and the solution XLearn (many samples) for the learned (i.e., extracted) governing model (9). This will provide us the “learned” drift b and “learned” diffusion σ, and thus, we also have the extracted stochastic model (9), together with the “learned” generator A.

Finally, with this learned generator A, we compute mean residence time u and escape probability p, by solving the deterministic partial differential Eqs. (4)(5) and (6)(8), respectively. We verify that our algorithm is effective by estimating the error of the maximal mean residence time

{error}(u)=maxDuLearnuTrue

and the error of the average escape probability

{error}(P)=PLearnPTrue

between the learned system (9) and the original system (1).

This paper is organized as follows. In Sec. II, we apply our method to learn a two-dimensional quasigeostrophic meandering jet model with additive noise and multiplicative noise and compute the mean residence time and escape probability. In Sec. III, we illustrate our method to learn two three-dimensional systems (a linear damped oscillator and the well-known Lorenz system). We summarize and conclude in Sec. IV.

Zonal shear flows occur naturally in both oceans and the atmosphere. The Gulf Stream is a well-known example. Based on RAFOS floats observations of the Gulf Stream, Bower et al.4,5 viewed the fluid motion as a steady and eastward meander propagation in the moving frame and divided the velocity field into three regimes: a central jet, exterior retrograde motion, and intermediate closed circulations above meander troughs and below crests. There was no exchange occurring between the three regimes in this model. Fluid particles in the intermediate regime execute periodic motion but never escape. However, float observations reveal that several particle trajectories pass through meander crests and troughs and then leave jet, which indicate that exchange does occur across some part of the Gulf Stream. Samelson15 modified the basic model and considered three different types of variability: a time-dependent spatially uniform meridional velocity superimposed on the basic flow, a time-dependent meander amplitude, and a propagating plane wave superimposed on the basic flow. In particular, he only considered time-periodic variabilities. Subsequently, Samelson,16 del-Castillo-Negrete et al.,8 Pratt et al.,11 Beigie et al.,3 Duan et al.,6 and Brannan et al.7 have obtained a series of outcomes from the point of view of dynamical systems.

As an approximate solution to the quasigeostrophic model for two-dimensional geophysical flows,13 the basic Bickley jet8 is

ψ(x,y)=tanh(y)+asech2(y)cos(kx)+cy,

where a=0.01, c=13(1+132β), and k=6c, 0β23. In this paper, we take β=13. Note β=2Ωrcosθ is the meridional derivative of the Coriolis parameter, where Ω is the rotation rate of the earth, r is the earth’s radius, and θ is the latitude. This stream function ψ(x,y) defines the basic meandering jet system x˙=ψy, y˙=ψx. The phase portrait for this deterministic system is in Fig. 1 (left). As it shows no fluid exchange between the eddies and the jet, it is not an appropriate model for a meandering jet.

FIG. 1.

Lagrangian fluid trajectories of the original meandering jet system (10)(13) with 20 initial positions. The fluid motion in the deterministic basic model (σ=0) (left) is steady, with no exchange of fluid between the eddies and the jet stream. Under additive noise (σ=0.3) (middle) and multiplicative noise (σ=0.3) (right), there are fluid exchanges between eddies and the jet stream.

FIG. 1.

Lagrangian fluid trajectories of the original meandering jet system (10)(13) with 20 initial positions. The fluid motion in the deterministic basic model (σ=0) (left) is steady, with no exchange of fluid between the eddies and the jet stream. Under additive noise (σ=0.3) (middle) and multiplicative noise (σ=0.3) (right), there are fluid exchanges between eddies and the jet stream.

Close modal

Thus, we incorporate random wind forcing and other fluctuations, either additive or multiplicative, in the model for the meandering jet, and consequently, we consider the following two stochastic dynamical systems:

dx=ψydt+σdB1,
(10)
dy=ψxdt+σdB2
(11)

and

dx=ψydt+σxdB1,
(12)
dy=ψxdt+σydB2,
(13)

where the noise intensity σ satisfies 0<σ<1, and B1(t), B2(t) are two independent Brownian motions. For both stochastic systems, a number of sample solution orbits are shown in Fig. 1 (middle, right). The central jet and two rows of recirculation eddies, which are called the northern and southern recirculation regions, are still visible. Outside the recirculation regimes are the exterior retrograde regimes. There is no exchange between regimes in unperturbed model, but exchange does occur in other two models with additive or multiplicative noises. These two stochastic models are more appropriate for modeling a meandering jet (as a simplified model for the Gulf Stream), as observations indicate fluid exchange.

A wide range of dynamical systems with intrinsic noise may be modeled with stochastic differential equations. However, identifying an appropriate SDE model from intermittent observations of the system is challenging, particularly if the dynamical process is nonlinear and the observations are noisy and indirect.10 Brunton et al.,20 Zhang et al.,18 and Dunker et al.21 presented some approaches to discover or learn governing physical laws from data, with underlying differential equations models.

We use our machine learning algorithm to data, from systems (10) and (11) or systems (12) and (13), to learn their governing laws in terms of the following model:

x˙=Θξ1,
(14)
y˙=Θξ2,
(15)

where Θ is a set of basis functions and Ξ=[ξ1,ξ2] defines the coefficients (or weights). To achieve this, the idea is as follows: First, we collect a time series of the system state (x(ti),y(ti)),i=1,,20000, from systems (10) and (11) or systems (12) and (13) with noise intensity σ=0.3, numerically by the Euler method [as our observation data (xTruei,yTruei)]. Then, we construct a library Θ consisting of basis functions with polynomials {1,x,y,,y5} and time derivative of Brownian motions {dB1/dt,dB2/dt}. We determine each column of coefficients ξk=[ξ1k,ξ2k,,ξ23k],k=1,2, by minimizing discrepancies Ei=120000xTrueixLearni2 and Ei=120000yTrueiyLearni2. Here, (xLearni,yLearni) is the numerical solution of the learned meandering jet model (14) and (15) by the Euler method (see Tables I and II).

TABLE I.

Identified coefficients for the learned meandering jet systems (14) and (15). Data are numerically generated via Euler method for (10) and (11) (additive noise) with initial position (− 0.2, 0.8)T in the eddy, for time t ∈ [0, 200] with stepsize 0.01.

Basisx˙y˙
2.211 465 0 × 10+00 1.135 687 6 × 10+00 
x 9.167 092 9 × 10−01 −3.105 340 6 × 10−01 
y −8.533 715 6 × 10+00 −7.883 008 9 × 10+00 
x2 2.158 193 3 × 10−01 −5.888 226 7 × 10−02 
xy −3.540 679 3 × 10+00 1.208 688 5 × 10+00 
y2 1.480 699 0 × 10+01 2.150 933 2 × 10+01 
x3 2.316 303 0 × 10−02 −2.068 737 6 × 10−03 
x2y −6.440 459 8 × 10−01 1.655 640 0 × 10−02 
xy2 5.216 326 3 × 10+00 −2.304 683 1 × 10+00 
y3 −1.356 879 2 × 10+01 −2.915 997 3 × 10+01 
x4 6.686 967 2 × 10−04 1.118 142 7 × 10−02 
x3y −7.380 868 5 × 10−02 6.441 070 0 × 10−02 
x2y2 5.297 775 4 × 10−01 2.256 811 7 × 10−01 
xy3 −3.578 395 4 × 10+00 2.294 521 4 × 10+00 
y4 5.631 493 3 × 10+00 1.977 609 2 × 10+01 
x5 2.851 040 0 × 10−04 1.238 256 6 × 10−03 
x4y −8.195 872 2 × 10−04 −8.028 851 6 × 10−04 
x3y2 4.227 503 3 × 10−02 −3.688 687 4 × 10−02 
x2y3 −1.161 596 8 × 10−01 −1.736 499 4 × 10−01 
xy4 9.944 926 4 × 10−01 −9.095 955 7 × 10−01 
y5 −6.810 078 8 × 10−01 −5.387 695 2 × 10+00 
dB1/dt 5.477 303 0 × 10−01 −9.739 278 8 × 10−06 
dB2/dt 2.707 655 1 × 10−07 5.477 183 7 × 10−01 
Basisx˙y˙
2.211 465 0 × 10+00 1.135 687 6 × 10+00 
x 9.167 092 9 × 10−01 −3.105 340 6 × 10−01 
y −8.533 715 6 × 10+00 −7.883 008 9 × 10+00 
x2 2.158 193 3 × 10−01 −5.888 226 7 × 10−02 
xy −3.540 679 3 × 10+00 1.208 688 5 × 10+00 
y2 1.480 699 0 × 10+01 2.150 933 2 × 10+01 
x3 2.316 303 0 × 10−02 −2.068 737 6 × 10−03 
x2y −6.440 459 8 × 10−01 1.655 640 0 × 10−02 
xy2 5.216 326 3 × 10+00 −2.304 683 1 × 10+00 
y3 −1.356 879 2 × 10+01 −2.915 997 3 × 10+01 
x4 6.686 967 2 × 10−04 1.118 142 7 × 10−02 
x3y −7.380 868 5 × 10−02 6.441 070 0 × 10−02 
x2y2 5.297 775 4 × 10−01 2.256 811 7 × 10−01 
xy3 −3.578 395 4 × 10+00 2.294 521 4 × 10+00 
y4 5.631 493 3 × 10+00 1.977 609 2 × 10+01 
x5 2.851 040 0 × 10−04 1.238 256 6 × 10−03 
x4y −8.195 872 2 × 10−04 −8.028 851 6 × 10−04 
x3y2 4.227 503 3 × 10−02 −3.688 687 4 × 10−02 
x2y3 −1.161 596 8 × 10−01 −1.736 499 4 × 10−01 
xy4 9.944 926 4 × 10−01 −9.095 955 7 × 10−01 
y5 −6.810 078 8 × 10−01 −5.387 695 2 × 10+00 
dB1/dt 5.477 303 0 × 10−01 −9.739 278 8 × 10−06 
dB2/dt 2.707 655 1 × 10−07 5.477 183 7 × 10−01 
TABLE II.

Identified coefficients for the learned meandering jet system (14) and (15). Data are numerically generated via the Euler method for (12) and (13) (multiplicative noise) with initial position (− 0.2, 0.8)T in the eddy, for time t ∈ [0, 200] with stepsize 0.01.

Basisx˙y˙
4.679 158 4 × 10+00 −7.586 793 5 × 10+00 
x 6.300 678 4 × 10−01 −2.108 732 6 × 10+00 
y −2.548 358 6 × 10+01 4.333 514 1 × 10+01 
x2 2.610 926 3 × 10−01 −1.595 241 1 × 10−01 
xy −1.742 074 1 × 10+00 1.056 439 5 × 10+01 
y2 6.116 905 8 × 10+01 −9.472 222 1 × 10+01 
x3 −3.182 419 0 × 10−03 −1.005 654 0 × 10−01 
x2y −9.661 535 4 × 10−01 −2.467 371 8 × 10−01 
xy2 8.853 546 5 × 10−01 −2.153 239 1 × 10+01 
y3 −7.682 962 4 × 10+01 9.686 362 9 × 10+01 
x4 −1.734 162 0 × 10−04 4.713 976 6 × 10−03 
x3y −6.630 334 9 × 10−03 2.637 329 7 × 10−01 
x2y2 1.143 603 0 × 10+00 1.300 516 0 × 10+00 
xy3 1.074 806 0 × 10+00 2.069 413 0 × 10+01 
y4 4.875 646 5 × 10+01 −4.404 292 2 × 10+01 
x5 −1.691 796 4 × 10−04 2.235 049 5 × 10−04 
x4y −3.953 219 7 × 10−03 −5.468 421 9 × 10−03 
x3y2 −1.339 436 5 × 10−02 −1.884 372 4 × 10−01 
x2y3 −4.833 223 9 × 10−01 −1.007 315 1 × 10+00 
xy4 −8.735 959 0 × 10−01 −7.798 103 5 × 10+00 
y5 −1.244 311 4 × 10+01 6.056 436 1 × 10+00 
xdB1/dt 5.477 237 4 × 10−01 4.326 505 7 × 10−05 
ydB2/dt 1.682 139 3 × 10−05 5.453 187 5 × 10−01 
Basisx˙y˙
4.679 158 4 × 10+00 −7.586 793 5 × 10+00 
x 6.300 678 4 × 10−01 −2.108 732 6 × 10+00 
y −2.548 358 6 × 10+01 4.333 514 1 × 10+01 
x2 2.610 926 3 × 10−01 −1.595 241 1 × 10−01 
xy −1.742 074 1 × 10+00 1.056 439 5 × 10+01 
y2 6.116 905 8 × 10+01 −9.472 222 1 × 10+01 
x3 −3.182 419 0 × 10−03 −1.005 654 0 × 10−01 
x2y −9.661 535 4 × 10−01 −2.467 371 8 × 10−01 
xy2 8.853 546 5 × 10−01 −2.153 239 1 × 10+01 
y3 −7.682 962 4 × 10+01 9.686 362 9 × 10+01 
x4 −1.734 162 0 × 10−04 4.713 976 6 × 10−03 
x3y −6.630 334 9 × 10−03 2.637 329 7 × 10−01 
x2y2 1.143 603 0 × 10+00 1.300 516 0 × 10+00 
xy3 1.074 806 0 × 10+00 2.069 413 0 × 10+01 
y4 4.875 646 5 × 10+01 −4.404 292 2 × 10+01 
x5 −1.691 796 4 × 10−04 2.235 049 5 × 10−04 
x4y −3.953 219 7 × 10−03 −5.468 421 9 × 10−03 
x3y2 −1.339 436 5 × 10−02 −1.884 372 4 × 10−01 
x2y3 −4.833 223 9 × 10−01 −1.007 315 1 × 10+00 
xy4 −8.735 959 0 × 10−01 −7.798 103 5 × 10+00 
y5 −1.244 311 4 × 10+01 6.056 436 1 × 10+00 
xdB1/dt 5.477 237 4 × 10−01 4.326 505 7 × 10−05 
ydB2/dt 1.682 139 3 × 10−05 5.453 187 5 × 10−01 

Figure 2 shows Lagrangian fluid trajectories of the learned meandering jet model (14) and (15), without noise, with additive noise, and with multiplicative noise.

FIG. 2.

Lagrangian fluid trajectories of the learned meandering jet model (14) and (15) with 20 initial positions: The fluid motion in the deterministic basic model when noise is absent (left) is steady, with no exchange of fluid between the eddies and the jet stream. Under additive noise (middle) and multiplicative noise (right), there are fluid exchanges between eddies and the jet stream.

FIG. 2.

Lagrangian fluid trajectories of the learned meandering jet model (14) and (15) with 20 initial positions: The fluid motion in the deterministic basic model when noise is absent (left) is steady, with no exchange of fluid between the eddies and the jet stream. Under additive noise (middle) and multiplicative noise (right), there are fluid exchanges between eddies and the jet stream.

Close modal

In this subsection, we will consider stochastic quasigeostrophic meandering jet models, (10)(11) and (12)(13), to discover mean residence time by our approach. We take an eddy (Fig. 3) as the bounded domain D (with boundary D composed with meander trough and crest).

FIG. 3.

An eddy: σ=0.

FIG. 3.

An eddy: σ=0.

Close modal

The mean residence time u(x,y) of the stochastic system (10) and (11) with additive noise for a trajectory starting in the eddy satisfies the following elliptic partial differential equation as in (4) and (5):

12σ2Δu+(ψy)ux+(ψx)uy=1,
(16)
u|D=0.
(17)

The mean residence time, as solution of this elliptic partial differential equation, is denoted by uTrue. For the corresponding learned meandering jet model (14) and (15), with learned coefficients in Table I, we can also set up a similar elliptic partial differential equation for learned mean residence time uLearn.

Similarly, the mean residence time u(x,y) of the stochastic system (12) and (13) with multiplicative noise, for a trajectory starting in the eddy, satisfies the following elliptic partial differential equation as in (4) and (5):

12σ2x2uxx+12σ2y2uyy+(ψy)ux+(ψx)uy=1,
(18)
u|D=0.
(19)

The mean residence time, as solution of this elliptic partial differential equation, is denoted by uTrue. For the corresponding learned meandering jet model (14) and (15), with learned coefficients in Table II, we can also set up a similar elliptic partial differential equation for learned mean residence time uLearn.

We use a finite element code to solve the preceding elliptic differential equations to get mean residence time uTrue and uLearn, for both additive and multiplicative noise cases. Figure 4 shows the mean residence time uLearn, as uTrue is barely distinguishable from uLearn when they are plotted together. Thus, we compare the error between the maximal values of the mean residence times uTrue and uLearn. In fact, the error of mean residence time between the learned system and the original system is {error}(u)=max(x,y)DuLearnuTrue=2.2142×104 for additive noise case, and {error}(u)=max(x,y)DuLearnuTrue=1.3×103 for multiplicative noise case.

FIG. 4.

Learned mean residence time uLearn for stochastic system (14) and (15) in an eddy. Additive noise case is on the left, and multiplicative noise case is on the right.

FIG. 4.

Learned mean residence time uLearn for stochastic system (14) and (15) in an eddy. Additive noise case is on the left, and multiplicative noise case is on the right.

Close modal

We again take the eddy (Fig. 3) to be the bounded domain D [with boundary D composed with meander trough (lower sub-boundary) and crest (upper sub-boundary)]. For the stochastic meandering jet system (10) and (11) with additive noise, the escape probability p(x,y) of a fluid particle, starting at (x,y) in an eddy and escaping through a boundary component Γ (either the trough or crest), satisfies the elliptic partial differential equations (6)(8),

12σ2Δp+(ψy)px+(ψx)py=0,
(20)
p|Γ=1,
(21)
p|DΓ=0.
(22)

For the corresponding learned meandering jet model (14) and (15), with learned coefficients in Table I, we can also set up a similar elliptic partial differential equation for the learned escape probability pLearn.

Similarly, for the stochastic meandering jet system (12) and (13) with multiplicative noise, the escape probability p(x,y) of a fluid particle, starting at (x,y) in an eddy and escaping through a boundary component Γ (either the trough or crest), satisfies the elliptic partial differential Eqs. (6)(8),

12σ2x2pxx+12σ2y2pyy+(ψy)px+(ψx)py=0,
(23)
p|Γ=1,
(24)
p|DΓ=0.
(25)

For the corresponding learned meandering jet model (14) and (15), with learned coefficients in Table II, we can also set up a similar elliptic partial differential equation for the learned escape probability pLearn.

Figures 5 and 6 show the finite element solution for the learned escape probability pLearn alone, as pTrue is barely distinguishable from pLearn when they are plotted together. Thus, we compare the error between the average escape probability values.

FIG. 5.

Learned escape probability PLearn for stochastic system (14) and (15) in an eddy, exiting from upper sub-boundary. Additive noise case is on the left, and multiplicative noise case is on the right.

FIG. 5.

Learned escape probability PLearn for stochastic system (14) and (15) in an eddy, exiting from upper sub-boundary. Additive noise case is on the left, and multiplicative noise case is on the right.

Close modal
FIG. 6.

Learned escape probability PLearn for stochastic system (14) and (15) in an eddy, exiting from lower sub-boundary. Additive noise case is on the left, and multiplicative noise case is on the right.

FIG. 6.

Learned escape probability PLearn for stochastic system (14) and (15) in an eddy, exiting from lower sub-boundary. Additive noise case is on the left, and multiplicative noise case is on the right.

Close modal

For fluid particles initially uniformly distributed in D (an eddy), we compute the average escape probability P for particles leaving D through the upper or lower sub-boundary Γ, given by P=1DDp(x,y)dxdy. The error of average escape probability, for fluid particles exiting from D through the upper sub-boundary, between the learned system and the original system is {error}(P)=PLearnPTrue=1.3949×105 both for additive noise and multiplicative noise cases. The error of the average escape probability, for particles leaving D through the lower sub-boundary, between the learned system and the original system is {error}(P)=PLearnPTrue=7.4133×104 both for additive noise and multiplicative noise cases.

In this section, we will illustrate our method in two systems with additive noise. The first system is for a three-dimensional stochastic linear oscillator. The second is the more complex stochastic Lorenz system. In these two examples, data from direct numerical simulations (as “observation data”) are used to discover mean residence time and escape probability.

We consider a linear stochastic system

dx=(0.1x2y)dt+ϵdB1,
(26)
dy=(2x0.1y)dt+ϵdB2,
(27)
dz=(0.3z)dt+ϵdB3,
(28)

where ϵ is the noise intensity (taken to be 0.9 here) and B1(t), B2(t), B3(t) are three independent Brownian motions. With sample path data for this system, we try to discover the governing equation,

x˙=Θξ1,
(29)
y˙=Θξ2,
(30)
z˙=Θξ3,
(31)

as in Brunton et al.,20 but we use a stochastic basis. Here, Θ is a set of basis functions and Ξ=[ξ1,ξ2,ξ3] are coefficients (or weights). As in Sec. II, we collect data (x(ti),y(ti),z(ti)),i=1,,20000, from samplewise simulations of the original system (26)(28), and construct a library Θ consisting of polynomial {1,x,y,z,,z4} and time derivatives of Brownian motions {dB1/dt,dB2/dt,dB3/dt}. We then solve a regression problem to determine the weights ξk=[ξ1k,,ξ38k],k=1,2,3 (see Table III and Fig. 7).

FIG. 7.

Trajectory of the three-dimensional stochastic linear system starting at (1,1,1)T: the trajectory for the original system in blue curve (left) and for the learned system in red curve (right).

FIG. 7.

Trajectory of the three-dimensional stochastic linear system starting at (1,1,1)T: the trajectory for the original system in blue curve (left) and for the learned system in red curve (right).

Close modal
TABLE III.

Identified coefficients for discovering a three-dimensional stochastic linear system with basis functions. Sample paths of (26)(28) with initial position (1, 1, 1)T are numerically generated, for t ∈ [0, 200] with stepsize 0.01.

Basisx˙y˙z˙
−3.901 872 6 × 10−06 −2.527 055 3 × 10−06 9.962 128 3 × 10−06 
x −9.997 259 3 × 10−02 2.000 017 8 × 10+00 −6.997 397 7 × 10−05 
y −1.999 946 3 × 10+00 −9.996 519 8 × 10−02 −1.371 972 5 × 10−04 
z 1.972 107 4 × 10−04 1.277 239 2 × 10−04 −3.005 035 1 × 10−01 
x2 4.195 535 5 × 10−05 2.717 246 7 × 10−05 −1.071 189 8 × 10−04 
xy −2.942 892 2 × 10−06 −1.905 969 8 × 10−06 7.513 692 0 × 10−06 
xz 1.225 505 7 × 10−04 7.937 011 3 × 10−05 −3.128 919 4 × 10−04 
y2 2.406 794 8 × 10−05 1.558 765 3 × 10−05 −6.144 946 4 × 10−05 
yz 3.387 555 3 × 10−04 2.193 956 7 × 10−04 −8.648 990 8 × 10−04 
z2 4.021 428 3 × 10−03 2.604 485 8 × 10−03 −1.026 737 4 × 10−02 
x3 −1.164 887 1 × 10−04 −7.544 413 7 × 10−05 2.974 150 0 × 10−04 
x2y −2.361 637 3 × 10−04 −1.529 519 0 × 10−04 6.029 651 9 × 10−04 
x2z −8.288 823 6 × 10−04 −5.368 272 6 × 10−04 2.116 274 2 × 10−03 
xy2 −1.177 680 6 × 10−04 −7.627 271 3 × 10−05 3.006 814 1 × 10−04 
xyz 3.517 902 6 × 10−05 2.278 376 4 × 10−05 −8.981 788 9 × 10−05 
xz2 −1.828 761 2 × 10−03 −1.184 400 7 × 10−03 4.669 130 7 × 10−03 
y3 −2.352 567 6 × 10−04 −1.523 645 0 × 10−04 6.006 495 4 × 10−04 
y2z −1.009 917 0 × 10−03 −6.540 746 8 × 10−04 2.578 485 6 × 10−03 
yz2 −6.201 471 5 × 10−03 −4.016 395 0 × 10−03 1.583 338 6 × 10−02 
z3 −1.242 848 7 × 10−02 −8.049 333 6 × 10−03 3.173 198 9 × 10−02 
x4 −9.374 085 3 × 10−05 −6.071 144 4 × 10−05 2.393 359 5 × 10−04 
x3y −2.102 290 2 × 10−06 −1.361 552 3 × 10−06 5.367 495 7 × 10−06 
x3z 7.505 776 6 × 10−04 4.861 130 7 × 10−04 −1.916 349 3 × 10−03 
x2y2 −1.271 197 6 × 10−04 −8.232 935 5 × 10−05 3.245 578 3 × 10−04 
x2yz 1.801 066 1 × 10−03 1.166 463 9 × 10−03 −4.598 420 6 × 10−03 
x2z2 2.218 039 3 × 10−03 1.436 517 4 × 10−03 −5.663 022 3 × 10−03 
xy3 −8.741 210 4 × 10−06 −5.661 261 8 × 10−06 2.231 776 1 × 10−05 
xy2z 7.652 722 2 × 10−04 4.956 300 3 × 10−04 −1.953 866 9 × 10−03 
xyz2 −3.053 939 8 × 10−05 −1.977 890 0 × 10−05 7.797 215 1 × 10−05 
xz3 4.838 899 2 × 10−05 3.133 922 4 × 10−05 −1.235 451 2 × 10−04 
y4 −3.262 586 8 × 10−05 −2.113 020 7 × 10−05 8.329 925 5 × 10−05 
y3z 1.799 062 5 × 10−03 1.165 166 3 × 10−03 −4.593 305 0 × 10−03 
y2z2 2.276 793 8 × 10−03 1.474 569 9 × 10−03 −5.813 032 3 × 10−03 
yz3 3.033 272 8 × 10−03 1.964 505 0 × 10−03 −7.744 448 8 × 10−03 
z4 5.742 310 1 × 10−03 3.719 018 2 × 10−03 −1.466 107 0 × 10−02 
dB1/dt 9.486 847 6 × 10−01 9.444 878 4 × 10−07 −3.723 349 0 × 10−06 
dB2/dt −1.808 827 7 × 10−05 9.486 715 8 × 10−01 4.618 237 2 × 10−05 
dB2/dt 4.246 885 2 × 10−05 2.750 503 4 × 10−05 9.485 748 7 × 10−01 
Basisx˙y˙z˙
−3.901 872 6 × 10−06 −2.527 055 3 × 10−06 9.962 128 3 × 10−06 
x −9.997 259 3 × 10−02 2.000 017 8 × 10+00 −6.997 397 7 × 10−05 
y −1.999 946 3 × 10+00 −9.996 519 8 × 10−02 −1.371 972 5 × 10−04 
z 1.972 107 4 × 10−04 1.277 239 2 × 10−04 −3.005 035 1 × 10−01 
x2 4.195 535 5 × 10−05 2.717 246 7 × 10−05 −1.071 189 8 × 10−04 
xy −2.942 892 2 × 10−06 −1.905 969 8 × 10−06 7.513 692 0 × 10−06 
xz 1.225 505 7 × 10−04 7.937 011 3 × 10−05 −3.128 919 4 × 10−04 
y2 2.406 794 8 × 10−05 1.558 765 3 × 10−05 −6.144 946 4 × 10−05 
yz 3.387 555 3 × 10−04 2.193 956 7 × 10−04 −8.648 990 8 × 10−04 
z2 4.021 428 3 × 10−03 2.604 485 8 × 10−03 −1.026 737 4 × 10−02 
x3 −1.164 887 1 × 10−04 −7.544 413 7 × 10−05 2.974 150 0 × 10−04 
x2y −2.361 637 3 × 10−04 −1.529 519 0 × 10−04 6.029 651 9 × 10−04 
x2z −8.288 823 6 × 10−04 −5.368 272 6 × 10−04 2.116 274 2 × 10−03 
xy2 −1.177 680 6 × 10−04 −7.627 271 3 × 10−05 3.006 814 1 × 10−04 
xyz 3.517 902 6 × 10−05 2.278 376 4 × 10−05 −8.981 788 9 × 10−05 
xz2 −1.828 761 2 × 10−03 −1.184 400 7 × 10−03 4.669 130 7 × 10−03 
y3 −2.352 567 6 × 10−04 −1.523 645 0 × 10−04 6.006 495 4 × 10−04 
y2z −1.009 917 0 × 10−03 −6.540 746 8 × 10−04 2.578 485 6 × 10−03 
yz2 −6.201 471 5 × 10−03 −4.016 395 0 × 10−03 1.583 338 6 × 10−02 
z3 −1.242 848 7 × 10−02 −8.049 333 6 × 10−03 3.173 198 9 × 10−02 
x4 −9.374 085 3 × 10−05 −6.071 144 4 × 10−05 2.393 359 5 × 10−04 
x3y −2.102 290 2 × 10−06 −1.361 552 3 × 10−06 5.367 495 7 × 10−06 
x3z 7.505 776 6 × 10−04 4.861 130 7 × 10−04 −1.916 349 3 × 10−03 
x2y2 −1.271 197 6 × 10−04 −8.232 935 5 × 10−05 3.245 578 3 × 10−04 
x2yz 1.801 066 1 × 10−03 1.166 463 9 × 10−03 −4.598 420 6 × 10−03 
x2z2 2.218 039 3 × 10−03 1.436 517 4 × 10−03 −5.663 022 3 × 10−03 
xy3 −8.741 210 4 × 10−06 −5.661 261 8 × 10−06 2.231 776 1 × 10−05 
xy2z 7.652 722 2 × 10−04 4.956 300 3 × 10−04 −1.953 866 9 × 10−03 
xyz2 −3.053 939 8 × 10−05 −1.977 890 0 × 10−05 7.797 215 1 × 10−05 
xz3 4.838 899 2 × 10−05 3.133 922 4 × 10−05 −1.235 451 2 × 10−04 
y4 −3.262 586 8 × 10−05 −2.113 020 7 × 10−05 8.329 925 5 × 10−05 
y3z 1.799 062 5 × 10−03 1.165 166 3 × 10−03 −4.593 305 0 × 10−03 
y2z2 2.276 793 8 × 10−03 1.474 569 9 × 10−03 −5.813 032 3 × 10−03 
yz3 3.033 272 8 × 10−03 1.964 505 0 × 10−03 −7.744 448 8 × 10−03 
z4 5.742 310 1 × 10−03 3.719 018 2 × 10−03 −1.466 107 0 × 10−02 
dB1/dt 9.486 847 6 × 10−01 9.444 878 4 × 10−07 −3.723 349 0 × 10−06 
dB2/dt −1.808 827 7 × 10−05 9.486 715 8 × 10−01 4.618 237 2 × 10−05 
dB2/dt 4.246 885 2 × 10−05 2.750 503 4 × 10−05 9.485 748 7 × 10−01 

We take cuboid D=[2,2]×[2,2]×[0,1], with boundary D, containing a sub-boundary Γ to be the surface z=1 (top boundary) or the surface z=0 (bottom boundary).

For the original system (26)(28), the mean residence time u (for solutions with initial points in D) and escape probability p (for solutions exiting through a sub-boundary Γ) satisfy the elliptic partial differential Eqs. (4)(5) and (6)(8), respectively. Hence,

12ϵΔu+(0.1x2y)ux+(2x0.1y)uy+(0.3z)uz=1,u|D=0.
12ϵΔp+(0.1x2y)px+(2x0.1y)py+(0.3z)pz=0,p|Γ=1,p|DΓ=0.

These equations can also be solved by a finite element method to get the mean residence time uTrue, the escape probability pTrue, and the average escape probability PTrue.

For the learned model (29)(31), we also have similar elliptic partial differential equations for the mean residence time uLearn, the escape probability pLearn and the average escape probability PLearn.

The mean residence time u(x,y,z) of the learned model is shown in Fig. 8. It is barely distinguishable from the mean residence time for the original system, and we thus do not show the latter. Instead, we compute the error (using 20 000 uniformly distributed points in D) between the maximal values of mean residence time for the learned and original systems: {error}(u)=max(x,y,z)DuLearnuTrue=3.7495×104.

FIG. 8.

Learned mean residence time uLearn for the learned linear stochastic system (29)(31), as viewed in two slices.

FIG. 8.

Learned mean residence time uLearn for the learned linear stochastic system (29)(31), as viewed in two slices.

Close modal

The escape probability p(x,y,z) of the learned system is shown in Fig. 9. It is barely distinguishable from the escape probability for the original system, and we thus do not show the latter. Instead, we calculate the error between the average escape probability values of the learned system and original system to be {error}(P)=PLearnPTrue=9.9478×107 (escaping from the top boundary) and {error}(P)=PLearnPTrue=1.0783×104 (escaping from the bottom boundary).

FIG. 9.

Learned escape probability for the learned linear stochastic system (29)–(31): escaping from the top surface (left) and the bottom surface (right). Two slices are shown for each subfigure.

FIG. 9.

Learned escape probability for the learned linear stochastic system (29)–(31): escaping from the top surface (left) and the bottom surface (right). Two slices are shown for each subfigure.

Close modal

We consider a stochastic Lorenz system,20 

dx=σ(yx)dt+ϵdB1,
(32)
dy=(ρxxzy)dt+ϵdB2,
(33)
dz=(xyβz)dt+ϵdB3,
(34)

where σ=10, β=83, ρ=28 are standard parameters and B1(t), B2(t), B3(t) are independent scalar Brownian motions. We take ϵ=0.9 in the following computations.

As before, a learned model is identified in the following form:

x˙=Θξ1,
(35)
y˙=Θξ2,
(36)
z˙=Θξ3,
(37)

where Θ is a set of basis function consisting of polynomials in (x,y,z) up to fourth order and time derivatives of three independent Brownian motions. Then, we solve a regression problem to determine the weights Ξ=[ξ1,ξ2,ξ3] (see Table IV).

TABLE IV.

Identified Lorenz system with basis functions. The sample paths data are numerically generated by solving (32)(34) with initial position (− 8, 7, 27)T, for t ∈ [0, 200] with stepsize 0.01.

Basisx˙y˙z˙
7.519 769 2 × 10−03 3.352 837 9 × 10−03 −1.369 622 8 × 10−02 
x −1.000 889 5 × 10+01 2.799 603 4 × 10+01 1.620 035 0 × 10−02 
y 1.000 469 0 × 10+01 −9.979 090 7 × 10−01 −8.541 396 7 × 10−03 
z −2.704 581 5 × 10−03 −1.205 891 2 × 10−03 −2.661 740 6 × 10+00 
x2 −4.563 336 2 × 10−04 −2.034 653 7 × 10−04 8.311 490 8 × 10−04 
xy −6.632 591 2 × 10−04 −2.957 272 1 × 10−04 1.001 208 0 × 10+00 
xz 1.013 984 2 × 10−03 −9.995 479 0 × 10−01 −1.846 833 1 × 10−03 
y2 4.180 401 0 × 10−04 1.863 914 5 × 10−04 −7.614 026 9 × 10−04 
yz −5.534 438 9 × 10−04 −2.467 639 1 × 10−04 1.008 022 1 × 10−03 
z2 3.086 039 3 × 10−04 1.375 971 7 × 10−04 −5.620 797 2 × 10−04 
x3 3.801 123 8 × 10−05 1.694 806 2 × 10−05 −6.923 225 5 × 10−05 
x2y −5.625 070 2 × 10−05 −2.508 048 9 × 10−05 1.024 529 4 × 10−04 
x2z 4.395 476 0 × 10−05 1.959 810 0 × 10−05 −8.005 756 6 × 10−05 
xy2 1.107 879 7 × 10−05 4.939 700 9 × 10−06 −2.017 850 9 × 10−05 
xyz 1.607 076 0 × 10−05 7.165 466 2 × 10−06 −2.927 068 4 × 10−05 
xz2 −3.438 403 0 × 10−05 −1.533 080 0 × 10−05 6.262 579 5 × 10−05 
y3 4.308 422 0 × 10−06 1.920 995 2 × 10−06 −7.847 199 7 × 10−06 
y2z −1.191 378 3 × 10−05 −5.311 996 1 × 10−06 2.169 932 2 × 10−05 
yz2 2.043 225 2 × 10−05 9.110 124 0 × 10−06 −3.721 454 4 × 10−05 
z3 −1.364 567 7 × 10−05 −6.084 195 1 × 10−06 2.485 372 9 × 10−05 
x4 2.471 056 3 × 10−07 1.101 769 3 × 10−07 −4.500 690 0 × 10−07 
x3y −1.531 861 5 × 10−06 −6.830 107 7 × 10−07 2.790 075 6 × 10−06 
x3z −6.061 635 5 × 10−07 −2.702 700 1 × 10−07 1.104 043 8 × 10−06 
x2y2 1.996 378 6 × 10−06 8.901 248 8 × 10−07 −3.636 129 7 × 10−06 
x2yz 1.158 815 1 × 10−06 5.166 806 5 × 10−07 −2.110 622 8 × 10−06 
x2z2 −1.034 200 2 × 10−06 −4.611 186 3 × 10−07 1.883 653 8 × 10−06 
xy3 −7.478 363 0 × 10−07 −3.334 376 1 × 10−07 1.362 081 2 × 10−06 
xy2z −2.273 687 1 × 10−07 −1.013 768 4 × 10−07 4.141 209 1 × 10−07 
xyz2 2.575 782 3 × 10−07 1.148 463 5 × 10−07 −4.691 434 1 × 10−07 
xz3 3.373 253 2 × 10−07 1.504 031 7 × 10−07 −6.143 917 9 × 10−07 
y4 −2.952 088 7 × 10−09 −1.316 247 2 × 10−09 5.376 824 5 × 10−09 
y3z −1.101 172 7 × 10−07 −4.909 796 3 × 10−08 2.005 635 0 × 10−07 
y2z2 −6.271 069 7 × 10−08 −2.796 080 5 × 10−08 1.142 189 3 × 10−07 
yz3 −2.384 655 6 × 10−07 −1.063 245 9 × 10−07 4.343 323 0 × 10−07 
z4 2.034 981 4 × 10−07 9.073 367 2 × 10−08 −3.706 439 4 × 10−07 
dB1/dt 9.458 149 5 × 10−01 −1.278 908 2 × 10−03 5.224 296 3 × 10−03 
dB2/dt −3.127 910 5 × 10−03 9.472 886 6 × 10−01 5.697 059 9 × 10−03 
dB3/dt 3.494 294 4 × 10−03 1.558 000 3 × 10−03 9.423 189 2 × 10−01 
Basisx˙y˙z˙
7.519 769 2 × 10−03 3.352 837 9 × 10−03 −1.369 622 8 × 10−02 
x −1.000 889 5 × 10+01 2.799 603 4 × 10+01 1.620 035 0 × 10−02 
y 1.000 469 0 × 10+01 −9.979 090 7 × 10−01 −8.541 396 7 × 10−03 
z −2.704 581 5 × 10−03 −1.205 891 2 × 10−03 −2.661 740 6 × 10+00 
x2 −4.563 336 2 × 10−04 −2.034 653 7 × 10−04 8.311 490 8 × 10−04 
xy −6.632 591 2 × 10−04 −2.957 272 1 × 10−04 1.001 208 0 × 10+00 
xz 1.013 984 2 × 10−03 −9.995 479 0 × 10−01 −1.846 833 1 × 10−03 
y2 4.180 401 0 × 10−04 1.863 914 5 × 10−04 −7.614 026 9 × 10−04 
yz −5.534 438 9 × 10−04 −2.467 639 1 × 10−04 1.008 022 1 × 10−03 
z2 3.086 039 3 × 10−04 1.375 971 7 × 10−04 −5.620 797 2 × 10−04 
x3 3.801 123 8 × 10−05 1.694 806 2 × 10−05 −6.923 225 5 × 10−05 
x2y −5.625 070 2 × 10−05 −2.508 048 9 × 10−05 1.024 529 4 × 10−04 
x2z 4.395 476 0 × 10−05 1.959 810 0 × 10−05 −8.005 756 6 × 10−05 
xy2 1.107 879 7 × 10−05 4.939 700 9 × 10−06 −2.017 850 9 × 10−05 
xyz 1.607 076 0 × 10−05 7.165 466 2 × 10−06 −2.927 068 4 × 10−05 
xz2 −3.438 403 0 × 10−05 −1.533 080 0 × 10−05 6.262 579 5 × 10−05 
y3 4.308 422 0 × 10−06 1.920 995 2 × 10−06 −7.847 199 7 × 10−06 
y2z −1.191 378 3 × 10−05 −5.311 996 1 × 10−06 2.169 932 2 × 10−05 
yz2 2.043 225 2 × 10−05 9.110 124 0 × 10−06 −3.721 454 4 × 10−05 
z3 −1.364 567 7 × 10−05 −6.084 195 1 × 10−06 2.485 372 9 × 10−05 
x4 2.471 056 3 × 10−07 1.101 769 3 × 10−07 −4.500 690 0 × 10−07 
x3y −1.531 861 5 × 10−06 −6.830 107 7 × 10−07 2.790 075 6 × 10−06 
x3z −6.061 635 5 × 10−07 −2.702 700 1 × 10−07 1.104 043 8 × 10−06 
x2y2 1.996 378 6 × 10−06 8.901 248 8 × 10−07 −3.636 129 7 × 10−06 
x2yz 1.158 815 1 × 10−06 5.166 806 5 × 10−07 −2.110 622 8 × 10−06 
x2z2 −1.034 200 2 × 10−06 −4.611 186 3 × 10−07 1.883 653 8 × 10−06 
xy3 −7.478 363 0 × 10−07 −3.334 376 1 × 10−07 1.362 081 2 × 10−06 
xy2z −2.273 687 1 × 10−07 −1.013 768 4 × 10−07 4.141 209 1 × 10−07 
xyz2 2.575 782 3 × 10−07 1.148 463 5 × 10−07 −4.691 434 1 × 10−07 
xz3 3.373 253 2 × 10−07 1.504 031 7 × 10−07 −6.143 917 9 × 10−07 
y4 −2.952 088 7 × 10−09 −1.316 247 2 × 10−09 5.376 824 5 × 10−09 
y3z −1.101 172 7 × 10−07 −4.909 796 3 × 10−08 2.005 635 0 × 10−07 
y2z2 −6.271 069 7 × 10−08 −2.796 080 5 × 10−08 1.142 189 3 × 10−07 
yz3 −2.384 655 6 × 10−07 −1.063 245 9 × 10−07 4.343 323 0 × 10−07 
z4 2.034 981 4 × 10−07 9.073 367 2 × 10−08 −3.706 439 4 × 10−07 
dB1/dt 9.458 149 5 × 10−01 −1.278 908 2 × 10−03 5.224 296 3 × 10−03 
dB2/dt −3.127 910 5 × 10−03 9.472 886 6 × 10−01 5.697 059 9 × 10−03 
dB3/dt 3.494 294 4 × 10−03 1.558 000 3 × 10−03 9.423 189 2 × 10−01 

A single trajectory of this stochastic system with initial condition (x,y,z)T=(8,7,27)T is shown in Fig. 10 (left), together with the same trajectory captured by the learned model in Fig. 10 (right).

FIG. 10.

One trajectory starting at (x,y,z)T=(8,7,27)T: for the original Lorenz system (left) and for the learned Lorenz system with basis in Table IV (right).

FIG. 10.

One trajectory starting at (x,y,z)T=(8,7,27)T: for the original Lorenz system (left) and for the learned Lorenz system with basis in Table IV (right).

Close modal

For the stochastic Lorenz system, we take a cuboid D to be either D1=[9,8]×[9,8]×[27,28] containing the left saddle as the left residence region or D2=[8,9]×[8,9]×[27,28] containing the right saddle as the right residence region. Then, mean residence time and escape probability satisfy the following elliptic partial differential equations, respectively:

12εΔu+(σ(yx))ux+(ρxxzy)uy+(xyβz)uz=1,u|Di=0,i=1,2.
12εΔp+(σ(yx))px+(ρxxzy)py+(xyβz)pz=0,p|Γ=1,p|DiΓ=0,i=1,2.

The finite element numerical solutions to these partial differential equations are denoted by uTrue and pTrue.

For the learned Lorenz system (35)(37), with basis and coefficients in Table IV, we can also set up the partial differential equations for the learned mean residence time uLearn, and the learned escape probability pLearn together with the average escape probability PLearn.

The learned mean residence time u of the learned Lorenz system (35)(37) is shown in Fig. 11. The maximal mean residence time in the left region D1 is uLearn(D1)=0.1238, while the maximal mean residence time in the right region D2 is uLearn(D2)=0.1236. Then, we compute the error of maximal mean residence time between the learned and the original systems, {error}(u(D1))=max(x,y,z)D1uLearnuTrue=1.0023×104 and {error}(u(D2))=max(x,y,z)D2uLearnuTrue=1.2284×105.

FIG. 11.

Learned mean residence time of the learned Lorenz system (35)(37) in D1 (left) and D2 (right): two slices are shown for each subfigure.

FIG. 11.

Learned mean residence time of the learned Lorenz system (35)(37) in D1 (left) and D2 (right): two slices are shown for each subfigure.

Close modal

We take Γ to be each surface of the region. The learned escape probability pLearn of particles exiting through each surface Γ of left and right regions are shown in Figs. 12 and 13, respectively, and the average escape probability of left region is PLearn(L1)=0.1065 (from x=9), PLearn(L2)=0.1167 (from x=8), PLearn(L3)=0.0955 (from y=9), PLearn(L4)=0.3485 (from y=8), PLearn(L5)=0.2343 (from z=27), PLearn(L6)=0.1244 (from z=28). Similarly, the average escape probability of the right region is PLearn(R1)=0.1161 (from x=8), PLearn(R2)=0.1061 (from x=9), PLearn(R3)=0.3469 (from y=8), PLearn(R4)=0.0947 (from y=9), PLearn(R5)=0.2367 (from z=27), PLearn(R6)=0.1264 (from z=28).

FIG. 12.

Learned escape probability of the learned Lorenz system (35)(37) in left region D1 exiting from a sub-boundary Γ: surface x=9 (upper left), x=8 (upper middle), y=9 (upper right), y=8 (lower left), z=27 (lower middle), and z=28 (lower right). Two slices are shown for each subfigure.

FIG. 12.

Learned escape probability of the learned Lorenz system (35)(37) in left region D1 exiting from a sub-boundary Γ: surface x=9 (upper left), x=8 (upper middle), y=9 (upper right), y=8 (lower left), z=27 (lower middle), and z=28 (lower right). Two slices are shown for each subfigure.

Close modal
FIG. 13.

Learned escape probability of the learned Lorenz system (35)(37) in right region D2 exiting from a sub-boundary Γ: surface x=8 (upper left), x=9 (upper middle), y=8 (upper right), y=9 (lower left), z=27 (lower middle), and z=28 (lower right). Two slices are shown for each subfigure.

FIG. 13.

Learned escape probability of the learned Lorenz system (35)(37) in right region D2 exiting from a sub-boundary Γ: surface x=8 (upper left), x=9 (upper middle), y=8 (upper right), y=9 (lower left), z=27 (lower middle), and z=28 (lower right). Two slices are shown for each subfigure.

Close modal

Moreover, for example, we compute the error of the average escape probability between the learned and the original systems only for two cases: {error}(P)(L4)=PLearn(L4)PTrue(L4)=6.5838×103 and {error}(P)(R3)=PLearn(R3)PTrue(R3)=6.431×103 (L4).

In summary, we have demonstrated a new data-driven method to determine dynamical quantities, mean residence time and escape probability, based on sample path data of complex systems, as long as these systems are modeled by stochastic differential equations.2,12,14,19,22,23

A novelty of this method is that it contains noise terms (in terms of Brownian motions) in the basis, in order to discover the governing stochastic differential equation. With the governing stochastic differential equation as the model for the sample path data, we can then compute the mean residence time and escape probability, which are deterministic quantities that carry significant dynamical information.

This method combines machine learning tools (i.e., extracting governing stochastic differential equations from data) and stochastic dynamical systems techniques (i.e., understanding random phenomena with deterministic quantities) to provide an efficient approach in examining stochastic dynamics from data.

To demonstrate that our method is effective, we have illustrated it on several stochastic systems with Brownian motions. It is expected to be applicable to data sets from other stochastic dynamical systems.

We would like to thank Min Dai and Jian Ren for helpful comments. The first author was supported by the China Scholarship Council under No. 201808220008, and the second author was supported by the China Scholarship Council under No. 201808220009. The second author was also supported by the Jilin Province Industrial Technology Research and Development Special Project under No. 2016C079. The third author was partly supported by the National Science Foundation (NSF) Grant No. 1620449.

1.
L.
Arnold
,
Random Dynamical Systems
(
Springer
,
New York
,
1998
).
2.
L.
Arnold
,
Stochastic Differential Equations: Theory and Applications
(
Wiley-Interscience
,
New York
,
1974
), pp. xvi+228.
3.
D.
Biegie
,
A.
Lenoad
, and
S.
Wiggins
, “
The dynamics associated with the chaotic tangles of two-dimensional quasiperiodic vector fields: Theory and applications
,”
Nonlinear Phenom. Atmos. Ocean. Sci.
40
,
47
138
(
1992
).
4.
A. S.
Bower
and
T.
Rossby
, “
Evidence of cross-frontal exchange processes in the Gulf stream based on isopycnal RAFOS float data
,”
J. Phys. Oceanogr.
19
,
1177
1190
(
1989
).
5.
A. S.
Bower
, “
A simple kinematic mechanism for mixing fluid parcels across a meandering jet
,”
J. Phys. Oceanogr.
21
,
173
180
(
1991
).
6.
J.
Duan
and
S.
Wiggins
, “
Fluid exchange across a meandering jet with quasiperiodic variability
,”
J. Phys. Oceanogr.
26
,
1176
1188
(
1995
).
7.
J. R.
Brannan
,
J.
Duan
, and
V. J.
Ervin
, “
Escape probability, mean residence time and geophysical fluid particle dynamics
,”
Physica D
133
,
23
33
(
1999
).
8.
D.
del-Castillo-Negrete
and
P. J.
Morris
, “
Chaotic transport by Rossby waves in shear flow
,”
Physica Fluid A
5
,
948
965
(
1993
).
9.
J.
Duan
,
An Introduction to Stochastic Dynamics
(
Cambridge University Press
,
New York
,
2015
).
10.
T.
Gao
and
J.
Duan
, “
Quantifying model uncertainty in dynamical systems driven by non-Gaussian Lévy stable noise with observations on mean exit time or escape probability
,”
Commun. Nonlinear Sci. Numer. Simulat.
39
,
1
6
(
2016
).
11.
L. J.
Pratt
,
M. S.
Lozier
, and
N.
Beliakova
, “
Parcel trajetories in quasigeostrophic jets: neutral models
,”
J. Phys. Oceanogr.
25
,
1451
1466
(
1995
).
12.
P.
Miller
,
C. K. R. T.
Jones
,
G.
Haller
, and
L.
Pratt
, “
Chaotic mixing across oceanic jets
,”
AIP Conf. Proc.
275
,
591
604
(
1996
).
13.
J.
Pedlosky
,
Geophysical Fluid Dynamics
, 2nd ed. (
Springer
,
New York
,
1987
).
14.
L.
Perko
,
Differential Equations and Dynamical Systems
, 3rd ed. (
Springer
,
New York
,
2001
).
15.
R. M.
Samelson
, “
Fluid exchange across a meandering jet
,”
J. Phys. Oceanogr.
22
(
4
),
431
440
(
1992
).
16.
R. M.
Samelson
, “
Stochastic forced current fluctuations in vertical shear and over topography
,”
J. Geographical Res.
94
,
8207
8215
(
1989
).
17.
B.
Oksendal
,
Stochastic Differential Equations
, 6th ed. (
Springer
,
New York
,
2005
).
18.
S.
Zhang
and
G.
Lin
, “
Robust data-driven discovery of governing physical laws with error bars
,”
Proc. R. Soc. A Math. Phys. Eng. Sci.
474
,
1
21
(
2018
).
19.
D. J.
Higham
, “
An algorithmic introduction to numerical simulation of stochastic differential equations
,”
SIAM Rev.
43
(
3
),
525
546
(
2001
).
20.
S. L.
Brunton
,
J. L.
Proctor
, and
J. N.
Kutz
, “
Discovering governing equations from data by sparse identification of nonlinear dynamical systems
,”
Proc. Natl. Acad. Sci. U.S.A.
113
(
15
),
3932
3937
(
2016
).
21.
L.
Dunker
,
G.
Bohner
,
J.
Boussard
, and
M.
Sahani
, “
Learning interpretable continuous-time models of latent stochastic dynamical systems
,”
Statistics
113
(
15
),
3932
3937
(
2019
).
22.
C.
Cotter
,
D.
Crisan
,
D. D.
Holm
,
W.
Pan
, and
I.
Shevchenko
, “Modeling uncertainty using circulation-preserving stochastic transport noise in a 2-layer quasi-geostrohic model,” e-print arXiv:1802.05711 (2018).
23.
M.
Crosskey
and
M.
Maggioni
, “
ATLAS: A geometric approach to learning high-dimensional stochastic systems near manifolds
,”
Multiscale Model. Simul.
15
(
1
),
110
156
(
2017
).