Positively Invariant Sets for ODEs and Numerical Integration

Peter Giesl

1 a

, Sigurdur Hafstein

2 b

and Iman Mehrabinezhad

2 c

Department of Mathematics, University of Sussex, Falmer, BN1 9QH, U.K.

Science Institute, University of Iceland, Dunhagi 3, 107 Reykjav

ık, Iceland

Keywords:

Positively Invariant Sets, Ordinary Differential Equations, Numerical Integration.

Abstract:

We show that for an ordinary differential equation (ODE) with an exponentially stable equilibrium and any

compact subset of its basin of attraction, we can ﬁnd a larger compact set that is positively invariant for both

the dynamics of the system and a numerical method to approximate its solution trajectories. We establish this

for both one-step numerical integrators and multi-step integrators using sufﬁciently small time-steps. Further,

we show how to localize such sets using continuously differentiable Lyapunov-like functions and numerically

computed continuous, piecewise afﬁne (CPA) Lyapunov-like functions.

1 INTRODUCTION

In this paper we study the ordinary differential equa-

tion (ODE)

x = f(x), f ∈C

), s ≥ 1, (1)

and numerical methods to approximate its solution for

given initial values. The solution x(t) to the initial

value problem (1) with x(0) = ξ

ξ is denoted by φ

φ(t,ξ

ξ).

Our ultimate goal is to characterize the long term

behavior of system (1) by computing Lyapunov func-

tions and contraction metrics. As many real-world

systems are modelled by (1), such a characterization

has numerous practical uses in engineering and sci-

ence. In this paper we assume that (1) possesses an

exponentially stable equilibrium at x

∈ R

, i.e. the

Jacobian Df(x

) ∈ R

n×n

of f at x

is Hurwitz (the real

part of all its eigenvalues is negative). We denote the

equilibrium’s basin of attraction by

A(x

) := {x ∈ R

: lim

t→∞

φ(t,x) = x

The Lyapunov stability theory is a generalization of

the concept of dissipative energy in physics. It is

the centerpiece of practical and theoretical stabil-

ity analysis and is treated in various detail in virtu-

ally all textbooks and monographs on ODE systems,

cf. e.g. (Hahn, 1967; Yoshizawa, 1966; Zubov, 1964;

Khalil, 2002; Sastry, 1999; Vidyasagar, 2002).

https://orcid.org/0000-0003-1421-6980

https://orcid.org/0000-0003-0073-2765

https://orcid.org/0000-0002-6346-9901

For a general ODE there is no obvious candi-

date for a Lyapunov function and there are no an-

alytical methods to compute a Lyapunov function.

Hence, various methods for the numerical compu-

tation of Lyapunov functions have been developed.

To name a few, in (Valmorbida and Anderson, 2017;

Vannelli and Vidyasagar, 1985) the computation of

rational Lyapunov functions was studied, in (An-

derson and Papachristodoulou, 2015; Parrilo, 2000)

sum-of-squared (SOS) polynomial Lyapunov func-

tions were computed using semi-deﬁnite optimization

(SOS method), see also (Ratschan and She, 2010;

Kamyar and Peet, 2015) for other approaches using

polynomials, and in (Giesl, 2007) a Zubov type PDE

was numerically solved using radial basis functions

(RBF method). For an overview of more methods see,

e.g., the review paper (Giesl and Hafstein, 2015b).

In (Julian et al., 1999; Marin

osson, 2002) linear pro-

gramming was used to compute continuous and piece-

wise afﬁne (CPA) Lyapunov functions; this approach

is referred to as the CPA method. In the CPA method

the domain, where the Lyapunov function is to be

computed, is triangulated, i.e. subdivided into sim-

plices, and a feasibility problem is derived, such that

its solution can be used to deﬁne a CPA Lyapunov

function for the system. In (Giesl and Hafstein, 2014;

Hafstein, 2004; Hafstein, 2005) it was shown that the

CPA method always succeeds in computing a Lya-

punov function for an ODE with an asymptotically

stable equilibrium, if the simplices in the triangulation

are sufﬁciently small and non-degenerate in a certain

sense.

Giesl, P., Hafstein, S. and Mehrabinezhad, I.

Positively Invariant Sets for ODEs and Numerical Integration.

DOI: 10.5220/0012189700003543

In Proceedings of the 20th International Conference on Informatics in Control, Automation and Robotics (ICINCO 2023) - Volume 1, pages 44-53

ISBN: 978-989-758-670-5; ISSN: 2184-2809

In (Giesl and Hafstein, 2015a), the CPA and the

RBF method were combined to deliver a method that

is as fast as the RBF method and delivers a veriﬁed

Lyapunov function as the CPA method. This was

achieved by solving a system of linear equations in the

RBF method rather than a linear optimization prob-

lem as in the CPA method. At the same time, the

obtained function is veriﬁed to be a Lyapunov func-

tion by checking that it satisﬁes the constraints of the

feasibility problem. Further, the authors proved that

this approach is constructive and one is always able

to compute a veriﬁed Lyapunov function in any com-

pact subset of an exponentially stable equilibrium’s

basin of attraction. A similar approach uses numer-

ical integration of solution trajectories of the ODE

to generate values for the variables of the feasibility

problem of the CPA method and then veriﬁes the con-

straints, see (Bj

ornsson et al., 2014; Bj

ornsson et al.,

2015; Bj

ornsson and Hafstein, 2017; Doban, 2016;

Doban and Lazar, 2016; Hafstein et al., 2014a; Haf-

stein et al., 2014b; Hafstein et al., 2015; Hafstein and

Valfells, 2017; Hafstein and Valfells, 2019; Li et al.,

2015) and also (Gudmundsson and Hafstein, 2015;

Hafstein, 2019a) for more implementation oriented

papers. This technique works well in practice and

in (Giesl and Hafstein, 2023) it is proved that it al-

ways works, assuming that one can ﬁnd a compact set

that is positively invariant for both the ODE and for a

numerical scheme to approximate its solution trajec-

tories. The main contribution of this paper is to prove

the existence of such a positively invariant set.

The results in this paper also make an important

contribution in the context of numerical methods for

contraction metrics, see (Giesl et al., 2023c), where

the results of this paper are used to derive a uniform

error estimate on compact sets. Contraction metrics

are Riemannian metrics such that the corresponding

distance of adjacent solution trajectories decreases ex-

ponentially over time. They have received consider-

able attention in the literature (Lewis, 1949; Lewis,

1951; Demidovi

c, 1961; Krasovski

i, 1963; Borg,

1960; Hartman, 1961; Hartman, 1964; Lohmiller and

Slotine, 1998; Aminzare and Sontag, 2014; Simpson-

Porco and Bullo, 2014; Forni and Sepulchre, 2014;

Giesl, 2015). The analytical computation of a con-

traction metric for an ODE is notoriously difﬁcult,

even more difﬁcult than the computation of a Lya-

punov function, as it requires the computation of a

matrix-valued function. Numerical methods for the

construction of contraction metrics include (Aylward

et al., 2008; Giesl and Hafstein, 2013; Giesl, 2019;

Giesl et al., 2023a), see also the recent review (Giesl

et al., 2023b).

Let us give an overview of the paper: In Section

2 we recall some facts about numerical integration

methods of ODEs and prove a theorem about approxi-

mations of solutions with multi-step methods. In Sec-

tion 3 we establish the existence of positively invari-

ant sets, both for the dynamics of the system (1) and a

numerical integration scheme to approximate the so-

lution trajectories in the basin of attraction of an expo-

nentially stable equilibrium; the main result is Theo-

rem 3.5. Such positively invariant sets are very useful,

in fact necessary, to prove that Lyapunov functions

and contraction metrics can be approximated arbitrar-

ily close on compact subsets of basins of attraction,

using numerical integration with subsequent numeri-

cal quadrature. We prove our results using the fourth-

order Adams-Bashforth (AB4) multi-step scheme ini-

tialized with fourth-order Runge-Kutta (RK4), but we

discuss how the results can be extended to AB-RK

numerical schemes of arbitrary order. Finally we con-

clude the paper in Section 4.

Notation: We deﬁne N

:= {0, 1, 2,...} as the set

of the natural numbers and N

:= N

\{0} as the

set of the positive natural numbers. We denote the

usual p-norms on R

and the corresponding induced

matrix norms by ∥·∥

, 1 ≤ p < ∞. For both vec-

tors in R

and matrices in R

n×n

we write ∥ ·∥

max

for the maximum absolute value norm, i.e. ∥x∥

max

i=1,2,...,n

| for a vector x ∈ R

n×n

and ∥A∥

max

i, j=1,2,...,n

i j

| for a matrix A =



i j



∈ R

n×n

Apart from the usual equivalence estimates for the p-

norms on R

, recall the norm equivalence ∥A∥

max

≤

∥A∥

≤ n∥A∥

max

for a matrix A ∈ R

n×n

and that

∥ · ∥

max

is not sub-multiplicative, but ∥Ab∥

max

≤

n∥A∥

max

∥b∥

max

for b := B ∈ R

n×n

or b := b ∈ R

We denote the closure of a set U ⊂ R

by U and its

boundary by ∂U.

2 NUMERICAL INTEGRATION

METHODS

For computing Lyapunov functions or contraction

metrics for the ODE (1), initial value problems can

be solved numerically for all vertices of a triangu-

lation of a relevant compact subset of R

, as dis-

cussed in the last section. As the number of vertices

can be very high, it is advantageous to use multi-step

methods rather than single-step methods, as these are

considerably faster for the same degree of precision.

Additionally, the examples in (Hafstein, 2019b) indi-

cate that the corrector step in the Adams-Bashforth-

Moulton predictor-corrector methods does not deliver

better results than the Adams-Bashforth method with-

out the corrector step. Therefore we will concentrate

Positively Invariant Sets for ODEs and Numerical Integration

in this paper on the Adams-Bashforth method, and

in order to be concrete, we will concentrate on the

Adams-Bashforth method of order four (AB4) initial-

ized with the usual Runge-Kutta method of the same

order (RK4). However, it is straight-forward to adapt

the proofs to the Adams-Bashforth method of any or-

der initialized with the Runge-Kutta of the same order

and we will elaborate at the appropriate places.

The true solution to the ODE (1) with initial-value

ξ is denoted by t 7→ φ

φ(t,ξ

ξ) and its numerical approx-

imations at t

= hi, i ∈ N

, by

(ξ

ξ) or

, i.e.

(ξ

ξ),

dependent on whether the initial-value ξ

ξ is clear from

the context or not. The time-step h > 0 is a parameter

of the numerical method. Thus, we set

= ξ

ξ and for

i = 0,1,2 we use the RK4 formulas

= hf(

) (2)

= hf(

+ k

/2)

= hf(

+ k

/2)

= hf(

+ k

)

i+1

+ 2k

+ k

For i ≥ 3 we use the AB4 formula and set

i+1

(3)



55f(

) −59f(

i−1

) + 37f(

i−2

) −9f(

i−3

)



Let S ⊂R

be a compact set and f ∈C

). It

is well known that one can ﬁnd a constant C

RK4

such

that for every ξ

ξ ∈ S we have

∥

(ξ

ξ) −φ

φ(h,ξ

ξ)∥

max

≤C

RK4

, (4)

where the constant C

RK4

depends on the derivatives of

f, up to and including the fourth order, in a compact,

convex set

S ⊃ S, fulﬁlling

(ξ

ξ),φ

φ(h,ξ

ξ) ∈

S for all

ξ ∈ S.

The existence of such a set

S is clear from the for-

mulas (2). The reason for this is simply that the er-

ror estimate (2) is derived by using Taylor polynomi-

als for t 7→ φ

φ(t,ξ

ξ) and using the fact that

φ(t,ξ

ξ) =

f(φ

φ(t,ξ

ξ)). From this it is obvious, that for τ

∗

> 0

there exists a constant C such that (4) holds true with

RK4

= C and all 0 < h ≤ τ

∗

and all ξ

ξ ∈ S. This is

the property we need for the assumptions in Theorem

3.5.

For multi-step methods like AB4 the error esti-

mate is usually formulated differently. Namely that

there exists a constant C

AB4

> 0 such that

∥

i+1

(ξ

ξ) −φ

φ((i + 1)h,ξ

ξ)∥

max

≤C

AB4

= φ

φ( jh, ξ

ξ) for j = i,i −1, i −2, i −3 in formula

(3), i.e. if the previous approximations

are exact.

Again, the constant C

AB4

depends on the derivatives

of f, up to and including the fourth order, in a com-

pact, convex set

S ⊃ S, fulﬁlling

i+1

(ξ

ξ),φ

φ( jh, ξ

ξ) ∈

for j = i,i −1, i −2,i −3 for all ξ

ξ ∈ S.

In this form the error estimate is not useful for

our application of computing Lyapunov functions and

contraction metrics. Therefore we prove the follow-

ing theorem.

Theorem 2.1. (Error estimate for AB4-RK4 method)

Consider the system (1) and assume that f ∈

). Then AB4 initialized with RK4 fulﬁlls the

property:

For any compact set S ⊂ R

there exist constants

C, h

∗

> 0 such that for all step-sizes 0 < h ≤ h

∗

have for any i ∈ N

that

∥

i+1

(ξ

ξ) −φ

φ(h,

(ξ

ξ))∥

≤Ch

whenever

(ξ

ξ),

(ξ

ξ),...,

(ξ

ξ) ∈ S.

Proof. Fix a constant 0 < h

≤1 and a compact, con-

vex set S

′

⊃ S such that φ

φ(h,ξ

ξ) ∈ S

′

for all ξ

ξ ∈ S and

i+1

(ξ

ξ) ∈S

′

, whenever

(ξ

ξ) ∈S, for j = 0,1,2,...,i,

i ∈N

, and when using step-size 0 < h ≤h

. Further-

more, let 0 < h

≤ h

be a constant and

S be a com-

pact, convex set such that φ([−3h

,0],S

′

) ⊂

S. Fix an

arbitrary, but constant ξ

ξ ∈ S for the rest of the proof.

For a ﬁxed step-size 0 < h ≤ h

denote φ

(y) :=

φ(ih,y), i ∈ Z, and by

, i ∈ N

, the approximation

to φ

(ξ

ξ) generated by the numerical method, i.e. AB4

initialized with RK4. There exists a constant C

RK4

0, independent of ξ

ξ ∈ S, such that for 0 < h ≤ h

have

∥

i+1

−φ

(

)∥

max

≤C

RK4

, (5)

as long as 0 < h ≤h

and

∈S, j = 0,1,...,i, for i =

0,1,2; this is just the classical local truncation error

estimate for RK4.

For the steps with the AB4 method it is advanta-

geous to deﬁne

AB4

j=i

) :=

[55f(x

) −59f(x

i−1

) + 37f(x

i−2

) −9f(x

i−3

)]

where either x

or x

= φ

(y). Note that x

refers

to the sequence (x

)

, and j in j = i refers to the index

j of the sequence (x

)

. i is the value of the last index

used in the sequence when AB4

j=i

) is computed

from x

i−1

i−2

i−3

. In the case x

we have

AB4

j=i

(

) =

i+1

if i ≥3

ICINCO 2023 - 20th International Conference on Informatics in Control, Automation and Robotics

and in the case x

= φ

(y) there exist a constant

AB4

> 0, such that

∥φ

i+1

(y) −AB4

j=i

(φ

(y))∥

max

≤C

AB4

(6)

as long as 0 < h ≤ h

and φ

(y) ∈

S for

j = i, i −1,i −2,i −3, in particular for all y ∈ S

′

. We

now prove the theorem by induction.

Denote by (I) the proposition:

There exist constants C,C

∗

> 0, such that for

every time-step 0 < h ≤ h

∗

we have for any i ∈ N

that

∥

−φ

(

i−1

)∥

max

≤Ch

, (7)

whenever

∈ S for k = 0,1,...,i −1, and addition-

ally we have for i ≥ 3 and with j = 0,1,2,3, that

∥

i−j

−φ

−j

(

)∥

max

≤C

∗

. (8)

The assertions of the theorem clearly follow from (I).

To prove (I) let us ﬁrst ﬁx the constants. Let L > 0

be a Lipschitz constant for f on

S and set

C := 2 max{C

RK4

AB4

}, (9)

∗

:= e

(1 + e

+ e

)C, (10)

and

∗

:= min



C −C

AB4

∗



. (11)

We now show that (I) holds true for i = 0,1,2,3.

Indeed, (7) follows immediately from (5), and (8) fol-

lows immediately from the standard error estimate for

an explicit numerical integrator with local truncation

error Ch

since 0 < h ≤ h

∗

, see e.g. (Sauer, 2012).

Now assume that (I) holds true for all natural num-

bers up to and including some i ≥ 3. We assume that

∈ S and show that (I) also holds true for i + 1.

Let us ﬁrst consider (7) with i replaced by i + 1.

Observe that

∥

i+1

−φ

(

)∥

max

= ∥AB4

j=i

(

) −φ

(

)∥

max

≤ ∥AB4

j=i

(

) −AB4

j=0

(φ

(

))∥

max

+ ∥AB4

j=0

(φ

(

)) −φ

(

)∥

max

(12)

and for the second term on the right-hand-side we

have the bound

∥AB4

j=0

(φ

(

)) −φ

(

)∥

max

≤C

AB4

(13)

by (6).

To bound the ﬁrst term on the right-hand-side of

(12) we use the formula for AB4, that φ

(

) =

the Lipschitz condition on f on

S, and induction hy-

pothesis (8), and we get

∥AB4

j=i

(

) −AB4

j=0

(φ

(

))∥

max

(14)

∥−59f(

i−1

) + 59f(φ

−1

(

))

+ 37f(

i−2

) −37f(φ

−2

(

))

−9f(

i−3

) + 9f(φ

−3

(

))∥

max

≤



59∥f(

i−1

) −f(φ

−1

(

))∥

max

+ 37∥f(

i−2

) −f(φ

−2

(

))∥

max

+ 9∥f(

i−3

) −f(φ

−3

(

))∥

max



≤



59∥

i−1

−φ

−1

(

)∥

max

+ 37∥

i−2

−φ

−2

(

)∥

max

+ 9∥

i−3

−φ

−3

(

)∥

max



≤

105hL

∗

< 5C

∗

Hence, (12), (13), and (14) deliver

∥

i+1

−φ

(

)∥

max

≤5C

∗

AB4

≤Ch

, (15)

because

∗

Lh +C

AB4

≤ 5C

∗

AB4

≤C

by (11).

Hence, the bound (7) in (I) holds true for i re-

placed by i + 1.

Let us now consider the bound (8) in (I) for i re-

placed by i + 1.

The case j = 0 is obvious and from

∥

i+1−j

−φ

−j

(

i+1

)∥

max

= ∥φ

−j

(φ

(

i+1−j

)) −φ

−j

(

i+1

)∥

max

≤ e

jLh

∥φ

(

i+1−j

) −

i+1

∥

max

the case j = 1, i.e.

∥

−φ

−1

(

i+1

)∥

max

≤e

∥φ

(

)−

i+1

∥

max

≤C

∗

follows from (15) and e

C ≤C

∗

. Here we used the

well known

∥φ

φ(t,a) −φ

φ(t,b)∥

max

≤ e

L|t|

∥a −b∥

max

The cases j = 2 and j = 3 now follow similarly from

(15) and the induction hypothesis (8). For j = 2 we

have

∥

i−1

−φ

−2

(

i+1

)∥

max

≤ ∥

i−1

−φ

−1

(

)∥

max

+ ∥φ

−1

(

) −φ

−2

(

i+1

)∥

max

≤ e

∥φ

(

i−1

) −

∥

max

+ e

2Lh

∥φ

(

) −

i+1

∥

max

≤ e

(1 + e

)Ch

Positively Invariant Sets for ODEs and Numerical Integration

and

(1 + e

)Ch

≤C

∗

For j = 3 we have

∥

i−2

−φ

−3

(

i+1

)∥

max

≤ ∥

i−2

−φ

−1

(

i−1

)∥

max

+ ∥φ

−1

(

i−1

) −φ

−2

(

)∥

max

+ ∥φ

−2

(

) −φ

−3

(

i+1

)∥

max

≤ e

∥φ

(

i−2

) −

i−1

∥

max

+ e

2Lh

∥φ

(

i−1

) −

∥

max

+ e

3Lh

∥φ

(

) −

i+1

∥

max

≤ e

(1 + e

+ e

2Lh

)Ch

and

(1 + e

+ e

2Lh

)C = C

∗

≤C

∗

Thus, we have proved the bound (8) of (I) for i

replaced by i + 1, which concludes the proof.

Deﬁnition 2.2. (Order of numerical methods)

A numerical method to solve

x = f(x), f ∈C

), (16)

is said to be of order p ∈ N, if for any compact set

S ⊂ R

there exist constants C

∗

> 0 such that for

all step-sizes 0 < h ≤ h

∗

we have for any i ∈ N

that

∥

i+1

(ξ

ξ) −φ

φ(h,

(ξ

ξ))∥

≤C

p+1

whenever

(ξ

ξ),

(ξ

ξ),...,

(ξ

ξ) ∈ S.

With Deﬁnition 2.2, Theorem 2.1 can be formu-

lated as: Assume f ∈ C

) in (1). Then AB4

initialized with RK4 is of order 4 in the sense of Def-

inition 2.2.

Remark 2.3. It is straight forward to adapt the proof

of Theorem 2.1, under the assumption that f in (1)

is in C

), to the Adams-Bashforth method of

order p initialized with Runge-Kutta of the same or-

der. Hence, an AB-RK pair of order p is a numerical

method of order p in the sense of Deﬁnition 2.2.

We are now ready to study positively invariant sets

for the ODE (1), that are also positively invariant for

the numerical method.

3 POSITIVELY INVARIANT SETS

A positively invariant set for system (1), i.e. a set

P ⊂ R

such that φ

φ(t,x) ∈ P for all t ≥ 0 whenever

x ∈P, is not necessarily positively invariant for a nu-

merical procedure to approximate its solution trajec-

tories. The following example is quite revealing for

the general situation; we show for a simple system

and Euler’s Method that

• no matter how small the time-step h > 0 is, the dis-

crete semi-dynamical system deﬁned by Euler’s

Method does not necessarily have the same basins

of attraction as the original system, and

• if we restrict Euler’s Method to certain compact

and positively invariant subsets of the basins of at-

traction of the system, then these sets are also pos-

itively invariant for the discrete semi-dynamical

system deﬁned by Euler’s Method for sufﬁciently

small step sizes.

Recall that semi-dynamical systems are dynamical

systems, with the exception that solution trajectories

are not deﬁned for negative times.

Example 3.1. Consider the system

θ = 1 and ˙r =

−r(1 − r

) in polar coordinates; the origin is an

asymptotically stable equilibrium and the circular

disc B

:= {x ∈ R

: ∥x∥

< 1} is its basin of attrac-

tion. In Cartesian coordinates the equations of mo-

tion are



˙x

˙y





−x(1 −x

−y

) −y

−y(1 −x

−y

) + x



=: f(x,y). (17)

Using Euler’s method at z

:= (0, y) and with step-

size h > 0 delivers the approximation

z :=





+ h



−y

−y(1 −y

)



at time h. Now, for y

< 1, we have

∥z∥

= h

+ (y −hy(1 −y

))

> 1

0 < g(y,h)

:= y

(1 + (1 −y

)

| {z }

=:a>0

+ 2y

−1)

| {z }

=:b<0

h + y

−1

|{z}

=:c<0

i.e. for

h >

−b +

√

−4ac

> 0.

Because g is continuous and lim

y→±1

g(y,h) = h

, for

a ﬁxed h > 0 we can always ﬁnd y close enough to 1

(or −1), such that ∥z∥

> 1. Hence, B

is not posi-

tively invariant for this system when Euler’s method

is used, no matter how small h > 0 is.

Now ﬁx 0 < r < 1 and consider the compact set

:= {x ∈ R

: ∥x∥

≤ r}. Note that B

is positively

ICINCO 2023 - 20th International Conference on Informatics in Control, Automation and Robotics

invariant for the dynamical system deﬁned by (17). If

h > 0 satisﬁes

(1 + (1 −r

)

)h + 2(r

−1) ≤ 0 (18)

i.e. h ≤

2(1 −r

)

1 + (1 −r

)

then for all z

= (x,y) with ∥z

∥

≤ r one has

∥z∥

= ∥z

+ hf(z

)∥

≤ ∥z

∥

i.e. B

is positively invariant for Euler’s method with h

fuﬁlling (18). Note that condition (18) can be derived

from considering the special case z

= (0,y).

We now show in Theorem 3.2 and Corollary 3.5

that for general systems, sublevel-sets S ⊂ R

of cer-

tain Lyapunov-like functions V are positively invari-

ant for both the system (1) and numerical methods of

order p in the sense of Deﬁnition 2.2, p ∈ N

and

1 ≤ p ≤ s, to approximate its solution trajectories in S

and sufﬁciently small step size h > 0.

Theorem 3.2. (Positively invariant sets) Consider the

system (1), let V ∈C

;R) and assume S is a com-

pact connected component of {x ∈ R

: V (x) ≤ m},

m ∈R. Further assume that ∇V (x) ·f(x) < 0 and that

∇V (x) points out of S for every x ∈ ∂S. Then S is

positively invariant for (1).

Further, assume that f ∈ C

) and that we

have a numerical method of order p in the sense of

Deﬁnition 2.2. Then there is an h

′

> 0 such that if

the time-step h of the numerical method fulﬁlls 0 <

h ≤h

′

, then

i+1

(ξ

ξ) ∈ S, whenever

(ξ

ξ) ∈ S for k =

0,1,2,...,i, i ∈ N

Proof. Deﬁne

δ := −

max

x∈∂S

∇V (x) ·f(x) > 0

and let ε > 0 be such that ∇V (x) ·f(x) ≤ −δ < 0 for

all

x ∈ W := {x ∈R

: d(x, ∂S) ≤ 2ε},

where

d(x, K) := min

y∈K

∥x −y∥

for a compact set K.

To see that S is positively invariant for (1) consider

that if it is not, then some solution trajectory starting

in S must intersect ∂S at some point x and then leave

S, that is, there exists an x ∈ ∂S and an τ

∗

> 0 such

that φ

φ(τ,x) ∈ W \S for all 0 < τ ≤ τ

∗

. Then

m < V (φ

φ(τ

∗

,x)) = V (x) +

∗

dτ

V (φ

φ(τ,x))dτ

= m +

∗

∇V (φ

φ(τ,x)) ·

dτ

φ(τ,x)dτ ≤ m −δτ

∗

a contradiction.

In order to prove the desired property for the nu-

merical method, set

V := {x ∈R

: d(x, ∂S) ≤ ε} ⊂ W ,

F := max{max

x∈S

∥f(x)∥

,1}, and

′

:= min{ε/(2max{F,C

}),1,τ

∗

Here and later in the proof C

∗

> 0 are the constants

for the numerical method from Deﬁnition 2.2.

Then, for x ∈ S \V and 0 ≤ h ≤ h

′

we have

∥φ

φ(h,x) −x∥

≤

∥f(φ

φ(s,x))∥

ds ≤ hF

and it follows that

d(φ

φ(h,x),S \V ) ≤ ε/2, ∀x ∈S \V . (19)

Note that from (19), Deﬁnition 2.2 and for

(ξ

ξ) ∈

S for k = 0, 1, 2, . . . , i, x :=

(ξ

ξ), we have

i+1

(ξ

ξ),S \V ) ≤ d(φ

φ(h,x),S \V )

+ ∥φ

φ(h,x) −

i+1

(ξ

ξ)∥

≤ ε/2 + ε/2 = ε.

Hence,

i+1

(x) ∈ S if the time-step of the numerical

method fulﬁlls 0 ≤ h ≤ h

′

, whenever

(ξ

ξ) ∈ S for

k = 0, 1, 2,...,i and

(ξ

ξ) = x ∈ S \V .

To ﬁnish the proof we need to show the statement

in the case

(ξ

ξ) ∈ S for k = 0, 1, 2, . . . , i and

x =

(ξ

ξ) ∈ S ∩V .

We assume on the contrary that there are sequences

∈ S and 0 < h

≤ h

′

, h

→ 0 as j → ∞, such that

(ξ

) /∈ S

for all j, although

(ξ

),...,

(ξ

) ∈S, x

(ξ

) ∈S∩V .

Here

(ξ

),...,

(ξ

)

is the sequence generated by the numerical method

with initial value ξ

∈ S and step-size h

Note that since S is positively invariant for (1)

we have φ

φ(h

) ∈ S and therefore V (x

) ≤ m and

V (φ

φ(h

)) ≤m for all j. Further, there exists I ∈N

such that for all j ≥ I we have φ

φ(θh

) ∈ W ∩S

for all θ ∈ [0,1] and V (

(ξ

)) > m. Moreover,

there is a convex and compact set

S ⊃ S such that

(ξ

) ∈

S for all j. Let L

be a Lipschitz con-

stant for V on

S and recall that by Deﬁnition 2.2 we

have

∥

(ξ

) −φ

φ(h

)∥

≤C

p+1

Positively Invariant Sets for ODEs and Numerical Integration

Now



V (

(ξ

)) −V (φ

φ(h

))



≤

∥

(ξ

) −φ

φ(h

)∥

≤

p+1

= L

V (

(ξ

)) −V (x

)

m −m

= 0,

and

V (

(ξ

)) −V (x

)

V (

(ξ

)) −V (φ

φ(h

))

V (φ

φ(h

)) −V (x

)

(20)

for all j ≥ I.

Deﬁne g

(t) := V (φ

φ(t,x

)). By the Mean-Value

theorem there exist θ

∈ (0,1) such that

) −g

(0) = g

′

(θ

≤ −δh

holds for all j ≥ I. From (20) it follows that

0 ≤ limsup

j→∞

V (

(ξ

)) −V (x

)

≤ limsup

j→∞

V (

(ξ

)) −V (φ

φ(h

))

+ limsup

j→∞

V (φ

φ(h

)) −V (x

)

≤ 0 −δ < 0,

which is a contradiction and thus the theorem is

proved.

Remark 3.3. Note that V (x,y) = x

+ y

is a Lya-

punov function for system (17) on B

and for 0 < r < 1

the set V

−1

([0,r

]) = {x ∈ R

: ∥x∥

≤ r} ⊂ B

is a

compact sublevel set that is positively invariant.

Further,

∇V (x,y) ·f(x,y) = −2r

(1 −r

) < 0

for x

+ y

= r

, i.e. ∇V (x, y) ·f(x,y) is less than a

negative constant at the boundary of V

−1

([0,r

]). The

last theorem then tells us that there must be an h

′

> 0

so small that V

−1

([0,r

]) is positively invariant for

Euler’s Method with step size 0 < h ≤h

′

too.

In applications it can be more convenient to have

the following version of Theorem 3.2, which can be

used for Lyapunov functions computed by ﬁrst ap-

proximately solving the Zubov’s PDE

∇V (x) ·f(x) = −

+ ∥f(x)∥

using generalized interpolation in reproducing kernel

Hilbert spaces and then interpolating the values over

the simplices of a triangulation T , see (Giesl and Haf-

stein, 2015a; Giesl et al., 2021) for more details. Note

that in the following theorem, CPA[T ] denotes the

set of these interpolating functions, called continuous

piecewise afﬁne (CPA) functions, which are afﬁne on

each simplex of the triangulation T and continuous

overall. Further, ∇V

∈ R

denotes the constant gra-

dient of V in the interior of a simplex S

∈ T .

Theorem 3.4 (CPA version of Thm. 3.2). Let V ∈

CPA[T ] and assume there is a compact connected

component S of {x ∈R

: V (x) ≤m}, m ∈ R. Assume

that ∇V

points out of S at x for every x ∈∂S∩S

and

every S

∈T , and that there is a constant c > 0 such

that ∇V

·f(x) ≤−c for every x in a neighbourhood of

∂S and every ν such that x ∈S

. Then S is positively

invariant for (1).

Further, assume that f ∈ C

) and that we

have a numerical method of order p in the sense of

Deﬁnition 2.2. Then there is an h

′

> 0 such that if

the time-step h of the numerical method fulﬁlls 0 <

h ≤h

′

, then

i+1

(ξ

ξ) ∈ S, whenever

(ξ

ξ) ∈ S for k =

0,1,2,...,i.

Proof. Essentially, the proof is the same as the proof

of Theorem 3.2; the existence of δ = c/2 and ε > 0

now follow directly from the assumptions. The only

reasoning that needs modiﬁcation is why

) −g

(0) ≤ −δh

For V ∈ CPA[T ] it follows because for

V (x) := lim sup

h→0+

V (φ

φ(h,x)) −V (x)

we have

V (x) ≤ max

ν: x∈S

∇V

·f(x) ≤ −δ,

see e.g. (Hafstein, 2020, Lem. 2.2), and by a gener-

alized Mean Value Theorem, see (Scheeffer, 1884) or

(Walter, 2004, Thm. 12.24).

Finally, we can state and prove the main result of

this paper.

Theorem 3.5. (Positively invariant sets for the system

and the numerical method) Let x

be an exponentially

stable equilibrium of (1), where f ∈C

) with

ICINCO 2023 - 20th International Conference on Informatics in Control, Automation and Robotics

p ∈ N, and let K ⊂ A(x

) be compact. Then there

exists a compact and connected set S, K ⊂S ⊂A(x

with the following property:

Assume we have a numerical method of order p

in the sense of Deﬁnition 2.2. Then there exists a con-

stant h

′

> 0, such that S is positively invariant both for

the original ﬂow φ

φ(0,ξ

ξ) = ξ

ξ ∈ S, t 7→ φ

φ(t,ξ

ξ), induced

by (1), and for the sequences

(ξ

ξ), i ∈ N

, generated

by the numerical method with step-size h, 0 < h ≤h

′

for the initial-values ξ

ξ ∈ S. In other words,

(ξ

ξ) ∈ S

for all ξ

ξ ∈ S and all i ∈N

Proof. By (Giesl, 2007, Thm. 2.46) there exists a

Lyapunov function V for the system (1) fulﬁlling

∇V (x) ·f(x) = −∥x −x

∥

1 + ∥f(x)∥

for all x ∈ A(x

). We set r := max

x∈K

V (x) and S :=

−1

([0,r]). Since V is also a Lyapunov function for

the system

x = f(x)(1 + ∥f(x)∥

)

−1/2

with a bounded right-hand-side, the set S ⊂ A(x

) is

compact. Since V (φ

φ(t,ξ

ξ)) ≤ r for all t ≥ 0 and

∈ φ

φ([0,∞),ξ

ξ) ⊂V

−1

([0,r]) = S

for all ξ

ξ ∈ S, the set S is also connected. Using this

Lyapunov function and Theorem 3.2 for the numeri-

cal method, the existence of h

′

> 0 with the claimed

properties follows.

Remark 3.6. By Theorem 2.1 the AB4 method initial-

ized by RK4 is a numerical method of order 4 in the

sense of Deﬁnition 2.2, and thus, fulﬁlls the assump-

tions in Theorem 3.5. By Remark 2.3 the same applies

to the Adams-Bashforth method of any order, initial-

ized with the Runge-Kutta method of the same order.

These results are used in (Giesl and Hafstein, 2023)

and (Giesl et al., 2023c) to prove that Lyapunov func-

tions and contraction metrics can be approximated

arbitrarily close on compact sets using numerical in-

tegration and quadrature.

4 CONCLUSIONS

We have shown that for an ODE with an exponen-

tially stable equilibrium x

and any compact subset K

of its basin of attraction A(x

), we can ﬁnd a compact

and connected set S, K ⊂S ⊂ A(x

), that is positively

invariant, both for the ODE and its numerical approx-

imation. We considered the concrete case of using the

Adams-Bashforth method of fourth order, initialized

with the usual Runge-Kutta method of fourth order,

but we also discussed obvious extensions to an arbi-

trary order. Finally, we demonstrated how such posi-

tively invariant sets can be computed in practice.

ACKNOWLEDGEMENT

The research done for this paper was partially sup-

ported by the Icelandic Research Fund in the grant

228725-051, Linear Programming as an Alternative

to LMIs for Stability Analysis of Switched Systems,

which is gratefully acknowledged.

REFERENCES

Aminzare, Z. and Sontag, E. (2014). Contraction methods

for nonlinear systems: A brief introduction and some

open problems. In Proceedings of the 53rd IEEE Con-

ference on Decision and Control, pages 3835–3847.

Anderson, J. and Papachristodoulou, A. (2015). Advances

in computational Lyapunov analysis using sum-of-

squares programming. Discrete Contin. Dyn. Syst. Ser.

B, 20(8):2361–2381.

Aylward, E., Parrilo, P., and Slotine, J.-J. (2008). Stability

and robustness analysis of nonlinear systems via con-

traction metrics and SOS programming. Automatica,

44(8):2163–2170.

ornsson, J., Giesl, P., Hafstein, S., Kellett, C., and Li,

H. (2014). Computation of continuous and piecewise

afﬁne Lyapunov functions by numerical approxima-

tions of the Massera construction. In Proceedings

of the CDC, 53rd IEEE Conference on Decision and

Control, pages 5506–5511, Los Angeles (CA), USA.

ornsson, J., Giesl, P., Hafstein, S., Kellett, C., and Li, H.

(2015). Computation of Lyapunov functions for sys-

tems with multiple attractors. Discrete Contin. Dyn.

Syst. Ser. A, 35(9):4019–4039.

ornsson, J. and Hafstein, S. (2017). Efﬁcient Lya-

punov function computation for systems with multiple

exponentially stable equilibria. Procedia Computer

Science, 108:655–664. Proceedings of the Interna-

tional Conference on Computational Science (ICCS),

Zurich, Switzerland, 2017.

Borg, G. (1960). A condition for the existence of orbitally

stable solutions of dynamical systems. Kungl. Tekn.

ogsk. Handl. 153.

Demidovi

c, B. (1961). On the dissipativity of a certain

non-linear system of differential equations. I. Vestnik

Moskov. Univ. Ser. I Mat. Meh., 1961(6):19–27.

Doban, A. (2016). Stability domains computation and sta-

bilization of nonlinear systems: implications for bio-

logical systems. PhD thesis: Eindhoven University of

Technology.

Doban, A. and Lazar, M. (2016). Computation of Lyapunov

functions for nonlinear differential equations via a

Yoshizawa-type construction. IFAC-PapersOnLine,

49(18):29 – 34.

Positively Invariant Sets for ODEs and Numerical Integration

Forni, F. and Sepulchre, R. (2014). A differential Lyapunov

framework for Contraction Analysis. IEEE Transac-

tions on Automatic Control, 59(3):614–628.

Giesl, P. (2007). Construction of Global Lyapunov Func-

tions Using Radial Basis Functions. Lecture Notes in

Math. 1904, Springer.

Giesl, P. (2015). Converse theorems on contraction metrics

for an equilibrium. J. Math. Anal. Appl., (424):1380–

1403.

Giesl, P. (2019). Computation of a contraction metric for

a periodic orbit using meshfree collocation. SIAM J.

Appl. Dyn. Syst., 18(3):1536–1564.

Giesl, P. and Hafstein, S. (2013). Construction of a CPA

contraction metric for periodic orbits using semideﬁ-

nite optimization. Nonlinear Anal., 86:114–134.

Giesl, P. and Hafstein, S. (2014). Revised CPA method to

compute Lyapunov functions for nonlinear systems. J.

Math. Anal. Appl., 410:292–306.

Giesl, P. and Hafstein, S. (2015a). Computation and veri-

ﬁcation of Lyapunov functions. SIAM J. Appl. Dyn.

Syst., 14(4):1663–1698.

Giesl, P. and Hafstein, S. (2015b). Review of computa-

tional methods for Lyapunov functions. Discrete Con-

tin. Dyn. Syst. Ser. B, 20(8):2291–2331.

Giesl, P. and Hafstein, S. (2023). Dynamical Systems in

Theoretical Perspective, chapter Lyapunov Functions

by Interpolating Numerical Quadrature: Proof of Con-

vergence, page (to appear). Springer. Springer Pro-

ceedings in Mathematics and Statistics.

Giesl, P., Hafstein, S., Haraldsdottir, M., Thorsteinsson, D.,

and Kawan, C. (2023a). Subgradient algorithm for

computing contraction metrics for equilibria. J. Com-

put. Dynamics, 10(2):281–303.

Giesl, P., Hafstein, S., and Kawan, C. (2023b). Review on

contraction analysis and computation of contraction

metrics. J. Comput. Dynamics, 10(1):1–47.

Giesl, P., Hafstein, S., and Mehrabinezhad, I. (2021). Com-

putation and veriﬁcation of contraction metrics for pe-

riodic orbits. J. Math. Anal. Appl., 503(2):Paper No.

125309, 32.

Giesl, P., Hafstein, S., and Mehrabinezhad, I. (2023c). Con-

traction metrics by numerical integration and quadra-

ture: Uniform error estimate. In Proceedings of the

20th International Conference on Informatics in Con-

trol, Automation and Robotics (ICINCO), page (sub-

mitted).

Gudmundsson, S. and Hafstein, S. (2015). Lyapunov func-

tion veriﬁcation: MATLAB implementation. In Pro-

ceedings of the 1st Conference on Modelling, Identiﬁ-

cation and Control of Nonlinear Systems (MICNON),

number 0235, pages 806–811.

Hafstein, S. (2004). A constructive converse Lyapunov the-

orem on exponential stability. Discrete Contin. Dyn.

Syst. Ser. A, 10(3):657–678.

Hafstein, S. (2005). A constructive converse Lyapunov

theorem on asymptotic stability for nonlinear au-

tonomous ordinary differential equations. Dynamical

Systems: An International Journal, 20(3):281–299.

Hafstein, S. (2019a). Computational Science - ICCS 2019:

19th International Conference, Faro, Portugal, June

12-14, 2019, Proceedings, Part V, chapter Numerical

Analysis Project in ODEs for Undergraduate Students,

pages 412–434. Springer.

Hafstein, S. (2019b). Numerical ODE solvers and integra-

tion methods in the computation of CPA Lyapunov

functions. In Proceedings of the 18th European Con-

trol Conference (ECC), pages 1136–1141.

Hafstein, S. (2020). CPA Lyapunov functions: Switched

systems vs. differential inclusions. In Proceedings of

the 17th International Conference on Informatics in

Control, Automation and Robotics (ICINCO), pages

745–753.

Hafstein, S., Kellett, C., and Li, H. (2014a). Computation

of Lyapunov functions for discrete-time systems using

the Yoshizawa construction. In Proceedings of 53rd

IEEE Conference on Decision and Control (CDC).

Hafstein, S., Kellett, C., and Li, H. (2014b). Continu-

ous and piecewise afﬁne Lyapunov functions using the

Yoshizawa construction. In Proceedings of the 2014

American Control Conference (ACC), pages 548–553

(no. 0170), Portland (OR), USA.

Hafstein, S., Kellett, C., and Li, H. (2015). Computing con-

tinuous and piecewise afﬁne Lyapunov functions for

nonlinear systems. J. Comput. Dyn, 2(2):227 – 246.

Hafstein, S. and Valfells, A. (2017). Study of dynamical

systems by fast numerical computation of Lyapunov

functions. In Proceedings of the 14th International

Conference on Dynamical Systems: Theory and Ap-

plications (DSTA), volume Mathematical and Numer-

ical Aspects of Dynamical System Analysis, pages

220–240.

Hafstein, S. and Valfells, A. (2019). Efﬁcient computation

of Lyapunov functions for nonlinear systems by in-

tegrating numerical solutions. Nonlinear Dynamics,

97(3):1895–1910.

Hahn, W. (1967). Stability of Motion. Springer, Berlin.

Hartman, P. (1961). On stability in the large for systems

of ordinary differential equations. Canadian J. Math.,

13:480–492.

Hartman, P. (1964). Ordinary Differential Equations. Wi-

ley, New York.

Julian, P., Guivant, J., and Desages, A. (1999). A

parametrization of piecewise linear Lyapunov func-

tions via linear programming. Int. J. Control, 72(7-

8):702–715.

Kamyar, R. and Peet, M. (2015). Polynomial optimization

with applications to stability analysis and control – an

alternative to sum of squares. Discrete Contin. Dyn.

Syst. Ser. B, 20(8):2383–2417.

Khalil, H. (2002). Nonlinear Systems. Pearson, 3. edition.

Krasovski

i, N. N. (1963). Problems of the Theory of Stabil-

ity of Motion. Mir, Moskow, 1959. English translation

by Stanford University Press.

Lewis, D. (1951). Differential equations referred to a vari-

able metric. Amer. J. Math., 73:48–58.

Lewis, D. C. (1949). Metric properties of differential equa-

tions. Amer. J. Math., 71:294–312.

Li, H., Hafstein, S., and Kellett, C. (2015). Computation of

continuous and piecewise afﬁne Lyapunov functions

ICINCO 2023 - 20th International Conference on Informatics in Control, Automation and Robotics

for discrete-time systems. J. Difference Equ. Appl.,

21(6):486–511.

Lohmiller, W. and Slotine, J.-J. (1998). On Contrac-

tion Analysis for Non-linear Systems. Automatica,

34:683–696.

Marin

osson, S. (2002). Lyapunov function construction for

ordinary differential equations with linear program-

ming. Dynamical Systems: An International Journal,

17:137–150.

Parrilo, P. (2000). Structured Semideﬁnite Programs and

Semialgebraic Geometry Methods in Robustness and

Optimization. PhD thesis: California Institute of

Technology Pasadena, California.

Ratschan, S. and She, Z. (2010). Providing a basin of attrac-

tion to a target region of polynomial systems by com-

putation of Lyapunov-like functions. SIAM J. Control

Optim., 48(7):4377–4394.

Sastry, S. (1999). Nonlinear Systems: Analysis, Stability,

and Control. Springer.

Sauer, T. (2012). Numerical Analysis. Pearson, 2nd edition.

Scheeffer, L. (1884). Zur Theorie der stetigen Funktionen

einer reellen Ver

anderlichen. Acta Math., 5(1):183–

194, 279–296.

Simpson-Porco, J. and Bullo, F. (2014). Contraction the-

ory on Riemannian manifolds. Systems Control Lett.,

65:74–80.

Valmorbida, G. and Anderson, J. (2017). Region of attrac-

tion estimation using invariant sets and rational Lya-

punov functions. Automatica, 75:37–45.

Vannelli, A. and Vidyasagar, M. (1985). Maximal Lya-

punov functions and domains of attraction for au-

tonomous nonlinear systems. Automatica, 21(1):69–

80.

Vidyasagar, M. (2002). Nonlinear System Analysis. Clas-

sics in Applied Mathematics. SIAM, 2. edition.

Walter, W. (2004). Analysis 1. Springer, 7th edition.

Yoshizawa, T. (1966). Stability theory by Liapunov’s sec-

ond method. Publications of the Mathematical Society

of Japan, No. 9. The Mathematical Society of Japan,

Tokyo.

Zubov, V. I. (1964). Methods of A. M. Lyapunov and their

application. Translation prepared under the auspices

of the United States Atomic Energy Commission;

edited by Leo F. Boron. P. Noordhoff Ltd, Groningen.

Positively Invariant Sets for ODEs and Numerical Integration