FACTA UNIVERSITATIS (NIŠ)

Ser. Math. Inform. 19 (2004), 73–92

AN EXTENSION OF LEVERRIER-FADDEEV ALGORITHM

USING A BASIS OF CLASSICAL ORTHOGONAL
POLYNOMIALS

J. Hernández and F. Marcellán

To our friend Giuseppe Mastroianni with occasion of his 65th anniversary

Abstract. In this paper, the Leverrier-Faddeev algorithm is extended in order

to compute both the determinant a(s) and the adjoint B(s) of the polynomial
matrix sp Ap + sp−1 Ap−1 + · · · + sA1 + A0 , where Ai ∈ Cn×n , i = 1, 2, . . . , p. B(s)
and a(s) are expressed in terms of a basis of classical orthogonal polynomials.

1. Introduction

When one considers a vibrating system consisting of N masses connected

by linear springs of stiffness (kr )N
r=1 and the whole lies in a straight line on
a smooth horizontal table and is excited by forces (Fr (t))N r=1 , the Newton’s
equations of motion for the system are

mr ür = Fr + θr+1 − θr , r = 1, 2, . . . , N − 1,
(1.1)
mN üN = FN − θN .

Hooke’s law states that the spring forces are given by

(1.2) θr = kr (ur − ur−1 ), r = 1, 2, . . . , N.

If the left hand end is pinned, then u0 = 0.

Received September 28, 2004.
2000 Mathematics Subject Classification. Primary 42C05; Secondary 33C45, 15A15.

73
74 J. Hernández and F. Marcellán

Forced vibration analysis concerns the solution of these equations for

given forcing functions Fr (t).
Free vibration analysis consists in finding solutions to the equations which
require no external excitation, i.e., Fr (t) ≡ 0, r = 1, 2, . . . , N , and which
satisfy the stated end conditions (see [5]).
In order to express equations (1.1)–(1.2) in matrix form, we get

(1.3) Aü + Cu = F,

where A = diag (m1 , m2 , . . . , mN ), F = (F1 , F2 , . . . , FN ),

 
k1 + k2 −k2 0 ··· 0
 . . .. 
 −k2 k 2 + k 3 −k 3 . . 
 
 
C=  .. .. .. 

 0 . . . 0 
 .. . .. 
 . −kN 
0 ··· 0 −kN kN

and u = (u1 , u2 , . . . , uN )T . A and C are called, respectively, the inertia and

stiffness matrices of the system.
If we assume initial conditions u(0) and u0 (0) are given, then using the
Laplace transform in (1.3) we get

(s2 A + C)ũ = F̃ (s) + sB + G,

where ũ denotes the Laplace transform of the vector function u. Thus,

(1.4) ũ = (s2 A + C)−1 F̃ (s) + (s2 A + C)−1 (sB + G).

In order to solve the problem, we need to find the inverse of the polyno-
mial matrix s2 A + C (see [1]). This question was analyzed in a more general
context by several authors. In particular, in [4] an algorithm has been given
to find the determinant of a general matrix polynomial as well as its adjoint
matrix. The connection with the Leverrier’s algorithm to find the determi-
nant of sIN − A is shown. The aim of our contribution is twofold. First,
we use an extension of the Leverrier’s algorithm stated in [6], [7] for polyno-
mial bases that is a more efficient computational method than the preceding
one taking into account the use of the canonical basis in the linear space
P of polynomials with complex coefficients. Second, we give an alternative
 An Extension of Leverrier-Faddeev Algorithm . . . 75

approach to the results of [4], [10] assuming the leading coefficient of the
polynomial matrix is any matrix, not necessary the unit matrix as there.
Notice that the inverse of the polynomial matrix in (1.4) is the transfer
function of the linear system (1.1)–(1.2).
The structure of the paper is the following. In Section 2, we give the basic
background about scalar orthogonal polynomials with a special emphasis
in the so called classical orthogonal polynomials. Then, we describe the
Leverrier’s algorithm taking into account a polynomial basis constituted by
classical orthogonal polynomials.
In Section 3 we introduce polynomial matrices as well as some special ex-
amples (orthogonal matrix polynomials), which are generated recursively. In
Section 4 we give an algorithm to find the determinant of a polynomial ma-
trix as well as the adjoint matrix. Furthermore, in Section 5, some examples
are presented.

2. Classical Orthogonal Polynomials

Given a linear functional u in the linear space P of polynomials with

complex coefficients, we call moments (un ) the complex numbers un =
hu, xn i, n ∈ N. Here h·, ·i denotes the duality bracket.

Definition 2.1. The linear functional u is said to be quasi-definite if the

principal submatrices of the infinite Hankel matrix H = (uj+k )+∞j,k=0 are
nonsingular.

Thus, we can introduce an inner product associated with the quasi-

definite linear functional u as follows:

φ(p, q) = hu, pqi, p, q ∈ P.

Proposition 2.1. ([2]) The following statements are equivalent:

(i) u is a quasi-definite linear functional.
(ii) There exists a sequence of monic polynomials (Pn ) such that

– deg Pn = n,
76 J. Hernández and F. Marcellán

– φ(Pn , Pm ) = 0 if n 6= m,

– φ(Pn , Pn ) 6= 0 for every n ∈ N.

(iii) There exist sequences of complex numbers (βn ) and γn with γn 6= 0

for every n ∈ N such that the sequence of monic polynomials (Pn ) defined by

(2.1) xPn (x) = Pn+1 (x) + βn Pn (x) + γn Pn−1 (x), n > 1

satisfy hu, Pn Pm i = 0, n 6= m, and hu, Pn 2 i 6= 0, for some linear functional u.

The matrix representation of the operator associated with the multipli-

cation by x in P is a tridiagonal matrix
 
β0 1 0 ···
 
 γ1 β1 1 . . . 
 
J =  0 γ 2 β2 1 .

 .. . . . . . . . . 
 . . . . . 

Notice that the characteristic polynomial of the principal submatrix Jn

of dimension n is the polynomial Pn . The recurrence relation (2.1) reads
as follows: the determinant of the matrix polynomial sIn+1 − Jn+1 can be
given recursively as a linear combination of the determinants of the matrix
polynomials sIn −Jn and sIn−1 −Jn−1 . Examples of orthogonal polynomials
with respect to a linear functional appear in Sturm-Liouville problems for
second and fourth order linear differential equations with polynomial coef-
ficients. More precisely, if u is a quasi-definite linear functional, then the
following statements are equivalent

Proposition 2.2. ([8]) (i) There exist polynomials φ and ψ with deg φ 6 2
and deg ψ = 1 such that the distributional equation

(2.2) D(φu) = ψu,

holds.
(ii) The sequence of monic polynomials (Pn ) orthogonal with respect to u
satisfies a second order linear differential equation

(2.3) φPn00 + ψPn0 + λn Pn = 0.

 An Extension of Leverrier-Faddeev Algorithm . . . 77

The equation (2.2) is called a Pearson differential equation. (2.3) means

essentially that (Pn ) is a hypergeometric family of polynomials. The linear
functional u satisfying a Pearson equation is said to be classical. The se-
quence of orthogonal polynomials with respect to u is said to be a classical
sequence of orthogonal polynomials.
For our purposes, we can use other two characterizations of classical or-
thogonal polynomials.

Proposition 2.3. ([8]) The following statements are equivalent:

(i) (Pn ) is a sequence of classical orthogonal polynomials.
(ii) The sequence of monic derivatives of (Pn ) is also a sequence of or-
P0
thogonal polynomials. We will denote Qn = n+1 .
n+1
(iii) In the expansion of Pn in terms of the polynomial basis (Qn ), the
coefficients for 0 6 k 6 n − 3 vanish. In other words,

(2.4) Pn (x) = Qn (x) + rn Qn−1 (x) + sn Qn−2 (x).

The sequences of classical orthogonal polynomials, up to a linear change

of variables, are

1. Hermite polynomials (Hn )

0
Hn+1 (x)
Hn (x) = .
n+1

2. Laguerre polynomials (Lαn ), α > −1,

¡ α ¢0
α Ln+1 (x) (Lαn )0 (x)
Ln (x) = +n .
n+1 n
³ ´
(α,β)
3. Jacobi polynomials Jn , α, β > −1,
à !0 !
Ã
(α,β) (α,β) 0
Jn+1 2n(α − β) Jn
Jn(α,β) (x) = (x) + (x)
n+1 (2n + α + β)(2n + α + β + 2)
n
à (α,β) !0
4n(n − 1)(n + α)(n + β) Jn−1
− (x).
(2n + α + β − 1)(2n + α + β)2 (2n + α + β + 1) n − 1
78 J. Hernández and F. Marcellán

4. Bessel polynomials (Bnα ), −α ∈ / N,

µ α ¶0 µ α ¶0
α Bn+1 4n Bn
Bn (x) = (x) + (x)
n+1 (2n + α)(2n + α + 2) n
µ α ¶0
4n(n − 1) Bn−1
+ 2
(x).
(2n + α − 1)(2n + α) (2n + α + 1) n − 1

Finally, we need an expression for xm Pn (x), n, m ∈ N, in terms of (Pn ). If

we use the three-term recurrence formula (2.1) m-times, then it is straight-
forward to verify that
n+m
X
(2.5) m
x Pn (x) = αkn,m Pk (x),
k=n−m

where αkn,m = hu, fm Pn Pk i/hu, Pk2 i and fm (x) = xm . Notice that αn+m
n,m
=1
and
hu, fm−1 Pn f1 Pk i
αkn,m =
hu, Pk2 i
hu, fm−1 Pn (Pk+1 + βk Pk + γk Pk−1 )i
=
hu, Pk2 i
hu, fm−1 Pn Pk+1 i hu, fm−1 Pn Pk i hu, fm−1 Pn Pk−1 i
= 2 + βk 2 + γk
hu, Pk i hu, Pk i hu, Pk2 i
2 i
hu, Pk+1 2 i
hu, Pk−1
n,m−1 n,m−1
= α + β α
k k + γ k αn,m−1 .
hu, Pk2 i k+1 hu, Pk2 i k−1

Now, using again the three-term recurrence formula (2.1), we get γk =

hu, Pk2 i/hu, Pk−1
2 i for k = 1, 2 . . . . Thus,

(2.6) αkn,m = γk+1 αk+1

n,m−1
+ βk αkn,m−1 + αk−1
n,m−1

k,1
and αk−1 = γk , αkk,1 = βk , αk+1
k,1
= 1, for k = 1, 2, . . . .

3. Polynomial Matrices

Let A(s) := sp Ap + sp−1 Ap−1 + · · · + sA1 + A0 be a matrix whose entries

are polynomials of degree at most p. The coefficients (Ak )pk=0 are matrices
in Cn×n .
 An Extension of Leverrier-Faddeev Algorithm . . . 79

Definition 3.1. α ∈ C is said to be a zero of A(s) if det A(α) = 0.

Thus it is an important question to find det A(s) using a finite number

of linear steps as an alternative to the natural approach based in terms
of determinants of polynomial matrices of dimension n − 1. On the other
hand, the adjoint of the matrix A(s) plays an important role in the transfer
function of a linear continuous system. Indeed
1
[A(s)]−1 = Adj A(s).
det A(s)
Notice that deg det A(s) 6 pn and thus the complexity of our problem in-
creases with n.
In the last years an increasing attention was paid to some cases of poly-
nomial matrices.

Definition 3.2. Let µ be a N × N positive definite matrix of measures

supported on the real line, and (Pn )+∞ n=0 a sequence of matrix polynomials
with deg Pn (s) = n. (Pn )+∞
n=0 is orthogonal with respect to µ if
Z
∗
(3.1) Pn (s)dµ(s)Pm (s) = 0, if n 6= m, n, m > 0,
R

and (Pn )+∞

n=0 is orthonormal with respect to µ if
Z
∗
(3.2) Pn (s)dµ(s)Pm (s) = δn,m IN , n, m > 0,
R

For a positive definite matrix of measures non-degenerate µ having mo-

ments of any order, there exists a sequence of orthonormal matrix polynomi-
als with respect to µ. This sequence satisfies a three-term recurrence relation
of the form
(3.3) sPn (s) = An+1 Pn+1 (s) + Bn Pn (s) + A∗n Pn−1 (s)
with the initial conditions P−1 (s) = 0 and P0 (s) = IN , where An is nonsin-
gular and Bn is hermitian, n > 0.
This three-term recurrence relation characterizes sequence of orthonor-
mal matrix polynomials. In other words, if a family of matrix polynomials
satisfies a three-term recurrence relation as (3.3), then this family is or-
thonormal with respect to some positive definite matrix of measures. This
result is known as Favard’s Theorem.
80 J. Hernández and F. Marcellán

Definition 3.3. Let (Pn )+∞ n=0 be a family of orthogonal matrix polynomials
with respect to µ, the matrix polynomial sequence of the second kind is the
+∞
family (Qn )n=0 , defined by
Z
Pn (t) − Pn (x)
(3.4) Qn (t) = dµ(x), n > 0.
t−x

The sequence (Qn )+∞n=0 satisfies the three-term recurrence relation (3.3),
with initial conditions Q0 (t) = 0 and Q1 (t) = A−1
1 .

We get from the three-term recurrence relation that

    
P0 (s) B0 A1 P0 (s)
 P1 (s)   A∗ B1 A2  P1 (s) 
   1  
(3.5) s  P (s)  =  A∗2 B2 A∗3  P2 (s) .
 2    
.. .. .. .. ..
. . . . .

The infinite dimensional matrix

 
B0 A1
 A∗ B1 A2 
 1 
J = A∗2 B2 A∗3 
 
.. .. ..
. . .

is said to be the N -Jacobi matrix associated with the family (Pn )+∞
n=0 . The
N -Jacobi matrix plays an important role in the study of the zeros and
quadrature formula for the matrix polynomials (Pn )+∞
n=0 .

Lemma 3.1. ([3]) For n ∈ N, the zeros of the matrix polynomials Pn (s)
are the same as those of the polynomial det (sInN − JnN ) with the same
multiplicity, where InN is the identity matrix of dimension nN and JnN is
the N -Jacobi matrix truncated of size nN :
 
B0 A1
 A∗ B1 A2 
 1 
 . .. . .. . .. 
(3.6) J = .
 
 An−2 Bn−2 An−1 
∗

A∗n−1 Bn−1

The next result is consequence of the previous lemma.

 An Extension of Leverrier-Faddeev Algorithm . . . 81

Theorem 3.1. ([3]) 1◦ The zeros of Pn have a multiplicity less than or

equal to N . Furthermore Pn has nN real zeros (taking into account their
multiplicities), n ∈ N.
2◦ If a is a zero of Pn with multiplicity p, then rank (Pn (a)) = N − p.
3◦ If we write xn,k (k = 1, . . . , nN ), the zeros of Pn in increasing order
(and taking into account their multiplicities), then

xn+1,k 6 xn,k 6 xn+1,k+N for k = 1, . . . , nN.

4◦ If a is both a zero of Pn and Pn+1 , then Pn (a) and Pn+1 (a) do not
have a common eigenvector associated to zero.
5◦ If xn,k is a zero of Pn with multiplicity N , then Pn (xn,k ) = 0. Fur-
1/N
thermore every complex value of xn,k is a zero of the N consecutive scalar
polynomials pnN +j (t) = det (PnN +j (t)) j = 0, 1, . . . , N − 1. In this case the
real number xn,k can not be a zero of the matrix polynomial Pn+1 .
6◦ Quadrature formula. Associated with every zero xn,k of the matrix
polynomial Pn (k = 1, . . . , nN ) there exists a N × N positive semidefinite
matrix Bn,k of rank one such that
nN
X
hP, Qi(Pn ) = P (xn,k )Bn,k Q∗ (xn,k ),
k=1

for P, Q matrix polynomials satisfying deg(P ) + deg(Q) 6 2n − 1, where

h·, ·i(Pn ) denotes the matrix inner product associated with (Pn )+∞
n=0 .

Let A be a N × N nonsingular matrix, and B be a N × N Hermitian

matrix. We define a sequence of orthonormal matrix polynomials (UnA,B )+∞
n=0 ,
by the three-term recurrence relation
A,B A,B
(3.7) tUnA,B (t) = A∗ Un+1 (t) + BUnA,B (t) + AUn−1 (t), n > 0,
A,B
with the initial conditions U−1 (t) = 0 and U0A,B (t) = IN . The elements of
this family are the Chebyshev matrix polynomials of the second kind.
If A is positive definite, then we can deduce an explicit form for the
corresponding matrix of measures. Indeed, we introduce the matrix
1
KA,B (z) = A−1/2 (B − zIN ) A−1/2 ,
2
82 J. Hernández and F. Marcellán

for z = x ∈ R. This matrix is hermitian and we can express KA,B (x) =

U (x)D(x)U ∗ (x), where D(x) is a diagonal matrix, with di,i (x) ∈ R, i =
1, . . . , N , and U (x)∗ U (x) = IN . Then
1 −1/2
(3.8) dµA,B (x) = A U (x)S(x)U ∗ (x)A−1/2 dx,
π
where S(x) = diag (s1 (x), . . . , sN (x)) and
( q
1 − d2i,i (x) if di,i (x) ∈ (−1, 1),
si (x) =
0 if di,i (x) ∈
/ (−1, 1),
for i = 1, . . . , N . The support of µA,B is
n ³ ´ o
supp (µA,B ) = x ∈ R : σ A−1/2 (xIN − B)A−1/2 ∩ [−2, 2] 6= ∅ .

Thus, we can express the matrix polynomial UnA,B (t) as

µ ¶
A,B −1/2 1 −1/2 −1/2
(3.9) Un (t) = A Un A (B − tIN )A A−1/2 .
2
+∞
where (Un )n=0 is the family of scalar Chebyshev polynomials of the second
kind.

4. The Algorithm

For a matrix A ∈ Cn×n , we presented in [6] an alternative version of

Leverrier algorithm for simultaneous computation of the characteristic poly-
nomial pA (λ) of A and the matrix Adj (λIn − A) in terms of a family of
classical orthogonal polynomials (Pk ), i. e.
n−1
X
(4.1) pA (λ) = Pn (λ) + ân−k Pk (λ),
k=0

n−2
X
(4.2) Ã(λ) = Pn−1 (λ)In + Pk (λ)B̂n−k−1 .
k=0

The algorithm presented in [6] is a follows:

DATA: {βk }n−1 n n−1 n
k=0 , {γk }k=1 , {rk }k=0 , {sk }k=1 .

Initial Condition: B̂−1 = 0, B̂0 = In .

 An Extension of Leverrier-Faddeev Algorithm . . . 83

FOR k = 1, 2, . . . , n − 1
1h ³ ´i
(4.3) âk = (βn−k − rn−k )trB̂k−1 + (γn−k+1 − sn−k+1 )trB̂k−2 − tr AB̂k−1 ,
k
(4.4) B̂k = AB̂k−1 + âk I − γn−k+1 B̂k−2 − βn−k B̂k−1 .

END (FOR)
1h ³ ´i
ân = (β0 − r0 )trB̂n−1 + (γ1 − s1 )trB̂n−2 − tr AB̂n−1 .
n

Now, we are interesting in the computation of a(s) := det (sp Ap + · · · +

sA1 + A0 ) and B(s) := Adj (sp Ap + · · · + sA1 + A0 ). Proceeding in the
same way as [7], that is, in the previous algorithm we replace A by A(s) :=
−(sp Ap + · · · + sA1 + A0 ), then we get
n−1
X
(4.5) ã(λ, s) := det (λIn − A(s)) = Pn (λ) + ân−k (s)Pk (λ),
k=0

as well as
n−2
X
(4.6) B̃(λ, s) := Adj (λIn − A(s)) = Pn−1 (λ)In + Pk (λ)B̂n−k−1 (s).
k=0

Thus, from (4.3) and (4.4) we get

³ ´
kâk (s) = (βn−k − rn−k )trB̂k−1 (s) − tr A(s)B̂k−1 (s)
(4.7)
+(γn−k+1 − sn−k+1 )trB̂k−2 (s), k = 1, . . . , n,

as well as

(4.8) B̂k (s) = âk (s)In − γn−k+1 B̂k−2 (s) − βn−k B̂k−1 (s) + A(s)B̂k−1 (s),

for k = 1, . . . , n − 1.
Thus, if λ = 0 in (4.5) and (4.6) then we get
n−1
X
(4.9) a(s) = ã(0, s) = Pn (0) + ân−k (s)Pk (0),
k=0
n−2
X
(4.10) B(s) = B̃(0, s) = Pn−1 (0)In + Pk (0)B̂n−k−1 (s).
k=0
84 J. Hernández and F. Marcellán

Taking into account deg (Pk (s)) = k for all k > 0, from (4.7) and (4.8)
we can assure that the degrees of the polynomial âk (s), k = 1, 2, . . . , n and
the polynomial matrix B̂k (s), k = 1, 2, . . . , n − 1 are at most equal to kp.
Thus for âk (s) and B̂k (s) we get the expansions
kp
X
âk (s) = ak,j Pj (s), ak,j ∈ C,
j=0
(4.11) kp
X
B̂k (s) = Pj (s)Bk,j , Bk,j ∈ Cn×n .
j=0

Substituting (4.11) in (4.7), we get


kp (k−1)p
X X
k ak,j Pj (s) = tr (βn−k − rn−k ) Pj (s)Bk−1,j
j=0 j=0
(k−2)p
X
+(γn−k+1 − sn−k+1 ) Pj (s)Bk−2,j
j=0 
(k−1)p
X
+(sp Ap + · · · + sA1 + A0 ) Pj (s)Bk−1,j 
j=0

(k−1)p
X
= tr (βn−k − rn−k ) Pj (s)Bk−1,j
j=0
(k−2)p
X
+(γn−k+1 − sn−k+1 ) Pj (s)Bk−2,j
j=0 
p (k−1)p
X X
+ sm Pj (s)Am Bk−1,j  .
m=0 j=0

If we replace (2.5) in the previous expression, we get

kp
X kp
X
k ak,j Pj (s) = αjj−p,p tr (Ap Bk−1,j−p ) Pj (s)
j=0 j=p
kp−1
à 1
!
X X
+ αjj−p+1,p−i tr (Ap−i Bk−1,j−p+1 ) Pj (s)
j=p−1 i=0
..
.
(k−1)p+1
Ãp−1 !
X X
+ αjj−1,p−i tr (Ap−i Bk−1,j−1 ) Pj (s)
j=1 i=0
 An Extension of Leverrier-Faddeev Algorithm . . . 85

(k−1)p
à p
!
X X
+ αjj,p−i tr (Ap−i Bk−1,j ) + (βn−k − rn−k )trBk−1,j Pj (s)
j=0 i=0
(k−1)p−1
Ãp−1 !
X X
+ αjj+1,p−i tr (Ap−i Bk−1,j+1 ) Pj (s)
j=0 i=0
..
.
(k−2)p+1
à 1 !
X X
+ αjj+p−1,p−i tr (Ap−i Bk−1,j+p−1 ) Pj (s)
j=0 i=0
(k−2)p ³ ´
X
+ αjj+p,p tr (Ap Bk−1,j+p ) + (γn−k+1 − sn−k+1 )trBk−2,j Pj (s).
j=0

Identifying coefficients we get

p
X
(4.12) kak,0 = δk trBk−1,0 + ηk trBk−2,0 + α00,p−i tr (Ap−i Bk−1,0 )
i=0
p−1 X
X l
+ α0p−l,p−i tr (Ap−i Bk−1,p−l ) ,
l=0 i=0
p
X l
X
kak,1 = δk trBk−1,1 + ηk trBk−2,1 + α1l−p+1,p−i tr (Ap−i Bk−1,l−p+1 )
l=p−1 i=0
p−1 X
X l
+ α1p−l+1,p−i tr (Ap−i Bk−1,p−l+1 ) ,
l=0 i=0
..
.
p X
X l
l−1,p−i
kak,p−1 = δk trBk−1,p−1 + ηk trBk−2,p−1 + αp−1 tr (Ap−i Bk−1,l−1 )
l=1 i=0
p−1 X
X l
2p−l−1,p−i
+ αp−1 tr (Ap−i Bk−1,2p−l+1 ) ,
l=0 i=0
p X
X l
kak,j = δk trBk−1,j + ηk trBk−2,j + αjj−p+l,p−i tr (Ap−i Bk−1,j−p+l )
l=0 i=0
p−1 X
X l
j+p−l,p−i
+ αp−1 tr (Ap−i Bk−1,j+p−l ) ,
l=0 i=0

for j = p, . . . , (k − 2)p, and

86 J. Hernández and F. Marcellán

kak,(k−2)p+1 = δk trBk−1,(k−2)p+1
p X
X l
(k−3)p+l+1,p−i ¡ ¢
+ α(k−2)p+1 tr Ap−i Bk−1,(k−3)p+l+1
l=0 i=0
p−1 X
X l
(k−1)p−l+1,p−i ¡ ¢
+ α(k−2)p+1 tr Ap−i Bk−1,(k−1)p−l+1 ,
l=1 i=0
..
.
kak,(k−1)p−1 = δk trBk−1,(k−1)p−1
p X
X l
(k−2)p+l−1,p−i ¡ ¢
+ α(k−1)p−1 tr Ap−i Bk−1,(k−2)p+l−1
l=0 i=0
p−1
X (k−1)p,p−i ¡ ¢
+ α(k−1)p−1 tr Ap−i Bk−1,(k−1)p ,
i=0
p X
X l
(k−2)p+l,p−i ¡ ¢
kak,(k−1)p = δk trBk−1,(k−1)p + α(k−1)p tr Ap−i Bk−1,(k−2)p+l ,
l=0 i=0
p X
X l
(k−2)p+l+1,p−i ¡ ¢
kak,(k−1)p+1 = α(k−1)p+1 tr Ap−i Bk−1,(k−2)p+l+1 ,
l=0 i=0
..
.
1 X
X l
(k−1)p+l−1,p−i ¡ ¢
kak,kp−1 = αkp−1 tr Ap−i Bk−1,(k−1)p+l−1 ,
l=0 i=0
(k−1)p,p ¡ ¢
kak,kp = αkp tr Ap Bk−1,(k−1)p ,

where δk = βn−k − rn−k , ηk = γn−k+1 − sn−k+1 .

In an analogous way, substituting (4.11) in (4.8), we deduce
kp
X kp
X p (k−1)p
X X
Pj (s)Bk,j = ak,j Pj (s)In − sm Pj (s)Am Bk−1,j
j=0 j=0 m=0 j=0
(k−1)p (k−2)p
X X
− βn−k Pj (s)Bk−1,j − γn−k+1 Pj (s)Bk−2,j .
j=0 j=0

Replacing (2.5) in the previous expression, we get

kp
X kp
X kp
X
Pj (s)Bk,j = ak,j Pj (s)In − Pj (s)αjj−p,p Ap Bk−1,j−p
j=0 j=0 j=p
 An Extension of Leverrier-Faddeev Algorithm . . . 87
kp−1
à 1
!
X X
− Pj (s) αjj−p+1,p−i Ap−i Bk−1,j−p+1
j=p−1 i=0
..
.
(k−1)p+1
Ãp−1 !
X X
− Pj (s) αjj−1,p−i Ap−i Bk−1,j−1
j=1 i=0
(k−1)p
à p
!
X X
− Pj (s) βn−k In + αjj,p−i Ap−i Bk−1,j
j=0 i=0
(k−1)p−1
Ãp−1 !
X X
− Pj (s) αjj+1,p−i Ap−i Bk−1,j+1
j=0 i=0
..
.
(k−2)p+1
à 1
!
X X
− Pj (s) αjj+p−1,p−i Ap−i Bk−1,j+p−1
j=0 i=0
(k−2)p
X ³ ´
− Pj (s) αjj+p,p Ap Bk−1,j+p + γn−k+1 Bk−2,j .
j=0

Identifying coefficients we obtain

(4.13) Bk,0 = ak,0 In − βn−k Bk−1,0 − γn−k+1 Bk−2,0
p−1 X
X l p
X
− α0p−l,p−i Ap−i Bk−1,p−l − α00,p−i Ap−i Bk−1,0 ,
l=0 i=0 i=0
Bk,1 = ak,1 In − βn−k Bk−1,1 − γn−k+1 Bk−2,1
p
X l
X
− α1l−p+1,p−i Ap−i Bk−1,l−p+1
l=p−1 i=0
p−1 X
X l
− α1p−l+1,p−i Ap−i Bk−1,p−l+1 ,
l=0 i=0
..
.
Bk,p−1 = ak,p−1 In − βn−k Bk−1,p−1 − γn−k+1 Bk−2,p−1
p X
X l
l−1,p−i
− αp−1 Ap−i Bk−1,l−1
l=1 i=0
p−1 X
X l
2p−l−1,p−i
− αp−1 Ap−i Bk−1,2p−l+1 ,
l=0 i=0
Bk,j = ak,j In − βn−k Bk−1,j − γn−k+1 Bk−2,j
88 J. Hernández and F. Marcellán
p X
X l
− αjj−p+l,p−i Ap−i Bk−1,j−p+l
l=0 i=0
p−1 X
X l
j+p−l,p−i
− αp−1 Ap−i Bk−1,j+p−l ,
l=0 i=0

for j = p, . . . , (k − 2)p, and

Bk,(k−2)p+1 = ak,(k−2)p+1 In − βn−k Bk−1,(k−2)p+1
p X
X l
(k−3)p+l+1,p−i
− α(k−2)p+1 Ap−i Bk−1,(k−3)p+l+1
l=0 i=0
p−1 X
X l
(k−1)p−l+1,p−i
− α(k−2)p+1 Ap−i Bk−1,(k−1)p−l+1 ,
l=1 i=0
..
.
Bk,(k−1)p−1 = ak,(k−1)p−1 In − βn−k Bk−1,(k−1)p−1
p X
X l
(k−2)p+l−1,p−i
− α(k−1)p−1 Ap−i Bk−1,(k−2)p+l−1
l=0 i=0
p−1
X (k−1)p,p−i
− α(k−1)p−1 Ap−i Bk−1,(k−1)p ,
i=0
Bk,(k−1)p = ak,(k−1)p In − βn−k Bk−1,(k−1)p
p X
X l
(k−2)p+l,p−i
− α(k−1)p Ap−i Bk−1,(k−2)p+l ,
l=0 i=0
Bk,(k−1)p+1 = ak,(k−1)p+1 In
p X
X l
(k−2)p+l+1,p−i
− α(k−1)p+1 Ap−i Bk−1,(k−2)p+l+1 ,
l=0 i=0
..
.
1 X
X l
(k−1)p+l−1,p−i
Bk,kp−1 = ak,kp−1 In − αkp−1 Ap−i Bk−1,(k−1)p+l−1 ,
l=0 i=0
(k−1)p,p
Bk,kp = ak,kp In − αkp Ap Bk−1,(k−1)p .

Thus, the algorithm can be described as follows

DATA: {βk }n−1 n n−1 n
k=0 , {γk }k=1 , {rk }k=0 , {sk }k=1 .

Initial Condition: â0 (s) = 1, B̂−1 (s) = 0, B̂0 (s) = In , and Bi,j = 0 if i, j < 0 or
i < pj.
 An Extension of Leverrier-Faddeev Algorithm . . . 89

FOR k = 1, . . . , n − 1

• From B̂k−2 (s) and B̂k−1 (s) and taking into account (4.12) we get âk (s).
• From (4.13) we get B̂k (s).

END (FOR)

From B̂n−2 (s) and B̂n−1 (s) taking into account (4.12) we get ân (s).

5. Examples

Example 5.1. Consider the polynomial matrix

 2 
s − s −s2 + s − 1 0
A(s) =  0 3 s2 + 2s  ,
s2 − 1 −s2 s2 + s

so that p = 2, n = 3 and
     
1 −1 0 −1 1 0 0 −1 0
A2 =  0 0 1  , A1 =  0 0 2  , A0 =  0 3 0 .
1 −1 1 0 0 1 −1 0 0

It is easy to verify that det A(s) = 3s4 − s3 − 4s2 + 2s,

 4 
s + 2s3 + 3s2 + 3s s4 + s −s4 − s3 + s2 − 2s
Adj A(s) =  s4 + 2s3 − s2 − 2s s4 − s2 −s4 − s3 + 2s2  .
2
s −1 −s2 s2 + s

We apply the above algorithm using the Chebyshev basis (Chebyshev

polynomials of first kind), (Tn )+∞
n=0 . The corresponding data are βn = rn =
0, n > 0 and γ1 = 1/2, s1 = 0, γn = 1/4, sn = −1/4, n > 2.
From (4.12) and (4.13) we get the following results.
For k = 1 :
a1,0 = 4, a1,1 = 0, a1,2 = 2;
     
7 3 0 1 −1 0 1 1 0
1
B1,0 =  0 2 −1  , B1,2 = 0 0 −2  , B1,2 =  0 2 −1  .
2
1 1 7 0 0 −1 −1 1 1
90 J. Hernández and F. Marcellán

For k = 2 :
3
a2,0 = 4, a2,1 = , a2,2 = 7, a2,3 = 2, a2,4 = 2;
2
   
19 3 1 18 4 −11
1 1
B2,0 =  −1 3 5  , B2,1 =  −2 0 −3  ,
8 4
12 8 16 0 −4 −12
     
4 1 0 2 0 −1 1 1 −1
B2,2 =  0 0 1  , B2,3 =  2 0 −1  , B2,4 =  1 1 −1  .
−3 0 3 0 0 0 0 0 0

For k = 3 :
9 5
a3,0 = , a3,1 = , a3,2 = 0, a3,3 = −1, a3,4 = 3, a3,5 = 0 = a3,6 .
8 4
Hence, from (4.9) we get
5 7
det A(s) = 3T4 (s) − T3 (s) − T2 (s) + T1 (s) − T0 (s),
4 8
and, from (4.10) we obtain
    
1 1 −1 2 0 −1 4 1 0
Adj A(s) = T4 (s)  1 1 −1  + T3 (s)  2 0 −1  + T2 (s)  0 0 1 
0 0 0 0 0 0 −3 0 3
   
18 4 −11 15 3 1
1 1
+T1 (s)  −2 0 −3  + T0 (s)  −1 −1 5  .
4 8
0 −4 −12 12 8 12

Example 5.2. Consider the matrices

" # " #
2 −1 −1 1
A= and B= .
−1 1 1 0
Taking into account the three-term recurrence relation (3.7), we get
" #
A,B
13t4 − 32t3 + 24t2 − 6t 21t4 − 52t3 + 39t2 − 8t
(5.1) U4 (t) = .
21t4 − 52t3 + 39t2 − 8t 34t4 − 84t3 + 63t2 − 14t
We need to compute det (U4A,B (t)) to find the zeros of U4A,B (t). For the
matrix polynomial U4A,B (t), p = 4, n = 2 and
" # "# " #
0 0 6 8 24 39
A0 = , A1 = − , A2 = ,
0 0 8 14 39 63
" # " #
32 52 13 21
A3 = − , A4 = .
52 84 21 34
 An Extension of Leverrier-Faddeev Algorithm . . . 91

We apply the algorithm of the previous section, using the Chebyshev basis
(Chebyshev polynomials of second kind), (Un )+∞
n=0 . The corresponding data
are βn = rn = 0, n > 0 and γn = 1/4, n > 1 and s1 = 0, sn = −1/4, n > 2.
From (4.12) and (4.13) we get:
For k = 1 :
221 489
a1,0 = , a1,1 = −78, a1,2 = , a1,3 = −116, a1,4 = 47,
8 4

and
" # " # " #
1 160 −99 −56 34 1 354 −219
B1,0 = , B1,1 = , B1,2 = ,
8 −99 61 34 −22 4 −219 135
" # " #
−84 52 34 −21
B1,3 = , B1,4 = .
52 −32 −21 13

For k = 2 :
2265 209 225 177
a2,0 = , a2,1 = − , a2,2 = 93, a2,3 = − , a2,4 = ,
128 4 2 2
69
a2,5 = −20, a2,6 = − , a2,7 = 4, a2,8 = 1.
4

Thus,

8
X 1
â2 (t) = a2,j Uj (t) = t8 + 4t7 − 19t6 − 26t5 + 111t4 − 90t3 + 20t2 + .
4
j=0

Notice that is not necessary to compute â1 (t) because U1 (0) = 0. Then
from (4.9), the determinant of U4A,B (t) is

det U4A,B (t) = U2 (0)+U0 (0)â2 (t) = t8 +4t7 −19t6 −26t5 +111t4 −90t3 +20t2 .

Acknowledgements. The work of the second author (F. Marcellán)

was supported by Dirección General de Investigación, Ministerio de Ciencia
y Tecnologı́a of Spain, under grant BFM2003-06335-C03-02. The work of the
first author (J. Hernández) has been supported by Fundación Universidad
Carlos III de Madrid.
92 J. Hernández and F. Marcellán

REFERENCES

1. S. Barnett: Leverrier’s algorithm: A new proof and extensions. SIAM

J. Matrix Anal. Appl. 10 (1989), 551–556.
2. T. S. Chihara: An Introduction to Orthogonal Polynomials. Gordon and
Breach, New York, 1978.
3. A. J. Duran and P. Lopez-Rodriguez: Orthogonal matrix polynomials:
Zeros and Blumenthal’s theorem. J. Approx. Theory 84 (1996), 96–118.
4. E. Emre and O. Hüseying: Generalization of Leverrier’s algorithm to
polynomial matrices of arbitrary degree. IEEE Trans. Automat. Control,
February 1975, 136.
5. G. M. L. Gladwell: Inverse Problems in Vibration. Kluwer Academic
Publishers, Dordrecht, 1986.
6. J. Hernández, F. Marcellán, and C. Rodrı́guez: Leverrier-Fadeev
algorithm and classical orthogonal polynomials. Rev. Acad. Colomb. Cienc.
28 (106) (2004), 39–47.
7. J. Hernández and F. Marcellán: Classical orthogonal polynomials and
Leverrier’s algorithm for singular pencils (submitted).
8. F. Marcellán, A. Branquinho, and J. C. Petronilho: Classical or-
thogonal polynomials: a functional approach. Acta Appl. Math. 34 (1994),
283–303.
9. B. G. Mertzios: Leverrier’s algorithm for singular systems. IEEE Trans.
Automat. Control AC-29 (1984), 652–653.
10. G. Wang and Y. Lin: A new extension of Leverrier’s algorithm. Linear
Algebra Appl. 180 (1993), 227–238.

Departamento de Matemáticas
Universidad Carlos III de Madrid
Avenida de la Universidad 30
28911, Leganés, Spain
e-mails: [email protected]
[email protected]