Journal of Information & Computational Science 11:17 (2014) 6175–6183
Available at http://www.joics.com
November 20, 2014
Constructing the Lyapunov Function by
Quasi-interpolation ?
Wenwu Gao a,b,∗
a School
of Economics, Anhui University, Hefei 230411, China
b Shanghai
Key Laboratory for Contemporary Applied Mathematics
School of Mathematical Sciences, Fudan University, Shanghai 200433, China
Abstract
The Lyapunov function plays a vital role in dynamical systems. One method of constructing the
Lyapunov function is the meshless collocation by radial basis functions. However, the meshless collocation
method gives a non-local Lyapunov function (with unwanted nonnegative orbital derivative in some
neighborhood of an equilibrium). To overcome the difficulty of the method, the paper proposes a scheme
for constructing the Lyapunov function by multiquadric trigonometric B-spline quasi-interpolation. The
scheme is efficient, simple and easy to compute. More importantly, it provides two additional shape
parameters for the constructed Lyapunov function. This implies that, by adjusting these two shape
parameters suitably, one can get different shapes of the Lyapunov function and thus obtain different
basins of attraction for a dynamical system. Moreover, by taking the union of these different basins, one
can even obtain a larger basin of attraction that provides better approximation to the exact one.
Keywords: Quasi-interpolation; Dynamical System; Lyapunov Function; Basin of Attraction; Meshless
Collocation
1
Introduction
Estimation of the basin of attraction is an important and hot topic in dynamical systems. A lot
of methods have been proposed in the literature (see [1, 2] and the reference therein). Among
them, the Lyapunov function method is widely used.
From the Lyapunov function point of view, the basin of attraction of an asymptotically stable
equilibrium for an autonomous dynamical system can be determined by sublevel sets of a Lyapunov function. Thus once a Lyapunov function has been constructed, the basin of attraction
can be obtained accordingly. There are several methods for constructing the Lyapunov function,
e.g., by numerical solutions of the Zubov’s equations [3], by optimizations [4], by radial basis
functions [5], etc.
?
Project supported by the National Nature Science Foundation of China Key Project (No. 91330201) and SGST
(No. 12DZ 2272800).
∗
Corresponding author.
Email address: [email protected] (Wenwu Gao).
1548–7741 / Copyright © 2014 Binary Information Press
DOI: 10.12733/jics20104898
6176
W. Gao / Journal of Information & Computational Science 11:17 (2014) 6175–6183
In particular, Giesl [6] constructed the Lyapunov function by the symmetric collocation method
using radial basis function, which was introduced in the context of Herimte-Birkhoff data interpolation [7]. For more details, we refer readers to [6].
However, the collocation method can not preserve the shapes. This may be the reason why
Giesl [8] pointed out that it would give a non-local Lyapunov function (with unwanted nonnegative
orbital derivative in some neighborhood of an equilibrium).
To overcome the difficulty, Giesl [8] proposed a mend by using the Taylor’s expansion in a small
neighborhood of an equilibrium coupled together with the symmetric collocation method. Giesl
also proposed some other techniques of modifying the method, see, e.g., Giesl [5], and Giesl and
Wendland [9] for instance. Nevertheless, all of the above methods require solving a large-scale
system of equations, which may be ill-conditioned for a large number of collocation points.
Motivated from Giesl’s works and the fair properties of quasi-interpolation, e.g., preserving
shapes better than interpolation (see [10] for instance), yielding solutions directly without the
need to solve any large-scale system of equations, etc, we use quasi-interpolation to construct the
Lyapunov function directly.
Quasi-interpolation has been widely discussed, for example, [11]-[15], [17]-[21] and so on. Particularly, some authors constructed quasi-interpolants for periodic data. Lyche, Schumaker and
Stanley [22] constructed a trigonometric B-spline quasi-interpolant with periodicity 2π. They also
discussed its approximation orders for high-order derivatives. Recently, by using a multiquadric
trigonometric function as the kernel, [23] provided a multiquadric trigonometric B-spline quasiinterpolant. Later, [24] derived its approximation orders for high-order derivatives and shape
preserving properties.
The paper goes further with the quasi-interpolant by discussing its applications in constructing
the Lyapunov functions and thus estimating basins of attraction for dynamical systems.
We focus on the two-dimensional autonomous dynamical system:
(
x˙ = f1 (x, y),
y˙ = f2 (x, y).
(1)
Here f1 (x, y), f2 (x, y) are some smooth functions such that (x0 , y0 ) is an asymptotically stable
equilibrium. According to Giesl [6], there exists a unique Lyapunov function V (x, y) satisfying
the equation
(
Vx (x, y)f1 (x, y) + Vy (x, y)f2 (x, y) = −(x − x0 )2 − (y − y0 )2 ,
V (x0 , y0 ) = V0 .
Moreover, by introducing the generalized polar coordinates transformations x = x0 + ar cos θ,
y = y0 + br sin θ, the above equation becomes
(
Vr (r, θ)F1 (r, θ) + Vθ (r, θ)F2 (r, θ) = −r2 ,
(2)
V (0, θ) = 0.
Here a, b, r ∈ R+ , θ ∈ [0, 2π), and
F1 (r, θ) =
1
(bf1 (r, θ) cos θ + af2 (r, θ) sin θ),
ab
W. Gao / Journal of Information & Computational Science 11:17 (2014) 6175–6183
F2 (r, θ) =
6177
1
(af2 (r, θ) cos θ − bf1 (r, θ) sin θ).
abr
To construct the Lyapunov function from the equation (2), we approximate the radial derivative with the first-order forward divided differences and use the first-order derivative of the multiquadric trigonometric B-spline quasi-interpolant to approximate the angular derivative.
The scheme is simple, efficient and easy to compute. Moreover, it overcomes the difficulty
of the symmetric collocation method (i.e., giving a non-local Lyapunov function with unwanted
nonnegative orbital derivative in some neighborhood of an equilibrium). More importantly, it
also provides two shape parameters (i.e., a, b) and thus gives different basins of attraction (by
adjusting these two parameters). Particularly, one can even get a larger basin of attraction by
taking the union of these different basins.
This paper is organized as follows.
To make the paper self-contained, some preliminaries about dynamical systems and multiquadric trigonometric B-spline quasi-interpolation are presented in Section 2. The main scheme
for constructing the Lyapunov function is provided in Section 3. Numerical examples are presented in Section 4. Conclusions and discussions are given in Section 5.
2
Prelinimaries
2.1
A Short Introduction to Dynamical Systems
This subsection will introduce some basic definitions of dynamical systems (see Giesl [5]).
Definition 1 (Flow operator) Let (x(t), y(t)) be the trajectory of the system (1) starting from
an initial point (ξ1 , ξ2 ), i.e., (x(0), y(0)) = (ξ1 , ξ2 ). Then the flow operator St is defined by
St (ξ1 , ξ2 ) = (x(t), y(t)).
Definition 2 (Equilibrium) A point (x0 , y0 ) ∈ R2 is called an equilibrium of the system (1), if
both f1 (x0 , y0 ) = 0 and f2 (x0 , y0 ) = 0.
Note that if (x0 , y0 ) is an equilibrium, then St (x0 , y0 ) = (x0 , y0 ) for all t ≥ 0, that is, the
trajectory starting from an equilibrium will stay at it for all the time.
Definition 3 (Equilibrium) A point (x0 , y0 ) ∈ R2 is called an equilibrium of the system (1), if
both f1 (x0 , y0 ) = 0 and f2 (x0 , y0 ) = 0.
Definition 4 (Stability and attractivity) Let (x0 , y0 ) be an equilibrium.
• (x0 , y0 ) is called stable, if for all ² > 0, there exists a δ > 0 such that
St (x, y) ∈ B² := {(x, y) ∈ R2 : k(x, y) − (x0 , y0 )k ≤ ²}
holds for all (x, y) ∈ Bδ ((x0 , y0 )) and for all t ≥ 0. Here and below, k · k denotes the
Euclidean norm in R2 .
6178
W. Gao / Journal of Information & Computational Science 11:17 (2014) 6175–6183
• (x0 , y0 ) is called attractive, if there exists a δ 0 > 0 such that for all (x, y) ∈ Bδ0 ((x0 , y0 )),
kSt (x, y) − (x0 , y0 )k converges to zero as t tends to infinity.
• (x0 , y0 ) is called asymptotically stable, if it is both stable and attractive.
Definition 5 (Basin of attraction) The basin of attraction of an asymptotically stable equilibrium
(x0 , y0 ) is defined as
A(x0 , y0 ) = {(x, y) ∈ R2 : lim St (x, y) = (x0 , y0 )}.
t→∞
From the definition, one can find that the basin of attraction consists of all the initial points such
that the trajectories starting from them approach the equilibrium as time tends to infinity.
Definition 6 (Orbital derivative) Given a function V (x, y) ∈ C 1 (R2 , R), its orbital derivative
(denoted as V 0 (x, y)) with respect to the system (1) is defined as the derivative along a trajectory
of the system, i.e., V 0 (x, y) = Vx (x, y)f1 (x, y) + Vy (x, y)f2 (x, y).
Definition 7 (Lyapunov function) Let (x0 , y0 ) be an equilibrium and let B be an open set containing (x0 , y0 ) as its interior. A function V (x, y) ∈ C 1 (B, R) is called the Lyapunov function if there exists a set K ⊂ B such that the orbital derivative V 0 (x, y) is negative for all
(x, y) ∈ K\{(x0 , y0 )}.
Based on the definitions of the orbital derivative and the Lyapunov function, one can find that
the Lyapunov function decreases strictly along the trajectories of a dynamical system and attains
its minimum at the equilibrium. Moreover, for an asymptotically stable equilibrium of an autonomous dynamical system, its basin of attraction can be determined through sublevel sets of a
Lyapunov function.
2.2
Multiquadric Trigonometric B-spline Quasi-interpolation
The construction of multiquadric trigonometric B-spline quasi-interpolation consists of three
steps [23].
Firstly, the kernel (multiquadric trigonometric function) is constructed:
q
φ(x) = c2 + sin2 x/2,
where c is a positive shape parameter.
Secondly, with respect to the kernel, the multiquadric trigonometric B-splines can be given in
the form
ψj (x) =
xj+1 −xj
φ(x −
2
xj+1 −xj
xj+1 −xj−1
sin
2
2
φ(x − xj+1 ) − cos
2 sin
xj )
−
cos
xj −xj−1
φ(x − xj ) − φ(x − xj−1 )
2
.
xj −xj−1
x
−x
2 sin 2 sin j+1 2 j−1
Finally, with these {ψj (x)} being the basis, the multiquadric trigonometric B-spline quasiinterpolation can be constructed as
Qf (x) =
N
X
j=1
sin
xj+1 − xj−1
f (xj )ψj (x),
2
(3)
W. Gao / Journal of Information & Computational Science 11:17 (2014) 6175–6183
6179
where 0 = x0 < · · · < xN = 2π < xN +1 = 2π + x1 .
Moreover, its approximation orders and optimal choices of the shape parameter for high-order
derivatives can be given in the following theorem [24].
Theorem 1 Let k be any given nonnegative integer and f (x) be a C k+2 continuous function
with periodicity 2π. Let Qf (x) be defined in (3) and h be the density of the ascending centers
1
{xj }, i.e., h = maxj {xj+1 − xj }. Then, by choosing the shape parameter c as c = O(h k+1 ), the
approximation orders of (Qf )(k) (x) to f (k) (x) can be given as
2
k(Qf )(k) − f (k) k∞ ≤ O(h k+1 ).
(4)
From the theorem, one can find that the quasi-interpolant approximates high-order derivatives
well and thus can be applied in numerical solutions of partial differential equations. In the
following section, we shall apply it in constructing the Lyapunov function for a dynamical system.
3
Construction of the Lyapunov Function
This section aims at providing a scheme for constructing the Lyapunov function by the multiquadric trigonometric B-spline quasi-interpolation and thus getting an estimation of the basin of
attraction for a dynamical system.
We consider the two-dimensional autonomous dynamical system:
(
x˙ = f1 (x, y),
y˙ = f2 (x, y),
where f1 (x, y), f2 (x, y) are some smooth functions such that (x0 , y0 ) is an asymptotically
stable equilibrium. Moreover, by introducing the generalized polar coordinates transformations
x = x0 + ar cos θ, y = y0 + br sin θ, where a, b, r ∈ R+ and θ ∈ [0, 2π), the dynamical system in
the polar form becomes
(
Vr (r, θ) = F1 (r, θ),
Vθ (r, θ) = F2 (r, θ),
where
1
(bf1 (r, θ) cos θ + af2 (r, θ) sin θ),
ab
1
F2 (r, θ) =
(af2 (r, θ) cos θ − bf1 (r, θ) sin θ).
abr
F1 (r, θ) =
Furthermore, based on Giesl [6], there exists a unique Lyapunov function V (r, θ) satisfying the
equation
(
Vr (r, θ)F1 (r, θ) + Vθ (r, θ)F2 (r, θ) = −r2 ,
(5)
V (0, θ) = 0.
We now provide a scheme for solving the equation numerically by multiquadric trigonometric
B-spline quasi-interpolation.
6180
W. Gao / Journal of Information & Computational Science 11:17 (2014) 6175–6183
Let M , N be two given positive integers. Denote by {(rn , θk )}, n = 0, 1, · · · , M , k =
0, 1, · · · , N , the scattered centers satisfying
0 = r0 < r1 < r2 < · · · < rM −1 < rM ,
0 = θ0 < θ1 < θ2 < · · · < θN −1 < θN = 2π, θN +1 = 2π + θ1 .
Let Vkn be approximations of V (r, θ) at the centers (rn , θk ).
The scheme for solving the equation (5) consists of three steps.
First, at the radius rn , n = 0, 1, 2, · · · , M − 1, discretizing the radial derivative of the equation
Vr (r, θ)F1 (r, θ) + Vθ (r, θ)F2 (r, θ) = −r2
with the first-order forward divided differences yields
Vkn+1 = Vkn −
(rn+1 − rn )(Vθ (rn , θk )F2 (rn , θk ) + rn2 )
.
F1 (rn , θk )
(6)
Second, using the first-order derivative of the multiquadric trigonometric B-spline quasi-interpolant
(3) to approximate Vθ (rn , θk ) in equation (6) gives an iterative formula
Vkn+1
=
Vkn
−
(rn+1 − rn )(F2 (rn , θk )
PN
Vjn sin
F1 (rn , θk )
j=1
θj+1 −θj−1 0
ψj (θk )
2
+ rn2 )
.
This implies that the approximations {Vkn+1 } at the radius rn+1 can be obtained by the ones at
the radius rn through the iterative formula.
Finally, solving the iterative formula step by step gives these approximations.
Remark 1 The scheme does not require solving any linear system of equations. Moreover, for a
given θ, it needs only to compute the first-order derivative ψj0 (θ), j = 1, 2, · · · , N, once for all.
Thus the scheme saves a lot of the CPU time for computation.
Up to now, we have proposed a scheme for constructing the Lyapunov function. The scheme is
simple, efficient and easy to compute. In the following section, we shall present some examples
to illustrate the validity of the scheme.
4
Numerical Examples
The section gives some numerical examples of constructing the Lyapunov functions for some
dynamical systems with the above scheme.
Example 1
The first example is the system (that is the Example 6.1) discussed in Giesl [8]:

 x˙ = −x + x3 ,
y
 y˙ = − + x2 .
2
W. Gao / Journal of Information & Computational Science 11:17 (2014) 6175–6183
6181
The system has an asymptotically stable equilibrium (0, 0) whose basin of attraction can be given
directly as
A(0, 0) = {(x, y) ∈ R2 : |x| < 1}.
Thus, we can compare our estimation with the exact one.
Fig. 1 shows the boundaries of exact basin of attraction (the two red lines) and sublevel sets of
the reconstructed Lyapunov function with the shape parameters a = 1, b = 4.
5
4
3
2
1
0
−1
−2
−3
−4
−5
−5 −4 −3 −2 −1
0
1
2
3
4
5
Fig. 1: Sublevel sets of the Lyapunov function and the boundaries of the exact basin of attraction
From the figure, we can find that our scheme gives a Lyapunov function with negative orbital
derivative in the basin of attraction except the equilibrium (at which the orbital derivative is
zero). Thus it overcomes the difficulty of the symmetric collocation method used in Giesl [6].
Moreover, compared with the technique provided in Giesl [8], our scheme gives a both local and
global Lyapunov function directly and thus is simple and easy to compute.
Example 2
As the second example, we consider the system
(
x˙ = −x(1 − x2 (y − 3)2 )
y˙ = −(y − 3)(1 − x2 (y − 3)2 ),
which is a special case of (cf. [25], pp. 122)
(
x˙ = xf1 (x, y)
y˙ = yf2 (x, y)
(7)
(8)
modeling the populations of two species that interact and effect each other.
The system (7) has an asymptotically stable equilibrium (0, 3) whose basin of attraction is
A(0, 3) = {(x, y) ∈ R2 : x2 (y − 3)2 − 1 < 0}.
Besides, to demonstrate the benefits of the shape parameters, we choose three pairs of the two
shape parameters: a = 1, b = 4; a = b = 1; a = 4, b = 1, respectively.
Fig. 2 shows the boundaries of exact basin of attraction (the blue curves) and sublevel sets of
the reconstructed Lyapunov function with different shape parameters.
6182
W. Gao / Journal of Information & Computational Science 11:17 (2014) 6175–6183
Remark 2 From the figure, we can find that our scheme provides three different basins of attraction corresponding to the three pairs of shape parameters. In addition, by taking the union
of these three basins of attraction, we can get a large one that gives better approximation to the
exact basin of attraction. Such a technique has been provided in [16].
8
6
4
y
2
0
−2
−4
−6
−8
−8
−6
−4
−2
0
x
2
4
6
8
Fig. 2: Sublevel sets of the Lyapunov function with different shape parameters and the boundaries of
the exact basin of attraction
5
Conclusions and Discussions
With the help of multiquadric B-spline quasi-interpolation, a scheme for constructing Lyapunov
functions for autonomous dynamical systems is provided in the paper. The scheme overcomes
the difficulty of the symmetric collocation method discussed in Giesl [6]. Moreover, compared
with other techniques provided in Giesl’s works, i.e., [5, 8, 9], the scheme is simple, efficient
and easy to compute. In addition, one can get different basins of attraction for a dynamical
system by adjusting the shape parameters of the constructed Lyapunov function properly. More
importantly, by taking the union of these basins, one can even get a larger one with more accurate
approximation to the exact basin of attraction.
Acknowledgement
We wish to express our great gratitude to the referees for their valuable comments and suggestions.
References
[1]
[2]
H. Chiang, M. Hirsch, F. Wu, Stability basins of nonlinear autonomous dynamical systems [J],
IEEE Transaction of Automatic Control, 30 (1988), 16-27
R. Genesio, M. Tartaglia, A. Vicino, On the estimation of asymptotic stability regions: State of
the art and new proposals [J], IEEE Transaction of Automatic Control, 30 (1985), 747-755
W. Gao / Journal of Information & Computational Science 11:17 (2014) 6175–6183
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
6183
F. Camilli, L. Gr¨
une, F. Wirth, A generalization of Zubov’s method to perturbed systems [J],
SIAM, Journal of Control Optimization, 40 (2001), 496-515
S. Hafstein, A constructive converse Lyapunov theorem in exponential stability [J], Discret and
Continuous Dynamical System, 10 (2004), 657-678
P. Giesl, Construction of global Lyapunov functions usig radial basis functions [C], Lecture Notes
in Mathematics, Vol. 1904, Springer, Berlin, 2007
P. Giesl, Construction of a global Lyapunov function using radial basis functions with a single
operator [J], Discrete and Continuou Dynamical System - Series B (DCDS-B), 7 (2007), 101-124
Z. M. Wu, Hermite-Birkhoff Interpolation for scattered data by radial basis functions [J], Approximation Theory and Its Application, 8 (1992), 1-10
P. Giesl, Construction of a local and global Lyapunov function using radial basis functions [J],
IMA Journal of Applied Mathematics, 5 (2008), 1-21
P. Giesl, H. Wendland, Meshless collocation: Error estimates with application to dynamical systems [J], SIAM Journal of Numerical Analysis, 45 (2007), 1723-1741
Z. M. Wu, R. Schaback, Shape preserving properties and convergence of univariate multiquadric
quasi-interpolation [J], Acta Mathematical Application Sinca, 10 (1994), 441-446
R. Beatson, M. Powell, Univariate multiquadric approximation: Quasi-interpolation to scattered
data [J], Constrctive Approximation, 8 (1992), 275-288
C. de Boor, G. Fix, Spline approximation by quasi-interpolants [J], Journal of Approximation
Theroy, 8 (1973), 19-45
Y. C. Hon, A quasi-radial basis functions method for american options pricing [J], Computers and
Mathematics with Application, 43 (2002), 513-524
R. Q. Jia, J. J. Lei, A new version of Strang-Fix conditions [J], Journal of Approximation Theory,
74 (1993), 221-225
Z. W. Jiang, R. H. Wang, An improved numerical solution of Burgers’ equation by cubic B-spline
quasi-interpolation [J], Journal of Information and Computation Science, 5 (2010), 1013-1021
W. W. Gao, Z. M. Wu, Quasi-interpolation for linear functional data [J], Journal of Computation
and Applied Mathematics, 236 (2012), 3256-3264
G. Strang, G. Fix, A Fourier analysis of the finite-element method [C], in: G. Geymonat (ed.),
Constructive Aspects of Functional Analysis, C. I. M. E., Rome, 1973, 793-840
R. H. Wang, C. J. Li, Bivariate quartic spline spaces and quasi-interpolation operators [J], Journal
of Computational and Applied Mathematics, 190 (2006), 325-338
R. H. Wang, X. Xu, Q. Fang, A kind of improved univariate multiquadric quasi-interpolation
operators [J], Computers and Mathematics with Application, 59 (2010), 451-456
Z. M. Wu, L. M. Ma, Generator, multiquadric generator, quasi-interpolation and multiquadric
quasi-interpolation [J], Applied Mathematics, Journal of Chinese University, 26 (2011), 390-400
R. G. Yu, R. H. Wang, C. G. Zhu, A numerical method for solving KdV equation with multilevel
B-spline quasi-interpolation [J], Applied Analysis, 92 (2013), 1682-1690
T. Lyche, L. Schumaker, S. Stanley, Quasi-interpolants based on trigonometric splines [J], Journal
of Approximation Theory, 95 (1998), 280-309
W. W. Gao, Z. M. Wu, A quasi-interpolation scheme for periodic data based on trigonometric
B-splines [J], Journal of Computational and Applied Mathematics, 271 (2014), 20-30
W. W. Gao, Z. M. Wu, Approximation orders and shape-preserving properties of multiquadric
trigonometric B-spline quasi-interpolation [J], Journal of Mathematical Analysis and Application,
to appear
C. Robinson, An Introduction to Dynamical Systems: Continuous and Discrete [M], Englewood
Cliffs, NJ: Prentice-Hall, 2004