# Probabilistic Basin of Attraction and Its Estimation Using Two Lyapunov Functions

Authors: Skuli Gudmundsson; Sigurdur Hafstein
Journal: Complexity
 (2018)

---

## Abstract

We study stability for dynamical systems specified by autonomous stochastic differential equations of the formdX(t)=f(X(t))dt+g(X(t))dW(t), with (X(t))t≥0 an Rd-valued Itô process and (W(t))t≥0 an RQ-valued Wiener process, and the functions f:Rd→Rd and g:Rd→Rd×Q are Lipschitz and vanish at the origin, making it an equilibrium for the system. The concept of asymptotic stability in probability of the null solution is well known and implies that solutions started arbitrarily close to the origin remain close and converge to it. The concept therefore pertains exclusively to system properties local to the origin. We wish to address the matter in a more practical manner: Allowing for a (small) probability that solutions escape from the origin, how far away can they then be started? To this end we define a probabilistic version of the basin of attraction, the γ-BOA, with the property that any solution started within it stays close and converges to the origin with probability at least γ. We then develop a method using a local Lyapunov function and a nonlocal one to obtain rigid lower bounds on γ-BOA.


---

## Body

## 1. Introduction

Lyapunov theory of stability [1] for deterministic systems has been very successful as a tool to study qualitative system behavior. For mathematically oriented studies of the Lyapunov theory, compare, for example, [2–4] or [5–7] for a more modern treatment including complete Lyapunov functions. For more application oriented surveys, compare, for example, [8–10]. The centerpiece of the theory is the Lyapunov function, a real function of state-space which is nonincreasing along solution trajectories. A deterministic Lyapunov function gives rigid estimates on the basin of attraction (BOA) of an equilibrium through its sublevel sets. The equilibrium under consideration can be assumed to be the origin without loss of generality and one then often speaks of the stability of the null solution. Its usefulness has sparked the development of methods to find Lyapunov functions. Breaking the system domain up and addressing each part separately can sometimes work in the deterministic case. With 0∈B⊂N⊂A⊂Rd, one starts by generating a local Lyapunov function W on the small neighbourhood N of the origin and then finding a nonlocal Lyapunov function V:A∖B→R+, using entirely different methods suited away from the equilibrium.In this paper (see Section2), we shall develop the theory necessary to enable such a splitting for a stochastic system. Advanced numerical methods taking advantage of the splitting are in the final stages of development by the authors and their collaborators; see [11] for first results on the local part at the equilibrium. An example using analytical solution at the origin and the shooting method further away is presented in Section 3 in this paper.To set the stage, we shall start by discussing the local and nonlocal Lyapunov functions in the deterministic setting. For definition of the notations used, see Section1.1. Consider a deterministic system given by an autonomous ordinary differential equation (ODE) x˙=f(x), where f:Rd→Rd is sufficiently smooth and such that f(0)=0. To characterize asymptotic stability of the null solution it is convenient to resort to the function classes K∞ and KL (defined in Section 1.1). The null solution is called(locally) asymptotically stable (AS) if there exists a KL function β and a neighbourhood A of the origin such that for every x∈A and t≥0 we have ϕ(t,x)≤β(x,t), where t↦ϕ(t,x) is the solution to the ODE started at x at t=0. In particular, A is a subset of the equilibrium’s BOA defined as {x∈Rd:limsupt→∞ϕ(t,x)=0}. If one can even take A=Rd the origin is calledglobally asymptotically stable (GAS). Obviously (local) AS is too weak a condition and GAS is an unnecessary strict condition for most applications. One is interested in AS where the set A is of reasonable size for the problem at hand. A rigorous lower bound on the BOA can be made if a Lyapunov function V for the system is known, that is, a continuously differentiable function V:U→R+ defined on an open neighbourhood U⊂Rd of the origin and such that μ1(x)≤V(x)≤μ2(x) and ∇V(x)·f(x)≤-μ3(x) for some μ1,μ2,μ3∈K∞ and all x∈U. For a Lyapunov function V all sublevel sets V-1([0,c]), c>0, that are compact subsets of U are subsets of the BOA.For the general deterministic setting, constructing a Lyapunov functionV on a reasonably sized set U is usually a very difficult problem. Often, however, the construction of a Lyapunov function on a small neighbourhood N of the origin is simple. Indeed, if the spectrum of the Jacobian J≔Df(0) of f at the origin does not have any purely imaginary points, then the origin is AS iff the equation J⊤P+PJ=-Q has a positive definite solution P∈Rd×d for a positive definite Q∈Rd×d. Indeed, then VP(x)=x⊤Px is a Lyapunov function for the ODE x˙=Jx on Rd and on some neighbourhood N of the origin it is a Lyapunov function for the nonlinear system x˙=f(x). The size of N, however, depends strongly on the nonlinearities of f; compare, for example, [12, p. 659] for an explicit estimate. By using this or a different method to compute a local Lyapunov function for the system x˙=f(x), one can compute a set B⊂N that is a subset of the equilibrium’s BOA. To assert that a larger set A⊃B is in the BOA, it now suffices to show that ϕ(t,x)≤β(x,t) whenever x,ϕ(t,x)∈A∖B, because the local Lyapunov function VP and the semigroup property ϕ(t,ϕ(s,x))=ϕ(t+s,x) of the flow deliver the missing part. For many methods to numerically compute Lyapunov functions for the system x˙=f(x) on a reasonably sized U⊂Rd it is either advantageous or even necessary to use the fact that a (small) set B is already known to be in the BOA; compare, for example, [12–18] or [19, Section 2.11] and the references therein for an overview.Now, consider a set of stochastic differential equations (SDEs) of the form:dX(t)=f(X(t))dt+g(X(t))dW(t), f:Rd→Rd and g:Rd→Rd×Q, with equilibrium at the origin; that is, f(0)=0 and g(0)=0. As in the deterministic case, a Lyapunov function is a function V defined on some neighbourhood U⊂Rd of the origin, such that μ1(x)≤V(x)≤μ2(x) for some μ1,μ2∈K∞. The condition for deterministic systems that V is decreasing along solution trajectories, that is, ∇V(x)·f(x)≤-α(x) for an α∈K∞, is replaced by the condition that the expectation of the process t↦V(Xy(t)) is a decreasing function of time. Here t↦Xy(t) is the stochastic process started at y at time t=0. By Itô’s formula(1)VXyt=Vy+∫0tLVXysds+∫0t∇VXys·gXysdWs,with (2)LVx≔∇Vx·fx+12∑i,j=1dgxgx⊤ij∂2V∂xi∂xjx,where [g(x)g(x)⊤]ij denotes the (i,j)-th element of the matrix g(x)g(x)⊤. The expectation of the second integral in the Itô formula is zero because the process t↦∇V(Xy(t))·g(Xy(t)) is adapted and therefore the decrease of the expectation, that is, E{V(Xy(t))}<V(y), is implied by LV(x)≤-α(x). In order to find such a Lyapunov function for a stochastic system, one could go about as in the deterministic case and solve the partial differential equation (PDE) LV(x)=-α(x) for an α∈K∞. However, in the deterministic case, the solution V to the PDE ∇V(x)·f(x)≤-α(x) will automatically have local minima at the AS equilibrium but this is not the case for a solution to LV(x)=-α(x). The PDE is markedly different, and the condition LV(x)≤-α(x) alone is no longer sufficient to ensure that the Lyapunov function V has a local minimum at the equilibrium at the origin. It is therefore a natural first approach to try and solve the PDE LV(x)=-α(x) imposing boundary conditions V(0)=0 and V(x)=V0>0 on the boundary ∂U. The problem with this approach is that the PDE cannot be strictly elliptic at the equilibrium at the origin and the standard theory for existence and uniqueness of solutions to PDE does not apply. Making acut-out of some small neighbourhood B containing the origin from the domain A and solving the PDE on A∖B instead can ensure strict ellipticity and hence existence and uniqueness of the solution; compare, for example, [20, Theorem 6.14], and better yet, solving LV(x)=0 so that solutions satisfy both a minimum and a maximum principle (of elliptic PDE) will guarantee that the minima and the maxima of V are taken on the boundary and with V=0 on ∂B and V=1 on ∂A we automatically have 0≤V≤1 on A∖B.Even if the PDE is not strictly elliptic, as will happen ifQ<d, the condition LV(x)<0 instead of LV(x)=0 means that one has more room to actually compute numerically a true Lyapunov function, rather than a numerical approximation of a function we know is a Lyapunov function, as done in [21] where solutions to the PDE are interpreted in the viscosity sense. This will be discussed briefly in Section 4.Despite these advantages, the cut-out of course removes the equilibrium itself from the domain of the Lyapunov function and one must take care off that the null solution is actually stable. In Section2, we describe the details of how this can be achieved by using a different (local) Lyapunov function on a superset of the cut-out set B. In Section 3, we show an example of our approach and in Section 4, we draw a roadmap of how the results in Section 2 can be used to develop a general method to give rigid estimates of γ-BOA of reasonable size for practical applications.


### 1.1. Notation

We denote a (column)-vectorx=(x1,x2,…,xd)⊤ in Rd in boldface and its transpose by x⊤. We also write the matrix-valued function g:Rd→Rd×Q in boldface. We denote the Euclidian norm of a vector x by x and for a matrix · denotes the induced matrix norm. We denote the scalar-product of two vectors x,y∈Rd by x·y. We write sets A⊂Rd in calligraphic and denote by A∘ its interior, by A¯ its closure, by ∂A its boundary, and by AC its complement. We define R+≔[0,∞) and denote by N the set of the natural numbers larger than zero.A functionα:R+→R+ is said to be of class K∞ if it is continuous, monotonically increasing, α(0)=0, and limx→∞α(x)=∞. A function γ:R+→R+ is said to be of class L if it is continuous, monotonically decreasing, and limx→∞γ(x)=0. A function β:R+2→R+ is said to be of class KL if r↦β(r,s) is of class K∞ for every fixed s∈R+ and s↦β(r,s) is of class L for every fixed r∈R+.We writeP and E for probability and expectation, respectively. The underlying probability spaces should always be clear from the context. Conditional probabilities and expectations are denoted by P(·∣·) and E(·∣·), respectively. We denote the characteristic function of a set A by 1A. The abbreviation a.s. stands foralmost surely, that is, with probability one, and =a.s. means equal a.s. By Lp(R+,Rd) and Lp(R+,Rd×Q), we denote the set of all Rd-valued and Rd×Q-valued adapted processes (Y(t))t∈R+ such that ∫0TYspds<∞ a.s. for all T∈R+. Finally, a∧b≔min⁡(a,b).


## 1.1. Notation

We denote a (column)-vectorx=(x1,x2,…,xd)⊤ in Rd in boldface and its transpose by x⊤. We also write the matrix-valued function g:Rd→Rd×Q in boldface. We denote the Euclidian norm of a vector x by x and for a matrix · denotes the induced matrix norm. We denote the scalar-product of two vectors x,y∈Rd by x·y. We write sets A⊂Rd in calligraphic and denote by A∘ its interior, by A¯ its closure, by ∂A its boundary, and by AC its complement. We define R+≔[0,∞) and denote by N the set of the natural numbers larger than zero.A functionα:R+→R+ is said to be of class K∞ if it is continuous, monotonically increasing, α(0)=0, and limx→∞α(x)=∞. A function γ:R+→R+ is said to be of class L if it is continuous, monotonically decreasing, and limx→∞γ(x)=0. A function β:R+2→R+ is said to be of class KL if r↦β(r,s) is of class K∞ for every fixed s∈R+ and s↦β(r,s) is of class L for every fixed r∈R+.We writeP and E for probability and expectation, respectively. The underlying probability spaces should always be clear from the context. Conditional probabilities and expectations are denoted by P(·∣·) and E(·∣·), respectively. We denote the characteristic function of a set A by 1A. The abbreviation a.s. stands foralmost surely, that is, with probability one, and =a.s. means equal a.s. By Lp(R+,Rd) and Lp(R+,Rd×Q), we denote the set of all Rd-valued and Rd×Q-valued adapted processes (Y(t))t∈R+ such that ∫0TYspds<∞ a.s. for all T∈R+. Finally, a∧b≔min⁡(a,b).


## 2. Stochastic Differential Equations and Stability

Stochastic differential equations (SDEs) often come up when modeling a physical system described by an ODE and perturbing this system with noise, corresponding to either uncertainty, measurement error, or unknown complications in the system. The stochastic integral as presented by Itô is in some sense the natural way to approach the situation when the underlying deterministic system responds causally to the noise. This is the case for example when describing financial markets when it is instrumental that the noise has no autocorrelation. In cases when the noise is an effective model of some unknown and/or complicated dynamic subsystem, then it can be more natural to represent this via the so-calledStratonovich stochastic integral. In this paper, we shall restrict our attention to the Itô stochastic integral since a Stratonovich approach is easily represented by a slightly modified Itô system. We study the stability of the origin of the SDE of Itô type(3)dXt=fXtdt+gXtdWt,where f:Rd→Rd and g:Rd→Rd×Q are (globally) Lipschitz continuous functions; that is, there exists a K>0 such that(4)fx-fy+gx-gy≤Kx-y∀x,y∈Rd,and such that f(0)=0 and g(0)=0.We considerstrong solutions to (3), that is, given a filtered probability space (Ω,F,(Ft)t≥0,P) satisfying theusual conditions (F complete, (Ft)t≥0 right continuous, F0 containing all F null sets), a RQ-valued Brownian motion W defined on [t0,∞), and an Ft0 measurable initial distribution Z fulfilling EZ2<∞, an Rd-valued stochastic process X defined on [t0,∞) is called a (unique) solution to (3) if it has the following properties (cf., e.g., [22, Section 2.3], [23, Section 21]):Existence (i)
t ↦ X ( t ) is continuous a.s. and Ft-adapted. (ii)
The processest↦f(X(t)) and t↦g(X(t)) belong to L1(R+,Rd) and L2(R+,Rd×Q), respectively. (iii)
For everyt≥t0, the equation(5)Xt=Z+∫t0tfXsds+∫t0tgXsdWs holds a.s. The second integral is interpreted in the Itô sense.Uniqueness. Any solution Y to (3) is indistinguishable from X; that is, for every t≥0, we have X(t)=a.s.Y(t).For deterministic initial solutions, that is,X(t0)=Z=x∈Rd a.s., we write Xx for the solution.The following theorem provides the standard solution theory, the backdrop to our study of stability later.Theorem 1.
The SDE (3) with initial distribution Z has a unique solution XZ. Further, t↦XZ(t) is a strong Markov process; that is, for every bounded Borel-measurable function φ:Rd→R, any Ft-stopping-time τ<∞ a.s., and any t≥0, we have (6)EφXZτ+t∣Fτ=EφXZτ+t∣XZτ.For a proof of the last theorem, compare, for example, [22, Section 2.3 and Theorem 2.9.3]. Note that since we have f(0)=0 and g(0)=0 the Lipschitz condition implies the so-calledlinear growth condition.A plethora of concepts are in use concerning the stability of SDEs. Here we will be concerned with the so-calledasymptotic stability in probability of the zero solution [24, (5.15)], also referred to asstochastic asymptotic stability [22, Definition 4.2.1]. For a more detailed discussion of the stability of SDEs, see the books by Khasminskii [24] or Mao [22]. For completeness, we recall the definitions of stability in probability and asymptotic stability in probability.Definition 2 (stability in probability).
The null solutionX(t)=a.s.0 to the SDE (3) is said to be stable in probability (SiP) if for every r>0 and every 0<ε<1 there exists a δ>0 such that (7)x≤δimpliesPsupt≥0Xxt≤r≥1-ε.Definition 3 (asymptotic stability in probability (ASiP)).
The null solutionX(t)=a.s.0 to the SDE (3) is said to be asymptotically stable in probability (ASiP) if it is SiP and in addition for every 0<ε<1 there exists a δ>0 such that (8)x≤δimpliesPlimt→∞Xxt=0≥1-ε.Our definitions of SiP and ASiP are equivalent to the more common(9)lim∥x∥→0Psupt>0Xxt≤r=1∀r>0forSiPandadditionallylim∥x∥→0Plimsupt→∞Xxt=0=1forASiP.The reason for our formulation is that we want to look at a more practical concept related to such stability, namely, stochastic analog of the BOA in the stability theory for deterministic systems. We therefore back from the limit x→0 and consider the following instead: Given someconfidence  0<γ≤1 how far from the origin the sample paths can start and still approach the equilibrium as t→∞ with probability greater than or equal to γ. This is the motivation for the next definition.Definition 4 (γ-basin of attraction (γ-BOA)).
Consider system (3) and let 0<γ≤1. One refers to the set (10)x∈Rd:Plimt→∞Xxt=0≥γγ-BOAas the γ-basin of attraction or short γ-BOA.Note that a 1-BOA corresponds to the usual BOA for deterministic systems.For the SDE (3), the associatedgenerator is given by(11)LVx≔∇Vx·fx+12∑i,j=1dgxgx⊤ij∂2V∂xi∂xjxfor some appropriately differentiable V:U→R with U⊂Rd. Notice that this is just the drift term in the expression for the stochastic differential of the process t↦V(X(t)); compare (1). As discussed in the Introduction, the generator for a stochastic system corresponds to the orbital derivative of a deterministic system.Definition 5 (local Lyapunov function).
Consider system (3). A function W∈C(N)∩C2(N∖{0}), where 0∈N⊂Rd is a domain, is called a (local) Lyapunov function for system (3) if there are functions μ1,μ2,μ3∈K∞, such that W fulfills the properties: (i)
μ1(x)≤W(x)≤μ2(x) for all x∈N.(ii)
LW(x)≤-μ3(x) for all x∈N∖{0}.Remark 6.
It is of vital importance thatW is not necessarily differentiable at the equilibrium, because otherwise a large number of systems with an ASiP null solution do not possess a Lyapunov function; compare [24, Remark 5.5].The following theorem provides the first centerpiece of Lyapunov stability theory for our application; compare [24, Theorem 5.5 and Corollary 5.1].Theorem 7.
If there exists a local Lyapunov functionW:N→R+ as in Definition 5 for system (3), then the null solution is ASiP. Further, let Wmax>0 and assume that W-1([0,Wmax]) is a compact subset of N. Then, for every 0<β<1, the set W-1([0,βWmax]) is a subset of the (1-β)-BOA of the origin.Proof.
That the null solution is ASiP follows directly by [24, Theorem 5.5 and Corollary 5.1].
To prove the second assertion fix0<β<1 and an arbitrary x∈W-1([0,βWmax]). We must show that (12)Plimt→∞Xxt=0≥1-β.Let ρ≔inf{t>0:Xx(t)∉W-1([0,Wmax])} be the first exit-time of the process t↦Xx(t) from W-1([0,Wmax]). Then ρ is a stopping-time and using identical argumentation as in the proof of [24, Theorem 5.5], we note that a.s. on {ρ=∞} we have limt→∞Xx(t)=0. It thus suffices to show that P{ρ=∞}≥1-β.
By [24, Lemma 5.4] the stochastic process t↦W(Xx(ρ∧t)) is a supermartingale; compare also Remark 8. Therefore (13)βWmax≥Wx≥EWXxρ∧t=E1ρ≥tWXxt+1ρ<tWXxρ≥E1ρ<tWXxρ=Pρ<tWmaxfor any t>0. Since {ρ<∞}=⋃n∈N{ρ<n}, we can take limits to get β≥P{ρ<∞}. Hence (14)Pρ=∞=1-Pρ<∞≥1-βwhich proves the theorem.Remark 8.
If the conditionLW(x)<0 is violated, even at only one point x∈N∖{0}, one cannot conclude that t↦W(Xx(t∧ρ)) is a supermartingale; compare [24, p. 149]. The reason why LW(x)≤0 does not have to hold true at the equilibrium at the origin is because it is so-calledinaccessible point of the SDE, compare again [24, p. 149], meaning that (15)P∃t>0:Xxt=0=0foreveryx≠0.Instead of trying to compute a local Lyapunov function for system (3) on a large domain N, we suggest to compute it on a small N, for example, by using linearization, and then to compute a nonlocal Lyapunov function, defined below, to compute a reasonably sized rigid estimate on the equilibrium’s γ-BOA. The advantages of this procedure were explained in the Introduction and are addressed again in Section 4, where we give a roadmap of our approach.Definition 9 (nonlocal Lyapunov function).
LetA,B⊂Rd, B⊂A∘, be simply connected compact neighbourhoods of the origin with C2 boundaries and set U≔A∖B∘. A function V∈C2(U) for system (3) such that (1)
0≤V(x)≤1 for all x∈U, V-1(0)=∂B, V-1(1)=∂A, and either(2a)
LV(x)<0 for all x∈U or(2b)
LV(x)≤0 and g(x)g(x)⊤ is positive definite for all x∈U,
is called a nonlocal Lyapunov function for system (3). We refer to ∂A as the outer boundary of U and ∂B as the inner boundary of U.Lemma 10.
Assume system (3) has a nonlocal Lyapunov function V∈C2(U) as in Definition 9 and let x∈U=A∖B∘. Let τ≔inf{t>0:Xx(t)∉U} be the exit-time from U. Then(16)Pτ<∞=1,PXxτ∈∂B≥1-Vx.Proof.
We first assume thatV fulfills (2a) in Definition 9 and set h≔-maxx∈ULV(x) and note that since the solution process remains in U for all s∈[0,t∧τ], we get by [24, Lemma 3.2] and the nonnegativity of V that (17)0≤EVXxτ∧t=Vx+E∫0τ∧tLVXxsds≤Vx-hEτ∧t≤Vx.Now (18)Eτ∧t=E1t>ττ+1t≤τt≥Pt≤τt,and then (19)0≤hPt≤τt≤hEτ∧t≤Vx∀t>0.Taking the limit t→∞ delivers P{∞≤τ}=0 and then P{τ<∞}=1.
If the condition (2b) in Definition9 is fulfilled (instead of condition (2a)), then, by [20, Theorem 6.14], it follows that the boundary value problem (20)LUx=-1∀x∈U,Ux=0∀x∈∂B,Ux=1∀x∈∂Ahas a solution U∈C2(U). Just as above with U substituted for V, it follows that P{τ<∞}=1.
For the second proposition of the lemma define the stopping-timesτ∂A≔inf{t>0:Xx(t)∈∂A} and τ∂B≔inf{t>0:Xx(t)∈∂B} and observe that τ=τ∂A∧τ∂B a.s. because t↦Xx(t) is continuous a.s. Now (21)VXxτ=1Xxτ∈∂B·VXxτB+1Xxτ∈∂A·VXxτAand then (22)EVXxτ=PXxτ∈∂B·0+PXxτ∈∂A·1=PXxτ∈∂A;that is, we have (23)1-PXxτ∈∂B=PXxτ∈∂A=EVXxτ≤Vxbecause P{τ<∞}=1, which completes the proof.We now show how to combine a local Lyapunov function and a nonlocal one to get a rigorous lower bound onγ-BOA. For the schematic representation of the interrelation between the (sub)level sets in Theorem 11 see Figure 1.Figure 1
A schematic figure of the (sub)level sets in Theorem11. Outward from the equilibrium (thick dot): ∂B=V-1(0)⊂W-1([0,βWmax])⊂W-1([0,Wmax]) and then V-1([0,λ])⊂V-1([0,α])⊂V-1([0,1])=∂A. In particular, W-1(Wmax)⊂V-1([0,λ]). The local Lyapunov function W asserts that a solution starting in the innermost area leaves the red area with probability no higher than β and the nonlocal Lyapunov function V asserts that a solution leaving the red area goes to the innermost area with probability no less than 1-λ.Theorem 11.
Consider system (3) and assume there exist a local Lyapunov function W:N→R+ as in Definition 5 and a nonlocal Lyapunov function V:U→R+ as in Definition 9 for the system. Assume further that the constants Wmax>0 and 0<β<1 from Theorem 7 and the constants 0<λ<α<1 and the set B from Lemma 10 are such that (24)W-1Wmax⊂V-10,λ,∂B=V-10⊂W-10,βWmax.Then the set V-1([0,α])∪B is a subset of the γ-BOA of the origin, where γ≔(1-α)(1-β)/(1-β1-λ).Proof.
Letx∈B∪V-1([0,α]) be arbitrary but fixed. We must show that (25)Plimt→∞Xxt=0≥1-α1-β1-β1-λ.Define the stopping-times τ∂A(y)≔inf{t>0:Xy(t)∈∂A}, τ∂B(y)≔inf{t>0:Xy(t)∈∂B}, and τ(y)≔τ∂A(y)∧τ∂B(y). In the following, we suppress the argument of stopping-times (as usual) if it can be seen from the context. Define for every y∈V-1([0,α])∪B the event(26)Ey≔∃t>0:Xyt∈∂A.Note that since ∂A is closed and the sample paths are continuous a.s., we have E(y)=a.s.{τ∂A<∞}.
Let us first assumex∈V-1([0,α]); we deal with the case x∈B below. Then (27)Ex=Xxτ∈∂A∪Xxτ∈∂B∩∃t>τ:Xxt∈∂A.Because the events {Xx(τ)∈∂A} and {Xx(τ)∈∂B} are disjoint we have(28)PEx=PXxτ∈∂A+PXxτ∈∂B∩∃t>τ:Xxt∈∂A=PXxτ∈∂A+P∃t>τ:Xxt∈∂A∣Xxτ∈∂BPXxτ∈∂B.The strong Markov property of t↦Xx(t) (cf. Theorem 1) delivers (cf., e.g., [25, Theorem 7.2.4] or [24, Theorem 2.9.5]) (29)P∃t>0:XXxτt∈∂A∣Xxτ∈∂B=P∃t>τ:Xxt∈∂A∣Xxτ∈∂B.Define (30)P≔supy∈∂BPEy.Now ν(G)≔P{Xx(τ)∈G}/P{Xx(τ)∈∂B}, G⊂∂B a Borel set with respect to the relative topology of ∂B, defines a probability measure on ∂B and, for example, by [23, Lemma 8.7], we get(31)P∃t>0:XXxτt∈∂A∣Xxτ∈∂B=∫∂BPEyνdy≤P∫∂Bνdy=P.Further, because P{τ<∞}=1 by Lemma 10, we have P{Xx(τ)∈∂B}=1-P{Xx(τ)∈∂A} and we get from (28) the estimate (32)PEx≤PXxτ∈∂A+1-PXxτ∈∂APwhich together with P{Xx(τ)∈∂A}≤α and the fact that P≤1 implies(33)PEx≤α+1-αP.Similarly, we get for every y∈V-1([0,λ]) that(34)PEy≤λ+1-λP.Now let z∈∂B be arbitrary and define the stopping-time ρ≔inf{t>0:Xz(t)∈W-1(Wmax)}. Then (35)PEz=Pρ<∞∩∃t>ρ:Xzt∈∂A=P∃t>ρ:Xzt∈∂A∣ρ<∞Pρ<∞.The assumption ∂B⊂W-1([0,βWmax]) delivers by Lemma 10 that P{ρ<∞}≤β. By using the assumption W-1(Wmax)⊂V-1([0,λ]) we can show just as above that (36)P∃t>ρ:Xzt∈∂A∣ρ<∞≤λ+1-λP.Thus(37)PEz≤λ+1-λPβand especially(38)P≤λ+1-λPβ,i.e.P≤βλ1-β1-λ.Plugging this estimate for P into (33) delivers(39)1-PEx≥1-α+1-α-P≥1-α1-β1-β1-λfor every x∈V-1([0,α]).
That estimate (39) also holds true if x∈B⊂W-1([0,βWmax]) can be seen from the fact that estimate (37) is clearly valid for any z∈W-1([0,βWmax]). Combining (37) and (38) and noting that 0<α we get (40)1-PEx≥1-β1-β1-λ>1-α1-β1-β1-λ.
We finish the proof by showing thatE(x)C=a.s.{τ∂A=∞} is equal to {limt→∞Xx(t)=0} a.s. Note first that t↦Xx(t) cannot stay in A∖B for an infinitely long contiguous time interval a.s. by Lemma 10. Second, as in the proof of Theorem 7, we have limt→∞Xx(t)=0 a.s. if t↦Xx(t) stays in W-1([0,Wmax]) for an infinitely long contiguous time interval. Thus {τ∂A=∞}∩{limt→∞Xxt≠0} implies the following event a.s.: There exists a sequence of increasing times(tn)n∈N, tn→∞, such that Xx(t2n-1)∈∂B and Xx(t2n)∈∂W-1(Wmax) for all n∈N.
Because∂B⊂W-1([0,βWmax]) and W-1(Wmax)⊂V-1([0,λ]) and by what we have showed above the probability of this event is no larger than β·(1-λ)·β·(1-λ)⋯=0, which completes the proof.Remark 12.
Note that in the proof of Theorem11 the conditions on f and g in (3) need only to hold within A, because we stop our solutions at the boundary ∂A, and therefore f and g can be altered outside of A without changing the conclusions of the theorem. It thus suffices that f and g fulfill the Lipschitz condition (4) on A, which is, for example, always the case if f and g are locally Lipschitz.


## 3. Worked Out Example

We consider an example to show how our approach works. Because we do not want to consider appropriate methods to solve the relevant PDEs in this paper we consider a 1-dimensional example, but a highly nonlinear one. The system is given by the SDE(41)dXt=a·sin⁡Xtdt+b·Xt1+Xt2dW~t,so that our coefficient functions in (3) are given by fx=asin⁡x and g(x)=bx/(1+x2) and we have only a single Wiener process, denoted by W~ because we reserve W for the local Lyapunov function, driving the noise.We first compute a local Lyapunov function for (41) as in Definition 5 using the ansatz W(x)=|x|p with 0<p<1. Direct calculations with x≠0 reveal (42)W′x=pxpx,W′′x=-p1-pxp-2,so (43)LWx=fxW′x+12g2xW′′x=asin⁡x·pxpx+b22x21+x22-p1-pxp-2=-pxpb21-p21+x22-asin⁡xx.We fix a=1, b=3, and p=1/2. This corresponds to the deterministic system x˙=sin⁡x with an unstable null solution. The noise, however, stabilizes the null solution, which can be seen from LW(x)<0 for all x∈[-2-1/2,2-1/2]∖{0}≈[-0.70711,0.70711]∖{0} because then (44)941+x22≥1>sin⁡xx.Thus {±2-1/2}=W-1(Wmax) with Wmax≔2-1/4. To compute a rigid lower estimate [-r1-β,r1-β] on the (1-β)-BOA of the equilibrium using W we can, by Theorem 7 and because W is an even function, solve W(r1-β)=βWmax. Thus r1-β=β22-1/2. In the following we use β=0.02 which delivers r0.98=4·10-4·2-1/2≈2.828·10-4.To enlarge our estimate we solve the PDELV(x)=0 (note that g(x)2>0 for x≠0) with the boundary conditions V(10-4)=0 and V(πk/2)=1 with k=5; compare Lemma 10 and Theorem 11. Because the system is 1-dimensional the PDE is also an ODE and can be conveniently solved with the shooting method. With the initial conditions V(10-4)=0 and V′(10-4)=0.00157296, we get V(2.5π)=0.99998 using a very conservative step-size of 10-5. This corresponds to B=[-10-4,10-4] and A=[-2.5π,2.5π]. From the results (cf. Figure 2), we can estimate V-1([0,0.09])⊃[-5.77,5.77]∖B and V-1([0,2·10-4])⊃[-0.75,0.75]∖B. This corresponds to α=0.09 and λ=2·10-4. Note that by extending V(x) canonically for negative x, we have {±10-4}=V-1(0)⊂[-r0.98,r0.98]=W-1([0,βWmax]) and {±2-1/2}=W-1(Wmax)⊂[-0.75,0.75]∖B⊂V-1([0,λ]). Thus Theorem 11 delivers that [-5.77,5.77] is affirmed to be contained in the γ-BOA of the equilibrium at the origin with (45)γ=1-α1-β1-β1-λ=1-0.091-0.021-0.021-2·10-4≈0.91.For comparison we get the rigid bound r1-0.1=0.12·2-1/2=0.00707, that is, [-0.00707,0.00707] on the 0.9-BOA, if we only use the local Lyapunov function W. Our combined method thus delivers a considerably larger set that is affirmed to be contained in the equilibrium’s 0.9-BOA, namely, [-5.77,5.77] instead of [-0.00707,0.00707] from the local Lyapunov function. For different γ the results are comparable. In Figure 2 we show nonlocal Lyapunov functions computed for the system, not only with k=5 but also for k=1,2,3,4 to give a better idea of what is happening.Figure 2
Nonlocal Lyapunov functions (solid blue) for system (41) computed on [10-4,πk/2] with k=1,2,3,4,5. Since the value V(x) is the probability that the solution started at x leaves the domain of the Lyapunov function V at the outer boundary, one obviously gets unnecessary conservative estimates on the γ-BOA if the domain is too small. Extending the domain beyond ca. 7.5 however does not lead to notably larger estimates. The local Lyapunov function W(x)=|x|1/2 (dotted red) is plotted on the interval [0,2-1/2].


## 4. Conclusions and Further Work

In this paper, we have worked out theory in Theorem11 for a method to give rigid estimates on the γ-basin of attraction (γ-BOA) of the null solution as defined in Definition 4 for stochastic systems given by SDEs (3). We showed how our method works for a highly nonlinear 1-dimensional system in Section 3. For higher dimensional systems, one has to consider appropriate methods for our setting to solve PDEs.Our roadmap for further advancing the method is as follows:(1)
Given system (3), we start by finding a local Lyapunov function W:N→R+ as in Definition 5 and an Wmax>0 as large as possible, such that W-1([0,Wmax]) is a compact subset of N.One method for doing this involves linearizing the system (i.e., the functionsf and g) via Taylor expansion around x=0. For a wide class of systems, we expect the linearized system to admit locally a Lyapunov function of the form W(x)=Qxp with 0<p<1 and Q∈Rd×d a positive definite matrix; compare [24, Theorem 5.12]. We have worked with our collaborators on a method using sum-of-squares programming for computing such Lyapunov functions for linear systems and the initial results are promising [11].(2)
Next we determine a nonlocal Lyapunov functionV:A∖B∘→R+ with A bounded but as large as practical. This must be done in such a way as to ensure B⊂W-1([0,βWmax]) and W-1(Wmax)⊂V-1([0,λ]) for appropriate 0<β, λ<1; compare Theorem 11. From these Lyapunov functions we know that for every α, λ<α<1, the set V-1([0,α])∪B is a subset of the γ-BOA of the null solution, where γ≔(1-α)(1-β)/(1-β1-λ).We propose to address this by solving the following PDE boundary value problem using an appropriate numerical method:(46)LV=-honA∖B∘,V=0on∂B,1on∂A.If the PDELV=-h is strictly elliptic, that is, the matrix g(x)g(x)⊤ is positive definite for all x∈A∖B∘, then the PDE possesses a unique solution; compare, for example, [20, Theorem 6.14]. In this case, we can set h=0 and from a rigid error estimate of the numerical approximate solution to the unique solution, we can compute rigid sublevel sets of V and thus give rigid lower bounds on γ-BOA.If the PDELV=-h is not strictly elliptic, then we suggest fixing h>0. In this case, we must be sure that the solution V we compute numerically really fulfills LV(x)<0 for all x∈A∖B∘ and that its sublevel sets are compatible with Theorem 11. That is, the computed function must be a true nonlocal Lyapunov function for the system and not just an approximation!We are optimistic that our suggested method to compute rigid estimates toγ-BOA can be developed into a useful tool for the study of stochastic systems, both in practice and in theory.



---

*Source: 2895658-2018-04-01.xml*