T Cluzeau
email: thomas.cluzeau@unilim.fr.

C Koutschan

A Quadrat
email: alban.quadrat@inria.fr.

M Tõnso

Effective algebraic analysis approach to linear systems over Ore algebras

Keywords: Linear systems theory, control theory, algebraic analysis, computer algebra, implementation T. Cluzeau, A. Quadrat

The purpose of this paper is to present a survey on the effective algebraic analysis approach to linear systems theory with applications to control theory and mathematical physics. In particular, we show how the combination of effective methods of computer algebrabased on Gröbner basis techniques over a class of noncommutative polynomial rings of functional operators called Ore algebrasand constructive aspects of module theory and homological algebra enables the characterization of structural properties of linear functional systems. Algorithms are given and a dedicated implementation, called OREALGEBRAIC-ANALYSIS, based on the Mathematica package HOLONOMICFUNCTIONS, is demonstrated.

Une approche par l'analyse algébrique effective des systèmes linéaires sur des algèbres de Ore Résumé : Le but de ce papier est de présenter un état de l'art d'une approche par l'analyse algébrique effective de la théorie des systèmes linéaires avec des applications à la théorie du contrôle et à la physique mathématique. En particulier, nous montrons comment la combinaison des méthodes effectives de calcul formelbasées sur les techniques de bases de Gröbner sur une classe d'algèbres polynomiales noncommutatives d'opérateurs fonctionnels appelée algèbres de Oreet d'aspects constructifs de théorie des modules et d'algèbre homologique permet la caractérisation de propriétés structurelles des systèmes linéaires fonctionnels. Des algorithmes sont donnés et une implémentation dédiée, appelée OREALGEBRAICANALYSIS, basée sur le package Mathematica HOLONOMIC-FUNCTIONS, est présentée.

Mots-clés : Théorie des systèmes linéaires, théorie du contrôle, analyse algébrique, calcul formel, implémentation

Introduction

To introduce the algebraic analysis approach to linear systems over Ore algebras, we use explicit examples. The model of a stirred tank studied in [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF] on page 7 is defined by the following mass balance equations

       d V (t) dt = -k V (t) S + F 1 (t) + F 2 (t), d (c(t) V (t)) dt = -c(t) k V (t) S + c 1 F 1 (t) + c 2 F 2 (t),
where F 1 and F 2 denote the flow rates of two incoming flows feeding the tank, c 1 and c 2 two constant concentrations of dissolved materials, c the concentration in the tank, V the volume, k an experimental constant, and S the constant cross-sectional area. The algebraic analysis approach can only handle linear systems. See [START_REF] Chakhar | An algebraic analysis approach to certain classes of nonlinear partial differential systems[END_REF] for a first attempt to extend the algebraic analysis approach to particular classes of nonlinear systems. We refer to [START_REF] Lange-Hegermann | Thomas Decomposition and Nonlinear Control Systems[END_REF] for the use of differential elimination techniques for studying this non linear system. If V 0 is a constant volume, c 0 a constant concentration, and

F 10 := (c 2 -c 0 ) (c 2 -c 1 ) k V 0 S , F 20 := (c 0 -c 1 ) (c 2 -c 1 ) k V 0 S , V (t) := V 0 + x 1 (t), c(t) := c 0 + x 2 (t), F 1 (t) := F 10 + u 1 (t), F 2 (t) := F 20 + u 2 (t),
then the linearized model around the steady-state equilibrium is defined by

       ẋ1 (t) = - 1 2 θ x 1 (t) + u 1 (t) + u 2 (t), ẋ2 (t) = - 1 θ x 2 (t) + c 1 -c 0 V 0 u 1 (t) + c 2 -c 0 V 0 u 2 (t), (1)
with the notation θ := V 0 /F 0 (the holdup time of the tank), where F 0 := k V 0 /S. See pages 8-9 of [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF]. The linear OD system (1) can then be studied by means of the standard analysis and synthesis techniques developed for linear OD systems. Now, if a transport delay of amplitude τ > 0 occurs in the pipe, then we obtain the following linear differential time-delay (DTD) system:

       ẋ1 (t) = - 1 2 θ x 1 (t) + u 1 (t) + u 2 (t), ẋ2 (t) = - 1 θ x 2 (t) + c 1 -c 0 V 0 u 1 (t -τ ) + c 2 -c 0 V 0 u 2 (t -τ ). (2)
For more details, see pages 449-451 of [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF]. Then, (2) can be studied by means of methods dedicated to linear DTD systems. Following [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF], if the valve settings are commanded by a process control computer which can only be changed at discrete instants and remain constant in between, the following discrete-time model of (1) can then be derived

     x 1 (n + 1) = e -∆ 2 θ x 1 (n) + 2 θ (1 -e -∆ 2 θ ) (u 1 (n) + u 2 (n)), x 2 (n + 1) = e -∆ θ x 2 (n) + θ (1 -e -∆ θ ) V 0 ((c 1 -c 0 ) u 1 (n) + (c 2 -c 0 ) u 2 (n)), (3)
where ∆ is the constant length of time intervals. For more details, see page 449 of [START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF]. Again, (3) can then be studied by means of standard techniques developed for linear discrete-time systems.

As shown above, a physical system can be modeled by means of different systems of functional equations, namely, systems whose unknowns are functions (e.g., OD systems, DTD systems, discrete-time systems). Moreover, the "same" system can be defined by means of different representations (e.g., state-space, input-output, polynomial, behaviors, geometric, systems over a ring, implicit, . . . representations). These representations are defined by different numbers of unknowns and equations. Linear systems are usually studied by means of dedicated mathematical methods which usually depend on the representations. The equivalences between different representations and different formulations of system-theoretic properties (e.g., controllability à la Kalman, controllability for polynomial systems, controllability à la Willems) are known for certain classes of linear functional systems.

We can wonder whether or not a unique mathematical approach to linear systems exists which satisfies the following two important requirements:

(a) The approach can handle the standard classes of linear functional systems studied in control theory by means of common mathematical concepts, methods, theorems, algorithms, and implementations.

(b) The approach does not depend on particular representations of the linear systems.

The goal of this paper is to show that the algebraic analysis approach satisfies these two points. Algebraic analysis (also called D-module theory) is a mathematical theory developed by B. Malgrange, J. Bernstein, M. Sato and his school in the sixties to study linear systems of partial differential (PD) equations by means of module theory, homological algebra, and sheaf theory (see [START_REF] Hotta | D-Modules, Perverse Sheaves, and Representation Theory[END_REF][START_REF] Kashiwara | Algebraic Study of Systems of Partial Differential Equations[END_REF][START_REF] Malgrange | Systèmes différentiels à coefficients constants[END_REF] and the references therein). In the nineties, algebraic analysis techniques were introduced in mathematical systems theory and control theory by U. Oberst, M. Fliess, and J.-F. Pommaret. For more details, see [START_REF] Fliess | Some basic structural properties of generalized linear systems[END_REF][START_REF] Fliess | Controllability and observability of linear delay systems: an algebraic approach[END_REF][START_REF] Oberst | Multidimensional constant linear systems[END_REF][START_REF] Pommaret | Partial Differential Control Theory[END_REF][START_REF] Pommaret | Algebraic analysis of linear multidimensional control systems[END_REF][START_REF] Zerz | An algebraic analysis approach to linear time-varying systems[END_REF] and the references therein.

Within the algebraic version of algebraic analysis, a linear system is studied by means of a finitely presented left module M [START_REF] Rotman | An Introduction to Homological Algebra[END_REF] over a ring D of functional operators, and its F-solutions are defined by the homomorphisms (namely, the left D-linear maps) from M to F, where F is a left D-module. We recall that a module is an algebraic structure which is defined by the same properties as the ones for a vector space but its scalars belong to a ring and not a field. Equivalent representations of a linear system yield isomorphic modules. These isomorphic modules are finitely presented by the different presentations, i.e., by the different matrices of functional operators defined by these representations. Hence, up to isomorphism, a linear system defines uniquely a finitely presented module. Structural (built-in) properties of linear systems, i.e., properties which do not depend on the representation of the system, then correspond to module properties (e.g., torsion elements, torsion-freeness, projectiveness, freeness). To study these module properties, we use homological algebra methods since they depend only on the underlying modules (up to isomorphism) and not on the presentations of the modules, i.e., not on the representations of the linear systems. Therefore, we have a way to study structural properties of linear systems independently of their representations. A second benefit of using homological algebra techniques is that large classes of linear functional systems can be studied by means of the same techniques, results, and algorithms since the standard rings of functional operators share the same properties. Only the "arithmetic" of the functional operators can be different. Based on Gröbner or Janet basis techniques for classes of noncommutative polynomial rings of functional operators, effective studies of module theory and homological algebra have recently been developed (see [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF][START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF][START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications. Les cours du CIRM[END_REF][START_REF] Robertz | Recent progress in an algebraic analysis approach to linear systems[END_REF] and the references therein). Dedicated symbolic packages such as OREMODULES, OREMORPHISMS and CLIPS have been developed [START_REF] Chyzak | OREMODULES: A symbolic package for the study of multidimensional linear systems[END_REF][START_REF] Cluzeau | OREMORPHISMS: A homological algebraic package for factoring, reducing and decomposing linear functional systems[END_REF][START_REF] Zerz | Algebraic systems theory and computer algebraic methods for some classes of linear control systems[END_REF].

The purpose of this paper is two-fold. We first give a brief overview of the algebraic analysis approach to linear systems defined over Ore algebras. We then show how a recent implementation of Gröbner bases for large classes of Ore algebras in a Mathematica package called HOLONOMICFUNCTIONS [START_REF] Koutschan | Advanced applications of the holonomic systems approach[END_REF][START_REF] Koutschan | HolonomicFunctions (user's guide)[END_REF] can be used to extend the classes of linear functional systems we can effectively study within the algebraic analysis approach. In particular, using the recent OREALGEBRAICANALYSIS package, we can now handle generic linearizations of explicit nonlinear functional systems or linear systems containing transcendental function (e.g., sin, cos, tanh) or special function coefficients (e.g., Airy or Bessel functions). These classes could not be studied by the OREMODULES, OREMORPHISMS or CLIP packages.

The paper is organized as follows. In Section 2, we explain that standard linear functional systems encountered in control theory can be studied by means of a polynomial approach over Ore algebras of functional operators, i.e., over a certain class of noncommutative polynomial rings. In Section 3, we shortly explain the concept of a Gröbner basis for left ideals and left modules over certain Ore algebras, and give algorithms to compute kernel and left/right inverses of matrices with entries in these Ore algebras. In Section 4, we introduce the algebraic analysis approach to linear systems theory and, using homological algebra techniques, we explain that this approach is an intrinsic polynomial approach to linear systems theory and we characterize standard system-theoretic properties in terms of module properties and homological algebra concepts that are shown to be computable. Finally, in Section 5, these results are illustrated on explicit examples which are studied by means of the OREALGEBRAICANALYSIS package.

Inria

This package is based on the Mathematica package HOLONOMICFUNCTIONS which contains Gröbner basis techniques for general classes of Ore algebras.

Linear systems over Ore algebras

In this section, we introduce the concept of a skew polynomial ring, an Ore extension and an Ore algebra [START_REF] Mcconnell | Noncommutative Noetherian Rings[END_REF] which will play important roles in what follows. To motivate the abstract definitions, let us start with standard examples of functional operators. In his treatises on differential equations, G. Boole used the idea of representing a linear OD equation r i=0 a i y (i) (t) = 0, where a i ∈ R, by means of the operator

P := r i=0 a i d i dt i , where d dt y(t) := y (1) (t) = ẏ(t)
∀ a 1 , a 2 ∈ A, d dt (a 1 + a 2 ) = da 1 dt + da 2 dt , d dt (a 1 a 2 ) = da 1 dt a 2 + a 1 da 2 dt ,
such as, for instance, the ring (resp., field) k[t] (resp., k(t)) of polynomials (resp., rational functions) in t with coefficients in a field k (e.g., k = Q, R, C) or C ∞ (R), then we can define the set of all the OD operators of the form r i=0 a i ∂ i with a i ∈ A. This set inherits a ring structure if the composition of OD operators is still an OD operator, i.e., if we have

  m j=0 b j ∂ j   n i=0 a i ∂ i = l k=0 c k ∂ k , (4)
for a certain l and for some c k ∈ A. In particular, such an identity should hold for m = 1 and n = 0, i.e., the composition of the two operators b 1 ∂ and a 0 has to be an OD operator. Since operators are understood by their actions on functions, we get

∀ y ∈ A, (b 1 ∂ a 0 ) y = b 1 ∂ (a 0 y) = b 1 d dt (a 0 y) = b 1 a 0 dy dt + da 0 dt y = b 1 a 0 ∂ + da 0 dt y.
Hence, on the OD operator level, we have the following commutation rule:

∀ a ∈ A, ∂ a = a ∂ + ȧ. (5)
It can be shown below that this commutation rule is enough to define a ring structure on the set of all the OD operators with coefficients in A. Note that the above commutation rule shows that this ring is usually noncommutative apart from the case where ȧ = 0 for all a ∈ A, i.e., the case where A is a ring of constants.

If we consider the case of a time-delay operator S defined by S y(t) = y(t -h), where h > 0, then to understand S a as an operator, where a is an element of a difference ring A of functions, namely, a commutative ring A of functions of t equipped with the endomorphism a(t) ∈ A -→ a(t -h) ∈ A, we have to apply it to a function y. We get

(S a(t)) y(t) = S (a(t) y(t)) = a(t -h) y(t -h) = (a(t -h) S) y(t),
i.e., on the operator level, we have the following commutation rule:

S a(t) = a(t -h) S. (6)
We note that in ( 5) and ( 6) the "degree" in ∂ or in S is 1 in both sides of the equalities. More generally, we can consider an operator ∂ which satisfies ∀ a ∈ A, ∂ a = σ ∂ + δ, where 0 = σ, δ ∈ A, so that both sides of the above expression have degree 1 in ∂. Clearly, σ and δ depend on a, i.e., σ(a) and δ(a). If we want to define a ring formed by elements which can uniquely be represented as

r i=0 a i ∂ i , we must have ∀ a 1 , a 2 ∈ A, ∂ (a 1 + a 2 ) = σ(a 1 + a 2 ) ∂ + δ(a 1 + a 2 ) = ∂ a 1 + ∂ a 2 = σ(a 1 ) ∂ + δ(a 1 ) + σ(a 2 ) ∂ + δ(a 2 ) = (σ(a 1 ) + σ(a 2 )) ∂ + δ(a 1 ) + δ(a 2 ),
which yields the following identities:

σ(a 1 + a 2 ) = σ(a 1 ) + σ(a 2 ), δ(a 1 + a 2 ) = δ(a 1 ) + δ(a 2 ).
Similarly, using the associativity of operators, we obtain

∀ a 1 , a 2 ∈ A, ∂ (a 1 a 2 ) = σ(a 1 a 2 ) ∂ + δ(a 1 a 2 ) = (∂ a 1 ) a 2 = (σ(a 1 ) ∂ + δ(a 1 )) a 2 = σ(a 1 ) (σ(a 2 ) ∂ + δ(a 2 )) + δ(a 1 ) a 2 = σ(a 1 ) σ(a 2 ) ∂ + σ(a 1 ) δ(a 2 ) + δ(a 1 ) a 2 ,
which yields the following identities:

σ(a 1 a 2 ) = σ(a 1 ) σ(a 2 ), δ(a 1 a 2 ) = σ(a 1 ) δ(a 2 ) + δ(a 1 ) a 2 .
We also have that ∂ = ∂ 1 = σ(1) ∂ + δ(1), which yields:

σ(1) = 1, δ(1) = 0.
The conditions on σ show that σ is an endomorphism of the ring A and δ is called a σ-derivation (if σ = id A , we find again the above definition of a derivation).

The concept of an Ore extension of a ring A was introduced by Ore [START_REF] Ore | Theory of non-commutative polynomials[END_REF] in 1933 to develop a unified mathematical framework to represent linear functional operators such as differential operators, difference and shift operators, q-shift and q-differential operators, and many more. Nowadays, this concept is widely used to state results and algorithms about linear functional operators in a concise and general form. For applications of this framework, for instance, to the problem of factoring operators or creative telescoping, see [START_REF] Bronstein | An introduction to pseudo-linear algebra[END_REF][START_REF] Chyzak | An extension of Zeilberger's fast algorithm to general holonomic functions[END_REF] and the references therein.

Definition 1 ([16]

) Let A be a ring. An Ore extension O := A[∂;σ,δ] of A is the noncommutative ring formed by all polynomials of the form n i=0 a i ∂ i , where n ∈ N and a i ∈ A, obeying the following commutation rule

∀ a ∈ A, ∂ a = σ(a) ∂ + δ(a), (7)
where σ is an endomorphism of A, namely, σ :

A -→ A satisfies ∀ a, b ∈ A,      σ(1) = 1, σ(a + b) = σ(a) + σ(b), σ(a b) = σ(a) σ(b),
and δ is a σ-derivation of A, namely, δ : A -→ A satisfies:

∀ a, b ∈ A, δ(a + b) = δ(a) + δ(b), δ(a b) = σ(a) δ(b) + δ(a) b. (8)
The Ore extension A[∂;σ,δ] is also called a skew polynomial ring.

Inria

Let O := A[∂;σ,δ] be a skew polynomial ring, P := n i=0 a i ∂ i ∈ O, where a n = 0, and

Q := m i=0 b i ∂ i ∈
O, where b m = 0. If A is a domain, i.e., A does not contain non-trivial zero divisors, then we have

P Q = (a n ∂ n + • • • ) (b m ∂ m + • • • ) = a n σ n (b m ) ∂ n+m + • • • ,
where • • • represents lower degree terms. Moreover, if σ is injective, we can define the degree of P to be n and the degree of Q to be m since we have:

∀ P, Q ∈ O, deg ∂ (P Q) = deg ∂ (P ) + deg ∂ (Q). A skew polynomial ring A[∂;σ,δ] has the structure of an A -A-bimodule, namely, O has a left module structure defined by ∀ a ∈ A, ∀ P = r i=0 a i ∂ i ∈ O : a P = r i=0 (a a i ) ∂ i ,
and a right A-module structure defined by

∀ a ∈ A, ∀ P = r i=0 a i ∂ i ∈ O : P a = r i=0 a i ∂ i a,
and they satisfy the following associativity condition: (a) If (A, δ) is a differential ring, i.e., A is a ring and δ is a derivation of A, i.e., δ satisfies (8) with σ = id A , then we can define the skew polynomial ring A[∂;id A , δ] of OD operators with coefficients in A. Then, (7) yields [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenrings nach einem nulldimensionalen Polynomideal[END_REF]. For instance, if we consider again (1), then we can define the algebra

∀ a 1 , a 2 ∈ A, ∀ P ∈ O, (a 1 P ) a 2 = a 1 (P a 2 ).
O := Q(θ,c 0 , c 1 , c 2 , V 0 )[∂; id A , δ]
of OD operators with coefficients in the field A := Q(θ,c 0 , c 1 , c 2 , V 0 ) of rational functions in the system parameters θ, c 0 , c 1 , c 2 , and V 0 , where δ := d dt is the trivial derivation of A, i.e., δ(a) = 0 for all a ∈ A. Thus, [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenrings nach einem nulldimensionalen Polynomideal[END_REF] implies that ∂ a = a ∂ for all a ∈ A, which shows that O is a commutative polynomial ring. Then, (1) can be rewritten as R η = 0, where:

R :=    ∂ + 1 2 θ 0 -1 -1 0 ∂ + 1 θ - c 1 -c 0 V 0 - c 2 -c 0 V 0    ∈ O 2×4 , η :=      x 1 (t) x 2 (t) u 1 (t) u 2 (t)      .
If one of the parameters is now a smooth function of t, then δ is no more the trivial derivation of A := C ∞ (R), and thus O is then a noncommutative polynomial ring in ∂ with coefficients in A.

A simple example of a noncommutative polynomial ring of OD operators is given by O

:= R[x][∂;id,δ],
where δ :

= d dx is the standard derivation on R[x]. The error function erf(x) := 2 √ π
x 0 e -t 2 dt satisfies the following ODE:

∂ 2 + 2 x ∂ erf(x) = 0.
(b) If we consider the algebra A := Q(θ,c 0 , c 1 , c 2 , V 0 , ∆, n) and the endomorphism σ(a(n)) := a(n + 1) of A, then we can define the skew polynomial ring O := A[S;σ,0] of forward shift operators, which encodes the commutation rule S a(n) = a(n + 1) S for a ∈ A. Then, (3) can be written as R η = 0, where:

R := S -e -∆ 2 θ 0 -2 θ (1 -e -∆ 2 θ ) -2 θ (1 -e -∆ 2 θ ) 0 S -e -∆ θ -α (c 1 -c 0 ) -α (c 2 -c 0 ) ∈ O 2×4 , α := θ (1 -e -∆ θ ) V 0 , η := (x 1 (n) x 2 (n) u 1 (n) u 2 (n)) T .
Since no entry of R is a (rational) function of n, we can only consider the algebra

A := Q(θ,c 0 , c 1 , c 2 , V 0 , ∆)
and σ = id A . We then get S a = a S for all a ∈ A, i.e., the ring of shift operators with constant coefficients is commutative.

A simple example of a noncommutative polynomial ring of shift operators is Q[n][S;σ,0], where σ(a(n)) = a(n + 1) for all a ∈ Q[n]. The Gamma function Γ(z) := +∞ 0 t z-1 e -t dt for (z) > 0 satisfies the following recurrence relation:

(S -n) Γ(n) = 0. (c)
(∂) = ∂ so that σ( r i=0 a i (t) ∂ i ) = r i=0 a i (t -h) ∂ i .
In particular, we have S ∂ = σ(∂) S = ∂ S, i.e., the two operators ∂ and S commute. This last identity encodes the following identity:

(∂ S)(y(t)) = ∂(y(t -h)) = ẏ(t -h) = (S ∂)(y(t)). (9)
Then, (2) can be rewritten as R η = 0, where:

R :=    ∂ + 1 2 θ 0 -1 -1 0 ∂ + 1 θ - (c 1 -c 0 ) V 0 S - (c 2 -c 0 ) V 0 S    ∈ O 2×4 , η := (x 1 (t) x 2 (t) u 1 (t) u 2 (t)) T .
(e) If we consider the difference (resp., divided difference) operator

a(t) -→ a(t + 1) -a(t) resp., a(t) -→ a(t) -a(t 0 ) t -t 0 ,
for a fixed t 0 ∈ R and for all a belonging to a field A of real-valued functions of t, then we can form the skew polynomial ring A[∂;σ,δ] of difference (resp., divided difference) operators with coefficients in A by respectively considering:

∀ a ∈ A, σ(a(t)) = a(t + 1), δ(a(t)) = a(t + 1) -a(t),    σ(a(t)) = a(t 0 ), δ(a(t)) = a(t) -a(t 0 ) t -t 0 .
If A is a (skew) field, then the right Euclidean division can be performed, i.e., the algebra O is a right Euclidean domain, and thus a principal left ideal domain, namely, every left ideal of O is finitely generated (see, e.g., [START_REF] Bronstein | An introduction to pseudo-linear algebra[END_REF][START_REF] Mcconnell | Noncommutative Noetherian Rings[END_REF]). Finally, if σ is also invertible, i.e., is an automorphism of A, then the left Euclidean division can also be performed, i.e., O is a left Euclidean domain, and thus a principal right ideal domain. More details on skew polynomial rings can be found in [START_REF] Mcconnell | Noncommutative Noetherian Rings[END_REF]. A left and right Euclidean domain is simply called a Euclidean domain.

A[∂ 1 ; σ 1 , δ 1 ] • • • [∂ m ; σ m , δ m ] := (• • • ((A[∂ 1 ; σ 1 , δ 1 ])[∂ 2 ; σ 2 , δ 2 ]) • • • )[∂ m ; σ m , δ m ]. If B := A[∂ 1 ; σ 1 , δ 1 ] • • • [∂ m-1 ; σ m-1 , δ m-1 ], then O := B[∂ m ; σ m , δ m ],
where σ m is an endomorphism of B and δ m is a σ m -derivation of B. In particular, we get:

∀ i = 1, . . . , m -1, ∀ a ∈ A, ∂ m ∂ i = σ m (∂ i ) ∂ m + δ m (∂ i ), ∂ m a = σ m (a) ∂ m + δ m (a).
Similarly, we have:

1 i < j m, ∀ a ∈ A, ∂ j ∂ i = σ j (∂ i ) ∂ j + δ j (∂ i ) ∂ j a = σ j (a) ∂ j + δ j (a). (10)
If we want that ∂ j commutes with ∂ i , σ j and δ j must satisfy the conditions:

1 i < j m, σ j (∂ i ) = ∂ i , δ j (∂ i ) = 0. (11)
Moreover, let us assume that σ j (A) ⊆ A and δ j (A) ⊆ A. Then, we have:

∂ j (∂ i a) = ∂ j (σ i (a) ∂ i + δ i (a)) = σ j (σ i (a) ∂ i ) ∂ j + δ j (σ i (a) ∂ i ) + σ j (δ i (a)) ∂ j + δ j (δ i (a)) = σ j (σ i (a)) σ j (∂ i ) ∂ j + σ j (σ i (a)) δ j (∂ i ) + δ j (σ i (a)) ∂ i + σ j (δ i (a)) ∂ j + δ j (δ i (a))
Using [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF], the above identity reduces to:

∂ j (∂ i a) = σ j (σ i (a)) ∂ i ∂ j + δ j (σ i (a)) ∂ i + σ j (δ i (a)) ∂ j + δ j (δ i (a)).
Since σ j (a) ∈ A and δ j (a) ∈ A, we also have:

∂ i (∂ j a) = ∂ i (σ j (a) ∂ j + δ j (a)) = σ i (σ j (a)) ∂ i ∂ j + δ i (σ j (a)) ∂ j + σ i (δ j (a)) ∂ i + δ i (δ j (a)).
If we have σ j (σ i (a)) = σ i (σ j (a)), δ j (σ i (a)) = σ i (δ j (a)), σ j (δ i (a)) = δ i (σ j (a)), and δ j (δ i (a)) = δ i (δ j (a)) for all a ∈ A, then we get

∂ j ∂ i a = ∂ i ∂ j a for all a ∈ A. Definition 2 Let k be a field. If A is a k-algebra, then an Ore extension of A of the form A[∂ 1 ; σ 1 , δ 1 ] • • • [∂ m ; σ m , δ m ]
is called an Ore algebra if σ j (A) ⊆ A and δ j (A) ⊆ A for j = 1, . . . , m, and:

1 i < j m, σ j (∂ i ) = ∂ i , δ j (∂ i ) = 0, 1 i, j m, i = j,      (σ j • σ i ) |A = (σ i • σ j ) |A , (δ j • σ i ) |A = (σ i • δ j ) |A , (δ j • δ i ) |A = (δ i • δ j ) |A .
We then have ∂ j ∂ i a = ∂ i ∂ j a for 1 i < j m and for all a ∈ A.

Finally, an Ore algebra

A[∂ 1 ; σ 1 , δ 1 ] • • • [∂ m ; σ m , δ m ] with coefficient ring A := k[x 1 , . . . ,
x n ] (resp., A := k(x 1 , . . . , x n )) is called a polynomial (resp., rational) Ore algebra.

If O := A[∂ 1 ; σ 1 , δ 1 ] • • • [∂ m ; σ m , δ m ]
is an Ore extension of a ring A, then P ∈ O can be expressed as P = 0 |ν| r p ν ∂ ν , where r ∈ N,

p ν ∈ A, ν := (ν 1 . . . ν m ) T ∈ N m , |ν| := ν 1 + • • • + ν m , and ∂ ν := ∂ ν1 1 • • • ∂ νm m .
Example 2 (a) If A := k[x 1 , . . . , x n ] (resp. A := k(x 1 , . . . , x n )), then the Ore algebra

O := A[∂ 1 ; σ 1 , δ 1 ] • • • [∂ n ; σ n , δ n ],
where σ i := id and

δ i := ∂ ∂xi for i = 1, . . . ,
d dt J n (t) = n t -1 J n (t) -J n+1 (t).
This equation can be rewritten as P J n (t) = 0, where

P := ∂ + S -n t -1 ∈ O. (c) If A is a k-algebra equipped with the following endomorphisms ∀ a ∈ A, σ k (a(i 1 , . . . , i n )) := a(i 1 , . . . , i k + 1, . . . , i n ), k = 1, . . . n, (e.g., A := k[i 1 , . . . , i n ], k(i 1 , . . . , i n ), or the algebra of real-valued sequences in (i 1 , . . . , i n ) ∈ Z n ), then A[S 1 ; σ 1 , 0] • • • [S n ; σ n , 0]
is the Ore algebra of multi-shift operators with coefficients in A.

(d) The ring of differential time-varying delay operators with S y(t) = y(t-h(t)), where h is a smooth function satisfying h(t) < t for all t larger than or equal to a certain T > 0, does not usually form an Ore algebra since we have

(∂ S)(y(t)) = ∂ y(t -h(t)) = (1 -ḣ(t)) ẏ(t -h(t)) = (1 -ḣ(t)) (S ∂)(y(t)),
i.e., ∂ S = (1 -ḣ) S ∂. It is an Ore algebra if and only if h is a constant function and we find then again [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF]. In [START_REF] Quadrat | Algebraic analysis for the Ore extension ring of differential time-varying delay operators[END_REF], it is shown that the ring of differential time-varying delay operators can be defined as an Ore extension and its properties are studied in terms of the function h.

For more examples of Ore algebras of functional operators and their uses in combinatorics and in the study of special functions, see [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF] and the references therein.

Theorem 1 can be used to prove that the Ore algebras defined in Example 2 are both left and right noetherian domains. We say that they are noetherian domains.

Finally, we introduce the concept of an involution of a ring which will be used in Sections 3.2 and 4.2.

Definition 3 Let O be an Ore algebra over a base field k. An involution of O is an anti-automorphism of order two of O, i.e., a k-linear map θ : O -→ O satisfying:

∀P 1 , P 2 ∈ O, θ(P 1 P 2 ) = θ(P 2 ) θ(P 1 ), θ • θ = id O .
Let us give a few examples of involutions.

Example 3 (a) If O is a commutative ring (e.g., O := k[x 1 , . . . , x n ]), then θ = id O is an involution of O. (b) Let O := A n (k) the polynomial Weyl algebra over k. Then, an involution of O is defined by θ(x i ) := x i and θ(∂ i ) := -∂ i for i = 1, . . . , n. More generally, if O := A[∂ 1 ; id, δ 1 ] • • • [∂ n ; id, δ n ] is a ring of PD
operators with coefficients in the differential ring (A, {δ 1 , . . . , δ n }), where δ i := ∂ ∂xi for i = 1, . . . , n, then an involution θ of O is defined by:

∀ a ∈ A, θ(a) := a, θ(∂ i ) := -∂ i , i = 1, . . . , n.

Gröbner basis techniques

In Section 2, we explain how standard linear functional systems can be defined by means of matrices of functional operators, i.e., by means of matrices with entries in noncommutative polynomial rings such as skew polynomial rings, Ore extensions, or Ore algebras. The idea of studying linear functional systems by means of the algebraic properties of their representations is well-developed in the polynomial approach [START_REF] Kailath | Linear Systems[END_REF]. If the ring of functional operators is a Euclidean domain, then Smith normal forms [START_REF] Kailath | Linear Systems[END_REF] can be extended to this noncommutative framework by considering the so-called Jacobson normal forms. For more details, implementations, and applications of Jacobson normal forms, see [START_REF] Levandovskyy | Computing diagonal form and Jacobson normal form of a matrix using Gröbner bases[END_REF] and the references therein. If the ring of functional operators is not a Euclidean domain (e.g., if the ring is usually defined by more than one functional operators), then such normal forms do not exist. But the Euclidean algorithm of multivariate (noncommutative) polynomials can still be used if the set of monomials appearing in the polynomials can be ordered in a particular way. This idea yields the concept of a Gröbner basis for a set of polynomials (i.e., for an ideal) or for a matrix (i.e., for a module).

In the next sections, we will state algorithms for the study of built-in properties of linear functional systems. These algorithms will be based on elimination techniques such as Gröbner basis techniques over noncommutative Ore algebras. Before doing so, we first motivate their uses by an explicit example.

-ν ∆ u + c u + ∇ p = 0, ∇ . u = 0,
where u ∈ R n is the velocity, p the pressure, ν the viscosity, and c the reaction coefficient. For simplicity reasons, let us consider the special case n = 2, i.e.

     E 1 := -ν (∂ 2 x u 1 + ∂ 2 y u 1 ) + c u 1 + ∂ x p = 0, E 2 := -ν (∂ 2 x u 2 + ∂ 2 y u 2 ) + c u 2 + ∂ y p = 0, E 3 := ∂ x u 1 + ∂ y u 2 = 0, (12)
with the standard notations ∂ x := ∂ ∂x and ∂ y := ∂ ∂y . We can wonder if the pressure p satisfies a system of PDEs by itself, i.e., if the components u 1 and u 2 of the speed can be eliminated from the equations of ( 12) to get PDEs only on p. Differentiating E 1 (resp., E 2 ) with respect to x (resp., y), we first obtain:

∂ x E 1 = -ν ∂ x (∂ 2 x u 1 + ∂ 2 y u 1 ) + c ∂ x u 1 + ∂ 2 x p = 0, ∂ y E 2 = -ν ∂ y (∂ 2 x u 2 + ∂ 2 y u 2 ) + c ∂ y u 2 + ∂ 2 y p = 0.
Similarly, we have:

ν (∂ 2 x E 3 + ∂ 2 y E 3 ) = ν ∂ x (∂ 2 x u 1 + ∂ 2 y u 1 ) + ν ∂ y (∂ 2 x u 2 + ∂ 2 y u 2 ) = 0, -c E 3 = -c (∂ x u 1 + ∂ y u 2 ) = 0.
Adding all the new differential consequences of the equations of (12), we get

∂ x E 1 + ∂ y E 2 + ν (∂ 2 x E 3 + ∂ 2 y E 3 ) -c E 3 = ∂ 2 x p + ∂ 2 y p = 0,
i.e., [START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF] yields ∆ p = 0, where ∆ := ∂ 2 x + ∂ 2 y is the Laplacian operator. This is an important result in hydrodynamics: the pressure must satisfy ∆ p = 0.

Gröbner basis techniques can be used for automatically eliminating (if possible) fixed unknowns. To do that, we first have to recast the above computations within a polynomial framework. Let us first consider the commutative ∈ D, where r ∈ N, a µ ∈ Q(ν,c), and µ := (µ x µ y ) T ∈ N 2 . Then, ( 12) can be rewritten as R η = 0, where:

polynomial ring D := Q(ν,c)[∂ x , ∂ y ] of PD operators in ∂ x and ∂ y with coefficients in Q(ν,c).
R :=    -ν ∆ + c 0 ∂ x 0 -ν ∆ + c ∂ y ∂ x ∂ y 0    ∈ D 3×3 , η :=    u 1 u 2 p    .
Then, the above computations correspond to the following matrix computations

(∂ x ∂ y ν ∆ -c)    E 1 E 2 E 3    = ((∂ x ∂ y ν ∆ -c) R)    u 1 u 2 p    = ∆ p
and using the fact that ∆ p = (0 0 ∆) η, we obtain:

(0 0 ∆) = (∂ x ∂ y ν ∆ -c) R ∈ D 1×3 R := {µ R | µ ∈ D 1×3 }.
We note that the D-submodule D 1×3 R of D 1×3 is formed by all the D-linear combinations of the rows of R. These combinations correspond to all the linear differential consequences of the equations of [START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF]. Within the operator framework, the fact that the pressure satisfies ∆ p = 0 can be rewritten as

(0 0 ∆) ∈ D 1×3 R. If R ∈ D q×p , then the (left) D-submodule L := D 1×q R of D 1×p is generated by the rows of R.
If D is a (noncommutative) polynomial ring, then a Gröbner basis of L is another set of generators of L, i.e., we have L = D 1×q R for a certain matrix R ∈ D q ×p , for which the so-called membership problem can easily be checked. The membership problem aims at deciding whether or not λ ∈ D 1×p belongs to D 1×q R. If D is a commutative polynomial ring with coefficients in a computable field, then Buchberger's algorithm [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenrings nach einem nulldimensionalen Polynomideal[END_REF] computes a Gröbner basis for a fixed monomial order. This result can be extended for some classes of noncommutative polynomial rings where the algorithm is proved to terminate. If a Gröbner basis R of L is known, then we can reduce any λ ∈ D 1×p with respect to this Gröbner basis in a unique way, i.e., there exists a unique λ ∈ D 1×p , called the normal form of λ, such that λ = λ + µ R for a certain µ ∈ D 1×q . Hence, we obtain that λ ∈ L if and only if we have λ = 0.

In the next sections, we first define the concept of a Gröbner basis for a finitely generated left ideal and then for a finitely generated left module.

Gröbner bases for ideals over Ore algebras

We first explain the basics of Gröbner bases using the standard commutative setting, i.e., for the case of a polynomial ring in several commuting variables, and then shortly explain how the theory can be extended to noncommutative Ore algebras.

Let x := x 1 , . . . , x n be a collection of variables, and let us denote by k[x] the ring of multivariate polynomials in x 1 , . . . , x n with coefficients in the field k. For α ∈ N n , we define the monomial

x α := x α1 1 • • • x αn n .
Unlike for univariate polynomials, there is no natural ordering of the monomials x α in a multivariate polynomial α∈N n c α x α . This is the reason for introducing the notion of monomial order, that is a total order ≺ on the set {x α | α ∈ N n } of x-monomials, namely an order ≺ which is total (i.e., we have either

x α ≺ x β or x β ≺ x α for all α, β ∈ N n , α = β). Definition 4 A monomial order on the set {x α | α ∈ N n } of x-monomials is called admissible if it satisfies the following conditions: (a) 1 ≺ x α ∀ α ∈ N n \ {(0, . . . , 0)}, (b) x α ≺ x β =⇒ x α x γ ≺ x β x γ ∀ α, β, γ ∈ N n .
Inria It follows that the set of monomials is well-founded with respect to any admissible monomial order, i.e., that each strictly decreasing sequence of monomials is finite. This is a crucial property for proving the termination of Buchberger's algorithm which computes a Gröbner basis of a polynomial ideal.

Example 5 We identify a monomial x α with the multi-index α ∈ N n .

(a) The lexicographic order on x-monomials is defined by α ≺ lex β whenever the first nonzero entry of β -α is positive. For instance, if we consider Q[x 1 , x 2 , x 3 ], then we have:

1 ≺ lex x 3 ≺ lex x 2 3 ≺ lex x 2 ≺ lex x 2 x 3 ≺ lex x 2 2 ≺ lex x 1 ≺ lex x 1 x 3 ≺ lex x 1 x 2 ≺ lex x 2 1 .
(b) The total degree order (also called degree reverse lexicographic order or graded reverse lexicographic order) on x-monomials is defined by α ≺ tdeg β whenever |α| < |β| or if we have |α| = |β|, then the last nonzero

entry of β -α is negative. It is also denoted ≺ degrevlex . For instance, if we consider Q[x 1 , x 2 , x 3 ], then we have: 1 ≺ tdeg x 3 ≺ tdeg x 2 ≺ tdeg x 1 ≺ tdeg x 2 3 ≺ tdeg x 2 x 3 ≺ tdeg x 1 x 3 ≺ tdeg x 2 2 ≺ tdeg x 1 x 2 ≺ tdeg x 2 1 .
(c) Let x := x 1 , . . . , x n and y := y 1 , . . . , y m be two collections of variables. Assume that an admissible monomial order ≺ X (resp., ≺ Y ) on x-monomials (resp., on y-monomials) is given. An elimination order is then defined by

u v ≺ w t ⇐⇒ u ≺ X w or u = w and v ≺ Y t,
where u, w (resp., v, t) are x-monomials (resp., y-monomials). An elimination order serves to eliminate the x i 's. The elimination order, which will be used in what follows, is the one induced by the total degree orders on x-monomials and y-monomials. This is a very common order called lexdeg. For instance, if we

consider Q[x 1 , x 2 , x 3 ], x = x 1 , x 2 , y = x 3 , ≺ X =≺ tdeg
and ≺ Y =≺ tdeg , then we have:

1 ≺ lexdeg x 3 ≺ lexdeg x 2 3 ≺ lexdeg x 2 ≺ lexdeg x 2 x 3 ≺ lexdeg x 1 ≺ lexdeg x 1 x 3 ≺ lexdeg x 2 2 ≺ lexdeg x 1 x 2 ≺ lexdeg x 2 1 .
Definition 5 Let P ∈ k[x] \ {0} and ≺ be an admissible monomial order. We can then define:

• The leading monomial lm ≺ (P ) of P to be the ≺-maximal monomial that appears in P with nonzero coefficient.

• The leading coefficient lc ≺ (P ) of P to be the coefficient of lm ≺ (P ).

• The leading term lt ≺ (P ) of P to be the product lc ≺ (P ) lm ≺ (P ).

When no confusion can arise, we skip the explicit mentioning of the monomial order in the subscripts. Hence, we can write P = lc(P ) lm(P ) + Q = lt(P ) + Q, where all monomials in the expanded expression of Q are strictly smaller (with respect to the chosen monomial order) than lm(P ).

Next, the concept of polynomial reduction is introduced, also called multivariate polynomial division, as it generalizes Euclidean division of univariate polynomials. For this purpose, we fix an admissible monomial order ≺ and use it in the following without any explicit mentioning. For nonzero polynomials P, Q ∈ k[x], one says that P is reducible by Q if lm(P ) is divisible by lm(Q). In other words, one can reduce P with respect to Q, and the result of the reduction is denoted by

red ≺ (P, Q) = red(P, Q) := P - lt(P ) lt(Q) Q.
It is important to notice that red(P, Q) = 0 or lm red(P, Q) ≺ lm(P ).

If G := {G 1 , . . . , G s } ⊆ k[x] \ {0} is
a set of polynomials, then red(P, G) denotes a polynomial obtained by iteratively reducing P with some elements of G until no such reduction is possible any more, i.e., the result is irreducible with respect to all elements of G.

Note that red(P, G) is usually not uniquely defined since it may depend on the choice of the polynomial G i that is used in a certain reduction step as demonstrated in the following example.

Example 6 Let us consider Q[x 1 , x 2 ] endowed with a total degree order (see (b) of Example 4). Choosing G :=

{G 1 , G 2 } with G 1 := x 1 x 2 -1 and G 2 := x 2 1 + x 2 + 1, the monomial x 2 1
x 2 can be reduced in two different ways yielding the two different irreducible polynomials

x 2 1 x 2 -x 1 G 1 = x 1 and x 2 1 x 2 -x 2 G 2 = -x 2 2 -x 2 . Definition 6 Let G denote the ideal generated by G 1 , . . . , G s ∈ k[x], i.e.: G = G 1 , . . . , G s := P 1 G 1 + • • • + P s G s | P 1 , . . . , P s ∈ k[x] .
Then, G is called a Gröbner basis with respect to the admissible monomial order ≺ if and only if one of the following equivalent statements holds:

(a) P ∈ G if and only if red ≺ (P, G) = 0. (b) red ≺ (P, G) is unique for any P ∈ k[x]. (c) If P ∈ G \ {0} then there exists G i ∈ G such that lm ≺ (G i ) divides lm ≺ (P ). (d) lm ≺ (P ) | P ∈ G \ {0} = lm ≺ (G 1 ), . . . ,
G := {G 1 , G 2 , G 3 }, with G 3 := x 2 2 + x 1 + x 2 . Since G 3 = x 2 G 2 -x 1 G 1 , we have G = G 1 , G 2 .
We claim that G is a Gröbner basis (see below for an algorithm which computes a Gröbner basis). Now, the monomial

x 2 1 x 2 reduces to x 1 since the polynomial -x 2 2 -x 2 = red(x 2 1 x 2 , G 2 ) is now reducible by G 3 yielding red(-x 2 2 -x 2 , G 3 ) = x 1 .
No further reductions can be done. Hence, we obtain

x 2 1 x 2 ∈ G .
Let us now shortly explain the principle of Buchberger's algorithm for computing a Gröbner basis of a polynomial ideal. Let lcm(m 1 , m 2 ) denote the least common multiple of the two monomials m 1 and m 2 . Buchberger's algorithm is based on the computation of the so-called S-polynomials. Definition 7 Given P, Q ∈ k[x] \ {0} and a monomial order ≺, we can define the S-polynomial S(P, Q) by:

S(P, Q) := lcm(lm(P ), lm(Q)) lt(P ) P - lcm(lm(P ), lm(Q)) lt(Q) Q.
Given a finite set {P 1 , . . . , P r } of elements of k[x] and an admissible monomial order ≺ on x-monomials, Buchberger's algorithm, which computes a Gröbner basis G := {Q 1 , . . . , Q s } of the ideal P 1 , . . . , P r of k[x], can be sketched as follows:

(a) Set G := {P 1 , . . . , P r } and let P be the set of pairs of distinct elements of G;

(b) While P = ∅, do:

• Choose (P i , P j ) ∈ P and remove it from P;

• Compute S(P i , P j ) and its reduction R ij := red ≺ (S(P i , P j ), G) by G;

• If R ij = 0, then: -Add {(P, R ij ) | P ∈ G} to P; -Add R ij to G; (c) Return G.
One can prove that the latter process terminates with a Gröbner basis G of the ideal P 1 , . . . , P r . For more details, we refer to [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenrings nach einem nulldimensionalen Polynomideal[END_REF][START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF][START_REF] Greuel | A Singular introduction to commutative algebra[END_REF][START_REF] Kandri-Rody | Non-commutative Gröbner bases in algebras of solvable type[END_REF]. While an ideal admits many different Gröbner bases with respect to the same monomial order, one can achieve uniqueness by means of the following definition. Inria Definition 8 A Gröbner basis G := {G 1 , . . . , G s } is said to be reduced if it satisfies the following two conditions:

• lc(G i ) = 1 for i = 1, . . . , s.

• Each monomial in G i is irreducible with respect to G \ {G i } for all i = 1, . . . , s.

Example 8

The Gröbner basis in Example 7 is a reduced one.

Remark 2 For an ideal of k[x] defined by a finite set of generators and a given monomial order ≺, one can compute a Gröbner basis, using, e.g., Buchberger's algorithm [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenrings nach einem nulldimensionalen Polynomideal[END_REF]. Algorithms for computing Gröbner bases are implemented in most of the computer algebra systems such as Maple, Mathematica, and Magma, or in dedicated computer algebra systems such as Singular and Macaulay2. However, in practice, such computations can be very costly, and it is still a topic of ongoing research to design faster algorithms for computing Gröbner bases. See the recent survey article [START_REF] Eder | A survey on signature-based Gröbner basis computations[END_REF] and the references therein.

Let us shortly state a few applications of Gröbner bases. Using the concept of a reduced Gröbner basis, we obtain a procedure to test whether or not two ideals of a commutative polynomial ring over a field, defined by different sets of generators, are equal: we check whether or not they have the same reduced Gröbner basis.

Solving a system of polynomial equations is an important application of Gröbner bases. For this purpose, we use the lexicographic order (see (a) of Example 5) which leads to a reduced Gröbner basis of a special form called "triangular" form. This means that some of the polynomials of the Gröbner basis depend only on certain variables, which simplifies the process of finding all solutions of the original system.

Example 9 Let G 1 , G 2 ∈ Q[x 1 , x 2 ]
be as in Example 6 but now endowed with the lexicographic order (see (a) of Example 5). Then, {x 3 2

+ x 2 2 + 1, x 1 + x 2 2 + x 2
} is a Gröbner basis with respect to this monomial order. Note that this Gröbner basis has a triangular form: the first element depends only on x 2 . The solutions of the polynomial system G 1 = G 2 = 0 can be obtained by first solving x 3 2 + x 2 2 + 1 = 0 and then plugging the solutions for

x 2 into x 1 = -(x 2 2 + x 2 ).
Gröbner basis techniques can also be used to develop an elimination theory. Let us state a standard problem for ideals: if I ⊆ k[x] is an ideal and y is a subset of x, then compute generators for the ideal I ∩ k[y]. To do that, we use the monomial order defined in (c) of Example 5. As explained in Section 4, elimination techniques play an important role in the effective study of module theory and homological algebra.

Example 10 If we consider again Example 9, we can check that we have:

G 1 , G 2 ∩ Q[x 2 ] = x 3 2 + x 2 2 + 1 .
The theory of Gröbner bases has been extended to noncommutative polynomial rings. See the work of Bergman [START_REF] Bergman | The Diamond lemma for ring theory[END_REF] for a very general and theoretic approach. A more algorithmically oriented but less general approach was presented in [START_REF] Kandri-Rody | Non-commutative Gröbner bases in algebras of solvable type[END_REF]. It only considers the so-called rings of solvable type (see also [START_REF] Kredel | Solvable Polynomial Rings[END_REF]). However, for our purposes, the latter suffices as most of the Ore algebras of interest are of solvable type. In this setting, again Buchberger's algorithm can be used to compute Gröbner bases, with only slight modifications due to noncommutativity.

Theorem 2 ( [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF][START_REF] Kredel | Solvable Polynomial Rings[END_REF]) Let k be a field, A := k[x 1 , . . . , x n ] the polynomial ring with coefficients in k, and

O := A[∂ 1 ; σ 1 , δ 1 ] • • • [∂ m ; σ m , δ m ]
a polynomial Ore algebra satisfying the following conditions

σ i (x j ) = a ij x j + b ij , δ i (x j ) = c ij , 1 i m, 1 j n, for certain a ij ∈ k \ {0}, b ij ∈ k, and c ij ∈ A.
Let ≺ be an admissible monomial order on the following set of monomials:

Mon(O) := {x α1 1 • • • x αm m ∂ ν1 1 • • • ∂ νn n | (α 1 , . . . , α m ) ∈ N m , (ν 1 , . . . , ν n ) ∈ N n }.
If the ≺-greatest term u in each non-zero c ij satisfies u ≺ x j ∂ i , then given a set of noncommutative polynomials in O, a noncommutative version of Buchberger's algorithm terminates for this admissible monomial order and its result is a Gröbner basis with respect to this order.

For more general results, we refer the reader to [START_REF] Kandri-Rody | Non-commutative Gröbner bases in algebras of solvable type[END_REF][START_REF] Kredel | Solvable Polynomial Rings[END_REF][START_REF] Levandovskyy | Non-commutative Computer Algebra for Polynomial Algebras: Gröbner Bases[END_REF]. In particular, for the Weyl algebra A n (Q) (see (c) of Example 3), the existence of Gröbner bases and the generalization of Buchberger's algorithm have been studied, e.g., in [START_REF] Levandovskyy | Plural -a computer algebra system for noncommutative polynomial algebras[END_REF][START_REF]Macaulay 2 project[END_REF][START_REF] Takayama | A system for computation in algebraic analysis[END_REF].

Example 11 Let us consider O := B 2 (Q) and the following linear PD system:

∂ 2 1 y = 0, x 1 ∂ 2 y + x 2 y = 0. (13)
Applying ∂ 1 to the second equation of ( 13), we get

x 1 ∂ 1 ∂ 2 y + ∂ 2 y + x 2 ∂ 1 y = 0. Applying again ∂ 1 to the equation then yields x 1 ∂ 2 1 ∂ 2 y + 2 ∂ 1 ∂ 2 y + x 2 ∂ 2
1 y = 0 and using (13), we get ∂ 1 ∂ 2 y = 0, and thus ∂ 2 y + x 2 ∂ 1 y = 0. Eliminating ∂ 2 y from the last equation by means of the second equation of ( 13), we obtain x 1 ∂ 1 y -y = 0. If we now apply ∂ 2 to the latter equation and use ∂ 1 ∂ 2 y = 0, we obtain ∂ 2 y = 0, which by substitution in the second equation of [START_REF] Cluzeau | OREMORPHISMS: A homological algebraic package for factoring, reducing and decomposing linear functional systems[END_REF] gives y = 0. The solution of ( 13) is then y = 0, a fact which is not obvious from [START_REF] Cluzeau | OREMORPHISMS: A homological algebraic package for factoring, reducing and decomposing linear functional systems[END_REF]. The computation of a Gröbner basis for the left O-ideal

I := O∂ 2 1 + O(x 1 ∂ 2 + x 2 )
for the total degree order follows the same line and yields I = O.

Gröbner bases for modules over Ore algebras

We now explain how we can extend the concept of a Gröbner basis from finitely generated left ideals to finitely generated left modules over an Ore algebra O. Let us first state again the definition of a module.

Definition 9 Let D be a noncommutative ring. A left D-module M is an abelian group (M, +) equipped with a scalar multiplication

D × M -→ M (d, m) -→ d m,
which satisfies the following properties

(a) d 1 (m 1 + m 2 ) = d 1 m 1 + d 1 m 2 , (b) (d 1 + d 2 ) m 1 = d 1 m 1 + d 2 m 1 , (c) (d 2 d 1 ) m 1 = d 2 (d 1 m 1 ), (d) 1 m 1 = m 1 , for all d 1 , d 2 ∈ D and for all m 1 , m 2 ∈ M .
Remark 3 The definition of a left D-module is similar to the one of a vector space but where the scalars belong to a noncommutative ring D and not to a (skew) field (e.g., Q, R, C) as for vector spaces.

A left D-module M is said to be finitely generated if M admits a finite set of generators, namely there exists a finite set S := {m i } i=1,...,r of elements of M such that for every m ∈ M , there exist d i ∈ D for i = 1, . . . , r such that:

m = r i=1 d i m i . S is
u f i ≺ v f j =⇒ w u f i ≺ w v f j . (b) u ≺ v =⇒ u f j ≺ v f j for j = 1, . . . , p. Inria Without loss of generality, we let f p ≺ f p-1 ≺ • • • ≺ f 1 .
There are two natural extensions of a monomial order to Mon(O 1×p ).

u f i ≺ v f j ⇐⇒ u ≺ v or u = v and f i ≺ f j .
(b) The position over term order on Mon(O 1×p ) induced by ≺ is defined by:

u f i ≺ v f j ⇐⇒ f i ≺ f j or f i = f j and u ≺ v.

Remark 4

The term over position order is of more computational value with regard to efficiency. The position over term order can be used to eliminate components.

If an admissible monomial order on Mon(O 1×p ) is fixed, then leading monomials and leading coefficients in O 1×p are defined similarly as in the case of ideals. Let R ∈ O q×p and L := O 1×q R be the left O-submodule of O 1×p . Buchberger's algorithm carries over to L. For more details, we refer, e.g., to [START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] Greuel | A Singular introduction to commutative algebra[END_REF].

:= Q(ν,c)[∂ x , ∂ y ]- submodule L := O 1×3 R of O 1×3 , i.e., (-ν ∆ + c) f 1 + ∂ x f 3 , (-ν ∆ + c) f 2 + ∂ y f 3 , ∂ x f 1 + ∂ y f 2 ,
for the position over term order induced by the monomial order ≺ tdeg (see (b) of Example 5). The Gröbner basis of L is then given by:

∂ x f 1 + ∂ y f 2 , (-ν ∂ 2 y + c) f 1 + ν ∂ x ∂ y f 2 + ∂ x f 3 , ∆ f 3 , (-ν ∆ + c) f 2 + ∂ y f 3 .
We find again that the pressure p satisfies ∆ p = 0 as shown in Example 4.

Let us shortly explain how Gröbner basis techniques can be used to compute left kernels (syzygy module computation), left factorizations and left inverses, . . . of matrices with entries in O. For more details, we refer to [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF].

P :=    p j=1 R ij η j -ζ i | i = 1, . . . , q    .
(b) Compute a Gröbner basis G of P in the free left O-module generated by the η j 's and the ζ i 's for j = 1, . . . , p and i = 1, . . . , q, namely, p j=1 Oη j ⊕ q i=1 Oζ i , with respect to a term order which eliminates the η j 's

(see (c) of Example 5). (c) Compute G ∩ ( q i=1 Oζ i ) = { q i=1 S ki ζ i | k = 1, . . . ,
P :=    p j=1 R ij η j -ζ i | i = 1, . . . , q    .
(b) Compute a Gröbner basis G of P in the free left O-module generated by the η j 's and the ζ i 's for j = 1, . . . , p and i = 1, . . . , q , namely, p j=1 Oη j ⊕ q i=1 Oζ i with respect to a term order which eliminates the η j 's (see (c) of Example 5).

(c) Define the following set:

Q :=    p j=1 R kj η j | k = 1, . . . , q    . (d) Compute the reduction H i of each element Q i of Q by G.
(e) If one of the H i 's contains η j , i.e., if the normal form of Q i contains not only ζ i 's, then return [ ], else return

R := (R ij ) ∈ O q×q , where H i = q j=1 R ij ζ j for i = 1, . . . , q.
P :=    p j=1 R ij η j -ζ i | i = 1, . . . , q    .
(b) Compute a Gröbner basis G of P in the free left O-module generated by the η j 's and the ζ i 's for j = 1, . . . , p and i = 1, . . . , q, namely, p j=1 Oη j ⊕ q i=1 Oζ i , with respect to a term order which eliminates the η j 's (see (c) of Example 5).

(c) Remove from G the elements which do not contain any η i and call H this new set.

(d) Write H in the form Q 1 (η 1 . . . η p ) T -Q 2 (ζ 1 . . . ζ q ) T , where Q 1 and Q 2 are two matrices with entries in O. (e) If Q 1 is invertible over O, then return S := Q -1 1 Q 2 ∈ O p×q , else return [ ].
Right analogues of the above algorithms (i. For an implementation of these algorithms in a computer algebra system, see [START_REF] Chyzak | OREMODULES: A symbolic package for the study of multidimensional linear systems[END_REF].

Inria 4 Algebraic analysis approach to linear systems theory

Linear functional systems and finitely presented left modules

As explained in Section 1, we study linear functional systems of the form R η = 0, where R ∈ D q×p , D is a noetherian domain (e.g., a noetherian Ore algebra O of functional operators (see Section 2)), and η is a vector of unknown functions. More precisely, if F is a left D-module (see Definition 9), then we can consider the following linear system or behavior:

ker F (R.) := {η ∈ F p | R η = 0}.
See Example 1 for the different models of the stirred tank considered in Section 1.

Remark 5

In this framework, we can consider the following classes of systems:

• State-space/input-output representation of 1-D linear systems. Considering, e.g.,

R := (∂ I n -A -B) ∈ O n×(n+m) , η := (x(t) T u(t) T ) T ∈ F n+m , R := (P (∂) -Q(∂)) ∈ O q×(q+r) , η := (y(t) T u(t) T ) T ∈ F q+r , where O := A ∂; id A , d
dt is a ring of OD operators with coefficients in a differential ring A and P has full row rank (i.e., ker O (.P ) = 0), we obtain the linear systems ẋ(t) = A(t) x(t) + B(t) u(t) and

P (∂) y(t) = Q(∂) u(t).
Similarly, we can consider the Ore algebra O := A[S;σ,0] of shift operators with coefficients in the difference ring A and S instead of ∂ in the above matrices to get the linear systems

x k+1 = A k x k + B k u k and P (S) y k = Q(S) u k .
• In the first above example, if we consider the Ore algebra O := A ∂; id A , d dt , where A := B[S;σ,0] and B is a difference ring, then we obtain the system ẋ(t) = A(t, S) x(t) + B(t, S) u(t), called in the literature a system over ring. Note that a general linear differential constant time-delay system is defined by R η = 0, where R ∈ O q×p , η ∈ F p and, e.g., F = C ∞ (R 0 ).

• General linear nD systems can be defined by R η = 0, where R ∈ O q×p and O is, for instance, one of the Ore algebras considered in Example 2. For instance, a simple discrete Roesser model can be defined by R η = 0, where Linear systems (e.g., a linearization of a nonlinear system around a given solution) can be studied within the algebraic analysis approach. The next example explains how the generic linearization of a nonlinear system can also be studied.

R := S 1 I r h -A 11 -A 12 -B 1 -A 21 S 2 I rv -A 22 -B 2 ∈ O (r h +rv)×(r h +rv+m) , η := (x T h x T v u T ) T , x h ∈ F r h , x v ∈ F rv , u ∈ F m ,

Example 13

We consider the nonlinear OD system defined by

ẋ(t) = f (x(t), u(t)), (14)
where we first suppose that f

= (f 1 • • • f n ) T
, where f i is a polynomial for i = 1, . . . , n. Let us denote X := X 1 , . . . , X n and U := U 1 , . . . , U m . Let k be a differential field (e.g., a field which is a differential ring), k{X, U } the differential ring formed by polynomials in a finite number of the X i 's, U j 's, and of their derivatives with coefficients in k, and p the differential ideal defined by the differential polynomials Ẋi -f i (X, U ) for i = 1, . . . , n, and their derivatives. Then, we can define the ring A := k{X, U }/p formed by the differential polynomials modulo the ideal p. If we denote by x i (resp., u j ) the residue class of X i (resp., U j ) in A, x := x 1 , . . . , x n , u := u 1 , . . . , u m , and A = k{x, u}, then these polynomials can be rewritten as polynomials in x i , u j , and the derivatives of the u j 's. Clearly, A is a differential ring with the derivation δ := d dt . It can be proved that p is a prime ideal, i.e., that A is an integral domain. Thus, we can define the quotient field K := Q(A) of A, i.e., the ring of fractions of A, which is a differential field for the derivation δ. Let O := B[∂;id A , δ] be the skew polynomial ring of OD operators with coefficients in B := A or K. The generic linearization of ( 14) is then defined by R η = 0, where R := ∂ I n -∂f ∂x -∂f ∂u ∈ D n×(n+m) and η := (dx T du T ) T , and can be studied by means of the finitely presented left O-module M := O 1×(n+m) /(O 1×n R). The cases of a rational, analytic or meromorphic function f can be studied similarly by considering the differential ring or field B formed by the rational/analytic/meromorphic functions which satisfy [START_REF] Cluzeau | Equivalences of linear functional systems[END_REF].

Within the algebraic analysis approach to linear systems theory [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF][START_REF] Oberst | Multidimensional constant linear systems[END_REF][START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications. Les cours du CIRM[END_REF][START_REF] Robertz | Recent progress in an algebraic analysis approach to linear systems[END_REF][START_REF] Zerz | An algebraic analysis approach to linear time-varying systems[END_REF], the linear system or behaviour is studied by means of the factor left D-module

M := D 1×p /(D 1×q R)
formed by the set of the residue classes π(λ) of λ ∈ D 1×p modulo the left D-submodule L := D 1×q R of D 1×p (i.e., π(λ) = π(λ ) if there exists µ ∈ D 1×q such that λ = λ + µ R) and equipped with the following left D-module structure: The left D-module M is said to be finitely presented and R is called a presentation matrix [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]. If {f j } j=1,...,p is the standard basis of D 1×p and y j := π(f j ) for j = 1, . . . , p, then {y j } j=1,...,p is a set of generators of the left D-module M (see Section 3.2). The generators y j 's of M satisfy non-trivial relations since we have:

∀ λ, λ ∈ D 1×p , ∀ d ∈ D, π(λ) + π(λ ) := π(λ + λ ), d π(λ) := π(d λ).
p j=1 R ij y j = 0, i = 1, . . . , q.
For the details of these results, see the chapter [START_REF] Cluzeau | Equivalences of linear functional systems[END_REF] of this book. Note that the y j 's do not belong to F but are just elements of M . To speak about F-solutions of R η = 0, we have to consider the homomorphisms from M to F, namely the maps f : M -→ F satisfying the following (left D-linear) condition:

∀ d 1 , d 2 ∈ D, ∀ m 1 , m 2 ∈ M, f (d 1 m 1 + d 2 m 2 ) = d 1 f (m 1 ) + d 2 f (m 2 ).
We recall that f ∈ hom D (M, F) is said to be an isomorphism if f is both injective and surjective [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]. If an isomorphism exists between M and F, then we say that M and F are isomorphic, which is denoted by M ∼ = F.

A standard result of homological algebra concerning the left exactness of the contravariant functor hom D ( • , F) [START_REF] Rotman | An Introduction to Homological Algebra[END_REF] yields the following fundamental result for the algebraic analysis approach of linear systems theory.

Theorem 3 We have the following isomorphism of abelian groups (i.e., Z-modules):

ker F (R.) ∼ = hom D (M, F). ( 15
)
For a direct proof of Theorem 3, see the chapter [START_REF] Cluzeau | Equivalences of linear functional systems[END_REF] of this book. ∈ ker F (R.), where F := C ∞ (R). But we have

R (∂ η) = ∂ 2 + ∂ (t -m) σ 2 - 1 σ 2 η = ∂ ∂ + t -m σ 2 η - 1 σ 2 η = - 1 σ 2 η, R (t η) = (t ∂ + 1) η + t (t -m) σ 2 η = t ∂ + t -m σ 2 η + η = η,
Inria which shows that neither η nor t η belongs to ker F (R.), i.e., ker F (R.) has no left D-module structure. However they are abelian groups (i.e., the Z-modules) and k-vector spaces if D is a k-algebra and k is a field included in the center of D:

Z(D) := {d ∈ D | d D = D d}.
If F := D, then hom D (M, D) inherits a right D-module structure [START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications. Les cours du CIRM[END_REF][START_REF] Rotman | An Introduction to Homological Algebra[END_REF].

Using the isomorphism (15), the linear system ker F (R.) depends only on M and F. Hence, we can study its built-in properties by means of those of the modules M and F. Note that the functional space F where the solutions are sought can be altered and the behaviour of the solutions highly depend on it (in a similar way as for the F-solutions of x 2 + 1 = 0 for F := R or C) [START_REF] Oberst | Multidimensional constant linear systems[END_REF]. In what follows, we shall suppose that F is a rich enough functional space (i.e., is an injective cogenerator left D-module [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]) so that F plays a similar role as the algebraic closure in algebraic geometry. Hence, we can study the properties of ker F (R.) by means of those of M . For the study of the role of F, we refer to [START_REF] Oberst | Multidimensional constant linear systems[END_REF][START_REF] Zerz | An algebraic analysis approach to linear time-varying systems[END_REF] and the references therein.

We also note that hom D (M, F) depends only on the isomorphism type of M , i.e., if M ∼ = M , then we have

hom D (M, F) ∼ = hom D (M , F). If M (resp., M ) is finitely presented by R ∈ D q×p (resp., R ∈ D q ×p ), then we get ker F (R.) ∼ = ker F (R .),
i.e., there is a 1-1 correspondence between the solutions of the first system and the solutions of the second one. For more details and applications of this result to Serre's reduction, Stafford's reduction, the decomposition problem, see [START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF][START_REF] Cluzeau | Equivalences of linear functional systems[END_REF][START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications. Les cours du CIRM[END_REF] and the references therein. Two different representations R ∈ D q×p and R ∈ D q ×p of the same linear system define two isomorphic modules:

M := D 1×p /(D 1×q R) ∼ = M := D 1×p /(D 1×q R ).
Homological algebra methods are developed to study modules up to isomorphism. In particular, even if a particular representation is used, fundamental theorems in homological algebra show that the results do not depend on it. For mathematical systems theory, it is a change of paradigm since systems are usually studied by means of their particular representations (e.g., state-space or polynomial representations). The equivalence between the different approaches is studied below. Within the algebraic analysis approach, we first define the equivalence of linear systems in terms of isomorphic left D-modules finitely presented by these representations, and then use mathematical methods which do only depend on the isomorphism type. For instance, if the concept of controllability is a built-in property of the linear system and not of its representation, then it should be a module property.

For standard classes of linear systems, it has been shown that certain definitions of controllability correspond to the concept of a torsion-free module (for more details, see Section 4.3). Let us introduce basic definitions of module theory [START_REF] Mcconnell | Noncommutative Noetherian Rings[END_REF][START_REF] Lam | Lectures on Modules and Rings[END_REF][START_REF] Rotman | An Introduction to Homological Algebra[END_REF].

Definition 11 Let D be a noetherian domain and M a finitely generated left D-module.

(a) M is free if there exists r ∈ N such that M ∼ = D 1×r . Then, r is called the rank of the free left D-module M and is denoted by rank D (M ).

(b) M is stably free if there exist r, s ∈ N such that M ⊕ D 1×s ∼ = D 1×r . Then, r -s is called the rank of the stably free left D-module M .

(c) M is projective if there exist r ∈ N and a left D-module N such that

M ⊕ N ∼ = D 1×r ,
where ⊕ denotes the direct sum of left D-modules.

(d) M is reflexive if the canonical left D-homomorphism ε : M -→ hom D (hom D (M, D), D), ε(m)(f ) = f (m),
for all f ∈ hom D (M, D) and all m ∈ M , is an isomorphism.

(e) M is torsion-free if the torsion left D-submodule t(M ) of M is 0, where:

t(M ) := {m ∈ M | ∃ d ∈ D \ {0} : d m = 0}.
The elements of t(M ) are called the torsion elements of M .

(f) M is torsion if t(M ) = M , i.e., if every element of M is a torsion element.

Considering s = 0 in (b) (resp., N := D 1×s in (c)) of Definition 11, a free (resp., stably free) module is stably free (resp., projective). A projective module is torsion-free since it can be embedded into a free, and thus into a torsion-free module. The converse of these results are not usually true. In some particular cases, they can hold.

(c) If D is the Weyl algebra A n (k) or B n (k),
where k is a field of characteristic 0, then every finitely generated projective left/right D-module is stably free and every finitely generated stably free left/right D-module of rank at least 2 is free (Stafford's theorem).

(d) If D := A ∂; id A , d
dt where A := k t is the ring of formal power series in t and k is a field of characteristic 0, or A := k{t} is the ring of locally convergent power series in t, where k := R or C, then every finitely generated projective left/right D-module is stably free and every finitely generated stably free left/right D-module of rank at least 2 is free.

In Section 4.3, we will give a dictionary between properties of a linear functional system and properties of the finitely presented left module associated with it.

Basic results of homological algebra

In this section, we briefly review how to effectively check whether or not a finitely presented left D-module M has torsion elements, is torsion-free, reflexive, or projective (see Definition 11), when D is a noetherian domain with finite global dimension [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]. To do that, let us introduce a few concepts of homological algebra [START_REF] Rotman | An Introduction to Homological Algebra[END_REF].

Let D be a noetherian domain, R ∈ D q×p , and

M := D 1×p /(D 1×q R) the left D-module finitely presented by R. If .R ∈ hom D (D 1×q , D 1×p ) is defined by .R : D 1×q -→ D 1×p λ -→ λ R, then we obtain coker D (.R) = D 1×p /im D (.R) = D 1×p /(D 1×q R) = M . Since D is a left noetherian ring, D 1×q
is a noetherian left D-module, i.e., every left D-submodule of D 1×q is finitely generated [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]. In particular, ker D (.R) is a finitely generated left D-module, i.e., there exists a finite generator set {λ i } i=1,...,r of ker D (.R).

Then, we have ker

D (.R) = im D (.R 2 ) = D 1×r R 2 , with the notation R 2 := (λ T 1 . . . λ T r ) T ∈ D r×q .
Let us introduce a few definitions.

Definition 12 (a) A complex of left D-modules is a sequence of left D-modules M i and D-homomorphisms

d i : M i -→ M i-1 for i ∈ Z such that d i • d i+1 = 0, i.e., im d i+1 ⊆ ker d i for all i ∈ Z. Such a complex is denoted by: . . . di+2 / / M i+1 di+1 / / M i di / / M i-1 di-1 / / M i-2 di-2 / / . . . (16)
(b) The defect of exactness of ( 16) at M i is the left D-module H(M i ) := ker d i / im d i+1 .

Inria (c) The complex ( 16) is said to be exact at M i if H(M i ) = 0, i.e., ker d i = im d i+1 , and exact if H(M i ) = 0 for all i ∈ Z.

(d) An exact sequence of the form

. . . .R3 / / D 1×p2 .R2 / / D 1×p1 .R1 / / D 1×p0 π / / M / / 0, (17)
where

R i ∈ D pi×pi-1 and .R i ∈ hom D (D 1×pi , D 1×pi-1 ) is defined by (.R i )λ = λ R i for all λ ∈ D 1×pi , is called a free resolution of M .
With R 1 = R, we can easily check that we have the following exact sequence 17)called a truncated free resolution of Mand using hom D (D 1×pi , F) ∼ = F pi , we then obtain the following complex . . .

0 / / ker D (.R 2 ) i / / D 1×r .R2 / / D 1×q .R1 / / D 1×p π / / M / /
F p2 R3. o o F p1 R2. o o F p0 R1. o o 0, o o (18)
where R i . :

F pi-1 -→ F pi is defined by (R i .) η = R i η for all η ∈ F pi-1 .
The extension Z-modules ext i D (M, F) are then the defects of exactness of [START_REF] Cox | Using Algebraic Geometry[END_REF].

Theorem 5 ([52])

The defects of exactness of [START_REF] Cox | Using Algebraic Geometry[END_REF] depend only on M and F, i.e., they do not depend on the choice of the free resolution [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF] of M . These abelian groups are denoted by:

ext 0 D (M, F) = hom D (M, F) = ker F (R 1 .), ext i D (M, F) = ker F (R i+1 .)/ im F (R i .), i 1.
Theorem 5 is a fundamental result of homological algebra. It shows that the ext i D (M, F)'s do not depend on a particular representation of the linear system.

Remark 8 Let us give an interpretation of the ext i D (M, F)'s. They define the obstructions of the solvability problem which aims at finding η ∈ F pi-1 which satisfies the inhomogeneous linear system R i η = ζ for a fixed ζ ∈ F pi . Indeed, if such an η exists, then we have

R i+1 ζ = R i+1 (R i η) = 0 since ker D (.R i ) = D 1×pi+1 R i+1 , i.e., ζ ∈ ker F (R i+1 .
). This condition is a necessary one for the solvability problem. This problem is solvable if and only if the residue class of The concept of a free resolution of a module can be extended to the concept of a projective resolution in which projective modules are used instead of (finitely generated) free left D-modules D 1×pi [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]. The length of a projective resolution is the number of non-zero projective modules defining this resolution. The minimal length of the projective resolutions of a left D-module M is called the left projective dimension of M and it is denoted by lpd D (M ). The left global dimension of a ring D is the supremum of lpd D (M ) for all left D-modules M and it is denoted by lgld D (M ). For more details, see [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]. Similar definitions can be given for right D-modules. If D is a noetherian ring, i.e., a left and a right noetherian ring, a result due to Kaplansky shows that the left and right global dimensions of D coincide [START_REF] Rotman | An Introduction to Homological Algebra[END_REF] and it is then denoted by gld(D). proved that ext i D (N, D) = 0 for i = 0, . . . , r -1 and ext r D (N, D) = 0, if and only if the algebraic variety defined by all the q × q minors of R has a strict complex dimension equal to n -r. This result is a generalization of the different concepts of coprimeness developed in the literature of multidimensional systems.

ζ in ext i D (M, F) is 0, i.e., if and only if ζ ∈ im F (R i .), which means that η ∈ F pi-1 exists such that ζ = R i η.
According to Theorem 6, certain module properties are characterized by the vanishing of some of the ext i D (N, D)'s. We point out that N := D q /(R D p ) is a right D-module, and thus, it does not define a linear system. To compute the ext i D (N, D)'s, we first have to compute a free resolution of right D-modules

0 N o o D q0 κ o o D q1 Q1. o o D q2 Q2. o o . . . Q3. o o (19)
where Q 1 := R, q 0 := q, and q 1 := p, then dualize it to get the following complex of left D-modules:

0 / / D 1×q0 .Q1 / / D 1×q1 .Q2 / / D 1×q2 .Q3 / / . . . Then, we have ext i D (N, D) = ker D (.Q i+1 )/im D (.Q i ) for i 0,
where we set im D (.Q 0 ) = 0. Since D is a left noetherian ring, the left D-module ker D (.Q i+1 ) is finitely generated, and thus there exists

Q i ∈ D q i-1 ×qi such that ker D (.Q i+1 ) = im D (.Q i ) = D 1×q i-1 Q i , which yields ext i D (N, D) = (D 1×q i-1 Q i )/(D 1×qi-1 Q i ).
If D := O is an Ore algebra satisfying the hypotheses of Theorem 2, then we can use Algorithm 1 to compute the matrix Q i and then Algorithm 2 to check whether or not there exists a matrix

Q i ∈ D q i-1 ×qi-1 such that Q i = Q i Q i , i.e.
, to check whether or not

D 1×q i-1 Q i = D 1×qi-1 Q i , i.
e., whether or not ext i D (N, D) is 0 for i 1. The only point that does not seem to be constructive is the use of Algorithm 1 to compute the free resolution of N since we have to compute right kernel and not left kernel. Moreover, the computation of Gröbner bases is usually not available in computer algebra systems for right ideals or right modules. To do that, we have to use an involution θ of D (see Definition 3). Indeed, it can be used to turn the right D-module structure into a left D-module structure as explained in the next lemma.

Inria Lemma 1 Let N be a right D-module and θ an involution of D. Then, we can define the left D-module N which is equal to N as a set, endowed with the same addition as N , and the left D-action on N is defined by:

∀ d ∈ D, ∀ n ∈ N , d n := n θ(d).
Let M := O 1×p /(O 1×q R) be a left O-module finitely presented by the matrix R ∈ O q×p and let θ be an involution of O. Then, we can define the matrix θ(R) := (θ(R ij )) T ∈ O p×q , i.e., the transpose of the matrix obtained by applying the involution θ to the matrix R component-wise. Note that we always have R = θ 2 (R), i.e., a matrix S can always be written as θ(T ) for a certain matrix T := θ(S). We now consider the left O-module finitely presented by θ(R), namely,

N := O 1×q /(O 1×p θ(R)). (20)
It is called the adjoint module of M . Then, one can prove that (b), (c), (d), and (e) of Theorem 6 hold where N is substituted by N . Hence, we can use Algorithm 1 to compute a free resolution of N

0 N o o D 1×q0 σ o o D 1×q1 .θ(Q1) o o D 1×q2 .θ(Q2) o o . . . , .θ(Q3)
o o then dualizing it by applying the involution θ to get the complex of left D-modules:

0 / / D 1×q0 .Q1 / / D 1×q1 .Q2 / / D 1×q2 .Q3 / / . . . Then, we have ext i D (N, D) = (D 1×q i-1 Q i )/D 1×qi-1 Q i ), where Q i ∈ D q i-1 ×qi is a matrix defined by ker D (.Q i+1 ) = im D (.Q i ).
Finally, as above, using Algorithm 2, we can check whether or not ext i D (N, D) = 0 for i 1. In this way, we can effectively check the conditions of Theorem 6, and thus whether or not M has torsion elements, is torsion-free, reflexive, or projective.

O (N, O) ∼ = ker 0/im O (.R) = O 1×3 /im O (.R) = M . If u (resp., v, p)
denotes the residue class of the first (resp., second, third) element of the standard basis of O 1×3 in M , then eliminating v and p (resp., u and p, resp., u and v) from [START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF] as shown in Example 12, we obtain:

     ∆ (ν ∆ -c) u = 0, ∆ (ν ∆ -c) v = 0, ∆ p = 0.
Hence, each generator u, v, p of M satisfies a PDE, i.e., is a torsion element.

Example 16 Let us illustrate Theorem 6 on a simple linear DTD system defined by:

ẋ1 (t) = x 1 (t) + x 2 (t -1) + u(t), ẋ2 (t) = x 1 (t -1) + x 2 (t) + u(t). (21)
Let

O := Q ∂; id Q , d
dt [δ; σ, 0] be the commutative Ore algebra of DTD operators, where σ is defined by σ(a(t)) = a(t -1) and M := O 1×3 /(O 1×2 R) the O-module finitely presented by the following matrix:

R := ∂ -1 -δ -1 -δ ∂ -1 -1 ∈ O 2×3 .
Let us introduce the Auslander transpose O ( N , O) = 0 for i 3, a fact that we will check again. Using Algorithm 1, we can check that N admits the free resolution

N := O 2 /(R O 3 ) of M . We note that we have N ∼ = N := O 1×2 /(O 1×3 R T ) because O is
0 N o o O 2 σ o o O 3 .R T o o O .Q T 2 o o 0, o o where Q T 2 := (1 1 ∂ -1 -δ).
Dualizing this exact sequence, we get the complex:

0 / / O 1×2 .R / / O 1×3 .Q2 / / O / / 0.
Then, we have

           hom O ( N , O) = ker O (.R) = 0, ext 1 O ( N , O) = ker O (.Q 2 )/im O (.R), ext 2 O ( N , O) = O/(O 1×3 Q 2 ) = 0, ext i O ( N , O) = 0, i 3, since R has full row rank and 1 ∈ O 1×3 Q 2 . Using Algorithm 1 again, we get ker O (.Q 2 ) = O 1×2 R , where: R := 1 -1 0 0 ∂ -1 -δ -1
.

By (b) of Theorem 6, we get t(M ) ∼ = (O 1×2 R )/(O 1×2 R).
It means that the rows of R modulo the system equations define a generating set of the torsion O-submodule t(M ) of M . The first (resp., second) row of R yields the torsion element z

1 := x 1 -x 2 (resp., z 2 := (∂ -1 -δ) x 2 -u = δ z 1 )
. Hence, t(M ) is generated by z 1 . If we consider the following inhomogeneous linear system

     x 1 -x 2 = z 1 , (∂ -1) x 1 -δ x 2 -u = 0, (∂ -1) x 2 -δ x 1 -u = 0,
then computing a Gröbner basis for a monomial order which eliminates x 1 , x 2 , and u, we obtain (∂ +δ -1) z 1 = 0.

Let us now study the torsion-free O-module M/t(M [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF][START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications. Les cours du CIRM[END_REF][START_REF] Robertz | Recent progress in an algebraic analysis approach to linear systems[END_REF]). By (e) of Theorem 6, the O-module M/t(M ) is projective. Using Algorithm 3, we can check that the following matrix

) := O 1×3 /(O 1×2 R ). One can show that M/t(M ) ∼ = O 1×3 Q 2 = O since ext 2 O ( N , O) = 0 (see, e.g.,
L := 1 0 0 0 0 -1
is a left inverse of R T , and thus S := L T is a right inverse of R , which shows again that M/t(M ) is a projective O-module by (f) of Theorem 6. By the Quillen-Suslin theorem (see (b) of Theorem 4), M/t(M ) is then a free O-module of rank 1. This result can be easily checked again by noticing that M/t(M ) is defined by

y 1 -y 2 = 0, (∂ -1 -δ) y 2 -v = 0, ⇐⇒ y 1 = y 2 , v = (∂ -1 -δ) y 2 ,
which shows that y 2 is a basis of M/t(M ). admits a left inverse S. This is equivalent to say that the adjoint system θ(R) λ = 0, i.e., λ + A T λ = 0,

B T λ = 0, (22)
has only the trivial solution λ = 0 since S θ(R) = I q yields λ = S (θ(R) λ) = 0. The above system is not a Gröbner basis for the total degree order since if we differentiate the zero-order equation, we get ḂT λ + B T λ = 0, i.e., using the first-order equation, we obtain the new zero-order equation ( ḂT -B T A T ) λ = 0. We can repeat the same procedure with this last equation. Hence, if we define the sequence of matrices B i defined by B 0 := B T and

B i+1 := Ḃi -B i A T for i 1, we obtain that (B T 0 B T 1 B T 2 .
. .) T λ = 0. Since A is supposed to be noetherian, the increasing sequence of A-submodules O k := k i=0 A 1×m B i of A 1×n stabilizes (see, e.g., [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]), i.e., there exists r ∈ N such that O s = O r for all s r. Then, we get: (see [START_REF] Kailath | Linear Systems[END_REF][START_REF] Kwakernaak | Linear Optimal Control Systems[END_REF]).

(22) ⇐⇒            λ + A T λ = 0,    B 0 . . . B r    λ = 0.

Dictionary between system properties and module properties

Let us introduce a few more definitions. It can be shown that a left D-module F is injective if and only if for every matrix R ∈ D q×p and ζ ∈ F q satisfying R 2 ζ = 0, where R 2 ∈ D r×q is any matrix such that ker D (.R) = im D (.R 2 ), there exists η ∈ F p solving the inhomogeneous linear system R η = ζ [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]. A standard result in homological algebra shows that there always exists an injective cogenerator left module for a ring D [START_REF] Rotman | An Introduction to Homological Algebra[END_REF].

Example 18

If Ω is an open convex subset of R n and F := C ∞ (Ω) or D (Ω) (i.e., the space of distributions with support in Ω), then F is an injective cogenerator

D := k ∂ 1 ; id k , ∂
Example 19 If F is the set of real-valued functions on R which are smooth except for a finite number of points, then F is an injective cogenerator left B 1 (R)-module [START_REF] Zerz | An algebraic analysis approach to linear time-varying systems[END_REF].

Based on the results of [21-23, 45, 46, 56], we can give the following definitions. (c) The linear system is said to be controllable if every observable is free.

(d) The linear system is said to be parametrizable if there exists a matrix Q ∈ D p×m such that ker F (R.) = Q F m , i.e., if for every η ∈ ker F (R.), there exists ξ ∈ F m such that η = Q ξ. Then, Q is called a parametrization and ξ a potential.

(e) The linear system is said to be flat if there exists a parametrization Q ∈ D p×m which admits a left inverse T ∈ D m×p , i.e, T Q = I p . In other words, a flat system is a parametrizable system such that every component ξ i of a potential ξ is an observable of the system. The potential ξ is then called a flat output.

We are now in position to state the correspondence between the properties of a linear system defined in Definition 14 and the properties of the associated finitely generated left module defined in Definition 11.

Theorem 7 ( [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF]) With the hypotheses and notations of Definition 14, we have:

(a) The observables of the linear system are in one-to-one correspondence with the elements of M .

(b) The autonomous elements of the linear system are in one-to-one correspondence with the torsion elements of M . Consequently, the linear system is controllable iff M is torsion-free.

(c) The linear system is parametrizable iff there exists a matrix Q ∈ D p×m such that we have

M := D 1×p /(D 1×q R) ∼ = D 1×p Q,
i.e., iff M is a torsion-free left D-module. Then, the matrix Q is a parametrization, i.e., ker

F (R.) = Q F m .
(d) The linear system is flat iff M is a free left D-module. Then, the bases of M are in one-to-one correspondence with the flat outputs of the linear system.

Example 20 Let us give the system interpretations of the results obtained in Example 16. First, we have the autonomous element z 1 (t) := x 1 (t) -x 2 (t) of ( 21) since it satisfies the autonomous DTD equation ż1 (t)z 1 (t) -z 1 (t -1) = 0. It is a non controllable element of (21) since its trajectory cannot be changed by means of u. Moreover, the controllable system associated with [START_REF] Fliess | Some basic structural properties of generalized linear systems[END_REF] is defined by M/t(M ) := O 1×3 /(O 1×2 R ), which is a free O-module of rank 1. Thus, if F is any O-module (e.g., F := C ∞ (R 0 )), then the corresponding linear system R η = 0, where η := (y 1 y 2 v) T , is flat and y 2 is a flat output. Finally, the matrix Q 2 defined in Example 16 is an injective parametrization of ker F (R .), i.e., we have R η = 0 if and only if there exists ξ ∈ F such that η = Q 2 ξ. Finally, we can check that Q 2 admits a left inverse S 2 := (0 1 0) (see Algorithm 3), which shows that ξ = (S 2 Q 2 ) ξ = S 2 η is uniquely defined by η.

Finally, we illustrate Theorems 6 and 7 with standard linear functional systems coming from control theory and mathematical physics.

Inria

Example 21 (a) Let us consider a wind tunnel model studied in [START_REF] Manitius | Feedback controllers for a wind tunnel model involving a delay: analytical design and numerical simulations[END_REF] and defined by the following linear DTD system: O (N, O) (see [START_REF] Eder | A survey on signature-based Gröbner basis computations[END_REF]) is then a parametrization of ( 23) and we have

     ẋ1 (t) + a x 1 (t) -k a x 2 (t -h) = 0, ẋ2 (t) -x 3 (t) = 0, ẋ3 (t) + ω 2 x 2 (t) + 2 ζ ω x 3 (t) -ω 2 u(t) = 0, (23)
(23) ⇐⇒          x 1 (t) = ω 2 k a z(t -h), x 2 (t) = ω 2 ż(t) -a ω 2 z(t), x 3 (t) = ω 2 z(t) + ω 2 a ż(t), u(t) = z (3) (t) + (2 ζ ω + a) z(t) + (ω 2 + 2 a ω ζ) ż(t) + a ω z(t),
for all z which belongs to an injective module F over the ring of DTD operators.

(b) Let us consider the first group of Maxwell equations defined by

     ∂ B ∂t + ∇ ∧ E = 0, ∇. B = 0, (24)
where B (resp., E) denotes the magnetic (resp., electric) field. We can prove that the differential module associated with [START_REF] Greuel | A Singular introduction to commutative algebra[END_REF] is reflexive (see [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF]). In particular, ( 24) is parametrizable and using the matrix Q 1 obtained in the computation of ext 1 O (N, O) (see ( 19)), we obtain

(24) ⇐⇒    E = - ∂ A ∂t -∇ V, B = ∇ ∧ A, (25)
where ( A, V ) is the so-called quadri-potential formed by smooth functions over R 3 . The second matrix Q 2 defining a free resolution of N (see [START_REF] Eder | A survey on signature-based Gröbner basis computations[END_REF]) then defines a parametrization of the inhomogeneous part of (25), i.e., we have

   - ∂ A ∂t -∇ V = 0, ∇ ∧ A = 0, ⇐⇒    A = ∇ ξ, V = - ∂ξ ∂t ,
where ξ is an arbitrary smooth function on R 3 (used, e.g., for the Lorenz gauge).

(c) Similarly as for the first group of Maxwell equations, we can prove that the equilibrium of the stress tensor defined by

                
defines a reflexive differential module (see [START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications. Les cours du CIRM[END_REF]) and we have the parametrization

(26) ⇐⇒                                                σ x = ∂ 2 χ 3 ∂y 2 + ∂ 2 χ 2 ∂z 2 + ∂ 2 ψ 1 ∂y ∂z , τ yz = - ∂ 2 χ 1 ∂y ∂z - 1 2 ∂ ∂x - ∂ψ 1 ∂x + ∂ψ 2 ∂y + ∂ψ 3 ∂z , σ y = ∂ 2 χ 1 ∂z 2 + ∂ 2 χ 3 ∂x 2 + ∂ 2 ψ 2 ∂z ∂x , τ zx = - ∂ 2 χ 2 ∂z ∂x - 1 2 ∂ ∂y ∂ψ 1 ∂x - ∂ψ 2 ∂y + ∂ψ 3 ∂z , σ z = ∂ 2 χ 2 ∂x 2 + ∂ 2 χ 1 ∂y 2 + ∂ 2 ψ 3 ∂x ∂y , τ xy = - ∂ 2 χ 3 ∂x ∂y - 1 2 ∂ ∂z ∂ψ 1 ∂x + ∂ψ 2 ∂y - ∂ψ 3 ∂z ,
where the ψ i 's and the χ j 's are smooth functions on R 3 . Finally, if we set ψ 1 = ψ 2 = ψ 3 = 0 (resp., χ 1 = χ 2 = χ 3 = 0), then we obtain the so-called Maxwell's parametrization (resp., Morera's parametrization).

For more details, see [START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications. Les cours du CIRM[END_REF].

(d) Let us consider the following time-varying linear OD system:

ẋ1 (t) -t u 1 (t) = 0, ẋ2 (t) -u 2 (t) = 0.
Using Example 17, we can easily check that this system is controllable, i.e., defines a stably free left module over O := A 1 (Q). By (c) of Theorem 4 (i.e., by Stafford's theorem), this module is then free, i.e., the timevarying linear system is flat by (d) of Theorem 7. The effective computation of an injective parametrization is usually a difficult task. To do that, following constructive versions of Stafford's theorems [START_REF] Quadrat | Computation of bases of free modules over the Weyl algebras[END_REF][START_REF] Quadrat | A constructive study of the module structure of rings of partial differential operators[END_REF] and their implementations in the STAFFORD package [START_REF] Quadrat | Computation of bases of free modules over the Weyl algebras[END_REF][START_REF] Quadrat | A constructive study of the module structure of rings of partial differential operators[END_REF], we obtain the following injective parametrization

           x 1 (t) = t 2 ξ 1 (t) -t ξ2 (t) + ξ 2 (t), x 2 (t) = t (t + 1) ξ 1 (t) -(t + 1) ξ2 (t) + ξ 2 (t), u 1 (t) = t ξ1 (t) + 2 ξ 1 (t) -ξ2 (t), u 2 (t) = t (t + 1) ξ1 (t) + (2 t + 1) ξ 1 (t) -(t + 1) ξ2 (t),
where ξ 1 and ξ 2 are arbitrary functions in a left O-module F, and:

ξ 1 (t) = (t + 1) u 1 (t) -u 2 (t), ξ 2 (t) = (t + 1) x 1 (t) -t x 2 (t).
In the language of module theory, {ξ 1 = (t + 1) u 1 -u 2 , ξ 2 = (t + 1) x 1 -t x 2 } is a basis of the free left A 1 (Q)-module M , where x 1 , x 2 , u 1 and u 2 denote here the generators of the module as explained in Section 4.1. Finally, we point out that the above injective parametrization does not contain singularities contrary to

u 1 (t) = t -1 ẋ1 (t), u 2 (t) = ẋ2 (t),
where x 1 and x 2 are arbitrary functions, which admits a singularity at t = 0.

(e) If we consider the following linear DTD system

ẏ1 (t) -y 1 (t -h) + 2 y 1 (t) + 2 y 2 (t) -2 u(t -h) = 0, ẏ1 (t) + ẏ2 (t) -u(t -h) -u(t) = 0, Inria Out[6]= Sn t + Dt t -n -1
The HOLONOMICFUNCTIONS package provides a rather general implementation of Ore algebras, which is advantageous for the applications in control theory discussed in Section 4. For instance, the coefficients of an Ore polynomial ring need not necessarily be polynomials or rational functions. The software also allows us to have, for example, elementary functions in the coefficients:

In[7]:= alg = OreAlgebra[Der[t]] Out[7]= K(t)[Dt; 1, Dt]
In [START_REF] Chyzak | An extension of Zeilberger's fast algorithm to general holonomic functions[END_REF]

:= op = ToOrePolynomial[Cos[t] * * Der[t] * * Sin[t], alg] Out[8]= sin(t) cos(t) Dt + cos 2 (t) In[9]:= op + Sin[t] ˆ2 Out[9]= sin(t) cos(t) Dt + sin 2 (t) + cos 2 (t)
Note that the obvious simplification in the last step is not carried out. By default, HOLONOMICFUNCTIONS keeps the coefficients of Ore polynomials in expanded form, without further simplifications. But there are options to specify a normal form for the coefficients and how to add and multiply them:

In[10]:= alg1 = OreAlgebra[Der[t], CoefficientNormal → Simplify, CoefficientPlus → (Simplify[#1 + #2]&), CoefficientTimes → (Simplify[#1 * #2]&)] Out[10]= K(t)[Dt; 1, Dt]
In [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF]:= op1 = ChangeOreAlgebra[op, alg1]

Out [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF]= sin(t) cos(t) Dt + cos 2 (t)

In [START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF]:= op1 + Sin[t]

ˆ2

Out [START_REF] Cluzeau | Factoring and decomposing a class of linear functional systems[END_REF]= sin(t) cos(t) Dt + 1

Ideally, these options are chosen in a way that expressions identically zero are actually simplified to 0. This is, for instance, not the case when dealing with rational function coefficients in expanded form (as we did above).

Apart from the coefficient domain, HOLONOMICFUNCTIONS provides also a lot of flexibility concerning Ore extensions. As we have seen already, the most common operator symbols are predefined, but there is also a way for the user to define own operator symbols. As an example, we can construct an Ore algebra with a generic Ore extension:

In [START_REF] Cluzeau | OREMORPHISMS: A homological algebraic package for factoring, reducing and decomposing linear functional systems[END_REF] In [START_REF] Mcconnell | Noncommutative Noetherian Rings[END_REF]

:= ToOrePolynomial[d ˆ2 * * t, alg] Out[16]= σ(σ(t))d 2 + (δ(σ(t)) + σ(δ(t))d + δ(δ(t))
Based on the arithmetic of Ore polynomials, an implementation of Buchberger's algorithm for computing Gröbner bases is part of the HOLONOMICFUNCTIONS package. In the following, we consider a family of orthogonal polynomials, namely the Legendre polynomials, which satisfy a second-order differential equation as well as a three-term recurrence. We represent these equations as operators in a suitable Ore algebra and show, by means of a Gröbner basis computation, that Buchberger's product criterion cannot be exploited in noncommutative domains (note that the two Ore polynomials have leading power products D 2 t and S 2 n , whose gcd is 1):

In[17]:= ode = (t ˆ2 -1) * D[f [n, t], t, t] + 2x * D[f [n, t], t] -n(n + 1) * f [n, t] Out[17]= t 2 -1 f (0,2) (n, t) + 2 t f (0,1) (n, t) -n (n + 1)f (n, t)
In [START_REF] Cox | Using Algebraic Geometry[END_REF]

:= rec = (n + 1) * f [n + 1, t] -t(2n + 1) * f [n, t] + n * f [n -1, t] Out[18]= n f (n -1, t) -(2 n + 1) t f (n, t) + (n + 1) f (n + 1, t)
In [START_REF] Eder | A survey on signature-based Gröbner basis computations[END_REF]

:= ops = ToOrePolynomial[{ode, rec}, f [n, x]] Out[19]= {(t 2 -1) D 2 t + 2 t Dt + (-n 2 -n), (n + 2) S 2 n + (-2 n t -3 t) Sn + (n + 1)} In[20]:= OreGroebnerBasis[ops] Out[20]= {(-n -1) Sn + (t 2 -1) Dt + (nt + t), (t 2 -1) D 2 t + 2 t Dt + (-n 2 -n)}
Although this paper is mostly about applications of the above-described methods in control theory, we want Inria to mention briefly the main application for which the HOLONOMICFUNCTIONS package has been developed. That is: proving special function identities, involving integrals and symbolic sums, in the spirit of Zeilberger's holonomic systems approach [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF]. Once the input functions are represented by their annihilators (together with initial conditions), one can use Gröbner basis techniques to compute the annihilator of an integral or sum, by employing the method of creative telescoping [START_REF] Zeilberger | The method of creative telescoping[END_REF]. An identity then is established, for example, by observing that both sides satisfy the same differential equation or recurrence. As an example, consider the following identity involving the Laguerre polynomials L a n (t) and the Bessel function J a (t):

e -t t a/2 n! L a n (t) = +∞ 0 e -τ τ a 2 +n J a 2 √ τ t dτ. (27)
By using closure properties of holonomic functions, the HOLONOMICFUNCTIONS package automatically computes the annihilator of the function on the left-hand side of [START_REF] Kailath | Linear Systems[END_REF]. The result is given as a Gröbner basis: For the right-hand side of ( 27), one computes the annihilator of the integrand, and then applies creative telescoping to it, in the form of Chyzak's algorithm [START_REF] Chyzak | An extension of Zeilberger's fast algorithm to general holonomic functions[END_REF]: Note that the first part of Out [START_REF] Fliess | Controllability and observability of linear delay systems: an algebraic approach[END_REF] agrees (up to sign) with Out [START_REF] Fliess | Some basic structural properties of generalized linear systems[END_REF], the annihilator of the left-hand side. In order to complete the proof of [START_REF] Kailath | Linear Systems[END_REF], one has to investigate whether the certificate (the second part of Out [START_REF] Fliess | Controllability and observability of linear delay systems: an algebraic approach[END_REF]) contributes an inhomogeneous part to the computed equations (this is not the case here), and one has to compare initial values. These steps are currently beyond the capabilities of the package and have to be done by hand; see the examples in [START_REF] Koutschan | Advanced applications of the holonomic systems approach[END_REF] where this is demonstrated in detail.

In[21]:= Annihilator[Exp[-t] * t ˆ(a/2) * n! * LaguerreL[n, a, t], {S[a], S[n], Der[t]}] Out[21]= {2 Sn -2 t Dt + (-a -2 n -2),
In[22]:= ann = Annihilator[Exp[-τ ] * τ ˆ(a/2 + n) * BesselJ[a, 2 Sqrt[τ t]], {S[a], S[n], Der[t], Der[τ ]}] Out[22]= {2 t Dt -2 τ Dτ + (a + 2 n -2 τ ), Sn -τ, τ 2 D 2 τ + (-a τ -2 n τ + 2 τ 2 + τ ) Dτ + (a n -a τ + n 2 -2 n τ + τ 2 + τ t + τ ), t S 2 a + (a τ + τ ) Dτ + (-a 2 -a n + a τ -a -n + τ t + τ )} In[23]:= CreativeTelescoping[ann, Der[τ ]] Out[23]= {{-2 Sn + 2 t Dt + (a + 2 n + 2),

The OREALGEBRAICANALYSIS package

A Mathematica package, called OREALGEBRAICANALYSIS, has been recently developed by the first, third, and fourth-named authors 1 . It is freely available with a library of examples (see [START_REF] Cluzeau | OREALGEBRAICANALYSIS: A Mathematica package for the algorithmic study of linear functional systems[END_REF]).

The OREALGEBRAICANALYSIS package can be used to study (determined/over-determined/underdetermined) linear functional systems appearing, e.g., in control theory and in mathematical physics. For instance, structural properties of linear functional systems can algorithmically be decided (e.g., existence and computation of autonomous elements, (injective, minimal, chain of) parametrizations, potentials, flat outputs, decide Willems' controllability and observability). We point out that the algorithms implemented in this package are generic in the sense that they do not depend on the Ore algebras.

To define, manipulate, and compute in Ore algebras of functional operators, we use the Mathematica package HOLONOMICFUNCTIONS described in the previous section. The package OREALGEBRAICANALYSIS extends these Gröbner basis techniques to finitely presented left modules over the same classes of Ore algebras. It also contains algorithms for module theory (e.g., test whether or not a module admits torsion elements, is torsion-free, reflexive, projective, stably free, free) and homological algebra (e.g., computation of free resolutions, projective dimension, extension modules with value in the underlying ring, invariants, . . . ).

The OREALGEBRAICANALYSIS package includes the main procedures implemented in the Maple packages OREMODULES [START_REF] Chyzak | OREMODULES: A symbolic package for the study of multidimensional linear systems[END_REF] and OREMORPHISMS [START_REF] Cluzeau | OREMORPHISMS: A homological algebraic package for factoring, reducing and decomposing linear functional systems[END_REF]. Since HOLONOMICFUNCTIONS can handle larger classes of Ore algebras than the Maple package ORE -ALGEBRA 2 , OREALGEBRAICANALYSIS can study larger classes of linear functional systems than the Maple packages OREMODULES and OREMORPHISMS. Moreover, the internal design of Mathematica can allow us to consider classes of systems which could not easily be considered in Maple such as generic linearizations of nonlinear functional systems defined by explicit equations and systems containing transcendental functions (e.g., trigonometric functions, special functions). See the following examples.

We will now shortly illustrate the main functions and applications of the OREALGEBRAICANALYSIS package with explicit examples. For more examples, see [START_REF] Cluzeau | OREALGEBRAICANALYSIS: A Mathematica package for the algorithmic study of linear functional systems[END_REF].

Example 22 Let us consider an example studied in [START_REF] Marquez-Martinez | A note on the accessibility for nonlinear time-delay systems[END_REF]. We start the Mathematica session by loading the package In [START_REF] Greuel | A Singular introduction to commutative algebra[END_REF]:= << OreAlgebraicAnalysisà nd then entering the system equations in the form:

In[25]:= eqs = {x 1 [t] → x1[t]u[t] + u[t -2],
x The matrix R of DTD operators which defines the generic linearization of the above nonlinear system is then given by:

In

Out[30]=      0 -D t 0 S -1 t + 1 -D t + u[t] 0 -D t S -1 t + x 1 [t] 0 0 D t S -1 t 2 -S -1 t 0 -S -1 t 2 + S -1 t -S -1 t -1 0     
The matrix Q is a parametrization of the controllable part (it is too large to be printed here; see [START_REF] Cluzeau | OREALGEBRAICANALYSIS: A Mathematica package for the algorithmic study of linear functional systems[END_REF]). Since Ann is not the identity matrix, we deduce that M admits nontrivial torsion elements and thus the corresponding system admits autonomous elements τ 1 , . . . , τ 4 , defined by:

In [START_REF] Kredel | Solvable Polynomial Rings[END_REF] We note that the first three autonomous elements τ 1 , τ 2 , τ 3 are trivial. The only nontrivial autonomous element is τ 4 = -dx 2 (t -2) + dx 2 (t -1) -dx 3 (t -1) -dx 3 (t), which satisfies τ4 = 0.

Example 23 Let us consider the following nonlinear DTD system considered in [START_REF] Califano | Controllability of driftless nonlinear time-delay systems[END_REF]:

In [START_REF] Lam | Lectures on Modules and Rings[END_REF] the matrix R of DTD operators which defines the generic linearization of the above nonlinear system

Out[39]=   0 -1 x 3 [t] 0 0 0 -D t 1 D t 0 -u[t] x 3 [t -1] S -1 t -x 3 [t -1]  
The matrix Q is too large to be printed here. For more details, see [START_REF] Cluzeau | OREALGEBRAICANALYSIS: A Mathematica package for the algorithmic study of linear functional systems[END_REF].

In [START_REF] Malgrange | Systèmes différentiels à coefficients constants[END_REF] The A-linear relations among the autonomous elements are given by:

In [START_REF] Marquez-Martinez | A note on the accessibility for nonlinear time-delay systems[END_REF]:= rels Out [START_REF] Marquez-Martinez | A note on the accessibility for nonlinear time-delay systems[END_REF]

= {u[t]τ [1][t -1] + τ [3][t] == 0, -x[3][t]τ [2][t] -τ [1] [t] == 0, τ [2][t] == 0}
Let us now prove that the set of autonomous elements can be generated by τ 1 . Let us introduce the matrix L defining rels: In [START_REF] Quadrat | An introduction to constructive algebraic analysis and its applications. Les cours du CIRM[END_REF]

:= Thread[Table[τ [i][t], {i, 3}] → ApplyMatrix[V, {τ [1][t]}]] Out[47]= {τ [1][t] → τ [1][t], τ [2][t] → 0, τ [3][t] → -u[t]τ [1][t -1]}
From this point, we will use some procedures which are not freely available (see [START_REF] Belikov | NLControl: A software project addressing Non-Linear Control systems[END_REF]). Finally, let us integrate the one-form defined by τ 1 :

In [START_REF] Quadrat | Computation of bases of free modules over the Weyl algebras[END_REF] Out [START_REF] Quadrat | A constructive study of the module structure of rings of partial differential operators[END_REF]= {x2[t] -

1 2 x3[t] 2 }
Thus, x 1 is an autonomous element of the nonlinear DTD system.

Example 1

 1 Let us give a few examples of skew polynomial rings.

Remark 1

 1 In Definition 2, the numbers m and n can be different. For instance, considering again (d) of Example 1, i.e., the Ore algebra O := A[∂;id A , δ][S; σ, 0], where, for instance, A := k[t], then we have m = 2 and n = 1.

  Let O := k(n)[S;σ,0] be the skew polynomial ring of forward shift operators considered in (b) of Example 1. Then, an involution of O can be defined by θ(n) := -n and θ(S) := -S. (d) Let O := k[t][∂;id,δ][S;σ,0] be the Ore algebra of differential time-delay operators defined by δ := d dt , and σ(a(t)) := a(t -1), where a ∈ k[t]. Then, an involution of O can be defined by θ(t) := -t, θ(∂) := ∂, and θ(S) := S.

Example 4

 4 In fluid mechanics, Stokes equations, which describe the flow of a viscous and incompressible fluid at low Reynolds number, are defined by

  The operators ∂ x and ∂ y commute, i.e., ∂ x ∂ y = ∂ y ∂ x , because of Schwarz's theorem and (12) has only constant coefficients. An element P ∈ D is of the form P = 0 µx+µy r a µ ∂ µx x ∂ µy y

  called a set of generators of M . Similar definitions hold for right D-modules. In what follows, we consider D to be a polynomial Ore algebra O. Let Mon(O) be the set of monomials of O and {f j } j=1,...,p the standard basis of the free finitely generated left O-module O 1×p := {(λ 1 . . . λ p ) | λ i ∈ O, i = 1, . . . , p}, namely the k th component of f j is 1 if k = j and 0 otherwise. First, we extend the monomial order ≺ from Mon(O) to the set of monomials of the form u f j , where u ∈ Mon(O) and j = 1, . . . , p, i.e., to Mon(O 1×p ) := p j=1 Mon(O) f j . This extension is also denoted by ≺ and it has to satisfy the following two conditions: (a) ∀ w ∈ Mon(O) :

Definition 10

 10 Let ≺ be an admissible monomial order on Mon(O), u, v ∈ Mon(O), and {f j } j=1,...,p the standard basis of the left O-module O 1×p . (a) The term over position order on Mon(O 1×p ) induced by ≺ is defined by:

Example 12

 12 We consider again Example 4. Let us compute a Gröbner basis of the O

Algorithm 1

 1 Computation of the left kernel of R ∈ O q×p , i.e., find S ∈ O r×q such that: ker O (.R) := {λ ∈ O 1×q | λ R = 0} = O 1×r S := {µ S | µ ∈ O 1×r }. • Input: An Ore algebra O satisfying the hypotheses of Theorem 2 and a finitely generated left O-submodule L := O 1×q R of O 1×p , where R ∈ O q×p . • Output: A matrix S ∈ O r×q such that ker O (.R) = O 1×r S. (a) Introduce the indeterminates η 1 , . . . , η p , ζ 1 , . . . , ζ q over O and define the set:

Algorithm 2

 2 r} by selecting the elements of G containing only the ζ i 's, and return S := (S ij ) ∈ O r×q . Computation of a left factorization: given two matrices R ∈ O q×p and R ∈ O q ×p , find a matrix R ∈ O q×q (if it exists) satisfying R = R R . • Input: An Ore algebra O satisfying the hypotheses of Theorem 2 and two matrices R ∈ O q×p and R ∈ O q ×p . • Output: R ∈ O q×q (if it exists) such that R = R R and [ ] otherwise. (a) Introduce the indeterminates η 1 , . . . , η p , ζ 1 , . . . , ζ q over O and define the set:

Algorithm 3

 3 Computation of a left inverse: given a matrix R ∈ O q×p , find (if it exists) a left inverse S ∈ O p×q of R over O, namely S R = I p . • Input: An Ore algebra O satisfying the hypotheses of Theorem 2 and R ∈ O q×p . • Output: A matrix S ∈ O p×q such that S R = I p if S exists and [ ] otherwise. (a) Introduce the indeterminates η 1 , . . . , η p , ζ 1 , . . . , ζ q over O and define the set:

  e., computation of right kernels, right factorizations, and right inverses) can be obtained by considering an involution of the Ore algebra O (see Definition 3). For instance, the computation of a right inverse of a matrix R ∈ O q×p over an Ore algebra can be done by applying Algorithm 3 to the matrix θ(R) := (θ(R ij )) T ∈ O p×q (obtained by applying an involution θ of O to each entry R ij of R and then transposing the result) and applying the involution to the left inverse T ∈ O q×p of θ(R) to get S := θ(T ) ∈ O p×q which then satisfies: R S = θ 2 (R) θ(T ) = θ(T θ(R)) = θ(I q ) = I q .

Remark 6

 6 If D := O is an Ore algebra satisfying the hypotheses of Theorem 2, then we can check if π(λ) = π(λ ) for λ, λ ∈ O 1×p since λ -λ ∈ L := O 1×q R if and only if red(λ -λ , G) = 0, where G is a Gröbner basis of L (see Section 3).

Remark 7 2 2 σ 2

 722 If D is not a commutative ring, then neither ker F (R.) nor hom D (M, F) are left D-modules. For instance, if we consider D := A 1 (Q), R := ∂ + t-m σ 2 ∈ D where t, m and σ are constants parameters (e.g., transcendental elements over Q), and M := D/(D R), then η := e -(t-m)

Theorem 4 (

 4 [START_REF] Mcconnell | Noncommutative Noetherian Rings[END_REF][START_REF] Lam | Lectures on Modules and Rings[END_REF][START_REF] Quadrat | A constructive study of the module structure of rings of partial differential operators[END_REF][START_REF] Rotman | An Introduction to Homological Algebra[END_REF]) We have the following results.(a) If D is a principal ideal domain, i.e., every left/right ideal of the domain D is principal (e.g., D :=A ∂; id A , d dt , where A := k, k(t), or k t [t -1] is the field of formal Laurent power series, where k is a field of characteristic 0, or A := k{t}[t -1 ] is the field of Laurent power series, where k := R, C), then every finitely generated torsion-free left/right D-module is free.(b) If D := k[x 1 , . . . ,x n ] is a commutative polynomial ring with coefficients in a field k, then every finitely generated projective D-module is free (Quillen-Suslin theorem).

Remark 9

 9 If F := D, then the ext i D (M, D)'s inherit a right D-module structure.

Example 14 (

 14 [START_REF] Mcconnell | Noncommutative Noetherian Rings[END_REF]) We have the following examples.(a) If A has finite left global dimension and σ is an automorphism of A, then we have:lgld(A) lgld(A[∂; σ, δ]) lgld(A) + 1. Moreover, if δ = 0, then we have lgld(A[∂; σ, δ]) = lgld(A) + 1. (b) If k is a field, then gld(k[x 1 , . . . , x n ]) = n. (c) Ifk is a field of characteristic 0 (e.g., k := Q, R, C), then gld(A n (k)) = n and gld(B n (k)) = n. Theorem 6 ([9]) Let D be a noetherian ring with finite global dimension gld(D) := n, M := D 1×p /(D 1×q R) the left D-module finitely presented by the matrix R ∈ D q×p , and N := D q /(R D p ) the so-called Auslander transpose of M . Then, we have the following results: (a) M is a torsion left D-module if and only if hom D (M, D) = 0. (b) t(M ) ∼ = ext 1 D (N, D). (c) M is a torsion-free left D-module if and only if ext 1 D (N, D) = 0. (d) M is a reflexive left D-module if and only if ext i D (N, D) = 0 for i = 1, 2. (e) M is a projective left D-module if and only if ext i D (N, D) = 0 for i = 1, . . . , n. (f) If R is a full row rank matrix, i.e., ker D (.R) = 0, then M is a projective left D-module if and only if N ∼ = ext 1 D (M, D) = 0, i.e., if and only if R admits a right inverse. Remark 10 If D := k[x 1 , . . . , x n ], then Theorem 6 and (a) of Example 14 show that the concepts of torsionfree, reflexive, and projective modules are instances of a sequence of n module properties characterized by the successive vanishing of the ext i D (N, D)'s for i = 1, . . . , n. If R has full row rank and k := R or C, it can be

Example 15

 15 Let us consider again Stokes equations defined in Example 4. With the notations of Example 12, using Theorem 6, we can easily prove that the finitely presented O-module M := O 1×3 /L is torsion. Indeed, since det R = 0, R has full row rank, i.e., ker O (R.) = 0, and we have the following free resolution of N duality, yields the following complex 0 / / O 1×3 .R / / O 1×3 / / 0, and thus we get t(M ) ∼ = ext 1

  a commutative ring and θ = id O . Let us explicitly compute the ext i O ( N , O)'s. Since gld(O) = 2 (see (b) of Example 14), one can prove that ext i

  Finally, since O is a commutative polynomial ring, we can use Remark 10 to prove again the results obtained above. Indeed, the ideal Fitt 0 (N ) defined by all the 2 × 2 minors of R is defined by the ideal I := (∂ + δ -1). The algebraic variety formed by the zeros of I is ∂ + δ -1 = 0, which is 1-dimensional. Using Remark 10, we then get ext 1 O (N, O) = 0, which proves again that t(M ) = 0. Similarly, ifN := O 2×1 /(R O 3×1 ) is the Auslander transpose of M/t(M ), then we have Fitt 0 (N ) = (∂ -1 -δ, 1) = O, which shows that ext i O (N , O) = 0 for i = 1,2, which proves again that M/t(M ) is a projective and thus a free O-module. Inria Finally, if R is a full row rank matrix, then (f) of Theorem 6 shows that M is a projective left D-module if and only if we have N ∼ = ext 1 D (M, D) = 0 or equivalently if and only if N = 0. In this case, we do not have to test the vanishing of all the ext i D (N, D)'s for i = 1, . . . , n as shown in (e) of Theorem 6. Example 17 Let O := A ∂; id A , d dt be a ring of OD operators with coefficients in a noetherian differential ring A, A ∈ A n×n , B ∈ A n×m , and the left O-module M := O 1×(n+m) /(O 1×n R) finitely presented by R := (∂ I n -A -B) ∈ O n×(n+m) which defines the linear system ẋ(t) = A x(t) + B u(t). Since R has full row rank, (f) of Theorem 6 shows that M is a projective left O-module if and only if R admits a right inverse, i.e., using the involution θ defined in (b) of Example 3, if and only if θ(R) := (-∂ I n -A T -B T ) T ∈ O (n+m)×n

Hence, the above 0 B T 1 B T 2 .

 012 system has the only solution λ = 0 if and only if the matrix O := (B T . . B T r ) T admits a left inverse with entries in A or equivalently if and only if O T admits a right inverse with entries in A. If A is a field, then we can take r = n -1 by the Cayley-Hamilton theorem and the condition on the existence of a right inverse for O T then becomes that O T has full row rank. Finally, if A := k is a field of constants, i.e., ȧ = 0 for all a ∈ A, as, e.g., A = Q or R, then we get the standard controllability condition rank k (B A B . . . A n-1 B) = n

Definition 13 (

 13 [START_REF] Rotman | An Introduction to Homological Algebra[END_REF]) We have the following definitions: (a) A left D-module F is said to be injective if for every left D-module M , we have ext i D (M, F) = 0 for i 1.(b) A left D-module F is said to be cogenerator if for every left D-module M and m ∈ M \ {0}, there exists ϕ ∈ hom D (M, F) such that ϕ(m) = 0.

Definition 14 (

 14 [START_REF] Chyzak | Effective algorithms for parametrizing linear control systems over Ore algebras[END_REF]) Let D be a noetherian domain, R ∈ D q×p , F an injective cogenerator left D-module, and the linear system (behaviour) defined by R and F:ker F (R.) := {η := (η 1 . . . η p ) T ∈ F p | R η = 0}.(a) An observable is a D-linear combination of the system variables η i .(b) An observable ψ(η) is called autonomous if it satisfies a D-linear relation by itself, i.e., d ψ(η) = 0 for some d ∈ D \ {0}. An observable is said to be free if it is not autonomous.

1 O

 1 where a, k, ω, and ζ are constant parameters. Checking that ext (N, O) = 0, (c) of Theorem 6 shows that (23) defines a torsion-free module over the ring of DTD operators, and thus (23) is parametrizable by (c) of Theorem 7. The matrix Q 1 obtained during the computation of ext 1

  := OreSigma[d] := σ; In[14]:= OreDelta[d] := δ; In[15]:= alg = OreAlgebra[d] Out[15]= K[d; σ, δ]

1 t 2 - 1 t- 1 S - 1 t 2 -

 121112 [START_REF] Kailath | Linear Systems[END_REF]:= MatrixForm[R = ToOrePolynomialD[ eqs, vars, A ]]Let M = A 1×4 /(A 1×3 R) be the left A-module finitely presented by the matrix R. The adjoint of R is then defined by:In[28]:= MatrixForm[ Radj = Involution[ R, A ]] x 1 [-t] -S -Let us check whether or not M is a torsion-free left A-module:In[START_REF] Koutschan | Advanced applications of the holonomic systems approach[END_REF]:= {Ann, Rp, Q} = Exti[ Radj, A, 1 ];MatrixForm[ Ann ]

In[ 36 ]

 36 := MatrixForm[R = ToOrePolynomialD[eqs, vars, A]] A-module M = A 1×4 /(A 1×3 R) finitely presented by R. Let us first compute the adjoint of R:In[37]:= MatrixForm[Radj = Involution[R, A]] -t] D t -x 2 [-1 -t] -x 3 [-t]Let us check whether or not M is a torsion-free left A-module:In[38]:= {Ann, Rp, Q} = Simplify[ Exti[ Radj, A, 1 ]]; MatrixForm[ Ann ] [t]D t -u [t]   In[39]:= MatrixForm[ Rp ]

In[ 43 ]

 43 := MatrixForm[L = ToOrePolynomialD[ rels, {τ [1][t], τ [2][t], τ [3][t]}, A]]Let us consider the following matrixIn[START_REF] Ore | Theory of non-commutative polynomials[END_REF]:= MatrixForm[γ = {{1, 0, 0}} ] Out[44]=1 0 0 which corresponds to the position of τ 1 . To express τ 2 and τ 3 in terms of τ 1 , we first check whether or not the matrix T , formed by stacking L with γ, admits a left inverse.In[START_REF] Pommaret | Partial Differential Control Theory[END_REF]:= U = LeftInverse[T = Join[L, γ], A] Out[45]= {{0, 0, 0, 1}, {0, 0, -1, 0}, {1, 0, 0, -u[t] S -1 t }}Hence, if we consider the last column of the left inverse U of T , i.e.In[START_REF] Pommaret | Algebraic analysis of linear multidimensional control systems[END_REF]:= MatrixForm[V = Take[U, All, -1]]

  := BookForm[sp = SpanK[{ApplyMatrixD[Rp[[1]], vars]}, t]] Out[48]= SpanK[-dx2[t] + x3[t]dx3[t]]In[START_REF] Quadrat | A constructive study of the module structure of rings of partial differential operators[END_REF]:= IntegrateOneForms[sp]

  is the first derivative of the function y. Note that d i dt i is the i th composition of the operator d dt . If the composition of operators is simply denoted by the standard product, we haved i dt i = dIt is important to note that the element a i ∈ R in the expression of P is seen as the multiplication operator y -→ a i y, and a i ∂ i stands for the composition of the two operators a i and ∂ i . As understood by G. Boole, the set of OD operators forms the commutative polynomial ring R[∂]. Algebraic techniques (e.g., Euclidean division) can then be used to study linear OD equations with constant coefficients. More generally, if A is a differential ring, namely a ring equipped with a derivation d

	dt	i . Hence, P
	can be rewritten as the polynomial P =

r i=0 a i ∂ i in ∂ := d dt with coefficients in R. dt : A -→ A satisfying the additivity condition and Leibniz's rule, namely,

  Similarly as the previous case, if h ∈ R 0 and A is a difference ring of functions of t with σ(a(t)) = a(t-h) for all a ∈ A as an endomorphism, then we can define the ring O := A[S;σ,0] of TD operators in S with coefficients in A. , δ] defined in (a) and then define the Ore extension O := B[S;σ,0] of B, where σ is the endomorphism of B defined by σ(a(t)) = a(t -h) for all a ∈ A and σ

	We have S a(t) = a(t -h) S, which is exactly (6).
	(d) If we want to reformulate (2) within the language of Ore extensions, we have to define the ring of DTD
	operators. To do that, we can first consider a difference-differential ring (A, σ, δ) and the skew polynomial ring B := A[∂;id A

  If A is a domain, i.e., A does not have non-trivial zero divisors, and σ is an injective endomorphism of A, If A is a left Ore domain, i.e., a domain A which satisfies the left Ore property which states that for a 1 , a 2 ∈ A\{0}, there exist b 1 , b 2 ∈ A\{0} such that b 1 a 1 = b 2 a 2 , and α is injective, then D is a left Ore domain. If A is a left (resp., right) noetherian ring, i.e., every left (resp., right) ideal of A is finitely generated, and α is an automorphism of A, then D is a left (right) noetherian ring. Moreover, if A is a domain, then D is a left Ore domain.

	(c) As shown in (d) of Example 1, we can iterate the construction of an Ore extension to obtain a multivariate
	noncommutative polynomial ring:
	Theorem 1 ([16]) Let D := A[∂;σ,δ] be a left skew polynomial ring over a ring A. Then, we have:
	(a) then D is a domain.
	(b) Inria

  n, is called the polynomial (resp., rational) Weyl algebra of PD operators with coefficients in A. It is denoted by A n (k) (resp., B n (k)).

	(b) We can combine the two skew polynomial algebras defined in (a) and (b) of Example 1 to obtain the Ore algebra O := Q(n,t)[∂;id,δ][S;σ,0] of differential-shift operators with coefficients in Q(n,t). The Bessel
	function of the first kind J n (t) satisfies the following functional equation:

  0, where i is the canonical injection and π the canonical projection. Repeating for R 2 what we did for R and so on, we get a free resolution (17) of M . If D := O is an Ore algebra satisfying the hypotheses of Theorem 2, then Algorithm 1 can be used to compute a free resolution of M .Applying the left exact contravariant functor hom

D (•, F)

[START_REF] Rotman | An Introduction to Homological Algebra[END_REF]

to the complex

. . . .R3 / / D 1×p2 .R2 / / D 1×p1 .R1 / / D 1×p0 / / 0,

obtained by removing M from (

  4 t 2 D 2 t + (4 t 2 + 4 t) Dt + (-a 2 + 2 a t + 4 n t + 4 t), 2 t S 2

a + (2 a t + 2 t 2 + 2 t) Dt + (-a 2 + a t -a + 2 n t + 2 t)}

  4 t 2 D 2

t + (4 t 2 + 4 t) Dt + (-a 2 + 2 a t + 4 n t + 4 t), 2 t S 2 a + (2 a t + 2 t 2 + 2 t) Dt + (-a 2 + a t -a + 2 n t + 2 t)}, {-2 τ, -4 τ t, -2 τ t}}

  2 [t] → u[t] + u[t -1], x 3 [t] → u[t -1] -u[t -2]}; vars = { x1[t], x2[t], x3[t], u[t] };Let us now introduce the following Ore algebra A of DTD operators:In[START_REF] Kandri-Rody | Non-commutative Gröbner bases in algebras of solvable type[END_REF]:= replA = ModelToReplacementRules[ eqs, t ]; A = OreAlgebraWithRelations[ Der[t], S[-1][t], replA ]

Out

[START_REF] Kandri-Rody | Non-commutative Gröbner bases in algebras of solvable type[END_REF]

= K(t)[Dt; 1, Dt][ S -1 t ; #1/ . t → t -1&, 0&]

  [START_REF] Bergman | The Diamond lemma for ring theory[END_REF][t] → du[-2 + t] -du[-1 + t] + dx 3 [t], τ [4][t] → -dx2[-2 + t] + dx2[-1 + t] -dx3[-1 + t] -dx3[t] } In[32]:= eqs Out[32]= { τ [1][t] == 0, τ [2][t] == 0, τ [3][t] == 0, τ [4] [t] == 0 } In[33]:= rels Out[33]= {-τ [2][t] -τ [3][t] == 0, -τ [1][t] == 0, τ [3][t] == 0, -τ [1][-2 + t] + τ [1][-1 + t] + τ [3][-1 + t] + τ [3][t] + τ [4] [t] == 0 }

:= {aut, eqs, rels} = AutonomousElements[ R, {dx1[t], dx2[t], dx3[t], du[t]}, τ, A, Relations → True]; aut Out[31]= {τ [1][t] → du[-1 + t] + du[t] -dx 2 [t], τ [2][t] → du[-1 + t] + u[t]dx1[t] + du[t]x1[t] -dx 1 [t] -dx 3 [t],

Inria τ

  := eqs = {x 1 [t] → x2[t -1]u[t], x 2 [t] → x3[t]u[t], x 3 [t] → u[t]}; vars = {x1[t], x2[t], x3[t], u[t]};Let us introduce the following Ore algebra A of DTD operatorsIn[START_REF] Lange-Hegermann | Thomas Decomposition and Nonlinear Control Systems[END_REF]:= replA = ModelToReplacementRules[ eqs, t ]; A = OreAlgebraWithRelations[ Der[t], S[-1][t], replA ]

Out

[START_REF] Lange-Hegermann | Thomas Decomposition and Nonlinear Control Systems[END_REF]

= K(t)[Dt; 1, Dt][ S -1 t ; #1/ . t → t -1&, 0&]

  := {aut, eqs, rels} = AutonomousElements[ R, {dx1[t], dx2[t], dx3[t], du[t]}, τ, A, Relations → True]; aut Out[40]= {τ [1][t] → -dx2[t] + dx3[t]x3[t], τ [2][t] → du[t] -dx 3 [t], Out[41]= {τ [1] [t] == 0, τ [2][t] == 0, -τ [3][t]u [t] + u[t]τ [3] [t] == 0}

τ [3][t] → -du[t]x2[t -1] -u[t]dx3[t -1]x3[t -1] + dx 1 [t]}

The autonomous elements τ 1 , τ 2 , τ 3 satisfy the following equations:

In

[START_REF] Manitius | Feedback controllers for a wind tunnel model involving a delay: analytical design and numerical simulations[END_REF]

:= eqs

RR n°8999

∂x1 • • • ∂ n ; id k , ∂ ∂xn -module,where k := R or C. For more details, see[START_REF] Oberst | Multidimensional constant linear systems[END_REF] and the references therein.RR n°8999

This work was supported by the PHC PARROT

29586NG between France and Estonia. 2 http://algo.inria.fr/chyzak/Mgfun/Sessions/Ore_algebra.html.RR n°8999

and M. Tõnso were supported by the PHC Parrot CASCAC (29586NG).

C. Koutschan was supported by the Austrian Science Fund (FWF): W1214.

using (e) or (f) of Theorem 6, then we can check that it defines a projective module over the commutative polynomial ring of DTD operators with constant coefficients. By the Quillen-Suslin theorem (see (b) of Theorem 4), this module is free. The computation of bases and injective parametrizations is usually intricate and requires an effective version of the Quillen-Suslin theorem [START_REF] Fabiańska | Applications of the Quillen-Suslin theorem to multidimensional systems theory[END_REF]. Using the QUILLENSUSLIN package [START_REF] Fabiańska | Applications of the Quillen-Suslin theorem to multidimensional systems theory[END_REF], we get the following injective parametrization

for all ξ belonging to a module over the ring of DTD operators. Finally, note that y 1 defines a basis of the free module defined by the above system.

5 Mathematica packages

The HOLONOMICFUNCTIONS package

The Mathematica package named HOLONOMICFUNCTIONS has been developed by the second-named author in the frame of his Ph.D. thesis [START_REF] Koutschan | Advanced applications of the holonomic systems approach[END_REF]. It can be downloaded for free from the website http://www.risc.jku. at/research/combinat/software/HolonomicFunctions/, and a complete documentation is given in the manual [START_REF] Koutschan | HolonomicFunctions (user's guide)[END_REF]. We start with a fresh Mathematica session and load the package with the following command:

In In its core, the package provides functionality to construct Ore algebras and to work with Ore polynomials. First, we demonstrate how this is done using a very standard application, namely the operator from (b) of Example 2. For this purpose, we define a multivariate Ore algebra with rational function coefficients, which is built up of the shift operator S n and the ordinary differential operator D t :

The symbol K in Out [START_REF] Belikov | NLControl: A software project addressing Non-Linear Control systems[END_REF] has no particular meaning, and just indicates that the constant field can be everything that covers the user's input; for example K could contain the rational numbers Q as a proper subfield.

We can now convert an input expression to an Ore polynomial that belongs to this Ore algebra and do some arithmetic (note the usage of the noncommutative multiplication * * in Mathematica):

To construct an Ore algebra with polynomial coefficients, we just have to include the variables n and t in the command. Note that each monomial is displayed according to the order in which the generators of the ring are given:

In [START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenrings nach einem nulldimensionalen Polynomideal[END_REF]