Numerical schemes of higher order for a class of nonlinear control systems (cid:3)

We extend a systematic method for the derivation of high order schemes for a(cid:14)nely controlled nonlinear systems to a larger class of systems in which the control variables are allowed to appear nonlinearly in multiplicative terms. Using an adaptation of the stochastic Taylor expansion to control systems we construct Taylor schemes of arbitrary high order and indicate how derivative free Runge-Kutta type schemes can be obtained.


Introduction
Traditional numerical schemes for ordinary differential equations, such as Runge-Kutta schemes, usually fail to attain their asserted order when applied to ordinary differential control equations due to the measurability of the control functions. A similar situation occurs with stochastic differential equations due to the nondifferentiability of the driving noise processes. To construct higher order numerical schemes for stochastic differential equations, one needs to start with an appropriate stochastic Taylor expansion to ensure consistency with the less robust stochastic calculus as well as a higher order of convergence. This is the opposite procedure to that used for numerical schemes for ordinary differential equations, where heuristic arguments are typically used to derive a scheme and the Taylor expansion is then used to establish its local discretization order.
In [9] it was shown that this approach for stochastic differential equations carries over to control systems with affine control (for these systems the stochastic Taylor expansion is essentially the same as the Fliess expansion [11]). In the present paper we will extend the results from [9] to a larger class of control systems allowing also nonlinearities in the control input. More precisely, we consider d-dimensional controlled nonlinear system with n-dimensional control functions of the form where t ∈ [t 0 , T ], x = (x 1 , . . . , x d ) ∈ IR d , the vector fields f j : IR ×IR d → IR d are sufficiently smooth in order to apply our expansion, the functions g j : IR × IR n → IR are continuous and the control functions u(t) are measurable and take values in a compact set U ⊂ IR n .
Numerical schemes for such systems play an important role in the numerical analysis of nonlinear control systems since in many algorithms the approximation of trajectories appears as a subproblem, see, e.g., the monographs [2] and [8].
The organization of this paper is as follows. We start with the introduction of the necessary notation in Section 2 and the precise statement of the Taylor expansion in Section 3. In Section 4 we explain how numerical Taylor and derivative free (i.e., Runge-Kutta type) schemes can be obtained, and finally in Section 5 we show a numerical example.

Setup and Notation
In the following sections we shall refer to the nonautonomous d-dimensional controlled differential equation (1), which we rewrite in the equivalent compact integral form where we set g 0 (t, u) ≡ 1 so that the first integral term can be included in the summation. We call a row vector α = (j 1 , j 2 , . . ., j l ), where j i ∈ {0, 1, . . ., m} for i = 1, . . ., l, a multi-index of length l := l(α) ≥ 1 and for completeness we write for the multi-index of length zero, that is, with l( ) = 0. We denote the set of all such multi-indices by M m .
For a multi-index α = (j 1 , j 2 , . . ., j l ) ∈ M m , some integrable control function u : IR → U m and an integrable function f : [t 0 , T ] → IR we define the multiple integral I α [f (·)] t 0 ,t recursively by We note that I α [f (·)] t 0 ,· : [t 0 , T ] → IR is continuous, hence integrable, so the integrals are well defined. For example, we obtain For simpler notation, we shall often abbreviate I α [f (·)] t 0 ,t to I α,t or just I α when f (t) ≡ 1 and shall explicitly write I α,u [f (·)] t 0 ,t , I α,u,t or I α,u when we want to emphasize a specific control function u.
For each α = (j 1 , . . ., j l ) ∈ M m and a function F : where the partial differential operators are defined by This definition requires the functions F , f 0 , f 1 , . . ., f m to be sufficiently smooth. For example, in the autonomous scalar case with d = m = 1 for the identity function where the dash denotes differentiation with respect to x. When the function F is not explicitly stated in the text we shall always take it to be the identity function Since different integrals can be expanded in forming a Taylor expansion, the terms with constant integrands cannot be written down completely arbitrarily. Rather, the sets of corresponding multi-indices must fo rm hierarchical and remainder sets. These sets can be defined in a very general way, see [13]. Here we only need the hierarchical and remainder sets defined by

Taylor expansions and approximations
We now formulate the Taylor expansion for the d-dimensional controlled system (2) using the terminology from the preceding section.

Theorem 1
Let F : IR + ×IR d → IR. Then for each N ≥ 0 the following Taylor expansion holds, provided all of the partial derivatives of F , f 0 , f 1 , . . ., f m and all of the multiple control integrals appearing here exist.
For the proof we refer to [9, Theorem 1], whose proof immediately carries over to our class of systems.
Based on Theorem 1 we can now construct Taylor approximations of arbitrary higher order. In the general multi-dimensional case d, m= 1, 2, . . . the Taylor approximation for N = 1, 2, 3, . . . is defined by with the coefficient functions F α corresponding to the function F (t, x) . When the function F (t, x) is N + 1 times continuously differentiable and the drift and control coefficients f 0 , f 1 , . . ., f m of the controlled differential equation (2) are N times continuously differentiable, then each of the integrals I α,t 0 ,t 0 +∆ (F α (·, x(·))) for α in the remainder set B(Γ N ) is of order ∆ N +1 . Since there are only finitely many, specifically (m + 1)!, remainder integrals, the truncation error here is where the constant K depends on N as well as on a compact set containing the initial value (t 0 , x(t 0 )) and the solution of the controlled differential equation.
For the function F (t, x) ≡ x k , the kth component of the vector x, and N = 1, 2 and 3, respectively, the solution x(t 0 + ∆) of the controlled differential equation (2) satisfies the componentwise approximations x and for k = 1, . . ., d, where we have written I (j) for I (j),t 0 ,t 0 +∆ and so on.

Numerical schemes
Using the Taylor approximation from the previous section we now construct numerical schemes by iterating Taylor approximations, or suitable derivative free approximations of those, over a partition of the time interval under interest. Schemes of arbitrary higher order N = 1, 2, . . . can be constructed by truncating the Taylor approximation corresponding to the the hierarchical set Γ N . Here we assume that the multiple control integrals I α are at our disposal, which is often feasible e.g. by using symbolic manipulators like maple. For a numerical approximation of such integrals see [9, Section 9]. Let {t 0 , t 1 , . . . , t n , . . . , } be a partition of the time interval [t 0 , T ] with stepsizes ∆ n = t n+1 −t n and maximal step size ∆ := max n ∆ n . In the general multi-dimensional case d, m = 1, 2, . . . for N = 1, 2, 3, . . . we define the Taylor scheme of order N for the controlled differential equation (2) is given componentwise by with the coefficient functions F k α corresponding to F (t, x) ≡ x k for k = 1, . . ., d and the multiple control integrals from (3). By standard arguments (see [12] or [10]) it follows from (6) that the global discretization error is of order N when the coefficients f j of the differential equation (2) are N times continuously differentiable.
Below, we write out the Taylor schemes for N = 1 and 2, where we distinguish the purely uncontrolled integrals, that is with multi-indices (0) and (0, 0) from the others, since no special effort is required for their evaluation.
The simplest nontrivial Taylor scheme is the Euler approximation with convergence order N = 1. It is given componentwise by for k = 1, . . ., d, where ∆ n = t n+1 − t n = I (0),tn,t n+1 . The kth component of the Taylor scheme of order N = 2 is given by L j 1 f j 2 ,k (t n , X n ) I (j 1 ,j 2 ),tn,t n+1 f or k = 1, . . ., d. For N = 3 we refer to [9]. A disadvantage of Taylor schemes is that the derivatives of various orders of the drift and control coefficients must be first derived and then evaluated at each step. Although nowadays symbolic manipulators [3] facilitate the computation of the derivatives in the schemes, it is useful to have approximations and schemes that avoid the use of derivatives of the drift and control coefficients in much the same way that Runge-Kutta schemes do in the more traditional setting since these often have other computational advantages.
Since the Euler or Taylor scheme of order 1 contains no derivatives of the f j , we illustrate this procedure for the second order Taylor scheme (11). In the autonomous case, from the ordinary Taylor expansion for f j we obtain Since O(∆ n ) I (i,j),tn,t n+1 = O(∆ 3 n ), the remainder is of the same order as the local discretization error if we replace the L i f j,k by this approximation.
In this way we obtain the second order derivative-free scheme in the autonomous case for k = 1, . . ., d. In the usual ODE case, that is with f j (x) ≡ 0 for j = 1, . . ., m, this is just the second order Runge-Kutta scheme known as the Heun scheme. This principle can be extended to obtain higher order derivative-free schemes. See [13] for analogous higher order derivative-free schemes for the stochastic case.
Note that all these schemes simplify considerable when the coefficients f j of the controlled differential equation (2) satisfy special properties. For example, if the control coefficients f 1 , . . ., f m are all constants or depend just on t, then all of the spatial derivatives of these control coefficients vanish and, hence, so do the corresponding higher order terms.
Another major simplification occurs under commutative control, that is when the f i satisfy L i f j,k (t, x) ≡ L j f i,k (t, x) for all i, j = 0, 1, . . ., m. Then, by the generalized integration-by-parts identities I (i,j),tn,t n+1 + I (j,i),tn,t n+1 = I (i),tn,t n+1 I (j),tn,t n+1 , i,j = 0, 1, . . ., m, we obtain L i f j,k (t n , X n ) I (i,j),tn,t n+1 + L j f i,k (t n , X n ) I (j,i),tn,t n+1 = L i f j,k (t n , X n ) I (i),tn,t n+1 I (j),tn,t n+1 , which involves more easily computed multiple control integrals of lower multiplicity. Note that this condition is similar to the one considered in [14], where the effect of time discretization of the control function is investigated and a second order scheme for the approximation of the reachable set is obtained.

A numerical example
We have tested the Euler (10) and Heun (12) Schemes from Section 4 with with control function u(t) = sin(100/t) and initial value x 0 = (0, 0) T . The resulting schemes have been simplified using the identity (13) such that the only remaining control integrals were I (1),0,t and I (0,1),0,t , which have been evaluated using maple. Note that the exact solution for this equation is easily verified to be x 1 (t) = I (1,0),0,t − I (1,1),0,t , x 2 (t) = I (1),0,t . The equation was solved on the interval [0, 1] with timestep ∆ = 1/N and N = 50, 100, . . ., 500. Figure 1 shows the resulting errors sup n=1,...,N x n − x(n∆) for the Heun and the Euler scheme. The left figure shows the error over N in a linear scale, the right figure shows the error over ∆ in a log-log scale.