Численное моделирование внутрикамерных нестационарных турбулентных течений. Часть 2

Let A(t) be a matrix function from the matrix Wiener algebra W ^pxp (T) that is invertible on the unit circle T. The representation

A(t) = A+ (t)D(t)A - (t), t G T,                             (1)

is called a left Wiener-Hopf factorization of A(t). Here A ± (t) belong to the group GW ±^xp (T) of invertible elements of the subalgebra W ±^xp (T), the middle factor D(t) is the diagonal matrix D(t) = diag [t ^{A 1} ,... ,t ^Ap] , where integers A 1 >  ... >  A _p are the left partial indices of A(t). The relation A 1 + ... + A p = к = ind _T det A(z) is valid. A similar representation in which the factors A ± are rearranged is called the right Wiener-Hopf factorization .

Mathematical modelling of wave diffraction, problems of dynamic elasticity and fracture mechanics, and geophysical problems are often reduced to the Wiener–Hopf factorization problem for matrix functions [1–4]. The factorization of matrix functions is also a powerful tool itself used in various areas of mathematics [5–7, 9].

Unfortunately, for the matrix case, there is no constrictive solution of the factorization problem in a general setting and it is very important to find cases when the problem can be solved effectively or explicitly. By the explicit (or constructive ) solution of the factorization problem we understand a clearly defined algorithmic procedure that should definitely terminate after a finite number of steps. There are not that many classes of matrix functions for which an explicit solution to factorization problem has been found.

The most important of them are classes of matrix polynomials [10, 11] and meromorphic matrix functions [12]. A detailed review of constructive methods for the factorization problem is presented in the works [13–15].

In addition to the aforementioned lack of availability of explicit solution to the factorization problem, in the general case, there is another obstacle to use the technique. This is possible instability of the factorization problem. Even if an explicit method for solving the particular factorization problem exists, each step of the respective algorithm can be executed exactly or approximately (numerically).

We say that the problem can be solved exactly if (i) the input data belonging to the Gaussian field Q(i) of complex rational numbers, and (ii) all steps of the explicit algorithm can be perform in the exact arithmetic . The instability of the problem leads to the fact that the explicit algorithm cannot be implemented numerically. As a rule, researchers developing a particular explicit factorization method usually ignore this issue. In fact, they implicitly assume that all steps of the proposed explicit algorithm can be carried out exactly that, unfortunately, is not always possible.

For the first time, the need to accurate study the way of numerical implementation of the explicit algorithm was highlighted in [16]. This has been done for matrix polynomials, where existence of an explicit solution of the factorization problem was proved in [11]. In [16], based on this work, a criterion for the exact factorizability of a matrix polynomial was obtained, and an exact algorithm for a solution of the factorization problem was developed. This algorithm was also implemented as the package ExactMPF in Maple. Thus, if the condition satisfies, the problem of an instability does not arise.

The package makes it easy to carry out numerical experiments with the Wiener–Hopf factorization for matrix polynomials. It can be used to construct an approximate canonical factorization with quaranteed accuracy for strictly nonsingular 2 x 2 matrix functions and to the integration of a discrete analog of the nonlinear Schro¨dinger equation by Inverse Scattering Transform. We hope that the application of the package will not be exhausted by these examples.

This paper is complimentary to [16], where the length was limited by the publisher rules. As a result, some crucial technical results have been omitted there. In particular, the algorithm for constructing essential polynomials was not described. In this work we fill this gap.

1.    Explicit Solution of the Factorization Problems for MatrixPolynomials

In this section, we present an explicit algorithm for the factorization of an arbitrary matrix polynomial. Our presentation is based on the results from [11, 12, 16].

N

Supposed that the matrix polynomial a(z) = £ a k z k , a k E C ^pxp , is invertible on the k =0

unit circle T. We will write its left and right Wiener–Hopf factorizations in the form

a(t) = l₊(t)d _L (t)l - (t), a(t) = r - (t)d _R (t)r₊ (t), t E T. (2)

Here d _L (t) = diag[t ^A 1 ,...,t ^{X p} ], and A 1 > ••• >  A _p ; d _R (t) = diag [t ^p 1 ,...,t ^{P p} ], where p1 <  ... <   p p . Note that left A j- and right p j indices are usually different sets of integers and constructions of the right and left factorizations are usually considered as two separate problems. For explicit construction of these factorizations we will use the method proposed in the work [11]. The method requires simultaneous considerations of the both factorizations.

Let A(z) = det a(z) and A(z) = A - (z)z ^к A ₊ (z), A _— ( to ) = 1, be the Wiener-Hopf factorization of A(z). The factorization is unique with the additional condition at infinity for the polynomial A - (z). In the sequel, we use, in fact, only one of the factors, namely, A _ (z) = 1 + A — z ^-1 + • • • + A — z ^-к , к = ind _T det a(z).

We expand the rational matrix function A -¹ (z)a(z) in the Laurent series at infinity: A -1 (z)a(z) = ^ N-^ C j z ^j . The coefficients C j G C ^pxp are computed recurrently in terms of matrix coefficients a j of the original matrix polynomial, a(z), and the coefficients A _j ^- of the scalar polynomial, A (z), (see [16, Eq.(2.5)]).

To construct the factorizations of a(z) we only need a finite number of the coefficients, c k , for k = -к, ..., 0,..., к. Denote с ^к _к := {с _-к , ... ,c₀,..., с _к } . The main tools for computations of the partial indices and factors in the factorizations of matrix functions are the so-called indices and essential polynomials of the sequence с — _к (see [11,17]). Let us define these notions.

Form a finite family of the block Toeplitz matrices of finite sizes:

Tk ||ci—j ^i=k,k+1,...,K , к A k A к, (3) j=0,1,...,K+k and study the structure of the right kerR Tk = {R G С(к+к+1)х1 |TkR = 0} and left kerL Tk = {L G C1x(K-k+1) lLTk = 0} kernels of Tk. Further it is more convenient to deal not with vectors R = (r0,r1,... ,гк+к)T G kerRTk, rj G Cpx1, but with their generating column-valued polynomials R(z) = r0 + r1 z + • • • + r k+K zk+к. We will use the spaces Nk of the generating polynomials instead of the spaces ker Tk .

By N _k ^R , - κ ≤ k ≤ κ, we denote the space of generating vector polynomials for vectors in ker R T k . Put N ^R _{K- 1} := { 0 } and let N R ₊₁ be (2к + 2)p-dimensional space of all column-valued polynomials whose degrees are not greater than 2к + 1.

Repeating the same line of reasoning, we can define the space N _k ^L , - κ ≤ k ≤ κ, of the row-valued generating polynomials in z ^-1 for the rows from ker _L T _k .

By d _k ^R , we denote a dimension of the right kernel N _k ^R and introduce the following integers: A R = d ^R — d ^R _{— 1} for —к A k A к +1. A sequence с — _к is called regular if A ^R _K = 0 and A ^R ₊₁ = 2p.

For a regular sequence, we have (see [11, 17])

R R     RR

^U     —к    — к +1      - к   к +1     ^p .

Since a monotone integer sequence is piecewise constant, then there are 2p integers Ц1 A ■■■ A M 2 p such that

R — к	= ••• = A R 1   =	0 ,
R A M i +1	= • • • = A R = µ i+1	i,
R A M 2p + 1	= • • • = A R+1 =	2 p.

The absence of the j -th row here means that M j +i = M j •

Definition 1. The integers Ц1, ... , ц_{2 р} defined by the relations (4) are the indices of the sequence с — _к .

Similarly, we can consider the sequence of the left kernel N _k ^L that will lead, however, to the same indices.

Furthermore, we define the right essential polynomials of the sequence с — _к . Note that N _k ^R and z N _k ^R are subspaces of N _k ^R ₊₁ as it follows from the definition of the spaces N _k ^R . The dimension, h ^R ₊₁ , of the complement H ^R ₊₁ of N R + zN R in N _k ^R ₊₁ is equal to A _k ^R ₊₁ — A _k ^R .

Then, Eqs. (4) imply that h ^R ₊₁ = 0 if and only if k = ^ j , j = 1,..., 2p. Moreover, in this case, h ^R ₊₁ is equal to the multiplicity, K j , of the index ^ j . Therefore, for k = ^ j we have

N R +1 = N R N , and for k = ^ j

N R +1 = N, +zN R ) ®H R +_V

Definition 2. Any polynomials R j (z), ..., R j _{+ K j} -1 (z) forming a basis for a complement H R +1 are called right essential polynomials of the sequence c - _K corresponding to the index µ j .

As a result, we have defined 2p indices ^1, ... , ^_{2 p} and 2p right essential polynomials R1(z),..., R_{2 p} (z) for any regular sequence c - _K . Similarly, we can define the left essential polynomials L1 (z),..., L_{2 p} (z) of the sequence c - _K .

In factorization problems, there are natural candidates for the role of indices and essential polynomials. To check that this is indeed the case, the following essentiality criterion can be used ( [17, Theorem 4.1], see also [11, p. 258],):

Theorem 1. The integers Ц1,...,ц2р are the indices and R1(z),..., R2p (z) essential polynomials of the regular sequence c-κκ if and only if the matrix are right

л ( aR {z ^-K-1 R1(z)} ••• a R    ^K R2 p (z)}

)

^A R =\       R1 (0)        •••         R2 _P (0)

is invertible. Here a R {z ^-K-1Rj (z)} = ^ ^ =1 ⁺¹ c _{K +1 -m} r m -

By Theorem 3.1 from [11], the sequence c _- ^κ _κ is regular and there exist respective essential polynomials R1 ( z ) , . . . , R2 p ( z ); L1 ( z ) , . . . , L2 p ( z ) such that

• (i)    the constant terms of the polynomials R1 ( z ) , . . . , R p ( z ) are equal to zero,

• (ii)    the leading terms of the polynomials L p +1 ( z ) , . . . , L2 p ( z ) are equal to zero.

Definition 3. The essential polynomials R1(z),..., R_{2 p} (z); L1 (z),..., L_{2 p} (z) satisfying the conditions (i), (ii) are called the factorization essential polynomials of the sequence.

Now we can formulate a final result on the explicit Wiener–Hopf factorization of a matrix polynomial a(z).

Theorem 2. [11, Theorem 3.2] Let ^1,...,^_{2 p} be the indices and R1 (z),..., R_{2 p} (z) L1 ( z ) , . . . , L2 p ( z ) are the right (left) factorization essential polynomials of the regular sequence c _- ^κ _κ . Let us introduce the p × p matrix functions

/^Lp +1 (z) \

R1(z) = (R1(z) ... R p (z)) , L2(z)= I . I

\ L2p (z) / and dL(z) = diag[z-11,..., z 1p], dR(z) = diag[z1p+1,... ,z12].

Then the left (X1 > • • • >  X _p ) and right (p1 < • • • <  p p ) partial indices and the factors (l ± (z), r ± (z)) of the respective factorizations of the matrix polynomial a(z) are defined by the formulas

X1 = —P1, . . . , Xp = —^p, P1 = Pp+1, . . . , Pp = ^2p, l-(z) = zK+1A—(z)d-1 (z)R-1 (z), l+ (t) = z-K-1 A—1 (z)a(z)R1 (z),(6)

r- (z) = A-(z)L2 1 (z), r+(z) = A-1(z)dR1 (z)L2(z)a(z).

In the statement of this theorem, we have corrected the misprints appeared in the formulas for the factors l₊ (z), r₊ (z) in [11, Theorem 3.2].

Let us list the basic steps of the presented factorization algorithm.

• 1.    Calculation of the Laurent coefficients Cj, —к < j< к, for the rational matrix functions A-1(z)a(z).

• 2.    Calculation of the indices for the sequence c^κ_κ.

• 3.    Calculation of the essential polynomials for the sequence c^κ_κ .

• 4.    Constructing the factorizations in accordance with Th.2.

2.    Exact Solution of the Factorization Problems for MatrixPolynomials

Here к = ind T det a(z) is a number of zeros of det a(z) in open disc | z | < 1.

Finding к and constructing the factorization of scalar polynomial A(z) can be considered as explicit procedures. Now calculation of the Laurent coefficient c j using recurrence relations requires a finite number of operations.

To calculate the indices µ 1 , . . . , µ 2p it is needed to find ranks of the matrices T k , -к <  j <  к. We can do it by means of linear algebra in a finite number of steps.

For this it is necessary to find bases of the ker R,_L T k , —к <  j <  к. We can do it by means of linear algebra in a finite number of steps.

Now this step can be done in an explicit form.

Thus, in accordance with our understanding of the explicit solution of the factorization problem given above, the presented algorithm indeed belongs to this class.

• 1.    The factorization of scalar polynomial A(z), in general case, can only be constructed approximately.

• 2.    Finding the indices and essential polynomials of the sequence c ^κ _κ requires calculating ranks and constructing bases of kernels for matrices T k . Unfortunately, those operations can be unstable.

Thus, in general, the proposed explicit factorization algorithm can not be implemented numerically.

Remark 1. A numerical implementation of the algorithm proposed in [10] meets into the same difficulties.

However, there is still a possibility to implement the algorithm exactly by utilising calculations in rational arithmetic. Obviously, we must demand that the coefficients a j of the original matrix polynomial a(z) must belong the Gaussian field Q(i) and the factorization of A(z) should be performed exactly. In this case the calculations of the Laurent coefficients c j and finding the indices µ ₁ , . . . , µ _2p can also be made exactly.

Now we have to make sure that finding the factorization essential polynomials can also be performed exactly . This was not done in [16] and it is the main goal of this work.

In the following theorem we describe the algorithm of finding these essential polynomials and prove that this algorithm can be implemented in the exact arithmetic if entries of the matrices C j , —к <  j <  к, belong to the field Q(i).

Theorem 3. Let скк := {с-к,... ,c0,..., ск} be a regular sequence of complex p x p matrices with entries from the field Q(i). Suppose that the indices Ц1,...,Ц2Р of the sequence satisfy the condition £p=1 ^j = -к, ^2=p+1 ^j = к, and the sequence has the factorization essential polynomials. Then these polynomials can be found by calculations in the exact arithmetic.

Proof. Let us restrict ourselves to considering only the right essential polynomials R1(z),..., R_{2 p} (z). Algorithm of finding the factorization essential polynomials is based on the criterion of essentiality (Th. 1).

By Definition 3, these polynomials have zero constant terms: R1 (0) = • • • = R _p (0) = 0. Hence, to construct the factorization essential polynomials R1(z),..., R _p (z) we must select first p vector polynomials R j (z) E N R j ₊₁ , j = 1,... ,p, such that the 2p x 2p matrix

A r =

( a R a "  R(z)}

{       0

••• S r {z- ^- Rp (z) }

•••           0

a R { z -^K-1 R _p+1 (z)}

R + 0'

• •• aR        R_{2 p} (z)}

~     ^R2 p ⁽⁰⁾

)

will be invertible or, in other words, the p x p submatrices

Л 11 = (a_R{z ^—K—1 R1 (z) }--- a_R ' R p (z) }) , Л ₂₂ = (R +1 (0) ••• R _2p (0))

will be invertible.

Now it will be convenient to introduce the distinct indices and to assign them the respective multiplicities. Moreover, it is necessary to highlight the border index µ p . Let v1 <  • • • < v s be the distinct indices of the sequence c — _K and K 1 ,..., K _s their multiplicities. Let the index µ _p coincides with ν _t .

We will use induction by a number of indices v 1 ,...,v _t . First we select the factorization essential polynomials R1 (z),... ,R _{K 1} (z) corresponding to the first index v 1 of K 1 -multiplicity. They are the generating polynomials of vectors forming a basis of ker T _v 1 ₊₁ . Since there are the factorization essential polynomials, there exist a basis of N _{V 1 +1} = ker T _{V 1 +1} such that R 1 (z) = zR1 (z), ... ,R _{K 1} (z) = zR _{K 1} (z) for some polynonials R 1 (z),... ,R _{K 1} (z) from the space N _{V 1 +1} = ker T _v 1 ₊₁ , where the matrix T _{V 1 +1} is obtained from T _{V 1 +1} by deleting of the first p columns, i.e. by deleting the first block column of the matrix. It is easily seen that corresponding vectors R1,..., R _K 1 form a basis ker T _{V 1 +1} . Thus, by virtue the essentiality criterion (Th.1) there exists a basis of ker T _{v 1 +1} such that the 2p x K1 submatrix

/ aR {z ^к 1R1 (z) }

a R {z ^к 1R_K1 (z)} A ⁰

of Л R has the rank is equal to к 1 . In fact, it is easy to show that this condition is fulfilled for any choice of a basis R1,..., R _{K 1} . Since the entries of the matrices C j belong to Q(i), the construction of this basis and calculation of the rank can be done in the exact arithmetic.

Thus, we can exactly construct the first K1 polynomials R1 (z),... ,R _{K 1} (z) such that R 1 (0) = 0,..., R _{K 1} (0) = 0, entries of the p x K1 matrix

⁽a R {z ^—K—1 R1(z)} • • • a R { z ^—к—1 R K 1 (z) } )

belong to Q(i), and this matrix has the rank equal to K1 .

Now we repeat these considerations for the other indices v2 , . . . , v t . Assume first that ^ _p <  M _p+1 . Recall that v _t = ^ _p and has the multiplicity K t . In this case K1 + • • • + K t coincides with the number of the indices Ц1,..., ^ _p , that is K1 + • • • + K t = p.

Suppose that we construct the polynomials

R1 (z),..., R k 1 (z); R k 1 +1 (z),..., R k 1 + k 2 (z);...; R K 1 +^+ K j - 1 +1 (z),..., R^^^ (z)

corresponding to the indices v₁,...,V j , 2 <  j <  v _t-1 , such that R ₁ (0) = 0,..., R k _{1 +} .. • _{+ K j} (0) = 0, entries of the p x ( k ₁ + • • • + K j ) matrix

(S r {z ^-K- R1 (z) }... S r       R .       (z) } )

belong to Q(i), and this matrix has the rank equal to k ₁ + • • • + K j .

Let us define the polynomials R _{K 1} . . _{K j +1}(z), ..., R _{K 1} . . _{K j+1} (z) corresponding to the index V j ₊i of the multiplicity K j₊₁ . These polynomials belong to the space N _{V j+1 +1} = ker T _{V j+1 +1} . Let n j ₊₁ is the dimension of the space N _{V j+1 +1} = ker T _{V j+1 +1} and the polynomials Q ₁ (z),..., Q n _j+1 (z) be a basis of this space. Here T _{V j+1 +1} is obtained from T V j +1 +1 by deleting the first block column of the matrix. Hence, zQ ₁ (z),... ,zQ n _j+1 (z) is a basis of the space N _{V j+1 +1} .

has the rank equal to k ₁ + • • • + K j + 1. This selection is always possible since the sequence c ^κ _{- κ} has the factorization essential polynomials. In a similar way, we select the other polynomials R _{K 1 +} ... _{+ K j +2} (z),..., R _{K 1} +-+ _{K j+1} (z) for which the matrix

(S r { z ^-к-1 R 1 (z) }...             R .     .      (z) } )

has the rank equal to k ₁ + • • • + K j₊₁ .

Hence, in the case of k ₁ + • • • + K t = p , we obtain, by induction, the polynomials R₁(z),..., R _p (z), for which the matrix

Л 11 = (S r { z ^-K-1 R 1 (z) }--- S r { z ^{- k - 1} R p ( z ) } )

over Q(i) has the rank equal to p. Thus, in this case, the first p factorization essential polynomials R ₁ (z),..., R _p (z) are exactly constructed.

Now, we build the polynomials R _p+1 (z),..., R _2p (z) such that the matrix Л ₂₂ is invertible. These polynomials must be sequentially chosen from the spaces N j ₊₁ = ker T j₊₁ , j = v _t+1 ,..., v s . It is clear that we can always choose polynomials R _p+1 (z),..., R _{p + K} t ₊₁ (z) from the basis N _{V t+1 +1} , such that vectors R _p+1 (0),... ,R _{p + K} t ₊₁ (0) are linear independent. Otherwise, the sequence c ^κ _{- κ} would not have factorization essential polynomials. Repeating these arguments for the indices V _{t .2} . . . , V _s we arrive to polynomials R _{p .1} ( z ) , . . . , R _{2 p} ( z ) for which the matrix Л ₂₂ is invertible. Therefore, in the case of R _p < R _{p +1}, the right factorization essential polynomials can always be found by the exact computation.

Let us consider now the case when the border index R p satisfies the equality R p = ^ _p+1 , or more precisely, when

R1 < . . . < Rp-l   Vt-1 < Rp-l+1    • • • Rp • • • Rp+m   vt < Rp+m+1 < . . . < R2p for some l > 0, m > 0. Then k1 + • • • + Kt-1 = p — l, Kt = l + m, Kt+1 + • • • + Ks = p — m and k1 + • • • + Kt = p + m > p. The right factorization polynomials R1 (z),..., R.^.- (z) corresponding to the indices v1,..., vt-1 we can construct as above. Recall that Kt is the number of the right essential polynomials RK1 .•••Kt_1+1(z),..., RK1+^Kt (z) corresponding to the index vt. They belong to the space NVt+1 = ker TVt+1. These polynomials are divided into two type. For the first l polynomials RK1+^Kt_1+1(z),..., Rp(z), the conditions

R _K 1 +-_{K t-1 +}i (0) = 0,..., R p (0) = 0 and the invertibility of the matrix Лц must be fulfilled. The remaining m polynomials R _{p +1}(z),... ,R _{K 1 +} + _K t (z) E N _v t ₊₁ must be chosen so that the vectors R _p+1 (0),..., R _{K 1 +} + _K t (0) are linearly independent.

The first type polynomials we can construct as above by choosing successively l polynomials from a basis Q ₁ (z), ... ,Q n _{t +1} of the space N _{V t +1} = ker T _{V t +1} . The existence of the factorization essential polynomials guarantees that this process can be carried out.

The remaining m polynomials R _{p +1}(z),..., R _{K 1 +} + _{K t} (z) must be chosen from the elements of a basis of the space N _{P p +1} = ker T _{V t +1} in a way that the rank of the matrix (R _p+1 (0) • • • R _{K 1 +}_____ + _{K t} (0)) is equal to m = k ₁ + • • • + K t — p. It is again possible since the sequence c - _K possesses the factorization essential polynomials.

By repeating this choice for the spaces N _{V j +1} = ker T _{V j +1} , j = v _{t +1} , ... ,v s , we obtain the polynomials R _{p +1}(z),..., R_{2 p} (z) for which the matrix Л₂₂ is invertible. Then, for the polynomials R₁ (z),...,R _p (z), R _{p +1}(z),..., R_{2 p} (z), the matrix Л R is invertible and these polynomials are the right factorization essential polynomials. To evaluate these polynomials, we have solved block Toeplitz systems with the coefficients belonging to Q(i) and have found the ranks of matrices with entries from this field. All such operations can be performed exactly.

To obtain the left factorization essential polynomials, we can carry out similar construction with the sequence of left kernels of matrices T _k , - κ ≤ k ≤ κ, or can apply a conformance procedure (see [17], Def. 5.3). This procedure can be also fulfilled exactly . The conformance procedure that is used to construct the left factorization essential polynomials, is described in [17].

□

After finding the indices and factorization essential polynomials, we can exactly construct the Wiener–Hopf factorizations using the formulas (5) – (7).

3.    Pseudo-Code for an Exact Constructing the Right Factorization Essential Polynomials

The full variant of the pseudo-code for the algorithm of simultaneous construction of the left and right factorizations is given in [18]. However, if only one type of the factorization is needed (for instance, the left factorization), using the full algorithm leads to a significant increasing in execution time. For this reason, in this section we give the pseudo-code for construction of the left factorization only. For simplicity, here we restrict ourselves to the case when ^ _p <  ^ _p+1 .

Algorithm. Indices and right factorization essential polynomials of a sequence

Input . The sequence c — _K := {c _-K , ... ,c₀ ..., c _K } , C j E Q p^xp (i).

Output . The indices ^ ₁ ,...,^₂ p and the matrix of the right factorization essential polynomials, R ₁ := RR(zz) • • • R p (z)),

• 1.    find the distinct indices v₁,...,v_s, their multiplicities k₁,...,k_s, form the indices ^₁,..., ц_2р, and the number t such that ц_р = v_t

• 2.    find the polynomials R₁(z),..., R_K1(z) forming a basis of the space N_V1+1= ker T_V1+1, define the matrix R₁:= zzR((z) • • • zR_K1 (z))

• 3.    form the matrix a_n := ^r_r{z^-KR₁(z)} • • • a_R{z^-KR_K1(z)})

• 4.    for j = 2,... ,t do

• 5.    find a basis Q1(z),... ,Qnj (z) of the space N_Vj+1 = ker T_Vj+1

• 6.    for k = 1,..., nj do

• 7.      form the matrix a₂:= (a11 a_R{z^-KQ_k(z)})

• 8.     if rank a₂= rank a11 + 1 then

• 9.         au. := a2

• 10.    Ri := (Ri zQk(z))

• 11.    end if

• 12.    end do

• 13.    end do

• 14.    if rank a11 = p then

• 15.    print “The factorization essential polynomials were not constructed. The factorization process is interrupted”

• 16.    stop

• 17.    end if

19. return ц1,...,ц_2р, R1 (z) 4.    Numerical Example

z ^µ1   ...     0

^. .         ^. .              ^. .                  .

0   ...  z ^µ2p

Now by formulas (5), (6) we can construct the left factorization of a matrix polynomial.

Based on the proposed algorithm, a procedure ExactFEP was developed, which is the main part of the ExactMPF package in Maple. The package is designed for the exact solution of the factorization problem for matrix polynomials. To access ExactMPF use the commands

• >    read("ExactMPF.txt");

• >    with(ExactMPF);

• >    with(LinearAlgebra);

To obtain the factorizations of a(z) we run the module SolverExactMPF with the argument a ( z ):

• >    lplus, dl, lminus, rminus, dr, rplus := SolverExactMPF(a):

The module SolverExactMPF returns the factors lplus, dl, lminus of the left factorization and the factors rminus, dr, rplus of the right factorization.

Let us give an example of using this package.

Example 1. Consider

36z 2 + 17z — 14 z ⁴ — z ² + 3z — 1 a(z) : =         0           z ² + 13z +15

z + 10 z ²

)

0                  0

The module SolverExactMPF gives in this case the following expression for the factors of a(z):

• >    lplus; dl; lminus;

  
  

  z + 10 z ²

  
  

  z ⁴ — z ² + 3z —

  z ² + 13z + 15 0

  
  

  1 z ² 0 0

  , 0 10

  0 0 1

  
  

  ^{1 +} 36z

  
  

  _-

  
  
  
  
  
  

  7 18z ²

  
  
  
  

The executing time is 0,500 seconds when computations were performed on a home desktop computer HP with Intel(R) Core(TM)i3-415T CPU, 3.00 GHz, 4G RAM, operating system Windows 10.

Acknowledgments . V.M. Adukov and N.V. Adukova were supported by funding from RFBR grant no. 20-41-740024. G. Mishuris was supported by funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement EffectFact no. 101008140.

• G. Mishuris is thankful to the Royal Society for the Wolfson Research Merit Award and to the Welsh Government for the Future Generation Industrial Fellowship.