Multivariate Canonical Polynomials in the Tau Method with Applications to Optimal Control Problems

The Tau method is a highly accurate technique that approximates differential equations efficiently. This paper discusses two approaches of the Tau Method: recursive and spectral. In the recursive Tau, the approximate solution of the differential equation is obtained in terms of a special polynomial basis called canonical polynomials. The present paper extends this concept to the multivariate canonical polynomial vectors and proposes a self starting algorithm to generate those vectors. In the spectral Tau, the approximate solution is obtained as a truncated series expansions in terms of a set of orthogonal polynomials where the coefficients of the expansions are obtained by forcing the defect of the differential equation to vanish at the some selected points. In this paper we use the spectral Tau to solve a class of optimal control problems associated with a nonlinear system of differential equations. Some numerical examples that confirm our method are given.


Introduction
The Tau method is a highly accurate technique that approximates differential equations without requiring the discretization of the given differential operator.Its basic idea is to perturb the right hand side of the differential equation in a way that an exact polynomial solution of the new equation can be found analytically.This method was devised in (Lanczos, 1956) to find polynomial approximations for simple linear ordinary differential equations (ODE) and it was extended later on to treat differential equations with different level of complexities, (see (Ortiz, 1969), (Ortiz & Samara, 1984), (El-Daou & Ortiz, 1992-1994), and (Liu & Pan, 1999)).
The Tau method has three equivalent approaches: Recursive, operational and spectral.The recursive approach, proposed in (Ortiz, 1969), permits to obtain an approximate polynomial solution expressed in terms of a special polynomials basis called canonical polynomials.This technique has been thoroughly investigated in a series of papers (see for example (Crisci & Russo, 1983), (Freilich & Ortiz, 1982) and (El-Daou & Ortiz, 1994-1998)).In the operational Tau, (see (Ortiz & Samara, 1981)), the ODE is transformed to a system of linear algebraic equations using some simple elementary matrices.The operational procedure was extended in (Liu & Pan, 1999) to solve mixed-order systems of linear ODEs with polynomial coefficients.In the same reference an automation of the operational approach has been reported.In (Ortiz & Samara, 1984) the operational Tau was shown effective in solving partial differential equations.In (Canuto, Hussaini, Quarteroni, & Zang, 2006) and (Gottlieb & Orszag, 1977 ), a spectral approach to the Tau method has been studied.This technique seeks an approximate solution in the form of a truncated series expansions of Chebyshev or Legendre polynomials.The coefficients of the series are computed by forcing the ODE to be exact at some selected points (called collocation points) and the supplementary conditions to be satisfied exactly.The spectral approach of the Tau Method guarantees spectral accuracy because the approximate solution is obtained in terms of orthogonal polynomials basis.
Although the three approaches of the Tau Method explained above appear to be different, it was shown in (El-Daou & Ortiz, 1992-1994) that they are equivalent in the sense that they yield the same approximate solution.However, the suitability of those approaches is judged by the problem under consideration.While the recursive and the operational Tau are more suitable from the computation point of view for ODE with polynomial coefficients, the spectral Tau enjoys a superiority when the ODE is nonlinear or its coefficients are not polynomials.
We point out that an important feature that distinguishes the Tau Method from the classical finite difference methods is that the Tau solution is obtained in a closed form on the whole domain of integration without discretization, while in the finite difference methods the domain is divided into small elements of stepsize h on which depends the accuracy of the method.
The present paper is to extend the recursive approach of the Tau Method to the numerical solution of systems of linear ODEs and to give a practical procedure that permits to construct approximate solutions.Further, we will show that the spectral Tau method is highly effective in tackling a class of optimal control problems (see (Flores Tlacuahuac, Terrazas Moreno, & Biegler, 2008) and (Jaddu & Majdalawi, 2014)).To this end, we generalize in Section 2 the concept of canonical polynomials to be adapted for system of ODEs.An algorithm to compute the Tau solution in terms of the canonical vectors will be given in Section 3. Section 4 is concerned with applying the spectral Tau to solve an optimal control problem whose the constraints are given as a system of nonlinear ODEs.Numerical examples illustrating the efficiency of our method are provided throughout the paper.

Canonical Polynomial Vectors
Let us consider a system of linear ODEs of dimension ν ≥ 1 written in the matrix form as where I ν is the ν identity matrix and A := (A i j (x)) ν i, j=1 is a ν × ν matrix with A i j (x) being functions of x.The superscript T means "Transpose".We shall assume for simplicity that all the A i j (x)'s are polynomials having the same degree d, and more compactly we can write where ; m = 0, 1, 2, . . ., d} are constant matrices.
Let us associate to (1) the initial conditions When ν = 1, system (1) reduces to a single equation (Dy)(x) = f (x).This case was fully discussed in Ortiz (1969) wherein a Tau solution ỹ(x) that approximates y(x) is obtained in the form with {q * k (x); k ≥ 0} being a sequence of functions, called canonical polynomials associated with D each one of which is an exact solution of the differential equation In (Ortiz, 1969, Theorem 3.3) a self starting recursive formula that generates the {q * k ; k ≥ 0} associated with a single ODE was developed.Next we extend the concept of canonical polynomial to systems of ODEs:

Definitions and Notation.
1. We call a vector Q (k)  i (x) an ith canonical vector of order k associated with the operator D if where e i is the ith column of The next theorem is a generalization of Theorem 3.3 given in (Ortiz, 1969): Theorem 1 Suppose that the matrix A d defined in ( 3) is non-singular.Then the canonical blocks {Q * k ; k ≥ 0} satisfy the recursion: In particular, if d = 0, then (5) becomes Proof.Let k ≥ 0 and i ∈ {1, 2, . . ., ν}.Let us apply the operator D, defined by (1), to the vector x k e i : The last identity is due to the linearity of D. Rearranging terms, and using the definition DQ (i) k = x k e i , we obtain: This implies (without loosing generality) that which gives Explicitly this means that for i = 1, 2, ..., ν . . .
In matricial form we have Since A d is non-singular, we obtain the desired recursion and this completes the proof of ( 5).
Although all the Q * k 's satisfy the self starting recursive formula (5), in practice only {Q * k ; k ≥ d} can be generated by this recursion while the remaining ones {Q * 0 , Q * 1 , ..., Q * d−1 } are not computable by the same formula.These are then called undefined.This point will be clarified further in the following example which illustrates Algorithm (5): Example 1.Consider the differential operator D defined by (1) with Since all the entries of A are of order 1, d = 1 and therefore (5) becomes where } is undefined because it can not be obtained from (6).However the execution of ( 6 2) 3) One observes that every This characteristic is not confined to the differential operator discussed in Example 1.In fact, for any operator of the form (1), every canonical block To justify this claim: Execute (5) when k = 0: for some constant square matrices {R j } depending on {A 0 , A 1 , . . ., A d−1 }.Proceeding this way we obtain by induction the following corollary: Corollary 2 Let {R ( j) k , k ≥ 0, s = 0, 1, ..., d − 1} be sequences of ν × ν matrices defined by the recursion: Then each Q * k , k ≥ 1, can be written as where {Q k , k ≥ 0, j = 1, 2, ..., ν} are defined by the recursion The proof this corollary follows from the fact that sequence {Q * k , k ≥ 0} is unique by construction and from comparing both sides of the identity:

Construction of the Tau Method Approximation
In the Tau method we associate to (1)-( 4) a perturbed problem of the form where H N is a free perturbation term adjusted in a way that Y N is a vector of polynomials that can be obtained analytically.Usually H N is chosen either of the form or with N being a prescribed positive integer, and are r vectors each one of which is formed of ν free parameters where the positive integer r will be fixed later.
V m (x) designates a polynomial of degree m that is usually chosen as the Chebyshev polynomial T m (x) or Legendre polynomial P m (x) shifted to the appropriate interval: The unknown vectors {τ j ; j = 1, 2, . . ., r − 1} are determined (i) by imposing the initial condition (9) on the desired approximate solution Y N , and (ii) by forcing Y N to be independent of the undefined terms {Q * j ; j = 0, 1, . . ., d − 1}.The latter can be realized by setting the coefficient of each Q * j ; j = 0, 1, . . ., d − 1 equal to zero.
Since ( 9) is the only initial condition to be satisfied, and since there are d undefined canonical blocks {Q * 0 , Q * 1 , ..., Q * d−1 }, we need d + 1 vectors τ j 's only i.e. we choose r = d + 1.
Proof.This follows once the right hand side of (8), f(x) + H N (x), is expressed in terms of {Q * i }.Let us consider first H N (x): , then H N (x) introduced in (10) can be written as Since DQ ( j) i = x i e j and D is linear, H N (x) becomes Hence, The same arguments apply to the function vector f(x) given by ( 2).Suppose that each entry in f(x) is a polynomial of degree α written as Therefore Adding ( 15) and ( 16) we get Thus the Tau problem ( 8) is written now in the form which implies that Y N is formally given by: Further, using ( 7) we can write Y N in terms of {Q i }: This expression holds for any choice of {τ k }.Since ( 19) contains undefined canonical blocks only we call it residual.
In order to eliminate this residual we set its coefficients equal to zero.That is for all s = 0, 1, . . ., d − 1, With this choice of the τ j 's, (18) reduces to (12): Imposing the initial condition ( 9) we get: ] T , we obtain ( 14) This completes the proof of the theorem.
The following corollary is a particular case of the previous theorem: Corollary 4 The assumptions of the previous theorem hold.If further d = 0 then where Proof.Since d = 0, all the canonical elements are defined and therefore the result is obtained from Theorem 1.
Example 2. Let us solve the initial value problem where A(x) is given in Example 1 and The exact solution is: Here d = 1 and therefore 0 } is the only undefined block.So we choose a perturbation term of the form: We applied the algorithm presented in Theorem 3 with N = 10 and V N = T N (x), the Chebyshev polynomial.We found that .
Figure 1 displays the exact errors in y 1 , y 2 and y 3 while Figure 2 shows the three components of the perturbation term H N (x).Throughout this section we considered H N as in (10).Following the same arguments we can derive analogous results when H N is of the form (11).It is worth to note that the Tau method with a perturbation of the form ( 11) is equivalent to the spectral collocation method at the zeros of V N (x), see (El-Daou & Ortiz, 1994).This equivalence permits to solve nonlinear differential equations using the spectral approach more effectively than the recursive Tau.In the next section we recall the main features of the spectral Tau method and we illustrate it by solving an optimal control problem with constraints given as a system of nonlinear differential equations.
To illustrate the spectral procedure we shall consider the Hicks-Ray reactor problem where it is desired to minimize the quadrature cost functional subject to the nonlinear constraints The two states are denoted by X(t) and Y(t), the control is denoted by U(t) and all other parameters are constants.
Our procedure starts with replacing the nonlinear constraints ( 24)-( 25) by an infinite sequence of linear systems of ODEs of the form (1).This is accomplished by employing an iterative procedure described as follows: Let ϕ(X, Y) := Xe −r/Y and W(t) := e −r/Y .The Taylor's series expansions of the bivariate function ϕ(X, Y) near a given point (X 0 , Y 0 ) allows to write where Dropping O 2 and using ( 27) in ( 24) yields: Again using ( 27)) and the approximation 28)-( 29) we obtain the sequence of linear systems which is of the form (1) with Clearly the coefficients {a ki , b ki , , i = 0, 1, . . ., N} of Xk (t) and Ỹk (t) will depend on {u ki ; i = 0, 1, . . ., N}.To determine the values of these u ki 's, we substitute Xk (t), Ỹk (t) and Ũk (t) in the objective function ( 23), and integrate it exactly to end up with the problem of minimizing a multivariate function of the form: min Ψ k (u k0 , u k1 , . . ., u kN ) that can be achieved using the direct approach.That is we solve the algebraic linear system formed by the gradient, ∂Ψ ∂u k j = 0, j = 0, 1, 2, . . ., N to obtain the unknowns {u k0 , u k1 , . . ., u kN }, and then we test the Hessian to verify the optimality.
For each k-cycle, we construct { Xk , Ỹk , Ũk } by the spectral Tau Algorithm, and repeat this process until a prescribed convergence tolerance ϵ is satisfied; that is until the iteration counter k reaches a certain k * such that We consider then { Xk * , Ỹk * , Ũk * } as the Tau approximation for {X, Y, U}.It is worth noting that the convergence of {( Xk , Ỹk , Ũk ), k ≥ 1} to the exact solution (X, Y, U) is guaranteed by Kantorovich Theorem that imposes conditions on the starting values { X0 , Ỹ0 , Ũ0 } and on the entries of the matrix A.
Figure 3 shows the profiles of the approximated states functions X(t) and Y(t) and the control U(t) computed by the spectral Tau method with N = 20.These were obtained when k * = 18 with tolerance ϵ = 10 −12 .This problem does not have an exact solution to compare but the same results can be obtained if this the problem is solved using the concept of Hamiltonian.
Table lists the coefficients of X20 , Ỹ20 and Ũ20 .The minimum value of the objective function = 2402.02746.