Document [original]

Technische Universit¨at Berlin

Institut f¨ur Mathematik

Linear Discrete-Time Descriptor

Systems

Tobias Br¨ull

Preprint 30/2007

Preprint-Reihe des Instituts f¨ur Mathematik

Technische Universit¨at Berlin

Report 30/2007 January 2008

Linear Discrete-Time Descriptor Systems

Diplomarbeit von

Tobias Br¨ull

betreut durch

Prof. Dr. V. Mehrmann

Technische Universit¨at Berlin,

Fakult¨at II f¨ur Mathematik und Naturwissenschaften,

Institut f¨ur Mathematik

Berlin, den 21. August 2007

Contents

1 Introduction 2

1.1 Applications.................................... 3

1.1.1 Discretization of linear differential algebraic equations . . . . . . . . . 4

1.1.2 Singular Leontief Systems . . . . . . . . . . . . . . . . . . . . . . . . 4

1.1.3 Backward Leslie Model . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.1.4 A self-excogitated example: The Bullwhip Effect . . . . . . . . . . . . 5

2 Linear discrete-time descriptor systems with constant coefficients 8

2.1 Solution with the Kronecker canonical form . . . . . . . . . . . . . . . . . . 8

2.2 Explicit representation of the solution . . . . . . . . . . . . . . . . . . . . . . 14

2.3 Conclusion..................................... 32

3 Linear discrete-time descriptor systems with variable coefficients 33

3.1 Canonicalforms.................................. 33

3.2 Forward global canonical form . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.2.1 Local and global invariants . . . . . . . . . . . . . . . . . . . . . . . . 66

3.2.2 Existence and uniqueness of solutions . . . . . . . . . . . . . . . . . . 74

3.3 Backward global canonical form . . . . . . . . . . . . . . . . . . . . . . . . . 75

4 Two-way global canonical form for linear descriptor systems 77

4.1 Existence and uniqueness of solutions . . . . . . . . . . . . . . . . . . . . . . 93

5 Algorithms for linear discrete-time descriptor systems 95

5.1 A method with the Drazin inverse . . . . . . . . . . . . . . . . . . . . . . . . 95

5.2 A method employing the GUPTRI-Algorithm . . . . . . . . . . . . . . . . . 98

5.3 A reduction method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

5.3.1 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

6 Summary 104

6.1 Outlook ...................................... 105

A Matlab Code 106

1

Chapter 1

Introduction

The purpose of this diploma thesis is to transfer the concepts and results presented in [8] into

the discrete-time case. In [8] continuous-time differential-algebraic equations of the form

F(t, x, ˙x) = 0,(1.1)

with F:I×Dx×D˙x→Cm,

where I⊂Rand Dx,D˙x⊂Cn,

are examined. Given such an equation one is looking for a differentiable function x:I→Cn

which solves (1.1) in the sense that F(t, x(t),˙x(t)) = 0 for all t∈I, where ˙xdenotes the

derivative of x. The analogous general form of a discrete-time descriptor system is

F(k, xk, xk+1) = 0, (1.2)

with F:K×Dxk×Dxk+1 →Cm,

where K:= {k∈Z:kb≤k≤kf},kb∈Z∪ {−∞},kf∈Z∪ {∞} and

Dxk,Dxk+1 ⊂Cn.

Such systems could also be called discrete-time systems of differential-algebraic equations,

since systems of the form (1.2) naturally arise by discretizing systems of the form (1.1)

through a difference quotient, e.g., ˙x(ti)≈x(ti+1)−x(ti)

ti+1−ti. Other more common names are

discrete-time singular systems (e.g., [16]), discrete-time semi-state systems and discrete-time

generalized state-space systems. One is looking for a sequence {yk}which solves (1.2) in the

sense that F(k, yk, yk+1) = 0 for all k∈K. Such a sequence is called a solution of (1.2).

Although not all concepts in [8] are reasonable in the discrete-time case (e.g., generalized

solutions) there are still many concepts which can be carried over into the discrete-time

case. To study all these concepts is a major task. For this reason, we only consider some

of the concepts in [8] for linear discrete-time systems, despite the importance of non-linear

discrete-time descriptor systems of the form (1.2).

We distinguish two types of linear discrete-time descriptor systems, namely systems with

constant and systems with variable coefficients. To introduce these two types we first define

the discrete interval

K:= {k∈Z:kb≤k≤kf},kb∈Z∪ {−∞},kf∈Z∪ {∞}. (1.3)

2

With this definition we call

Exk+1 =Axk+fk,k∈K, (1.4)

where E, A ∈Cm,n, fk∈Cm,

alinear discrete-time descriptor system with constant coefficients. Such systems represent a

special case of the second type

Ekxk+1 =Akxk+fk,k∈K, (1.5)

where Ek, Ak∈Cm,n, fk∈Cm,

which is called a linear discrete-time descriptor system with variable coefficients. Together

with one of the systems (1.4) or (1.5) we also often require a solution to satisfy an initial

condition

xk0= ˆx, where k0∈K. (1.6)

Other notation is listed in the following table.

diag (A1,...,Ak) Denotes the block diagonal matrix with A1,...,Akon

the block diagonal.

N{1,2,3,...}

N0N∪ {0}={0,1,2,3,...}

R+R∩]0; ∞] = ]0; ∞]

Zis a basis of the vector

space V

When talking of a matrix Zbeing a basis of some vector

space V⊂Rlthen this means that the columns of Z

form a basis of V.

xk

iFor xk∈Cn,xk

imeans the i-th component of xk.

{Ak}kf

k=kbA sequence of Ak∈Cm,n for k=kb,...,kfwhere

kb, kf∈Z.

{Ak}k≥kbA sequence of Ak∈Cm,n for k≥kb∈Z.

{Ak}k≤kfA sequence of Ak∈Cm,n for k≤kf∈Z.

{Ak}k∈KFor a subset K⊂Z, this expression means the sequence

of the Ak∈Cm,n for k∈K.

Table 1.1: Notation used in this text

1.1 Applications

In this section we will discuss some applications where linear discrete-time descriptor systems

are used.

3

1.1.1 Discretization of linear differential algebraic equations

Consider the general linear continuous-time differential algebraic equation

E(t) ˙x(t) = A(t)x(t) + f(t), t∈[t0, tf]. (1.7)

We define a grid t0< t1< . . . < tN=tfand introduce xi:= x(ti), fi:= f(ti), Ei:=

1

ti+1−tiE(ti) and Ai:= 1

ti+1−tiE(ti) + A(ti). Approximating

˙x(ti)≈xi+1 −xi

ti+1 −ti

,

in (1.7) yields

E(ti)xi+1 −xi

ti+1 −ti

=A(ti)x(ti) + f(ti), for all i∈ {0,...,N −1}.

This system is equivalent to

Eixi+1 =Aixi+fi, for all i∈ {0,...,N −1},

which is a linear discrete-time descriptor system with variable coefficients.

1.1.2 Singular Leontief Systems

Leontief systems [3, 14] have the form

xk=Axk+B

|{z}

:= ˜

E

(xk+1 −xk) + dk

|{z}

:=−fk

(1.8)

⇔˜

Exk+1 = (I−A+B)

|{z }

:= ˜

A

xk+fk

⇔˜

Exk+1 =˜

Axk+fk

where

A,B∈Rn,n,xk,dk∈Rn, n ∈N.

Such systems describe the production of an economy with ndistinct sectors. A widely used

example is n= 3 with the sectors agriculture, manufacturing, and service, i.e., the primary,

secondary, and tertiary sector of industry. Here xk

iis the (monetary) output of the sector

iin the time period k, whereas dk

iis the (monetary) customer demand for products of the

sector iin the time period kand is prescribed.

The term Axkis there to consider inter-sector relations in the economy, i.e., it takes into

account that any sector may need output from all the other sectors to produce its output.

For example, it seems reasonable that any produced good of the secondary sector requires

4

some service (like telecommunications) from the tertiary sector. In this case the entry a3,2

of the matrix Awould have to be positive. In this way the (also often considered) equation

xk=Axk+dk

arises, which describes the output of the sectors in dependence of the customer demand.

This is a linear equation and the matrix Ais called consumption matrix.

Finally, the term B(xk+1 −xk) may describe investments. When the output level in sector

iincreases from period kto k+ 1 there may be investments necessary to accomplish this

increase. For example, a growth in the primary sector demands for an increase of the

agricultural machinery. So when the primary sector grows from time period kto k+ 1, it

is necessary to produce this machinery a priori in period k. In this case the entry b2,1of

the matrix Bwould have to be positive. In this way equation (1.8) arises. This is a linear

discrete-time descriptor system. The matrix Bis called the capital coefficients matrix.

The matrix B(and thus ˜

E) may be singular, as stated in [3, 14], since some rows of Bmay

only contain zero elements. For example, a growth in any of the three sectors of industry

does not really demand for an increase in the production of the primary sector (i.e., of food).

At least the increase is (monetarily) inferior which means that the first row of Bis zero.

1.1.3 Backward Leslie Model

A Leslie Model [2] has the form

xk+1 =Txk,

where

T∈Rn,n,xk∈Rn,n∈N.

The model describes the age distribution of a given population in time. The population is

divided into ndistinct age classes. Here xk

iis the number of individuals in age class iin time

period k. Considering the individual birth and death rates of the age classes, the matrix Tis

constructed. As shown in [2] this matrix may be singular. If one wants to determine an age

distribution in the past, given a present age distribution, one has to solve the Leslie Model

backwards, i.e., one has to solve a system of the form

Txl+1 =xl.

Obviously the same situation occurs every time one wants to determine a state in the past

given a present state if the behavior of the system is given by a difference equation with

singular T.

1.1.4 A self-excogitated example: The Bullwhip Effect

Consider a chain of four warehouses where each warehouse may send and receive goods to

and from the preceding warehouse. The last warehouse accepts its orders from a customer.

The first warehouse issues orders to a manufacturer and may also send goods back to the

5

Figure 1.1: Sketch of the supply chain

manufacturer. Let xk

i,i∈ {1,2,3,4}be the stock in warehouse iin time period k. Further,

let yk

i,i∈ {0,1,2,3}be the amount of goods moved from warehouse ito warehouse i+ 1

(where warehouse 0 is the manufacturer) after the time period k(and thus before time period

k+1). A negative value of yk

imeans that goods were actually sent back from warehouse i+1

to warehouse i. Finally, let yk

4be the demand of the customer which is satisfied after time

period k(and thus before time period k+ 1). Assume that this demand is given externally

by a (non-negative) sequence fk∈R. From this setup one obviously gets the five equations

xk+1

i=xk

i+yk

i−1−yk

ifor i∈ {1,2,3,4},

yk

4=fk.

Since there are nine variables to be described there are four equations missing. Therefore

suppose that the manager of every warehouse iin every time period ktries to have a stock

in his warehouse that is equal to the amount of goods moved from warehouse ito warehouse

i+ 1 (=yk

i) after the time period k(and thus before time period k+ 1) plus a safety stock,

which he chooses proportional to yk

i. This assumption then leads to the four additional

equations

xk

i=ayk

ifor i∈ {1,2,3,4},

with a > 1.

Using the state vector xk=yk

0xk

1yk

1xk

2yk

2xk

3yk

3xk

4yk

4Tand the matrices

E:=







0

1

0

1

0

1

0

1

0







A:=







0−1a0 0 0 0 0 0

1 1 −1 0 0 0 0 0 0

0 0 0 −1a0 0 0 0

0 0 1 1 −1 0 0 0 0

0 0 0 0 0 −1a0 0

0 0 0 0 1 1 −1 0 0

0 0 0 0 0 0 0 −1a

0 0 0 0 0 0 1 1 −1

0 0 0 0 0 0 0 0 −1







,

one can write the system as

Exk+1 =Axk+00000000fkT.

Since the matrix Ais regular (which can be seen by exchanging adjacent rows), the pencil

(E, A) is regular and one can compute the unique solution, given a sequence of fk. In Figure

6

1.2 the sequence fkand the corresponding xk

ifor i∈ {1,2,3,4}of such a solution are plotted.

One sees that the further the warehouses are away from the customer the bigger the stock

fluctuations are. This is called the ”Bullwhip Effect” by economists.

Note that this example makes the assumption of complete information. (It has a high index

(I calculated an index of 5) and only infinite eigenvalues. Also it is possible for stocks to

become negative. All this shows that the example is of no practical use.)

0 5 10 15 20 25 30 35 40

4

5

6

7

8

9

10

−x−0−0−0−

−0−x−0−0−

−0−0−x−0−

−0−0−0−x−

Demand

Figure 1.2: Stock development and demand

7

Chapter 2

Linear discrete-time descriptor

systems with constant coefficients

It would be ideal to develop a theory for the most general type of discrete-time descriptor

systems (1.2). Nevertheless, such a theory is unlikely to exist. Thus, restrictions have to be

imposed on (1.2) in order to get a nice theory. Obviously, the theory gets nicer the more

restrictions one imposes on the problem.

It seems reasonable to start by imposing very strong restrictions on (1.2) and then loosen

the restrictions as one proceeds. The advantage in doing so is that one may first derive some

basic results, which than can be generalized.

Here we start with the very simple case of linear discrete-time descriptor systems with

constant coefficients (i.e., systems of the type (1.4) ), analogous to [8].

2.1 Solution with the Kronecker canonical form

We begin by recapitulating some basic results of linear algebra, that are needed afterwards.

First let us review the Kronecker canonical form.

Theorem 2.1. [8] Let E, A ∈Cm,n. Then there exist nonsingular matrices P∈Cm,m and

Q∈Cn,n such that for all λ∈C

P(λE −A)Q= diag Lǫ1,...,Lǫp,Mη1,...,Mηq,Jρ1,...,Jρr,Nσ1,...,Nσs,(2.1)

where the diagonal blocks have the following forms:

1. Every Lǫiblock is of size ǫi×(ǫi+ 1), ǫi∈N0and has the form

Lǫj=λ





0 1

......

0 1





−





1 0

......

1 0





.(2.2)

8

2. Every Mηjblock is of size (ηj+ 1) ×ηj, ηj∈N0and has the form

Mηj=λ





1

0...

...1

0







−





0

1...

...0

1







.(2.3)

3. Every Jρkblock is of size ρk×ρk, ρk∈Nand has the form

Jρk=λ





1

...

1







−





λk1

......

...1

λk







.(2.4)

4. Every Nσlblock is of size σl×σl, σl∈Nand has the form

Nσl=λ





0 1

......

...1

0







−





1

...

1







.(2.5)

The Kronecker canonical form is unique up to permutation of the blocks, i.e., the kind, size

and number of the blocks are invariant for the matrix pair (E, A).

Definition 2.2. Let E, A ∈Cn,n. Then the matrix pair (E, A) is called regular if and only

if det λE −Adoes not vanish identically.

Theorem 2.3. [8] Let E, A ∈Cn,n and (E, A)be regular. Then there exist nonsingular

matrices P, Q ∈Cn,n such that for all λ∈Cwe have

P(λE −A)Q=λI0

0N−J0

0I,(2.6)

where J is a matrix in Jordan canonical form and N is a nilpotent matrix also in Jordan

canonical form. Moreover, it is allowed that one or the other block is not present.

Thus, the Kronecker canonical form of a regular matrix pencil only has blocks of type (2.4)

and (2.5).

Definition 2.4. [8] Let E, A ∈Cn,n, let the matrix pair (E, A) be regular and let the

Kronecker canonical form of (E, A) be given by (2.6). Then the quantity νdefined by

Nν= 0, Nν−16= 0, i.e., by the index of nilpotency of N in (2.6), if the nilpotent block in

(2.6) is present and by ν= 0 if it is absent, is called the index of the matrix pair (E, A),

denoted by ν= ind(E, A).

9

Definition 2.5. Let E∈Cn,n. Further, let νbe the index of the matrix pair (E, I). Then

νis also called the index of E and denoted by ind(E) = ν.

Consider an arbitrary matrix pencil λE −Awith the Kronecker canonical from (2.1). When

we are interested in the solution of the associated discrete-time descriptor system consisting

of (1.4) and (1.6) with kb=k0∈Zwe can study the problem in the coordinates of the

Kronecker canonical form, i.e., we can look at the equivalent problem

PEQ˜xk+1 =PAQ˜xk+P fk, for k∈K,

˜xk0=Q−1ˆx,

with ˜xk=Q−1xk. Since the pencil (PEQ, PAQ) is block diagonal one can compute the

solution for each block separately. This is done in the following.

1. Consider a block of type (2.2), i.e., let

λEL−AL=λ





0 1

......

0 1





−





1 0

......

1 0





∈Cǫ,ǫ+1.

Then the system (1.4) corresponding to this pencil, together with (1.6), i.e.,

(ELxk+1 =ALxk+fk,k≥k0

xk0= ˆx,(2.7)

is equivalent to (xk+1

i+1 =xk

i+fk

i,i= 1,...,ǫ,

xk0= ˆx.(2.8)

If we have ǫ= 0 we have no equations but one variable. Thus, in this case every

sequence satisfying the initial condition is a solution. Therefore let us assume in the

following that ǫ > 0. In this case recursion of (2.8) yields

xk

i=xk−1

i−1+fk−1

i−1

=xk−2

i−2+fk−2

i−2+fk−1

i−1

=...=



min(i−1,k−k0)

X

j=1

fk−j

i−j

+xk−min(i−1,k−k0)

i−min(i−1,k−k0)

=









Pi−1

j=1 fk−j

i−j+xk−i+1

1, if k−k0> i −1,

Pk−k0

j=1 fk−j

i−j+ ˆxi−k+k0, if k−k0≤i−1,

for k=k0,...,kf. Thus, once the initial condition ˆx, values for xk0+1

1, xk0+2

1,...and the

right hand sides fk0, fk0+1,... are given one can uniquely determine the solution. For

10

a given initial condition and right hand side, there are still kf−k0degrees of freedom,

i.e., the solution space has dimension kf−k0.

Written with matrices the following result is obtained. Multiplying (2.7) with ET

Lfrom

the left results in

ET

LELxk+1 =ET

LALxk+ET

Lfk, (2.9)

where

ET

LEL=0 0

0Iǫ1

ǫis a projector and

ET

LAL=0 0

Iǫ01

ǫis nilpotent with nilpotency index ǫ+ 1.

Partitioning the vector xk=xk+1

1

xk+1

21

ǫin (2.9) results in

0

xk+1

2=ET

LALxk+ET

Lfk

⇔xk+1 =ET

LALxk+ET

Lfk+





xk+1

1

0

.

0







.

Recursion of this equation leads to

xk=ET

LALk−k0ˆx+

k−k0−1

X

j=0 ET

LALj





ET

Lfk−1−j+





xk−j

1

0

.

0













.

2. Consider a block of type (2.3), i.e., let

λEM−AM=λ





1

0...

...1

0







−





0

1...

...0

1







∈Cη+1,η.

In the case η= 0 this means that we are looking at sequences in which each iterate

is a vector of length 0. Such sequences are of no interest. Therefore let us assume in

the following that η > 0. In this case the system (1.4) corresponding to this pencil,

together with (1.6), i.e.,

(EMxk+1 =AMxk+fk,k≥k0

xk0= ˆx,(2.10)

11

is equivalent to 









xk+1

1=fk

1,

xk+1

i=xk

i−1+fk

i,i= 2,...,η,

0 = xk

η+fk

η+1,

xk0= ˆx.

Recursion of the second equation and the first equation yields

xk

i=xk−1

i−1+fk−1

i

=xk−2

i−2+fk−2

i−1+fk−1

i

=...=



min(i−1,k−k0)

X

j=1

fk−j

i+1−j

+xk−min(i−1,k−k0)

i−min(i−1,k−k0)

=









Pi−1

j=1 fk−j

i+1−j+xk−i+1

1

|{z}

=fk−i

1

, if k−k0> i −1,

Pk−k0

j=1 fk−j

i+1−j+xk0

i−k+k0, if k−k0≤i−1,

=









Pi

j=1 fk−j

i+1−j, if k−k0> i −1,

Pk−k0

j=1 fk−j

i+1−j+ ˆxk0

i−k+k0, if k−k0≤i−1.

Because of the additional equation xk

η=−fk

η+1 the previous equation implies that

−fk

η+1 =xk

η=









Pη

j=1 fk−j

η+1−j, if k−k0> η −1,

Pk−k0

j=1 fk−j

η+1−j+ ˆxk0

η−k+k0, if k−k0≤η−1.

This provides the consistency condition for the inhomogeneity

0 =

η

X

j=0

fk−j

η+1−jfor all k > k0+η−1,

and the consistency condition for the initial condition

ˆxk0

η−k+k0=−

k−k0

X

j=0

fk−j

η+1−jfor all k0≤k≤k0+η−1.

Replacing i=η−k+k0in the last equation yields the nicer form

ˆxk0

i=−

η−i

X

j=0

fη−i+k0−j

η+1−jfor all 1 ≤i≤η.

12

Again one can get the same result in terms of matrices. Multiplying (2.10) with ET

M

from the left leads to

ET

MEM

|{z }

=I

xk+1 =ET

MAMxk+ET

Mfk

⇒xk+1 =ET

MAMxk+ET

Mfk

= (ET

MAM)2xk−1+ (ET

MAM)ET

Mfk−1+ET

Mfk

=...= (ET

MAM)k−k0+1 ˆx+

k−k0

X

j=0

(ET

MAM)jET

Mfk−j. (2.11)

One should notice, that

(ET

MAM) = 





1 0

......

1 0











0

1...

...0

1







=





0

1...

......

1 0







∈Cη,η,

is nilpotent with nilpotency index η. Thus, (2.11) can also be written as

xk+1 = (ET

MAM)k−k0+1 ˆx+

min(k−k0,η−1)

X

j=0

(ET

MAM)jET

Mfk−j.

3. Consider a block of type (2.4), i.e., let

λEJ−AJ=λ





1

...

1







−





λk1

......

...1

λk







∈Cρ,ρ.

Then the system (1.4) corresponding to this pencil, together with (1.6) is

(xk+1 =AJxk+fk,k≥k0

xk0= ˆx.(2.12)

Simple recursion yields

xk=AJxk−1+fk−1

=A2

Jxk−2+AJfk−2+fk−1

=...=Ak−k0

Jˆx+

k−k0

X

i=1

Ai−1

Jfk−i.

13

4. Consider a block of type (2.5), i.e., let

λEN−AN=λ





0 1

......

...1

0







−





1

...

1







∈Cσ,σ.

Then the system (1.4) corresponding to this pencil, together with (1.6) is

(ENxk+1 =xk+fk,k≥k0

xk0= ˆx.(2.13)

It follows that

xk=ENxk+1 −fk

=E2

Nxk+2 −ENfk+1 −fk

=...=−

σ−1

X

i=0

Ei

Nfk+i,

because the matrix ENis nilpotent with index σ. Thus, in this case we see that the

solution only depends on present and future right hand sides.

2.2 Explicit representation of the solution

In order to determine an explicit solution of (1.4) one can employ the Drazin inverse.

Definition 2.6. Let E∈Cn,n have the index ν. A matrix X∈Cn,n satisfying

EX =XE, (2.14)

XEX =X, (2.15)

XEν+1 =Eν, (2.16)

is called a Drazin inverse of Eand denoted by ED.

From this definition some basic results can be derived.

Lemma 2.7. Consider matrices E, A ∈Cn,n with EA =AE. Then

EAD=ADE,

EDA=AED,(2.17)

EDAD=ADED

where EDdenotes the Drazin inverse of E.

14

Proof. [8], Theorem (2.21), p. 25.

Also, recall the following Theorem.

Theorem 2.8. Let E∈Cn,n with ν= ind(E). There is one and only one decomposition

E=˜

C+˜

N(2.18)

with the properties

˜

C˜

N=˜

N˜

C= 0,˜

Nν= 0,˜

Nν−16= 0,ind( ˜

C)≤1.(2.19)

For this decomposition the following statements hold:

˜

CD˜

N=˜

N˜

CD= 0,(2.20a)

ED=˜

CD,(2.20b)

˜

C˜

CD˜

C=˜

C,(2.20c)

˜

CD˜

C=EDE,(2.20d)

˜

C=EEDE,˜

N=EI−EDE.(2.20e)

Proof. [8], Theorem (2.22), p. 25.

With these preliminaries one can start like in [8] to find solutions in the special case that

the matrices Eand Acommute.

Lemma 2.9. Let E, A ∈Cn,n with EA =AE and E=˜

C+˜

Nbe the decomposition (2.18).

Then the following propositions hold.

1. Let {xk}kf+1

k=k0be a solution of (1.4). Set

xk

1:= EDExk,xk

2:= I−EDExk,(2.21)

fk

1:= EDEfk,fk

2:= I−EDEfk.

Then, we have

˜

Cxk+1

1=Axk

1+fk

1,(2.22)

˜

Nxk+1

2=Axk

2+fk

2,(2.23)

for k=k0,...,kf.

2. Let {xk

1}kf+1

k=k0and {xk

2}kf+1

k=k0be solutions of

˜

Cxk+1

1=Axk

1+EDEfk,(2.24)

˜

Nxk+1

2=Axk

2+ (I−EDE)fk,

for k=k0,...,kf, respectively. Then {xk}kf+1

k=k0given by

xk:= EDExk

1+ (I−EDE)xk

2

is a solution of (1.4).

15

3. Let {xk

1}kf+1

k=k0have the form

xk

1=EDExk,(2.25)

for some {xk}kf+1

k=k0and let {fk

1}kf+1

k=k0have the form

fk

1=EDEfk,(2.26)

for some {fk}kf+1

k=k0. Then {xk

1}kf+1

k=k0is a solution of

˜

Cxk+1

1=Axk

1+fk

1(2.27)

if and only if {xk

1}kf+1

k=k0is a solution of

xk+1

1=EDAxk

1+EDfk

1.(2.28)

Proof. 1. From (2.20) and (2.17) it follows that

˜

CD˜

CA (2.20d)

=EDEA (2.17)

=AEDE=A˜

CD˜

C(2.29)

and

˜

Nxk

1=˜

NEDExk(2.20d)

=˜

N˜

CD˜

Cxk(2.20a)

= 0, (2.30)

˜

Cxk

2=˜

C−˜

C˜

CD˜

Cxk(2.20c)

= 0, (2.31)

˜

Nfk

1=˜

NEDEfk(2.20d)

=˜

N˜

CD˜

Cfk(2.20a)

= 0,

˜

Cfk

2=˜

C−˜

C˜

CD˜

Cfk(2.20c)

= 0. (2.32)

for all k=k0,...,kf+ 1. Further, from (1.4) we obtain

˜

C+˜

Nxk+1

1+xk+1

2=Axk

1+xk

2+fk

1+fk

2(2.33)

(2.29)

⇒(˜

CD˜

C˜

C+˜

CD˜

C˜

N

|{z}

(2.19)

= 0

)xk+1

1+xk+1

2=A˜

CD˜

Cxk

1+xk

2+˜

CD˜

Cfk

1+fk

2

(2.31),(2.32)

⇒˜

CD˜

C˜

Cxk+1

1

(2.21)

=A˜

CD˜

Cxk

1+fk

1=AEDExk

1+fk

1

(2.21)

=Axk

1+fk

1,

for all k=k0,...,kf, which shows (2.22), since

˜

CD˜

C˜

C(2.14)

=˜

C˜

CD˜

C(2.20c)

=˜

C.

Subtracting (2.22) from (2.33) yields

˜

Cxk+1

2

|{z}

(2.31)

= 0

+˜

Nxk+1

1

|{z }

(2.30)

= 0

+˜

Nxk+1

2=Axk

2+fk

2,

for all k=k0,...,kf, which shows (2.23).

16

2. Applying the Definition of xk+1 leads to

Exk+1 =EEDExk+1

1+E(I−EDE)xk+1

2

(2.20e)

=˜

Cxk+1

1+˜

Nxk+1

2

(2.24)

=A(xk

1+xk

2) + EDEfk+ (I−EDE)fk

=Axk+fk.

3. Multiplying (2.27) with ˜

CD(2.20b)

=EDfrom the left leads to

˜

CD˜

Cxk+1

1=EDAxk

1+EDfk

1. (2.34)

From (2.25) one can also obtain

(I−˜

CD˜

C)xk+1

1= 0. (2.35)

Adding (2.34) and (2.35) then immediately shows (2.28). Conversely, multiplying

(2.28) by ˜

Cfrom the left gives

˜

Cxk+1

1=˜

C˜

CDAxk

1+˜

C˜

CDfk

1

=AEDExk

1+EDEfk

1

(2.25),(2.26)

=Axk

1+fk

1,

since

˜

C˜

CD(2.20e),(2.20b)

=EEDEED(2.17)

=EDEEDE(2.15)

=EDE.

Lemma 2.10. Let E, A ∈Cn,n with EA =AE,k0∈Zand v∈Cn. Then the following

statements hold.

1. Let ˆv=EDEv. Then

xk:= (EDA)k−k0ˆv,k=k0, k0+ 1,... (2.36)

solves the homogeneous linear discrete-time descriptor system

Exk+1 =Axk,k=k0, k0+ 1,... (2.37)

2. Let ˆv=ADAv. Then

xk:= (ADE)k0−kˆv,k=k0, k0−1,... (2.38)

solves the homogeneous linear discrete-time descriptor system

Exk+1 =Axk,k=k0−1, k0−2,... (2.39)

17

3. Let ˆv∈range ADA∩range EDE. Then

xk:= ((EDA)k−k0ˆv,k=k0, k0−1,...

(ADE)k0−kˆv,k=k0−1, k0−2,... (2.40)

solves the homogeneous linear discrete-time descriptor system

Exk+1 =Axk,k∈Z.(2.41)

Proof. 1. We have

Exk+1 =E(EDA)(EDA)k−k0EDEv

(2.17)

=A(EDA)k−k0EDEEDEv

=A(EDA)k−k0EDEv

=Axkfor all k=k0, k0+ 1,... .

2. In this case we obtain

Axk=A(ADE)k0−kADAv

=A(ADE)(ADE)k0−k−1ADAv

(2.17)

=E(ADE)k0−k−1ADAADAv

=E(ADE)k0−k−1ADAv

=Exk+1 for all k=k0−1, k0−2,... .

3. This follows from 1. and 2., since the definitions of xk0from 1. and 2. coincide.

Since

(EDA)k−k0EDEv =EDE(EDA)k−k0v,

it is clear, that the solution xkstays in the subspace range EDEfor all k≥k0. An

analogous conclusion is possible for the case 2. in Lemma 2.10. In case 3. of Lemma 2.10

the solution even stays in range ADA∩range EDEall the time.

Theorem 2.11. Let E, A ∈Cn,n with EA =AE and suppose that

kernel (E)∩kernel (A) = {0}.(2.42)

Then,

(I−EDE)ADA= (I−EDE).(2.43)

Proof. [8], Theorem (2.28), p. 30.

18

Theorem 2.12. Let E, A ∈Cn,n with EA =AE satisfy (2.42). Also, let k0∈Z. Then the

following statements hold.

1. Let {xk}k≥k0be any solution of (2.37). Then {xk}k≥k0has the form (2.36) for some

ˆv∈range EDE.

2. Let {xk}k≤k0be any solution of (2.39). Then {xk}k≤k0has the form (2.38) for some

ˆv∈range ADA.

3. Let {xk}k∈Zbe any solution of (2.41). Then {xk}k∈Zhas the form (2.40) for some

ˆv∈range ADA∩range EDE.

Proof. Using the decomposition (2.18) we get the following results.

1. We have

A˜

N(2.20e)

=AE(I−EDE)(2.17)

=E(I−EDE)A(2.20e)

=˜

NA. (2.44)

Furthermore for any x∈Cnone has

A˜

Nx = 0 ⇒ADA˜

Nx = 0

⇒(I−EDE)ADA˜

Nx = 0

(2.43)

⇒(I−EDE)˜

Nx = 0 (2.45)

(2.20e)

⇒˜

Nx = 0.

Let {xk}k∈Zbe any solution of (2.37). From Lemma 2.9 part 1. we get {xk

1}k≥k0,

{xk

2}k≥k0with xk=xk

1+xk

2which solve (2.22) and (2.23), respectively. With ν= ind(E)

one then obtains

0(2.19)

=˜

Nνxk+1

2

(2.23)

=˜

Nν−1Axk

2

(2.44)

=A˜

Nν−1xk

2,k≥k0

(2.45)

⇒˜

Nν−1xk

2= 0, k≥k0.

Discarding the identity for k=k0then yields

˜

Nν−1xk

2= 0, k≥k0+ 1.

Shift

⇒˜

Nν−1xk+1

2= 0, k+ 1 ≥k0+ 1

⇒˜

Nν−1xk+1

2= 0, k≥k0

⇒...⇒˜

Nxk

2= 0, k≥k0

(2.23)

⇒Axk

2= 0, k≥k0(2.46)

⇒xk

2

(2.21)

= (I−EDE)xk

2

(2.43)

= (I−EDE)ADAxk

2

|{z}

(2.46)

= 0

= 0, k≥k0

19

⇒xk=xk

1,k≥k0

Lemma 2.9 3.

⇒xk

1solves xk+1

1= (EDA)xk

1,k≥k0

Recursion

⇒xk

1= (EDA)k−k0xk0

1,k≥k0

⇒xk=xk

1= (EDA)k−k0xk0

1

(2.21)

= (EDA)k−k0EDExk0. (2.47)

2. Let {xk}k≤k0be any solution of (2.39). Set

l0:= −k0and yl:= x−l,l≥l0.

By replacing k=−lin (2.39) one obtains

Ex−l+1 =Ax−l,−l=−l0−1,−l0−2,...

⇒Ex−(l−1) =Ax−l,l=l0+ 1, l0+ 2,...

⇒yl,l≥l0is a solution of Eyl−1=Ayl,l≥l0+ 1

⇒yl,l≥l0is a solution of Ayl+1 =Eyl,l≥l0

(2.47)

⇒yl= (ADE)l−l0ADAyl0,l≥l0.

Undoing the replacements then yields

x−l= (ADE)l−l0ADAx−l0,l≥l0,

xk= (ADE)−k+k0ADAxk0,−k≥ −k0,

xk= (ADE)k0−kADAxk0,k≤k0. (2.48)

3. Let {xk}k∈Zbe any solution of (2.41). Then from (2.47) we have

xk= (EDA)k−k0EDExk0,k≥k0

k=k0

⇒xk0=EDExk0⇒xk0∈range EDE.

Also we know from (2.48) that

xk= (ADE)k0−kADAxk0,k≤k0

k=k0

⇒xk0=ADAxk0⇒xk0∈range ADA.

Thus, the claim of the Theorem follows with ˆv=xk0.

Remark 2.13. One may think that it is not meaningful to look at case 3. of the previous

Theorem, since in most cases one starts at some time point and then calculates into the

future. But as shown by the following Lemma, also those solutions (where one starts at k0∈

Zand calculates a solution for k≥k0) are almost completely in the subspace range ADA∩

range EDE.

20

Lemma 2.14. Let E, A ∈Cn,n with EA =AE satisfy (2.42). Also, let k0∈Zand let

νE= ind(E),νA= ind(A). Then the following statements hold.

1. Let {xk}k≥k0be any solution of (2.37). Then for all k≥k0+νAit holds that xk∈

range ADA∩range EDE.

2. Let {xk}k≤k0be any solution of (2.39). Then for all k≤k0−νEit holds that xk∈

range EDE∩range ADA.

Proof. 1. Since k≥k0+νAit follows that k=ˆ

k+k0+νAwith ˆ

k≥0. From Theorem

2.12 we then know that for some v∈Cnwe have

ADAxk(2.36)

=ADA(EDA)k−k0

|{z }

=(ED)k−k0Ak−k0

EDEv

=ADA(ED)k−k0AνAAˆ

kEDEv

(2.16)

=ADA(ED)k−k0ADAνA+1Aˆ

kEDEv

= (ED)k−k0ADAAD

|{z }

(2.15)

=AD

AνA+1Aˆ

kEDEv

= (ED)k−k0ADAνA+1Aˆ

kEDEv

(2.16)

= (EDA)k−k0EDEv

=xk.

Also, we naturally get

EDExk(2.36)

=EDE(EDA)k−k0EDEv

= (EDA)k−k0EDEED

|{z }

(2.15)

=ED

Ev (2.49)

=xk,

and thus the assertion follows.

2. As in (2.49) one gets that ADAxk=xk. Let k=−ˆ

k+k0−νEwith ˆ

k≥0. Then again

for some v∈Cnit follows that

EDExk(2.38)

=EDE(ADE)k0−k

|{z }

=(AD)k0−kEk0−k

ADAv

=EDE(AD)k0−kEνEEˆ

kADAv

(2.16)

=EDE(AD)k0−kEDEνE+1Eˆ

kADAv

= (AD)k0−kEDEED

|{z }

(2.15)

=ED

EνE+1Eˆ

kADAv

21

= (AD)k0−kEDEνE+1Eˆ

kADAv

(2.16)

= (ADE)k0−kADAv

=xk.

Remark 2.15. From Lemma 2.14 one might presume that it is meaningful to require that

the initial condition satisfies

xk0∈range ADA∩range EDE, (2.50)

since only in this case it is possible to calculate the solution into the future (i.e., calculate

xkfor k≥k0) and into the past (i.e., calculate xkfor k≤k0).

Also, only in case that (2.50) holds, we get something like an invertibility of the operator

that calculates xk+1 from xk. To understand this, imagine that a fixed xk0is given. From

this we calculate a finite number of steps κinto the future. Thus, we have xk0+κ. From this

state we then calculate κsteps back into the past to obtain ˜xk0. We then have xk0= ˜xk0if

condition (2.50) holds. Otherwise we cannot be sure that xk0= ˜xk0holds, as shown in the

following example.

Example 2.16. Consider the homogeneous linear discrete-time descriptor system defined

by





1 0 0

0 1 0

0 0 0



|{z }

:=E

xk+1 =



000

010

001



|{z }

:=A

xk,k≥0, x0=



1

0

. (2.51)

Clearly, we have EA =AE,ED=E,AD=Aand condition (2.42) holds. Thus, the pencil

(E, A), corresponding to system (2.51), satisfies all assumptions from Lemma 2.14 which

means that the iterate x1has to be in range ADA. Indeed,

Ax0=



0

1

0

⇒x1=



0

1

0

∈range ADA. (2.52)

Now let us calculate back one step from (2.52), i.e., let us consider the reversed system

A˜xl+1 =E˜xl,l≤0, ˜x−1=



0

1

0

.

We see that

E˜x−1=



0

1

0

⇒˜x0=



0

1

0

,

22

and thus

˜x0=



0

1

0

6=x0.

Theorem 2.17. Let E, A ∈Cn,n with EA =AE satisfy (2.42). Also, let νE= ind(E),

νA= ind(A),{fk}k∈Zwith fk∈Cnand k0∈Z. Then the following statements hold.

1. The linear discrete-time descriptor system

Exk+1 =Axk+fk,k≥k0

has the particular solution {xk}k≥k0with

xk:=

k−1

X

j=k0

(EDA)k−j−1EDfj

|{z }

:=xk

1

−(I−EDE)

νE−1

X

i=0

(ADE)iADfk+i

|{z }

:=xk

2

for k≥k0.(2.53)

For the construction of the iterate xkonly the fkwith k≥k0have to be employed.

2. The linear discrete-time descriptor system

Exk+1 =Axk+fk,k≤k0−1 (2.54)

has the particular solution {xk}k≤k0with

xk:= (I−ADA)

νA−1

X

i=0

(EDA)iEDfk−i−1−

k0

X

j=k+1

(ADE)j−k−1ADfj−1for k≤k0.

(2.55)

For the construction of the iterate xkonly the fkwith k≤k0−1have to be employed.

Proof. 1. Let E=˜

C+˜

Nbe the decomposition (2.18). Then

EDExk

1

(2.17)

=

k−1

X

j=k0

(EDA)k−j−1EDEED

|{z }

=ED

fj=xk

1,

(I−EDE)xk

2=−(I−EDE)(I−EDE)

νE−1

X

i=0

(ADE)iADfk+i=xk

2.

One can also conclude, that for all k≥k0it follows that

˜

Cxk+1

1=˜

C

k

X

j=k0

(EDA)k+1−1−jEDfj

23

=˜

C k−1

X

j=k0

(EDA)k−jEDfj+EDfk!

=˜

C (EDA)

k−1

X

j=k0

(EDA)k−j−1EDfj+EDfk!

(2.20e),(2.15)

=AEDExk

1+EDEfk

=Axk

1+EDEfk,

and with

(I−EDE)EνE=((I−EDE)EDEνE−1= 0 , if νE≥1,

(I−EDE) = (I−I) = 0 , if νE= 0, (2.56)

we obtain

˜

Nxk+1

2

(2.20e)

=E(I−EDE)xk+1

2

=−(I−EDE)

νE−1

X

i=0

(ADE)i+1fk+i+1

(2.56)

=−(I−EDE)

νE−2

X

i=0

(ADE)i+1fk+i+1

(2.43)

=−(I−EDE)ADA

νE−1

X

i=1

(ADE)ifk+i

=−A(I−EDE)

νE−1

X

i=0

(ADE)iADfk+i+ (I−EDE)fk

=Axk

2+ (I−EDE)fk.

With these results and Lemma 2.9 part 2 one immediately gets that

xk=EDExk

1+ (I−EDE)xk

2=xk

1+xk

2

is a solution and thus the assertion follows.

2. By replacing l:= −kand l0:= −k0in (2.54) one gets the system

Ex−l+1 =Ax−l+f−l,−l≤ −l0−1

⇒Ex−(l−1) =Ax−l+f−l,l≥l0+ 1.

By further replacing yl:= x−l, gl:= −f−l−1for l≥l0one gets

Eyl−1=Ayl+f−l,l≥l0+ 1

⇒Eyl=Ayl+1 +f−l−1,l+ 1 ≥l0+ 1

⇒Ayl+1 =Eyl−f−l−1,l≥l0

⇒Ayl+1 =Eyl+gl,l≥l0.

24

We then get a solution of this last system as

1.

⇒yl=

l−1

X

j=l0

(ADE)l−j−1ADgj−(I−ADA)

νA−1

X

i=0

(EDA)iEDgl+i

⇒x−l=−

l−1

X

j=l0

(ADE)l−j−1ADf−j−1+ (I−ADA)

νA−1

X

i=0

(EDA)iEDf−(l+i)−1

⇒xk= (I−ADA)

νA−1

X

i=0

(EDA)iEDfk−i−1−

−k−1

X

j=−k0

(ADE)−k−j−1ADf−j−1

= (I−ADA)

νA−1

X

i=0

(EDA)iEDfk−i−1−

k0

X

j=k+1

(ADE)j−k−1ADfj−1.

Finding a particular solution for the general case (i.e., (1.4) with K=Z) is more complicated.

Similar to Lemma 2.9 we obtain the following result.

Lemma 2.18. Let E, A ∈Cn,n with EA =AE satisfy (2.42). Further, let E=˜

C+˜

Nand

analogously A=˜

D+˜

Mbe decompositions as in (2.18). Let {xk

1}k∈Z,{xk

2}k∈Z,{xk

3}k∈Zbe

solutions of

˜

Cxk+1

1=˜

Mxk

1+ (I−ADA)fk,(2.57)

˜

Cxk+1

2=˜

Dxk

2+ADAEDEfk,(2.58)

˜

Nxk+1

3=˜

Dxk

3+ (I−EDE)fk,(2.59)

respectively. Then {xk}k∈Zwith

xk:= (I−ADA)xk

1+ADAEDExk

2+ (I−EDE)xk

3,

is a solution of

Exk+1 =Ak+fk.

Proof. First of all, we see that

(I−ADA) + (I−EDE) + ADAEDE

=I−ADA+I−(I−ADA)EDE

|{z }

(2.43)

= (I−ADA)

=I−ADA+I−(I−ADA) = I.

(2.60)

Furthermore, we have

˜

D(I−ADA) = AADA(I−ADA) = 0 , (2.61)

˜

M(ADAEDE) = A(I−ADA)(ADAEDE) = 0 , (2.62)

˜

M(I−EDE) = A(I−ADA)(I−EDE) = A(I−EDE)−(I−EDE)ADA= 0. (2.63)

25

With these identities we get

Exk+1 =E(I−ADA)

|{z }

(2.43)

= (I−ADA)EDE

xk+1

1+EADA EDE

|{z}

(2.15)

=EDEEDE

xk+1

2+E(I−EDE)

|{z }

=(I−EDE)(I−EDE)

xk+1

3

= (I−ADA)˜

Cxk+1

1+ADAEDE˜

Cxk+1

2+ (I−EDE)˜

Nxk+1

3

= (I−ADA)˜

Mxk

1+ (I−ADA)fk+

ADAEDE˜

Dxk

2+ADAEDEfk+

(I−EDE)˜

Dxk

3+ (I−EDE)fk

(2.60)

=fk+ (I−ADA)˜

Mxk

1+ADAEDE˜

Dxk

2+ (I−EDE)˜

Dxk

3

=fk+

˜

M(I−ADA)xk

1+˜

D(I−ADA)xk

1

|{z }

(2.61)

= 0

˜

MADAEDExk

2

|{z }

(2.62)

= 0

+˜

DADAEDExk

2

˜

M(I−EDE)xk

3

|{z }

(2.63)

= 0

+˜

D(I−EDE)xk

3

=fk+Axk.

Here we have used, that ˜

D=AADAand thus ˜

Dcommutes with the matrices Eand A.

Using Lemma 2.18 we can construct a particular solution for the general case (i.e., (1.4) with

K=Z).

Lemma 2.19. Let E, A ∈Cn,n with EA =AE satisfying (2.42). Also, let νE= ind(E),

νA= ind(A),{fk}k∈Z⊂Cnand k0∈Z. Then a solution {xk}k∈Zof

Exk+1 =Axk+fk

is given by

xk:= (I−ADA)

νA−1

X

i=0

(EDA)iEDfk−i−1

|{z }

=:xk

1

+(ADAPk−1

j=k0(EDA)k−j−1EDfjfor k≥k0

−EDEPk0

j=k+1(ADE)j−k−1ADfj−1for k≤k0

|{z }

:=xk

2

−(I−EDE)

νE−1

X

i=0

(ADE)iADfk+i

|{z }

:=xk

3

, for k∈Z.

(2.64)

26

Proof. Considering the decompositions E=˜

C+˜

Nand A=˜

D+˜

Mas in (2.18) we have

˜

Mxk

1=A(I−ADA)

νA−1

X

i=0

(EDA)iEDfk−i−1

(2.17)

= (I−ADA)

νA−1

X

i=0

(EDA)i+1fk−i−1

see (2.56)

= (I−ADA)

νA−2

X

i=0

(EDA)i+1fk−i−1

= (I−ADA)

νA−1

X

i=1

(EDA)ifk−i

= (I−ADA) νA−1

X

i=0

(EDA)ifk−i−fk!

(2.43)

=−(I−ADA)fk+ (I−ADA)EDE

νA−1

X

i=0

(EDA)ifk−i

(2.17)

=−(I−ADA)fk+E(I−ADA)

νA−1

X

i=0

(EDA)iEDfk−i

=−(I−ADA)fk+Exk+1

1

=−(I−ADA)fk+ ( ˜

C+˜

N)xk+1

1

=−(I−ADA)fk+˜

Cxk+1

1,

where the last identity holds, since xk

1has the form xk

1= (I−ADA)yk

1for some yk

1and

˜

Nxk

1=E(I−EDE)(I−ADA)yk

1=E



I−ADA−EDE(I−ADA)

|{z }

(2.43)

= (I−ADA)





yk

1= 0. (2.65)

As in Theorem 2.17, part 1. one obtains

˜

Nxk+1

3=Axk

3+ (I−EDE)fk= ( ˜

D+˜

M)xk

3+ (I−EDE)fk.

Again as in (2.65) it follows that ˜

Mxk

3= 0,

and thus ˜

Nxk+1

3=˜

Dxk

3+ (I−EDE)fk.

27

Finally, for k≥k0one has

˜

Cxk+1

2=˜

CADA

k

X

j=k0

(EDA)k−jEDfj

=˜

C

|{z}

=EEDE

EDADA k−1

X

j=k0

(EDA)k−jfj+fk!

=EDEADA k−1

X

j=k0

(EDA)k−jfj+fk!

=EDEADA EDA

k−1

X

j=k0

(EDA)k−j−1fj+fk!

=AADA

k−1

X

j=k0

(EDA)k−j−1EDfj+ADAEDEfk

=AADAADA

k−1

X

j=k0

(EDA)k−j−1EDfj+ADAEDEfk

=˜

DADA

k−1

X

j=k0

(EDA)k−j−1EDfj+ADAEDEfk

=˜

Dxk

2+ADAEDEfk,

and for k < k0analogously,

˜

Dxk

2=−˜

DEDE

k0

X

j=k+1

(ADE)j−k−1ADfj−1

=−AADAEDEAD

k0

X

j=k+1

(ADE)j−k−1fj−1

=−AADEDE

k0

X

j=k+1

(ADE)j−k−1fj−1

=−AADEDE k0

X

j=k+2

(ADE)j−k−1fj−1+fk!

=−AADEDEfk−ADAADEEDE

k0

X

j=k+2

(ADE)j−k−2fj−1

=−AADEDEfk−ADEEDEEDE

k0

X

j=k+2

(ADE)j−k−2fj−1

28

=−AADEDEfk+EEDE −EDE

k0

X

j=k+2

(ADE)j−k−2ADfj−1!

=−AADEDEfk+˜

Cxk+1

2.

Lemma 2.18 then implies the assertion.

We finally combine Lemmas 2.10, 2.14, 2.19 and Theorem 2.17.

Theorem 2.20. Let E, A ∈Cn,n with EA =AE satisfy (2.42). Also, let νE= ind(E),

νA= ind(A),{fk}k∈Zwith fk∈Cnand k0∈Z. Then the following statements hold.

1. Every solution {xk}k≥k0of

Exk+1 =Axk+fk,k≥k0(2.66)

satisfies

xk=EDAk−k0EDEv +

k−1

X

j=k0EDAk−j−1EDfj

−I−EDEνE−1

X

i=0 ADEiADfk+i

(2.67)

for k≥k0and for some v∈Cn.

2. Every solution {xk}k≤k0of

Exk+1 =Axk+fk,k≤k0−1 (2.68)

satisfies

xk=ADEk0−kADAv + (I−ADA)

νA−1

X

i=0

(EDA)iEDfk−i−1

−

k0

X

j=k+1

(ADE)j−k−1ADfj−1

(2.69)

for k≤k0and for some v∈Cn.

3. Every solution {xk}k∈Zof

Exk+1 =Axk+fk,k∈Z(2.70)

29

satisfies

xk=(I−ADA)

νA−1

X

i=0

(EDA)iEDfk−i−1

+





ADA(EDA)k−k0EDEv +Pk−1

j=k0(EDA)k−j−1EDfjfor k≥k0

EDE(ADE)k0−kADAv −Pk0

j=k+1(ADE)j−k−1ADfj−1for k≤k0

−(I−EDE)

νE−1

X

i=0

(ADE)iADfk+i

(2.71)

for k∈Zand for some v∈Cn.

Proof. Since the problem is linear any solution may be written as a particular solution of

the inhomogeneous problem plus a solution of the homogeneous problem.

Corollary 2.21. Let the assumptions of Theorem 2.20 hold. Then the following statements

hold.

1. The initial value problem consisting of (2.66) and (1.6) possesses a solution if and only

if there exists a v∈Cnwith

ˆx=EDEv −I−EDEνE−1

X

i=0 ADEiADfk0+i.(2.72)

If this is the case, then the solution is unique.

2. The initial value problem consisting of (2.68) and (1.6) possesses a solution if and only

if there exists a v∈Cnwith

ˆx=ADAv + (I−ADA)

νA−1

X

i=0

(EDA)iEDfk0−i−1.(2.73)

If this is the case, then the solution is unique.

3. The problem consisting of (2.70) and (1.6) possesses a solution if and only if there

exists a v∈Cnwith

ˆx= (I−ADA)

νA−1

X

i=0

(EDA)iEDfk0−i−1

+ADAEDEv (2.74)

−(I−EDE)

νE−1

X

i=0

(ADE)iADfk0+i.

If this is the case, then the solution is unique.

30

Remark 2.22. Analogously to Remark 2.15 one might presume that it is meaningful to

require that the initial condition satisfies (2.74), even when only calculating into the future.

Following the approach in [8] we now consider the case that (E, A) is regular but Eand A

do not commute. In this case we need the following Lemma.

Lemma 2.23. [1] Let E, A ∈Cn,n with (E, A)regular. Let ˜

λ∈Cbe chosen such that

˜

λE −Ais nonsingular. Then the matrices

˜

E=˜

λE −A−1E,˜

A=˜

λE −A−1A

commute.

Proof. We have

˜

λ˜

E−˜

A=˜

λ˜

λE −A−1E−˜

λE −A−1A

=˜

λE −A−1˜

λE −A

=I

and thus

˜

A=˜

λ˜

E−I.

Therefore we finally obtain

˜

E˜

A=˜

E˜

λ˜

E−I

=˜

λ˜

E−I˜

E

=˜

A˜

E.

Remark 2.24. [8] Since the factor ˜

λE −A−1represents a simple scaling of the descrip-

tor system, results similar to Theorem 2.20 and Corollary 2.21 hold for the general case

provided that the coefficient matrices form a regular matrix pair. We only need to perform

the replacements

E←˜

λE −A−1E,A←˜

λE −A−1A,f←˜

λE −A−1f,

in Theorem 2.20 and Corollary 2.21.

Also note that condition (2.42) is equivalent to the regularity of a matrix pair (E, A), which

satisfies EA =AE. Thus, the assumptions of Theorem 2.20 and Corollary 2.21 can essen-

tially be reduced to the regularity of the original matrix pair (E, A).

31

2.3 Conclusion

In this chapter we have first examined the behavior of a general linear discrete-time descriptor

system with constant coefficients, i.e., of a system of the form (1.4), with the help of the

Kronecker canonical form. Here we have only considered the case where one has an initial

condition at some point k0∈Zand wants to get a solution for all k≥k0.

Then we concentrated on regular systems, i.e., on systems of the form (1.4) where the matrix

pair (E, A) is regular. For such systems we wrote down the explicit solution with the help

of the Drazin inverse. In contrast to the continuous-time case one can distinguish between

three different cases for such systems. The first case is where one has an initial condition

given at point k0∈Zand only wants to get a solution for indices k≥k0. The second case is

where one has an initial condition given at point k0∈Zand only wants to get a solution for

indices k≤k0. These first two cases are closely related, since the first case can be transferred

into the second one by a variable substitution. The third case is really different from the first

two cases. Here also an initial condition is given at some point k0∈Zbut one is looking for

a solution for indices k≥k0as well as for indices k≤k0. This puts stronger restrictions on

the initial condition, i.e., the set of consistent initial conditions in the third case is smaller

than in the first or second case.

32

Chapter 3

Linear discrete-time descriptor

systems with variable coefficients

In this chapter we consider the more general case, where the matrices Eand Aare allowed

to change with time, i.e., system (1.5) is considered. Thus, we give up one restriction and

expect the theory to become more complicated. This analysis is done analogously as in [8].

Note that there are other approaches to systems of the form (1.5), for example [12].

3.1 Canonical forms

As in the continuous case (see [8]) the unique solvability of (1.5) and the regularity of all

(Ek, Ak) for k∈Kare completely independent in the discrete case. This can be seen by the

following three examples.

Example 3.1. Let K=Nand let the matrix sequences {Ek}k∈K,{Ak}k∈K⊂C2,2be given

by

(E1, A1) = 0 0

0 0,1 0

0 1,

(Ei, Ai) = 1 0

0 1,0 0

0 0for i≥2.

Then (1.5) is equivalent to

x1=−f1,

xi+1 =fifor i≥2.

Thus, there is no equation for x2and it can be chosen arbitrarily.

Example 3.2. From [6] pp.24-25. Let K=Zand let the matrix sequences {Ek}k∈K,

{Ak}k∈K⊂C2,2be given by

(Ek, Ak) = 0 0

−1k,−1k−1

0 0 for k∈K.

33

Since

det(λEk−Ak) = 

1 1 −k

−λ λk =λk −(−λ)(1 −k) = λ,

the matrix pair (Ek, Ak) is regular for all k∈K. By defining a sequence {xk}k∈Kthrough

xk:= ckk−1

1,

where the scalar sequence {ck}may be chosen arbitrarily, one gets

Ekxk+1 =0 0

−1k(k+ 1) −1

1ck+1 = 0,

Akxk=−1k−1

0 0 k−1

1ck= 0.

Thus, the solution is not unique, even if one defines an initial condition like x0= 0.

Example 3.3. From [8] pp.56-57. Let K=Zand let the matrix sequences {Ek}k∈K,

{Ak}k∈K⊂C2,2be given by

(Ek, Ak) = 0 0

1−k,−1k

0 0 for k∈K.

Let fk= [fk

1, fk

2]Tbe partitioned conformably with Ekand Akand let xk

1and xk

2be defined

according to Table 1.1. Then (1.5) implies

0 = −xk

1+kxk

2+fk

1,

xk+1

1−kxk+1

2=fk

2,

which implies that

−xk+1

1+ (k+ 1)xk+1

2=−fk+1

1

⇒ −xk+1

1+kxk+1

2=−xk+1

2−fk+1

1

⇒fk

2=xk+1

2+fk+1

1

⇒xk+1

2=fk

2−fk+1

1

⇒xk+1

1=fk

2+kfk

2−kfk+1

1.

Thus, the solution is uniquely determined by the sequence {fk}k∈Kalthough all (Ek, Ak) are

singular.

Considering an arbitrary discrete-time descriptor system with variable coefficients and an

initial condition we see that

Ekxk+1 =Akxk+fk,k∈Z,x0= ˆx

⇔PkEkQk+1Q−1

k+1xk+1 =PkAkQkQ−1

kxk

|{z}

=:˜xk

+Pkfk,k∈Z,Q−1

0x0=Q0ˆx

⇔PkEkQk+1˜xk+1 =PkAkQk˜xk+Pkfk,k∈Z, ˜x0=Q0ˆx,

as long as all Pkand Qkare invertible. This leads to the following definition.

34

Definition 3.4. Two sequences of matrix pairs {(Ek, Ak)}k∈Kand {(˜

Ek,˜

Ak)}k∈Kwith

Ek, Ak,˜

Ek,˜

Ak∈Cm,n are called globally equivalent (on K) if there exist two pointwise

nonsingular matrix sequences

{Pk}k∈Kwith Pk∈Cm,m,

{Qk}k∈K∪{kf+1}with Qk∈Cn,n,

such that

PkEkQk+1 =˜

Ekand PkAkQk=˜

Ak, (3.1)

for all k∈K. We denote this by {(Ek, Ak)}k∈K∼ {(˜

Ek,˜

Ak)}k∈K.

Lemma 3.5. The relation introduced in Definition 3.4 is an equivalence relation.

Proof. Reflexivity: We have {(Ek, Ak)}k∈K∼ {(Ek, Ak)}k∈Kwith Pk=Imand Qk=Infor

all k∈Kand Qkf+1 =In.

Symmetry: From {(Ek, Ak)} ∼ {(˜

Ek,˜

Ak)}, it follows that

PkEkQk+1 =˜

Ekand PkAkQk=˜

Ak,

with nonsingular Pk,Qk. Hence,

Ek=P−1

k˜

EkQ−1

k+1 and Ak=P−1

k˜

AkQ−1

k.

Transitivity: From {(Ek, Ak)} ∼ {(˜

Ek,˜

Ak)}and {(˜

Ek,˜

Ak)} ∼ {(ˆ

Ek,ˆ

Ak)}it follows that

PkEkQk+1 =˜

Ekand PkAkQk=˜

Ak,

˜

Pk˜

Ek˜

Qk+1 =ˆ

Ekand ˜

Pk˜

Ak˜

Qk=ˆ

Ak,

for all k∈K, where all Qk,˜

Qk,Pk,˜

Pk, are nonsingular. From this one immediately sees

that {(Ek, Ak)} ∼ {(ˆ

Ek,ˆ

Ak)}using the matrix sequences ˜

PkPkand Qk˜

Qk.

Definition 3.6. Two pairs of matrices (E, A), ( ˜

E, ˜

A)∈Cm,n are called locally equivalent

if there exist matrices P∈Cm,m and Q, R ∈Cn,n that are all nonsingular, such that

˜

E=PEQ and ˜

A=PAR. (3.2)

Again, we denote this by (E, A)∼(˜

E, ˜

A).

Lemma 3.7. The relation introduced in Definition 3.6 is an equivalence relation.

Proof. By using Lemma 3.5 in the special case K={1}we immediately obtain the assump-

tion.

For the proof of Theorem 3.9 we recall the notion of the echelon form used here.

35

Lemma 3.8. Let A∈Cm,n be a matrix with rank (A) = r. Then there exist invertible

matrices P∈Cm,m and Q∈Cn,n such that

PAQ =Ir0

0 0(3.3)

is in echelon form.

For convenience, we say in the following that a matrix is a basis of a vector space if this is

valid for its columns. For matrix pairs of block matrices we also use the convention that

corresponding blocks (i.e., blocks in the same block row and block column) have the same

number of rows and columns.

Theorem 3.9. Let E, A ∈Cm,n and introduce the following matrices:

Zbasis of corange (E) = kernel EH,(3.4)

Ybasis of corange (A) = kernel AH.(3.5)

Then, the quantities

rf= rank (E), (rank of E; corresponds to forward direction) (3.6a)

rb= rank (A), (rank of A; corresponds to backward direction) (3.6b)

hf= rank ZHA, (rank of ZHA; forward direction) (3.6c)

hb= rank YHE

=rf+hf−rb, (rank of YHE; backward direction) (3.6d)

c=rb−hf, (common part) (3.6e)

a= min(hf, n −rf), (algebraic part) (3.6f)

s=hf−a, (strangeness) (3.6g)

d=rf−c−s, (differential part) (3.6h)

u=n−rf−a, (undetermined variables) (3.6i)

v=m−rf−hf, (vanishing equations) (3.6j)

are invariant under (3.2), and (E, A)is locally equivalent to the canonical form

(˜

E, ˜

A) =













Is0 0 0 0

0Id0 0 0

0 0 Ic0 0

0 0 0 0 0







,







0 0 0 0 0

0 0 Ic0 0

0 0 0 Ia0

Is0 0 0 0

0 0 0 0 0













.(3.7)

We have that either s= 0,u= 0 or s=u= 0. The quantities (3.6) are called local

characteristics or local invariants of the matrix pair (E, A).

36

Proof. Let (Ei, Ai), i= 1,2, be locally equivalent, i.e., let P, Q, R be invertible matrices of

appropriate size such that

E2=PE1Qand A2=PA1R. (3.8)

Since

rank (E2) = rank (PE1Q) = rank (E1) ,

rank (A2) = rank (PA1R) = rank (A1) ,

it follows that rfand rbare invariant under local equivalence. First note that hfis indepen-

dent of the choice of the basis of Z. To see this let Zand ˜

Zbe two bases of corange (E).

Then there exists a regular matrix MZwith

˜

Z=ZMZ.

Then, from

rank ˜

ZHA= rank MZHZHA= rank ZHA

the statement follows.

Let Z2be a basis of corange (E2), i.e kernel E2H= range (Z2). Then Z1:= PHZ2is a basis

of corange (E1), since Z1has full column rank and

kernel E1H= range (Z1) .

To prove this, note that

x∈kernel E1H

⇒0 = E1Hx=Q−HE2HP−Hx

⇒0 = E2HP−Hx

⇒P−Hx∈kernel E2H= range (Z2)

⇒there exists a zsuch that P−Hx=Z2z

⇒x=PHZ2z=Z1z∈range (Z1) .

Conversely, we have

x∈range (Z1)

⇒there exists a zsuch that x=Z1z=PHZ2z

⇒P−Hx=Z2z∈range (Z2) = kernel E2H

⇒0 = E2HP−Hx=QHE1HPHP−Hx=QHE1Hx

⇒0 = E1Hx

⇒x∈kernel E1H.

37

This implies

rank Z2HA2= rank Z2HPA1Q= rank Z1HA1Q= rank Z1HA1,

which shows that also hfis invariant under local equivalence. By exchanging the roles of A

and Eas well as the roles of Zand Yin the previous argument we also see that hb(defined

as rank YHE) is invariant under local equivalence. Since all other quantities are functions

of those first three invariant quantities all quantities are invariant under local equivalence.

Let Z′be a matrix, such that the composed matrix Z′Zis invertible. Also let X, Y be

matrices, such that Y(Z′)HEX=Irf0is in echelon form (3.3). Note that (Z′)HE

has full row rank, since Z′spans the vector space range (E). Thus, the following equivalence

transformations in the sense of Definition 3.6 can be applied to the matrix pair (E, A) (where

all corresponding blocks of both matrices have the same size).

(E, A)∼(Z′)HE

0,(Z′)HA

ZHA

∼Irf0

0,Y(Z′)HA

ZHA

∼











Is0 0 0

0Id+c0 0

0 0 0 0







,





∗ ∗ ∗ ∗

0 0 Ia0

Is0 0 0

0 0 0 0













. (3.9)

Looking at matrix pair (3.9) one can (because of the equivalence relation (3.2)) permute the

columns of the second matrix without permuting the columns of the first matrix. This shows

that the size of the Iablock can be increased by decreasing the size of the Isblock (as long

as the last block column does not vanish) and, conversely, that the size of the Isblock can

be increased at the expense of the size of the Iablock (as long as the second block column

does not vanish). To get a canonical form we choose the Iablock to be of maximum size.

Hence, it is clear that either the last block column vanishes, i.e., u= 0, or that there is no

Is, i.e., s= 0 (or even s=u= 0), as stated in the assertion. These cases are considered

separately in the following.

If u= 0, then (3.9) reads as

(E, A)∼











Is0 0

0Id+c0

0 0 0







,





∗ ∗ ∗

0 0 Ia

Is0 0

0 0 0













∼











Is0 0

0Id+c0

0 0 0







,





0∗0

0 0 Ia

Is0 0

0 0 0













38

∼













Is0 0 0

0Id0 0

0 0 Ic0

0 0 0 0







,







0 0 0 0

0 0 Ic0

0 0 0 Ia

Is0 0 0

0 0 0 0













.

In this case we have hf=s+a≥a=n−rf, which is consistent with (3.6f).

If s= 0, then (3.9) reads as

(E, A)∼





Id+c0 0

0 0 0

0 0 0

,



∗ ∗ ∗

0Ia0

0 0 0





∼





Id+c0 0

0 0 0

0 0 0

,



∗0∗

0Ia0

0 0 0





∼









Id0 0 0

0Ic0 0

0 0 0 0





,





0 0 0 0

0Ic0 0

0 0 Ia0

0 0 0 0









.

In this case we have hf=a=n−rf−u≤n−rfwhich is also consistent with (3.6f).

Finally, the identity (3.6d) can be derived from the canonical form (3.7). Therefore, let ˜

Z

and ˜

Ybe bases of corange ˜

Eand corange ˜

A, where ˜

Eand ˜

Aare the matrices in the

canonical form (3.7), respectively. Then we have

˜

Y=







Is0 0

0Id0

0 0 0

0 0 Iv







.

From this we see, that

rank ˜

YH˜

E= rank











Is0 0 0 0 0

0Id0 0 0 0

0 0 0 0 0 Iv









Is0 0 0 0

0Id0 0 0

0 0 Ic0 0

0 0 0 0 0













=s+d(3.6h)

=s+rf−c−s

=rf−c(3.6e)

=rf−rb+hf,

which shows (3.6d).

39

Note that the form (3.7) can also be obtained by first reducing (E, A) to Kronecker canoni-

cal form and then applying further local equivalence transformations. Doing so shows that

if s > 0, then there has to be at least one real Kronecker block of the form (2.3), i.e., a

Kronecker block with dimension greater or equal to 2×1. Anyway, smay be zero even when

there is any number of Kronecker blocks of the form (2.3), i.e., the presence of Kronecker

blocks of the form (2.3) is necessary for s > 0 but not sufficient.

Comparing this result to the analogous result from [8] (Theorem 3.7) one notices the addi-

tional ”common” part. This part cannot be eliminated, since with the equivalence relation

(3.2) the matrix Amay not be changed by means of the matrix E.

3.2 Forward global canonical form

Like in the constant rank case we first study the case where one starts at some time point

(here this time point is always k= 0) and calculates into the future, i.e., one tries to

get a solution for k≥0. In order to derive a global canonical form, some constant rank

assumptions are introduced. Milder assumptions are necessary in this case, than in the case

where one wants to get a solution for all k∈Z. Despite the issue that we only want to get

a solution for k≥0 we may also consider linear descriptor systems with equations for all

k∈Z(i.e., systems of the form {(Ek, Ak)}k∈Z), since this simplifies moving to the case where

we want to get a solution for all k∈Z. This is no restriction, since every linear descriptor

system of the form {(Ek, Ak)}k∈N0can be extended to one of the form {(Ek, Ak)}k∈Zby

choosing Ek=E0and Ak=A0for all k < 0.

Note that we use here the term canonical form in a way that differs from the terminology of

abstract algebra.

Lemma 3.10. Consider system (1.5) and introduce the matrix sequence {Zk}k∈Kwhere

Zkis a basis of corange (Ek) = kernel EH

kfor all k∈K.(3.10)

Let

rk

f= rank (Ek),k∈K,(3.11a)

rk

b= rank (Ak),k∈K,(3.11b)

hk

f= rank ZkHAk,k∈K,(3.11c)

be the local characteristics of each matrix pencil (Ek, Ak)with k∈K. Then, these charac-

teristic sequences are invariant under global equivalence (3.1).

Assume further that the two local characteristic sequences

rf≡rk

fand hf≡hk

f(3.12)

are constant for all k∈K. Then the sequence of matrix pairs {(Ek, Ak)}k∈Kis globally

equivalent to the sequence











E(1)

kE(2)

k

0 0

0 0 

,



A(1)

k0

0Ihf

0 0 







k∈K

,(3.13)

40

where all matrices [E(1)

kE(2)

k]have full row rank, i.e., they all are of rank rf.

Proof. The invariance of the local characteristics follows directly from Theorem 3.9. Let

Z′

kbe a basis of range (Ek) for all k∈K.

Then [Z′

kZk] is invertible for all k∈Kand Z′

k

HEkhas full row rank rf. Transforming with

the transpose of this sequence from the left yields

{(Ek, Ak)}k∈K∼Z′

k

HEk

0,Z′

k

HAk

ZkHAkk∈K

∼









E(1)

kE(2)

k

0 0

0 0 

,



A(1)

kA(2)

k

0Ihf

0 0 







k∈K

∼









E(1)

kE(2)

k

0 0

0 0 

,



A(1)

k0

0Ihf

0 0 







k∈K

.

Writing down the equations from (1.5) connected with the form (3.13) one obtains

E(1)

kxk+1

1+E(2)

kxk+1

2=A(1)

kxk

1+fk

1,

0 = xk

2+fk

2,

0 = fk

3.

Assuming that K=N0, this system is equivalent to

E(1)

kxk+1

1=A(1)

kxk

1+˜

fk

1,

0 = xk

2+fk

2,

0 = fk

3.

(where ˜

fk

1=fk

1+E(2)

kfk+1

2) which is connected with the sequence of matrix pairs











E(1)

k0

0 0

0 0

,



A(1)

k0

0Ihf

0 0 









. (3.14)

Since this step is reversible the set of solution sequences is not altered. One also may notice

that the new right hand side ˜

fkcan depend on the right hand side of the former next right

hand side fk+1. Thus, one can view this step as an index reduction.

Analogous to [8] Theorem 3.14, in the following Theorem it is shown, that the so obtained

reduced sequences of matrix pairs are still globally equivalent, if the original sequences have

been.

41

Theorem 3.11. Assume that the sequences of matrix pairs

{(Ek, Ak)}k∈Z=









E(1)

kE(2)

k

0 0

0 0 

,



A(1)

k0

0Ihf

0 0 







k∈Z

and

{(˜

Ek,˜

Ak)}k∈Z=









˜

E(1)

k˜

E(2)

k

0 0

0 0 

,



˜

A(1)

k0

0Ihf

0 0 







k∈Z

are globally equivalent on Zand in the form (3.13). In particular, suppose that (3.12) holds

and that all hE(1)

kE(2)

kiand all h˜

E(1)

k˜

E(2)

kihave full row rank rf. Then the sequences of matrix

pairs {(E(1)

k, A(1)

k)}and {(˜

E(1)

k,˜

A(1)

k)}are also globally equivalent on Z.

Proof. By assumption, there exist two pointwise nonsingular matrix sequences

{Pk}k∈Z⊂Cn,n,

{Qk}k∈Z⊂Cm,m,

such that

PkEk=˜

EkQk+1, (3.15a)

PkAk=˜

AkQk, (3.15b)

for all k∈Z. By partitioning the transforming matrices appropriately we get

˜

AkQk=



˜

A(1)

k0

0Ihf

0 0 

"Q(1,1)

kQ(1,2)

k

Q(2,1)

kQ(2,2)

k#

=



˜

A(1)

kQ(1,1)

k˜

A(1)

kQ(1,2)

k

Q(2,1)

kQ(2,2)

k

0 0 

=

PkAk=





P(1,1)

kP(1,2)

kP(1,3)

k

P(2,1)

kP(2,2)

kP(2,3)

k

P(3,1)

kP(3,2)

kP(3,3)

k









A(1)

k0

0Ihf

0 0 

(3.16)

=





P(1,1)

kA(1)

kP(1,2)

k

P(2,1)

kA(1)

kP(2,2)

k

P(3,1)

kA(1)

kP(3,2)

k





,

and

˜

EkQk+1 =



˜

E(1)

k˜

E(2)

k

0 0

0 0 

"Q(1,1)

k+1 Q(1,2)

k+1

Q(2,1)

k+1 Q(2,2)

k+1 #

42

=



˜

E(1)

kQ(1,1)

k+1 +˜

E(2)

kQ(2,1)

k+1 ˜

E(1)

kQ(1,2)

k+1 +˜

E(2)

kQ(2,2)

k+1

0 0

0 0 

=

PkEk=





P(1,1)

kP(1,2)

kP(1,3)

k

P(2,1)

kP(2,2)

kP(2,3)

k

P(3,1)

kP(3,2)

kP(3,3)

k









E(1)

kE(2)

k

0 0

0 0 

(3.17)

=





P(1,1)

kE(1)

kP(1,1)

kE(2)

k

P(2,1)

kE(1)

kP(2,1)

kE(2)

k

P(3,1)

kE(1)

kP(3,1)

kE(2)

k





.

From (3.16) we obtain that

P(3,2)

k= 0 for all k∈Z.

Let p∈C1,rfbe any row of the matrices P(2,1)

k, P(3,1)

k. Then from (3.17) we get

p[E(1)

kE(2)

k] = 0.

But since all matrices [E(1)

kE(2)

k] have full row rank as stated in Theorem 3.10, it follows

that also p= 0. Thus, we get that

P(2,1)

k= 0, P(3,1)

k= 0 for all k∈Z,

which means that the left transforming matrices take the form

Pk=





P(1,1)

kP(1,2)

kP(1,3)

k

0P(2,2)

kP(2,3)

k

0 0 P(3,3)

k





.

Hence, the diagonal matrices P(1,1)

k,P(2,2)

k,P(3,3)

khave to be nonsingular. Since we also get

from (3.16) that

Q(2,1)

k=P(2,1)

kA(1)

k= 0,

it follows that all matrices Q(1,1)

k,Q(2,2)

kare invertible. With this, from (3.17) and (3.16) we

finally get that

P(1,1)

kE(1)

k=˜

E(1)

kQ(1,1)

k+1 +˜

E(2)

kQ(2,1)

k+1

|{z}

=0

=˜

E(1)

kQ(1,1)

k+1

and

P(1,1)

kA(1)

k=˜

A(1)

kQ(1,1)

k,

which proves the claim by employing {P(1,1)

k}k∈Zand {Q(1,1)

k}k∈Zas transforming matrix

sequences.

43

Corollary 3.12. Under the assumptions of Theorem 3.11 the sequences of matrix pairs











E(1)

k0

0 0

0 0

,



A(1)

k0

0Ihf

0 0 









and 









˜

E(1)

k0

0 0

0 0

,



˜

A(1)

k0

0Ihf

0 0 









are globally equivalent on Z.

Proof. Using the matrices from the proof of Theorem 3.11 one immediately sees that global

equivalence is achieved by employing the sequences of transforming matrices









P(1,1)

k0 0

0I0

0 0 I







and Q(1,1) 0

0I.

Remark 3.13. The preceding results allow for an inductive procedure closely related to

the corresponding procedure for continuous-time systems [8]. For an original sequence of

matrix pairs {(Ek, Ak)}k∈Z=: {(Ek,0, Ak,0)}k∈Zwe define a sequence (of sequences of matrix

pairs) {{(Ek,i, Ak,i)}k∈Z}i∈N0by the following procedure. First we reduce {(Ek,i, Ak,i)}k∈Zby

Lemma 3.10 to the from (3.13) assuming that the local invariants rf=: rf,i and hf=: hf,i

are constant for all matrix pairs on the whole interval Z. Then we reduce the so obtained

sequence of matrix pairs to the form (3.14) which yields the next sequence of matrix pairs

{(Ek,i+1, Ak,i+1)}k∈Z.

This whole iterative process (although derived from [8]) is very similar to Luenberger’s

shuffle algorithm, which is described in [13] for discrete-time discrete descriptor systems

with constant coefficients.

Observe that we have to have the constant rank assumptions of the form (3.12) for every

step of the procedure. The so obtained sequence of local invariants (which are also global

invariants) {(rf,i, hf,i)}i∈N0is characteristic for a given equivalence class of sequences of

matrix pairs, due to Corollary 3.12. Several properties of this sequence are summed up in

the following Lemma.

Lemma 3.14. Let the sequences {(rf,i, hf,i)}i∈N0and {{(Ek,i, Ak,i)}k∈Z}i∈N0be defined as

in Remark 3.13. In particular, let the constant rank assumptions (3.12) hold. Defining the

quantities

hf,−1:= 0,(3.18a)

ai:= hf,i −hf,i−1∀i∈N0,(3.18b)

44

di:= rf,i +hf,i ∀i∈N0,(3.18c)

si:= rf,i −rf,i+1 ∀i∈N0,(3.18d)

w0:= m−rf,0−hf,0,(3.18e)

wi:= si−1−ai∀i∈N0,(3.18f)

vi:= m−rf,i −hf,i ∀i∈N0,(3.18g)

there exist ξ, µ ∈N0so that:

rf,i ≥rf,i+1 ∀i∈N0,(3.19a)

di≥di+1 ∀i∈N0,(3.19b)

sj= 0 ∀j≥µ,(3.19c)

hf,i ≤hf,i+1 ∀i∈N0,(3.19d)

ai≥ai+1 ∀i∈N0,(3.19e)

di≥hf,i ∀i∈N0,(3.19f)

aj= 0 ∀j≥ξ,(3.19g)

vi=v0+w1+...+wi∀i∈N0,(3.19h)

si≤ai∀i∈N0,(3.19i)

ai+1 ≤si∀i∈N0,(3.19j)

si≥si+1 ∀i∈N0,(3.19k)

ai−wi+1 ≤ai−1−wi∀i∈N,(3.19l)

ai, di, si, wi, vi≥0∀i∈N0.(3.19m)

For any µ∈N0with the property (3.19c) we also have

ai+1 = 0,rf,µ =rf,i,dµ=di∀i≥µ.(3.20)

Proof. Let i∈N0be any non-negative integer. Then we know from Lemma 3.10 that

{(Ek,i, Ak,i)}k∈Z∼









E(1)

k,i E(2)

k,i

0 0

0 0 

,



A(1)

k,i 0

0Ihf,i

0 0 







k∈Z

Remark 3.13

⇒ {(Ek,i+1, Ak,i+1)}k∈Z∼









E(1)

k,i 0

0 0

0 0

,



A(1)

k,i 0

0Ihf,i

0 0 







k∈Z

⇒rf,i = rank [E(1)

k,i E(2)

k,i ]≥rank hE(1)

k,i 0i=rf,i+1 ⇒(3.19a).

Since [E(1)

k,i E(2)

k,i ] has full row rank rf,i, we get

dim range E(1)

k,i + dim corange E(1)

k,i =rf,i

45

and independent of this

dim range E(1)

k,i =rf,i+1.

For k∈Zlet Zkbe a basis of corange E(1)

k,i . Then we know that

hf,i+1 =hf,i + rank ZH

kA(1)

k,i ≤hf,i + dim corange E(1)

k,i . (3.21)

Combining these equations yields

rf,i +hf,i = dim range E(1)

k,i + dim corange E(1)

k,i +hf,i

|{z }

(3.21)

≥hf,i+1

≥rf,i+1 +hf,i+1 ⇒(3.19b).

Since the sequence {rf,i}is non-increasing and bounded by zero, it becomes stationary at

some point µ, which implies (3.19c). From (3.21) one can also see (3.19d). (3.19f) follows

from (3.18c).

Now we show via induction that for all i∈N0we have

{(Ek,i, Ak,i)}k∈Z∼



















E(1)

k,i E(2)

k,i 0

0 0 0





,





A(1)

k,i 0 0

0Iai0

0 0 Ihf,i−1

0 0 0

















k∈Z

, (3.22)

where all [E(1)

k,i E(2)

k,i ] have full row rank.

Induction basis: i= 0

Here hf,−1= 0 and because of Lemma 3.10 we get

{(Ek,0, Ak,0)}k∈Z={(Ek, Ak)}k∈Z∼









E(1)

k,0E(2)

k,0

0 0

0 0 

,



A(1)

k,00

0Ia0

0 0 







k∈Z

,

since a0=hf,0.

Induction step: i→i+ 1 with the help of (3.22).

Since (3.22) is a special form of (3.13), one can immediately perform the reduction to (3.14)

by

{(Ek,i+1, Ak,i+1)}k∈Z∼



















E(1)

k,i 0 0

0 0 0





,





A(1)

k,i 0 0

0Iai0

0 0 Ihf,i−1

0 0 0

















k∈Z

∼









E(1)

k,i 0

0 0

0 0

,



A(1)

k,i 0

0Iai+hf,i−1

0 0 







k∈Z

46

∼





















E(1)

k,i+1 0

0 0





,





A(1)

k,i+1 0

ˆ

A(2)

k,i+1 0

0Ihf,i

0 0





















k∈Z

,

where all E(1)

k,i+1 have full row rank rf,i+1. Adapting the indexing we can proceed with

{(Ek,i+1, Ak,i+1)}k∈Z∼





















E(1)

k,i+1 E(2)

k,i+1 0

0 0 0







,





A(1)

k,i+1 A(2)

k,i+1 0

0Iai+1 0

0 0 0

0 0 Ihf,i

0 0 0





















k∈Z

∼



















E(1)

k,i+1 E(2)

k,i+1 0

0 0 0





,





A(1)

k,i+1 0 0

0Iai+1 0

0 0 Ihf,i

0 0 0

















k∈Z

,

which completes the induction. From the form (3.22) we obtain rf,i −rf,i+1 ≤aifor all

i∈N0. From the induction step we can further see that ai+1 = rank ˆ

A(2)

k,i+1≤rf,i −rf,i+1

for all i∈N0, since ˆ

A(2)

k,i+1 only has rf,i −rf,i+1 rows. This shows (3.19e).

Let µbe as in (3.19c). Then (3.22) implies that

{(Ek,µ, Ak,µ)}k∈Z∼



















E(1)

k,µ E(2)

k,µ 0

0 0 0





,





A(1)

k,µ 0 0

0Iaµ0

0 0 Ihf,µ−1

0 0 0

















k∈Z

.

Further reduction steps show that for all j≥1 we have

{(Ek,µ+j, Ak,µ+j)}k∈Z∼









E(1)

k,µ 0

0 0

0 0

,



A(1)

k,µ 0

0Ihf,µ

0 0 







k∈Z

, (3.23)

since rf,µ+j=rf,µ which means that all E(1)

k,µ have full row rank. Thus, applying a reduction

step to (3.23) does not change anything. This shows that hf,µ+j=hf,µ which is equivalent

to aµ+j= 0. From this (3.20) follows.

Note that there exists a positive integer ξfor which the sequence hf,i gets stationary, i.e.,

hf,i =hf∈N0for all i≥ξ, which follows from the boundedness of the sequence (hf,i ≤m)

and because of (3.19d). This implies (3.19g).

(3.19h) follows since on the one hand we have

vi−v0=m−rf,i −hf,i −(m−rf,0−hf,0)

= (rf,0+hf,0)−(rf,i +hf,i),

47

and on the other hand we have

w1+...+wi

(3.18f)

=s0+...+si−1−a1−...−ai

(3.18d)

=rf,0−rf,i −(a1+...+ai)

(3.18b)

=rf,0−rf,i −(hf,i −hf,0).

(3.19i) can again be seen from the form (3.22). For this, note that we have

si

(3.18d)

=rf,i −rf,i+1

= rank hE(1)

k,i E(2)

k,i i−rank E(1)

k,i 

≤ai,

where the last inequality holds, since the E(2)

k,i matrices only have aicolumns as one can see

from the block structure of (3.22). (3.19l) can be obtained in the following way. Let i∈N

be arbitrarily. Then we have that

wi+ai

(3.18f)

=si−1

(3.19i)

≤ai−1≤ai−1+wi+1

⇒ai−wi+1 ≤ai−1−wi,

where we have used that all wiare non-negative, which is shown below.

To prove the non-negativity of all constants defined in (3.18) without using (3.19l), first

observe that ai≥0, since the sequence {hf,i}i≥−1is non-decreasing. All diare non-negative,

since all rf,i and all hf,i are non-negative. All siare non-negative because the sequence {rf,i}

is non-increasing.

We then see that vi=m−diand thus {vi}is a non-decreasing sequence. Since v0=

m−rf,0−hf,0=m−rank (Ek)−rank ZH

kAk≥0 (where Zkis a basis of corange (Ek)) all

vihave to be non-negative.

Also note, that for all i≥1 we have

vi−vi−1=m−rf,i −hf,i −(m−rf,i−1−hf,i−1)

=rf,i−1−rf,i −(hf,i −hf,i−1) (3.24)

=si−1−ai=wi,

which shows the non-negativity of all wi, since {vi}is non-decreasing.

To finally show (3.19j) and (3.19k) note that for all i∈N0we have

si+1

(3.19i)

≤ai+1

(3.19m)

≤ai+1 +wi+1

(3.18f)

=si.

The previous Lemma 3.14 leads to the following Definition.

48

Definition 3.15. Let {(Ek, Ak)}k∈Zbe a sequence of matrix pairs. Let the sequence

{(rf,i, hf,i)}i∈N0(as described in Remark 3.13) be well defined. In particular, let (3.12) hold

for every entry {(Ek,i, Ak,i)}k∈Zof the sequence (of sequences of matrix pairs) in Remark

3.13. Then, with the definitions (3.18) we call

µ= min{i∈N0|si= 0}(3.25)

the strangeness index of the sequence of matrix pairs {(Ek, Ak)}k∈Zand of the associated

descriptor system (1.5). In the case that µ= 0 we call {(Ek, Ak)}k∈Zand (1.5) strangeness-

free.

In the proof of Lemma 3.14 we were able to see that (under some constant rank assumptions)

every sequence of matrix pairs {(Ek,i, Ak,i)}k∈Zis equivalent to a sequence of the form (3.22).

With i=µwe have sµ= 0 and thus, from (3.23) we see that

{(Ek,µ+1, Ak,µ+1)}k∈Z∼









E(1)

k,µ 0

0 0

0 0

,



A(1)

k,µ 0

0Ihf,µ

0 0 







k∈Z

,

with all E(1)

k,µ having full row rank rf,µ. Finally, one can further reduce all the matrices E(1)

k,µ

to echelon form (3.3) [ Irf,µ 0 ] by global equivalence achieving (with adapted indexing)

{(Ek,µ+1, Ak,µ+1)}k∈Z∼









Irf,µ 0 0

0 0 0

0 0 0

,



A(1)

k,µ 0A(2)

k,µ

0Ihf,µ 0

000







k∈Z

, (3.26)

which can be regarded as a canonical from. One notices that in general not only µbut

µ+ 1 reduction steps are necessary to get to the canonical form, although after µreduction

steps a strangeness-free sequence has already been reached. This situation can be avoided

by introducing a further constant rank assumption in every step of the reduction process,

which was described in Remark 3.13 (see [6]).

Further reductions of (3.26) are possible under additional assumptions.

Theorem 3.16. Let {(Ek, Ak)}k∈Z⊂Cm,n be a strangeness-free sequence of matrix pairs in

the form (3.26), i.e., let

{(Ek, Ak)}k∈Z=









Irf0 0

0 0 0

0 0 0

,



A(1)

k0A(2)

k

0Ihf0

000







k∈Z

.(3.27)

Then the following statements hold.

1. Let k0, kf∈Zwith k0≤kf. Set K:= [k0, kf]∩Zand for k∈Klet all A(1)

kbe

invertible. Then we have

{(Ek, Ak)}k∈Z∼









Irf0 0

0 0 0

0 0 0

,



˜

A(1)

k0˜

A(2)

k

0Ihf0

000







k∈Z

,

where ˜

A(1)

k=Irffor k∈K.

49

2. Let all A(1)

khave the same constant rank, i.e., let rank A(1)

k=pfor all k∈Z. Then

we have

{(Ek, Ak)}k∈Z∼



















Ip0 0 0

0Irf−p0 0

0 0 0 0





,





˜

A(1)

k˜

A(2)

k0˜

A(3)

k

0 0 0 ˜

A(4)

k

0 0 Ihf0

0 0 0 0

















k∈Z

,

where all [˜

A(1)

k,˜

A(2)

k]have full row rank pfor all k∈Z.

3. For all k∈Zlet rank A(2)

k=rf. Then we have

{(Ek, Ak)}k∈Z∼









Irf0 0

0 0 0

0 0 0

,



0 0 A(2)

k

0Ihf0

0 0 0 







k∈Z

.

Proof. In all cases the global equivalence transformations are only applied to the first block

row and the first and third block column. Thus, it is sufficient to look at the sequence of

matrix pairs nIrf0,hA(1)

kA(2)

kiok∈Z. (3.28)

1. Since the identity matrix Irfin (3.28) has to be kept, the allowed global equivalence

transformations are limited. Effectively, the transformations on A(1)

kthat are allowed

are given by

A(1)

k∼PkA(1)

kP−1

k−1. (3.29)

In the following we show by induction that a finite sequence of invertible matrices can

be transformed to identity matrices by the equivalence relation (3.29), i.e., we show

that for ki≤kfthere exist invertible matrices Pk∈Cn,n, for k0−1≤k≤kiwith

PkA(1)

kP−1

k−1=Irffor all k0≤k≤ki. (3.30)

Induction basis: ki=k0

Choose Pk0the inverse of A(1)

k0and Pk0−1=Irf.

Induction step: ki→ki+ 1 ≤kf

Because of (3.30) there exist invertible matrices Pk0−1,...,Pkiwith

Pk0A(1)

k0P−1

k0−1=Irf,

.

PkiA(1)

kiP−1

ki−1=Irf.

50

By setting Pki+1 =Pki(A(1)

ki+1)−1we obtain the equation

Pki+1A(1)

ki+1P−1

ki=Irf.

Thus, by induction the assertion follows.

2. First reduce all A(1)

kto echelon form (3.3) PkA(1)

kQk. Then we have

PkA(1)

k=ˆ

A(1)

k

0

with ˆ

A(1)

khaving full row rank p, since all Qkwere invertible. From this we get with

global equivalence

Irf0,hA(1)

kA(2)

ki∼Pk0,hPkA(1)

kPkA(2)

ki

∼PkPk−10,hPkA(1)

kPk−1−1PkA(2)

ki

= Ip0 0

0Irf−p0,"˜

A(1)

k˜

A(2)

k˜

A(3)

k

0 0 ˜

A(4)

k#!,

with ˆ

A(1)

kP−1

k−1= [ ˜

A(1)

k,˜

A(2)

k] having full row rank pfor all k∈Z.

3. Without loss of generality we may assume that for all k∈Zwe have A(2)

k= [ ˜

A(2)

k,˜

A(3)

k]

with ˜

A(2)

kbeing invertible (otherwise permute columns, which is a global equivalence

transformation). Hence it follows that

(Ek, Ak) = 





Irf0 0 0

0 0 0 0

0 0 0 0

,



A(1)

k0˜

A(2)

k˜

A(3)

k

0Ihf0 0

0 0 0 0 





∼







Irf0 0 0

0 0 0 0

0 0 0 0

,



A(1)

k0˜

A(2)

k˜

A(3)

k

0Ihf0 0

0 0 0 0 







Irf0 0 0

0Ihf0 0

−(˜

A(2)

k)−1˜

A(1)

k0Irf0

0 0 0 In−rf−a−rf











=





Irf0 0

0 0 0

0 0 0

,



0 0 A(2)

k

0Ihf0

0 0 0 





and thus the assertion is proved.

Note that in the proof of part 1 of Theorem 3.16, Pki+1 may be changed, while all other

elements of the matrix sequence {Pk}are not altered. Thus, part 1 of Theorem 3.16 can be

extended to the case where either k0=−∞ or kf=∞.

51

Example 3.17. Consider the constant sequence of matrix pairs consisting of a regular

matrix pencil in Kronecker canonical form

{(Ek, Ak)}k∈Z:= 









0 1 0

0 0 1

0 0 0

,



1 0 0

0 1 0

0 0 1







k∈Z

. (3.31)

This pencil only has one Kronecker block of size 3 corresponding to the infinite eigenvalue

with nilpotency index 3. Using Definition 2.4 pencil (3.31) has Kronecker index 3. To

determine the strangeness index we follow the procedure leading to Definition 3.15. Using

the sign red

∼to denote reducing a sequence of matrix pairs from (3.13) to (3.14) we have that











0 1 0

0 0 1

0 0 0

,



1 0 0

0 1 0

0 0 1







k∈Z

red

∼









010

000

000

,



100

010

001







k∈Z

red

∼









0 0 0

0 0 0

,



1 0 0

0 1 0

0 0 1







k∈Z

.

The sequence of characteristic values (along with some of those defined in (3.18)) is thereby

given as

(rf,0, hf,0, a0, s0) = (2,1,1,1),

(rf,1, hf,1, a1, s1) = (1,2,1,1),

(rf,2, hf,2, a2, s2) = (0,3,1,0),

(rf,3, hf,3, a3, s3) = (0,3,0,0),

(rf,4, hf,4, a4, s4) = (0,3,0,0),

.

which shows that the strangeness index of this sequence of matrix pairs is 2. Thus, this

strangeness index shows the same behavior as the one defined in [8] for constant coefficient

matrix pairs.

One also can guess that increasing the order of the nilpotent block will increase the strange-

ness index accordingly.

Example 3.18. To determine the strangeness index of Example 3.2 we observe that

0 0

−1k,−1k−1

0 0 k∈Z

∼−1k

0 0,0 0

−1k−1k∈Z

∼1 0

0 0,0 0

1 0k∈Z

∼0 1

0 0,0 0

0 1k∈Z

red

∼0 0

0 0,0 0

0 1k∈Z

.

52

The sequence of characteristic values is thereby given as

(rf,0, hf,0, a0, s0) = (1,1,1,1),

(rf,1, hf,1, a1, s1) = (0,1,0,0),

(rf,2, hf,2, a2, s2) = (0,1,0,0),

.

showing that the strangeness index of this sequence of matrix pairs is 1.

Example 3.19. To determine the strangeness index of Example 3.3 we observe that

0 0

1−k,−1k

0 0k∈Z

∼1−k

0 0 ,0 0

−1kk∈Z

∼1 0

0 0,0 0

−1 1k∈Z

∼1 0

0 0,0 0

0 1k∈Z

red

∼1 0

0 0,0 0

0 1k∈Z

.

The sequence of characteristic values is thereby given as

(rf,0, hf,0, a0, s0) = (1,1,1,0),

(rf,1, hf,1, a1, s1) = (1,1,0,0),

(rf,2, hf,2, a2, s2) = (1,1,0,0),

.

showing that the strangeness index of this sequence of matrix pairs is 0 although the matrix

pairs in this example are all singular while the solution is uniquely determined by the right

hand side.

Analogously to [8] one can derive a canonical form for sequences of matrix pairs with well

defined strangeness index without the use of shifts. For convenience, we denote unspecified

blocks in a matrix by ∗.

Theorem 3.20. Let the strangeness index µof the sequence of matrix pairs {(Ek, Ak)}k∈Z

as in (3.25) be well defined. Then, with the definitions from (3.18), {(Ek, Ak)}k∈Zis globally

equivalent to a sequence of the form











Irf,µ 0Wk

0 0 Fk

0 0 Gk

,



∗ ∗ 0

0 0 0

0 0 Ihf,µ









k∈Z

,(3.32)

53

with

Fk=





0F(µ)

k∗

......

...F(1)

k

0







,Gk=





0G(µ)

k∗

......

...G(1)

k

0







,(3.33)

where all F(i)

kand G(i)

khave sizes wi×ai−1and ai×ai−1, respectively, and all Wk= [ ∗ · · · ∗ ]

are partitioned accordingly. In particular, all F(i)

kand G(i)

ktogether have full row rank, i.e.,

rank "F(i)

k

G(i)

k#!=ai+wi=si−1

(3.19i)

≤ai−1∀k∈Z.(3.34)

Proof. With an inductive argument we show that

{(Ek, Ak)}k∈Z∼

{(˜

Ek,i,˜

Ak,i)}k∈Z:=























E(1,1)

k,i E(1,2)

k,i ∗ · · · ∗

0 0 Fk,i ∗

.

..

.

.......

.

..

.

....Fk,1

0 0 ··· ··· 0

0 0 Gk,i ∗

.

..

.

.......

.

..

.

....Gk,1

0 0 ··· ··· 0







,







A(1)

k,i 0 0 ··· 0

0 0 ··· ··· 0

.

..

.

..

.

..

.

..

.

0 0 ··· ··· 0

0Iai0··· 0

.

. 0 .......

.

..

.

.......0

0 0 ··· 0Ia0























,

(3.35)

where:

1. All [ E(1,1)

k,i , E(1,2)

k,i ] are of full row rank rf,i and rank E(1,1)

k,i =rf,i+1.

2. With Zk,i being bases of corange E(1,1)

k,i , we have that rank ZH

k,iA(1)

k,i =ai+1 for all

k∈Z.

3. rank Fk,j

Gk,jis full for all j∈ {1,...,i}and for all k∈Z.

Induction basis: i= 0

From Lemma 3.10 we know that

{(Ek, Ak)} ∼ 









E(1,1)

k,0E(1,2)

k,0

0 0

0 0 

,



A(1)

k,00

0Ihf,0

0 0 









∼









E(1,1)

k,0E(1,2)

k,0

0 0

0 0 

,



A(1)

k,00

0 0

0Ia0











,

54

with [ E(1,1

k,0, E(1,2)

k,0] = rf,0from which the first part of 1. immediately follows. If we perform

the step from (3.13) to (3.14), then we find that

rf,1= rank E(1,1)

k,0

holds and thus 1. is shown. To see 2. let Zk,0be bases of corange E(1,1)

k,0for all k∈Z.

Then 



Zk,00 0

0Iw00

0 0 Ia0



is a basis of corange 





E(1,1)

k,00

0 0

0 0



,

whereas the latter matrix is the one in (3.14). Thus,

hf,1= rank 





ZH

k,00 0

0Iw00

0 0 Ia0







A(1)

k,00

0 0

0Ia0





= rank ZH

k,0A(1)

k,0+a0.

This gives

rank ZH

k,0A(1)

k,0=hf,1−a0=hf,1−hf,0=a1.

Induction step: i→i+ 1 with the help of (3.35)

Applying Lemma 3.10 to the sequence of matrix pairs {(E(1,1)

k,i , A(1)

k,i )}k∈Z(in (3.35)) one

obtains that

{(E(1,1)

k,i , A(1)

k,i )}k∈Z∼







E(1,1)

k,i+1 E(1,2)

k,i+1

0 0

0 0 

,



A(1)

k,i+1 0

0 0

0Iai+1









, (3.36)

where [ E(1,1)

k,i+1 , E(1,2)

k,i+1 ] is of full row rank rf,i+1 due to part 1. and 2. of the inductive

assumption. Applying the transformation corresponding to (3.36) to the original sequence

of matrix pairs (3.35) yields that (Ek, Ak) is globally equivalent to























E(1,1)

k,i+1 E(1,2)

k,i+1 E(1,3)

k,i+1 ∗ · · · ∗

0 0 E(2,3)

k,i+1 ∗ ∗

0 0 E(3,3)

k,i+1 ∗ ∗

0 0 0 Fk,i ∗

.

..

.

..

.

.......

.

..

.

..

.

....Fk,1

0 0 0 ··· ··· 0

0 0 0 Gk,i ∗

.

..

.

..

.

.......

.

..

.

..

.

....Gk,1

0 0 0 ··· ··· 0







,







A(1)

k,i+1 0 0 0 ··· 0

0 0 0 0 ··· 0

0Iai+1 0 0 ··· 0

0 0 0 ··· ··· 0

.

..

.

..

.

..

.

..

.

..

.

..

.

0 0 0 ··· ··· 0

0 0 Iai0··· 0

.

. 0 0 .......

.

..

.

..

.

.......0

0 0 0 ··· 0Ia0























55

∼























E(1,1)

k,i+1 E(1,2)

k,i+1 E(1,3)

k,i+1 ∗ · · · ∗

0 0 Fk,i+1 ∗ ∗

0 0 0 Fk,i ∗

.

..

.

..

.

.......

.

..

.

..

.

....Fk,1

0 0 0 ··· ··· 0

0 0 Gk,i+1 ∗ ∗

0 0 0 Gk,i ∗

.

..

.

..

.

.......

.

..

.

..

.

....Gk,1

0 0 0 ··· ··· 0







,







A(1)

k,i+1 0 0 0 ··· 0

0 0 0 0 ··· 0

0 0 0 ··· ··· 0

.

..

.

..

.

..

.

..

.

..

.

..

.

0 0 0 ··· ··· 0

0Iai+1 0 0 ··· 0

0 0 Iai0··· 0

.

. 0 0 .......

.

..

.

..

.

.......0

0 0 0 ··· 0Ia0























,

(3.37)

by defining Fk,i+1 := E(2,3)

k,i+1 and Gk,i+1 := E(3,3)

k,i+1. Due to the nature of global equivalence it

follows that

rf,i = rank hE(1,1)

k,i E(1,2)

k,i i

= rank 









E(1,1)

k,i+1 E(1,2)

k,i+1 E(1,3)

k,i+1

0 0 E(2,3)

k,i+1

0 0 E(3,3)

k,i+1









.

With regard to the fact that all [ E(1,1)

k,i+1 , E(1,2)

k,i+1 ] also have full rank rf,i+1 one observes that

rank "E(2,3)

k,i+1

E(3,3)

k,i+1#!= rank Fk,i+1

Gk,i+1=rf,i −rf,i+1

(3.18d)

=si,

which means that part 3 of the inductive assumption is shown. Further, we notice that

Gk,i+1 is in the same block row as Iai+1 and in the same block column as Iaiwhich means

that all Gk,i+1 are of size ai+1 ×ai. Since "E(2,3)

k,i+1

E(3,3)

k,i+1#has (rf,i −rf,i+1) rows, it follows that all

Fk,i+1 are of size (rf,i −rf,i+1 −ai+1)×ai

(3.18f)

=wi+1 ×hf,i.

Performing i+ 1 reductions from (3.13) to (3.14) for the sequence (3.37) gives

rf,i+1 = rank hE(1,1)

k,i+1 E(1,2)

k,i+1i.

Analogously performing i+ 2 reductions from (3.13) to (3.14) on the sequence of matrix

pairs (3.37) one finds that

rf,i+2 = rank hE(1,1)

k,i+1i.

56

This shows part 1. To prove part 2. let Zk,i+1 be bases of corange E(1,1)

k,i+1. Performing again

i+ 2 reductions from (3.13) to (3.14) on the sequence (3.37) and denoting the so obtained

sequence by {(ˆ

Ek,ˆ

Ak)}k∈Z, we have that

ˆ

Zk:=







Zk,i+1 0··· ··· ··· ··· 0

0Iwi+1

....

.

...........

.

....Iw0

....

.

....Iai+1

....

.

.......0

0··· ··· ··· ··· 0Ia0







are bases of corange ˆ

Ek. Since all ˆ

Akthen only contain the A(1)

k,i+1 and Iajblock entries,

it is clear that

hf,i+2 = rank ˆ

ZH

kˆ

Ak= rank Zk,i+1A(1)

k,i+1+a0+...+ai+1 (3.38)

⇒rank Zk,i+1A(1)

k,i+1=hf,i+2 −ai+1 −...−a0(3.39)

(3.18b)

⇒rank Zk,i+1A(1)

k,i+1=ai+2. (3.40)

This proves (3.35). Selecting i=µin (3.35) yields

{(Ek, Ak)} ∼























E(1,1)

k,µ E(1,2)

k,µ ∗ · · · ∗

0 0 Fk,µ ∗

.

..

.

.......

.

..

.

....Fk,1

0 0 ··· ··· 0

0 0 Gk,µ ∗

.

..

.

.......

.

..

.

....Gk,1

0 0 ··· ··· 0







,







A(1)

k,µ 0 0 ··· 0

0 0 ··· ··· 0

.

..

.

..

.

..

.

..

.

0 0 ··· ··· 0

0Iaµ0··· 0

.

. 0 .......

.

..

.

.......0

0 0 ··· 0Ia0























(3.18b)

=









E(1,1)

k,µ ˜

Wk

0Fk

0Gk



,



∗0

0 0

0Ihf,µ











.

By induction we obtain that all hE(1,1)

k,µ E(1,2)

k,µ ihave full row rank rf,µ and that all E(1,1)

k,µ

have rank rf,µ+1. Further, from (3.20) we know that rf,µ =rf,µ+1. Thus, all E(1,1)

k,µ have to

have full row rank rf,µ and one may reduce each of these matrices to echelon form (3.3) to

obtain the form (3.32).

57

Further reduction of the Fkand Gkblocks in (3.32) are possible.

Corollary 3.21. Let the strangeness index µof the sequence of matrix pairs {(Ek, Ak)}k∈Z

as in (3.25) be well defined. Then, with the definitions from (3.18), {(Ek, Ak)}k∈Zis globally

equivalent to a sequence of the form











Irf,µ 0Wk

0 0 Ck

0 0 Dk

,



∗ ∗ 0

0 0 0

0 0 Ihf,µ









k∈Z

,(3.41)

with

Ck=





0 0 Iwµ0

Iwµ−10

...

Iw10

0 0







,(3.42)

Dk=







0 0 0 D(µ)

k0∗ · · · 0∗

0 0 0 E(µ)

k0∗ · · · 0∗

0 0 0 D(µ−1)

k··· 0∗

0 0 0 E(µ−1)

k··· 0∗

....

.

..

.

0D(1)

k

0E(1)

k

0 0







,(3.43)

where all D(i)

kand E(i)

khave sizes wi+1 ×(ai−1−wi)and (ai−wi+1)×(ai−1−wi), respectively,

and all Wk= [ ∗ · · · ∗ ]are partitioned accordingly. In particular, all D(i)

kand E(i)

ktogether

have full row rank, i.e.,

rank "D(i)

k

E(i)

k#!=wi+1 + (ai−wi+1) = ai∀k∈Z.(3.44)

The diagonal blocks in Dkall are square matrices.

Proof. The only differences between the forms (3.32) and (3.41) are in the (2,3) and (3,3)

blocks. We can transform these blocks without affecting the other structure (the Wkblock

is affected, but this does not matter, since its structure is not used). Thus, it is sufficient

to only consider the sequence of matrix pairs which is built of the blocks mentioned before.

58

From Theorem 3.20 we know that

Fk

Gk,0

Ihf,µ =























0F(µ)

k∗

......

...F(1)

k

0

0G(µ)

k∗

......

...G(1)

k

0







,







0· · · · · · 0

.

..

.

..

.

0· · · · · · 0

Iaµ...

...

Ia0























,

where all F(i)

khave full row rank wi, due to (3.34). Fixing any i∈ {1,...,µ}we see that

all F(i)

kcan be reduced to echelon form by transforming with invertible matrices in block

row µ+ 1 −iand in block column µ+ 2 −i. This indeed affects the matrices G(i)

kand the

identity matrices Iai−1(by a multiplication with an invertible matrix from the right) but

the G(i)

kmatrices still have full row rank, together with the new F(i)

kmatrices (which now

are in echelon form), i.e., (3.34) still holds. Also the identity matrices Iai−1can be restored

by transforming with invertible matrices from the left. This alters the G(i)

konce more but,

again, condition (3.34) is not affected. From this we see that

Fk

Gk,0

Ihf,µ 

∼























0Iwµ0∗ ∗ · · · · · · ∗

Iwµ−10.

.

...∗ ∗

Iw10

0

0G(µ,1)

kG(µ,2)

k∗ ∗ · · · · · · ∗

G(µ−1,1)

kG(µ−1,2)

k

.

...∗ ∗

G(1,1)

kG(1,2)

k

0







,0

Ihf,µ 

















∼























0Iwµ0 0 ∗0∗

Iwµ−10.

.

..

.

...0∗

Iw10

0

0 0 G(µ,2)

k0∗0∗

0G(µ−1,2)

k

.

..

.

...0∗

0G(1,2)

k

0







,0

Ihf,µ 

















59

∼























0Iwµ0 0 0 ··· 0 0

Iwµ−10.

.

..

.

...0 0

Iw10

0

0 0 G(µ,2)

k0∗0∗

0G(µ−1,2)

k

.

..

.

...0∗

0G(1,2)

k

0







,0

Ihf,µ +Bk

















,

where all the matrices Bkare upper triangular and nilpotent. Thus, these Bkmatrices can

again be eliminated by adding multiples of a row kto a row l, where always k > l (i.e.,

by transforming from the left). By splitting the block rows in the lower part (i.e., the part

corresponding to the Gkblock) and using

G(i,2) ="D(i)

k

E(i)

k#,

we finally obtain the assertion, where (3.44) follows from (3.34), since (3.34) now reads

rank 





Iwi0

0D(i)

k

0E(i)

k





=si−1.

The converse of Theorem 3.20 clearly is easier to understand than Theorem 3.20 itself.

Lemma 3.22. Let µ∈N0and let {ˆai}i∈N0be a non-increasing sequence with ˆaµ+1 = 0

and ˆaµ6= 0. Further let {ˆwi}i∈N0be a non-negative sequence with ˆwi+ ˆai≤ˆai−1. Set

ˆ

hf= ˆa0+...+ ˆaµand let ˆrf∈N0be an integer. Assume that Ek, Ak∈Cm,n are matrices

such that the sequence {(Ek, Ak)}k∈Zis globally equivalent to











Iˆrf0Wk

0 0 Fk

0 0 Gk

,



∗ ∗ 0

0 0 0

0 0 Iˆ

hf









k∈Z

,(3.45)

with

Fk=





0F(µ)

k∗

......

...F(1)

k

0







,Gk=





0G(µ)

k∗

......

...G(1)

k

0







,(3.46)

60

where all F(i)

kand G(i)

khave sizes ˆwi×ˆai−1and ˆai×ˆai−1, respectively, and all Wk= [ ∗ · · · ∗ ]

are partitioned accordingly. Also assume that all F(i)

kand G(i)

ktogether have full row rank,

i.e.,

rank "F(i)

k

G(i)

k#!= ˆai+ ˆwi∀k∈Z.(3.47)

Then {(Ek, Ak)}k∈Zhas a well defined strangeness index µas defined in (3.25). The sequence

of characteristic values is thereby given as

rf,i = ˆrf+

µ

X

j=i+1

ˆwj+

µ

X

j=i+1

ˆaj,(3.48)

hf,i =

i

X

j=0

ˆaj.(3.49)

Proof. Performing ireductions from (3.13) to (3.14) on the sequence of pencils (3.46) shows

that for the resulting sequence we have

{(Ek,i, Ak,i)}k∈Z∼























Iˆrf0∗ ∗ · · · ∗ 0

0 0 0 F(µ)

k∗0

.

..

.

..

.

.......0

.

..

.

..

.

....F(i+1)

k0

0 0 0 ··· ··· 0 0

0 0 0 G(µ)

k∗0

.

..

.

..

.

.......0

.

..

.

..

.

....G(i+1)

k0

0 0 0 ··· ··· 0 0







,







∗ ∗ 0· · · · · · · · · 0

0 0 .

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

0 0 0 · · · · · · · · · 0

0 0 Iˆaµ

.

..

.

....

.

..

.

....

0 0 Iˆai

0 0 Iˆ

hf,i−1























,

where ˆ

hf,i−1:= ˆa0+...+ ˆai−1. Thus, one can see, that the rank of all Ek,i is constant in k

and given by (3.48). Once Zkare bases of corange (Ek,i) one also notices, that the rank of

ZH

kAk,i is also constantly equal to ˆai+ˆ

hf,i−1, which shows (3.49).

We derive another Lemma from Theorem 3.20 which will be needed in the next section.

Lemma 3.23. With the assumptions and notation of Theorem 3.20 we have

rank (FiGi+1 . . . Gl) =

µ

X

j=l−i+1

wjand (3.50)

rank (GiGi+1 . . . Gl) =

µ

X

j=l−i+1

aj(3.51)

61

for all l≥i. Also, we know that all non-zero rows of FiGi+1 . . . Gland all non-zeros rows of

GiGi+1 . . . Glare linear independent.

Proof. Since all matrices Fkand Gkare strictly block-upper-triangular with µ+1 block-rows

and block-columns it is clear that all matrices Fkand Gkare nilpotent with nilpotency index

µ+ 1. Also we have that

Fk1Gk2. . . Gkµ+1 = 0 and

Gk1Gk2. . . Gkµ+1 = 0

for any integer sequences k1,...,kµ+1. For l=ithe claim is that all Fihave rank Pµ

j=1 wj

and that all Gihave rank Pµ

j=1 aj, which follows directly from condition (3.34). Thus, it is

sufficient to show the statement of the Lemma for i < l < µ +i.

For i < l < µ +iwe have

Gi+1 . . . Gl=







0··· 0H(µ)

i+1,l ∗

.

.......

.

....H(l−i)

i+1,l

.

. 0

.

..

.

0··· ··· 0







,

where we have used the definition H(j)

p,q := G(j)

p. . . G(j−q+p)

q. Since l−i < µ, the ma-

trix Gi+1 . . . Glmay have two block-rows and two block-columns that are non zero. Pre-

multiplying with Fior Giyields

FiGi+1 . . . Gl=







0··· 0F(µ)

iH(µ−1)

i+1,l ∗

.

.......

.

....F(l−i+1)

iH(l−i)

i+1,l

.

. 0

.

..

.

0··· ··· 0







and

GiGi+1 . . . Gl=







0··· 0H(µ)

i,l ∗

.

.......

.

....H(l−i+1)

i,l

.

. 0

.

..

.

0··· ··· 0







.

62

In the following we show that rank H(j)

i,l =ajand rank F(j)

iH(j−1)

i+1,l =wjand also that

H(j)

i,l and F(j)

iH(j−1)

i+1,l are both of full row rank. This then immediately proves the claim.

We first observe, that H(j)

i,l is a product of a aj×aj−1matrix and a aj−1×aj−2matrix

and so on. Since {ai}i∈N0is a non-increasing sequence, due to (3.19e), this means that H(j)

i,l

is a product of full row rank matrices (G(t)

s), where the column dimension increases with

every further post-multiplication. Thus, H(j)

i,l has to have full row rank aj. With the very

same argument we see that F(j)

iH(j−1)

i+1,l is of full row rank wj, since from (3.34) we know that

wj≤aj−1.

In [6] index concepts for discrete-time descriptor systems are examined. There an approach

similar to the following is taken.

Theorem 3.24. Let {(Ek, Ak)}k∈Zbe a sequence of matrix pairs and introduce the following

matrices:

Zkbasis of corange (Ek) = kernel EkH,

Z′

kbasis of range (Ek),

Tkbasis of kernel (Ek),

T′

kbasis of cokernel (Ek) = range EkH,

Vkbasis of corange ZH

kAkTk−1,

V′

kbasis of range ZH

kAkTk−1,

for all k∈Z. With this, define the following quantities:

rf,k := rank (Ek),

˜ak:= rank ZH

kAkTk−1,

˜sk:= rank VH

kZH

kAkT′

k−1,

for all k∈Z. Under the assumption that

rf≡rf,k,˜a≡˜ak,˜s≡˜sk, k ∈Z,(3.52)

the sequence {(Ek, Ak)}k∈Zis globally equivalent to























I˜s0 0 0

0I˜

d0 0

0 0 0 0







,





0A(1,2)

k0A(1,4)

k

0A(2,2)

k0A(2,4)

k

0 0 I˜a0

I˜s000

0 0 0 0























,(3.53)

where the abbreviation ˜

d:= rf−˜shas been used.

63

Proof. We have

{(Ek, Ak)}k∈Z∼Z′

k

HEkT′

kZ′

k

HEkTk

ZkHEkT′

kZkHEkTk,Z′

k

HAkT′

k−1Z′

k

HAkTk−1

ZkHAkT′

k−1ZkHAkTk−1k∈Z

=Z′

k

HEkT′

k0

0 0,Z′

k

HAkT′

k−1Z′

k

HAkTk−1

ZkHAkT′

k−1ZkHAkTk−1k∈Z

∼









Irf0

0 0

0 0

,



˜

A(1,1)

k˜

A(1,2)

k

V′

k

HZkHAkT′

k−1V′

k

HZkHAkTk−1

VkHZkHAkT′

k−1VkHZkHAkTk−1









k∈Z

=









Irf0

0 0

0 0

,



˜

A(1,1)

k˜

A(1,2)

k

V′

k

HZkHAkT′

k−1V′

k

HZkHAkTk−1

VkHZkHAkT′

k−10







k∈Z

∼









Irf0 0

0 0 0

0 0 0

,



˜

A(1,1)

k˜

A(1,2)

k˜

A(1,3)

k

V′

k

HZkHAkT′

k−1I˜a0

VkHZkHAkT′

k−10 0 







k∈Z

∼









Irf0 0

0 0 0

0 0 0

,



˜

A(1,1)

k0˜

A(1,3)

k

0I˜a0

VkHZkHAkT′

k−10 0 







k∈Z

∼





















I˜s0 0 0

0I˜

d0 0

0 0 0 0







,





˜

A(1,1)

k˜

A(1,2)

k0˜

A(1,4)

k

˜

A(2,1)

k˜

A(2,2)

k0˜

A(2,4)

k

0 0 I˜a0

I˜s000

0 0 0 0





















k∈Z

∼





















I˜s0 0 0

0I˜

d0 0

0 0 0 0







,





0˜

A(1,2)

k0˜

A(1,4)

k

0˜

A(2,2)

k0˜

A(2,4)

k

0 0 I˜a0

I˜s000

0 0 0 0





















k∈Z

.

By further reducing (3.53) to the form























0 0 0 0

0I˜

d0 0

0 0 0 0







,





0A(1,2)

k0A(1,4)

k

0A(2,2)

k0A(2,4)

k

0 0 I˜a0

I˜s000

0 0 0 0























, (3.54)

[6] then defines a sequence of characteristic values which we denote by

{(r(i)

f,˜a(i),˜s(i))}i∈N0(3.55)

64

and a sequence of sequences of matrix pairs similar to Remark 3.13. In the following we

compare the approach from [6] (as in Theorem 3.24) to our approach.

Lemma 3.25. With the assumptions and symbols from the preceding Theorem 3.24 and hf,k

defined as in (3.11c) we have that

hf,k = ˜ak+ ˜sk= ˜a+ ˜s=: hf.

In particular, under the assumptions of Theorem 3.24, assumption (3.12) is satisfied (and

thus all assumptions of Lemma 3.10).

Proof. Using the notation of Theorem 3.24, we see that Vkis a basis of kernel TH

k−1AH

kZk,

from which it is clear that TH

k−1AH

kZkVk= 0. Transposing the very last identity yields

VH

kZH

kAkTk−1= 0. (3.56)

Also we know that

V′

k

HZH

kAkTk−1has full row rank ˜ak, (3.57)

since V′

kis a basis of range ZH

kAkTk−1. Hence, we finally obtain

hf,k = rank ZH

kAk

= rank V′

k

HZH

kAkTk−1V′

k

HZH

kAkT′

k−1

VH

kZH

kAkTk−1VH

kZH

kAkT′

k−1

(3.56)

= rank V′

k

HZH

kAkTk−1V′

k

HZH

kAkT′

k−1

0VH

kZH

kAkT′

k−1

(3.57)

= rank V′

k

HZH

kAkTk−10

0VH

kZH

kAkT′

k−1

= rank V′

k

HZH

kAkTk−1+ rank VH

kZH

kAkT′

k−1

= ˜ak+ ˜sk.

The result then follows, since rf,k, ˜ak, and ˜skare constant for all k∈Z.

Remark 3.26. Note that (3.53) can be transformed to a sequence of matrix pairs of the

form (3.13) by simple block column permutations. From this we see that reducing (3.13)

to (3.14) is quite the same as reducing (3.53) to (3.54). Hence, Lemma 3.25 suggests that

our approach is more general than the one in [6], since it only requires hf,k = ˜ak+ ˜skto be

constant in every step of the reduction procedure.

Remark 3.27. Suppose that a sequence of matrix pairs {(Ek, Ak)}k∈Zsatisfies the assump-

tions of Theorem 3.24. Further suppose, that the associated reduced sequence (3.54) also

satisfies the constant rank assumptions of Theorem (3.24). Then we know that the first two

elements of the sequence {(rf,i)}i∈N0(as defined in Remark 3.13) are given by rf,0=˜

d+ ˜s

and rf,1=˜

d. Thus, from (3.18d) it follows that s0=rf,0−rf,1= ˜s= ˜s(0), where ˜s(0) is as

in (3.55). We conjecture that it is possible to show that for the sequence of characteristic

values (3.55) we have ˜s(i)=sifor all i∈N0, where siis as defined in (3.18d).

65

3.2.1 Local and global invariants

Analogously to [8] in this subsection the connections between the global invariants of a given

sequence of matrix pairs and the local invariants of corresponding inflated descriptor systems

are investigated. The inflated descriptor system of a sequence of matrix pairs {(Ek, Ak)}k∈Z

of stage l∈Nis given by

M(l)

k(z(l))k=N(l)

k(z(l))k+ (g(l))k, (3.58)

where

M(l)

k=





Ek

−Ak+1 Ek+1

......

−Ak+lEk+l







,

N(l)

k=





Ak0··· 0

0 0 ··· 0

.

..

.

..

.

0 0 ··· 0







,

(z(l))k=





xk

.

xk+l



,

(g(l))k=





fk

.

fk+l



.

This corresponds to combining lsteps of the original descriptor system in one step. For

every l∈Nand every k∈Z, we can determine the local characteristic values (3.6) of the

matrix pair (M(l)

k, N(l)

k). The next Theorem shows that those local characteristic values are

invariant under global equivalence transformations of the original sequence of matrix pairs.

Theorem 3.28. Let {(Ek, Ak)}k∈Zand {(˜

Ek,˜

Ak)}k∈Zbe two globally equivalent sequences

of matrix pairs, i.e., suppose there exist pointwise invertible matrix sequences {Pk}k∈Zand

{Qk}k∈Zwith

Ek=Pk˜

EkQk+1,

Ak=Pk˜

AkQk.

Then for any k∈Zand for any l∈Nthe matrix pairs (M(l)

k, N(l)

k)and (˜

M(l)

k,˜

N(l)

k)(corre-

sponding to (3.58)) are locally equivalent as in Definition (3.6).

66

Proof. Let k∈Zand l∈Nbe fixed. Define the matrices

Π := 





Pk

...

Pk+l





,

Θ := 





Qk+1 ...

Qk+l+1





,

Ψ := 





Qk

...

Qk+l





.

Then we have

Π˜

M(l)

kΘ = 





Pk˜

Ek

−Pk+1 ˜

Ak+1 Pk+1 ˜

Ek+1

......

−Pk+l˜

Ak+lPk+l˜

Ek+l







Θ

=





Pk˜

EkQk+1

−Pk+1 ˜

Ak+1Qk+1 Pk+1 ˜

Ek+1Qk+2

......

−Pk+l˜

Ak+lQk+lPk+l˜

Ek+lQk+l+1







=





Ek

−Ak+1 Ek+1

......

−Ak+lEk+l







=M(l)

k.

For the other matrix we also get

Π˜

N(l)

kΨ = 





Pk˜

AkQk0··· 0

0 0 ··· 0

.

..

.

..

.

0 0 ··· 0







=N(l)

k

from which the local equivalence follows.

To make a statement about the relation between the global characteristics of the original

sequence of matrix pairs and the local characteristics of its inflated descriptor systems as in

[8] we need the following Lemma.

67

Lemma 3.29. With the assumptions and notation of Theorem 3.20 we have

rank 











Gk. . . Gk+l

FkGk+1 . . . Gk+l

.

Fk+l−1Gk+l

Fk+l













=

µ

X

j=l+1

aj+

l+1

X

i=1

µ

X

j=i

wj,(3.59)

for all k∈Zand all l∈N0.

Proof. We use induction to prove the claim.

Induction basis: l= 0

Here (3.59) reads

rank Gk

Fk=

µ

X

j=1

aj+

µ

X

j=1

wj,

which follows directly from Theorem 3.20.

Induction step: l→l+ 1 with the help of (3.59).

Since all non-zero rows of Gk+l+1

Fk+l+1 are linear independent we know, that

rank 











Gk. . . Gk+l+1

FkGk+1 ...Gk+l+1

.

Fk+lGk+l+1

Fk+l+1













= rank



















Gk...Gk+l

FkGk+1 . . . Gk+l

.

Fk+l−1Gk+l

Fk+l







0

0I







Gk+l+1

Fk+l+1 





= rank 











Gk. . . Gk+l

FkGk+1 . . . Gk+l

.

Fk+l−1Gk+l

Fk+l







Gk+l+1







+ rank (Fk+l+1) . (3.60)

To determine the rank of







Gk. . . Gk+l

FkGk+1 . . . Gk+l

.

Fk+l−1Gk+l

Fk+l







|{z }

=:H1

Gk+l+1

|{z }

=:H2

we first notice the special structure of the block entries of the matrices H1and H2, which

was presented in the proof of Lemma 3.23. Lemma 3.23 also shows, that all non-zero rows

68

of H1are linear independent, since due to the inductive assumption we have

rank (H1) =rank (Gk. . . Gk+l) +

rank (FkGk+1 . . . Gk+l) + ...+

rank (Fk+l−1Gk+l) +

rank (Fk+l) .

By multiplying H1with H2from the right more rows that only contain zeros are generated.

However, the non-zero rows of H1H2are still linear independent. To see this, note that

range (GkGk+1 . . . Gk+lGk+l+1)⊂range (GkGk+1 . . . Gk+l) ,

range (FiGi+1 . . . Gk+lGk+l+1)⊂range (FiGi+1 . . . Gk+l) , for all i=k,...,l,

and remember the special structure presented in the proof of Lemma 3.23. This proves that

rank (H1H2) =rank (Gk. . . Gk+lGk+l+1) +

rank (FkGk+1 . . . Gk+lGk+l+1) + ...+

rank (Fk+l−1Gk+lGk+l+1) +

rank (Fk+lGk+l+1)

=

µ

X

j=l+2

aj+

µ

X

j=l+2

wj+...+

µ

X

j=3

wj+

µ

X

j=2

wj. (3.61)

Thus, we have

rank 











Gk...Gk+l+1

FkGk+1 . . . Gk+l+1

.

Fk+lGk+l+1

Fk+l+1













(3.60)

= rank (H1H2) + rank (Fk+l+1)

(3.61)

=

µ

X

j=l+2

aj+

µ

X

j=l+2

wj+...+

µ

X

j=3

wj+

µ

X

j=2

wj+

µ

X

j=1

wj,

which proves the claim.

Lemma 3.29 now will be used to show the following Theorem.

Theorem 3.30. With Ek, Ak∈Cm,n let the strangeness index µof {(Ek, Ak)}k∈Zas in

(3.25) be well defined with global characteristic values (rf,i, hf,i), i ∈N0. Fix any ˆ

k∈Zand

for any l∈N0, let (M(l)

ˆ

k, N(l)

ˆ

k)be the inflated descriptor system corresponding to (3.58) with

local characteristic values (˜rf,l,˜

hf,l)as in (3.6). Then,

˜rf,l = (l+ 1)m−

l

X

i=0

ai−

l

X

i=0

vi

(3.18)

=

l

X

i=0

rf,i +

l−1

X

i=0

hf,i,(3.62)

˜

hf,l =

l

X

i=0

ai

(3.18b)

=hf,l,(3.63)

with videfined in (3.18g).

69

Proof. By Theorem 3.28, we may assume that the sequence {(Ek, Ak)}k∈Zis already in

canonical form (3.32). With this and by using k=ˆ

k, from (3.58) we get that

(M(l), N(l)) := (M(l)

ˆ

k, N(l)

ˆ

k)

=











Ek

−Ak+1 Ek+1

......

−Ak+lEk+l







,





Ak0··· 0

0 0 ··· 0

.

..

.

.....

.

0 0 ··· 0













=













Irf,µ 0Wk

0 0 Fk

0 0 Gk

∗ ∗ 0Irf,µ 0Wk+1

0 0 0 0 0 Fk+1

0 0 −Ihf,µ 0 0 Gk+1

......

∗ ∗ 0Irf,µ 0Wk+l

0 0 0 0 0 Fk+l

0 0 −Ihf,µ 0 0 Gk+l







, N(l)







∼













Irf,µ 0 0

0 0 Fk

0 0 Gk

0 0 0 Irf,µ 0 0

0 0 0 0 0 Fk+1

0 0 −Ihf,µ 0 0 Gk+1

......

0 0 0 Irf,µ 0 0

0 0 0 0 0 Fk+l

0 0 −Ihf,µ 0 0 Gk+l







, N(l)







,

where we have to start eliminating the ∗and Wkblocks from the bottom of the matrix.

Again starting from the bottom, we eliminate the Fk+l−1,...,Fkand Gk+l−1,...,Gkblocks

70

with help of the −Ihf,µ blocks, which are below these blocks. By doing so we obtain that

(M(l), N(l)) is locally equivalent to













Irf,µ 0 0

0 0 Fk

0 0 Gk

0 0 0 Irf,µ 0 0

0 0 0 0 0 Fk+1

0 0 −Ihf,µ 0 0 Gk+1

0 0 0 ...

0 0 −Ihf,µ

...

...Irf,µ 0 0

...0 0 Fk+l−1

...0 0 Gk+l−1

0 0 0 Irf,µ 0 0

0 0 0 0 0 Fk+l

0 0 −Ihf,µ 0 0 Gk+l







, N(l)







∼













Irf,µ 0 0

0 0 Fk

0 0 Gk

0 0 0 Irf,µ 0 0

0 0 0 0 0 Fk+1

0 0 −Ihf,µ 0 0 Gk+1

0 0 0 ...

0 0 −Ihf,µ

...

...Irf,µ 0 0 0

...0 0 0 Fk+l−1Gk+l

...0 0 0 Gk+l−1Gk+l

0 0 0 Irf,µ 0 0

0 0 0 0 0 Fk+l

0 0 −Ihf,µ 0 0 0







, N(l)







71

∼...∼













Irf,µ 0 0 0

0 0 0 FkR(1,l)

k

0 0 0 GkR(1,l)

k

0 0 0 Irf,µ 0 0 0

0 0 0 0 0 0 Fk+1R(2,l)

k

0 0 −Ihf,µ 0 0 0 0

0 0 0 ...

0 0 0 ....

.

0 0 −Ihf,µ

...

...Irf,µ 0 0 0

...0 0 0 Fk+l−1R(l,l)

k

...0 0 0 0

0 0 0 Irf,µ 0 0

0 0 0 0 0 Fk+l

0 0 −Ihf,µ 0 0 0







, N(l)







,

(3.64)

by using the definition R(i,j)

k:= Qj

l=iGk+l. By Lemma 3.29 we can determine the rank of

M(l)as the sum of the ranks of the block columns. Thus, we have

˜rf,l =(l+ 1)rf,µ +lhf,µ + rank 











R(0,l)

k

FkR(1,l)

k

.

Fk+lR(l+1,l)

k













Lemma 3.29

= (l+ 1)rf,µ +lhf,µ +

µ

X

j=l+1

aj+

µ

X

j=l+1

wj+...+

µ

X

j=1

wj. (3.65)

With (3.65) we show (3.62) by induction. For l= 0 the inductive argument follows directly

from (3.18g), since (Eˆ

k, Aˆ

k) = (M(0)

ˆ

k, N(0)

ˆ

k) and a0=hf,0because of (3.18b). Assume that

(3.62) holds for any l < µ. Then we obtain

˜rf,l+1

(3.65)

= (l+ 2)rf,µ + (l+ 1)hf,µ +

µ

X

j=l+2

aj+

µ

X

j=l+2

wj+...+

µ

X

j=1

wj

= (l+ 1)rf,µ +lhf,µ +

µ

X

j=l+1

aj+

µ

X

j=l+2

wj+

µ

X

j=l+1

wj+...+

µ

X

j=1

wj

+rf,µ +hf,µ −al+1

72

(3.65)

= ˜rf,l +rf,µ +hf,µ −al+1 +

µ

X

j=l+2

wj

= ˜rf,l +rf,µ +hf,µ −al+1 +

µ

X

j=1

wj−

l+1

X

j=1

wj

(3.19h)

= ˜rf,l +rf,µ +hf,µ −al+1 + (vµ−v0)−(vl+1 −v0)

= ˜rf,l −al+1 −vl+1 +rf,µ +hf,µ +vµ

|{z }

(3.18g)

=m

,

which, with help of the inductive assumption, yields (3.62).

Let Zbe a basis of corange M(l). We want to determine rank ZHN(l)in order to prove

the statement about ˜

hf,l. Since only the first (three) block rows of N(l)contain non-zero

entries, we only need to consider the first (three) block columns of ZH, i.e., the first (three)

block rows of Z. Let us denote the first (three) block rows by the matrix 



Z1

Z2

Z3

. Then due

to the structure of (3.64) we have Z1= 0. Further, we have that

rank (Z2) =vµ−

µ

X

j=l+1

wj,

rank (Z3) =hf,µ −

µ

X

j=l+1

aj,

because of Lemma 3.23. Therefore, we know that

˜

hf,l = rank ZHN(l)= rank 

0ZH

2ZH

3



∗ ∗ 0

0 0 0

0 0 Ihf,µ







= rank ZH

3=hf,µ −

µ

X

j=l+1

aj

(3.18b)

=hf,µ −(hf,µ −hf,µ−1)−...−(hf,l+1 −hf,l)

=hf,l =

l

X

j=0

aj.

This proves the claim.

Obviously, the formulas (3.62) and (3.63) correspond to the formulas in [8] when one adds

up the algebraic part and the strangeness to one quantity. Also, as in [8], we get the converse

73

of the result, i.e., the knowledge of the sequence {(˜rf,l,˜

hf,l)}allows for the determination of

the sequence {(rf,i, hf,i)}of the global characteristic values of the sequence of matrix pairs

{(Ek, Ak)}k∈K.

Corollary 3.31. Let the strangeness index µof {(Ek, Ak)}k∈Zas in (3.25) be well defined

and let (˜rf,l,˜

hf,l),l= 0, . . . , µ be the local characteristic values of (M(l)

ˆ

k, N(l)

ˆ

k)for some ˆ

k∈Z

(as in Theorem 3.30). Then the sequence {(rf,i, hf,i)}i∈N0of the global characteristic values

of the sequence of matrix pairs {(Ek, Ak)}k∈Zcan be obtained from

hf,l =˜

hf,l,(3.66a)

rf,l = ˜rf,l −˜rf,l−1−˜

hf,l−1,(3.66b)

where we have used ˜

hf,−1= 0 and ˜rf,−1= 0.

Proof. Using the formulas from Theorem 3.30 we immediately get (3.66a). Also, from The-

orem 3.30 we see that

˜rf,l −˜rf,l−1=

l

X

i=0

rf,i +

l−1

X

i=0

hf,i −

l−1

X

i=0

rf,i −

l−2

X

i=0

hf,i =rf,l +hf,l−1,

which proves (3.66b).

3.2.2 Existence and uniqueness of solutions

Concerning existence and uniqueness of sequences of matrix pairs with well defined strange-

ness index we get similar results as in [8].

Theorem 3.32. Let the strangeness index µof the sequence {(Ek, Ak)}k∈Zas in (3.25) be

well defined. Then the discrete descriptor system (1.5) is equivalent (in the sense that there

is a one-to-one correspondence between the solution/sequence spaces) to a discrete descriptor

system of the form

xk+1

1=A(1)

kxk

1+A(3)

kxk

3+fk

1,rf,µ (3.67a)

0 = xk

2+fk

2,hf,µ (3.67b)

0 = fk

3,vµ(3.67c)

where with uµ:= n−rf,µ −hf,µ we have xk

3∈Cuµand each of the inhomogeneities fk

1,fk

2,

fk

3is determined by the original right hand sides fk,...,fk+µ+1 as in (1.5) for all k∈Z.

For the associated forward problem

Ekxk+1 =Akxk+fk, for all k∈N0,(3.68)

we also have the following results:

74

1. (3.68) is solvable if and only if the vµconsistency conditions conditions

fk

3= 0

are fulfilled for all k∈N0.

2. An initial condition x0= ˆxtogether with (3.68) is consistent if and only if in addition

the hf,µ conditions

x0

2= ˆx2=−f0

2

are satisfied.

3. The corresponding initial value problem is uniquely solvable if and only if in addition

uµ= 0

holds.

Proof. Under the assumptions of the Theorem the original sequence {(Ek, Ak)}k∈Zcan be

transformed to the form (3.26) by µ+ 1 reduction steps and proper global equivalence

transformations. Both of these operations generate a one-to-one correspondence of solutions.

3.3 Backward global canonical form

Looking back at Theorem 2.20 we recall that one can look at three different cases of descriptor

systems. The first case is where one starts at a point in time k0, and calculates into the

future, i.e., one calculates a solution {xk}k≥k0. This case has been considered in the previous

section.

The second case is where one starts at a point in time k0, and calculates into the past, i.e.,

one calculates a solution {xk}k≤k0. This case is closely related to the first case. To see this

suppose that a descriptor system of the form

Ekxk+1 =Akxk+fk,k∈Z, (3.69)

xk0= ˆx,

is given and we are looking for a solution for all k≤k0. Then, by defining yk:= x−k+1 and

gk:= f−k, (3.69) is equivalent to

E−kx−k+1 =A−kx−k+f−k,k∈Z

⇔E−kyk=A−kyk+1 −gk,k∈Z

⇔A−kyk+1 =E−kyk+gk,k∈Z.

By calculating the solution of the very last system into the future with the initial condition

y−k0+1 = ˆx, i.e., by calculating {yk}k≥−k0+1, we see through resubstitution, that we got a so-

lution {yk}k≥−k0+1 ={x−k+1}k≥−k0+1 ={xk+1}−k≥−k0+1 ={xk+1}k≤k0−1={xk}k−1≤k0−1=

{xk}k≤k0, i.e., a solution corresponding to the second case of Theorem 2.20. Thus, we do not

have to consider the second case separately. We make the following definition.

75

Definition 3.33. Let {(Ek, Ak)}k∈Zbe a sequence of matrix pairs. Then

{(A−k, E−k)}k∈Z(3.70)

is called the reversed sequence of matrix pairs. Analogously, the descriptor system corre-

sponding to (3.70) is called the reversed descriptor system. Also, the strangeness index of

(3.70) is call reversed strangeness index and is denoted by µb(for backwards). In contrast to

this, the strangeness index of the original sequence is also called ordinary strangeness index

or forward strangeness index and denoted by µf.

Clearly the reversed system is very similar to the original system and obviously the following

Lemma may be derived.

Lemma 3.34. Let {(Ek, Ak)}k∈Zand {(˜

Ek,˜

Ak)}k∈Zbe two globally equivalent sequences of

matrix pairs. Then the reversed sequences are also globally equivalent.

Proof. By assumption we know that there exist invertible matrices Pk, Qksuch that

Ek=Pk˜

EkQk+1,

Ak=Pk˜

AkQk,

for all k∈Z. Substituting kby −kthen yields

E−k=P−k˜

E−kQ−k+1,

A−k=P−k˜

A−kQ−k,

for all k∈Z. Setting Rk:= P−kand Sk:= Q−k+1, this condition then is finally equivalent

to

A−k=Rk˜

A−kSk+1,

E−k=Rk˜

E−kSk,

for all k∈Z, which proves the claim with the transformation matrix sequences {Rk}and

{Sk}.

76

Chapter 4

Two-way global canonical form for

linear descriptor systems

In the previous chapter we generalized the first and second case of Theorem 2.20 to linear

discrete-time descriptor systems with variable coefficients. This leaves the third case of

Theorem 2.20, i.e., the case where one wants to get a solution on the whole Z. This case is

really different from the first two cases and it will therefore be studied in this chapter. To see

that the third case of Theorem 2.20 really is different from the first two cases consider the

form (3.26). The problem is that in this (strangeness-free) form (3.26) the A(1)

kare allowed

to be arbitrary. Consider a descriptor system which only consists of the (1,1) block in (3.26).

For such a system one can easily compute the unique value of xk0+1 once the value of xk0

is given. In contrast, if the value for xk0is given there may be many choices of appropriate

xk0−1values (e.g., xk0= 0xk0−1,xk0−1=xk0−2,xk0−2=xk0−3, ...) or even no possible

choice of an appropriate xk0−1value (e.g., xk0=xk0−1,xk0−1= 0xk0−2, given that xk06= 0),

depending on the sequence of the A(1)

kmatrices. Also, the solvability may vary from iterate

to iterate. This shows that additional rank assumptions are necessary to generalize the third

case of Theorem 2.20 to linear discrete-time descriptor systems with variable coefficients.

One approach that suggests itself is to not only demand the system itself to have well defined

strangeness index but to also demand the reversed system to have well defined strangeness

index.

Lemma 4.1. For k∈Zlet Ek, Ak∈Cm,n be such matrices, that the strangeness index µf

and the reversed strangeness index µbof {(Ek, Ak)}k∈Zare both well defined. Perform one

step of index reduction from (3.13) to (3.14) on the reversed sequence {(A−k, E−k)}k∈Zand

denote the so obtained sequence by {(˜

A−k,˜

E−k)}k∈Z. Then, not only the reversed strange-

ness index ˜µb(i.e., the strangeness index of {(˜

A−k,˜

E−k)}k∈Z) but also the strangeness index

˜µfof {(˜

Ek,˜

Ak)}k∈Zis well defined. We have ˜µf≤µfand ˜µb≤µb.

Proof. Since the strangeness index of {(Ek, Ak)}k∈Zis well defined and all Akhave the

same constant rank (since the strangeness index of {(A−k, E−k)}k∈Zis also well defined) we

know from Corollary 3.21 and from Theorem 3.16 part 2. that with p∈ {1,...,rf,0}, the

definitions from (3.18) and with bi:= ai−wi+1 for i∈N0the sequence {(Ek, Ak)}k∈Zis

77

globally equivalent to a sequence of matrix pairs of the form























Ip0 0 ∗ · · · · · · · · · · · · · · · · · · · · · ∗

0Irf,0−p0∗ · · · · · · · · · · · · · · · · · · · · · ∗

0··· 0 0 0 Iwµf0

.

..

.

.Iwµf−10

.

..

.

....

.

..

.

.Iw10

0··· 0 0

0··· 0 0wµf+1 0 0 D(µf)

k0∗ · · · 0∗

.

..

.

. 0 0bµf0E(µf)

k0∗ · · · 0∗

.

..

.

. 0wµf0 0 D(µf−1)

k··· 0∗

.

..

.

. 0 0bµf−10E(µf−1)

k··· 0∗

.

..

.

.....

.

..

.

..

.

. 0 D(1)

k

.

..

.

. 0 E(1)

k

.

..

.

. 0w10

0··· 0 0 0b0







,







A(1)

kA(2)

kA(3)

k0· · · · · · · · · · · · · · · · · · · · · 0

000 .

.

..

.

0··· 0.

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

0··· 0 0 · · · · · · · · · · · · · · · · · · · · · 0

0··· 0Iwµf+1

.

..

.

.Ibµf

.

..

.

.Iwµf

.

..

.

.Ibµf−1

.

..

.

.Iwµf−1

.

..

.

.Ibµf−2

.

..

.

....

.

..

.

.Iw1

0··· 0Ib0





















k∈Z

,

78

where all hA(1)

kA(2)

kA(3)

kihave full row rank. Then, with Irf,0−pall ∗blocks in the same row

can be eliminated. This may yield ∗blocks in the first row of the left matrices (next to

hA(1)

kA(2)

kA(3)

ki), but these ∗blocks can again be eliminated by the identity matrices below.

Thus, the system is also equivalent to























Ip0 0 ∗ · · · · · · · · · · · · · · · · · · · · · ∗

0Irf,0−p0 0 · · · · · · · · · · · · · · · · · · · · · 0

0··· 0 0 0 Iwµf0

.

..

.

.Iwµf−10

.

..

.

....

.

..

.

.Iw10

0··· 0 0

0··· 0 0wµf+1 0 0 D(µf)

k0∗ · · · 0∗

.

..

.

. 0 0bµf0E(µf)

k0∗ · · · 0∗

.

..

.

. 0wµf0 0 D(µf−1)

k··· 0∗

.

..

.

. 0 0bµf−10E(µf−1)

k··· 0∗

.

..

.

.....

.

..

.

..

.

. 0 D(1)

k

.

..

.

. 0 E(1)

k

.

..

.

. 0w10

0··· 0 0 0b0







,(4.1)







A(1)

kA(2)

kA(3)

k0· · · · · · · · · · · · · · · · · · · · · 0

0 0 0 .

.

..

.

0··· 0.

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

0··· 0 0 · · · · · · · · · · · · · · · · · · · · · 0

0··· 0Iwµf+1

.

..

.

.Ibµf

.

..

.

.Iwµf

.

..

.

.Ibµf−1

.

..

.

.Iwµf−1

.

..

.

.Ibµf−2

.

..

.

....

.

..

.

.Iw1

0··· 0Ib0





















k∈Z

.

79

Performing one reduction step on the reversed of system (4.1) and reversing the reduced

system back shows that {(˜

Ek,˜

Ak)}k∈Zis globally equivalent to























Ip0 0 ∗ · · · · · · · · · · · · · · · · · · · · · ∗

0Irf,0−p0 0 · · · · · · · · · · · · · · · · · · · · · 0

0··· 0 0 0 Iwµf0

.

..

.

.Iwµf−10

.

..

.

....

.

..

.

.Iw10

0··· 0 0

0··· 0 0wµf+1 0 0 D(µf)

k0∗ · · · 0∗

.

..

.

. 0 0bµf0E(µf)

k0∗ · · · 0∗

.

..

.

. 0wµf0 0 D(µf−1)

k··· 0∗

.

..

.

. 0 0bµf−10E(µf−1)

k··· 0∗

.

..

.

.....

.

..

.

..

.

. 0 D(1)

k

.

..

.

. 0 E(1)

k

.

..

.

. 0w10

0··· 0 0 0b0







,







A(1)

k0A(3)

k0· · · · · · · · · · · · · · · · · · · · · 0

0 0 0 .

.

..

.

0··· 0.

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

0··· 0 0 · · · · · · · · · · · · · · · · · · · · · 0

0··· 0Iwµf+1

.

..

.

.Ibµf

.

..

.

. 0wµf

.

..

.

.Ibµf−1

.

..

.

. 0wµf−1

.

..

.

.Ibµf−2

.

..

.

....

.

..

.

. 0w1

0··· 0Ib0





















k∈Z

.

80

Reordering of rows shows that {(˜

Ek,˜

Ak)}k∈Zis globally equivalent to























Ip0 0 ∗ · · · · · · · · · · · · · · · · · · · · · ∗

0Irf,0−p0 0 · · · · · · · · · · · · · · · · · · · · · 0

0··· 0 0 0 Iwµf0

.

..

.

.Iwµf−10

.

..

.

....

.

..

.

.Iw10

0··· 0 0

0··· 0 0wµf0 0 D(µf−1)

k··· 0∗

.

..

.

.....

.

..

.

..

.

. 0 D(1)

k

.

..

.

. 0w10

.

..

.

. 0wµf+1 0 0 D(µf)

k0∗ · · · 0∗

.

..

.

. 0 0bµf0E(µf)

k0∗ · · · 0∗

.

..

.

. 0 0bµf−10E(µf−1)

k··· 0∗

.

..

.

....0E(1)

k

0··· 0 0 0b0







,







A(1)

k0A(3)

k0· · · · · · · · · · · · · · · · · · · · · 0

000 .

.

..

.

0··· 0.

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

0··· 0 0 · · · · · · · · · · · · · · · · · · · · · 0

0··· 0 0wµf

.

..

.

. 0wµf−1

.

..

.

....

.

..

.

. 0w1

.

..

.

.Iwµf+1

.

..

.

.Ibµf

.

..

.

.Ibµf−1

.

..

.

....

0··· 0Ib0





















k∈Z

.

81

Reordering of columns shows that {(˜

Ek,˜

Ak)}k∈Zis globally equivalent to























Ip0∗ · · · ∗ 0∗ · · · · · · · · · · · · ∗

0Irf,0−p0 0 · · · · · · · · · · · · · · · · · · · · · · · ·

Iwµf

...

Iw1

0··· ··· ··· ··· 0

D(µf−1)

k...

D(1)

k

0

0wµf+1 0D(µf)

k∗ · · · ∗

0 0bµfE(µf)

k∗....

.

0bµf−1E(µf−1)

k∗

...E(1)

k

0b0







,







A(1)

k0· · · · · · 0A(3)

k

Iwµf+1

Ibµf

Ibµf−1

Ibµf−2

...

Ib0





















k∈Z

∼









Iˆrf0∗

0 0 ˜

Dk

0 0 ˜

Ek

,



˜

A(1)

kA(3)

k0

0 0 0

0 0 Iˆ

hf









k∈Z

,

by setting ˆ

hf=bµf+...+b0, ˆrf=rf,0+wµf+...+w1and

˜

Dk=





0D(µf)

k0

......

...D(1)

k

0







,

82

˜

Ek=





0bµfE(µf)

k0

......

...E(1)

k

0b0







,

since wµf+1 = 0 (which can be seen from (3.24), (3.18g), (3.18c) and (3.20)).

By setting ˆai:= biand ˆwi:= wi+1 for i∈N0we finally see, that the sequence {ˆai}i∈N0is

non-increasing due to (3.19l). Also it is clear that

ˆwi+ ˆai=wi+1 +bi=wi+1 +ai−wi+1 =ai

(3.34)

≤ai−1−wi=bi= ˆai−1.

It may happen that some ˆaµf=... = ˆa˜µf+1 = 0 in which case the strangeness index ˜µfof

{(˜

Ek,˜

Ak)}k∈Zhas been decreased. Anyway, with (3.44), all assumptions of Lemma 3.22 are

fulfilled, which shows that the strangeness index ˜µfreally is well defined.

That the reversed strangeness index ˜µbis still well defined follows from the fact that we

have performed the reduction step on the reversed system. To understand this, recall that

Definition 3.15 demands the constant rank assumptions (3.12) to hold for every step of the

reduction procedure. Actually performing one reduction step uses only the constant rank

assumptions (3.12) of the first step and the so obtained system will still satisfy the constant

rank assumptions (3.12) in every further reduction step.

The index reduction performed in the previous Lemma 3.22 will be used frequently in the

following, which is why we introduce the following Definition.

Definition 4.2. Let {(Ek, Ak)}k∈Zbe a sequence of matrix pairs. Then performing one

step of index reduction from form (3.13) to (3.14) on the reversed sequence {(A−k, E−k)}k∈Z

and re-reversing the so obtained sequence is called one step of reversed index reduction. In

contrast to this, the index reduction from form (3.13) to (3.14) on the original sequence

{(Ek, Ak)}k∈Zis also called ordinary index reduction or forward index reduction.

Lemma 4.1 shows that under the assumption that both the strangeness index and the reversed

strangeness index are well-defined one can perform ordinary and reversed index reduction

steps at will. One may conjecture, that the order in which one performs the index reduction

steps is of no meaning, as long as one performs the same number of reduction steps, but this

is false as the following example shows.

Example 4.3. Consider the constant sequence of matrix pairs

1

0,0

1k∈Z

. (4.2)

First performing one ordinary step of index reduction on (4.2) yields the sequence

0

0,0

1k∈Z

. (4.3)

83

Performing one step of reversed index reduction on this sequence does not alter the sequence

any more. First performing one reversed step of index reduction on (4.2), however, yields

the sequence

0

1,0

0k∈Z

, (4.4)

which again is not altered anymore by one further step of ordinary index reduction. Com-

paring (4.3) with (4.4) clearly shows that these two sequences are not globally equivalent,

since (corresponding to Lemma 3.10) those matrix pairs do not have the same characteristic

values.

Remark 4.4. Example 4.3 also shows that a step of reversed index reduction may really

change the ordinary index. Therefore note that (4.2) has an ordinary strangeness index of

1. Nonetheless, performing one step of reversed index reduction on (4.2) yields (4.4), which

surely has an ordinary strangeness index of 0.

Let us derive the canonical form under the assumption that both the strangeness index and

the reversed strangeness index are well defined.

Theorem 4.5. For k∈Zlet Ek, Ak∈Cm,n be matrices, such that the strangeness index and

the reversed strangeness index of {(Ek, Ak)}k∈Zare both well defined. Define the matrices

Zkas a basis of corange (Ek)for k∈Z,

Ykas a basis of corange (Ak)for k∈Z.

Then, there exist hf, hb, q ∈N0such that for all k∈Zwe have

hf=rank ZH

kAk, (forward direction)

hb=rank YH

kEk, (backward direction)

q=hf+hb−rank  YH

kEk

ZH

k+1Ak+1.(4.5)

These quantities in (4.5) are invariant under global equivalence and we have

{(Ek, Ak)}k∈Z∼























E(1)

k0 0 E(2)

k

0Ihb−q0 0

0 0 Iq0

0 0 0 0







,







A(1)

kA(2)

k0 0

0 0 0 0

0 0 Iq0

0 0 0 Ihf−q

0 0 0 0























,(4.6)

where for all k∈Zthe matrices hE(1)

kE(2)

kiand hA(1)

kA(2)

kihave full row rank.

84

Proof. That rank ZH

kAkand rank YH

kEkare invariant under global equivalence follows

from Lemma 3.9. That they are constant for all k∈Zfollows from the fact that the

strangeness index and the reversed strangeness index are both well defined.

To show that qis independent of the choice of the bases, let Ykand ˜

Ykbe bases of corange (Ak)

and let Zkand ˜

Zkbe bases of corange (Ek) for all k∈Z. Then for all k∈Zthere exists

invertible matrices MYkand MZksuch that Yk=˜

YkMYkand Zk=˜

ZkMZk. This shows that

rank  YH

kEk

ZH

k+1Ak+1= rank M−H

Yk0

0M−H

Zk+1  YH

kEk

ZH

k+1Ak+1= rank  ˜

YH

kEk

˜

ZH

k+1Ak+1,

and thus, that qis independent of the choice of the bases. To show the invariance under

global equivalence, let {(˜

Ek,˜

Ak)}k∈Zbe globally equivalent to {(Ek, Ak)}k∈Z, i.e., let Qkand

Pkbe invertible matrices, such that for all k∈Zwe have

Ek=Pk˜

EkQk+1,

Ak=Pk˜

AkQk.

Since

0 = YH

kAk=YH

kPk˜

AkQk=PH

kYkH˜

AkQk,

0 = ZH

kEk=ZH

kPk˜

EkQk+1 =PH

kZkH˜

EkQk+1,

it is clear that ˆ

Yk:= PH

kYkis a basis of corange ˜

Akand that ˆ

Zk:= PH

kZkis a basis of

corange ˜

Ek. With this we see that

rank  ˆ

YH

k˜

Ek

ˆ

ZH

k+1 ˜

Ak+1=rank  YH

kPk˜

Ek

ZH

k+1Pk+1 ˜

Ak+1

=rank  YH

kPk˜

EkQk+1

ZH

k+1Pk+1 ˜

Ak+1Qk+1

=rank  YH

kEk

ZH

k+1Ak+1,

which means that qdoes only depend on the equivalence class.

Since the strangeness index of {(Ek, Ak)}k∈Zis well defined, it is clear that the sequence can

be transformed to the form (3.13). Since the reversed strangeness index is also well defined,

we also know that all Akhave constant rank. Thus, in (3.13) all A(1)

kmatrices also have to

have constant rank. Thus, by transforming the first block row of (3.13) from the left we have

that {(Ek, Ak)}k∈Zis equivalent to





















E(1,1)

kE(1,2)

k

E(2,1)

kE(2,2)

k

0 0





,





A(1,1)

k0

0 0

0Ihf

0 0

















k∈Z

, (4.7)

85

with all A(1,1)

khaving full (constant) row rank. Performing one (ordinary) reduction step

from (3.13) to (3.14) on (4.7) yields the sequence





















E(1,1)

k0

E(2,1)

k0

0 0





,





A(1,1)

k0

0 0

0Ihf

0 0

















k∈Z

. (4.8)

Then it follows from Lemma 4.1 that (4.8) still has well-defined reversed strangeness index.

Let ˜

Ykbe bases of corange 









A(1,1)

k0

0 0

0Ihf

0 0









. Then clearly ˜

Yk=





0 0

I0

0 0

0I





, since all A(1,1)

khave

full row rank. Thus, since the reversed strangeness index of (4.8) is well defined, we know

that for every k∈Zthe matrix

˜

YH

k





E(1,1)

k0

E(2,1)

k0

0 0





=E(2,1)

k0

0 0

has to have constant rank, which means that all E(2,1)

khave to have constant rank. Let us

say all E(2,1)

kmatrices have constant rank ˆg. By reducing all E(2,1)

kin (4.7) to echelon form

and adapting the indexing we then see that {(Ek, Ak)}k∈Zis globally equivalent to























E(1,1)

kE(1,2)

kE(1,3)

k

0IˆgE(2,3)

k

0 0 E(3,3)

k

0 0 0







,





A(1,1)

kA(1,2)

k0

0 0 0

0 0 Ihf

0 0 0





















k∈Z

∼





















E(1,1)

k0˜

E(1,3)

k

0IˆgE(2,3)

k

0 0 E(3,3)

k

0 0 0







,





A(1,1)

kA(1,2)

k0

0 0 0

0 0 Ihf

0 0 0





















k∈Z

∼





















E(1,1)

k0˜

E(1,3)

k

0Iˆg0

0 0 E(3,3)

k

0 0 0







,





A(1,1)

kA(1,2)

k∗

0 0 0

0 0 Ihf

0 0 0





















k∈Z

86

∼





















E(1,1)

k0˜

E(1,3)

k

0Iˆg0

0 0 E(3,3)

k

0 0 0







,





A(1,1)

kA(1,2)

k0

0 0 0

0 0 Ihf

0 0 0





















k∈Z

,

where all hA(1,1)

kA(1,2)

kihave full row rank. Also, all 



E(1,1)

k0˜

E(1,3)

k

0Iˆg0

0 0 E(3,3)

k



have full row rank,

since those matrices are equivalent to the matrices "E(1,1)

kE(1,2)

k

E(2,1)

kE(2,2)

k#as in (4.7), which have

full row rank. So all E(3,3)

kalso have full row rank. Reducing all E(3,3)

kto echelon form then

finally shows that {(Ek, Ak)}k∈Zis globally equivalent to























E(1,1)

k0E(1,3)

kE(1,4)

k

0Iˆg0 0

0 0 Iˆq0

0 0 0 0







,







A(1,1)

kA(1,2)

k0 0

0 0 0 0

0 0 Iˆq0

0 0 0 Ihf−ˆq

0 0 0 0





















k∈Z

∼























E(1,1)

k0 0 E(1,4)

k

0Ihb−ˆq0 0

0 0 Iˆq0

0 0 0 0







,







A(1,1)

kA(1,2)

k0 0

0 0 0 0

0 0 Iˆq0

0 0 0 Ihf−ˆq

0 0 0 0





















k∈Z

,

where hb:= ˆg+ ˆqhas been used. Finally, we have q= ˆq, since (as shown above) the quantity

defined in (4.5) is invariant under global equivalence and qcan directly be calculated from

the last sequence of matrix pairs.

From the form (4.6) one may conjecture that it is also possible to show Theorem 4.5 by

defining

q=hf+hb−rank YH

kEk

ZH

kAk, (4.9)

instead of (4.5). This is not the case. If one would do so, qwould not be invariant under

global equivalence any more as shown by the following example.

Example 4.6. Define the (constant) sequence of matrix pairs

{(Ek, Ak)}k∈Z:= 0 1 0

0 0 0,0 0 0

0 1 0k∈Z

,

87

which has hf= 1, hb= 1 and with both (4.9) or (4.5) q= 1. Transforming this sequence

from the right by the sequence {Qk}k∈Zdefined through

Q2k=



100

010

001

and Q2k+1 =



0 1 0

1 0 0

0 0 1

for all k∈Z,

will yield a sequence {(˜

Ek,˜

Ak)}k∈Z={(EkQk+1, AkQk)}k∈Zwhich satisfies

(˜

E2k,˜

A2k) = 1 0 0

0 0 0,0 0 0

0 1 0 and

(˜

E2k+1,˜

A2k+1) = 0 1 0

0 0 0,0 0 0

1 0 0 for all k∈Z.

This sequence would have q= 0 if one would apply definition (4.9).

The same result as in Theorem 4.5 can be obtained under a weaker assumption.

Corollary 4.7. Let {(Ek, Ak)}k∈Zbe a sequence and define the matrices

Zkas a basis of corange (Ek)for k∈Z,

Ykas a basis of corange (Ak)for k∈Z.

Assume that the quantities

rf=rf,k ≡rank (Ek),(4.10a)

hf=hf,k ≡rank ZH

kAk,(4.10b)

hb=hb,k ≡rank YH

kEk,(4.10c)

q=qk≡hf,k +hb,k −rank  YH

kEk

ZH

k+1Ak+1,(4.10d)

(which are invariant under global equivalence as shown in Theorem 4.5) are constant for all

k∈Z. Then we also have the relation (4.6), where for all k∈Zthe matrices hE(1)

kE(2)

ki

and hA(1)

kA(2)

kihave full row rank.

Proof. First we note that under the given assumptions the sequence {(Ek, Ak)}k∈Zis globally

equivalent to a sequence of the form





















E(1,1)

kE(1,2)

k

E(2,1)

kE(2,2)

k

0 0





,





A(1,1)

k0

0 0

0Ihf

0 0

















k∈Z

,

88

with all "E(1,1)

kE(1,2)

k

E(2,1)

kE(2,2)

k#and all A(1,1)

khaving full row rank. Since qis invariant under global

equivalence, it is clear that

rank E(2,1)

kE(2,2)

k

0Ihf= rank E(2,1)

k0

0Ihf

has to be constant for all k∈Z. Thus, also all E(2,1)

khave to have constant rank. The

remainder of the proof can then be carried out analogously to the proof of Theorem 4.5.

Applying one step of ordinary and one step of reversed index reduction to the form (4.6)

yields the form























E(1,1)

k0 0 0

0Ihb−q0 0

0 0 0 0







,







A(1,1)

k0 0 0

0 0 0 0

0 0 Iq0

0 0 0 Ihf−q

0 0 0 0





















k∈Z

. (4.11)

It is clear that first applying one step of reversed and then one step of ordinary index

reduction will yield another form (i.e., the Isblock then stays in the left matrices and is

therefore missing in the right matrices), as in Example 4.3.

Remark 4.8. The preceding results allow for an inductive procedure similar to Remark

3.13. For an original sequence of matrix pairs {(Ek, Ak)}k∈Z=: {(Ek,0, Ak,0)}k∈Zwe define a

sequence (of matrix pair sequences) {{(Ek,i, Ak,i)}k∈Z}i∈Zby the following procedure. First

we reduce {(Ek,i, Ak,i)}k∈Zby Corollary 4.7 to the from (4.6) assuming that the local invari-

ants rf=: rf,i,hf=: hf,i,hb=: hb,i and q=: qiare constant for all matrix pairs on the whole

interval Z. Then we reduce the so obtained sequence of matrix pairs by one step of ordinary

and one step of reversed index reduction to the form (4.11), which yields the next sequence of

matrix pairs {(Ek,i+1, Ak,i+1)}k∈Z. The so obtained sequence of values {(rf,i, hf,i, hb,i, qi)}i∈N0

is characteristic for a given equivalence class of sequences of matrix pairs, due to Corollary

3.12, Corollary 4.7 and Lemma 3.34.

Remark 4.9. Under the assumption that the strangeness index and the reversed strange-

ness index of the sequence {(Ek, Ak)}k∈Zare both well defined, all constant rank assumptions

which are required in Remark 4.8 are satisfied, because of Lemma 4.1 and Theorem 4.5.

To define a strangeness index under the assumptions of Remark 4.8 we need a Lemma similar

to Lemma 3.14.

Lemma 4.10. Let the sequences {(rf,i, hf,i, hb,i, qi)}i∈N0and {{(Ek,i, Ak,i)}k∈Z}i∈N0be defined

as in Remark 4.8. In particular, let the constant rank assumptions (4.10) hold for every step

89

of the reduction process in Remark 4.8. Defining the quantities

rb,i := rf,i −hb,i +hf,i ∀i∈N0,(4.12a)

sE,i := rf,i −rf,i+1 ∀i∈N0,(4.12b)

sA,i := rb,i −rb,i+1 ∀i∈N0,(4.12c)

si:= sE,i +sA,i ∀i∈N0,(4.12d)

there exists a µ∈N0so that

rb,i = rank (Ak,i)∀i∈N0,k∈Z,(4.13a)

rf,i+1 ≤rf,i ∀i∈N0,(4.13b)

rb,i+1 ≤rb,i ∀i∈N0,(4.13c)

sE,i =sA,i =si= 0 ∀i≥µ.(4.13d)

Proof. (4.13a) follows directly from (3.6d). Let i∈N0be any non-negative integer. Then

we know from Corollary 4.7 that

{(Ek,i, Ak,i)}k∈Z∼























E(1)

k,i 0 0 E(2)

k,i

0Ihb,i−qi0 0

0 0 Iqi0

0 0 0 0







,







A(1)

k,i A(2)

k,i 0 0

0 0 0 0

0 0 Iqi0

0 0 0 Ihf,i−qi

0 0 0 0





















k∈Z

⇒ {(Ek,i+1, Ak,i+1)}k∈Z∼























E(1)

k,i 0 0 0

0Ihb,i−qi0 0

0 0 0 0







,







A(1)

k,i 0 0 0

0 0 0 0

0 0 Iqi0

0 0 0 Ihf,i−qi

0 0 0 0





















k∈Z

.

This clearly shows that rank (Ek,i+1)≤rank (Ek,i) and rank (Ak,i+1)≤rank (Ak,i), which

implies (4.13b) and (4.13c). Since we know that both of the sequences {rf,i}i∈N0and {rb,i}i∈N0

are non-increasing and bounded by zero, they have to become stationary a some point µ,

which shows (4.13d).

The previous Lemma 4.10 leads to the following Definition.

Definition 4.11. Let {(Ek, Ak)}k∈Zbe a sequence of matrix pairs. Let the sequence

{(rf,i, hf,i, hb,i, qi)}i∈N0(as described in Remark 4.8) be well defined. In particular, let (4.10)

hold for every entry {(Ek,i, Ak,i)}k∈Zof the sequence (of sequences of matrix pairs) in Remark

4.8. Then, with the definitions (4.12) we call

µ= min{i∈N0|si= 0}(4.14)

90

the two-way strangeness index of the sequence of matrix pairs {(Ek, Ak)}k∈Zand of the

associated descriptor system (1.5). In the case that µ= 0 we call {(Ek, Ak)}k∈Zand (1.5)

two-way strangeness-free.

Since one step of the iterative procedure described in Remark 4.8 involves one step of ordinary

and one step of reversed index reduction it may happen that the two-way strangeness index

is smaller than the ordinary strangeness index, as shown in the following Example.

Example 4.12. Consider the sequence of (constant) matrix pairs

{(Ek, Ak)}k∈Z=









1 0

0 1

0 0

,



0 0

1 0

0 1







k∈Z

.

With Definition 3.15 we get the sequences

(rf,0, hf,0, a0, s0) = (2,1,1,1),

(rf,1, hf,1, a1, s1) = (1,2,1,1),

(rf,2, hf,2, a2, s2) = (0,3,1,0),

(rf,3, hf,3, a3, s3) = (0,3,0,0),

(rf,4, hf,4, a4, s4) = (0,3,0,0),

.

and thus an ordinary strangeness index of 2. With Definition 4.11, however, we face the

reduction process











1 0

0 1

0 0

,



0 0

1 0

0 1







k∈Z

∼









0 1

1 0

0 0

,



1 0

0 0

0 1







k∈Z

reduction

∼









0 0

1 0

0 0

,



0 0

0 1







k∈Z

and thus the sequences

(rf,0, hf,0, hb,0, q0, rb,0, sE,0, sA,0, s0) = (2,1,1,0,2,1,1,2),

(rf,1, hf,1, hb,1, q1, rb,1, sE,1, sA,1, s1) = (1,1,1,0,1,0,0,0),

(rf,2, hf,2, hb,2, q2, rb,2, sE,2, sA,2, s2) = (1,1,1,0,1,0,0,0),

.

which shows that the two-way strangeness index is 1.

On the other hand, the ordinary strangeness index and the two-way strangeness index can

also coincide as shown in the following example.

91

Example 4.13. Consider the sequence of (constant) matrix pairs

{(Ek, Ak)}k∈Z=









0 1 0

0 0 1

0 0 0

,



1 0 0

0 1 0

0 0 1







k∈Z

. (4.15)

With Definition 3.15 we get the sequences

(rf,0, hf,0, a0, s0) = (2,1,1,1),

(rf,1, hf,1, a1, s1) = (1,2,1,1),

(rf,2, hf,2, a2, s2) = (0,3,1,0),

(rf,3, hf,3, a3, s3) = (0,3,0,0),

(rf,4, hf,4, a4, s4) = (0,3,0,0),

.

and thus an ordinary strangeness index of 2 as shown in Example 3.17. With Definition 4.11

however we face the same reduction process. This is due to the fact that the second matrix

in the matrix pair (4.15) has full row rank, and thus every step of reversed index reduction

has no effect. Thus we get the sequences

(rf,0, hf,0, hb,0, q0, rb,0, sE,0, sA,0, s0) = (2,1,0,0,3,1,0,1),

(rf,1, hf,1, hb,1, q1, rb,1, sE,1, sA,1, s1) = (1,2,0,0,3,1,0,1),

(rf,2, hf,2, hb,2, q2, rb,2, sE,2, sA,2, s2) = (0,3,0,0,3,0,0,0),

(rf,3, hf,3, hb,3, q3, rb,3, sE,3, sA,3, s3) = (0,3,0,0,3,0,0,0),

.

which shows that the two-way strangeness index is also 2.

In Example 4.13 we have seen that the ordinary strangeness index can be equal to the two-

way strangeness index. However, in Example 4.12 we have seen that the ordinary strangeness

index can be bigger as the two-way strangeness index. Also, we observe that one step of

two-way index reduction involves one step of ordinary index reduction. Thus, it should be

possible to show the following Conjecture.

Conjecture 4.14. 1. Every sequence of matrix pairs with well defined two-way strange-

ness index has well defined ordinary strangeness index.

2. Let the strangeness index µf, the reversed strangeness index µband the two-way strange-

ness index µall be well defined. Then we have 2µ≥max(µf, µb)≥µ.

92

4.1 Existence and uniqueness of solutions

With the notation of Remark 4.8 we know that for the two-way strangeness index µwe have

{(Ek,µ, Ak,µ)} ∼























E(1)

k,µ 0 0 E(2)

k,µ

0Ihb,µ−qµ0 0

0 0 Iqµ0

0 0 0 0







,







A(1)

k,µ A(2)

k,µ 0 0

0 0 0 0

0 0 Iqµ0

0 0 0 Ihf,µ−qµ

0 0 0 0





















k∈Z

⇒ {(Ek,µ+1, Ak,µ+1)} ∼























E(1)

k,µ 0 0 0

0Ihb,µ−qµ0 0

0 0 0 0







,







A(1)

k,µ 0 0 0

0 0 0 0

0 0 Iqµ0

0 0 0 Ihf,µ−qµ

0 0 0 0





















k∈Z

.

But we also know from the definitions (4.12) that rank (Ek,µ) = rank (Ek,µ+1) from which

we see that qµ= 0 and that E(1)

k,µ is a matrix with full row rank for all k∈Z. From

rank (Ak,µ) = rank (Ak,µ+1), we analogously see that all A(1)

k,µ already have full row rank.

Thus, every sequence with well defined two-way strangeness index can be transformed by µ+1

reduction steps and appropriate global equivalence transformations to a two-way strangeness-

free sequence of the form

{(Ek,µ+1, Ak,µ+1)} ∼ 



















E(1)

k,µ 0 0

0Ihb,µ 0

0 0 0





,





A(1)

k,µ 0 0

0 0 0

0 0 Ihf,µ

0 0 0

















k∈Z

, (4.16)

where all A(1)

k,µ and all E(1)

k,µ have full row rank. By transformations of the (1,1)-block in (4.16)

one can also achieve that

E(1)

k,µ =I0for all k∈N0

and (4.17)

A(1)

k,µ =I0for all k≤ −1.

From (4.16) with (4.17) one can derive a statement similar to Theorem 3.32 for the case

where one wants to get a solution for all k∈Z.

Theorem 4.15. Let the two-way strangeness index µof the sequence {(Ek, Ak)}k∈Zas in

(4.14) be well defined. Then the discrete descriptor system (1.5) is equivalent (in the sense

that there is a one-to-one correspondence between the solution/sequence spaces) to a discrete

93

descriptor system of the form

xk+1

1=A(1)

kxk

1+A(2)

kxk

4+fk

1,k≥0,rf,µ −hb,µ (4.18a)

xk−1

1=E(1)

k−1xk

1+E(2)

k−1xk

4−fk−1

1,k≤0,rb,µ −hf,µ (4.18b)

xk+1

2=fk

2,hb,µ (4.18c)

0 = xk

3+fk

3,hf,µ (4.18d)

0 = fk

4,m−rf,µ −hf,µ (4.18e)

where with uµ:= n−rf,µ −hf,µ we have xk

4∈Cuµand each of the inhomogeneities fk

1,fk

2,

fk

3,fk

4is determined by the original right hand sides fk−µ−1,...,fk,...,fk+µ+1 as in (1.5)

for all k∈Z. For the problem

Ekxk+1 =Akxk+fk,k∈Z,(4.19)

we also have the following results:

1. (4.19) is solvable if and only if the vµ:= m−rf,µ−hf,µ consistency conditions conditions

fk

4= 0

are fulfilled for all k∈Z.

2. An initial condition x0= ˆxtogether with (4.19) is consistent if and only if in addition

the hf,µ +hb,µ conditions

x0

2= ˆx2=f−1

2,

x0

3= ˆx3=−f0

3,

are satisfied.

3. The corresponding initial value problem is uniquely solvable if and only if in addition

uµ= 0

holds.

Proof. The proof can be carried out as the proof of Theorem 3.32 by using (4.16) and

(4.17).

94

Chapter 5

Algorithms for linear discrete-time

descriptor systems

In the following some algorithms are proposed, which may be used to solve linear discrete-

time descriptor systems (as described in the Chapters 2, 3 and 4). Some of the algorithms

are only suited for the constant coefficient case, while others may also be used in the variable

coefficient case.

5.1 A method with the Drazin inverse

Theorem 2.20 can be used to compute the solution of a regular linear discrete-time descriptor

system (E, A)∈Cn,n ×Cn,n with constant coefficients. There are two problems that arise

in this case.

1. Typically Eand Awill not commute. If they do not commute, one first has to deter-

mine an

invertible matrix Psuch that PEPA =PAPE, (5.1)

i.e., that PE and PA commute. One of the following approaches can be chosen to

determine such a matrix P:

(a) Lemma 2.23 can be used to compute such a matrix. The problem with this

approach is to determine a proper ˜

λ(as in Lemma 2.23), such that ˜

λE −Ais

good conditioned. It is a good idea to choose a ˜

λwhich is nowhere close to any

eigenvalue of λE −A, so that ˜

λE −Ais well conditioned.

(b) One may also employ the Kronecker product to determine such a matrix. For this

note that (5.1) is equivalent to solving the Sylvester equation

EPA =APE, (5.2)

for an invertible P. (5.2) can be rewritten as a system of linear equations, Gy = 0,

using the Kronecker product (see [11]). By calculating a singular value decompo-

sition of Gwe can then determine a basis of kernel (G). Each of the basis vectors

95

can then be reshaped to a matrix P, which solves (5.2). However, we still have to

check if that matrix is invertible. Note that a SVD of Gcan be obtained by com-

puting an SVD of Eand Aseparately (see [11]), which makes the computation

considerably faster.

(c) The GUPTRI-Algorithm (see [4]) together with the solution of a generalized

Sylvester equation (to decouple the eigenvalue infinity from the finite eigenval-

ues; compare [19]) could be employed to calculate a matrix Pwith (5.1). To see

this in the real case, let Uand Vbe orthogonal matrices such that

E=VEfEu

0E∞UT,A=VAfAu

0A∞UT,

is in GUPTRI form, i.e., Efis regular and triangular, Afis quasi-triangular,

E∞is strict upper triangular and A∞is regular and triangular. Following the

approach in [19] we can solve the generalized Sylvester equation

EfY−ZE∞=−Eu,

AfY−ZA∞=−Au.

This equation can be solved uniquely and we obtain matrices Yand Zsuch that

with

W:= VI Z

0I,T:= I−Y

0IUT,

we get

λE −A=WλEf−Af0

0λE∞−A∞T.

Thus, we can rewrite equation (5.2) as

W−1ET−1TPW

|{z}

:= ˜

P

W−1AT−1=W−1AT−1TP WW−1ET−1,

which is equivalent to

Ef0

0E∞˜

PAf0

0A∞=Af0

0A∞˜

PEf0

0E∞. (5.3)

For example, choosing

˜

P=E−1

f0

0A−1

∞,

solves equation (5.3), which means that T−1˜

PW−1solves (5.2). Anyway, the

solution is not unique as one can see by assuming that Af= 0 and E∞= 0 (this

means that we only have the infinite and the zero eigenvalue) in which case every

˜

P= diag (P1, P2)

solves (5.3). Thus, the question arises whether there is a solution which is better

conditioned than T−1˜

PW−1.

96

2. To compute the solution with Theorem 2.20 one needs to compute the products ADx

and EDxfor arbitrary xand also the index of Eand Ais needed.

Some effort has been made to numerically compute the Drazin inverse [5, 18, 20, 21].

One could also compute the Schur form of Eor A, respectively and then solve a

Sylvester equation (to decouple the zero eigenvalues from the other non-zero eigen-

values, similar to the approach in [19]) to calculate the Drazin inverse. To see this,

let

E=USnSu

0SzUH,

be the Schur decomposition of a matrix E in such a way that Szonly has zero eigen-

values and Snonly has non-zero eigenvalues. Following the approach in [19] we can

solve the Sylvester equation

ZSz−SnZ=−Su.

This equation can be solved uniquely, since Szand Snhave no common eigenvalues

and we obtain a matrix Zsuch that with

W:= UI−Z

0I,

we get

E=WSn0

0SzW−1.

Computing the Jordan form of the blocks Sn=VnJnV−1

nand Sz=VzJzV−1

zshows

that

E=WVn0

0Vz

|{z }

:=V

Jn0

0Jz

|{z }

:=J

V−1W−1,

which proves that Jis the Jordan canonical form of E. Thus, we know that

ED=WV JDV−1W−1=WV J−1

n0

0 0V−1W−1

=WV V−1

nS−1

nVn0

0 0V−1W−1=WS−1

n0

0 0W−1,

which means that we can compute the Drazin inverse of E without computing the

Jordan canonical form.

Anyway, one should remember, that (like with the ordinary inverse) computing EDx

without explicitly computing EDhas better numerical properties.

97

Looking at Theorem 2.20 we see that for one element of the solution sequence there are many

calculations of the form EDxnecessary, which makes the whole procedure very expensive.

Also this method (in general) involves the transformation with (regular but) non-unitary

matrices and it can only be used for regular matrix pencils.

All this shows that this method will be of no practical use.

5.2 A method employing the GUPTRI-Algorithm

The GUPTRI-Algorithm [4] can also be used directly to compute the solution of a regu-

lar linear discrete-time descriptor system (E, A)∈Cn,n ×Cn,n with constant coefficients.

Once we have computed the GUPTRI form we can determine the solution by backward

substitution. Let us therefore assume that

(E, A) = EfEu

0E∞,AfAu

0A∞, (5.4)

already is in GUPTRI form. Also assume for simplicity that (E, A) has no complex eigen-

values, i.e., Efis regular and triangular, Afis triangular, E∞is strict upper triangular and

A∞is regular and triangular. Let us further denote the entries of Eand Aby [ei,j]i,j=1,...,n

and [ai,j]i,j=1,...,n, respectively, and let rbe the dimension of the Efand Afblocks. Then

we know that ei,j =ai,j = 0 for all i > j. Additionally, we know that ei,i = 0 for all

i=r+ 1,...,n and that ai,i 6= 0 for all i=r+ 1,...,n. The equations of the descriptor

system (1.4) belonging to (5.4) can thus be written as

n

X

j=i

ei,jxk+1

j=

n

X

j=i

ai,jxk

j+fk

i, for i= 1,...,r, (5.5)

n

X

j=i+1

ei,jxk+1

j=

n

X

j=i

ai,jxk

j+fk

i, for i=r+ 1,...,n. (5.6)

To solve the equations (5.6) we start with the last equation (i.e., for i=n). This equation

can be written as 0 = an,nxk

n+fk

n, which allows for the determination of all xk

n, once all right

hand sides are available. More precisely, we have xk

n=−fk

n

an,n . We further observe that the

last but one of the equations (5.6) can be written as

en−1,nxk+1

n=an−1,n−1xk

n−1+an−1,nxk

n+fk

n−1

⇔xk

n−1=1

an−1,n−1en−1,nxk+1

n−an−1,nxk

n−fk

n−1

=1

an−1,n−1an−1,n

fk

n

an,n

−en−1,n

fk+1

n

an,n

−fk

n−1.

Proceeding in this way we see that we can obtain all xk

ifor i=r+ 1,...,n, once all right

hand sides are available. Note, that by doing so we need the right hand side fk+n−r−1to

98

compute xk

r+1. Thus, there may be more right hand sides necessary than demanded by the

index, because the λE∞−A∞block does not reveal the index.

Given an initial condition for the regular part (5.5) we can also compute the forward solution,

since Efis regular. A backward solution may also be computed, if we divide the λEf−Af

block further into a part belonging to the non-zero eigenvalues and a part belonging to the

zero eigenvalues, which is automatically done by the GUPTRI algorithm. The backward

solution can then be found by performing the very same kind of process described above (for

the equations (5.6)) for the zero eigenvalues of λEf−Af.

Finally, note that the GUPTRI form can also be computed for singular pencils. One could

examine if this GUPTRI form also allows for an easy determination of the solution of such

descriptor systems.

This method will be considerably faster than the one introduced in section 5.1 since we

only have to compute the GUPTRI form and perform backward substitution. Nevertheless,

this method has the drawback that there may have to be n−r−1 future right hand sides

available, where n−r−1 is the number of infinite eigenvalues.

5.3 A reduction method

In this section an algorithm is presented, which can also be applied to discrete-time de-

scriptor systems with variable coefficients, as long as the system has a well defined forward

strangeness index and/or a well defined backward strangeness index. The matlab code of

the algorithm can be found in Appendix A and is called solve_dds. We will use the line

numbers on the left side of the algorithm to refer to specific parts of the code and we will use

true-type fonts to reference to parameters, variables and other interna of the code solve_dds

(e.g., Efun means the first parameter of the matlab function in Appendix A).

The code starts in lines 91-101 by first setting up a few variables and then determining if a

forward computation is necessary (result saved in do_forward) and if a backward computa-

tion is necessary (result saved in do_backward).

If the forward computation is necessary, solve_dds first computes the quantities r_f(i+1)

=rf,i and h_f(i+1) =hf,i as in Remark 3.13 and mu_f = max(µf,1), where µfis the forward

strangeness index (in lines 103-128). If the backward computation is necessary solve_dds

also computes the quantities r_b and h_b as in Remark 3.13 for the reversed descriptor

system and mu_b = max(µb,1), where µbis the reversed strangeness index (in lines 130-155).

mu_f and mu_b are then later updated by the real strangeness indices in lines 225 and 269.

In lines 158-184 all constraints on the initial condition are determined and the initial condi-

tion xhat is altered so that new initial condition is the consistent initial condition which is

closest to the old initial condition. As shown in Corollary 2.21, there are more constraints

on the initial condition if one is looking for a two-way solution. This is represented by lines

160-169 of the code. Later, in lines 186-272 the actual forward and backward solutions are

determined. Finally, in lines 274-281 some terminal computations are done.

The actual core of the algorithm is the function advance_inflated_system, which is needed

for the index determination as well as the computation of the solution. In the following we

99

will shortly illustrate the functioning of advance_inflated_system in the forward case.

Note that in every call to advance_inflated_system the matlab cell arrays Ek,Ak and fk

are inserted into the function and are also returned by the function. These cell arrays are

assumed to have µ:=mu (parameter to the function advance_inflated_system) entries on

entry and they have µ+ 1 entries on exit. These cell arrays are interpreted as the entries of

the matrix

D(1) := 





−A1E1−f1

.......

.

−Aµ−1Eµ−1−fµ−1

−AµEµ−fµ







,

where Ai:= Ak{i},Ei:= Ek{i} and fi:= fk{i}. Note that this matrix can be generated

with the help of the function dumparray (in line 415). Clearly, any vector

˜xµ=





x1

.

xµ

xµ+1

1







(5.7)

with D(1) ˜xµ= 0 gives a sequence {xi}µ+1

i=1 , which solves the equations Eixi+1 =Aixi+fifor

i∈ {1,...,µ}. In lines 294-297 matrix D(1) is first extended to the matrix

D(2) := 





−A1E1−f1

.......

.

−AµEµ−fµ

−Aµ+1 Eµ+1 −fµ+1







,

where Aµ+1 := feval(Afun,kact+mu),Eµ+1 := feval(Efun,kact+mu) and fµ+1 := feval

( ffun,kact+mu), since we only consider the forward case here. Then advance_inflated

_system performs a singular value decomposition of Eµ+1 to discover corange (Eµ+1) (in line

301). This yields a matrix Uµ+1 with

Uµ+1Eµ+1 =E(1)

µ+1

0,

where E(1)

µ+1 has full row rank rf,0. Setting D(3) := diag (Imµ, Uµ+1)D(2) shows that

D(2) ˜xµ+1 = 0 ⇔D(3) ˜xµ+1 = 0,

100

where ˜xµ+1 is as in (5.7). D(3) has the form

D(3) =





−A1E1−f1

.......

.

−AµEµ−fµ

−A(1)

µ+1 E(1)

µ+1 −fµ+1

1

−A(2)

µ+1 0−fµ+1

2







.

In line 319 we then perform a singular value decomposition of A(2)

µ+1 =U(2)

µ+1S(2)

µ+1V(2)

µ+1

T.

Setting D(4) := diag Imµ+rf,0, U(2)

µ+1TD(3)diag Inµ, V (2)

µ+1, In+1shows that the solution sets

of D(3) ˜xµ+1 = 0 and D(4) ˜xµ+1 = 0 differ, but there is a one-to-one correspondence between

the sets through a linear mapping, which can be represented by V(2)

µ+1. To keep track of this

linear mapping this transformation is stored in line 335, so that we can compute the original

solution later. D(4) has the form

D(4) =







−A1E1−f1

.......

.

−AµE(3)

µE(4)

µ−fµ

−A(3)

µ+1 −A(4)

µ+1 E(1)

µ+1 −fµ+1

1

−˜

S(2)

µ+1 0 0 −˜

fµ+1

2

0 0 0 −˜

fµ+1

3







,

where ˜

S(2)

µ+1 has full row and full column rank hf,0and EµV(2)

µ+1 =hE(3)

µE(4)

µi. Thus, we can

eliminate the entries above −˜

S(2)

µ+1 which yields the matrix

D(5) =







−A1E1−f1

.......

.

−Aµ0E(4)

µ−˜

fµ

0−A(4)

µ+1 E(1)

µ+1 −˜

fµ+1

1

−˜

S(2)

µ+1 0 0 −˜

fµ+1

2

0 0 0 −˜

fµ+1

3







.

This is done in lines 337-345. One can see that we have computed the form (3.13) for the

single matrix pair (Eµ+1, Aµ+1). This process, which leads from D(2) to D(5), is then repeated

for the next matrix pair (h0E(4)

µi, Aµ). This is done over and over again until the final

upper left matrix pair is reached. This final matrix pair is then handled in the lines 347 -

396. Call this final matrix D(6).

We see that due to the structure of the main program all matrix pairs in D(1) already have

been reduced at least once. The process leading from D(1) to D(6) reduces all matrix pairs

once again. Thus, the upper left matrix pair in D(6) has been reduced µ+ 1 times. As

101

we know from (3.26) µ+ 1 reductions suffice to obtain a strangeness-free from, once µis

the forward strangeness index. Thus, the upper left matrix pair in D(6) can be used to

compute one iterate of the solution, which is done in lines 237 - 252. In these lines the

consistency conditions on the inhomogeneity are ignored and the undetermined components

of the solution are set to zero.

5.3.1 Numerical results

To present numerical results, we consider a linear differential-algebraic equation from [8],

namely

0 0

1−t˙x(t) = −1t

0 0x(t) + tsin(t)

t+ cos(t), for t∈[−7; 7]. (5.8)

As shown in [8] on page 57,

x(t) = t2+tcos(t)−t2cos(t)

t+ cos(t)−sin(t)−tcos(t)(5.9)

is the only solution of (5.8). This means that there is also only one consistent initial condition.

Note that (5.8) is a special case of (1.7). Thus, we know that by discretizing (5.8) with the

explicit Euler method and the equidistant grid

. . . < −2h < −h < 0< h < 2h < . . . , (5.10)

for some step size h∈R+yields the discrete-time descriptor system

0 0

1/h −kxk+1 =−1kh

1/h −kxk+kh sin(kh)

kh + cos(kh),k∈Z. (5.11)

Solving this discrete-time descriptor system with the algorithm from Appendix A for a given

hon the discrete interval Z∩[−7

h;7

h] yields an approximation to the actual solution (5.9).

Figures 5.1 and 5.2 show the exact solution (5.9) and the approximation for several step sizes

h. One can clearly see how the approximation gets better as the step sizes gets smaller.

In any case only one solution was found and although an inconsistent initial condition was

given to the algorithm from Appendix A (parameter xhat) the algorithm found the only

consistent initial condition.

In table 5.1 the runtime and the approximation errors of the algorithm form Appendix A

are shown for different step sizes h. The runtime in seconds relates to a 64bit processor

from Intel with 2.13GHz. The ∅solution difference value is the average error between the

iterates of the solution computed from (5.11) with the help of the algorithm from Appendix

A and the exact solution of (5.8), i.e., (5.9). The max. solution difference value is the

maximum error between the iterates of the solution computed from (5.11) with the help of

the algorithm from Appendix A and the exact solution of (5.8), i.e., (5.9).

One clearly sees how the errors are proportional to the step size h. Also the runtimes seem

to be inverse proportional to the step size h. Thus, the explicit Euler method is convergent

of order one, when applied to the differential-algebraic equation (5.8).

102

−8 −6 −4 −2 0 2 4 6 8

−10

−5

0

5

10

15

20

25

30

35

40

h = 1

computed solution

real solution

−8 −6 −4 −2 0 2 4 6 8

−10

−5

0

5

10

15

20

25

30

h = 0.5

computed solution

real solution

Figure 5.1: Discretization with h = 1 and h = 0.5

−8 −6 −4 −2 0 2 4 6 8

−10

−5

0

5

10

15

20

25

30

h = 0.1

computed solution

real solution

−8 −6 −4 −2 0 2 4 6 8

−10

−5

0

5

10

15

20

25

30

h = 0.01

computed solution

real solution

Figure 5.2: Discretization with h = 0.1 and h = 0.01

h runtime in

seconds

∅solution difference max. solution difference

1 0.01305 6.7431 31.945

0.5 0.02230 2.9053 13.926

0.1 0.09910 0.51967 2.3753

0.05 0.1944 0.2565 1.1601

0.01 0.9677 0.050795 0.22757

0.001 10.96 0.0050684 0.022657

0.0001 362.75 0.00050673 0.002265

Table 5.1: Runtime and approximation errors for different step sizes h

103

Chapter 6

Summary

In this diploma thesis some concepts from the book [8] are transferred to the discrete-

time case. In chapter 2 we discuss linear discrete-time descriptor systems with constant

coefficients. This is done in terms of the Kronecker canonical form and afterwards in explicit

form in analogy to [8]. In chapter 2 we already note that one may distinguish three cases

of (linear) discrete time descriptor systems. The first case is the forward case, i.e., the case

where one has an initial condition given at point k0∈Zand only wants to get a solution for

indices k≥k0. The second case is the backward case, i.e., the case where one has an initial

condition given at point k0∈Zand only wants to get a solution for indices k≤k0. These

first two cases are closely related, since the first case can be transferred into the second one

by a variable substitution. The third case is the two-way case, i.e., the case where an initial

condition is given at some point k0∈Zbut one is looking for a solution for indices k≥k0as

well as for indices k≤k0. This case puts stronger restrictions on the initial condition, i.e.,

the set of consistent initial conditions in the third case is smaller than in the first or second

case.

In chapter 3 we then move on to descriptor systems with variable coefficients. We first

identify the local characteristics of such systems. Then, with the help of some constant rank

assumptions, we derive a canonical form, which allows for statements about the existence

and uniqueness of solutions for variable coefficient descriptor systems in the forward case. To

achieve this canonical form one has to employ global equivalence transformations and shifts

of equations. The shift is the discrete-time analogon to the differentiation in the continuous-

time case. A strangeness index is defined, which counts the number of shifts necessary to

obtain the canonical form. In section 3.2.1 the inflated descriptor system is introduced and it

is shown (analogous to [8]) that, under some constant rank assumptions, one can determine

the characteristics of the original systems once one knows the local invariants of the inflated

systems. Section 3.3 discusses how to carry over the results from the forward case to the

backward case.

Chapter 4 is devoted to the combination of the forward and backward base. Here even

stronger constant rank assumptions are introduced, which make it possible to obtain a two-

way canonical form. With this two-way canonical form one can make statements about

the existence and uniqueness of solutions in the two-way case. A new strangeness index is

104

defined for systems, which satisfy this stronger constant rank assumptions.

Finally, in chapter 5 two algorithms for descriptor systems with constant coefficients and

one algorithm for descriptor systems with variable coefficients are introduced. The first

algorithms will never be used. The second one may be fast enough but it does not respect

the actual index of the descriptor system. The third algorithm seems to be the appropriate

one.

6.1 Outlook

We have only considered the general linear discrete-time descriptor system, but other special

cases might also be interesting, for example the periodic time-variant case, which is consid-

ered in [17]. Another interesting topic would be polynomial descriptor systems, i.e., systems

of the form

A(l)

kxk+l+...+A(1)

kxk+1 +A(0)

kxk= 0 for k∈Z,

although the theory might get very complicated for such systems. For such systems one

could also investigate the associated control problem, which has already been done for non-

polynomial descriptor systems (e.g., [10, 15, 22]).

The constant rank assumptions (3.12) make it possible to obtain the canonical form (3.26).

One can relax these assumptions, similar to Hypothesis 3.48 in [8]. Then, it should be possible

to create an algorithm that computes the solution of a linear discrete-time descriptor system

under these relaxed assumptions.

It would be interesting to investigate the discretization of linear continuous-time descriptor

systems (i.e., differential-algebraic equations) and connections between properties of the

original continuous-time descriptor system and the discretized descriptor system. Another

thing one could do is the derivation of a canonical form similar to the one in Theorem 3.20

for the two-way reduction process. Having obtained such a canonical form one could try to

prove Conjecture 4.14.

Other topics, which are commonly examined for linear systems of the form xk+1 =Akxk+fk,

also make sense for discrete-time descriptor systems. Such topics include Stability and

Stochastic Systems (see [7, 9]).

Finally, of course, non-linear descriptor systems also have to be considered.

Acknowledgements

I want to thank Volker Mehrmann for numerous revisions of the manuscript and continuous

support during its creation and Peter Kunkel for the discussion in which he noticed that

the form (3.26) only allows for statements about the existence and uniqueness of forward

solutions.

105

Appendix A

Matlab Code

Note that although the following matlab code works, there are major improvements neces-

sary. For example, it is not always necessary to compute the singular value decompositions of

the whole matrix in the function advance inflated system. Also, there may occur an division

by zero error in line 243 and line 199 for underdetermined systems, which can be avoided.

The file can be downloaded from http://www.math.tu-berlin.de/~bruell/matlab/.

1% solve_dds Compute a solution of a discrete - time descriptor system .

2% solve_dds (Efun ,Afun ,ffun ,kb ,k0 ,kf , xhat , tol) tries to compute the

3% iterates x_kb , ... , x_kf of a solution of one of the discrete -time

4% descriptor systems

5%

6% (1) E_k x_{k +1} = A_k x_{k} + f_k for kb <= k <= infty

7% (2) E_k x_{k +1} = A_k x_{k} + f_k for -infty <= k <= kf

8% (3) E_k x_{k +1} = A_k x_{k} + f_k for -infty <= k <= infty

9%

10 % along with xhat = x_{k0},

11 % where the matrices E_k and A_k and the vector f_k are obtained by

12 % evaluation the function Efun , Afun and ffun , respectively , i.e.

13 % E_k = feval( Efun , k ), A_k = feval ( Afun , k ),

14 % f_k = feval( ffun , k ).

15 % If no such solution exists , the right hand side f_k is changed so

16 % that a solution exists . If there are multiple solutions , one

17 % solution is selected . If the initial value ’xhat ’ is inconsistent ,

18 % the consistent initial value is chosen which is closest to

19 % ’xhat ’ ( in the 2- norm ).

20 % If the constant rank assumptions form [1] are not satisfied an

21 % error message is issued and the algorithm is aborted , although

22 % there might still exist a solution .

23 % If kb == k0 and k0 < kf equation (1) is considered .

24 % If kb < k0 and k0 == kf equation (2) is considered .

25 % If kb < k0 and k0 < kf equation (3) is considered .

26 % If kb == k0 and k0 == kf no equation is considered and no real

27 % computations are performed .

28 %

29 % Input parameters:

30 %

106

31 % Efun : A function handle or name of a function that takes

32 % one input parameter k and returns the matrix E_k .

33 %

34 % Afun : A function handle or name of a function that takes

35 % one input parameter k and returns the matrix A_k .

36 %

37 % ffun : A function handle or name of a function that takes

38 % one input parameter k and returns the vector f_k .

39 %

40 % kb : The index of the first iterate to compute .

41 % k0 : The index , where to apply the initial condition .

42 % k0 has to satisfy kb <= k0 <= kf.

43 % kf : The index of the last iterate to compute .

44 %

45 % xhat : The desired value of the solution at iterate x_k0 .

46 %

47 % tol : Tolerance , below which a singular value is considered

48 % zero .

49 %

50 % Output parameters

51 %

52 % x : The solution iterates x_kb , ... , x_kf side by side

53 % in a matrix . Thus , x (: ,k- kb +1) represents the x_k

54 % iterate of the solution . x(: ,1) represents the x_kb

55 % iterate and x(kf -kb +1) represents the x_kf iterate .

56 %

57 % isunique : This is set to 1 if there is only one unique

58 % solution and to 0 otherwise .

59 %

60 % r_f : The sequence of the r_f as in [1]. Note that r_f(i)

61 % corresponds to r_{f,i -1} in the forward - only

62 % reduction process as described in [1].

63 %

64 % h_f : The sequence of the h_f as in [1]. Note that h_f(i)

65 % corresponds to h_{f,i -1} in the forward - only

66 % reduction process as described in [1].

67 %

68 % mu_f : The forward strangeness index of the system as

69 % defined in [1].

70 %

71 % r_b : The sequence of the r_b as in [1]. Note that r_b(i)

72 % corresponds to r_{b,i -1} in the backward - only

73 % reduction process as described in [1].

74 %

75 % h_b : The sequence of the h_b as in [1]. Note that h_b(i)

76 % corresponds to h_{b,i -1} in the backward - only

77 % reduction process as described in [1].

78 %

79 % mu_b : The backward strangeness index of the system as

80 % defined in [1].

81 %

82 % Reference:

107

83 % [1] Bruell , T.

84 % Linear discrete - time descriptor systems ;

85 % Diploma - Thesis (2007);

86 % http :// www . math .tu - berlin .de /~ bruell

87 %

88 function [x, isunique ,r_f ,h_f ,mu_f ,r_b ,h_b ,mu_b ]=...

89 solve_dds(Efun ,Afun ,ffun ,kb ,k0 ,kf ,xhat ,tol )

90

91 r_f =[]; h_f =[];

92 mu_f = 0;

93 r_b =[]; h_b =[];

94 mu_b = 0;

95

96 dimmat = feval( Efun , k0 );

97 [m,n] = size ( dimmat );

98 clear dimmat;

99

100 do_forward = ( k0 < kf );

101 do_backward = ( k0 > kb );

102

103 if( do_forward )

104 Ek_f = {};

105 Ak_f = {};

106 Qk_f = {};

107 fk_f = {};

108

109 kact = k0;

110

111 % determine the forward index

112 haveindex = 0;

113 while ( ~ haveindex )

114 [Ek_f ,Ak_f ,Qk_f ,fk_f ,r_f , h_f ]= advance_inflated_system ...

115 (Efun ,Afun ,ffun ,m,n ,Ek_f ,Ak_f ,Qk_f ,fk_f ,kact ,mu_f ,r_f ,h_f ,tol ,0);

116

117 if( ( mu_f >= 1 && r_f ( mu_f +1) == r_f(mu_f ) ) )

118 haveindex = 1;

119 else

120 mu_f = mu_f + 1;

121 end

122 end

123

124 r_mu_f = r_f ( mu_f +1);

125 h_mu_f = h_f ( mu_f +1);

126 else

127 r_mu_f = -1; h_mu_f = -1;

128 end

129

130 if( do_backward )

131 Ek_b = {};

132 Ak_b = {};

133 Qk_b = {};

134 fk_b = {};

108

135

136 kact = k0 -1;

137

138 % determine the backward index

139 haveindex = 0;

140 while ( ~ haveindex )

141 [Ek_b ,Ak_b ,Qk_b ,fk_b ,r_b , h_b ]= advance_inflated_system ...

142 (Afun ,Efun ,ffun ,m,n ,Ek_b ,Ak_b ,Qk_b ,fk_b ,kact ,mu_b ,r_b ,h_b ,tol ,1);

143

144 if( ( mu_b >= 1 && r_b ( mu_b +1) == r_b(mu_b ) ) )

145 haveindex = 1;

146 else

147 mu_b = mu_b + 1;

148 end

149 end

150

151 r_mu_b = r_b ( mu_b +1);

152 h_mu_b = h_b ( mu_b +1);

153 else

154 r_mu_b = -1; h_mu_b = -1;

155 end

156

157 % find all constraints that the initial condition has to fulfill

158 constraintA = zeros (0 ,n );

159 constraintf = zeros (0 ,1);

160 if( do_forward )

161 constraintA = [ constraintA ;- Ak_f {1}( r_mu_f +1: r_mu_f +h_mu_f ,:)...

162 *Qk_f {1} ’];

163 constraintf = [ constraintf ; fk_f {1}( r_mu_f +1: r_mu_f + h_mu_f )];

164 end

165 if( do_backward )

166 constraintA = [ constraintA ; Ak_b {1}( r_mu_b +1: r_mu_b +h_mu_b ,:)...

167 *Qk_b {1} ’];

168 constraintf = [ constraintf ; fk_b {1}( r_mu_b +1: r_mu_b + h_mu_b )];

169 end

170

171 % find the consistent initial condition that is closest to the given

172 % one

173 [U,S ,V] = svd ( constraintA );

174 cnum = rank(S); % compute the number of constraints

175 cnum

176

177 ytilde = zeros(n,1);

178 xtilde = V’ * xhat ;

179 ftilde = U’ * constraintf ;

180

181 ytilde (1: cnum) = S(1: cnum ,1: cnum) \ ftilde (1: cnum );

182 ytilde ( cnum +1: n) = xtilde ( cnum +1: n );

183

184 xhat = V * ytilde;

185

186 if( do_forward )

109

187 kact = k0;

188

189 xtt_old = (Qk_f {1}) ’ * xhat ;

190

191 % start forward solver

192 for i =1:(kf -k0)

193 xtt = zeros (n ,1);

194

195 xtt (1: h_mu_f ) = -Ak_f {2}( r_mu_f +1: r_mu_f +h_mu_f ,1: h_mu_f )\...

196 fk_f {2}( r_mu_f +1 : ( r_mu_f + h_mu_f ) );

197

198 xtt ( h_mu_f +1: h_mu_f + r_mu_f ) = ...

199 Ek_f {1}(1: r_mu_f , h_mu_f +1: h_mu_f + r_mu_f )\(...

200 Ak_f {1}(1: r_mu_f ,:) * xtt_old ...

201 + fk_f {1}(1: r_mu_f ) ...

202 - Ek_f {1}(1: r_mu_f ,:) * xtt );

203

204 % choose a solution

205 xtt ( h_mu_f + r_mu_f +1: n) = zeros (n- r_mu_f - h_mu_f ,1);

206

207 x(1:n, k0 -kb+i+1) = ( Qk_f {2}) * xtt ;

208 xtt_old = xtt;

209

210 % proceed one step

211 kact = kact +1;

212 for i =1: prod (size ( Ek_f )) -1

213 Ek_f {i }= Ek_f {i +1};

214 Ak_f {i }= Ak_f {i +1};

215 fk_f {i }= fk_f {i +1};

216 Qk_f {i }= Qk_f {i +1};

217 end

218 [Ek_f ,Ak_f ,Qk_f ,fk_f ,r_f , h_f ]= advance_inflated_system ...

219 (Efun ,Afun ,ffun ,m,n ,Ek_f ,Ak_f ,Qk_f ,fk_f ,kact ,mu_f ,r_f ,h_f ,tol ,0);

220 end

221

222 % compute the _real_ forward strangeness index (i.e. as in [1])

223 for i= length ( r_f ) -1: -1:1

224 if( r_f (i) == r_f(i +1) )

225 mu_f = i -1;

226 end

227 end

228 end

229

230 if( do_backward )

231 kact = k0 -1;

232

233 xtt_old = (Qk_b {1}) ’ * xhat ;

234

235 % start backward solver

236 for i =0:(k0 -kb -1)

237 xtt = zeros (n ,1);

238

110

239 xtt (1: h_mu_b ) = Ak_b {2}( r_mu_b +1: r_mu_b +h_mu_b ,1: h_mu_b )\...

240 fk_b {2}( r_mu_b +1 : ( r_mu_b + h_mu_b ) );

241

242 xtt ( h_mu_b +1: h_mu_b + r_mu_b ) = ...

243 Ek_b {1}(1: r_mu_b , h_mu_b +1: h_mu_b + r_mu_b )\( ...

244 Ak_b {1}(1: r_mu_b ,:) * xtt_old ...

245 - fk_b {1}(1: r_mu_b ) ...

246 - Ek_b {1}(1: r_mu_b ,:) * xtt );

247

248 % choose a solution

249 xtt ( h_mu_b + r_mu_b +1: n) = zeros (n- r_mu_b - h_mu_b ,1);

250

251 x(1:n, k0 -kb -i) = (Qk_b {2}) * xtt;

252 xtt_old = xtt;

253

254 % proceed one step

255 kact = kact -1;

256 for i =1: prod (size ( Ek_b )) -1

257 Ek_b {i }= Ek_b {i +1};

258 Ak_b {i }= Ak_b {i +1};

259 fk_b {i }= fk_b {i +1};

260 Qk_b {i }= Qk_b {i +1};

261 end

262 [Ek_b ,Ak_b ,Qk_b ,fk_b ,r_b , h_b ]= advance_inflated_system ...

263 (Afun ,Efun ,ffun ,m,n ,Ek_b ,Ak_b ,Qk_b ,fk_b ,kact ,mu_b ,r_b ,h_b ,tol ,1);

264 end

265

266 % compute the _real_ backward strangeness index (i.e. as in [1])

267 for i= length ( r_b ) -1: -1:1

268 if( r_b (i) == r_b(i +1) )

269 mu_b = i -1;

270 end

271 end

272 end

273

274 if( ( ~ do_forward || r_mu_f + h_mu_f == n ) && ...

275 ( ~ do_backward || r_mu_b + h_mu_b == n ) )

276 isunique = 1;

277 else

278 isunique = 0;

279 end

280

281 x(:,k0 -kb +1) = xhat ;

282

283

284

285 function [Ek ,Ak ,Qk ,fk ,r,h]= advance_inflated_system ...

286 (Efun , Afun , ffun ,m ,n,Ek ,Ak ,Qk ,fk ,kact ,mu ,r,h,tol , backward )

287

288 if( backward )

289 signed_mu = -mu;

290 else

111

291 signed_mu = mu;

292 end

293

294 Ek {mu +1} = feval ( Efun , kact + signed_mu );

295 Ak {mu +1} = feval ( Afun , kact + signed_mu );

296 Qk{mu +1} = eye (n);

297 fk {mu +1} = feval ( ffun , kact + signed_mu );

298

299 for j=mu : -1:1

300 % calculate a svd of E

301 [U,S,V] = svd (Ek {1+j }); % Ek == U * S * V’

302

303 % determine the rank ( only for sure )

304 loc_r = sum( getdiag (S) > tol );

305 if( loc_r ~= r(mu -j+1) )

306 error ([’r not invariant ! (k=’, num2str ( kact),’,loc_r=’ ,...

307 num2str ( loc_r ),’,r(’, num2str (mu -j +1) , ’)= ’ ,...

308 num2str (r(mu -j+1)) , ’)’]);

309 end

310

311 % apply transformation

312 S( ( loc_r +1):m, : ) = zeros (m- loc_r ,n);

313 Ek {1+ j} = S * V’; % = U ’* Ekact

314 Ak {1+j} = U ’* Ak {1+j};

315 fk {1+j} = U ’* fk {1+j};

316

317 % calculate a svd of A

318 A2 = Ak {1+ j}( loc_r +1:m, : );

319 [U,S,V] = svd(A2 ); % Ak == U * S * V’

320

321 % determine the rank ( only for sure )

322 loc_h = sum( getdiag (S) > tol );

323 if( loc_h ~= h(mu -j+1) )

324 error ([’h not invariant ! (k=’, num2str ( kact),’,loc_h=’ ,...

325 num2str ( loc_h ),’,h(’, num2str (mu -j +1) , ’)= ’ ,...

326 num2str (h(mu -j+1)) , ’)’]);

327 end

328

329 % apply transformation

330 S( loc_h +1:m-loc_r , : ) = zeros (m- loc_r -loc_h ,n);

331 Ek{j} = Ek{j}*V;

332 Ak {1+ j}( loc_r +1:m, : ) = S;

333 Ak {1+ j}( 1: loc_r , : ) = Ak{1+j}( 1: loc_r , : ) * V;

334 fk {1+ j}( loc_r +1:m, : ) = U’* fk {1+ j}( loc_r +1:m, : );

335 Qk {1+ j} = Qk {1+ j}* V;

336

337 % do block elimination (in A)

338 fk {1+j}(1: loc_r ) = fk {1+j}(1: loc_r) - Ak {1+j}(1: loc_r ,1: loc_h )*...

339 (S (1: loc_h ,1: loc_h )\ fk {1+ j}(( loc_r +1):( loc_r+loc_h )));

340 Ak {1+j }(1: loc_r ,1: loc_h ) = zeros (loc_r , loc_h );

341

342 % do block elimination (in E)

112

343 fk{j }(1: loc_r ) = fk{j}(1: loc_r) + Ek{j }(1: loc_r ,1: loc_h )*...

344 (S (1: loc_h ,1: loc_h )\ fk {1+ j}(( loc_r +1):( loc_r+loc_h )));

345 Ek {j }(1: loc_r ,1: loc_h) = zeros (loc_r , loc_h );

346

347 end

348

349 [U,S,V] = svd (Ek {1}); % Ek == U * S * V’

350

351 % determine the rank

352 loc_r = sum ( getdiag (S) > tol );

353 if( length (r) >= mu+1 )

354 if ( loc_r ~= r( mu +1) )

355 error ([’r not invariant ! (k=’, num2str ( kact),’,loc_r=’ ,...

356 num2str ( loc_r ),’,r(’, num2str ( mu +1), ’)=’ ,...

357 num2str (r( mu +1)) , ’) ’]);

358 end

359 else

360 r( mu +1) = loc_r ;

361 end

362

363 % apply transformation

364 S( r(mu +1)+1:m, : ) = zeros(m-r(mu+1),n);

365 Ek {1} = S * V ’; % = U ’* Ekact

366 Ak {1} = U ’* Ak {1};

367 fk {1} = U ’* fk {1};

368

369 A2 = Ak {1}( r(mu +1)+1:m, : );

370 [U,S,V] = svd (A2 ); % Ak == U * S * V’

371

372 % determine the rank

373 loc_h = sum ( getdiag (S) > tol );

374 if( length (h) >= mu+1 )

375 if ( loc_h ~= h( mu +1) )

376 error ([’h not invariant ! (k=’, num2str ( kact),’,loc_h=’ ,...

377 num2str ( loc_h ),’,h(’, num2str ( mu +1), ’)=’ ,...

378 num2str (h( mu +1)) , ’) ’]);

379 end

380 else

381 h( mu +1) = loc_h ;

382 end

383

384 % apply transformation

385 S( h( mu +1)+1:m-r(mu +1) , : ) = zeros(m-r(mu+1)- h(mu+1), n);

386 Ek {1} = Ek {1};

387 Ak {1}( r(mu +1)+1:m, : ) = S;

388 Ak {1}( 1:r(mu +1) , : ) = Ak {1}( 1:r(mu+1), : ) * V;

389 fk {1}( r(mu +1)+1:m, : ) = U’* fk {1}( r(mu +1)+1:m, : );

390 Qk {1} = Qk {1} * V;

391

392 % eliminate in A

393 fk {1}(1: r(mu +1)) = fk {1}(1: r(mu +1)) - Ak {1}(1: r( mu +1) ,1: h(mu +1))*...

394 (S (1: h( mu +1) ,1: h(mu +1))\ fk {1}(( r(mu +1)+1):( r( mu +1)+ h(mu +1))));

113

395 Ak {1}(1: r( mu +1) ,1: h( mu +1)) = zeros (r( mu +1),h( mu +1));

396

397

398

399 % Returns the diagonal elements of a matrix in a vector .

400 function x=getdiag(A)

401

402 if( size(A,1) > 1 && size (A ,2) > 1 )

403 x = diag (A);

404 else

405 if( prod ( size(A)) > 0 )

406 x = A(1 ,1);

407 else

408 x = [];

409 end

410 end

411

412

413 % Dumps the inflated descriptor system represented by the arrays

414 % Ek , Ak and fk.

415 function dumparray (Ek ,Ak ,fk ,m,n)

416 outsys = zeros (0 ,n );

417 outf = [];

418 for i=1: prod ( size(Ek))

419 outsys = [ outsys , zeros ((i -1)*m,n ); ...

420 zeros(m,(i -1)* n) , -Ak{i}, Ek{i} ];

421 outf = [ outf ;-fk{i}];

422 end

423 [ outsys , outf ]

114

Bibliography

[1] S. Campbell,Singular Systems of Differential Equations I, Pitman, San Francisco,

1980.

[2] S. Campbell and J. C.D. Meyer,Generalized Inverses of Linear Transformations,

General Publishing Company, 1979, ch. 9.3 and 9.4, pp. 181–187.

[3] S. L. Campbell,Nonregular singular dynamic Leontief systems, Econometrica, Vol.

47 (1979), pp. 1565–1568.

[4] J. Demmel and B. K˚

agstr¨

om,The generalized Schur decomposition of an arbitrary

pencil A - λB: Robust software with error bounds and applications. part I and II, ACM

Transactions on Mathematical Software, Vol. 19 (1993), pp. 160–201.

[5] T. N. E. Greville,The Souriau-Frame algorithm and the Drazin pseudoinverse,

Linear Algebra and its Applications, (1973), pp. 205–208.

[6] D. D. Hai,Phuong trinh sai phan tuyen tinh voi he so ca co hang thay doi (in viet-

namese; translates: On linear implicit nonautonomous difference equations), master’s

thesis, College of Natural Science, Vietnam National University, 2006.

[7] C. Heij, A. Ran, and F. van Schagen,Introduction to Mathematical Systems

Theory - Linear Systems, Identification and Control, Birkh¨auser Verlag AG, Basel,

2006.

[8] P. Kunkel and V. Mehrmann,Differential-Algebraic Equations - Analysis and Nu-

merical Solution, European Mathematical Society, Z¨urich, 2006.

[9] V. Kuˇ

cera,Analysis and Design of Discrete Linear Control Systems, Prentice Hall,

Upper Saddle River, NJ, USA, 1991.

[10] F. L. Lewis,Fundamental, reachability and observability matrices for discrete descrip-

tor systems, IEEE Transactions on automatic control, Vol. 30 (1985), pp. 502–505.

[11] C. F. V. Loan,The ubiquitous Kronecker product, Journal of computational and

applied mathematics, (2000), pp. 85–100.

[12] D. G. Luenberger,Dynamic equations in descriptor form, IEEE Transactions on

Automatic Control, Vol. 22 (1977), pp. 312–321.

115

[13] , Time-invariant descriptor systems, Automatica, Vol. 14 (1978), pp. 473 – 480.

[14] D. G. Luenberger and A. Arbel,Singular dynamic Leontief systems, Econome-

trica, Vol. 45 (1977), pp. 991–995.

[15] V. Mehrmann,The Autonomous Linear Quadratic Control Problem, Springer-Verlag,

Berlin, 1991.

[16] B. Mertzios and F. Lewis,Fundamental matrix of discrete singular systems, Cir-

cuits, Systems, and Signal Processing, Vol. 8 (1989), pp. 341–355.

[17] J. Sreedhar and P. V. Dooren,Periodic descriptor systems: Solvability and con-

ditionability, IEEE Transactions on Automatic Control, Vol. 44 (1999), pp. 310–313.

[18] P. S. Stanimirovic and M. D. Petkovic,Computing generalized inverse of poly-

nomial matrices by interpolation, Applied Mathematics and Computation, 172 (2006),

pp. 508–523.

[19] T. Stykel,Numerical solution and perturbation theory for generalized Lyapunov equa-

tions, Linear Algebra and its Applications, (2002), pp. 155–185.

[20] Y. Wei and H. Wu,The representation and approximation for Drazin inverse, Journal

of computational and applied mathematics, (2000), pp. 417–432.

[21] J. H. Wilkinson,Note on the practical significance of the Drazin inverse, tech. report,

School of Humanities and Sciences, Stanford University, April 1979.

[22] L. Zhang, J. Lam, and Q. Zhang,Lyapunov and Riccati equations of discrete-time

descriptor systems, IEEE Transactions on automatic control, Vol. 44 (1999), pp. 2134–

2139.

116