Theorem 1 (Bolzano-Weierstrass) Every bounded sequence in ℝ^{p} contains a convergent subsequence.

Proof We prove this theorem only for sequences in ℝ, that is, for the case when p = 1. Let a sequence a_{1}, a_{2}, …, a_{n}, …, denoted {a_{n}}, be bounded. Let Z be the set of all numbers x such that x ≤ a_{n} for infinitely many n. The set Z is non-empty. Indeed, the number −M belongs to the set Z, where |a_{n}| ≤ M, that is, the inequality −M ≤ a_{n} holds for all n. The set Z is bounded above, for if x ∈ Z, then x ≤ M because if x > M, then the inequality x ≤ a_{n} would not be satisfied by any a_{n}. Since the set Z is non-empty and bounded above, it has a least upper bound. Denote this bound by g. Thus

g = sup {x : x ≤ a_{n}} .

We now show that g is the limit of a certain subsequence of the sequence {a_{n}}. This means that we have to define a sequence of positive integers m_{1} < m_{2} < … < m_{n} < … such that

lim_{n→∞} a_{m_{n}} = g.

By Theorem ??, it follows that for any ε > 0 there exists infinitely many a_{n}’s such that

g − ε < a_{n} < g + ε.

Let ε = 1 in g − ε < a_{n} < g + ε. Then,

g − 1 < a_{n} < g + 1.

We choose one of n’s for which the above relation is satisfied and denote it by m_{1}. We thus have

g − 1 < a_{m_{1}} < g + 1.

Now let ε = 1/2. Then, there exist infinitely many a_{n}’s such that

g − \frac{1}{2} < a_{n} < g + \frac{1}{2}.

We choose one a_{n} for which the above relation hold and denote it by a_{m_{2}}, where m_{1} < m_{2}. In general, if m_{n} is defined, we choose m_{n+1} in such a way that

g − \frac{1}{n + 1} < a_{m_{n+1}} < g + \frac{1}{n + 1}.

Since lim_{n→∞} \frac{1}{n} = 0, we have

lim_{n→∞} (g − \frac{1}{n}) = g = lim_{n→∞} (g + \frac{1}{n}) .

The above two relations imply

lim_{n→∞} a_{m_{n}} = g.

The inequality m_{n} < m_{n+1}, on the other hand, implies that the sequence {a_{m_{n}}} is a subsequence of the sequence {a_{n}}. Thus, a bounded sequence contains a convergent subsequence, and the proof is complete.

□

1 The Inverse Function Theorem

Let S be an open set in ℝ^{n} and let f : S → ℝ^{n} be a map, where

f(x) = [ f_{1}(x) f_{2}(x) ⋯ f_{n}(x) ]^{T} .

Recall that if all partial derivatives of all functions f_{i} exist and are continuous in S, then we say that f is a C^{1} map in S. We say that f is a C^{p} map in S if each component f_{i} is continuous and has continuous partial derivatives in S for all k = 1, 2, …, p.

Let f : ℝ^{n} → ℝ^{n} be a C^{p}-map, p > 0, on an open subset S of ℝ^{n}, and let a be a point of S. We say that f is locally C^{p}-invertible at a if there is an open set U ⊂ S containing a and an open set V containing f(a) such that f : U → V has a C^{p} inverse h : V → U so that for all x ∈ U,

h(f(x)) = x.

Figure 1 illustrates the sets U and V for one variable case.

Figure 1: An illustration of the sets U and V that appear in the inverse function theorem.

The Jacobian matrix of f : ℝ^{n} → ℝ^{n}, denoted J(x), is defined as

J(x) = Df(x) = [ ∇^{T} f_{1}(x) ; ∇^{T} f_{2}(x) ; ⋮ ; ∇^{T} f_{n}(x) ] = [ Df_{1}(x) ; Df_{2}(x) ; ⋮ ; Df_{n}(x) ] .

A set S in ℝ^{n} is convex if for every x and y in S, the line segment joining x and y is in S. The line segment joining x and y, denoted [x, y], can be described as

[x, y] = {λx + (1 − λ)y : 0 ≤ λ ≤ 1}.  (1)

Therefore, a set S in ℝ^{n} is convex, if for every x and y in S and every λ ∈ [0, 1],

λx + (1 − λ)y ∈ S.

In the proof of the Inverse Function Theorem, we will use the following Mean Value Theorem for maps.

Theorem 2 Let S ⊆ ℝ^{n} be an open set and let f : S → ℝ^{q}. Suppose that S contains the points x, y and the line segment [x, y], and that f is differentiable at every point of [x, y]. Then, there exists a point ξ on [x, y] such that

∥f(y) − f(x)∥ ≤ ∥Df(ξ)(y − x)∥  (2)

Proof If f(y) − f(x) = 0 in ℝ^{q}, then (2) holds and we are done. Suppose that f(y) − f(x) ≠ 0. Let

z = (f(y) − f(x)) / ∥f(y) − f(x)∥ .

Note that ∥z∥ = 1. Define

g(x) = f^{T} (x)z.  (3)

We have

g(y) − g(x) = (f(y) − f(x))^{T} z = ∥f(y) − f(x)∥ .  (4)

Using the Mean Value Theorem for real-valued functions, we obtain

g(y) − g(x) = Dg(ξ) (y − x)  (5)

for some ξ ∈ [x, y], where

Dg(ξ) = z^{T} Df(ξ).  (6)

Substituting (6) and (4) into (5) yields

∥f(y) − f(x)∥ = z^{T} Df(ξ) (y − x) .  (7)

Using Schwarz’s Inequality and the fact that ∥z∥ = 1 gives

∥f(y) − f(x)∥ ≤ ‖z^{T} ‖ ∥Df(ξ) (y − x)∥ ≤ ∥Df(ξ) (y − x)∥ ,

which completes the proof of the theorem.

□

We now state and prove the Inverse Function Theorem for maps.

Theorem 3 Let f : ℝ^{n} → ℝ^{n} be a C^{p}-map, p > 0, on an open subset S of ℝ^{n}, and let a be a point of S such that the Jacobian matrix J(a) is invertible, and b = f (a). Then, there exist open sets U and V in ℝ^{n} such that a ∈ U, b ∈ V, f is one-to-one on U, f(U) = V, and f is locally C^{p}-invertible at a.

Proof We first show that there exists an open subset U ⊆ S in which f is one-to-one. Since J^{−1} (a) exists, we can choose s so that

‖J^{−1} (a)‖ = \frac{1}{2s}.  (8)

By assumption Df is continuous at a. Hence, there is a ball U ⊆ S centered at a such that

∥J (x) − J (a)∥ < s, x ∈ U.  (9)

For each y ∈ ℝ^{n}, we associate

F (x) = x + J^{−1} (a) (y − f(x)) .  (10)

We have f(x) = y if and only if x is a fixed point of the map F given above. The Jacobian matrix of F is

DF (x) = I_{n} − J^{−1} (a) J(x) = J^{−1} (a) (J (a) − J(x)) .  (11)

Applying (8) and (9) to the above gives

∥DF (x)∥ < \frac{1}{2}, x ∈ U.  (12)

By the Mean Value Theorem for maps, stated on page 3, for any x^{(1)}, x^{(2)} ∈ U there exists ξ ∈ [x^{(1)}, x^{(2)}] such that

‖F (x^{(2)}) − F (x^{(1)})‖ ≤ ‖DF (ξ) (x^{(2)} − x^{(1)})‖ .  (13)

The point ξ is in U because the set U, which is a sphere, is convex. Using (12) gives

‖F (x^{(2)}) − F (x^{(1)})‖ ≤ ∥DF (ξ)∥ ‖x^{(2)} − x^{(1)}‖
≤ \frac{1}{2} ‖x^{(2)} − x^{(1)}‖ .  (14)

Therefore, F is a contraction map and can have at most one fixed point in U for a given y, which means that f(x) = y for at most one x ∈ U. This proves that f is one-to-one in U.

Let V = f(U). We will show that for any y^{(0)} ∈ V there is a neighborhood contained in V, which means that the set V open. If y^{(0)} ∈ V, then there is x^{(0)} ∈ U such that y^{(0)} = f (x^{(0)}). Denote by B an open ball with center at x^{(0)} and radius r > 0. Choose r so small that the closed ball B̄ is in U, that is,

B̄ = {x : ‖x − x^{(0)}‖ ≤ r} ⊆ U.

We will now show that for any y^{(0)} ∈ V there is an open ball with center at y^{(0)} and radius rs that is contained in V. For a given y^{(0)} fix y such that

‖y − y^{(0)}‖ < rs.

Then, using (8), we obtain

‖F (x) − x^{(0)}‖ = ‖J^{−1} (a) (y − y^{(0)})‖
≤ ‖J^{−1} (a)‖ rs
= \frac{r}{2}.  (15)

If x is in B̄, then using (14) and (15), we get

‖F (x) − x^{(0)}‖ ≤ ‖F (x) − F (x^{(0)})‖ + ‖F (x^{(0)}) − x^{(0)}‖
< \frac{1}{2} ‖x − x^{(0)}‖ + \frac{r}{2}
≤ r.  (16)

Hence, if x ∈ B̄, then F (x) ∈ B, which means that F is a contraction of B̄ into B̄. Since B̄ is a closed subset of ℝ^{n}, it is complete. By Theorem ??, the map F has a fixed point x ∈ B̄. Furthermore, y = f(x). Therefore,

y ∈ f \binom{}{B̄} ⊆ f(U) = V.

This completes the proof of the fact that the set V is open. In summary, an inverse function h : V → U exists such that x_{i} = h_{i}(y) for y ∈ V.

We will now prove that the inverse function h is a C^{1} map. Applying the Mean Value Theorem to the i-th component f_{i} : ℝ^{n} → ℝ of f, we obtain

f_{i} (x^{(2)}) − f_{i} (x^{(1)}) = Df_{i} (ξ^{(i)}) (x^{(2)} − x^{(1)}) , i = 1, 2, …, n,

where Df_{i} = [ \frac{∂f_{i}}{∂x_{1}} \frac{∂f_{i}}{∂x_{2}} ⋯ \frac{∂f_{i}}{∂x_{n}} ]. The above equations can be represented in matrix form as

f (x^{(2)}) − f (x^{(1)}) = [ Df_{1} (ξ^{(1)}) ; Df_{2} (ξ^{(2)}) ; ⋮ ; Df_{n} (ξ^{(n)}) ] (x^{(2)} − x^{(1)}) ,  (17)

By assumption, the Jacobian matrix of f is invertible at a, that is,

det Df (a) = det [ Df_{1} (a) ; Df_{2} (a) ; ⋮ ; Df_{n} (a) ] ≠ 0.

Since Df is continuous, there is a subset U of S within which det Df(x) ≠ 0 and

det [ Df_{1} (ξ^{(1)}) ; Df_{2} (ξ^{(2)}) ; ⋮ ; Df_{n} (ξ^{(n)}) ] ≠ 0.  (18)

Therefore, we can represent (17) as

h (y^{(2)}) − h (y^{(1)}) = [ Df_{1} (ξ^{(1)}) ; Df_{2} (ξ^{(2)}) ; ⋮ ; Df_{n} (ξ^{(n)}) ]^{−1} (y^{(2)} − y^{(1)}) .  (19)

Let

[ Df_{1} (ξ^{(1)}) ; Df_{2} (ξ^{(2)}) ; ⋮ ; Df_{n} (ξ^{(n)}) ]^{−1} = [ D^{−1}f_{1} (ξ^{(1)}) ; D^{−1}f_{2} (ξ^{(2)}) ; ⋮ ; D^{−1}f_{n} (ξ^{(n)}) ] ,

then we can write

h_{i} (y^{(2)}) − h_{i} (y^{(1)}) = D^{−1}f_{i}(ξ^{(i)}) (y^{(2)} − y^{(1)}) , i = 1, 2, …, n.  (20)

Consider the limit

lim_{y^{(2)}→y^{(1)}} \frac{|h_{i} (y^{(2)}) − h_{i} (y^{(1)}) − D^{−1}f_{i}(x^{(1)}) (y^{(2)} − y^{(1)})|}{∥y^{(2)} − y^{(1)}∥},  (21)

where D^{−1}f_{i}(x^{(1)}) is the ith row of D^{−1}f (x^{(1)}). Substituting (20) into (21) gives

lim_{y^{(2)}→y^{(1)}} \frac{|(D^{−1}f_{i} (ξ^{(1)}) − D^{−1}f_{i} (x^{(1)})) (y^{(2)} − y^{(1)})|}{∥y^{(2)} − y^{(1)}∥} = 0

because

lim_{x^{(2)}→x^{(1)}} D^{−1}f_{i} (ξ^{(1)}) = D^{−1}f_{i} (x^{(1)}) .

Thus,

Dh_{i}(y) = D^{−1}f_{i}(h(y)),

which is continuous, and so h is of class C^{1}. Furthermore, if f is of class C^{p}, so is h. □

1.1 The Implicit Function Theorem

The implicit function theorem is closely related to the inverse function theorem. We illustrate the concept behind the implicit function theorem on a function on ℝ^{2}. Let U be an open set in ℝ^{2}, and let

f : U → ℝ

be a C^{1}-function. Let [ x_{0} y_{0} ]^{T} be a point in U, and let

f(x_{0}, y_{0}) = 0.

Note that the right-hand side of the above equation can be an arbitrary constant. Without loss of generality, we can assume that this constant is 0, as in the above equation. Our goal will be to determine whether there is a differentiable function y = ϕ(x) defined near x = x_{0} such that ϕ(x_{0}) = y_{0} and

f(x, ϕ(x)) = 0

for all x near x_{0}. If such a function ϕ exists, then we say that y = ϕ(x) is the function determined implicitly by f.

As another illustrating example, consider a system of linear equations written as

[ D B ] [ \begin{matrix}x \\ y \end{matrix} ] − b = 0,

where D ∈ ℝ^{m×n}, B ∈ ℝ^{m×m}, x ∈ ℝ^{n}, and y ∈ ℝ^{m}. A unique solution for y exists if det B ≠ 0. We now generalize this result to a nonlinear case.

Theorem 4 Let

f : ℝ^{n+m} → ℝ^{m}

be a C^{p}, p ≥ 1, function, where the ith component of f is

f_{i}(x_{1}, x_{2}, …, x_{n}, y_{1}, y_{2}, …, y_{m}), i = 1, 2, …, m.

Let x_{0} ∈ ℝ^{n} and y_{0} ∈ ℝ^{m}, and let f (x_{0}, y_{0}) = 0. Assume that at (x_{0}, y_{0})

det \frac{∂f}{∂y} ≠ 0.

Then, there exists an open set U ⊂ ℝ^{n} containing x_{0} and open set V ⊂ ℝ^{m} containing y_{0} such that there exists a C^{p} function

ϕ : U → V, y = ϕ(x),

with ϕ(x_{0}) = y_{0} and f(x, ϕ(x)) = 0 for all x in U.

Proof Consider the following system of equations:

z_{1} = x_{1}
⋮
z_{n} = x_{n}
z_{n+1} = f_{1}(x, y)
⋮
z_{n+m} = f_{m}(x, y),

where

z_{10} = x_{10}
⋮
z_{n0} = x_{n0}
z_{n+1,0} = f_{1}(x_{0}, y_{0}) = 0
⋮
z_{n+m,0} = f_{m}(x_{0}, y_{0}) = 0.

The above system defines a map F : ℝ^{n+m} → ℝ^{n+m}. We claim that F is locally C^{p}-invertible at [ x_{0}^{T} , y_{0}^{T} ]^{T} ∈ ℝ^{n+m}. We compute the Jacobian matrix of F. We have

DF = [ I_{n} 0_{\frac{∂f}{∂x} \frac{∂f}{∂y}} ] .

The determinant of the above Jacobian matrix evaluated at (x_{0}, y_{0}) is nonzero because

det DF = det \frac{∂f}{∂y} .

Hence, the inverse mapping theorem applies, and we have the inverse system

x_{1} = z_{1}
⋮

x_{n} = z_{n}
y_{1} = h_{1}(z_{1}, z_{2}, …, z_{n+m})
⋮
y_{m} = h_{m}(z_{1}, z_{2}, …, z_{n+m}).

such that

z_{n+1} = f_{1}(x, h_{1}(x, z_{n+1}, …, z_{n+m}), …, h_{m}(x, z_{n+1}, …, z_{n+m}))
⋮
z_{n+m} = f_{m}(x, h_{1}(x, z_{n+1}, …, z_{n+m}), …, h_{m}(x, z_{n+1}, …, z_{n+m}))

for all z_{n+1}, z_{n+2}, …, z_{n+m}. In particular, for

z_{n+1} = z_{n+2} = … = z_{n+m} = 0

we have

f (x, h(x, 0)) = 0,

and

h (x_{0}, 0) = y_{0}.  (22)

Thus,

y = ϕ(x) = h(x, 0).

Note that so defined function ϕ is indeed the required function because by (22),

y_{0} = h(x_{0}, 0) = ϕ(x_{0}).

This proves the implicit function theorem.

□

=============================================

We now discuss a numerical approach to computing the largest Lyapunov exponent. One approach is to numerically integrate (??) to obtain a reference trajectory ϕ_{t}(x_{0}), starting from a given initial point x_{0}. Then, at each time step t_{k} integrate (??) using the initial condition ϕ_{t_{k}}(x_{0}) + εη arbitrarily close to ϕ_{t_{k}}(x_{0}), where ∥η∥ = 1. Hence,

ϕ_{t_{k+1}} (ϕ_{t_{k}}(x_{0}) + εη) ≈ ϕ_{t_{k+1}}(x_{0}) + ε Φ_{t_{k+1}}(x_{0})η.

Let d_{0}(t_{k−1}) = ε and d(t_{k}) = ε ∥Φ_{t_{k+1}}(x_{0})η∥. Then, we could use (??) to evaluate the largest Lyapunov exponent.

We now describe a more efficient numerical method for evaluating the largest Lyapunov exponent, that can be found, for example, in the paper by Parker and Chua [4]. First, select the time interval ΔT and iteration count N. Then, integrate (??) subject to a given initial condition x(t_{0}) = x_{0} to obtain a reference trajectory ϕ_{t}(x_{0}) of the nonlinear system over the time interval [t_{0}, NΔT]. Select an arbitrary Δx_{0} such that ∥Δx_{0}∥ = 1. Find the solution Δx(t_{1}) = Δx_{1} to

\frac{d}{dt}Δx = Df(ϕ_{t}(x_{0}))Δx  (23)

subject to the initial condition Δx_{0}. Observe that we can set

d_{0}(t_{0}) = ∥Δx_{0}∥ = 1,

and

d(t_{1}) = ∥Φ_{t_{1}}(x_{0})Δx_{0}∥ = ∥Δx_{1}∥ .

Next, normalize the solution Δx_{1} and use this normalized solution as an initial condition to solve the differential equation (23) again. Repeat the process until t_{N} = NΔT. Denote Δx(t_{N}) = Δx_{N}. Then,

Δx_{N} = Φ_{t_{N}} (x_{N−1})Δx_{N−1}/∥Δx_{N−1}∥
= Φ_{t_{N}} (x_{N−1})Φ_{t_{N−1}}(x_{N−2})Δx_{N−2}/ (∥Δx_{N−1}∥∥Δx_{N−2}∥)
⋮
= Φ_{t_{N}} (x_{N−1})Φ_{t_{N−1}}(x_{N−2})
⋯ Φ_{t_{1}}(x_{0})Δx_{0}/ (∥Δx_{N−1}∥∥Δx_{N−2}∥ ⋯ ∥Δx_{1}∥) ,  (24)

Observe that we can set d_{0}(t_{k−1}) = 1 and d(t_{k}) = ∥Δx_{k}∥. We also have

Φ_{t_{N}} (x_{0})Δx_{0} = Φ_{t_{N}} (x_{N−1})Φ_{t_{N−1}}(x_{N−2}) ⋯ Φ_{t_{1}}(x_{0})Δx_{0}.  (25)

Comparing (24) with (25) gives

Φ_{t_{N}} (x_{0})Δx_{0} = Δx_{N}∥Δx_{N−1}∥∥Δx_{N−2}∥ ⋯ ∥Δx_{1}∥.  (26)

Therefore, for large N,

μ_{1} ≈ \frac{1}{NΔT} ln ∥Φ_{t_{N}} (x_{0})Δx_{0}∥
= \frac{1}{NΔT} ln ∏_{k=1}^{N} ∥Δx_{k}∥
= \frac{1}{NΔT} ∑_{k=1}^{N} ln ∥Δx_{k}∥  (27)
= \frac{1}{NΔT} ∑_{k=1}^{N} ln \frac{d(t_{k})}{d_{0}(t_{k−1})},  (28)

which is the defining formula (??).

======================

LaSalle’s Invariance Principle

In this section, we generalize the Lyapunov function concept using LaSalle’s approach as reported in [3]. This generalization is based on the idea of the so called invariant set. We first present this generalization for nonlinear dynamical systems modeled by the difference equation

x(k + 1) = F (x(k)),  (29)

where k ∈ Z_{+}, Z_{+} denotes the set of nonnegative integers, x : Z_{+} → ℝ^{n}, that is, x(k) is the state of the system model (29) at time k, and F : ℝ^{n} → ℝ^{n}. We assume that the vector-valued function F is continuous on ℝ^{n}, that is, if x^{j} → y then F (x^{j}) → F (y). If x is a vector, then F (x) is a vector and is the value of the function F at x.

The solution to the initial value problem

x(k + 1) = F (x(k)), x(0) = x^{0}

is

x(k) = x^{k} = F^{k} (x^{0}) ,

where F^{k} is the k-th iterate of F, that is, F^{0} = I_{n}, and F^{k} = F F^{k−1}.

A state x_{e} ∈ ℝ^{n} is a fixed point or an equilibrium state of the system x(k + 1) = F (x(k)) if F (x_{e}) = x_{e}. Thus the solution of the system that starts at x_{e} stays at x_{e} for all k ≥ 0.

A subset A ⊂ ℝ^{n} is called positively invariant if F (A) ⊂ A, that is, if x ∈ A then F (x) ∈ A. It is clear that an equilibrium solution is invariant. The adjective “positive” is added to emphasize the fact that we will be interested in the behavior of the system for k ≥ 0.

The solution through x^{0} is bounded if there is a constant M such that ∥F^{k} (x^{0}) ∥ ≤ M for all k ≥ 0.

Consider again a particular solution F^{k} (x^{0}). A point y ∈ ℝ^{n} is a positive limit point or a positive point of accumulation of F^{k} (x^{0}) if there is a subsequence k_{i} such that F^{k_{i}} (x^{0}) → y. The set of all positive limit points of F^{k} (x^{0}) is called the positive limit set and denoted Ω (x^{0}).

A sequence {x^{k}} is said to approach a set A if

ρ (x^{k}, A) = inf {∥x^{k} − y∥ : y ∈ A} → 0.

Lemma 1 If F is continuous and F^{k} (x^{0}) is bounded for all k ≥ 0, then Ω (x^{0}) is nonempty, compact positively invariant set. Furthermore, F^{k} (x^{0}) → Ω (x^{0}) as k → ∞. Proof Because the sequence {x^{k}} is bounded, by the Bolzano-Weierstrass theorem, stated on page 1, there is at least one positive limit point of {x^{k}} and so Ω (x^{0}) is non-empty and bounded. For a proof that the set Ω (x^{0}) is closed, we refer to Khalil [page 652][2]. Because the set Ω (x^{0}) is also bounded, it is therefore compact.

We now show that Ω (x^{0}) is positively invariant. For this let y ∈ Ω (x^{0}). Hence, there is a subsequence k_{i} such that x^{k_{i}} → y. Because F is continuous, x^{k_{i}+1} = F (x^{k_{i}}) → F (y) as i → ∞. Therefore F (y) ∈ Ω (x^{0}) which implies that Ω (x^{0}) is invariant. Note that we do not have to have y = F (y).

Finally, we will show that F^{k} (x^{0}) → Ω (x^{0}) as k → ∞. We shall prove the statement by contradiction. First observe that ρ (x^{k}, Ω (x^{0})) is bounded because both {x^{k}} and Ω (x^{0}) are. Assume that ρ (x^{k}, Ω (x^{0})) does not converge to zero. Then, there is a subsequence k_{i} such that x^{k_{i}} → y and ρ (x^{k_{i}}, Ω (x^{0})) → a > 0. But since y ∈ Ω (x^{0}),

ρ (x^{k_{i}}, Ω (x^{0})) ≤ ∥x^{k_{i}} − y∥ → 0,

and thus ρ (x^{k_{i}}, Ω (x^{0})) → 0 which is a contradiction. The proof of the lemma is now complete.

□

Let V : ℝ^{n} → ℝ. Define the Lyapunov difference of V relative to x(k + 1) = F (x(k)) as

ΔV (x) = V (F (x)) − V (x),

that is,

ΔV (x(k)) = V (x(k + 1)) − V (x(k)).

If ΔV (x(k)) ≤ 0, then V is non-increasing along solutions of the system.

We say that V is a Lyapunov function for x(k + 1) = F (x(k)) on a set G ⊆ ℝ^{n} if 1. V is continuous on ℝ^{n}, and

2. ΔV (x(k)) ≤ 0 for all x(k) ∈ G.

A c-level surface of V is the set

V^{−1}(c) = {x : V (x) = c, x ∈ ℝ^{n}}.

Theorem 5 (LaSalle’s Invariance Principle) If

1. V is a Lyapunov function for x(k + 1) = F (x(k)) in G, and

2. the solution F^{k} (x^{0}) is bounded and in G,

then, there is a number c such that F^{k} (x^{0}) → M∩V^{−1}(c), where M is the largest positively invariant set contained in the set E = {x ∈ ℝ^{n} : ΔV = 0} ∩ Ḡ, where Ḡ denotes the closure of G.

Proof Because x^{k} = F^{k} (x^{0}) is bounded and in G, by Lemma 1 the positive limit set Ω (x^{0}) is non-empty and in G. Furthermore, x^{k} tends to Ω (x^{0}). Now {V (x^{k})} is nonincreasing and bounded below, so V (x^{k}) → c for some c. If y ∈ Ω (x^{0}), then there is a subsequence k_{i} such that x^{k_{i}} → y, and hence V (x^{k_{i}}) → V (y), that is, V (y) = c. Thus V (Ω (x^{0})) = c, and hence Ω (x^{0}) ⊂ V^{−1}(c). Because V (Ω (x^{0})) = c and Ω (x^{0}) is positively invariant, we have ΔV (Ω (x^{0})) = 0. In summary

x^{k} → Ω (x^{0}) ⊂ {x ∈ ℝ^{n} : ΔV = 0} ∩ Ḡ ∩ V^{−1}(c).

Since Ω (x^{0}) is positively invariant we have Ω (x^{0}) ⊂ M, and the proof of the theorem is complete.

□

Using the preceding theorem we can study stability properties of difference equations.

LaSalle’s Invariance Principle can be used to prove a result concerning the instability of equilibrium states.

Theorem 6 Let x_{e} be an equilibrium state of x(k + 1) = F (x(k)). Let ΔV be positive definite with respect to x_{e}, and let V take positive values arbitrarily close to x_{e}, and V (x_{e}) = 0, then x_{e} is unstable.

Proof We prove the theorem by contradiction. First observe that −V is a Lyapunov function because Δ(−V) < 0. Since x^{k} = F^{k} (x^{0}) is bounded we can apply the Invariance Principle to conclude that x^{k} tends to {x : Δ(−V (x)) = 0} = {x_{e}}, that is, x^{k} → x_{e}. Using the above, we conclude that V (x^{k}) → V (x_{e}) = 0. However, ΔV (x^{k}) > 0 and so

V (x^{k}) ≥ V (x^{k−1}) ≥ ⋯ ≥ V (x^{0}) > 0.

Hence, V (x^{k}) > 0 and thus we cannot have V (x_{e}) = 0. The above contradiction completes the proof.

□

1.2 Properties of M-Matrices

We start with the definition of the M-matrix.

Definition 1 A real n × n matrix A is an M-matrix if all its off-diagonal elements are non-positive, that is, a_{ij} ≤ 0, i ≠ j, and all its principal minors are positive.

Note that if A is an M-matrix, then its diagonal elements are positive since its principal minors are positive. In establishing the needed properties of the M-matrices, we use the following technical lemma due to Fiedler and Pták [1].

Lemma 2 The following properties of an arbitrary matrix A ∈ ℝ^{n×n} are equivalent:

(i) The leading principal minors of A are positive,

(ii) There exists a lower triangular matrix L ∈ ℝ^{n×n} and an upper triangular matrix U ∈ ℝ^{n×n} both with positive diagonal elements such that A = LU

Proof We first prove the implication (i)⇒(ii). The proof is by induction. Obviously the implication holds for n = 1. Suppose the statement holds for all real matrices of order n − 1. Let us represent A ∈ ℝ^{n×n} as follows:

A = [ \begin{matrix}A_{n−1} & a \\ b^{T} & a_{nn} \end{matrix} ] ,

where A_{n−1} ∈ ℝ^{(n−1)×(n−1)}, a ∈ ℝ^{(n−1)×1}, b^{T} ∈ ℝ^{1×(n−1)}, and a_{nn} ∈ ℝ. From the inductive hypothesis, the matrix A_{n−1} has positive leading principal minors, A_{n−1} = L_{n−1}U_{n−1}, where L_{n−1} ∈ ℝ^{(n−1)×(n−1)} is a lower triangular matrix with positive diagonal elements, and U_{n−1} ∈ ℝ^{(n−1)×(n−1)} is an upper triangular matrix with positive diagonal elements. Observe that

[ \begin{matrix}I_{n−1} & 0 \\ −b^{T} A_{n−1}^{−1} & 1 \end{matrix} ] A = [ \begin{matrix}A_{n−1} & a \\ 0 & a_{nn} − b^{T} A_{n−1}^{−1}a \end{matrix} ] ,

where from the assumptions det A_{n−1} > 0, and det A > 0. Hence,

0 < det A = (a_{nn} − b^{T} A_{n−1}^{−1}a) det A_{n−1}.

Thus, we have

a_{nn} − b^{T} A_{n−1}^{−1}a > 0.

Let

L = [ \begin{matrix}L_{n−1} & 0 \\ b^{T} U_{n−1}^{−1} & 1 \end{matrix} ] ,

and

U = [ \begin{matrix}U_{n−1} & L_{n−1}^{−1}a \\ 0 & a_{nn} − b^{T} A_{n−1}^{−1}a \end{matrix} ] .

One can check that

A = LU,

and that the matrices L and U satisfy (ii).

We now prove the implication (i)⇐(ii). From the generalized Cauchy’s Theorem it follows that the leading principal minors of A = LU are just the products of the corresponding leading principal minors of L and U. This is because the matrix L is lower triangular with positive diagonal elements and the matrix U is upper triangular with positive diagonal elements. This completes our proof of the lemma.

□

The following theorem, also due to Fiedler and Pták [1], establishes some properties of the M-matrices in terms of their LU-decomposition.

Theorem 7 Let A ∈ ℝ^{n×n} be an M-matrix. Then, there is a lower triangular matrix L ∈ ℝ^{n×n} and an upper triangular matrix U ∈ ℝ^{n×n} both with positive diagonal elements. Their off-diagonal elements are non-positive. Furthermore, both L^{−1} and U^{−1} are nonnegative which implies that A^{−1} = U^{−1}L^{−1} is non-negative.

Proof We prove this theorem by induction. The theorem holds for n = 1. Suppose the theorem holds for all matrices of order n − 1. Consider now the factorization from the proof of Lemma 2. We have

A = LU = [ \begin{matrix}L_{n−1} & 0 \\ b^{T} U_{n−1}^{−1} & 1 \end{matrix} ] [ \begin{matrix}U_{n−1} & L_{n−1}^{−1}a \\ 0 & a_{nn} − b^{T} A_{n−1}^{−1}a \end{matrix} ] ,

where L_{n−1}, U_{n−1}, L_{n−1}^{−1}, U_{n−1}^{−1} and A_{n−1}^{−1} = U_{n−1}^{−1}L_{n−1}^{−1} satisfy the inductive hypotheses. We see that l_{ij} ≤ 0 for i > j and u_{kl} ≤ 0 for k < l which means that the off-diagonal elements of L and U are non-positive. Their diagonal elements are positive. The inverse of L is

L^{−1} = [ \begin{matrix}L_{n−1}^{−1} & 0 \\ −b^{T} L_{n−1}^{−1} & 1 \end{matrix} ] .

From the inductive assumption L_{n−1}^{−1} and A_{n−1}^{−1} are non-negative and hence L_{n−1}^{−1} is nonnegative. The inverse of U is

U^{−1} = [ \begin{matrix}U_{n−1}^{−1} & −A_{n−1}^{−1} \frac{a}{a_{nn}−b^{T}A_{n−1}^{−1}a} \\ 0 & (a_{nn} − b^{T} A_{n−1}^{−1}a)^{−1} \end{matrix} ] .

By the inductive assumptions and the fact that a_{nn} − b^{T} A_{n−1}^{−1}a > 0, the matrix U^{−1} is non-negative, and the proof is complete.

□

We now consider a more general class of matrices called Z-matrices.

Definition 2 A real n × n matrix A whose off-diagonal elements are non-positive is called a Z-matrix.

Note that M-matrices are Z-matrices. To simplify our presentation, we use the notation B > 0, to mean that all the elements of B are positive, where B can be a matrix or a vector. Respectively, B ≥ 0 means that all the elements of B are non-negative. We are now ready to state and prove the following theorem concerning properties of Z-matrices. The theorem can be found in Fiedler and Pták [1].

Theorem 8 Let A ∈ ℝ^{n×n} be a Z-matrix. Then, the following conditions are equivalent:

(i) There exists a vector x ≥ 0 such that Ax > 0; (ii) There exists a vector y > 0 such that Ay > 0;

(iii) There exists a diagonal matrix D with positive diagonal elements such that ADe > 0, where e = [ 1, …, 1 ]^{T} ∈ ℝ^{n};

(iv) There exists a diagonal matrix D with positive diagonal elements such that the matrix W = AD has a dominant positive main diagonal, that is for i = 1, 2, …, n, we have

w_{ii} > ∑n ; j≠i^{j=1} w_{ij}.

Proof We prove this theorem using the following scheme:

(i) ⇒ (ii) ⇒ (iii) ⇒ (iv) ⇒ (i).

We begin with the implication (i)⇒(ii). Thus suppose that there exists a vector x ≥ 0 such that Ax > 0. Let e = [ 1, …, 1 ]^{T} ∈ ℝ^{n}, ϵ > 0, and y = x + ϵe > 0. Then

Ay = A(x + ϵe) = Ax + ϵAe > 0

for a sufficiently small ϵ > 0.

We next prove that (ii)⇒(iii). For this let y = [ y_{1}, …, y_{n} ]^{T} > 0 be such that Ay > 0. Let

D = [ \begin{matrix}y_{1} & 0 & 0 & … & 0 \\ 0 & y_{2} & 0 & … & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & 0 & … & y_{n} \end{matrix} ] .

It then follows that ADe = Ay > 0.

To prove (iii)⇒(iv), let D > 0 be a diagonal matrix with positive diagonal elements such that ADe > 0. We denote

AD = W = [ \begin{matrix}w_{11} & −w_{12} & … & −w_{1n} \\ −w_{21} & w_{22} & … & −w_{2n} \\ ⋮ & ⋮ & ... & ⋮ \\ −w_{n1} & −w_{n2} & … & w_{nn} \end{matrix} ] ,

where w_{ij} ≥ 0 for i ≠ j. Note that ADe = W e > 0 which is equivalent to

w_{ii} > ∑_{j≠i^{j=1}}^{n} w_{ij}, i = 1, 2, …, n.

Since w_{ij} ≥ 0, we have

w_{ii} > ∑_{j≠i^{j=1}}^{n} |w_{ij}|.

The above means that W is a row diagonal dominant matrix with positive main diagonal. Note a_{ii} > 0 since 0 < w_{ii} = a_{ii}d_{i} and d_{i} > 0.

To complete the proof of the theorem we now show that (iv)⇒(i). For this let D be a diagonal matrix with positive diagonal elements such that W = AD has a positive dominant main diagonal. Then

0 < W e = ADe = Ax,

where

x = [ \begin{matrix}d_{1} & … & 0 \\ ⋮ & ... & ⋮ \\ 0 & … & d_{n} \end{matrix} ] [ 1 ⋮ 1 ] = [ d_{1} ; ⋮ ; d_{n} ] > 0.

Hence (i) is satisfied.

□

To proceed further, we need the following two technical results.

Lemma 3 An n × n matrix A with complex coefficients such that

|a_{ii}| > ∑_{j≠i^{j=1}}^{n} |a_{ij}| for i = 1, 2, …, n,

is nonsingular, that is, det A ≠ 0.

Proof We prove the above lemma by contraposition. Suppose that det A = 0. Then, there exists a nontrivial solution x ≠ 0 to Ax = 0, that is for k = 1, 2, …, n,

∑_{j=1}^{n} a_{kj}x_{j} = 0.

Let

x_{i} = max_{1≤k≤n} |x_{k}|.

Rearranging ∑_{j=1}^{n} a_{ij}x_{j} = 0 yields

−a_{ii}x_{i} = a_{i1}x_{1} + ⋯ + a_{i,i−1}x_{i−1} + a_{i,i+1}x_{i+1} + ⋯ + a_{in}x_{n}.

Taking the absolute values of the both sides of the above equation gives

|a_{ii}||x_{i}| = |a_{i1}x_{1} + ⋯ + a_{i,i−1}x_{i−1} + a_{i,i+1}x_{i+1} + ⋯ + a_{in}x_{n}| ≤ ∑_{j≠i^{j=1}}^{n} |a_{ij}||x_{j}|.

Thus

|a_{ii}| ≤ ∑_{j≠i^{j=1}}^{n} |a_{ij}|\frac{|x_{j}|}{|x_{i}|} ≤ ∑_{j≠i^{j=1}}^{n} |a_{ij}|

which contradicts the assumption that |a_{ii}| > ∑n ; j≠i^{j=1} |a_{ij}|.

□

The above lemma says that if a matrix A is row diagonal dominant then it is nonsingular. Following the arguments of the proof of the lemma, we can show that if a square, real or complex, matrix is column diagonal dominant then it is also nonsingular. We state this result in the following corollary.

Corollary 1 A real or complex n × n matrix A such that for j = 1, 2, …, n,

|a_{jj}| > ∑_{k≠j^{k=1}}^{n} |a_{kj}|

is nonsingular.

We now use Lemma 3 to prove the following theorem due to Gershgorin.

Theorem 9 Given a real or complex n × n matrix A. Let D_{i}, i = 1, 2, …, n, denotes the i-th closed disk in the complex plane centered at a_{ii} whose radius is r_{i} = ∑n ; j≠i^{j=1} |a_{ij}|. Then the eigenvalues of A are located in the union ⋃_{i=1}^{n} D_{i}.

Proof The eigenvalues of A are the roots of the characteristic equation

det[λI − A] = det [ \begin{matrix}λ − a_{11} & −a_{12} & … & −a_{1n} \\ −a_{21} & λ − a_{22} & … & −a_{2n} \\ ⋮ & ⋮ & ... & ⋮ \\ −a_{n1} & −a_{n2} & … & λ − a_{nn} \end{matrix} ] = 0.

By Lemma 3, the characteristic polynomial is nonzero if for all i = 1, 2, …, n,

|λ − a_{ii}| > ∑_{j≠i^{j=1}}^{n} |a_{ij}|,

that is, for all i = 1, 2, …, n

|λ − a_{ii}| > r_{i}.

Thus the eigenvalues of A must be in the set ⋃_{i=1}^{n} D_{i}.

□

In the following theorem, due to Fiedler and Pták [1], we collect a number of equivalent conditions characterizing the M-matrices. Theorem 10 Let A be a Z-matrix. Then the following conditions are equivalent

(i) There exists a vector x ≥ 0 such that Ax > 0; (ii) There exists a diagonal matrix D with positive diagonal such that the matrix W = AD

has a positive main diagonal and is row diagonal dominant; (iii) Eigenvalues of W have positive real parts and principal minors of W are positive; (iv) All principal minors of A are positive, that is, A is an M-matrix;

(v) All leading principal minors of A are positive; (vi) There exists a lower triangular matrix L and an upper triangular matrix U both with

positive diagonals such that A = LU. Furthermore, l_{ij} ≤ 0 for i > j, and u_{kl} ≤ 0 for k < l; (vii) A^{−1} exists and A^{−1} ≥ 0. Proof The implication (i)⇒(ii) follows from Theorem 8. To prove (ii)⇒(iii), we represent the matrix W = AD as follows

W = AD = [ \begin{matrix}w_{11} & −w_{12} & … & −w_{1n} \\ −w_{21} & w_{22} & … & −w_{2n} \\ ⋮ & ⋮ & ... & ⋮ \\ −w_{n1} & −w_{n2} & … & w_{nn} \end{matrix} ] ,

where w_{ij} ≥ 0 for i ≠ j, and w_{ii} > 0, as follows from Theorem 8. Furthermore

w_{ii} > ∑_{j≠i^{j=1}}^{n} w_{ij} for i = 1, 2, …, n.

The above implies that each principal minor matrix of W is also row diagonal dominant. By the Theorem of Gershgorin, the eigenvalues of W as well as its principal minor matrices are located on the disks with centers at w_{ii} and radii r_{i} < w_{ii}. Hence they must have real positive parts.

We shall now show that (iii)⇒(iv) holds. First note that since the determinant of a square matrix is the product of its eigenvalues and the product of complex eigenvalues of a real matrix is positive, it follows that all principal minors of the matrix W are positive. Next, the matrix W has the form

W = AD = [ \begin{matrix}a_{11} & −a_{12} & … & −a_{1n} \\ −a_{21} & a_{22} & … & −a_{2n} \\ ⋮ & ⋮ & ... & ⋮ \\ −a_{n1} & −a_{n2} & … & a_{nn} \end{matrix} ] [ \begin{matrix}d_{11} & … & & 0 \\ & d_{22} \\ ⋮ & & ... & ⋮ \\ 0 & … & & d_{nn} \end{matrix} ] ,

where a_{ij} ≥ 0 for i ≠ j, a_{ii} > 0, and d_{ii} > 0. Thus if the principal minors of W are positive, the principal minors of A must be positive, and vice versa. This is because

W ( \begin{matrix}i_{1} & i_{2} & … & i_{k} \\ i_{1} & i_{2} & … & i_{k} \end{matrix} ) = A [ \begin{matrix}i_{1} & i_{2} & … & i_{k} \\ i_{1} & i_{2} & … & i_{k} \end{matrix} ] [ \begin{matrix}d_{i_{1}i_{1}} & … & & & & 0 \\ & d_{i_{2}i_{2}} \\ & & . & . & . \\ 0 & … & & & & d_{i_{k}i_{k}} \end{matrix} ] ,

where

W ( \begin{matrix}i_{1} & i_{2} & … & i_{k} \\ i_{1} & i_{2} & … & i_{k} \end{matrix} )

is a matrix formed from i_{1}, i_{2}, …, i_{k} columns and i_{1}, i_{2}, …, i_{k} rows of the matrix W.

The implication (iv)⇒(v) follows from the observation that the leading principal minors are also principal minors.

The implications (v)⇒ (vi) ⇒(vii) follow from Theorem 7.

Finally, to complete the proof, we show that (vii)⇒(i). For this note that if A^{−1} exists and A^{−1} ≥ 0, then

x = A^{−1}e ≥ 0,

where e = [1, 1, …, 1]^{T} ∈ ℝ^{n}. Hence

Ax = e > 0,

which implies that (i) is satisfied.

□

References

[1] M. Fiedler and V. Pták. On matrices with non-positive off-diagonal elements and positive principal minors. Czechoslovak Journal of Mathematics, 12(87):382–400, 1962.

[2] H. K. Khalil. Nonlinear Systems. Prentice Hall, Upper Saddle River, NJ 07458, second edition, 1996.

[3] J. P. LaSalle. The Stability and Control of Discrete Processes. Springer-Verlag, New York, 1986.

[4] T. S. Parker and L. O. Chua. Chaos: A tutorial for engineers. Proceedings of the IEEE, 75(8):982–1008, August 1987.