# On Solving System of Linear Differential-Algebraic Equations Using Reduction Algorithm

Authors: Srinivasarao Thota
Journal: International Journal of Mathematics and Mathematical Sciences
 (2020)

---

## Abstract

In this paper, we present a new reduction algorithm for solving system of linear differential-algebraic equations with power series coefficients. In the proposed algorithm, we transform the given system of differential-algebraic equations into another simple equivalent system using the elementary algebraic techniques. This algorithm would help to implement the manual calculations in commercial packages such as Mathematica, Maple, MATLAB, Singular, and Scilab. Maple implementation of the proposed algorithm is discussed, and sample computations are presented to illustrate the proposed algorithm.


---

## Body

## 1. Introduction

In many applications of science and engineering, for example, simulation of electric circuits [1–4], mechanical systems [5, 6], and chemical reactions subject to invariants [7–13], the systems of differential-algebraic equations (DAEs) arise naturally, and these systems of DAEs consist of algebraic equations and differential operations. Many engineers and scientists have studied the system of DAEs from a theoretical as well as a numerical point of view and created many new approaches to solve the system of linear differential-algebraic equations; see, for example, [14–25].In this paper, we are concerned with a linear system of differential-algebraic equations of the following form:(1)Ax∂ux+ℬxux=fx,where x is a complex variable, Ax and ℬx are m×n matrices with analytic functions entries, fx is an m-dimensional column matrix with analytic functions entries, ux is an n-dimensional unknown column matrix which is going to be determined, and ∂=d/dx is a differential operator. In this paper, we focus on creating a new reduction algorithm using elementary algebraic techniques as well as the implementation of the proposed algorithm in Maple. Using this algorithm, we can transform a given system of DAEs into another equivalent system, where we can solve the reduced system easily.The rest of the paper is organized as follows: in Section2, we present a new reduction algorithm to solve the given system of DAEs with certain examples to illustrate the proposed reduction algorithm and Section 3 discusses the Maple implementation of the proposed algorithm with sample computations.


## 2. A New Reduction Algorithm

LetK be a subfield of the field of complex numbers ℂ. Note that (ℚ⊆K⊆ℂ). We denote the ring of formal power series by Kx in the variable x and Kx denotes its quotient field, that is, Kx=Kxx−1. The ring of differential operators is denoted by Kx∂ with coefficients in Kx, that is, the set of finite sums ∑ai∂i with ai∈Kx is equipped with the addition and the multiplication defined by(2)∂i∂j=∂i+j,i,j∈ℕ,∂f=f∂+dfdx,where f∈Kx. Recall the system of DAEs in equation (1):(3)Ax∂ux+ℬxuz=fx,where Ax,ℬx∈Kxm×n, fx∈Kxm, and ux∈Kxn. The corresponding matrix differential operator of system (3) is L=A∂+ℬ∈Kx∂m×n. We recall the basic concepts of the matrix differential operators; see [23, 26–29] for further details.Definition 1.
A matrix differential operatorT∈Kx∂n×n is said to be unimodular matrix if there exists V∈Kx∂n×n such that VT=TV=In.Definition 2.
Two matrix differential operatorsL,L˜∈Kx∂m×n are said to be equivalent if there exist two unimodular matrices S∈Kx∂m×m and T∈Kx∂n×n such that L˜=SLT.
The following section presents a new reduction algorithm using elementary algebraic techniques.


### 2.1. Reduction Algorithm

The following lemma is one of the essential steps to create a new reduction algorithm. It shows that any matrix of formal power series centered at origin can be transformed into a block matrix.Lemma 1 (see [23, 26, 29]).
LetA∈Kxn×n. Then there exist two unimodular matrices S,T∈Kxn×n such that(4)SA=A110,AT=A110,SAT=A11000,whereAij denote the i-th row and j-th column block of matrix A, A11∈Kxr×r is a block matrix, and r is the rank of matrix A.Suppose thatL∈Kx∂m×n is a given matrix differential operator. Using Lemma 1, we can construct two unimodular matrix differential operators S1 and T1, by finding the basis of left null space and right null space of the matrix differential operator L, such that(5)L1=S1LT1=A1∂+ℬ1,where A1=A11000 and ℬ1=ℬ11000. Now using Lemma 1 to matrix A1 of (5), we can get an unimodular matrix S2 such that(6)L2=S2L1=A2∂+ℬ2,where A2=S2A1=A^11000 and A^11 is invertible matrix, and ℬ2=S2ℬ1=ℬ^11000. Again, using Lemma 1 to matrix ℬ2 of (6), we can construct an unimodular matrix T2 such that(7)L3=L2T2=A3∂+ℬ3,where A3=A2T2=A^11000 and ℬ3=ℬ2T2=ℬ^11000. If we denote S=S2S1 and T=T1T2, then we have two unimodular matrix differential operators S∈GLmKx and T∈GLnKx such that the given system is in reduced form as follows:(8)L˜=A˜∂+ℬ˜,where L˜=SLT has the form L˜11000 and L˜11 is invertible matrix differential operator; A˜ has the form A˜11000; and ℬ˜ has the form ℬ˜11000. Note that the ranks of A˜11 and ℬ˜11 may not be the same, but the ranks of A˜11 and ℬ˜11 are the same only when the coefficient matrices A and ℬ have the same rank. We generalize the above reduction algorithm in the following theorem.Theorem 1.
LetL=AD+ℬ∈Kx∂m×n be a matrix differential operator. Then we can construct two unimodular matrices S∈GLmKx and T∈GLnKx such that the given system L=A∂+ℬ∈Kx∂m×n is in reduced form:(9)L˜=A˜∂+ℬ˜,where L˜=SLT has the form L˜11000 and L˜11 is invertible matrix differential operator; A˜ has the form A˜11000; and ℬ˜ has the form ℬ˜11000.Proof.
Using Lemma1 to L, we can construct two unimodular matrices S∈GLmKx and T∈GLnKx. If we use the substitution ux=Tvx in (3) and left multiplying the resultant equation with S, we have the following reduced form:(10)SLTvx=SAT∂vx+SℬTvx=Sfx,orL˜v=A˜∂v+ℬ˜v=f˜,where L˜=SLT∈Kx∂m×n, A˜=SAT∈Kxm×n, ℬ˜=SℬT∈Kxm×n, and f˜=Sf∈Kxm×1.
Indeed, (i) if rank (L) < rank (A) and rank (L) < rank (ℬ), then the reduced system of DAEs (10) has the following form:(11)L˜v=A˜∂v+ℬ˜v=f˜,where L˜=L¯11∂L¯120L¯21L¯220000,v=v1v2v3,f˜=f¯1f¯2f¯3. Hence, the system of DAEs in (3) is decomposed into two systems as follows:(12)L¯11∂v1+L¯12v2=f¯1,L¯21v1+L¯22v2=f¯2,with some necessary conditions on the right-hand side expressed by f¯3=0.(ii)
If rank (L) = rank (A) = rank (ℬ), then the reduced DAS (10) has the following form:(13)L˜v=A˜∂v+ℬ˜v=f˜,where L˜=L¯11000,v=v1v2,f˜=f¯1f¯2. Hence, the system of DAEs in (3) is decomposed as L¯11v1=f¯1 with some necessary conditions on the right-hand side expressed by f¯2=0.
In the following section, we present certain examples to illustrate the proposed method presented in Theorem1.


### 2.2. Examples

Example 1.
Consider a matrix differential operator of DAS.(14)L=A∂+ℬ=0x1−x11100101x1−x211x−11021−xx−10−1∂+1001000−100000000x+101000000=1x∂1−x∂∂+1∂∂0−1∂0∂x∂1−x∂2∂∂∂x−1∂+x+1∂12∂∂−x∂x−1∂0−∂.
Using Lemma1 to L, one can construct two unimodular matrices S1 and T1 (obtained using a basis of left null space and right null space of L) as follows:(15)left null space ofL=2∂∂−10∂+11−∂01,right null space ofL=x+1−x+1∂2x−1∂+11−3∂2x−1∂+10x+1∂−x−12x−1∂+11.
Now the unimodular matricesS1 and T1 are(16)S1=100000100000100000102∂∂−10∂+11−∂01,T1=1000x+1−x+1∂2x−1∂+101001−3∂2x−1∂+1001000001x+1∂−x−12x−1∂+100001.
Thus, multiplying operatorL on the left and right by S1 and T1 yields the operator(17)L1=S1LT1=1x∂1−x∂∂+1⋮0∂0−1∂⋮0∂x∂1−x∂2∂⋮0∂x−1∂+x+1∂1⋮0………………………………0000⋮0,where(18)A1=A11000=0x1−x1⋮01001⋮01x1−x2⋮01x−110⋮0……………0000⋮0,ℬ1=ℬ11000=1001⋮000−10⋮00000⋮00x+101⋮0……………0000⋮0.
Now, using Lemma1 to matrix A1, we can construct an unimodular matrix S2, using a basis of left null space of A1, as follows:(19)S2=100000100000010−1−110000001.
We have(20)L3=S2L2=A2∂+ℬ2=1x∂1−x∂∂+1⋮0∂0−1∂⋮0∂x−1∂+x+1∂1⋮0−101−1⋮0……………………………0000⋮0,where(21)A2=A^11000=0x1−x1⋮01001⋮01x−110⋮0……………0000⋮00000⋮0,ℬ2=ℬ^11000=1001⋮000−10⋮00x+101⋮0−101−1⋮0……………0000⋮0.
Again, using Lemma1 to matrix ℬ2 of the matrix differential operator L2, we can construct an unimodular matrix T2 using a basis of right null space of ℬ2 as follows:(22)T2=100−10010−1x+10001000001000001.
We have(23)L3=L2T2=A3∂+ℬ3=1x∂1−x∂∂1+x⋮0∂0−10⋮0∂x−1∂+x+1∂−2x∂1+x⋮0−1010⋮0……………………………0000⋮0,where(24)A3=A2T2=A˜11000=0x1−x1x+1⋮01000⋮01x−11−2x1+x⋮0……………0000⋮00000⋮0,ℬ3=ℬ2T2=ℬ˜11000=1000⋮000−10⋮00x+101⋮0−1010⋮0……………0000⋮0.
If we denoteS=S2S1 and T=T1T2, then we have two unimodular matrix differential operators:(25)S=S2S1=100000100000010−1−11002∂∂−10∂+11−∂01,T=T1T2=100−1x+1−x+1∂2x−1∂+1010−1x+11−3∂2x−1∂+1001000001x+1∂−x−12x−1∂+100001.
We have that the given system is in reduced form as follows:(26)L˜=A˜∂+ℬ˜=1x∂1−x∂∂1+x⋮0∂0−10⋮0∂x−1∂+x+1∂−2x∂1+x⋮0−1010⋮0……………………………0000⋮0,where L˜=SLT has the form L˜11000 and(27)L˜11=1x∂1−x∂∂1+x∂0−10∂x−1∂+x+1∂−2x∂1+x−1010,is invertible matrix differential operator.Example 2.
Consider a matrix differential operator as given below:(28)L=∂1+x∂1∂+11−x∂∂−1∂+1∂∂+1∂−∂+1−2x−1∂∂−1−∂−2=1x+1010−x101111−1−2x−11−1∂+0011100−1101010−1−2.
Applying the proposed algorithm in Theorem1 to matrix differential operator (28) similar to example 1, one can construct two unimodular matrix differential operators S and T as(29)S=10000100−1−1101−101,T=101−x∂x∂+x+1−x∂+∂+xx∂+x+101−2∂+1∂x∂+x+1∂2−1∂x∂+x+100010010,such that the given system is in reduced form as follows:(30)L˜=∂1+x∂001−x∂0000000000,where L˜=SLT has the form L˜11000 and(31)L˜11=∂1+x∂1−x∂,is invertible matrix differential operator.Example 3.
Consider the following system of differential-algebraic equations to verify that the reduced system and the given system of DAEs have the same solution:(32)u1′+u3′+u3=0,u1+2u2′+u3=x,2u1′+2u3′+2u3=0.
The solution of the given system (32) is u1=c1,u2=x2/2+c2e−x/2−c1x/2+c3 and u3=c2e−x. In particular, if we take u1=0, then the solution becomes u1=0,u2=x2/2+c2e−x/2+c2 and u3=c1e−x.
The operator notationLu=f of the given system (14) is given by(33)∂0∂+112∂12∂02∂+1u1u2u3=0x0,where L=∂0∂+112∂12∂02∂+1. Now, apply the proposed algorithm to matrix differential operator L to get a reduced operator L˜ with two unimodular matrix differential operators S and T. We get(34)L˜=∂0012∂0000,S=100010−201,T=10−1+∂∂0112∂2001.
System (33) is reduced to L˜zx=f˜x, where ux=Tzx and f˜x=Sfx. On simplification, we can get zx=T−1ux and f˜x=Sfx.
Therefore, we have(35)L˜zx=f˜x⟹L˜T−1ux=Sfx⇒∂0012∂000010−1+∂∂0112∂2001−1u1u2u3=100010−2010x0⇒∂0012∂0000u1+∂+1∂u3u2−u32∂2u3=0x0⇒∂u1+∂+1u3u1+2∂u2+u30=0x0.
Now the reduced system of DAEs is(36)u1′+u3′+u3=0,u1+2u2′+u3=x.
Solution of the reduced system (36), for u1=0, is u1=0,u2=x2/2+c1e−x/2+c2 and u3=c1e−x.
One can observe that the solution of the given system of DAEs (32) and the reduced system of DAEs (36) have the same solution. We can also observe that solving the reduced system (36) (contains two equations only) is simple compared to solving the given system (32) (contains three equations).


## 2.1. Reduction Algorithm

The following lemma is one of the essential steps to create a new reduction algorithm. It shows that any matrix of formal power series centered at origin can be transformed into a block matrix.Lemma 1 (see [23, 26, 29]).
LetA∈Kxn×n. Then there exist two unimodular matrices S,T∈Kxn×n such that(4)SA=A110,AT=A110,SAT=A11000,whereAij denote the i-th row and j-th column block of matrix A, A11∈Kxr×r is a block matrix, and r is the rank of matrix A.Suppose thatL∈Kx∂m×n is a given matrix differential operator. Using Lemma 1, we can construct two unimodular matrix differential operators S1 and T1, by finding the basis of left null space and right null space of the matrix differential operator L, such that(5)L1=S1LT1=A1∂+ℬ1,where A1=A11000 and ℬ1=ℬ11000. Now using Lemma 1 to matrix A1 of (5), we can get an unimodular matrix S2 such that(6)L2=S2L1=A2∂+ℬ2,where A2=S2A1=A^11000 and A^11 is invertible matrix, and ℬ2=S2ℬ1=ℬ^11000. Again, using Lemma 1 to matrix ℬ2 of (6), we can construct an unimodular matrix T2 such that(7)L3=L2T2=A3∂+ℬ3,where A3=A2T2=A^11000 and ℬ3=ℬ2T2=ℬ^11000. If we denote S=S2S1 and T=T1T2, then we have two unimodular matrix differential operators S∈GLmKx and T∈GLnKx such that the given system is in reduced form as follows:(8)L˜=A˜∂+ℬ˜,where L˜=SLT has the form L˜11000 and L˜11 is invertible matrix differential operator; A˜ has the form A˜11000; and ℬ˜ has the form ℬ˜11000. Note that the ranks of A˜11 and ℬ˜11 may not be the same, but the ranks of A˜11 and ℬ˜11 are the same only when the coefficient matrices A and ℬ have the same rank. We generalize the above reduction algorithm in the following theorem.Theorem 1.
LetL=AD+ℬ∈Kx∂m×n be a matrix differential operator. Then we can construct two unimodular matrices S∈GLmKx and T∈GLnKx such that the given system L=A∂+ℬ∈Kx∂m×n is in reduced form:(9)L˜=A˜∂+ℬ˜,where L˜=SLT has the form L˜11000 and L˜11 is invertible matrix differential operator; A˜ has the form A˜11000; and ℬ˜ has the form ℬ˜11000.Proof.
Using Lemma1 to L, we can construct two unimodular matrices S∈GLmKx and T∈GLnKx. If we use the substitution ux=Tvx in (3) and left multiplying the resultant equation with S, we have the following reduced form:(10)SLTvx=SAT∂vx+SℬTvx=Sfx,orL˜v=A˜∂v+ℬ˜v=f˜,where L˜=SLT∈Kx∂m×n, A˜=SAT∈Kxm×n, ℬ˜=SℬT∈Kxm×n, and f˜=Sf∈Kxm×1.
Indeed, (i) if rank (L) < rank (A) and rank (L) < rank (ℬ), then the reduced system of DAEs (10) has the following form:(11)L˜v=A˜∂v+ℬ˜v=f˜,where L˜=L¯11∂L¯120L¯21L¯220000,v=v1v2v3,f˜=f¯1f¯2f¯3. Hence, the system of DAEs in (3) is decomposed into two systems as follows:(12)L¯11∂v1+L¯12v2=f¯1,L¯21v1+L¯22v2=f¯2,with some necessary conditions on the right-hand side expressed by f¯3=0.(ii)
If rank (L) = rank (A) = rank (ℬ), then the reduced DAS (10) has the following form:(13)L˜v=A˜∂v+ℬ˜v=f˜,where L˜=L¯11000,v=v1v2,f˜=f¯1f¯2. Hence, the system of DAEs in (3) is decomposed as L¯11v1=f¯1 with some necessary conditions on the right-hand side expressed by f¯2=0.
In the following section, we present certain examples to illustrate the proposed method presented in Theorem1.


## 2.2. Examples

Example 1.
Consider a matrix differential operator of DAS.(14)L=A∂+ℬ=0x1−x11100101x1−x211x−11021−xx−10−1∂+1001000−100000000x+101000000=1x∂1−x∂∂+1∂∂0−1∂0∂x∂1−x∂2∂∂∂x−1∂+x+1∂12∂∂−x∂x−1∂0−∂.
Using Lemma1 to L, one can construct two unimodular matrices S1 and T1 (obtained using a basis of left null space and right null space of L) as follows:(15)left null space ofL=2∂∂−10∂+11−∂01,right null space ofL=x+1−x+1∂2x−1∂+11−3∂2x−1∂+10x+1∂−x−12x−1∂+11.
Now the unimodular matricesS1 and T1 are(16)S1=100000100000100000102∂∂−10∂+11−∂01,T1=1000x+1−x+1∂2x−1∂+101001−3∂2x−1∂+1001000001x+1∂−x−12x−1∂+100001.
Thus, multiplying operatorL on the left and right by S1 and T1 yields the operator(17)L1=S1LT1=1x∂1−x∂∂+1⋮0∂0−1∂⋮0∂x∂1−x∂2∂⋮0∂x−1∂+x+1∂1⋮0………………………………0000⋮0,where(18)A1=A11000=0x1−x1⋮01001⋮01x1−x2⋮01x−110⋮0……………0000⋮0,ℬ1=ℬ11000=1001⋮000−10⋮00000⋮00x+101⋮0……………0000⋮0.
Now, using Lemma1 to matrix A1, we can construct an unimodular matrix S2, using a basis of left null space of A1, as follows:(19)S2=100000100000010−1−110000001.
We have(20)L3=S2L2=A2∂+ℬ2=1x∂1−x∂∂+1⋮0∂0−1∂⋮0∂x−1∂+x+1∂1⋮0−101−1⋮0……………………………0000⋮0,where(21)A2=A^11000=0x1−x1⋮01001⋮01x−110⋮0……………0000⋮00000⋮0,ℬ2=ℬ^11000=1001⋮000−10⋮00x+101⋮0−101−1⋮0……………0000⋮0.
Again, using Lemma1 to matrix ℬ2 of the matrix differential operator L2, we can construct an unimodular matrix T2 using a basis of right null space of ℬ2 as follows:(22)T2=100−10010−1x+10001000001000001.
We have(23)L3=L2T2=A3∂+ℬ3=1x∂1−x∂∂1+x⋮0∂0−10⋮0∂x−1∂+x+1∂−2x∂1+x⋮0−1010⋮0……………………………0000⋮0,where(24)A3=A2T2=A˜11000=0x1−x1x+1⋮01000⋮01x−11−2x1+x⋮0……………0000⋮00000⋮0,ℬ3=ℬ2T2=ℬ˜11000=1000⋮000−10⋮00x+101⋮0−1010⋮0……………0000⋮0.
If we denoteS=S2S1 and T=T1T2, then we have two unimodular matrix differential operators:(25)S=S2S1=100000100000010−1−11002∂∂−10∂+11−∂01,T=T1T2=100−1x+1−x+1∂2x−1∂+1010−1x+11−3∂2x−1∂+1001000001x+1∂−x−12x−1∂+100001.
We have that the given system is in reduced form as follows:(26)L˜=A˜∂+ℬ˜=1x∂1−x∂∂1+x⋮0∂0−10⋮0∂x−1∂+x+1∂−2x∂1+x⋮0−1010⋮0……………………………0000⋮0,where L˜=SLT has the form L˜11000 and(27)L˜11=1x∂1−x∂∂1+x∂0−10∂x−1∂+x+1∂−2x∂1+x−1010,is invertible matrix differential operator.Example 2.
Consider a matrix differential operator as given below:(28)L=∂1+x∂1∂+11−x∂∂−1∂+1∂∂+1∂−∂+1−2x−1∂∂−1−∂−2=1x+1010−x101111−1−2x−11−1∂+0011100−1101010−1−2.
Applying the proposed algorithm in Theorem1 to matrix differential operator (28) similar to example 1, one can construct two unimodular matrix differential operators S and T as(29)S=10000100−1−1101−101,T=101−x∂x∂+x+1−x∂+∂+xx∂+x+101−2∂+1∂x∂+x+1∂2−1∂x∂+x+100010010,such that the given system is in reduced form as follows:(30)L˜=∂1+x∂001−x∂0000000000,where L˜=SLT has the form L˜11000 and(31)L˜11=∂1+x∂1−x∂,is invertible matrix differential operator.Example 3.
Consider the following system of differential-algebraic equations to verify that the reduced system and the given system of DAEs have the same solution:(32)u1′+u3′+u3=0,u1+2u2′+u3=x,2u1′+2u3′+2u3=0.
The solution of the given system (32) is u1=c1,u2=x2/2+c2e−x/2−c1x/2+c3 and u3=c2e−x. In particular, if we take u1=0, then the solution becomes u1=0,u2=x2/2+c2e−x/2+c2 and u3=c1e−x.
The operator notationLu=f of the given system (14) is given by(33)∂0∂+112∂12∂02∂+1u1u2u3=0x0,where L=∂0∂+112∂12∂02∂+1. Now, apply the proposed algorithm to matrix differential operator L to get a reduced operator L˜ with two unimodular matrix differential operators S and T. We get(34)L˜=∂0012∂0000,S=100010−201,T=10−1+∂∂0112∂2001.
System (33) is reduced to L˜zx=f˜x, where ux=Tzx and f˜x=Sfx. On simplification, we can get zx=T−1ux and f˜x=Sfx.
Therefore, we have(35)L˜zx=f˜x⟹L˜T−1ux=Sfx⇒∂0012∂000010−1+∂∂0112∂2001−1u1u2u3=100010−2010x0⇒∂0012∂0000u1+∂+1∂u3u2−u32∂2u3=0x0⇒∂u1+∂+1u3u1+2∂u2+u30=0x0.
Now the reduced system of DAEs is(36)u1′+u3′+u3=0,u1+2u2′+u3=x.
Solution of the reduced system (36), for u1=0, is u1=0,u2=x2/2+c1e−x/2+c2 and u3=c1e−x.
One can observe that the solution of the given system of DAEs (32) and the reduced system of DAEs (36) have the same solution. We can also observe that solving the reduced system (36) (contains two equations only) is simple compared to solving the given system (32) (contains three equations).


## 3. Maple Implementation

In this section, we discuss the Maple implementation of the algorithm by creating different data types. Using the Maple package, one can obtain the two unimodular matrix differential operatorsS,T and the reduced matrix differential operator of the given system. In Maple implementation, x is complex variable and δ=d/dx is the differential operator.Input: A and ℬ, the coefficient matrices of a given matrix differential operator L=A∂+ℬ.Output:L˜, S, and T, the reduced matrix differential operator L˜ of a given matrix differential operator L and two unimodular matrix differential operators S and T.


### 3.1. Pseudocode

(1)
A,B←coefficientmatricies(2)
m←rowdimensionofL(3)
n←columndimensionofL(4)
r←rankofA(5)
k←rankofB(6)
NS_Ls←leftnullspaceofL(7)
S1←identitymatrixwithLNS_Lsasbottomblockmatrix(8)
NS_Lt←rightnullspaceofL(9)
T1←identitymatrixwithNS_Ltasrightblockmatrix(10)
L1←S1.L.T1(11)
A1,B1←coefficientmatriciesofL1(12)
NS_A1←leftnullspaceofA1(13)
S2←identitymatrixwithLNS_A1asbottomblockmatrix(14)
L2←S2.L1(15)
A2,B2←coefficientmatriciesofL2(16)
NS_B2←rightnullspaceofB2(17)
T2←identitymatrixwithNS_B2asrightblockmatrix(18)
L3←L3.T2(19)
S←S2.S1(20)
T←T1.T2


### 3.2. Maple Code

Using the procedure ArrangeZeroLines, we can replace the zero rows at the bottom and zero columns at right side of the matrix.ArrangeZeroLines := proc (M::Matrix)local m,n,ZR,L1,DR,m1,ZM,A1,ZC,L2,DC,m2,ZM2,A2;uses LinearAlgebra;m,n := op(1, M);ZR := LinearAlgebra: −ZeroVector[row](n);L1 := [seq(‘if‘(LinearAlgebra: −Equal(M[i],ZR),i,NULL), i = 1..m)];DR := LinearAlgebra: −DeleteRow(M, L1);m1 := op([1, 1], DR);ZM := LinearAlgebra: −ZeroMatrix(m − m1, n);A1:= convert(linalg[blockmatrix](2, 1, [DR, ZM]), Matrix);ZC:= LinearAlgebra: −ZeroVector(m);L2 := [seq(‘if‘(LinearAlgebra: −Equal(A1[()..(),i],ZC),i,NULL), i = 1..n)];DC := LinearAlgebra: −DeleteColumn(A1, L2);m2 := op([1, 2], DC);ZM2 := LinearAlgebra: −ZeroMatrix(m, n − m2);A2 := linalg[blockmatrix](1, 2, [DC, ZM2]);return convert(A2, Matrix);end proc:The following procedure is DAEs_Reduction: the reduced matrix differential operator of a given matrix differential operator with two unimodular matrix differential operators. In this procedure,δ=d/dz is differential operator and x is complex variable.DAEs_Reduction := proc (A::Matrix, B::Matrix)local L,m,n,r,k,Id,transpose_L,NS_Ls,Id_partA, NS_partA,S1,NS_Lt,Id_part3,NS_part3,T1,L1na,L1,A1,B1,transpose_A1,NS_A1,Id_partA1,NS_partA1,S2,L2na,L2,A2,B2,NS_B2,Id_partB2,NS_partB2,T2,L3na,L3,A3,B3,S,T;uses MatrixPolynomialAlgebra;L := A∗delta + B;m := LinearAlgebra: −RowDimension(L);n := LinearAlgebra: −ColumnDimension(L);r := MTM: −rank(A);k := MTM: −rank(B);Id := LinearAlgebra: −IdentityMatrix(n);transpose_L := LinearAlgebra: −Transpose(L);NS_Ls := LinearAlgebra: −NullSpace(transpose_L);Id_partA := LinearAlgebra: −Transpose(Matrix(‘∼‘[LinearAlgebra: −Transpose](‘∼‘[convert]([seq(Id[i],i = 1..m-nops(NS_Ls))],Matrix))));NS_partA := LinearAlgebra: −Transpose (Matrix(‘∼‘[convert]([seq(NS_Ls[i], i = 1 .. Nops(NS_Ls))], Matrix)));S1 := convert(linalg: −blockmatrix(2,1,[Id_partA, NS_partA]),Matrix);NS_Lt := ‘∼‘[simplify](LinearAlgebra: −NullSpace(L));Id_part3 := Matrix(‘∼‘[LinearAlgebra: −Transpose](‘∼‘[convert]([seq(Id[i], i = 1 .. n-nops(NS_Lt))], Matrix)));NS_part3 := Matrix(‘∼‘[convert]([seq(NS_Lt[i], i = 1 .. Nops(NS_Lt))], Matrix));T1 := convert(linalg: −blockmatrix(1,2,[Id_part3,NS_part3]),Matrix);L1na := ‘∼‘[simplify](S1.L.T1);L1 := ArrangeZeroLines(L1na);A1:= seq(MatrixPolynomialAlgebra: −Coeff(L1,delta,i),i = 0..1)[2];B1 := seq(MatrixPolynomialAlgebra: −Coeff(L1,delta,i),i = 0..1)[1];transpose_A1:= LinearAlgebra: −Transpose(A1);NS_A1:= LinearAlgebra: −NullSpace(transpose_A1);Id_partA1:=LinearAlgebra: −Transpose(Matrix(‘∼‘[LinearAlgebra: −Transpose](‘∼‘[convert]([seq(Id[i], i = 1 .. m-nops(NS_A1))], Matrix))));NS_partA1 := LinearAlgebra: −Transpose(Matrix(‘∼‘[convert]([seq(NS_A1[i], i = 1 .. Nops(NS_A1))], Matrix)));S2 := convert(linalg: −blockmatrix(2,1,[Id_partA1,NS_partA1]),Matrix);L2na := ‘∼‘[simplify](S2.L1);L2 := ArrangeZeroLines(L2na);A2 := seq(MatrixPolynomialAlgebra: −Coeff(L2,delta,i),i = 0..1)[2];B2 := seq(MatrixPolynomialAlgebra: −Coeff(L2,delta,i),i = 0..1)[1];NS_B2 := ‘∼‘[simplify](LinearAlgebra: −NullSpace(B2));Id_partB2 := Matrix(‘∼‘[LinearAlgebra: −Transpose](‘∼‘[convert]([seq(Id[i],i = nops(NS_B2)+1..n)],Matrix)));NS_partB2 := Matrix(‘∼‘[convert]([seq(NS_B2[i], i = 1..nops(NS_B2))], Matrix));T2 := LinearAlgebra: −IdentityMatrix(n);L3na := ‘∼‘[simplify](L2.T2);L3 := ArrangeZeroLines(L3na);A3 := seq(MatrixPolynomialAlgebra: −Coeff(L3,delta,i),i = 0..1)[2];B3 := seq(MatrixPolynomialAlgebra: −Coeff(L3,delta,i),i = 0..1)[1];S := ‘∼‘[simplify](S2.S1);T := ‘∼‘[simplify](T1.T2);return S, T, L3end proc:


### 3.3. Sample Computations

Example 4.
Consider the following matrix differential operator as given in Example2:(37)L=∂1+x∂1∂+11−x∂∂−1∂+1∂∂+1∂−∂+1−2x−1∂∂−1−∂−2,where A=1x+1010−x101111−1−2x−11−1 and ℬ=0011100−1101010−1−2. Using Maple implementation, with DAEs_Reduction in the proposed algorithm, we have the following:>A := Matrix([[1, x + 1, 0, 1], [0, −x, 1, 0], [1, 1, 1, 1], [−1, −2x − 1, 1, −1]]).>B := Matrix([[0, 0, 1, 1], [1, 0, 0, −1], [1, 0, 1, 0], [1, 0,−1, −2]]).(38)A≔1x+1010−x101111−1−2x−11−1,B≔0011100−1101010−1−2.>S, T, Lred := DAEs_Reduction(A, B).(39)10000100−1−1101−101,10−δx−1δx+x+1−δx+δ+xδx+x+101−2δ+1δδx+x+1δ2−1δδx+x+100010010,δδ1+x001−δx0000000000.
From Maple implementation, we have that two unimodular matrices differential operatorsS,T are(40)S=10000100−1−1101−101,T=101−x∂x∂+x+1−x∂+∂+xx∂+x+101−2∂+1∂x∂+x+1∂2−1∂x∂+x+100010010,and the reduced matrix differential operator of the given DAS (37) is(41)L˜=∂1+x∂001−x∂0000000000.
One can also verify in Maple thatL˜=SLT as follows:>L := Matrix([[delta, delta∗(1 + x), 1, delta + 1],[1, −delta∗x, delta, −1], [delta + 1, delta, delta + 1, delta],[−delta + 1, −delta∗(1 + 2∗x), delta − 1, −delta − 2]]);(42)δδ1+x1δ+11−δxδ−1δ+1δδ+1δ−δ+1−δ1+2xδ−1−δ−2.>simplify∼(S.L.T);(43)δδ1+x001−δx0000000000.


## 3.1. Pseudocode

(1)
A,B←coefficientmatricies(2)
m←rowdimensionofL(3)
n←columndimensionofL(4)
r←rankofA(5)
k←rankofB(6)
NS_Ls←leftnullspaceofL(7)
S1←identitymatrixwithLNS_Lsasbottomblockmatrix(8)
NS_Lt←rightnullspaceofL(9)
T1←identitymatrixwithNS_Ltasrightblockmatrix(10)
L1←S1.L.T1(11)
A1,B1←coefficientmatriciesofL1(12)
NS_A1←leftnullspaceofA1(13)
S2←identitymatrixwithLNS_A1asbottomblockmatrix(14)
L2←S2.L1(15)
A2,B2←coefficientmatriciesofL2(16)
NS_B2←rightnullspaceofB2(17)
T2←identitymatrixwithNS_B2asrightblockmatrix(18)
L3←L3.T2(19)
S←S2.S1(20)
T←T1.T2


## 3.2. Maple Code

Using the procedure ArrangeZeroLines, we can replace the zero rows at the bottom and zero columns at right side of the matrix.ArrangeZeroLines := proc (M::Matrix)local m,n,ZR,L1,DR,m1,ZM,A1,ZC,L2,DC,m2,ZM2,A2;uses LinearAlgebra;m,n := op(1, M);ZR := LinearAlgebra: −ZeroVector[row](n);L1 := [seq(‘if‘(LinearAlgebra: −Equal(M[i],ZR),i,NULL), i = 1..m)];DR := LinearAlgebra: −DeleteRow(M, L1);m1 := op([1, 1], DR);ZM := LinearAlgebra: −ZeroMatrix(m − m1, n);A1:= convert(linalg[blockmatrix](2, 1, [DR, ZM]), Matrix);ZC:= LinearAlgebra: −ZeroVector(m);L2 := [seq(‘if‘(LinearAlgebra: −Equal(A1[()..(),i],ZC),i,NULL), i = 1..n)];DC := LinearAlgebra: −DeleteColumn(A1, L2);m2 := op([1, 2], DC);ZM2 := LinearAlgebra: −ZeroMatrix(m, n − m2);A2 := linalg[blockmatrix](1, 2, [DC, ZM2]);return convert(A2, Matrix);end proc:The following procedure is DAEs_Reduction: the reduced matrix differential operator of a given matrix differential operator with two unimodular matrix differential operators. In this procedure,δ=d/dz is differential operator and x is complex variable.DAEs_Reduction := proc (A::Matrix, B::Matrix)local L,m,n,r,k,Id,transpose_L,NS_Ls,Id_partA, NS_partA,S1,NS_Lt,Id_part3,NS_part3,T1,L1na,L1,A1,B1,transpose_A1,NS_A1,Id_partA1,NS_partA1,S2,L2na,L2,A2,B2,NS_B2,Id_partB2,NS_partB2,T2,L3na,L3,A3,B3,S,T;uses MatrixPolynomialAlgebra;L := A∗delta + B;m := LinearAlgebra: −RowDimension(L);n := LinearAlgebra: −ColumnDimension(L);r := MTM: −rank(A);k := MTM: −rank(B);Id := LinearAlgebra: −IdentityMatrix(n);transpose_L := LinearAlgebra: −Transpose(L);NS_Ls := LinearAlgebra: −NullSpace(transpose_L);Id_partA := LinearAlgebra: −Transpose(Matrix(‘∼‘[LinearAlgebra: −Transpose](‘∼‘[convert]([seq(Id[i],i = 1..m-nops(NS_Ls))],Matrix))));NS_partA := LinearAlgebra: −Transpose (Matrix(‘∼‘[convert]([seq(NS_Ls[i], i = 1 .. Nops(NS_Ls))], Matrix)));S1 := convert(linalg: −blockmatrix(2,1,[Id_partA, NS_partA]),Matrix);NS_Lt := ‘∼‘[simplify](LinearAlgebra: −NullSpace(L));Id_part3 := Matrix(‘∼‘[LinearAlgebra: −Transpose](‘∼‘[convert]([seq(Id[i], i = 1 .. n-nops(NS_Lt))], Matrix)));NS_part3 := Matrix(‘∼‘[convert]([seq(NS_Lt[i], i = 1 .. Nops(NS_Lt))], Matrix));T1 := convert(linalg: −blockmatrix(1,2,[Id_part3,NS_part3]),Matrix);L1na := ‘∼‘[simplify](S1.L.T1);L1 := ArrangeZeroLines(L1na);A1:= seq(MatrixPolynomialAlgebra: −Coeff(L1,delta,i),i = 0..1)[2];B1 := seq(MatrixPolynomialAlgebra: −Coeff(L1,delta,i),i = 0..1)[1];transpose_A1:= LinearAlgebra: −Transpose(A1);NS_A1:= LinearAlgebra: −NullSpace(transpose_A1);Id_partA1:=LinearAlgebra: −Transpose(Matrix(‘∼‘[LinearAlgebra: −Transpose](‘∼‘[convert]([seq(Id[i], i = 1 .. m-nops(NS_A1))], Matrix))));NS_partA1 := LinearAlgebra: −Transpose(Matrix(‘∼‘[convert]([seq(NS_A1[i], i = 1 .. Nops(NS_A1))], Matrix)));S2 := convert(linalg: −blockmatrix(2,1,[Id_partA1,NS_partA1]),Matrix);L2na := ‘∼‘[simplify](S2.L1);L2 := ArrangeZeroLines(L2na);A2 := seq(MatrixPolynomialAlgebra: −Coeff(L2,delta,i),i = 0..1)[2];B2 := seq(MatrixPolynomialAlgebra: −Coeff(L2,delta,i),i = 0..1)[1];NS_B2 := ‘∼‘[simplify](LinearAlgebra: −NullSpace(B2));Id_partB2 := Matrix(‘∼‘[LinearAlgebra: −Transpose](‘∼‘[convert]([seq(Id[i],i = nops(NS_B2)+1..n)],Matrix)));NS_partB2 := Matrix(‘∼‘[convert]([seq(NS_B2[i], i = 1..nops(NS_B2))], Matrix));T2 := LinearAlgebra: −IdentityMatrix(n);L3na := ‘∼‘[simplify](L2.T2);L3 := ArrangeZeroLines(L3na);A3 := seq(MatrixPolynomialAlgebra: −Coeff(L3,delta,i),i = 0..1)[2];B3 := seq(MatrixPolynomialAlgebra: −Coeff(L3,delta,i),i = 0..1)[1];S := ‘∼‘[simplify](S2.S1);T := ‘∼‘[simplify](T1.T2);return S, T, L3end proc:


## 3.3. Sample Computations

Example 4.
Consider the following matrix differential operator as given in Example2:(37)L=∂1+x∂1∂+11−x∂∂−1∂+1∂∂+1∂−∂+1−2x−1∂∂−1−∂−2,where A=1x+1010−x101111−1−2x−11−1 and ℬ=0011100−1101010−1−2. Using Maple implementation, with DAEs_Reduction in the proposed algorithm, we have the following:>A := Matrix([[1, x + 1, 0, 1], [0, −x, 1, 0], [1, 1, 1, 1], [−1, −2x − 1, 1, −1]]).>B := Matrix([[0, 0, 1, 1], [1, 0, 0, −1], [1, 0, 1, 0], [1, 0,−1, −2]]).(38)A≔1x+1010−x101111−1−2x−11−1,B≔0011100−1101010−1−2.>S, T, Lred := DAEs_Reduction(A, B).(39)10000100−1−1101−101,10−δx−1δx+x+1−δx+δ+xδx+x+101−2δ+1δδx+x+1δ2−1δδx+x+100010010,δδ1+x001−δx0000000000.
From Maple implementation, we have that two unimodular matrices differential operatorsS,T are(40)S=10000100−1−1101−101,T=101−x∂x∂+x+1−x∂+∂+xx∂+x+101−2∂+1∂x∂+x+1∂2−1∂x∂+x+100010010,and the reduced matrix differential operator of the given DAS (37) is(41)L˜=∂1+x∂001−x∂0000000000.
One can also verify in Maple thatL˜=SLT as follows:>L := Matrix([[delta, delta∗(1 + x), 1, delta + 1],[1, −delta∗x, delta, −1], [delta + 1, delta, delta + 1, delta],[−delta + 1, −delta∗(1 + 2∗x), delta − 1, −delta − 2]]);(42)δδ1+x1δ+11−δxδ−1δ+1δδ+1δ−δ+1−δ1+2xδ−1−δ−2.>simplify∼(S.L.T);(43)δδ1+x001−δx0000000000.


## 4. Conclusion

In this paper, we discussed a new reduction algorithm to solve a given system of linear differential-algebraic equations with power series coefficients. With the proposed algorithm, one can transform the given matrix differential operator into another simple equivalent matrix differential operator using the elementary algebraic techniques. Certain examples are presented to illustrate the algorithm. The implemented Maple package is discussed and sample computations are presented.



---

*Source: 6671926-2020-12-16.xml*