Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

HTML conversions sometimes display errors due to content that did not convert correctly from the source. This paper uses the following packages that are not yet supported by the HTML conversion tool. Feedback on these issues are not necessary; they are known and are being worked on.

  • failed: extarrows

Authors: achieve the best HTML results from your LaTeX submissions by following these best practices.

License: arXiv.org perpetual non-exclusive license
arXiv:2402.18315v1 [math.DS] 28 Feb 2024

Rare events in a stochastic vegetation-water dynamical system based on machine learning

Yang Li11footnotemark: 1 liyangbx5433@163.com School of Automation, Nanjing University of Science and Technology, 200 Xiaolingwei Street, Nanjing 210094, Jiangsu Province, China Shenglan Yuan22footnotemark: 233footnotemark: 3 shenglanyuan@gbu.edu.cn School of Sciences, Great Bay University, Songshan Lake International Innovation Entrepreneurship Community A5, Dongguan 523000, Guangdong Province, China Great Bay Institute for Advanced Study, Songshan Lake International Innovation Entrepreneurship Community A5, Dongguan 523000, Guangdong Province, China Shengyuan Xu44footnotemark: 4 syxu@njust.edu.cn
Abstract

Stochastic vegetation-water dynamical systems play a pivotal role in ecological stability, biodiversity, water resource management, and adaptation to climate change. This research proposes a machine learning-based method for analyzing rare events in stochastic vegetation-water dynamical systems with multiplicative Gaussian noise. Utilizing the Freidlin-Wentzell large deviation theory, we derive the asymptotic expressions for the quasipotential and the mean first exit time. Based on the decomposition of vector field, we design a neural network architecture to compute the most probable transition paths and the mean first exit time for both non-characteristic and characteristic boundary scenarios. The results indicate that this method can effectively predict early warnings of vegetation degradation, providing new theoretical foundations and mathematical tools for ecological management and conservation. Moreover, the method offers new possibilities for exploring more complex and higher-dimensional stochastic dynamical systems.

keywords:
Stochastic vegetation-water system, Rare events, Machine learning, Most probable path, Mean first exit time
journal: Applied Mathematics and Computation

1 Introduction

The stochastic noise is indeed very common in ecosystems, which refers to random fluctuations in ecological systems that are not predictable or explainable by deterministic factors QY ; YLZ ; TTYBD . These fluctuations can arise from a variety of sources, including environmental conditions, dispersal patterns and interactions among species, etc. Under random fluctuations, rare events SSCA , i.e., transition from one stable state to another, frequently occur, even for weak noise. Thus they are well worth investigating, such as vegetation degradation and species extinction YW .

The Freidlin-Wentzell large deviation theory FW is actually a principle that is mainly applied to study rare events of dynamical systems with small random perturbations. As the noise intensity decreases, the convergence rate at which the sample trajectories converge to the reference orbit is exponential in terms of noise. Large deviation techniques are widely used, which have become an extremely active branch of applied probability. They can estimate the escape probability of stochastic systems, reckon the probability of deviation from the reference orbit, and quantify the asymptotic probability of errors in hypothesis testing.

In large deviation theory, the action functional is an important concept that builds the relationship between the stationary probability distribution and the path distribution of stochastic dynamical systems FW . It is defined as the integral of a certain Lagrangian, similar to classical mechanics. The minimization of the action functional leads to the most probable path connecting given initial and final states, along which the stochastic system moves with highest probability than other paths YD ; LDLZ . Although variables in statistical mechanics do not move along a definite trajectory, the most probable path offers us a very intuitive way to comprehend the stochastic behaviors of the dynamical systems and predict their evolution. Furthermore, the most probable path can also be used in optimization and control problems to find optimal solutions that maximize or minimize the objective functions in stochastic dynamical systems, which provides an appropriate and effective method for exploring the properties of stochastic models in practical applications LKBMM . The calculation of the most probable path typically involves statistical methods and numerical simulation techniques. For instance, the Monte Carlo method is a random sampling technique, which can be utilized to simulate the behavior of a system under diverse conditions and parameters. Through extensive simulation and statistics, we can more accurately find the most probable path and its associated probability DMSSS .

In addition, quasipotential is commonly used in physics and engineering to generalize the concept of potential function to nongradient dynamical systems FW ; ZAAH ; Ao . It is defined as the global minimal value of action functional about both possible paths and time length, which can be used to describe rare transition events in various physical systems, such as fluid mechanics, electrodynamics and ecosystems. The quasipotential function exponentially dominates the magnitude of mean first exit time and stationary probability density.

Mean first exit time is a statistical quantity that represents the expected value of the time required for a stochastic system to escape from a confined state or region for the first time. In stochastic processes and diffusion theory, it provides a measure of how long a stochastic system can transition from one state to another, which can be valuable in grasping and optimizing the performance of complex dynamical system related to the properties of the medium, such as its diffusivity, geometry, and boundary conditions. In practical applications, it is a critical component used to assess the engineering reliability of first passage failure CZ and to describe the activation process in neuron systems FTPVB . Apart from the quasipotential, more accurate perturbation expression of mean first exit time also depends on the exponential prefactor function of WKB approximation NKMS ; MST ; MS ; LZXDL .

The traditional numerical methods for calculating these previously mentioned quantities include the geometric minimum action method (GMAM) HVE and ordered upwind method (OUM) C . The former method is grounded in geometric action, aiming to iterate the most probable paths with fixed endpoints in the path space and obtain their corresponding action functionals. It converts the time parameterization of the action functional into arc length, simplifying the integration over infinite time to a finite-length integration, thereby finding the path corresponding to the minimum action. OUM is a numerical method for discretizing phase space by considering the influence of each node in an ordered manner and using an upwind difference scheme to compute the quasipotential at each node. This method facilitates efficient computation within the discretized phase space. GMAM and OUM are two different numerical methods based on the concepts of geometric action and discretization, respectively, for solving different types of problems. However, both GMAM and OUM have nonnegligible shortcomings, such as multiple minima for GMAM and great computational cost for OUM, especially for the high-dimensional case. Due to their limitations, there is a need to explore new methods.

Machine learning is an important direction in the field of artificial intelligence A . It analyzes large amounts of data and improves existing algorithms to make them more intelligent and enhance their generalization capabilities SWS . The powerful features of machine learning are mainly manifested in the following aspects. Firstly, machine learning is data-driven. By learning from a large amount of data, it can discover patterns and rules in the data, enabling more accurate predictions and decisions B . Secondly, machine learning algorithms can automatically adjust parameters and models, reducing human intervention and errors; thus, making the results more objective and precise HLWFKB . Thirdly, machine learning algorithms can quickly adapt to different datasets and tasks, exhibiting strong adaptability, facilitating learning and application across diverse environments AR . Fourthly, as algorithms continuously improve and optimize, the interpretability of machine learning models is also steadily improving, enabling a better understanding of the internal mechanisms of data and models CCC . Fifthly, machine learning has widespread applications in various fields, such as image recognition, speech recognition, natural language processing, recommendation systems, financial risk control, etc., bringing about tremendous changes and adding significant business value across different sectors BB . In summary, machine learning is a powerful technology that facilitates more accurate and effective data processing through automatic learning and optimization algorithms JM . With the continuous development of technology, the application prospects of machine learning will become even broader.

At present, machine learning has been widely applied in the study of stochastic dynamics. For instance, researchers designed some data-driven methods to discover stochastic dynamical systems with Gaussian or non-Gaussian noise LD ; CYDK . Xu et al. developed a novel deep learning method to compute the probability density by solving the Fokker-Planck equation XZLZLK . There are also some scholars devoted to exploring data-driven machine learning methods for computing the most probable paths of stochastic dynamical systems LDL ; WGCD . Based on large deviation theory and machine learning, Li et al. LYX calculated most probable transition path and designed a control strategy to control the mean first exit time of stochastic dynamical systems to achieve a desired value. Machine learning can also be used to compute the large deviation prefactors in the highly complex nonlinear stochastic systemsLYLL . These methods to study rare events are mainly used for stochastic systems with additive noise.

The investigation of stochastic vegetation-water ecosystems is of great significance due to their notable impact on ecological stability, biodiversity, water resource management, adaptation to climate change, and soil health. Considering their complex structures with multiplicative noise ZXLQ , this paper aims to develop a machine learning method to handle rare transition events of stochastic dynamical systems with multiplicative Gaussian noise, rather than the case with additive noise LZXDL ; LYLL ; LLR . Furthermore, this method will be utilized to calculate the most probable path, quasipotential, and mean first exit time of the stochastic vegetation-water system ZXLQ . The analysis of the dynamic behavior of this system provides a theoretical basis and mathematical methods for understanding and controlling the phenomenon of vegetation degradation.

The structure of this work is as follows. In Section 2, we present the framework of the stochastic vegetation-water system and explain the meaning of the involved parameters. We reveal the dynamical structures of the corresponding deterministic system and simulate the transition phenomenon between the metastable states of stochastic system. In Section 3, we describe the concepts of large deviation theory, quasipotential and mean first exit time, and derive their expressions asymptotically. Then a machine learning method is proposed to compute these quantities for the vegetation-water system in both non-characteristic and characteristic boundary cases. In Section 4, we use the machine learning algorithm to compute the rare transition events of the stochastic vegetation-water system and provide a mathematical basis for early warning of vegetation degradation. In Section 5, we summarize the results of the paper and discuss some important future prospects.

2 Stochastic vegetation-water system

The vegetation-water dynamical system is a complex ecosystem where the vegetation and water interact to maintain the balance. However, due to various influencing factors such as climate change, terrain, soil type, and the hydrological cycle, the dynamics of the vegetation-water system often exhibit uncertain and unpredictable stochasticity that manifests in several ways. Firstly, meteorological processes, such as rainfall intensity and evaporation, are stochastic, leading to stochasticity in hydrological processes (e.g., river discharge, groundwater levels). These hydrological processes directly influence the state of the vegetation. Secondly, the growth and distribution of vegetation are also affected by many stochastic factors. For example, seed dispersal, plant growth rates, and the occurrence of pests and diseases are stochastic processes that lead to stochasticity in vegetation spatial distribution and density. In addition, the interaction between vegetation and water in the system is also stochastic. For instance, vegetation consumes water through transpiration, and water returns to the system through soil infiltration, surface runoff, etc. To understand and predict the stochasticity of the vegetation-water dynamical system, mathematical tools such as probability theory and stochastic processes are required to establish quantitative models of the system’s state and behavior. Furthermore, modern technological means, like machine leaning techniques, are employed to obtain extensive observational data of the system, supporting model validation and refinement. The stochasticity of vegetation-water systems is a crucial characteristic that plays a key role in predicting and responding to practical issues such as water resource management and ecological restoration. Studying and comprehending this stochasticity is essential for effective management and conservation efforts.

More specifically, we consider a stochastic vegetation-water system

x˙1=ρx1(x2x1K)βx1x1+x0+σ1x12ξ1(t),x˙2=Rαx2λx1x2+σ2ξ2(t),subscript˙𝑥1𝜌subscript𝑥1subscript𝑥2subscript𝑥1𝐾𝛽subscript𝑥1subscript𝑥1subscript𝑥0subscript𝜎1superscriptsubscript𝑥12subscript𝜉1𝑡subscript˙𝑥2𝑅𝛼subscript𝑥2𝜆subscript𝑥1subscript𝑥2subscript𝜎2subscript𝜉2𝑡\begin{array}[]{l}\dot{x}_{1}=\rho x_{1}\big{(}x_{2}-\frac{x_{1}}{K}\big{)}-% \beta\frac{x_{1}}{x_{1}+x_{0}}+\sigma_{1}x_{1}^{2}\xi_{1}(t),\\ \dot{x}_{2}=R-\alpha x_{2}-\lambda x_{1}x_{2}+\sigma_{2}\xi_{2}(t),\end{array}start_ARRAY start_ROW start_CELL over˙ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_ρ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - divide start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_K end_ARG ) - italic_β divide start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG + italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ξ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_t ) , end_CELL end_ROW start_ROW start_CELL over˙ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_R - italic_α italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_λ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_ξ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( italic_t ) , end_CELL end_ROW end_ARRAY (2.1)

where the variable x1subscript𝑥1x_{1}italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT represents the biomass of vegetation, while x2subscript𝑥2x_{2}italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT signifies the moisture level of the soil. The interaction terms ρx1x2𝜌subscript𝑥1subscript𝑥2\rho x_{1}x_{2}italic_ρ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT and λx1x2𝜆subscript𝑥1subscript𝑥2-\lambda x_{1}x_{2}- italic_λ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT depict the relationship between vegetation biomass and water. The term ρx12/K𝜌superscriptsubscript𝑥12𝐾-\rho x_{1}^{2}/K- italic_ρ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_K limits the growth of biomass due to competition for shared resources, such as water or soil nutrients. The term βx1/(x1+x0)𝛽subscript𝑥1subscript𝑥1subscript𝑥0-\beta x_{1}/(x_{1}+x_{0})- italic_β italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) illustrates the impact of herbivores and other influencing factors. The parameter R𝑅Ritalic_R denotes the average rainfall, and the term αx2𝛼subscript𝑥2-\alpha x_{2}- italic_α italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT stands for the loss of water from the soil, which could be caused by percolation or evaporation. Taking into account random environmental disturbances affecting vegetation competition for shared resources and rainfall, the term ρx12/K𝜌superscriptsubscript𝑥12𝐾-\rho x_{1}^{2}/K- italic_ρ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_K transforms into ρx12(1+σ~1ξ1(t))/K𝜌superscriptsubscript𝑥121subscript~𝜎1subscript𝜉1𝑡𝐾-\rho x_{1}^{2}(1+\tilde{\sigma}_{1}\xi_{1}(t))/K- italic_ρ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( 1 + over~ start_ARG italic_σ end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_ξ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_t ) ) / italic_K and R𝑅Ritalic_R becomes R(1+σ~2ξ2(t))𝑅1subscript~𝜎2subscript𝜉2𝑡R(1+\tilde{\sigma}_{2}\xi_{2}(t))italic_R ( 1 + over~ start_ARG italic_σ end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_ξ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( italic_t ) ). Consequently, these disturbances are captured by multiplicative noise σ1x12ξ1(t)subscript𝜎1superscriptsubscript𝑥12subscript𝜉1𝑡\sigma_{1}x_{1}^{2}\xi_{1}(t)italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ξ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_t ) and additive noise σ2ξ2(t)subscript𝜎2subscript𝜉2𝑡\sigma_{2}\xi_{2}(t)italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_ξ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( italic_t ), where σ1=ρσ~1/Ksubscript𝜎1𝜌subscript~𝜎1𝐾\sigma_{1}=-\rho\tilde{\sigma}_{1}/Kitalic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = - italic_ρ over~ start_ARG italic_σ end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_K and σ2=Rσ~2subscript𝜎2𝑅subscript~𝜎2\sigma_{2}=R\tilde{\sigma}_{2}italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_R over~ start_ARG italic_σ end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT. The driving terms ξ1(t)subscript𝜉1𝑡\xi_{1}(t)italic_ξ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_t ) and ξ2(t)subscript𝜉2𝑡\xi_{2}(t)italic_ξ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( italic_t ) are independent Gaussian white noises with

𝔼[ξi(t)]=0,𝔼[ξi(t)ξi(t+τ)]=εδ(τ),i=1,2.formulae-sequence𝔼delimited-[]subscript𝜉𝑖𝑡0formulae-sequence𝔼delimited-[]subscript𝜉𝑖𝑡subscript𝜉𝑖𝑡𝜏𝜀𝛿𝜏𝑖12\mathbb{E}[\xi_{i}(t)]=0,\quad\mathbb{E}[\xi_{i}(t)\xi_{i}(t+\tau)]=% \varepsilon\delta(\tau),\quad i=1,2.blackboard_E [ italic_ξ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_t ) ] = 0 , blackboard_E [ italic_ξ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_t ) italic_ξ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_t + italic_τ ) ] = italic_ε italic_δ ( italic_τ ) , italic_i = 1 , 2 .

Here, ε𝜀\varepsilonitalic_ε is a small parameter, implying that the noise intensity is weak. The drift coefficient and diffusion matrix can be written as

b(x)=(ρx1(x2x1K)βx1x1+x0Rαx2λx1x2),a(x)=σ(x)σ(x)=(σ12x1400σ22),formulae-sequence𝑏𝑥𝜌subscript𝑥1subscript𝑥2subscript𝑥1𝐾𝛽subscript𝑥1subscript𝑥1subscript𝑥0𝑅𝛼subscript𝑥2𝜆subscript𝑥1subscript𝑥2𝑎𝑥𝜎𝑥superscript𝜎top𝑥superscriptsubscript𝜎12superscriptsubscript𝑥1400superscriptsubscript𝜎22b(x)=\left(\begin{array}[]{c}\rho x_{1}\big{(}x_{2}-\frac{x_{1}}{K}\big{)}-% \beta\frac{x_{1}}{x_{1}+x_{0}}\\ R-\alpha x_{2}-\lambda x_{1}x_{2}\\ \end{array}\right),\quad a(x)=\sigma(x)\sigma^{\top}(x)=\left(\begin{array}[]{% cc}\sigma_{1}^{2}x_{1}^{4}&0\\ 0&\sigma_{2}^{2}\\ \end{array}\right),italic_b ( italic_x ) = ( start_ARRAY start_ROW start_CELL italic_ρ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - divide start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_K end_ARG ) - italic_β divide start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG end_CELL end_ROW start_ROW start_CELL italic_R - italic_α italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_λ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_CELL end_ROW end_ARRAY ) , italic_a ( italic_x ) = italic_σ ( italic_x ) italic_σ start_POSTSUPERSCRIPT ⊤ end_POSTSUPERSCRIPT ( italic_x ) = ( start_ARRAY start_ROW start_CELL italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_CELL end_ROW end_ARRAY ) ,

where

σ(x)=(σ1x1200σ2).𝜎𝑥subscript𝜎1superscriptsubscript𝑥1200subscript𝜎2\sigma(x)=\left(\begin{array}[]{cc}\sigma_{1}x_{1}^{2}&0\\ 0&\sigma_{2}\\ \end{array}\right).italic_σ ( italic_x ) = ( start_ARRAY start_ROW start_CELL italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_CELL end_ROW end_ARRAY ) .

We first consider the corresponding deterministic system with σ1=σ2=0subscript𝜎1subscript𝜎20\sigma_{1}=\sigma_{2}=0italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 0:

x˙1=ρx1(x2x1K)βx1x1+x0,x˙2=Rαx2λx1x2.subscript˙𝑥1𝜌subscript𝑥1subscript𝑥2subscript𝑥1𝐾𝛽subscript𝑥1subscript𝑥1subscript𝑥0subscript˙𝑥2𝑅𝛼subscript𝑥2𝜆subscript𝑥1subscript𝑥2\begin{array}[]{l}\dot{x}_{1}=\rho x_{1}\big{(}x_{2}-\frac{x_{1}}{K}\big{)}-% \beta\frac{x_{1}}{x_{1}+x_{0}},\\ \dot{x}_{2}=R-\alpha x_{2}-\lambda x_{1}x_{2}.\end{array}start_ARRAY start_ROW start_CELL over˙ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_ρ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - divide start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_K end_ARG ) - italic_β divide start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG , end_CELL end_ROW start_ROW start_CELL over˙ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_R - italic_α italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_λ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT . end_CELL end_ROW end_ARRAY (2.2)

In this paper, we fix these system parameters as ρ=1𝜌1\rho=1italic_ρ = 1, K=10𝐾10K=10italic_K = 10, β=3𝛽3\beta=3italic_β = 3, x0=1subscript𝑥01x_{0}=1italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 1, α=1𝛼1\alpha=1italic_α = 1, λ=0.12𝜆0.12\lambda=0.12italic_λ = 0.12.

Let b(x)=𝟎𝑏𝑥𝟎b(x)=\textbf{0}italic_b ( italic_x ) = 0 and we can derive the fixed points of the system (2.2). It is observed that SN1=(0,Rα)SN10𝑅𝛼\text{SN1}=\big{(}0,\frac{R}{\alpha}\big{)}SN1 = ( 0 , divide start_ARG italic_R end_ARG start_ARG italic_α end_ARG ) is a trivial fixed point for arbitrary R𝑅Ritalic_R. When x10subscript𝑥10x_{1}\neq 0italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ≠ 0, b(x)=𝟎𝑏𝑥𝟎b(x)=\textbf{0}italic_b ( italic_x ) = 0 implies that

f(x1)=ρλx13+ρ(λx0+α)x12+(ραx0+KβλρRK)x1+KβαρRKx0=0,x2=Rλx1+α.𝑓subscript𝑥1𝜌𝜆superscriptsubscript𝑥13𝜌𝜆subscript𝑥0𝛼superscriptsubscript𝑥12𝜌𝛼subscript𝑥0𝐾𝛽𝜆𝜌𝑅𝐾subscript𝑥1𝐾𝛽𝛼𝜌𝑅𝐾subscript𝑥00subscript𝑥2𝑅𝜆subscript𝑥1𝛼\begin{array}[]{l}f(x_{1})=\rho\lambda x_{1}^{3}+\rho(\lambda x_{0}+\alpha)x_{% 1}^{2}+(\rho\alpha x_{0}+K\beta\lambda-\rho RK)x_{1}+K\beta\alpha-\rho RKx_{0}% =0,\\ x_{2}=\frac{R}{\lambda x_{1}+\alpha}.\end{array}start_ARRAY start_ROW start_CELL italic_f ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) = italic_ρ italic_λ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + italic_ρ ( italic_λ italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_α ) italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( italic_ρ italic_α italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_K italic_β italic_λ - italic_ρ italic_R italic_K ) italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_K italic_β italic_α - italic_ρ italic_R italic_K italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0 , end_CELL end_ROW start_ROW start_CELL italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = divide start_ARG italic_R end_ARG start_ARG italic_λ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_α end_ARG . end_CELL end_ROW end_ARRAY

As shown in Fig. 1, a saddle-node bifurcation occurs at R=Rc𝑅subscript𝑅𝑐R=R_{c}italic_R = italic_R start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT. Therefore, a node SN2 and a saddle US emerge. Let f(x1)=0𝑓subscript𝑥10f(x_{1})=0italic_f ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) = 0 and f(x1)=0superscript𝑓subscript𝑥10f^{\prime}(x_{1})=0italic_f start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) = 0. We have Rc=1.4278subscript𝑅𝑐1.4278R_{c}=1.4278italic_R start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 1.4278. SN1 is stable if R<2.998𝑅2.998R<2.998italic_R < 2.998, otherwise it is unstable due to the collision of US. Above all, there are three states for deterministic vegetation-water system (2.2) depending on the value range of R𝑅Ritalic_R: bare for R<Rc𝑅subscript𝑅𝑐R<R_{c}italic_R < italic_R start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT, bistable for Rc<R<2.998subscript𝑅𝑐𝑅2.998R_{c}<R<2.998italic_R start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT < italic_R < 2.998, and vegetated for R>2.998𝑅2.998R>2.998italic_R > 2.998.

Refer to caption
Figure 1: Bifurcation digram of the vegetation system about the parameter R𝑅Ritalic_R. It exhibits two branches, the blue curve representing the stable equilibria, while the red dashed branch representing the unstable equilibria. Upon varying the control parameter R𝑅Ritalic_R, the two branches approach each other until they meet at the critical value of Rc=1.4278subscript𝑅𝑐1.4278R_{c}=1.4278italic_R start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 1.4278. At this point, they annihilate each other in a saddle-node bifurcation.

It can be seen from the saddle-node bifurcation diagram in Fig. 1 that the biomass of vegetation depends greatly on the average rainfall R𝑅Ritalic_R. This finding aligns with the actual ecological phenomena. When R<Rc𝑅subscript𝑅𝑐R<R_{c}italic_R < italic_R start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT, the rainfall is insufficient for vegetation to survive, which corresponds to the bare state; when R>2.998𝑅2.998R>2.998italic_R > 2.998, the rainfall is ample, and thus the vegetation can grow fully without the phenomenon of vegetation disappearance, corresponding to the vegetated state. Meanwhile, for rainfall values within the range Rc<R<2.998subscript𝑅𝑐𝑅2.998R_{c}<R<2.998italic_R start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT < italic_R < 2.998, the bistable phenomenon emerges as the rainfall is neither too little nor too much. Investigating the vegetation state within this range of rainfall has practical significance for applications.

In this paper, we choose R=1.55𝑅1.55R=1.55italic_R = 1.55, i.e., bistable state, for investigation. As demonstrated in Fig. 2, the system (2.2) has two stable fixed points SN1=(0,1.55)SN101.55\text{SN1}=(0,1.55)SN1 = ( 0 , 1.55 ) and SN2=(4.6366,0.9959)SN24.63660.9959\text{SN2}=(4.6366,0.9959)SN2 = ( 4.6366 , 0.9959 ). The basins of attraction for these fixed points are separated by the stable manifold of the saddle US=(1.6667,1.2917)US1.66671.2917\text{US}=(1.6667,1.2917)US = ( 1.6667 , 1.2917 ), denoted by a purple curve. The unstable manifold of US is indicated by a green curve.

Refer to caption
Figure 2: For the selected value of R=1.55𝑅1.55R=1.55italic_R = 1.55, the system (2.2) possesses two stable fixed points SN1=(0,1.55)𝑆𝑁101.55SN1=(0,1.55)italic_S italic_N 1 = ( 0 , 1.55 ) and SN2=(4.6366,0.9959)𝑆𝑁24.63660.9959SN2=(4.6366,0.9959)italic_S italic_N 2 = ( 4.6366 , 0.9959 ), which are separated by the stable manifold of the saddle point US=(1.6667,1.2917). This stable manifold is delineated by a purple curve. Additionally, the unstable manifold of the saddle point US is denoted by a green curve.

Based on the ecological significance of stochastic vegetation-water model (2.1), we investigate the transition phenomena of the system initially located at SN2, i.e., with lush vegetation, under random perturbations approaching SN1. Since the noise intensity ε𝜀\varepsilonitalic_ε is small, these phenomena are referred to as rare events. Indeed, when the system crosses the boundary, i.e., the stable manifold of US, it will flow along the unstable manifold of the US to the point SN1. Therefore, we mainly focus on the process of the system escaping from the attractor domain of SN2 driven by noise perturbations, which can provide important information for early warning and intervention of vegetation loss. A typical noise induced transition trajectory simulated by Monte Carlo is illustrated in Fig. 3.

Refer to caption
Figure 3: A representative transition trajectory of the stochastic vegetation-water model (2.1) depicted through Monte Carlo simulation.

3 Theory and method

Under the condition of weak Gaussian noise, rare exit events can be analyzed by utilizing Freidlin-Wentzell large deviation theory. In this section, we first review the concepts of the most probable exit path, quasipotential, and mean first exit time, and derive their asymptotic expressions. Then we design a machine learning method to compute these quantities numerically.

3.1 Large deviation theory

Now we reformulate the stochastic vegetation-water system (2.1) into the following form

dx(t)=b(x)dt+σ(x)dBε(t),𝑑𝑥𝑡𝑏𝑥𝑑𝑡𝜎𝑥𝑑superscript𝐵𝜀𝑡dx(t)=b(x)dt+\sigma(x)dB^{\varepsilon}(t),italic_d italic_x ( italic_t ) = italic_b ( italic_x ) italic_d italic_t + italic_σ ( italic_x ) italic_d italic_B start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT ( italic_t ) ,

where Bε(t)=(B1ε(t),B2ε(t))superscript𝐵𝜀𝑡superscriptsuperscriptsubscript𝐵1𝜀𝑡superscriptsubscript𝐵2𝜀𝑡topB^{\varepsilon}(t)=(B_{1}^{\varepsilon}(t),B_{2}^{\varepsilon}(t))^{\top}italic_B start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT ( italic_t ) = ( italic_B start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT ( italic_t ) , italic_B start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT ( italic_t ) ) start_POSTSUPERSCRIPT ⊤ end_POSTSUPERSCRIPT is a two-dimensional Brownian motion. Given the scenario of weak noise intensity, the stationary distribution ps(x)subscript𝑝𝑠𝑥p_{s}(x)italic_p start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ( italic_x ) of the stochastic system can be assumed to have the following WKB form

ps(x)C(x)exp{ε1V(x)},similar-tosubscript𝑝𝑠𝑥𝐶𝑥superscript𝜀1𝑉𝑥p_{s}(x)\sim C(x)\exp\{-\varepsilon^{-1}V(x)\},italic_p start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ( italic_x ) ∼ italic_C ( italic_x ) roman_exp { - italic_ε start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT italic_V ( italic_x ) } ,

where C(x)𝐶𝑥C(x)italic_C ( italic_x ) is referred to as exponential prefactor, and V(x)𝑉𝑥V(x)italic_V ( italic_x ) is called quasipotential, which characterizes the possibility of the stochastic state fluctuating into the vicinity of the specific point x𝑥xitalic_x. In addition to the stationary distribution, the quasipotential also exponentially dominates the magnitude of mean first exit time and exit location distribution. According to Freidlin-Wentzell large deviation theory, the quasipotential is defined as the minimum of the action functional along the absolutely continuous path connecting the fixed point x¯¯𝑥\bar{x}over¯ start_ARG italic_x end_ARG and the specific point x𝑥xitalic_x, in the sense that

V(x):=infT>0infφC[0,T]{S(φ):φ(0)=x¯,φ(T)=x},assign𝑉𝑥subscriptinfimum𝑇0subscriptinfimum𝜑𝐶0𝑇conditional-set𝑆𝜑formulae-sequence𝜑0¯𝑥𝜑𝑇𝑥V(x):=\inf_{T>0}\inf_{\varphi\in C[0,T]}\{S(\varphi):\varphi(0)=\bar{x},\,% \varphi(T)=x\},italic_V ( italic_x ) := roman_inf start_POSTSUBSCRIPT italic_T > 0 end_POSTSUBSCRIPT roman_inf start_POSTSUBSCRIPT italic_φ ∈ italic_C [ 0 , italic_T ] end_POSTSUBSCRIPT { italic_S ( italic_φ ) : italic_φ ( 0 ) = over¯ start_ARG italic_x end_ARG , italic_φ ( italic_T ) = italic_x } ,

where the action functional S(φ)𝑆𝜑S(\varphi)italic_S ( italic_φ ) has the following form

S(φ)=120T(φ˙b(φ))a1(φ˙b(φ))𝑑t.𝑆𝜑12superscriptsubscript0𝑇superscript˙𝜑𝑏𝜑topsuperscript𝑎1˙𝜑𝑏𝜑differential-d𝑡S(\varphi)=\frac{1}{2}\int_{0}^{T}(\dot{\varphi}-b(\varphi))^{\top}a^{-1}(\dot% {\varphi}-b(\varphi))dt.italic_S ( italic_φ ) = divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT ( over˙ start_ARG italic_φ end_ARG - italic_b ( italic_φ ) ) start_POSTSUPERSCRIPT ⊤ end_POSTSUPERSCRIPT italic_a start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ( over˙ start_ARG italic_φ end_ARG - italic_b ( italic_φ ) ) italic_d italic_t .

By substituting the WKB approximation into the stationary Fokker-Planck equation and combining the lowest order terms of ε𝜀\varepsilonitalic_ε, we obtain Hamilton-Jacobi equation

V(x),b(x)+12V(x),a(x)V(x)=0.𝑉𝑥𝑏𝑥12𝑉𝑥𝑎𝑥𝑉𝑥0\langle\nabla V(x),b(x)\rangle+\frac{1}{2}\langle\nabla V(x),a(x)\nabla V(x)% \rangle=0.⟨ ∇ italic_V ( italic_x ) , italic_b ( italic_x ) ⟩ + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ⟨ ∇ italic_V ( italic_x ) , italic_a ( italic_x ) ∇ italic_V ( italic_x ) ⟩ = 0 .

Note that the above equation has a geometric interpretation. By combining the two terms on the left, we gain

V(x),b(x)+12a(x)V(x)=0.𝑉𝑥𝑏𝑥12𝑎𝑥𝑉𝑥0\langle\nabla V(x),b(x)+\frac{1}{2}a(x)\nabla V(x)\rangle=0.⟨ ∇ italic_V ( italic_x ) , italic_b ( italic_x ) + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x ) ∇ italic_V ( italic_x ) ⟩ = 0 .

Therefore, we get the following orthogonal relationship

V(x)b(x)+12a(x)V(x).perpendicular-to𝑉𝑥𝑏𝑥12𝑎𝑥𝑉𝑥\nabla V(x)\perp b(x)+\frac{1}{2}a(x)\nabla V(x).∇ italic_V ( italic_x ) ⟂ italic_b ( italic_x ) + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x ) ∇ italic_V ( italic_x ) .

Define l(x):=b(x)+12a(x)V(x)assign𝑙𝑥𝑏𝑥12𝑎𝑥𝑉𝑥l(x):=b(x)+\frac{1}{2}a(x)\nabla V(x)italic_l ( italic_x ) := italic_b ( italic_x ) + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x ) ∇ italic_V ( italic_x ). The vector field b(x)𝑏𝑥b(x)italic_b ( italic_x ) has the following decomposition

b(x)=12a(x)V(x)+l(x).𝑏𝑥12𝑎𝑥𝑉𝑥𝑙𝑥b(x)=-\frac{1}{2}a(x)\nabla V(x)+l(x).italic_b ( italic_x ) = - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x ) ∇ italic_V ( italic_x ) + italic_l ( italic_x ) .

The next-to-leading-order approximation of the WKB approximation to the stationary Fokker-Planck equation leads to the following transport equation about the exponential prefactor C(x)𝐶𝑥C(x)italic_C ( italic_x ):

C,b+aV+C(divb+12a:H(V)+A,V)=0,\langle\nabla C,b+a\nabla V\rangle+C\big{(}\text{div}b+\frac{1}{2}a:H(V)+% \langle A,\nabla V\rangle\big{)}=0,⟨ ∇ italic_C , italic_b + italic_a ∇ italic_V ⟩ + italic_C ( div italic_b + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a : italic_H ( italic_V ) + ⟨ italic_A , ∇ italic_V ⟩ ) = 0 ,

where H𝐻Hitalic_H denotes Hessian matrix, a:H=i,jaijHij:𝑎𝐻subscript𝑖𝑗subscript𝑎𝑖𝑗subscript𝐻𝑖𝑗a:H=\sum_{i,j}a_{ij}H_{ij}italic_a : italic_H = ∑ start_POSTSUBSCRIPT italic_i , italic_j end_POSTSUBSCRIPT italic_a start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT, and the vector A(x)𝐴𝑥A(x)italic_A ( italic_x ) is defined by Ai(x)=j=12aijxjsubscript𝐴𝑖𝑥superscriptsubscript𝑗12subscript𝑎𝑖𝑗subscript𝑥𝑗A_{i}(x)=\sum\limits_{j=1}^{2}\frac{\partial a_{ij}}{\partial x_{j}}italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_x ) = ∑ start_POSTSUBSCRIPT italic_j = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT divide start_ARG ∂ italic_a start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_x start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_ARG, i.e.,

A(x)=(4σ12x130).𝐴𝑥4superscriptsubscript𝜎12superscriptsubscript𝑥130A(x)=\left(\begin{array}[]{c}4\sigma_{1}^{2}x_{1}^{3}\\ 0\\ \end{array}\right).italic_A ( italic_x ) = ( start_ARRAY start_ROW start_CELL 4 italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL end_ROW end_ARRAY ) .

Putting the decomposition b(x)=12a(x)V(x)+l(x)𝑏𝑥12𝑎𝑥𝑉𝑥𝑙𝑥b(x)=-\frac{1}{2}a(x)\nabla V(x)+l(x)italic_b ( italic_x ) = - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x ) ∇ italic_V ( italic_x ) + italic_l ( italic_x ) into the transport equation yields the following result

C,12aV+l+C(divl+12A,V)=0.𝐶12𝑎𝑉𝑙𝐶div𝑙12𝐴𝑉0\langle\nabla C,\frac{1}{2}a\nabla V+l\rangle+C(\text{div}l+\frac{1}{2}\langle A% ,\nabla V\rangle)=0.⟨ ∇ italic_C , divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ∇ italic_V + italic_l ⟩ + italic_C ( div italic_l + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ⟨ italic_A , ∇ italic_V ⟩ ) = 0 .

The above equation can be rewritten as

lnC,12aV+l=F,𝐶12𝑎𝑉𝑙𝐹\langle\nabla\ln C,\frac{1}{2}a\nabla V+l\rangle=-F,⟨ ∇ roman_ln italic_C , divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ∇ italic_V + italic_l ⟩ = - italic_F , (3.3)

where

F(x)=divl(x)+12A(x),V(x).𝐹𝑥div𝑙𝑥12𝐴𝑥𝑉𝑥F(x)=\text{div}l(x)+\frac{1}{2}\langle A(x),\nabla V(x)\rangle.italic_F ( italic_x ) = div italic_l ( italic_x ) + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ⟨ italic_A ( italic_x ) , ∇ italic_V ( italic_x ) ⟩ .

According to Freidlin-Wentzell large deviation theory, the most probable path of fluctuation dynamics satisfies the following equation

φ˙x(t)=b(φx(t))+a(φx(t))V(φx(t))=12a(φx(t))V(φx(t))+l(φx(t)).superscript˙𝜑𝑥𝑡𝑏superscript𝜑𝑥𝑡𝑎superscript𝜑𝑥𝑡𝑉superscript𝜑𝑥𝑡12𝑎superscript𝜑𝑥𝑡𝑉superscript𝜑𝑥𝑡𝑙superscript𝜑𝑥𝑡\dot{\varphi}^{x}(t)=b(\varphi^{x}(t))+a(\varphi^{x}(t))\nabla V(\varphi^{x}(t% ))=\frac{1}{2}a(\varphi^{x}(t))\nabla V(\varphi^{x}(t))+l(\varphi^{x}(t)).over˙ start_ARG italic_φ end_ARG start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) = italic_b ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) + italic_a ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) ∇ italic_V ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) = divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) ∇ italic_V ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) + italic_l ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) .

This can also be confirmed by the results of the method of characteristics applied to the Hamilton-Jacobi equation. Therefore, along the most probable path, the left-hand side of equation (3.3) is transformed into the complete differential of lnC(x)𝐶𝑥\ln C(x)roman_ln italic_C ( italic_x ), i.e.,

ddtlnC(φx(t))=F(φx(t)).𝑑𝑑𝑡𝐶superscript𝜑𝑥𝑡𝐹superscript𝜑𝑥𝑡\frac{d}{dt}\ln C(\varphi^{x}(t))=-F(\varphi^{x}(t)).divide start_ARG italic_d end_ARG start_ARG italic_d italic_t end_ARG roman_ln italic_C ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) = - italic_F ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) .

So that we can integrate the prefactor function as

C(x)exp{0F(φx(t))𝑑t}.similar-to𝐶𝑥superscriptsubscript0𝐹superscript𝜑𝑥𝑡differential-d𝑡C(x)\sim\exp\{-\int_{-\infty}^{0}F(\varphi^{x}(t))dt\}.italic_C ( italic_x ) ∼ roman_exp { - ∫ start_POSTSUBSCRIPT - ∞ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT italic_F ( italic_φ start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ( italic_t ) ) italic_d italic_t } . (3.4)

If x𝑥xitalic_x is a saddle point, then the upper limit of the above integral is infinite.

Assume that D𝐷Ditalic_D is the attraction domain of SN2. Define the first exit time of the stochastic state from D𝐷Ditalic_D as

τDε=inf{t0:x(t)D,x(0)=SN2}.subscriptsuperscript𝜏𝜀𝐷infimumconditional-set𝑡0formulae-sequence𝑥𝑡𝐷𝑥0SN2\tau^{\varepsilon}_{D}=\inf\{t\geq 0:x(t)\notin D,x(0)=\text{SN2}\}.italic_τ start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT = roman_inf { italic_t ≥ 0 : italic_x ( italic_t ) ∉ italic_D , italic_x ( 0 ) = SN2 } .

It is a random time, and its average quantity is called the mean first exit time. This can provide important quantitative information regarding the disappearance of vegetation. According to large deviation theory, the mean first exit time is exponentially dominated by the minimal value of the quasipotential along the boundary D𝐷\partial D∂ italic_D:

limε0εln𝔼τDε=infxDV(x).subscript𝜀0𝜀𝔼subscriptsuperscript𝜏𝜀𝐷subscriptinfimum𝑥𝐷𝑉𝑥\lim_{\varepsilon\rightarrow 0}\varepsilon\ln\mathbb{E}\tau^{\varepsilon}_{D}=% \inf_{x\in\partial D}V(x).roman_lim start_POSTSUBSCRIPT italic_ε → 0 end_POSTSUBSCRIPT italic_ε roman_ln blackboard_E italic_τ start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT = roman_inf start_POSTSUBSCRIPT italic_x ∈ ∂ italic_D end_POSTSUBSCRIPT italic_V ( italic_x ) .

Usually, this minimization is achieved at the saddle US. Then we have

𝔼τDε=LDεexp{V(US)}.𝔼subscriptsuperscript𝜏𝜀𝐷subscriptsuperscript𝐿𝜀𝐷𝑉US\mathbb{E}\tau^{\varepsilon}_{D}=L^{\varepsilon}_{D}\exp\{V(\text{US})\}.blackboard_E italic_τ start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT = italic_L start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT roman_exp { italic_V ( US ) } .

In general, We can calculate the prefactor LDεsuperscriptsubscript𝐿𝐷𝜀L_{D}^{\varepsilon}italic_L start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT in two distinct scenarios.

Case A. non-characteristic boundary (see BR22 ; BR for reference).

We assume that the domain D𝐷Ditalic_D is an open, smooth and connected subset of nsuperscript𝑛\mathbb{R}^{n}blackboard_R start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT that satisfies the following conditions, where n(y)𝑛𝑦n(y)italic_n ( italic_y ) represents the exterior normal vector at yD𝑦𝐷y\in\partial Ditalic_y ∈ ∂ italic_D.

(A1)

The deterministic system x˙=b(x)˙𝑥𝑏𝑥\dot{x}=b(x)over˙ start_ARG italic_x end_ARG = italic_b ( italic_x ) has a unique fixed point x¯¯𝑥\bar{x}over¯ start_ARG italic_x end_ARG within D𝐷Ditalic_D that attracts all the trajectories originating from D𝐷Ditalic_D. Additionally, the inner product between b(y)𝑏𝑦b(y)italic_b ( italic_y ) and n(y)𝑛𝑦n(y)italic_n ( italic_y ) is negative for all y𝑦yitalic_y on the boundary D𝐷\partial D∂ italic_D, i.e., b(y),n(y)<0𝑏𝑦𝑛𝑦0\langle b(y),n(y)\rangle<0⟨ italic_b ( italic_y ) , italic_n ( italic_y ) ⟩ < 0 for all yD𝑦𝐷y\in\partial Ditalic_y ∈ ∂ italic_D.

(A2)

The function V𝑉Vitalic_V is continuously differentiable (C1superscript𝐶1C^{1}italic_C start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT) in D𝐷Ditalic_D; for any xD¯𝑥¯𝐷x\in\bar{D}italic_x ∈ over¯ start_ARG italic_D end_ARG, the most probable path φtxsuperscriptsubscript𝜑𝑡𝑥\varphi_{t}^{x}italic_φ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT approaches x¯¯𝑥\bar{x}over¯ start_ARG italic_x end_ARG as t𝑡t\rightarrow-\inftyitalic_t → - ∞; and 12a(y)V(y)+l(y),n(y)>012𝑎𝑦𝑉𝑦𝑙𝑦𝑛𝑦0\langle\frac{1}{2}a(y)\nabla V(y)+l(y),n(y)\rangle>0⟨ divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_y ) ∇ italic_V ( italic_y ) + italic_l ( italic_y ) , italic_n ( italic_y ) ⟩ > 0 for all yD𝑦𝐷y\in\partial Ditalic_y ∈ ∂ italic_D.

(A3)

The minimum of V𝑉Vitalic_V over D𝐷\partial D∂ italic_D is attained at a single point x*superscript𝑥x^{*}italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT. At this point,

μ*=12a(x*)V(x*)+l(x*),n(x*)>0,superscript𝜇12𝑎superscript𝑥𝑉superscript𝑥𝑙superscript𝑥𝑛superscript𝑥0\mu^{*}=\langle\frac{1}{2}a(x^{*})\nabla V(x^{*})+l(x^{*}),n(x^{*})\rangle>0,italic_μ start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT = ⟨ divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ) ∇ italic_V ( italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ) + italic_l ( italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ) , italic_n ( italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ) ⟩ > 0 ,

and the quadratic form h*:ξξ,2V(x*)ξ:superscriptmaps-to𝜉𝜉superscript2𝑉superscript𝑥𝜉h^{*}:\xi\mapsto\langle\xi,\nabla^{2}V(x^{*})\xi\rangleitalic_h start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT : italic_ξ ↦ ⟨ italic_ξ , ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_V ( italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ) italic_ξ ⟩ has positive eigenvalues on the hyperplane n(x*)={ξn:ξ,n(x*)=0}𝑛superscriptsuperscript𝑥perpendicular-toconditional-set𝜉superscript𝑛𝜉𝑛superscript𝑥0n(x^{*})^{\perp}=\{\xi\in\mathbb{R}^{n}:\langle\xi,n(x^{*})\rangle=0\}italic_n ( italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT ⟂ end_POSTSUPERSCRIPT = { italic_ξ ∈ blackboard_R start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT : ⟨ italic_ξ , italic_n ( italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ) ⟩ = 0 }.

If b(y),n(y)<0𝑏𝑦𝑛𝑦0\langle b(y),n(y)\rangle<0⟨ italic_b ( italic_y ) , italic_n ( italic_y ) ⟩ < 0 for all yD𝑦𝐷y\in\partial Ditalic_y ∈ ∂ italic_D, the boundary is designated as non-characteristic. This condition guarantees that the dynamical trajectories, originating from the closure D¯¯𝐷\bar{D}over¯ start_ARG italic_D end_ARG, will remain confined within D𝐷Ditalic_D, and that the vector field is perpendicular to the boundary. Leveraging Assumption (A1), we derive the integral formula

λDε=xD12a(x)V(x)+l(x),n(x)C(x)exp{ε1V(x)}𝑑xsuperscriptsubscript𝜆𝐷𝜀subscript𝑥𝐷12𝑎𝑥𝑉𝑥𝑙𝑥𝑛𝑥𝐶𝑥superscript𝜀1𝑉𝑥differential-d𝑥\lambda_{D}^{\varepsilon}=\int\limits_{x\in\partial D}\langle\frac{1}{2}a(x)% \nabla V(x)+l(x),n(x)\rangle C(x)\exp\{-\varepsilon^{-1}V(x)\}dxitalic_λ start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT = ∫ start_POSTSUBSCRIPT italic_x ∈ ∂ italic_D end_POSTSUBSCRIPT ⟨ divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x ) ∇ italic_V ( italic_x ) + italic_l ( italic_x ) , italic_n ( italic_x ) ⟩ italic_C ( italic_x ) roman_exp { - italic_ε start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT italic_V ( italic_x ) } italic_d italic_x

for the exit rate λDε=[𝔼τDε]1superscriptsubscript𝜆𝐷𝜀superscriptdelimited-[]𝔼superscriptsubscript𝜏𝐷𝜀1\lambda_{D}^{\varepsilon}=[\mathbb{E}\tau_{D}^{\varepsilon}]^{-1}italic_λ start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT = [ blackboard_E italic_τ start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT ] start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. By employing the second-order expansion of the potential V𝑉Vitalic_V in the vicinity of the point x*superscript𝑥x^{*}italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT, we obtain an equivalent relation for the prefactor

LDεsuperscriptsubscript𝐿𝐷𝜀\displaystyle L_{D}^{\varepsilon}italic_L start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT 1C(x)μdeth(2πε)n1similar-toabsent1𝐶superscript𝑥superscript𝜇superscriptsuperscript2𝜋𝜀𝑛1\displaystyle\sim\frac{1}{C(x^{\ast})\mu^{\ast}}\sqrt{\frac{\det h^{\ast}}{(2% \pi\varepsilon)^{n-1}}}∼ divide start_ARG 1 end_ARG start_ARG italic_C ( italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) italic_μ start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT end_ARG square-root start_ARG divide start_ARG roman_det italic_h start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT end_ARG start_ARG ( 2 italic_π italic_ε ) start_POSTSUPERSCRIPT italic_n - 1 end_POSTSUPERSCRIPT end_ARG end_ARG (3.5)
1μdeth(2πε)n1exp{0F(φx(t))𝑑t},similar-toabsent1superscript𝜇superscriptsuperscript2𝜋𝜀𝑛1superscriptsubscript0𝐹superscript𝜑superscript𝑥𝑡differential-d𝑡\displaystyle\sim\frac{1}{\mu^{\ast}}\sqrt{\frac{\det h^{\ast}}{(2\pi% \varepsilon)^{n-1}}}\exp\{\int_{\infty}^{0}F(\varphi^{x^{\ast}}(t))dt\},∼ divide start_ARG 1 end_ARG start_ARG italic_μ start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT end_ARG square-root start_ARG divide start_ARG roman_det italic_h start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT end_ARG start_ARG ( 2 italic_π italic_ε ) start_POSTSUPERSCRIPT italic_n - 1 end_POSTSUPERSCRIPT end_ARG end_ARG roman_exp { ∫ start_POSTSUBSCRIPT ∞ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT italic_F ( italic_φ start_POSTSUPERSCRIPT italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT end_POSTSUPERSCRIPT ( italic_t ) ) italic_d italic_t } ,

where we utilize the approximation in (3.4).

Case B. characteristic boundary BR22 ; BR .

The basin domain D𝐷Ditalic_D is characteristic in the sense that the inner product between the vector field b(y)𝑏𝑦b(y)italic_b ( italic_y ) and the exterior normal vector n(y)𝑛𝑦n(y)italic_n ( italic_y ) vanishes for all points y𝑦yitalic_y within the domain D𝐷Ditalic_D, i.e., b(y),n(y)=0𝑏𝑦𝑛𝑦0\langle b(y),n(y)\rangle=0⟨ italic_b ( italic_y ) , italic_n ( italic_y ) ⟩ = 0 for all yD𝑦𝐷y\in Ditalic_y ∈ italic_D. We consider the metastable scenario that the deterministic system x˙=b(x)˙𝑥𝑏𝑥\dot{x}=b(x)over˙ start_ARG italic_x end_ARG = italic_b ( italic_x ) possesses two stable fixed points x¯1subscript¯𝑥1\bar{x}_{1}over¯ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and x¯2subscript¯𝑥2\bar{x}_{2}over¯ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, whose respective basins of attraction are divided by a smooth hypersurface S𝑆Sitalic_S. We focus on exit events from the basin of attraction D𝐷Ditalic_D associated with x¯1subscript¯𝑥1\bar{x}_{1}over¯ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and introduce the following set of assumptions.

(B1)

All trajectories of the deterministic system x˙=b(x)˙𝑥𝑏𝑥\dot{x}=b(x)over˙ start_ARG italic_x end_ARG = italic_b ( italic_x ) initiated on the hypersurface S𝑆Sitalic_S remain confined to S𝑆Sitalic_S and ultimately converge to a single fixed point xSsuperscript𝑥𝑆x^{\ast}\in Sitalic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ∈ italic_S; furthermore, the Jacobi matrix b(x)𝑏superscript𝑥\nabla b(x^{\ast})∇ italic_b ( italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) possesses n1𝑛1n-1italic_n - 1 eigenvalues with negative real part and a single positive eigenvalue denoted by λ*superscript𝜆\lambda^{*}italic_λ start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT.

(B2)

With respect to the quasipotential V𝑉Vitalic_V associated with x¯1subscript¯𝑥1\bar{x}_{1}over¯ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT, there exists a unique (up to time shift) trajectory ρ=(ρt)tD𝜌subscriptsubscript𝜌𝑡𝑡𝐷\rho=(\rho_{t})_{t\in\mathbb{R}}\subset Ditalic_ρ = ( italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) start_POSTSUBSCRIPT italic_t ∈ blackboard_R end_POSTSUBSCRIPT ⊂ italic_D such that

limtρt=x¯1,limt+ρt=x*,andV(x*)=𝒮,+[ρ].formulae-sequencesubscript𝑡subscript𝜌𝑡subscript¯𝑥1formulae-sequencesubscript𝑡subscript𝜌𝑡superscript𝑥and𝑉superscript𝑥subscript𝒮delimited-[]𝜌\lim_{t\rightarrow-\infty}\rho_{t}=\bar{x}_{1},\quad\lim_{t\rightarrow+\infty}% \rho_{t}=x^{*},\quad\text{and}\quad V(x^{*})=\mathcal{S}_{-\infty,+\infty}[% \rho].roman_lim start_POSTSUBSCRIPT italic_t → - ∞ end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = over¯ start_ARG italic_x end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , roman_lim start_POSTSUBSCRIPT italic_t → + ∞ end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT , and italic_V ( italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT ) = caligraphic_S start_POSTSUBSCRIPT - ∞ , + ∞ end_POSTSUBSCRIPT [ italic_ρ ] .
(B3)

The quasipotential V𝑉Vitalic_V is smooth in a neighborhood of ρ=(ρt)t𝜌subscriptsubscript𝜌𝑡𝑡\rho=(\rho_{t})_{t\in\mathbb{R}}italic_ρ = ( italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) start_POSTSUBSCRIPT italic_t ∈ blackboard_R end_POSTSUBSCRIPT. Additionally, the vector field l𝑙litalic_l defined by l(x)=b(x)+12a(x)V(x)𝑙𝑥𝑏𝑥12𝑎𝑥𝑉𝑥l(x)=b(x)+\frac{1}{2}a(x)\nabla V(x)italic_l ( italic_x ) = italic_b ( italic_x ) + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x ) ∇ italic_V ( italic_x ) satisfies the orthogonality relation V(x),l(x)=0𝑉𝑥𝑙𝑥0\langle\nabla V(x),l(x)\rangle=0⟨ ∇ italic_V ( italic_x ) , italic_l ( italic_x ) ⟩ = 0.

In this context, the quasipotential V𝑉Vitalic_V attains its minimum on the hypersurface S𝑆Sitalic_S precisely at the point x*superscript𝑥x^{*}italic_x start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT. Moreover, the trajectory ρ𝜌\rhoitalic_ρ is designated as the most probable exit path, satisfying the differential equation

ρ˙t=12a(ρt)V(ρt)+l(ρt),t.formulae-sequencesubscript˙𝜌𝑡12𝑎subscript𝜌𝑡𝑉subscript𝜌𝑡𝑙subscript𝜌𝑡for-all𝑡\dot{\rho}_{t}=\frac{1}{2}a(\rho_{t})\nabla V(\rho_{t})+l(\rho_{t}),\quad% \forall t\in\mathbb{R}.over˙ start_ARG italic_ρ end_ARG start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) ∇ italic_V ( italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) + italic_l ( italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) , ∀ italic_t ∈ blackboard_R .

For any given t𝑡t\in\mathbb{R}italic_t ∈ blackboard_R, this path coincides with the trajectory (φsx)s0subscriptsuperscriptsubscript𝜑𝑠𝑥𝑠0(\varphi_{s}^{x})_{s\leq 0}( italic_φ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT ) start_POSTSUBSCRIPT italic_s ≤ 0 end_POSTSUBSCRIPT that connects x¯¯𝑥\bar{x}over¯ start_ARG italic_x end_ARG to x=ρt𝑥subscript𝜌𝑡x=\rho_{t}italic_x = italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT, according to the relation

φsx=ρs+t,s0.formulae-sequencesuperscriptsubscript𝜑𝑠𝑥subscript𝜌𝑠𝑡for-all𝑠0\varphi_{s}^{x}=\rho_{s+t},\quad\forall s\leq 0.italic_φ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_x end_POSTSUPERSCRIPT = italic_ρ start_POSTSUBSCRIPT italic_s + italic_t end_POSTSUBSCRIPT , ∀ italic_s ≤ 0 .

To describe the prefactor LDϵsuperscriptsubscript𝐿𝐷italic-ϵL_{D}^{\epsilon}italic_L start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ϵ end_POSTSUPERSCRIPT in this scenario, we formulate the following supplementary assumption.

(B4)

The matrix H=limt+2V(ρt)superscript𝐻subscript𝑡superscript2𝑉subscript𝜌𝑡H^{\ast}=\lim_{t\rightarrow+\infty}\nabla^{2}V(\rho_{t})italic_H start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = roman_lim start_POSTSUBSCRIPT italic_t → + ∞ end_POSTSUBSCRIPT ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_V ( italic_ρ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) exists and possesses n1𝑛1n-1italic_n - 1 positive eigenvalues and a single negative eigenvalue.

Relying on these four assumptions, an asymptotic formula for estimating the expected time taken by the process to exit the domain D𝐷Ditalic_D is given by

LDεπλ|detH|det2V(x¯)exp{F(φ(t))𝑑t}.similar-tosuperscriptsubscript𝐿𝐷𝜀𝜋superscript𝜆superscript𝐻superscript2𝑉¯𝑥superscriptsubscript𝐹𝜑𝑡differential-d𝑡L_{D}^{\varepsilon}\sim\frac{\pi}{\lambda^{\ast}}\sqrt{\frac{|\det H^{\ast}|}{% \det\nabla^{2}V(\bar{x})}}\exp\Big{\{}\int_{-\infty}^{\infty}F(\varphi(t))dt% \Big{\}}.italic_L start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ε end_POSTSUPERSCRIPT ∼ divide start_ARG italic_π end_ARG start_ARG italic_λ start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT end_ARG square-root start_ARG divide start_ARG | roman_det italic_H start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT | end_ARG start_ARG roman_det ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_V ( over¯ start_ARG italic_x end_ARG ) end_ARG end_ARG roman_exp { ∫ start_POSTSUBSCRIPT - ∞ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_F ( italic_φ ( italic_t ) ) italic_d italic_t } . (3.6)

This expression provides an approximation for the expected exit time from the domain D𝐷Ditalic_D under the given conditions.

3.2 Machine learning algorithm

It is seen in subsection 3.1 that the computations of the most probable path and the mean first exit time depend on the results of quasipotential and rotational component, i.e., the decomposition of the vector field. In this subsection, we aim to propose a machine learning method to compute these quantities of the stochastic vegetation-water system (2.1) based on this decomposition.

We design a neural network architecture to achieve this goal. The input of the network is set as the coordinate x=(x1,x2)𝑥superscriptsubscript𝑥1subscript𝑥2topx=(x_{1},x_{2})^{\top}italic_x = ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT ⊤ end_POSTSUPERSCRIPT. The output of the network is (V^θ,lθ)n+1subscript^𝑉𝜃subscript𝑙𝜃superscript𝑛1(\hat{V}_{\theta},l_{\theta})\in\mathbb{R}^{n+1}( over^ start_ARG italic_V end_ARG start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , italic_l start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ) ∈ blackboard_R start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT. The quasipotential function is defined as Vθ(x)=V^θ(x)+|xx¯|2subscript𝑉𝜃𝑥subscript^𝑉𝜃𝑥superscript𝑥¯𝑥2V_{\theta}(x)=\hat{V}_{\theta}(x)+|x-\bar{x}|^{2}italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x ) = over^ start_ARG italic_V end_ARG start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x ) + | italic_x - over¯ start_ARG italic_x end_ARG | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT to guarantee the unboundedness of V(x)𝑉𝑥V(x)italic_V ( italic_x ) and |V(x)|𝑉𝑥|\nabla V(x)|| ∇ italic_V ( italic_x ) |. Here, x¯¯𝑥\bar{x}over¯ start_ARG italic_x end_ARG is the stable fixed point SN2, and θ𝜃\thetaitalic_θ denotes the training parameters of the neural network.

In order to train the neural network, we choose N𝑁Nitalic_N points randomly in the attraction domain of SN2 and construct a loss function as follows:

L=Ldyn+λ1Lorth+λ2L0.𝐿subscript𝐿dynsubscript𝜆1subscript𝐿orthsubscript𝜆2subscript𝐿0L=L_{\text{dyn}}+\lambda_{1}L_{\text{orth}}+\lambda_{2}L_{0}.italic_L = italic_L start_POSTSUBSCRIPT dyn end_POSTSUBSCRIPT + italic_λ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_L start_POSTSUBSCRIPT orth end_POSTSUBSCRIPT + italic_λ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT .

Since the vector field has the decomposition b(x)=12a(x)V(x)+l(x)𝑏𝑥12𝑎𝑥𝑉𝑥𝑙𝑥b(x)=-\frac{1}{2}a(x)\nabla V(x)+l(x)italic_b ( italic_x ) = - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x ) ∇ italic_V ( italic_x ) + italic_l ( italic_x ), the first part Ldynsubscript𝐿dynL_{\text{dyn}}italic_L start_POSTSUBSCRIPT dyn end_POSTSUBSCRIPT of the loss function can be set as

Ldyn=1Ni=1N[b(xi)+12a(xi)Vθ(xi)lθ(xi)]2,subscript𝐿dyn1𝑁superscriptsubscript𝑖1𝑁superscriptdelimited-[]𝑏subscript𝑥𝑖12𝑎subscript𝑥𝑖subscript𝑉𝜃subscript𝑥𝑖subscript𝑙𝜃subscript𝑥𝑖2L_{\text{dyn}}=\frac{1}{N}\sum_{i=1}^{N}[b(x_{i})+\frac{1}{2}a(x_{i})\nabla V_% {\theta}(x_{i})-l_{\theta}(x_{i})]^{2},italic_L start_POSTSUBSCRIPT dyn end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG italic_N end_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT [ italic_b ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_a ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ∇ italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) - italic_l start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ,

where the gradient of the quasipotential is realized by automatic differentiation technique. Due to the orthogonal relation V(x)l(x)bottom𝑉𝑥𝑙𝑥\nabla V(x)\bot l(x)∇ italic_V ( italic_x ) ⊥ italic_l ( italic_x ), the second part Lorthsubscript𝐿orthL_{\text{orth}}italic_L start_POSTSUBSCRIPT orth end_POSTSUBSCRIPT of the loss function is assigned as

Lorth=1Ni=1N[Vθ(xi)lθ(xi)]2|Vθ(xi)|2|lθ(xi)|2+δ.subscript𝐿orth1𝑁superscriptsubscript𝑖1𝑁superscriptdelimited-[]subscript𝑉𝜃subscript𝑥𝑖subscript𝑙𝜃subscript𝑥𝑖2superscriptsubscript𝑉𝜃subscript𝑥𝑖2superscriptsubscript𝑙𝜃subscript𝑥𝑖2𝛿L_{\text{orth}}=\frac{1}{N}\sum_{i=1}^{N}\frac{[\nabla V_{\theta}(x_{i})\cdot l% _{\theta}(x_{i})]^{2}}{|\nabla V_{\theta}(x_{i})|^{2}|l_{\theta}(x_{i})|^{2}+% \delta}.italic_L start_POSTSUBSCRIPT orth end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG italic_N end_ARG ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT divide start_ARG [ ∇ italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ⋅ italic_l start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG | ∇ italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_l start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_δ end_ARG .

Here, the small parameter δ1much-less-than𝛿1\delta\ll 1italic_δ ≪ 1 is chosen to avoid a zero denominator, ensuring numerical stability. Besides, the third part L0subscript𝐿0L_{0}italic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT of the loss function is set as L0=Vθ(x¯)2subscript𝐿0subscript𝑉𝜃superscript¯𝑥2L_{0}=V_{\theta}(\bar{x})^{2}italic_L start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( over¯ start_ARG italic_x end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT to guarantee the fact that the quasipotential of the stable fixed point x¯=SN2¯𝑥SN2\bar{x}=\text{SN2}over¯ start_ARG italic_x end_ARG = SN2 is zero. In this paper, we choose the weight parameters λ1=1subscript𝜆11\lambda_{1}=1italic_λ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 1 and λ2=0.1subscript𝜆20.1\lambda_{2}=0.1italic_λ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 0.1, which are used to balance the three parts of loss function.

After training the neural network, we obtain the quasipotential function Vθ(x)subscript𝑉𝜃𝑥V_{\theta}(x)italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x ) and the rotational component l(x)𝑙𝑥l(x)italic_l ( italic_x ). Thus the most probable path can be integrated by the equation

x˙=b(x)+a(x)Vθ(x)˙𝑥𝑏𝑥𝑎𝑥subscript𝑉𝜃𝑥\dot{x}=b(x)+a(x)\nabla V_{\theta}(x)over˙ start_ARG italic_x end_ARG = italic_b ( italic_x ) + italic_a ( italic_x ) ∇ italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x )

in reverse time, starting from the end point. Additionally, the mean first exit time can be computed using the asymptotic expression provided in subsection 3.1.

4 Results

In this section, we present the results of applying the proposed machine learning algorithm to the stochastic vegetation-water system (2.1). The hyperparameters we set are as follows: The neural network is configured with 6 hidden layers, each having 20 nodes. The activation function in hidden layers is chosen as tanh\tanhroman_tanh, while the output layer uses the identify function. We utilize the Adam optimizer with a learning rate of 0.001. The small parameter in the loss function is set to δ=0.001𝛿0.001\delta=0.001italic_δ = 0.001. The neural network is trained for 1000000 epochs.

In Section 2, upon observing the expressions of the noise intensities σ1=ρσ~1/Ksubscript𝜎1𝜌subscript~𝜎1𝐾\sigma_{1}=-\rho\tilde{\sigma}_{1}/Kitalic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = - italic_ρ over~ start_ARG italic_σ end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_K and σ2=Rσ~2subscript𝜎2𝑅subscript~𝜎2\sigma_{2}=R\tilde{\sigma}_{2}italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_R over~ start_ARG italic_σ end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT with ρ=1𝜌1\rho=1italic_ρ = 1, K=10𝐾10K=10italic_K = 10, and R=1.55𝑅1.55R=1.55italic_R = 1.55, it is worth noting that the order of magnitude of σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT should be significantly smaller than that of σ2subscript𝜎2\sigma_{2}italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT in practical system. Consequently, for this context, we assume the values of σ1=0.1subscript𝜎10.1\sigma_{1}=0.1italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 0.1 and σ2=1subscript𝜎21\sigma_{2}=1italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 1.

Refer to caption
Figure 4: The loss function is decreasing with respect to the increasing number of epochs during the training of the neural network.

Now we apply our proposed method to the exit problem of the system (2.1). We randomly and uniformly select 10000 points in the domain [1,7]×[0,2]1702[1,7]\times[0,2][ 1 , 7 ] × [ 0 , 2 ], within which N=8018𝑁8018N=8018italic_N = 8018 collocation points located on the right-hand side of the stable manifold of US are used to train the neural network. As depicted in Fig. 4, the loss function is reduced to the magnitude of 105superscript10510^{-5}10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT, indicating good convergence of the algorithm. The quasipotential function and rotational components of the training results are exhibited in Fig. 5.

Refer to caption
(a) Learned Vθ(x)subscript𝑉𝜃𝑥V_{\theta}(x)italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x )
Refer to caption
(b) Learned l1θ(x)subscript𝑙1𝜃𝑥l_{1\theta}(x)italic_l start_POSTSUBSCRIPT 1 italic_θ end_POSTSUBSCRIPT ( italic_x )
Refer to caption
(c) Learned l2θ(x)subscript𝑙2𝜃𝑥l_{2\theta}(x)italic_l start_POSTSUBSCRIPT 2 italic_θ end_POSTSUBSCRIPT ( italic_x )
Figure 5: The quasipotential function Vθ(x)subscript𝑉𝜃𝑥V_{\theta}(x)italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x ) and the two rotational components l1θ(x)subscript𝑙1𝜃𝑥l_{1\theta}(x)italic_l start_POSTSUBSCRIPT 1 italic_θ end_POSTSUBSCRIPT ( italic_x ), l2θ(x)subscript𝑙2𝜃𝑥l_{2\theta}(x)italic_l start_POSTSUBSCRIPT 2 italic_θ end_POSTSUBSCRIPT ( italic_x ) are learned by machine learning.

We first consider the case of a non-characteristic boundary. We take x1=3subscript𝑥13x_{1}=3italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 3 as the boundary, as it has practical ecological significance. On the one hand, this boundary is located on the right side of the natural boundary, i.e., the stable manifold of the saddle point, and thus it can be used for early warning of vegetation degradation, leaving sufficient time for manual intervention. On the other hand, this boundary is relatively simple, only requiring measurement of vegetation biomass.

By employing the gradient descent method, we can locate the point x=(3,1.0632)superscript𝑥31.0632x^{\ast}=(3,1.0632)italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = ( 3 , 1.0632 ) with the minimal quasipotential on the boundary x1=3subscript𝑥13x_{1}=3italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 3. Starting from xsuperscript𝑥x^{\ast}italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT, we can use the inverse time integration to determine the most probable path connecting xsuperscript𝑥x^{\ast}italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT and SN2, as demonstrated by the red curve in Fig. 6. The blue dashed line represents the path obtained by the shooting method BMLSM . The consistency of these results with those obtained through machine learning is evident.

Refer to caption
Figure 6: Along with the identification of the point x=(3,1.0632)superscript𝑥31.0632x^{\ast}=(3,1.0632)italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = ( 3 , 1.0632 ) as the location of minimal quasipotential, the results of the most probable path connecting xsuperscript𝑥x^{\ast}italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT and SN2 are consistent with the outcomes of the machine learning.

Denote H¯=2Vθ(SN2)¯𝐻superscript2subscript𝑉𝜃SN2\bar{H}=\nabla^{2}V_{\theta}(\text{SN2})over¯ start_ARG italic_H end_ARG = ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( SN2 ). We expand the Hamilton-Jacobi equation near the fixed point x¯=SN2¯𝑥SN2\bar{x}=\text{SN2}over¯ start_ARG italic_x end_ARG = SN2 to obtain the algebraic Riccati equation

H¯1Q¯+Q¯H¯1=a¯,superscript¯𝐻1superscript¯𝑄top¯𝑄superscript¯𝐻1¯𝑎\bar{H}^{-1}\bar{Q}^{\top}+\bar{Q}\bar{H}^{-1}=\bar{a},over¯ start_ARG italic_H end_ARG start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT over¯ start_ARG italic_Q end_ARG start_POSTSUPERSCRIPT ⊤ end_POSTSUPERSCRIPT + over¯ start_ARG italic_Q end_ARG over¯ start_ARG italic_H end_ARG start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT = over¯ start_ARG italic_a end_ARG ,

where

Q¯=b(x¯),a¯=a(x¯).formulae-sequence¯𝑄𝑏¯𝑥¯𝑎𝑎¯𝑥\bar{Q}=-\nabla b(\bar{x}),\quad\bar{a}=a(\bar{x}).over¯ start_ARG italic_Q end_ARG = - ∇ italic_b ( over¯ start_ARG italic_x end_ARG ) , over¯ start_ARG italic_a end_ARG = italic_a ( over¯ start_ARG italic_x end_ARG ) .

Therefore, we can get

H¯=(0.05430.06080.06082.9133).¯𝐻0.05430.06080.06082.9133\bar{H}=\left(\begin{array}[]{cc}0.0543&0.0608\\ 0.0608&2.9133\\ \end{array}\right).over¯ start_ARG italic_H end_ARG = ( start_ARRAY start_ROW start_CELL 0.0543 end_CELL start_CELL 0.0608 end_CELL end_ROW start_ROW start_CELL 0.0608 end_CELL start_CELL 2.9133 end_CELL end_ROW end_ARRAY ) .

Then det(H¯)=0.1546¯𝐻0.1546\det(\bar{H})=0.1546roman_det ( over¯ start_ARG italic_H end_ARG ) = 0.1546. Besides,

μ=(12σ12(x1)4Vθx1(x)+l1θ(x))=0.022,det(h)=2Vθx22(x)=2.2602,Vθ(x)=0.0691.superscript𝜇12superscriptsubscript𝜎12superscriptsuperscriptsubscript𝑥14subscript𝑉𝜃subscript𝑥1superscript𝑥subscript𝑙1𝜃superscript𝑥0.022formulae-sequencesuperscriptsuperscript2subscript𝑉𝜃superscriptsubscript𝑥22superscript𝑥2.2602subscript𝑉𝜃superscript𝑥0.0691\begin{array}[]{l}\mu^{\ast}=-\big{(}\frac{1}{2}\sigma_{1}^{2}(x_{1}^{\ast})^{% 4}\frac{\partial V_{\theta}}{\partial x_{1}}(x^{\ast})+l_{1\theta}(x^{\ast})% \big{)}=0.022,\\ \det(h^{\ast})=\frac{\partial^{2}V_{\theta}}{\partial x_{2}^{2}}(x^{\ast})=2.2% 602,\quad V_{\theta}(x^{\ast})=0.0691.\end{array}start_ARRAY start_ROW start_CELL italic_μ start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = - ( divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT divide start_ARG ∂ italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG ( italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) + italic_l start_POSTSUBSCRIPT 1 italic_θ end_POSTSUBSCRIPT ( italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) ) = 0.022 , end_CELL end_ROW start_ROW start_CELL roman_det ( italic_h start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) = divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ( italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) = 2.2602 , italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) = 0.0691 . end_CELL end_ROW end_ARRAY

Due to the asymptotic expression of mean first exit time, we can compute the functional relation between mean first exit time and the noise intensity ε𝜀\varepsilonitalic_ε, as illustrated in Fig. 7. The Monte Carlo simulations confirm the machine learning results, wherein the small error mainly stems from the asymptotic expression itself.

Refer to caption
Figure 7: The machine learning result, represented by the blue curve, is in agreement with the Monte Carlo simulation, indicated by the red stars. This agreement underscores the validity of both methods and their consistency in predicting the relationship between mean first exit time and noise intensity, despite the presence of a small error in the non-characteristic boundary case.

Next, we consider the characteristic boundary, namely, the stable manifold of US. This serves as the natural boundary of the system (2.1). When noise perturbation causes the system to reach this boundary, it tends to move towards SN1 along the unstable manifold of US. Therefore, once the system escapes from this boundary, it is not far from vegetation degradation, and effective measures need to be taken to change the situation.

Refer to caption
Figure 8: The machine learning-computed most probable path from SN2 to US closely agrees with the result obtained from the shooting method.

According to large deviation theory, the saddle US represents the point with minimal quasipotential on the boundary i.e., the exit point. To obtain the most probable path, we integrate the equation x˙=b(x)+a(x)Vθ(x)˙𝑥𝑏𝑥𝑎𝑥subscript𝑉𝜃𝑥\dot{x}=b(x)+a(x)\nabla V_{\theta}(x)over˙ start_ARG italic_x end_ARG = italic_b ( italic_x ) + italic_a ( italic_x ) ∇ italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( italic_x ), commencing from the neighborhood of US with time reversal. As seen in Fig. 8, the most probable path computed by machine learning aligns well with the result obtained via the shooting method. It should be noted that on the most probable path from SN2 to US, the soil moisture level x2subscript𝑥2x_{2}italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT undergoes a process of first decreasing and then increasing. Although the water content at US is much higher than that of the steady state SN2, vegetation inevitably degrades. Hence, in scenarios where vegetation is inadequate, enhancing soil moisture alone is insufficient to halt desertification. Under such circumstances, the effect of planting plants artificially is significantly more beneficial, which is far greater than merely increasing soil moisture.

Denote H=2Vθ(US)superscript𝐻superscript2subscript𝑉𝜃USH^{\ast}=\nabla^{2}V_{\theta}(\text{US})italic_H start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( US ). We solve the equation

(H)1(Q)+Q(H)1=a,superscriptsuperscript𝐻1superscriptsuperscript𝑄topsuperscript𝑄superscriptsuperscript𝐻1superscript𝑎(H^{\ast})^{-1}(Q^{\ast})^{\top}+Q^{\ast}(H^{\ast})^{-1}=a^{\ast},( italic_H start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ( italic_Q start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT ⊤ end_POSTSUPERSCRIPT + italic_Q start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ( italic_H start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT = italic_a start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ,

where

Q=b(US),a*=a(US).formulae-sequencesuperscript𝑄𝑏USsuperscript𝑎𝑎USQ^{\ast}=-\nabla b(\text{US}),\quad a^{*}=a(\text{US}).italic_Q start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = - ∇ italic_b ( US ) , italic_a start_POSTSUPERSCRIPT * end_POSTSUPERSCRIPT = italic_a ( US ) .

We can obtain

H=(0.04460.42280.42281.3305).superscript𝐻0.04460.42280.42281.3305H^{\ast}=\left(\begin{array}[]{cc}-0.0446&0.4228\\ 0.4228&1.3305\\ \end{array}\right).italic_H start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = ( start_ARRAY start_ROW start_CELL - 0.0446 end_CELL start_CELL 0.4228 end_CELL end_ROW start_ROW start_CELL 0.4228 end_CELL start_CELL 1.3305 end_CELL end_ROW end_ARRAY ) .

Then det(H)=0.238superscript𝐻0.238\det(H^{\ast})=-0.238roman_det ( italic_H start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) = - 0.238. Specifically, λ=0.3721superscript𝜆0.3721\lambda^{\ast}=0.3721italic_λ start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = 0.3721, Vθ(US)=0.1643subscript𝑉𝜃US0.1643V_{\theta}(\text{US})=0.1643italic_V start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT ( US ) = 0.1643. Thus the mean first exit time crossing this boundary can be computed via its asymptotic approximation, as plotted in Fig. 9. The machine learning results are also validated by Monte Carlo method. It is observed that the mean first exit time in the characteristic boundary case is much greater than the one in the non-characteristic boundary case, given the same noise intensity. Therefore, these two boundaries can be used for two-level early warning.

Refer to caption
Figure 9: In the characteristic boundary case, the machine learning result, denoted by the blue curve, aligns well with the Monte Carlo simulation marked by red stars. This concordance reinforces the consistency of both methods and their reliability in predicting the correlation between mean first exit time and noise intensity.

Finally, we want to investigate how different noise combinations affect the escape behavior of the system. We consider three cases:

(i)

σ1=0.1subscript𝜎10.1\sigma_{1}=0.1italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 0.1, σ2=1subscript𝜎21\sigma_{2}=1italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 1;

(ii)

σ1=0.08subscript𝜎10.08\sigma_{1}=0.08italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 0.08, σ2=0.1subscript𝜎20.1\sigma_{2}=0.1italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 0.1;

(iii)

σ1=0.1subscript𝜎10.1\sigma_{1}=0.1italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = 0.1, σ2=0.8subscript𝜎20.8\sigma_{2}=0.8italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 0.8.

We take the characteristic boundary as an illustrative example. Fig. 10 illustrates the most probable exit paths in the above three cases. Taking case (i) as a benchmark, we find that the path moves downward as σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT decreases, while the path moves upward when σ2subscript𝜎2\sigma_{2}italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT decreases. This phenomenon is understandable. When σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT decreases, the effect of noise in the x1subscript𝑥1x_{1}italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT direction becomes smaller, so the path tends to be vertical, allowing the noise in the x2subscript𝑥2x_{2}italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT direction to exert a more significant influence and causing the path to bend downward. Conversely, as σ2subscript𝜎2\sigma_{2}italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT decreases, the path becomes more horizontal, causing it to move upward.

Refer to caption
Figure 10: The most probable exit paths from the metastable state SN2 to the saddle point US for the three cases (i)-(iii).

In addition, Fig. 11 depicts the mean first exit time of the three cases (i)-(iii), along with a comparison of Monte Carlo simulations. Obviously, irrespective of the direction in which the noise is reduced, the mean first exit time will increase. Reducing σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT from 0.1 to 0.08 and reducing σ2subscript𝜎2\sigma_{2}italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT from 1 to 0.8 result in the same proportional reduction. However, the growth of mean first exit time after reducing σ2subscript𝜎2\sigma_{2}italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT is significantly greater than the result of reducing σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. Therefore, the noise in the x2subscript𝑥2x_{2}italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT direction plays a dominant role in the escape process of the system.

Refer to caption
Figure 11: The graph compares the mean first exit time among the three cases (i)-(iii), including a side-by-side comparison with Monte Carlo simulations.

5 Conclusion and future perspective

In this paper, we proposed a machine learning method to investigate rare events of stochastic dynamical systems with multiplicative Gaussian noise. We computed the most probable paths, quasipotential, and mean first exit time of a stochastic vegetation-water system in both non-characteristic and characteristic boundaries via the machine learning algorithm. We analyzed the dynamics of the system and explored the feasibility of using the exit phenomenon to establish early warnings for vegetation degradation.

We can innovatively apply machine learning methods to more complex and high-dimensional stochastic dynamical systems driven by Gaussian multiplicative noise. Furthermore, rare events, most probable paths, quasipotential, transition rates, and mean first exit time of the dynamical system with non-Gaussian Lévy noise are also worthy of further study.

Acknowledgement

The authors acknowledge support from the National Natural Science Foundation of China (Grant Nos. 12302035, 62073166, 62221004), the Natural Science Foundation of Jiangsu Province (Grant No. BK20220917), the Key Laboratory of Jiangsu Province, the Shandong Provincial Natural Science Foundation under Grant ZR2021ZD13, and the Project on the Technological Leading Talent Teams Led by Frontiers Science Center for Complex Equipment System Dynamics (FSCCESD220401).

Data availability

Numerical algorithms source code associated with this article can be found, in the online version, at https://github.com/liyangnuaa/rare-events-in-stochastic-vegetation-system.

References

  • (1) M. Qiao, S. Yuan, Analysis of a stochastic predator-prey model with prey subject to disease and Lévy noise, Stochastics and Dynamics 19(05) (2019) 1950038.
  • (2) S. Yuan, Y. Li, Z. Zeng, Stochastic bifurcations and tipping phenomena of insect outbreak systems driven by α𝛼\alphaitalic_α-stable Lévy processes, Mathematical Modelling of Natural Phenomena 17 (2022) 34.
  • (3) A. Tesfay, D. Tesfay, S. Yuan, J. Brannan, J. Duan, Stochastic bifurcation in single-species model induced by α𝛼\alphaitalic_α-stable Lévy noise, Journal of Statistical Mechanics: Theory and Experiment 2021(10) (2021) 103403.
  • (4) S.E. Selvan, M.S.P. Subathra, A.H. Christinal, U. Amato, On the benefits of Laplace samples in solving a rare event problem using cross-entropy method, Applied Mathematics and Computation 225 (2013) 843-859.
  • (5) S. Yuan, Z. Wang, Bifurcation and chaotic behavior in stochastic Rosenzweig-MacArthur prey-predator model with non-Gaussian stable Lévy noise, International Journal of Non-Linear Mechanics 150 (2023) 104339.
  • (6) M.I. Freidlin, A.D. Wentzell, Random perturbations of dynamical systems, Springer, Berlin, Germany, 2012.
  • (7) S. Yuan, J. Duan, Action functionals for stochastic differential equations with Lévy noise, Communications on Stochastic Analysis 13(3) (2019) 10.
  • (8) Y. Li, J. Duan, X. Liu, Y. Zhang, Most probable dynamics of stochastic dynamical systems with exponentially light jump fluctuations, Chaos: An Interdisciplinary Journal of Nonlinear Science 30(6) (2020) 063142.
  • (9) D.G. Luchinsky, I.A. Khovanov, S. Berri, R. Mannella, P.V.E. McClintock, Optimal fluctuations and the control of chaos, International Journal of Bifurcation and Chaos 12(3) (2002) 583-604.
  • (10) M.I. Dykman, P.V.E. McClintock, V.N. Smelyanski, N.D. Stein, N.G. Stocks, Optimal paths and the prehistory problem for large fluctuations in noise-driven systems, Physical Review Letters 68(18) (1992) 2718-2721.
  • (11) J.X. Zhou, M.D.S. Aliyu, E. Aurell, S. Huang, Quasi-potential landscape in complex multi-stable systems, Journal of the Royal Society Interface 9(77) (2012) 3539-3553.
  • (12) P. Ao, Potential in stochastic differential equations: novel construction, Journal of Physics A: Mathematical and General 37 (2004) L25-L30.
  • (13) L. Chen, W. Zhu, First passage failure of quasi non-integrable generalized Hamiltonian systems, Archive of Applied Mechanics 80 (2010) 883-893.
  • (14) I. Franovic, K. Todorovic, M. Perc, N. Vasovic, N. Buric, Activation process in excitable systems with multiple noise sources: One and two interacting units, Physical Review E 92 (2015) 062911.
  • (15) T. Naeh, M.M. Kłosek, B.J. Matkowsky, Z. Schuss, A Direct Approach to the Exit Problem, SIAM Journal on Applied Mathematics 50(2) (1990) 595-627.
  • (16) B. Matkowsky, Z. Schuss, C. Tier, Diffusion across characteristic boundaries with critical points, SIAM Journal on Applied Mathematics 43(4) (1983) 673–695.
  • (17) R.S. Maier, D.L. Stein, A scaling theory of bifurcations in the symmetric weak-noise escape problem, Journal of Statistical Physics 83(3) (1996) 291–357.
  • (18) Y. Li, F. Zhao, S. Xu, J. Duan, X. Liu, A deep learning method for computing mean exit time excited by weak Gaussian noise, Nonlinear Dynamics https://doi.org/10.1007/s11071-024-09280-w.
  • (19) M. Heymann, E. Vanden-Eijnden, The geometric minimum action method: A least action principle on the space of curves, Communications on Pure and Applied Mathematics 61(8) (2008) 1052-1117.
  • (20) M.K. Cameron, Finding the quasipotential for nongradient SDEs, Physica D: Nonlinear Phenomena 241(18) (2012) 1532-1550.
  • (21) E. Alpaydin, Machine learning, Mit Press, London, England, 2021.
  • (22) S.N. Steinmann, Q. Wang, Z.W. Seh, How machine learning can accelerate electrocatalysis discovery and optimization, Materials Horizons 10(2) (2023) 393-406.
  • (23) G. Bonaccorso, Machine learning algorithms. Packt Publishing Ltd, 2017.
  • (24) K. Hippalgaonkar, Q. Li, X. Wang, J.W. Fisher III, J. Kirkpatrick, T. Buonassisi, Knowledge-integrated machine learning for materials: lessons from gameplaying and robotics, Nature Reviews Materials 8(4) (2023) 241-260.
  • (25) M. Amini, A. Rahmani, Agricultural databases evaluation with machine learning procedure, Australian Journal of Engineering and Applied Science 8(2023) (2023) 39-50.
  • (26) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, L. Zdeborová, Machine learning and the physical sciences, Reviews of Modern Physics 91(4) (2019) 045002.
  • (27) J.P. Bharadiya, Machine learning and AI in business intelligence: Trends and opportunities, International Journal of Computer 48(1) (2023) 123-134.
  • (28) M.I. Jordan, T.M. Mitchell, Machine learning: Trends, perspectives, and prospects, Science 349(6245) (2015) 255-260.
  • (29) Y. Li, J. Duan, A data-driven approach for discovering stochastic dynamical systems with non-Gaussian Lévy noise, Physica D: Nonlinear Phenomena 417 (2021) 132830
  • (30) X. Chen, L. Yang, J. Duan, G.E. Karniadakis, Solving inverse stochastic problems from discrete particle observations using the Fokker–Planck equation and physics-informed neural networks, SIAM Journal on Scientific Computing 43(3) (2021) B811-B830
  • (31) Y. Xu, H. Zhang, Y. Li, K. Zhou, Q. Liu, J. Kurths, Solving Fokker-Planck equation using deep learning, Chaos: An Interdisciplinary Journal of Nonlinear Science 30(1) (2020) 013133.
  • (32) Y. Li, J. Duan, X. Liu, Machine learning framework for computing the most probable paths of stochastic dynamical systems, Physical Review E 103(1) (2021) 012124.
  • (33) W. Wei, T. Gao, X. Chen, J. Duan, An optimal control method to compute the most likely transition path for stochastic dynamical systems with jumps, Chaos: An Interdisciplinary Journal of Nonlinear Science 32 (2022) 051102.
  • (34) Y. Li, S. Yuan, S. Xu, Controlling mean exit time of stochastic dynamical systems based on quasipotential and machine learning, Communications in Nonlinear Science and Numerical Simulation 126 (2023) 107425.
  • (35) Y. Li, S. Yuan, L. Lu, X. Liu, Computing large deviation prefactors of stochastic dynamical systems based on machine learning, Chinese Physics B in press https://doi.org/10.1088/1674-1056/ad12a8 (2023).
  • (36) H. Zhang, W. Xu, Y. Lei, Y. Qiao, Early warning and basin stability in a stochastic vegetation-water dynamical system, Communications in Nonlinear Science and Numerical Simulation 77 (2019) 258-270.
  • (37) B. Lin, Q. Li, W. Ren, A data driven method for computing quasipotentials, In Mathematical and Scientific Machine Learning PMLR 145 (2022) 652–670
  • (38) F. Bouchet, J. Reygner, Path integral derivation and numerical computation of large deviation prefactors for non-equilibrium dynamics through matrix riccati equations, Journal of Statistical Physics 189 (2) (2022) 21.
  • (39) F. Bouchet, J. Reygner, Generalisation of the Eyring-Kramers transition rate formula to irreversible diffusion processes, Annales Henri Poincaré 17(12) (2016) 3499-3532.
  • (40) S. Beri, R. Mannella, D.G. Luchinsky, A. Silchenko, P.V. McClintock, Solution of the boundary value problem for optimal escape in continuous stochastic systems and maps, Physical Review E 72(3) (2005) 036131.