Feedback Stabilization Methods for the Numerical Solution of Systems of Ordinary Differential Equations

In this work we study the problem of step size selection for numerical schemes, which guarantees that the numerical solution presents the same qualitative behavior as the original system of ordinary differential equations, by means of tools from nonlinear control theory. Lyapunov-based and Small-Gain feedback stabilization methods are exploited and numerous illustrating applications are presented for systems with a globally asymptotically stable equilibrium point. The obtained results can be used for the control of the global discretization error as well.


Introduction
It is well-known that step size control can enhance the performance of a numerical scheme when applied to a system of Ordinary Differential Equations (ODEs). In fact the use of the word "control" suggests that methods and techniques used in Mathematical Control Theory can be (in principle) used in order to achieve certain objectives for the numerical solution of systems of ODEs. For example, in [6] the authors use a "Proportional-Integral" technique which is similar to the "Proportional-Integral" controller used in Linear Systems Theory in order to keep the local discretization error within certain bounds (see also [7,8,13]). Theoretical results on the behavior of adaptive timestepping methods have been presented in [20,22] and the control theoretic notion of input-to-state stability (ISS) has been successfully used in [10,11] in order to explain the behavior of attractors under discretization.
In this work, we develop tools for nonlinear systems which are similar to methods used in Nonlinear Control Theory. We consider the problem of selecting the step size for numerical schemes so that the numerical solution presents the same qualitative behavior as the original system of ODEs. It is well-known that any consistent and stable numerical scheme for ODEs inherits the asymptotic stability of the original equation in a practical sense, even for more general attractors than equilibria see for instance [10,11] and [28] (Chapter 7). Practical Asymptotic stability means that the system exhibits an asymptotically stable set close to the original attractor, i.e., in our case a small neighbourhood around the equilibrium point, which shrinks down to the attractor as the time step h tends to 0. In contrast to these results, in this paper we investigate the case in which the numerical approximation is asymptotically stable in the usual sense, i.e., not only practically.
Here, we concentrate on nonlinear systems for which an equilibrium point is the global attractor. In Section 2 of the present work, it is shown how the problem of appropriate step size selection can be converted to a rigorous abstract feedback stabilization problem for a particular hybrid system (see also [17]-the reader should notice that the standard stability analysis of numerical schemes uses the study of a discrete-time system e.g., [13,14,16,21,28]; not a hybrid system). Therefore, we are in a position to use all methods of feedback design for nonlinear systems. Specifically, we consider methods based on Small-Gain Theorems, methods based on Lyapunov functions. Both methods have been used widely in Nonlinear Systems Theory for the solution of feedback stabilization problems (see [1,3,15,18,19,25,27] and references therein). In the present work, the above methods are used for the step size selection for numerical schemes for ODEs (see Section 3 and Section 4). General results are developed for arbitrary consistent Runge-Kutta schemes (see Theorem 4.5, Theorem 4.9 and Theorem 4.12 below) and specific results are given for specific Runge-Kutta schemes (see Corollary 4.7 and Theorem 4.15 below). The analysis is based on the Control Lyapunov Function methodology. The obtained results constitute nonlinear extensions of well-known properties of numerical schemes (e.g., A-stability, cf. Corollary 4.16). On the other hand small-gain methods allow the proposal of novel numerical schemes (see Theorem 3.1 below) for certain classes of nonlinear systems.
A number of applications of the obtained results is developed in Sections 5, 6 and 7. More specifically, we consider the possibility of using explicit schemes for stiff linear systems of ODEs: our results constitute rigorous nonlinear extensions of the ideas presented in [2] (see Section 6). Furthermore, we consider the problem of controlling the global discretization error: our results show that rigorous global discretization error control is possible only after stabilization of the numerical scheme (by appropriate step size selection-see Section 7). Another application of stabilization methods based on Small-Gain analysis to systems described by Partial Differential Equations is presented in Section 5.
Thus, the contribution of the paper is twofold. On the one hand, our control theoretic approach yields new insight into the stability properties of numerical schemes and as such it adds another means to the toolbox for stability investigations of numerical schemes. On the other hand, our method leads to the design of new discretization schemes and step size control algorithms, which beyond the mere control of the local discretization error also take care of the global qualitative behaviour.
Notations Throughout this paper we adopt the following notations: * Let n A ℜ ⊆ be a set. By ) ; ( , we denote the class of continuous functions on A , which take values in Ω . is an integer, we denote the class of differentiable functions on A with continuous derivatives up to order k , which take values in Ω . By ) ; ( i.e., We say that a function . We say that a continuous function is radially unbounded if the following property holds: "for every 0 is compact". * For a sufficiently smooth function ℜ → ℜ n V : we denote by ) ( ) ( : the Lie derivative of V along f and we define recursively ( )

Description of the Problem
Consider the autonomous system is a locally Lipschitz vector field with 0 ) 0 ( = f . For every n z ℜ ∈ 0 and 0 ≥ t , the solution of (2.1) with initial condition 0 ) 0 ( z z = will be denoted by ) , ( 0 z t z . The numerical approximation of system (2.1) will be the hybrid system:  [17]): Step i : 1) Given i τ and We will further assume that there exists a continuous, non-decreasing function It should be noticed that the hybrid system (2.2) under hypothesis (2.3) is an autonomous system, which satisfies the "Boundedness-Implies-Continuation" property and for each locally bounded input + + ℜ → ℜ : u and n x ℜ ∈ 0 there exists a unique absolutely continuous function which satisfies (2.2) (see [17]). Some remarks are needed in order to justify the name "numerical approximation of system (2.1)" for the hybrid system (2.2): is a consistency condition for the numerical scheme applied to obtained by using the Runge-Kutta numerical scheme (2.4), (2.5) and using the discretization step sizes )) ( exp( )) ( ( for some continuous, non-decreasing function is the function involved in (2.3). The reader should notice that appropriate step size restriction can always guarantee that (2.7) holds for as defined by (2.6). Notice that the implicit function theorem for (2.4) guarantees for each fixed . However, it must be emphasized again that the above step size restriction may be conservative (e.g., for explicit schemes).
Using Theorem II.3.1 in [12], (2.7), the fact that ) ; ( If we further assume that there exists a neighborhood satisfying the following properties: there exists a constant 0 > Λ and an integer 1 then it follows from (2.8a) that there exists a neighborhood N N ⊆ with Ñ 0 ∈ and a constant 0 is Uniformly Globally Asymptotically Stable (UGAS) for (2.1) (in the sense described in [24]; see also [17,19]). Our goal is to be able to produce numerical solutions using the numerical approximation (2.2) which have the correct qualitative behavior. More specifically, we would like to be in a position to know a continuous function ] , 0 ( : r n → ℜ ϕ so that the numerical solution produced by (2.2) has the correct qualitative behavior (e.g., ). However, we would like to be able to guarantee that the correct behavior for the numerical solution can be obtained by using arbitrary discretization step sizes smaller than )) ( ( i x τ ϕ (i.e., if we obtain the correct qualitative behavior using the discretization step sizes )) ( , we would like to obtain the correct qualitative behavior using the discretization step sizes )) ( exp( )) ( ( is an arbitrary locally bounded function). This is equivalent by requiring that n ℜ ∈ 0 is Uniformly Robustly Globally Asymptotically Stable (URGAS) for (2.2) (in the sense described in [17]).
The reader should notice that continuity for the function ] , 0 ( : r n → ℜ ϕ is essential: without assuming continuity it may happen that 0 ) ( inf lim 0 = → x x ϕ and this would require discretization step sizes of vanishing magnitude as +∞ → t . Moreover, since we want to be able to determine a continuous function ] , 0 ( : r n → ℜ ϕ , which "stabilizes" the hybrid system (2.2), we are essentially studying a feedback stabilization problem for the hybrid system (2.2). Hence, we are in a position to pose the problem rigorously. We consider the following feedback stabilization problems: It is well known that any consistent and stable numerical scheme for ODEs inherits the asymptotic stability of the original equation in a practical sense, even for more general attractors than equilibria see for instance [10,11] or [28] (Chapter 7). Practical Asymptotic stability means that the system exhibits an asymptotically stable set close to the original attractor, i.e., in our case a small neighbourhood around the equilibrium point, which shrinks down to the attractor as the time step h tends to 0.
Here, the property we are looking for, i.e., "real" asymptotic stability, is a stronger property which cannot in general be deduced from practical stability. In [28] (Chapter 5), several results for our problem for specific classes of ODEs are derived using classical numerical stability concepts like A-stability, B-stability and the like. In contrast to this reference, in the sequel we use nonlinear control theoretic analysis and feedback design techniques; more precisely Small-Gain and Lyapunov function techniques in Sections 3 and 4, respectively for solving Problems (P1) and (P2). This allows us to obtain asymptotic stability results under different structural assumptions and for more general classes of systems as in [28] (Chapter 5).

Small-Gain Methodology
One of the tools used in mathematical control theory for nonlinear feedback design is the methodology based on small-gain results. The methodology was first used in [15] where a nonlinear small-gain result based on the notion of Input-to-State Stability (ISS-see [26]) was presented. Since then it has been applied successfully to many feedback stabilization problems. Recently, the small-gain theorem was extended to general control systems including hybrid systems (see [18]). Consequently, the small-gain methodology for feedback design appears to be applicable for the solution of problem (P2) for certain classes of nonlinear systems (2.1). Consider the following system: We also assume that m ℜ ∈ 0 is UGAS for (3.1a). Under the previous assumptions, using the fact that system (3.1) has a structure of systems in cascade, we may prove by induction that for every n j ,..., The proof for 1 = j is based on the fact that for every ℜ ∈ 10 x and for every measurable = satisfies the following estimate: Consequently, the solution of ) ( ) ( satisfies the uniform ISS property from the input is GAS for (3.1a), a well-known corollary of the small-gain theorem (systems in cascade) guarantees UGAS for the composite system. The proof is similar for all n j ,..., 1 = .
Suppose that a stable numerical scheme is available for (3.1a), i.e., there exist ]) , is URGAS for the hybrid system: We propose the following first order numerical scheme for (3.1b): The above scheme is a partitioned scheme which treats differently the state i x and the states 1 1 ,..., , for each differential equation. The resulting hybrid system is system (3.5) with: We have: is uniformly RGAS for system (3.5), (3.7).
The proof of the above theorem is based on the following technical lemma.

Lemma 3.2: Let
Proof: Notice that for every 0 ≥ i it holds that: and using the definition 0 , we obtain the following bound from (3.10): Using the above inequalities in conjunction with (3.11) we obtain for all 0 ≥ i : Notice that for every 0 it holds that: . Combining (3.12) and (3.13) we obtain (3.13). The proof is complete. < The proof of the Theorem 3.1 follows induction: we show that for every n k ,..., . Remark 3.2(b) in [18] (systems in cascade) guarantees URGAS for the composite system.

Lyapunov function based Step Selection
While the small-gain methodology can be applied to certain systems of differential equations and can yield numerical methods which guarantee numerical stability, it cannot be applied in a systematic way. On the other hand the Lyapunov-based feedback design methods can be applied to general nonlinear systems of differential equations and yield explicit formulas for the feedback law (see [25]). In this section we apply the Lyapunov-based feedback design methodology for the solution of Problems (P1) and (P2). It is well known that Lyapunov functions exist for every asymptotically stable ODE system and in many applications one can even give explicit formulas for these functions (some examples can be found in Section 6). However, even if a Lyapunov function is not exactly known, under suitable assumptions on the ODE system, certain structural properties of the Lyapunov function can be obtained (cf., e.g., Proposition 4.4, below) and used in our context. Hence, the main task of this section is to derive conditions under which the Lyapunov function for the ODE system can be used in order to conclude stability for the hybrid system (2.2), i.e., for the numerical approximation of system (2.1).
The results will be developed in the following way: first (subsection 4.I) we provide some background material needed for the derivation of the main results. In Subsection 4.II we consider general consistent Runge-Kutta schemes and provide sufficient conditions for the solvability of Problem (P1) and (P2). The results are specialized for the explicit Euler method. Finally, in Subsection 4.III, we present special results for the implicit Euler numerical scheme.

4.I. Background Material
The crucial technical result that allows the use of Lyapunov functions for hybrid systems of the form (2.2) is the following lemma. Here it should be recalled that and a constant 0 > σ such that for every locally bounded (an extension of the corresponding notion for systems described by ODEs, see [23]).

Lemma 4.1: Consider system (2.2) and suppose that there exist a continuous, positive definite and radially unbounded function
and a continuous, positive definite function such that for every n x ℜ ∈ the following inequality holds for all )] is robustly K-exponentially stable for (2.2).
Proof: Notice first that by virtue of (2.3) the following claim holds: Let 0 ≥ R (arbitrary) and consider the solution ) ; , ( 0 u x t x of (2.2) corresponding to arbitrary locally bounded is continuous, positive definite and radially unbounded, it follows from Lemma 3.5 in [19] that there exist functions ∞ ∈ K a a 2 1 , such that the following inequality holds: Using induction and (4. Inequality (4.3) in conjunction with (4.2) and the above claim shows that It is clear that, by virtue of (4.2), (4.5) and the above claim it follows that ε . The previous estimate implies uniform robust global attractivity.
We next establish (4.5) by contradiction. Let 0 > ε (arbitrary). Suppose on the contrary that there exists . The previous inequality in conjunction with inequalities (4.1), (4.4) and definition (4.6) implies ) , ( )) ; , Furthermore, notice that by virtue of (2.3) there exists a continuous non-decreasing mapping . By virtue of the above claim and previous inequalities we get . The previous inequality implies that n ℜ ∈ 0 is robustly K-exponentially stable for (2.2). The proof is complete. < The essential problem with the use of Lemma 4.1 is the knowledge of the Lyapunov function V . In the sequel, we will use a Lyapunov function for the continuous-time system in order to construct a Lyapunov function for its hybrid numerical approximation. To this end we use the following definition.

Definition 4.2: A positive definite, radially unbounded function
is called a Lyapunov function for system (2

.1) if the following inequality holds for all
In the following subsections, we show that under certain assumptions a Lyapunov function V for (2.1) can be used as a Control Lyapunov Function (see [1,25,27]) in order to design the step function ] , 0 ( : involved in problems (P1) and (P2). Therefore, the Lyapunov function for the original system (2.1) will be used as a Lyapunov function for its numerical approximation (2.2). The following technical results will be used in the subsequent subsections and their proofs are provided at the Appendix.

Lemma 4.3:
be a Lyapunov function for system (2.1). Then the following statements hold: (i) There exists a locally Lipschitz, positive definite function such that the following inequality holds for all n x ℜ ∈ : . Then for every , there exists a continuous, positive definite function such that the following inequality holds for all be the function from statement (i) above and let (4.14)

4.II. General Runge-Kutta Schemes
The following theorem provides sufficient conditions for the solvability of problem (P2) based on the Lyapunov function for the dynamical system (2.1).

Theorem 4.5:
Suppose that there exists an integer 1 ≥ p and a Lyapunov function ) ; ( for system (2.1). Consider system (2.2) that corresponds to a Runge-Kutta scheme for (2.1) and suppose that: is robustly K-exponentially stable for (2.2).
Proof: Since for each fixed n x ℜ ∈ the mapping )) , times continuously differentiable, we have from Taylor's theorem for all p j ,..., Since the Runge-Kutta numerical scheme is of order Consequently, we obtain for all p j ,..., or equivalently for all )] robustly K-exponentially stable for (2.2). The proof is complete. < Remark 4.6: (a) Theorem 4.5 shows that given a Runge-Kutta numerical scheme with order 1 ≥ p satisfying (2.7) and a system of , then for every ) 1 , 0 ( ∈ λ and every This fact follows from (2.7) and the observation that ϕ for x close to zero. Consequently, the numerical solution of (2.1) with sufficiently small step size will give the correct dynamic behavior. x h Consequently, we obtain the following corollary.

Corollary 4.7-Explicit Euler method: Suppose that there exists a Lyapunov function
is simply locally Lipschitz and that there exist constants , then for every ) 1 , 0 ( ∈ λ and every This fact follows from (2.7) the observation that ).
The following theorem provides sufficient conditions for the solvability of problem (P2) based on the Lyapunov function for the dynamical system (2.1). It should be noticed that the hypotheses of the following theorem are of different nature from the hypotheses (i)-(iii) of Theorem 4.5. In particular, the conditions in the following theorem do not require the Lyapunov function to be smoother than 1 C . ii) There exists iii) There exists a constant for all n x ℜ ∈ and is a continuous positive is the continuous positive definite function with . The above inequality in conjunction with (4.21) implies: then for every ) 1 , 0 ( ∈ λ and every compact This fact follows from (2.9) and the observation that for x close to zero. The reader should notice the differences with Remark 4.6(a) (less regularity requirements). Example 4.14 below will show that in practice it is not necessary to compute )) , 0 ( ; ( 0 +∞ ℜ ∈ n C ϕ in order to be able to produce a qualitatively correct numerical solution; you only need to know a Lyapunov function for (2.1) with the properties described above.
The following example illustrates Remark 4.6 and Remark 4.10.

Example 4.11: Consider three planar systems with
) For each of the systems we can use the Lyapunov function Then we obtain: . More specifically, we have 1 2 1 = = q q and 3 3 = q . Remark 4.6(a) shows that for 1 = k and 3 = k we are in a position to apply any consistent Runge-Kutta numerical scheme with sufficiently small step size and produce a qualitatively correct numerical solution. The same conclusion is derived by Remark 4.10 (notice that (4.11) holds for each of the systems with ) ( ) ( : ).
On the other hand, the requirements presented in Remark 4.6(a) or Remark 4.10 are not fulfilled for 2 = k . Similarly, the requirements presented in Remark 4.8 are not fulfilled for 2 = k . Consequently, we cannot conclude that the application of any consistent Runge-Kutta numerical scheme with sufficiently small step size will produce a qualitatively correct numerical solution. Numerical solutions with the explicit Euler and the Heun scheme confirm these results (see Figure 1 and Figure 2 for both numerical schemes show the correct behavior (see Figure 3 and Figure 4).
We would like to emphasize that Theorem 4.5 and Theorem 4.9 do not state that there does not exist a Runge-Kutta scheme which produces an asymptotically stable approximation for system ) , since it gives a merely sufficient but not necessary condition. In fact, for instance the implicit Euler scheme produces an asymptotically stable approximation, which we will rigorously show in Example 4.19 below. < Based on the general Theorem 4.9, the following theorem shows that for the special case of a locally exponentially stable ODE system, problem (P1) is always solvable.
Moreover, the proof of Lemma 4.3 shows that (see (A3)) for all . We notice that the inequality is any function satisfying with is a continuous positive definite function with (4.13). Clearly, (4.13) implies that is the symmetric, positive definite matrix involved in (4.13). Notice that without loss of generality we may assume that there exists constant 0 > K such that (2.9) holds with is a consequence of local exponential stability). Consequently, by virtue of (2.9), (4.26), we can guarantee that (4.25) holds for every ) . Therefore, from all the above we conclude that we may define , so that all requirements of Theorem 4.9 are fulfilled. The proof is complete. < Remark 4.13: Theorem 4.12 is an existence result which does not give an explicit estimate for ) (x ϕ , i.e., it answers (P1) but does not solve (P2). However, similar to Remark 4.8 and 4.10 we can conclude that the numerical approximation is asymptotically stable on each compact set S for sufficiently small step size h. Note that local exponential stability is not a necessary condition for asymptotic stability of explicit Runge-Kutta schemes, as Example 4.11 shows ( 2 0 ℜ ∈ is GAS but not locally exponentially stable for system ) ).
Due to the complexity of calculations needed in order to satisfy the requirements of Theorem 4.5 and Theorem 4.9, we suggest the following algorithm for Runge-Kutta schemes and a given ) Otherwise halve the step size (i.e., 2 / h h = ) and check again" is a Lyapunov function for (2.1) for which there exist a constant 0 > K and a neighborhood . Therefore, in practice we do not have to compute the step function ) (x ϕ that guarantees robust global asymptotic stability of the numerical approximation. The following example illustrates this point.

Example 4.14:
In this example we consider four different explicit numerical schemes: the explicit Euler scheme, Heun's scheme, the improved polygonal scheme (a second order method) and Kutta's 3 rd order scheme. The numerical schemes are applied to the planar system: . For all numerical schemes (except the explicit Euler method) the calculation of the maximum allowed time-step by using Theorem 4.5 or Theorem 4.9 is very complicated. However, holds. The following figure shows the graph of the maximum allowable time-step for the four numerical methods with It should be noticed that all higher order schemes allow greater time-steps for 1 x close to zero than the time-steps allowed by the explicit Euler method (notice that for holds for the explicit Euler method is 2 ). However, for large values of 1 x the maximum allowable time-step for higher order schemes are considerably smaller than the time-step allowed by the explicit Euler method. It is clear that a higher order method does not necessarily allow higher values for the maximum allowable time step for which the inequality )

4.III. Implicit Runge-Kutta Schemes
In this subsection we show how Lyapunov function based arguments can be used for implicit schemes. In order to keep the presentation technically simple, we restrict ourselves to the implicit Euler scheme for which we can prove the following result.  is URGAS for the corresponding system (2.2) Proof: Define the functions By virtue of (4.7) both functions are continuous and positive definite.
is convex the following inequality holds for all n x x ℜ ∈ 2 1 , : We apply (4.29) with . We get: By virtue of (4.7), (4.30) implies that ) We distinguish the following cases: In this case (4.30) in conjunction with definition (4.28) of 1 W we obtain: Consequently, in any case we obtain for all )] ( , 0 [ x h ϕ ∈ and n x ℜ ∈ : = is a positive definite function. Lemma 4.1 yields the assertion. The proof is complete. < The following corollary shows that Theorem 4.15 can be seen as a nonlinear generalization of the well-known Astability property of the implicit Euler method.

Corollary 4.16: Consider the system of ODEs
is a Hurwitz matrix. Then the implicit Euler method is URGAS for arbitrary step size 0 > h .
Proof: As pointed out in Section 2, the implicit Euler method is well defined for each step size 0 > h . Furthermore, the system Ax is a symmetric, positive definite matrix (see Theorem 5.7.18 in [27]). This Lyapunov function is obviously convex and Theorem 4.15 yields the assertion. <

Remark 4.17:
The main result in [9] guarantees that if 5 , 4 ≠ n then there exists a homeomorphism n n , being a diffeomorphism on }) 0 { \ ( n ℜ and 1 C on n ℜ such that the transformed system (2.1) admits the convex Lyapunov function Consequently, the implicit Euler can be applied to the transformed system. However, for numerical purposes, the method is not practical, since the homeomorphism is usually not available. On the other hand, for certain classes of systems Theorem 4.15 is directly applicable. One such class are the so called gradient systems, as shown in the following example.

Example 4.18:
Consider the following class of systems: The class of systems of the form (4.34) is a generalization of the class of the so-called gradient systems (see [28]).
It follows from Theorem 4.15 that the implicit Euler scheme can be applied successfully for the numerical approximation of the solutions of (4.34) for every 0

Application of the Small-Gain Step Selection Methodology
Consider the infinite-dimensional dynamical system: being locally Lipschitz, 0 > c and initial condition ) , under the following hypothesis Using the method of characteristics and hypothesis (H), it can be shown that the infinite-dimensional dynamical system (5.1) admits a unique classical solution ) ]; Moreover, the zero solution is globally asymptotically stable, since for every ) ]; (uniform global attractivity).
Using a uniform space grid of 1 + n points with space discretization step and approximating the spatial derivative by the backward difference scheme , we obtain the following set of ordinary differential equations: It is clear that system (5.2) has the structure of system (3.1) with ) ( ) ( is the constant involved in Hypothesis (H), then inequalities (3.2) hold as well with K z . Theorem 3.1 allows us to conclude that for every 0 > h the numerical scheme: will give the correct qualitative behavior. The reader should notice that for the case 0 ) and the numerical scheme (5.4) is related to the so-called implicit upwind numerical scheme for the advection equation, which is unconditionally stable.

Applications of the Lyapunov-based Step Selection Methodology
In this section we present some applications for the Lyapunov-based step selection methodology that was provided in the Section 4. It should be emphasized that the Lyapunov-based step selection methodology can be (in principle) applied to all dynamical systems for which a Lyapunov function is known with a globally asymptotically stable and locally exponentially stable equilibrium (Theorem 4.7). However, as the following applications show there are certain classes of systems that we can guarantee more properties or require special attention.

Application 1: Solution of Nonlinear Programming Problems
There are many nonlinear programming problems which can be solved by constructing a dynamical system with a globally asymptotically stable equilibrium point which coincides with the minimizer of the nonlinear programming problem (see [5,29,30,31]). A special feature for such methods is that a Lyapunov function is available; however the position of the equilibrium point is not known (this is what we seek). Consider the following nonlinear programming problem: Problem (P) may be solved by means of differential equations if we further assume that the function ( ) has the unique equilibrium point , which is UGAS for (6.2). This fact can be proved by using the Lyapunov function (notice that it is radially unbounded). Notice that ( ) . Thus the dynamical system (6.2) can be solved by means of Runge-Kutta methods with a Lyapunov-based step selection methodology: each Runge-Kutta method applied to the dynamical system (6.2) will yield a method for the solution of the nonlinear programming problem (P).
Here we will discuss the explicit Euler method.  (6.4) and it is clear that (6.4) will achieve convergence to the (unknown) equilibrium point of (6.2), provided that (the

Application 2: Control Systems under Feedback Control
One class of dynamical systems for which a Lyapunov function is known is the class of control systems for which a continuous feedback stabilizer is designed by using a Lyapunov-based methodology (see [1,3,19,25]). This is evident for the class of so-called "triangular control systems" (see [3]). Consider the following triangular control system: ).
Using backstepping (see [3]), we are in a position to construct a smooth function ℜ → ℜ n k : with 0 ) 0 ( = k and a positive definite and radially unbounded smooth function is locally exponentially stable for the closed-loop system (6.6) with ) (x k u = and for every 0 ≥ Δ there exist constants 0 , 2 such that the following inequality holds: Consequently, Corollary 4.7 guarantees that the explicit Euler method can be used for the numerical approximation of the closed-loop system (6.6) with ) (x k u = . Furthermore, Corollary 4.7 can be used in order to obtain explicit estimate of the allowable discretization time step for the explicit Euler method. Indeed, notice that all requirements of Corollary 4.7 hold (notice that (6.7), in conjunction with (6.8) show that for every bounded neighborhood of the origin). Formula (4.20) (combined with (6.7)) provides an explicit upper bound for the function ]) , 0 ( Other Runge-Kutta numerical schemes can be used as well. Notice that the backstepping procedure achieves the construction of the Lyapunov function by constructing a diffeomorphism n n is a symmetric, positive definite matrix. Then Theorem 4.15 guarantees that the implicit Euler can be used as well for the transformed closed-loop system (6.6) with ) (x k u = : It follows that for every 0 > r , ) 1 , 0 ( ∈ λ , the implicit Euler scheme can be applied for (6.11) with . This fact was observed in [17].
is the stability function of the scheme and h is the (constant) discretization step size. The possibility of using larger discretization step size for explicit Runge-Kutta methods than the one allowed by the classical analysis must be considered. This possibility was considered in [2] where it was shown that a sequence of "small" time steps can allow a "big" time step for explicit ODE solvers.
Here for simplicity, we consider the Explicit Euler scheme. The fact that guarantees that the numerical solution produced by the explicit Euler scheme has the correct qualitative behavior.
Notice that the quantity depends heavily on the direction of the vector n x ℜ ∈ and can allow greater discretization step sizes than the one produced by classical stability analysis.
Example 5.1: Consider the singularly perturbed and stiff linear system Classical analysis (Dahlquist test) for the explicit Euler scheme with constant step size requires that the discretization step must satisfy 500 (a value for the discretization step size comparable to the value provided by the limit of the classical analysis). Figure 6 shows the numerical solution obtained by using the explicit Euler scheme for (6.14) with step size Figure 7 shows the sequence of discretization time steps for the explicit Euler scheme with step size . It can be seen that the numerical solution presents the correct qualitative behavior. Moreover, Figure 7 shows a repeated pattern for the step selection: after a number of "small" time-steps the scheme allows a "big" time step. This is in agreement with the results in [2]. However, the "big" time step creates errors in the evaluation of the state variable ) ( 1 t x and this explains the behavior presented in Figure 6. After 500 explicit Euler steps the value of time is 71372 . 12 = t : if we had applied 500 Euler steps with constant step size we would have reached at most 1 = t . The accuracy of the numerical solution can be improved by imposing a higher value for λ . Figure 8 shows the numerical solution obtained by using the explicit Euler scheme for (5.14) with step size . Figure 9 shows the sequence of discretization time steps for the explicit Euler scheme with step size ⎪ ⎭  Figure 9: The sequence of time steps for the Explicit Euler scheme for (5.14) with It is clear that a trade-off between the allowable time steps and the accuracy of the numerical solution exists, as expected. <

Stabilization of Numerical Schemes and Control of the Global Discretization Error
It is well-known that the numerical solution produced by a Runge-Kutta numerical scheme of order The global discretization error ) ; , for the numerical solution with initial condition Clearly we have such that   Therefore, by virtue of (7.7), the demand ε τ

Conclusions
In this work, we considered the problem of step size selection for numerical schemes such that the numerical solution presents the same qualitative behavior as the original system of ODEs. Specifically, we developed tools for nonlinear systems with a globally asymptotically stable equilibrium point which are similar to methods used in Nonlinear Control Theory. It is shown how the problem of appropriate step size selection can be converted to a rigorous abstract feedback stabilization problem for a particular hybrid system. Feedback stabilization methods based on Lyapunov functions and Small-Gain results were employed. The obtained results have been applied to several examples and have been shown to be efficient for controlling the global discretization error.
The methodology presented in the present work (transformation of the step selection problem to a feedback stabilization problem and use of modern nonlinear control theory for the solution of the problem) can be used for more complicated numerical problems such as the step size selection problem for: the numerical approximation of the solution of infinite-dimensional systems (systems described by partial differential equations or systems described by retarded functional differential equations), systems with more complicated attractors, time-varying systems, systems with inputs.
Future work will address these problems.  [4] is the unique function satisfying (A15) with 0 ) 0 ( = V . An inspection of the proof yields that if equation (A15) holds on a forward invariant set for (2.1) then uniqueness holds on this set (because uniqueness is established by looking at trajectories in forward time). Notice that by virtue of (A15) and definition (A13), it follows that (4.14) holds.
Picking a forward invariant neighborhood ) 0 ( δ B N ⊂ of zero (such a forward invariant neighborhood of zero exists because n ℜ ∈ 0 is asymptotically stable), we observe by (A13) that the function Px It follows that (4.13) holds. The proof is complete. <