Analytic Methods for Differential Algebraic Equations
aa r X i v : . [ m a t h . G M ] J a n Analytic Methods for Differential Algebraic Equations
Samiya A Alkhairy Earth, Atmospheric, and Planetary Sciences DepartmentMassachusetts Institute of Technology77 Mass Ave Office 54-317, Cambridge, MA 02139
Abstract
We introduce methods for deriving analytic solutions from differential-algebraic systems of equations(DAEs), as well as methods for deriving governing equations for analytic characterization which is currentlylimited to very small systems as it is carried out by hand. Analytic solutions to the system and analyticcharacterization through governing equations provide insights into the behaviors of DAEs as well as theparametric regions of operation for each potential behavior. DAEs are a mixture of algebraic, ordinarydifferential, and/or partial differential equations which are multivariate constitutive equations that arisefrom fundamental laws of a discipline which relate multiple dependent variables. For each system (DAEs),and choice of dependent variable, there is a corresponding governing equation which is univariate ODEor PDE that is typically higher order than the constitutive equations of the system. We first introducea direct formulation for representing systems of linear differential-algebraic scalar constitutive equationsthat are homogeneous or nonhomogeneous with constant or nonconstant coefficients. Unlike state spaceformulations, our formulation follows very directly from the system of constitutive equations without theneed for introducing state variables or singular matrices. Using this formulation for the system of constitutiveequations (DAEs), we develop methods for deriving analytic expressions for the full solution (complementaryand particular) for all dependent variables of systems that consist of constant coefficient ordinary-DAEs andspecial cases of partial-DAEs. We also develop methods for deriving the governing equation for a chosendependent variable for the constant coefficient ordinary-DAEs and partial-DAEs as well as special casesof variable coefficient DAEs. The methods can be automated with symbolic coding environments therebyallowing for dealing with systems of any size while retaining analytic nature. This is relevant for interpretablemodeling, analytic characterization and estimation, and engineering design in which the objective is to tuneparameter values to achieve specific behavior. Such insights cannot directly be obtained using numericalsimulations.
Keywords:
Differential algebraic equations, computational engineering, dynamic systems, governingequations, symbolic computation, differential operator theory Correspondence to [email protected], [email protected]
Preprint submitted to Applied Mathematics and Computation February 5, 2021 . Introduction
Many systems in engineering, science, and social science are described by fundamental principles and laws.These laws are expressed as simultaneous equations in terms of multiple variables dependent on time and / orspace. We refer to these equations as constitutive equations. These multivariate equations relate the variousdependent variables of a system and, together, fully describe said system. By the nature of the laws resultingin constitutive equations, the problem is typically perfectly determinable: There are as many equations asthere are unknowns (dependent variables) and all equations are independent of one another. An equationis not a combination of the other equations but rather arises from a distinct law. The set of constitutiveequations that describe a system is referred to as the system of constitutive equations. The systems we areinterested in consist of ordinary differential, partial differential and/or algebraic equations. The differentialequations in these systems may be constant or nonconstant coefficient, homogeneous or nonhomogeneous,and of any order. Some such systems are linear, and others are nonlinear (potentially linearized). The systemmay have one or multiple independent variables. Constitutive systems of differential-algebraic equations arisein a variety of problems such as multibody mechanics [1, 2], membrane distillation [3, 4], chemical kinetics[5, 6] and circuit modified node analysis and other network equation methods [7, 8, 9, 10]. Therefore, itis important to develop analytical (and numerical) methods for studying these systems. Domain expertsmay be interested in: determining solutions via analytic methods or numerical simulation, interpretationsand characterization through analytic study such as determining possible behaviors or parametric regionsof operation, studying control aspects of these systems, or parameter estimation problems.We are particularly interested in developing analytic methods for: (a) obtaining analytic expressions forsolutions to the system, as well as (b) deriving governing equations - which may be thought of as decouplingor reducing to a single dependent variable.Reduction of a system of constitutive equations into a governing equation for a given dependent variableproduces a single univariate ODE or PDE (uni-ODE or uni-PDE) that is of equal or higher order compared tothe constitutive equations . Examples of governing equations derived from constitutive equations include thewave equation, Helmholtz equation, and the governing equation describing the displacement of a particularmass in a vibrating system of masses.The structure of the governing equation is an extremely powerful analytic tool for studying a systemsystems. Consider, for instance, the classification of governing equations that are partial differential equa-tions into elliptic/parabolic/hyperbolic. A system-specific example for the benefit of studying the governingequation is the well-known Webster horn equation which is derived from the constitutive equations in orderto study a horn’s properties which is of relevance for a variety of systems. Another particularly simple exam-ple is that combining equations to arrive at the governing equation for ventilator-respiratory system is keyto developing a noninvasive estimation method for determining patient respiratory system parameters andflow variables from measurements at the ventilator [11]. These examples motivate the need for developing asystematic method for deriving the governing equation from the system of constitutive equations. Despitethis, the derivation of the governing equation from differential-algebraic equations is typically done by handand is therefore only derived if the number of constitutive equations is very small (generally consisting of2-3 dependent variables).The form of the analytic solutions to constitutive systems of differential-algebraic equations (sys-DAEs)is also key to studying these systems. Current methods to obtain analytic solutions exist for certain classesof sys-DAEs. These require first formulating the system using the general state-space formulation . Theformulation does not retain the natural structure of the system but requires introducing state variables.Current analytic solution methods are based on the state-space formulation, and a survey of these solution A governing equation is specific to a system of constitutive equations and a single chosen dependent variable. In the caseof homogeneous constant coefficient systems, the governing equation is the same for all dependent variables Alternatively, the solution for a single dependent variable of the system is determined from the governing equation after it isderived from the constitutive equations by hand. In this case, single equation (uni-ODE) methods such as those of Appendix Bare used.
We use the following abbreviations for univariate equations (uni) and multivariate systems of equations(sys) in this paper. All equations may be homogeneous or nonhomogeneous. As mentioned previously, weare only concerned with linear (L) equations and cases with scalar differential operators and do not considervector differential operators. All systems we consider consist of independent equations which are equal innumber to the number of dependent variables.
Abbreviation Definition uni-LAE Algebraic equationuni-LODE Ordinary differential equationuni-LPDE Partial differential equationsys-LAE System of algebraic equationssys-LODEs System of ordinary differential equationssys-LPDEs System of partial differential equationssys-LODAEs System of ordinary differential-algebraic equationssys-LPDAEs System of partial differential-algebraic equationssys-LDAEs sys-LODAEs and sys-LPDAEscc (prefix, not needed for LAEs) Constant coefficient - e.g. uni-ccLODEsvc (prefix, does not apply to LAEs) Variable coefficient - e.g. sys-vcLODAEs
Table 1: Abbreviations. Note that in the literature the abbreviation DAE usually refers to sys-ODAEs.
In Section 1.1, we motivated the need for methods for analytically solving a system of constitutiveequations as well as methods for deriving governing equations from systems of constitutive equations. Moreexplicitly, our objectives are to develop the following:1. Objective I: Formulation for systems of constitutive equations that are sys-LDAEs.2. Objective II: Systematic methods for deriving governing equations (uni-LODE / uni-LPDE) fromsys-LDAEs for a chosen dependent variable.3. Objective III: Systematic methods for deriving analytic expressions for full solutions of dependentvariables from sys-LDAEs. 3he constant coefficient class of systems is the natural starting point to develop any analytic methodsas evidenced by the literature. Hence, over the course of this paper, we focus on developing the formulation(objective I) and simple analytic methods (objectives II, III) for homogeneous or nonhomogeneous sys-ccLODAEs that include multivariate algebraic equations and ODEs of any order. We then discuss methodsand limitations for other types of sys-LDAEs (sys-ccLPDAEs, sys-vcLODAEs, and sys-vcLPDAEs).The constitutive equations can unambiguously and directly be represented using the proposed formu-lation - as opposed to having to define state variables. The formulation further directly admits algebriacequations and allows for partial differential equations. We use the formulation to develop methods for de-riving the governing equation and obtaining analytic solutions to the system in its native domains. Themethods are suitable for systems of any size and can be programmed into an symbolic system such as Mapleand Mathematica using existing linear algebraic algorithms. The simplicity and the ability to encode theproposed analytic methods in symbolic toolboxes enables them to be used by domain experts who encountersys-LDAEs in their systems (for example in certain distillation equations or linear circuit modified nodalanalysis) rather than being limited in use to those specialized in mathematics or systems theory.
In developing the formulation and methods towards the objectives in Section 1.3, we build on conceptsfrom linear algebra and differential operator theory. Due to this, and the fact that the methods introducedin this paper are intended to be used in a variety of fields by engineers, scientists, and social scientists,we review differential operator theory concepts we use in Appendix B. We refer to equations from theappendix throughout the paper and the reader is encouraged to start with Appendix B if unfamiliar withthese concepts. Also included in Appendix B, are some modifications we included in differential operatortheory that allow our methods to be more direct and general than they otherwise would be.We introduce our formulation and methods starting with sys-ccLODAEs in Sections 2, 3, and 4. Wethen introduce the formulation, methods, and limitations for sys-ccLPDAEs in Section 5, and finally forsys-vcLODAEs and sys-vcLPDAEs in Section 6. Throughout this paper, we provide many examples so thatthe formulation and methods can be better understood and applied. The examples include transmissionlines and vibrating mass systems.
2. Formulating Constitutive sys-ccLODAEs as a Linear System of Equations
We introduce a formulation for the sys-LDAE system of constitutive equations (objective I), and focushere, in particular, on the sys-ccLODAEs. The formulation is specifically constructed for developing theanalytic methods for objectives II and III (deriving the governing equation and obtaining analytic solu-tions) rather than, for instance, numerical simulations, study of properties from controls perspectives, ordetermining eigenvalues and eigenvectors for certain classes of equations.
Consider a sys-LAE, Ax = b (1) l · · · l n ... . . . ... l n · · · l nn x ... x n = b ... b n (2)where we are only concerned with the case of a nonsingular square matrix A which is the categoryrelevant for constitutive equations.Our formulation for the general class of sys-LDAEs involves expressing the system of constitutive equa-tions as if it were a sys-LAE with each row corresponding to a constitutive equation. For sys-ccLODAEs,this is, 4 ( D ) x ( t ) = f ( t ) (3)or, more explicitly, P ( D ) · · · P n ( D )... . . . ... P n ( D ) · · · P nn ( D ) x ( t )... x n ( t ) = f ( t )... f n ( t ) (4)with t as the independent variable. M ( D ) is a matrix where the entries are polynomials of differentialoperator D ddt ; x ( t ) is a vector of time-dependent variables, x i ( t ); and f ( t ) is a vector of forcing terms f i ( t ).By nature of the set of the inherent laws / constitutive equations for a typical system of interest, thematrix of operators will always be square and the rows are all independent. Note that AEs are included inthis form by assigning all polynomials in the corresponding row to be of zeroth order. We also note thatmany integral and integro-differential constitutive equations may be formulated as equivalent differentialconstitutive equations by differentiation (though this may lead to some loss of information). Therefore, ourformulation as a system of linear equations - objective I, and methods for objectives II and III also extendto systems that include constitutive integro-differential equations that can be formulated accordingly (bymultiplication of the corresponding row with some D k ). Any df i dt terms are absorbed into f and there is noneed for defining ant new intermediate variables.Our formulation is appropriate for, and is the first step for, deriving the governing equation (objectiveII) for constitutive sys-ccLODAEs, sys-ccLPDAEs and certain classes of sys-vcLODAEs and sys-vcLPDAEs;deriving an analytic expression for the particular and full solutions that does not require simulation (objectiveIII) as we discuss for ccLODAEs and certain classes of ccLPDAEs; and analytically characterizing the systemand obtaining physical insights. The methods are appropriate for the simple, systematic and automatableanalytic study of systems of any size. In this section, we exemplify our formulation for sys-ccLODAEs. We include the equivalent state-spaceformulation for these examples.For our first example, consider the following set of constitutive equations, ¨ x ( t ) − x ( t ) = f ( t )˙ x ( t ) − x ( t ) = f ( t )˙ x ( t ) + x ( t ) − x ( t ) + x ( t ) = f ( t ) (5)Our formulation (Equation 3) follows directly from the above and is, D − D − − D + 1 x x x = f f f (6)Such a class of systems also has an easily obtainable corresponding state-space formulation of equationA.1. If the system contains an equation with a second order derivative in one of the dependent variables,then, for the state-space formulation, a new dependent variable is defined and the defining equation addedto the system of equations. For the above example, this leads to,˙ x x x x = − − x x x x + f f f (7)5or our second example, consider the following set of constitutive equations, ¨ x ( t ) − x ( t ) + x ( t ) = f ( t ) − x ( t ) + ˙ x ( t ) + ˙ x ( t ) = f ( t )˙ x ( t ) + x ( t ) + x ( t ) = f ( t ) x ( t ) − x ( t ) − x ( t ) + x ( t ) = f ( t ) (8)Our formulation (Equation 3) which follows directly from the above system is, D − − D DD + 1 11 − − x x x x = f f f f (9)The formulation for the same set of equations using singular state-space (Equation A.2) is, ˙ x x x x x = − − − − − x x x x x + f f f f (10)The matrix multiplying ˙ x ( t ) encodes the linear algebraic equation in the final row, and the constantcoefficient ODE that contain derivatives in more than one dependent variable in the third row.The choice of formulation is the basis for developing methods for objectives II and III. Notice thatthe proposed formulation follows unambiguously from the constitutive equations. There is no introductionof intermediate variables (state variables) to study a system with higher order derivatives. The proposedformulation does not involve singular matrices which are traditionally needed for algebraic equations andequations with derivatives in more than one state variable. As discussed in Section 5, the formulation alsoadmits partial differential equations.In Appendix A, we review existing methods for objective III that are for the general state-space formu-lation. The rest of the paper is concerned with developing methods for objectives II and III based on ourdirect formulation.
3. Methods for Deriving Governing uni-ccLODE from Constitutive sys-ccLODAEs
In this section, we develop methods for objective II - systematically deriving the higher (or equal)order governing uni-ccLODE for a particular dependent variable from a system of constitutive equationssys-ccLODAEs formulated as in Equation 3.
In order to describe our method for objective II, consider Gaussian elimination for sys-LAEs (Equation2). Gaussian elimination involves a series of row operations which leads to a final row, h ( l ) x n = n X i =1 g i ( l ) b i (11)where l is a vector containing the entries l ij of matrix A . And h and g i are functions of their argumentsin a form dictated by Gaussian elimination row operations. Row permutation is occasionally needed duringGaussian elimination. 6o derive higher order governing uni-LODEs/uni-LPDEs in a specified dependent variable, we use Gaus-sian elimination. Let x n ( t ) be the variable of interest for the governing equation. We formulate the systemas in Equation 3 with the variable of interest as the final variable .We then carry out steps of Gaussian elimination on the system, which yields a parallel to Equation 11, h ( ˜ P ) x n ( t ) = n X i =1 g i ( ˜ P ) f i ( t ) (12)where ˜ P is a vector of the polynomial operators P ij ( D ) which are entries of matrix M ( D ) of Equation 3. h and g i have forms arising from row operations. Note that we do not need to keep track of the order ofoperations of Gaussian elimination to reach Equation 12 for the case of sys-ccLODAEs and they may beconsidered commutative because the equations are constant coefficient. This arises because Dax = aDx .However, this will not generally be the case for sys-vcLODAEs and sys-vcLPDAEs as we discuss in Section6.2.We then manipulate this expression into the governing uni-ccLODE by multiplication of both sides by alldenominators of the n th rows of the elimination matrices used in Gaussian elimination to obtain Equation12 to yield the following form. This manipulation via multiplication is possible owing to the constantcoefficients nature. P ( D ) x n ( t ) = n X i =1 Q i ( D ) f i ( t ) (13)Where P ( D ) and all Q i ( D ) are polynomials in the differential operator, D .Equivalently, we may simply write, P ( D ) x n ( t ) = φ ( t ) (14)where P ( D ) is a polynomial of differential operators of equal or higher degree than the individualconstitutive equations and φ ( t ) is an effective forcing term. This equation is clearly in the form of a higherorder governing uni-ccLODE, thereby achieving objective II for sys-ccLODAEs.Note that if we only need to determine P ( D ) - which is the case of the homogeneous sys-ccLODAEs,then we may utilize the constant coefficient nature and alternatively derive it using P ( D ) = det M ( D ). Therelationship to the the method described above can be clearly inferred from the fact that the determinantis the multiplication of pivots which are in turn also involved in Gaussian elimination and multiplication ofdenominators of the n th rows of the elimination matrices.Because the proposed method for objective II builds on fundamental concepts of linear algebra, existingalgorithms may be utilized for symbolic code implementations. The methods for sys-ccLODAEs are suitable for systems of any size. However, many systems, includingclassically studied systems, have two or three dependent variables and corresponding constitutive equations.We therefore include the specific form for the corresponding higher order governing equations (Equation 14)here.
For a 2x2 sys-ccLODAE, the result of Gaussian elimination (Equation 12) is, (cid:18) P ( D ) − P ( D ) P ( D ) P ( D ) (cid:19) x ( t ) = − P ( D ) P ( D ) f ( t ) + f ( t ) (15) Related to this is the fact that the higher order governing equation for homogeneous sys-ccLODAEs is the same for anychosen dependent variable. However, we explicitly state choosing the dependent variable to be the final variable in the vector x here as it is necessary to specify for the nonhomogeneous case, for the case in which M ( D ) of Equation 3 can be written asblock diagonal, and for the variable coefficient case. P ( D ), the denominator of last (2 nd ) row of the eliminationmatrix involved in Gaussian elimination to generate h ( ˜ P ), yields the higher order governing equation for x . Equation 14 for 2x2 system is, (cid:18) P ( D ) P ( D ) − P ( D ) P ( D ) (cid:19)| {z } P ( D ) x ( t ) = − P ( D ) f ( t ) + P ( D ) f ( t ) | {z } φ ( t ) (16)The proposed method is in fact a systemization of the elimination methods done by hand to derive thegoverning equation for 2x2 systems. For the three variable system, we get for Equation 12. h ( ˜ P ) = P − P P P − P − P P P P − P P P (cid:0) P − P P P (cid:1) (17)with P ij being polynomials P ij ( D ). We have dropped explicitly writing D for compactness.and, n X i =1 g i ( ˜ P ) f i ( t ) = − P P − P − P P P P − P P P (cid:18) − P P (cid:19)! f − P − P P P P − P P P f + f (18)where f i = f i ( t ) and we have again dropped explicitly writing the dependence on t for compactness.By multiplication by denominators of the third rows of the elimination matrices used in Gaussian elimi-nation, P (cid:0) P − P P P (cid:1) , we get for Equation 14, the (potentially higher order) governing uni-ccLODE, P ( D ) = ( P P − P P ) P + P P P + P ( P P − P P ) − P P P (19)which may also have been obtained directly using the determinant approach.And, φ ( t ) = (cid:0) P P − P P (cid:1) f ( t ) − ( P P − P P ) f ( t ) + (cid:0) P P − P P (cid:1) f ( t ) (20) In this section, we illustrate the method for deriving the governing equation from sys-ccLODAEs usingan example motivated by classic electronic traveling-wave amplifiers from microwave theory and an earlyphenomenological model of cochlear amplification [12, 13]. This system consists of two coupled continuoustransmission lines and we use the constitutive equations in the frequency domain as an example. Derivingthe governing equation is useful in determining the class of possible behaviors, characterizing the form, andparametric regions of solutions of interest for such a system e.g. regarding power transfer between the linesand determining which parameter choices provide desired amplification patterns.We choose the circuit configuration drawn in [12] and express constitutive equations relating the depen-dent variables accordingly. In this system, the top line has variables I t ( x, ω ) , V t ( x, ω ), and the bottom linehas variables, I b ( x, ω ) , V b ( x, ω ) . We then express these in the form of Equation 3 as, The identity of the specific parameters and variables are not important - as opposed to the form of these relationshipsto illustrate the systematic method of deriving the governing equation. In the case of the chosen configuration, Z st ( ω ) = jωL ot , Z sb ( ω ) = jωL ob , Y pt ( ω ) = jωL t + R t + jωC t , Z pb ( ω ) = jωC b + R b . The effective sources from both lines are due tothe voltage across the capacitance. Z st ( ω ) DZ pt ( ω ) D − αZ sb ( ω ) Dβ ( ω ) D D Y bp ( ω ) I t ( x, ω ) V t ( x, ω ) I b ( x, ω ) V b ( x, ω ) = 0 (21)where the coupling between the two transmission lines is due to α and β ( ω ), and D ddx is the spatialderivative operator. Notice that setting the coupling parameters to zero leads to separation into two blockmatrices and, in that case, carrying out the methods of Section 3.1 leads to two independent wave equations(one for the dependent variables of the top line and another for the dependent variables of the bottom line).Applying the method described in Section 3.1 with V b ( ω ) as the dependent variable of interest , givesafter Gaussian elimination (Equation 12), (cid:18) Y bp − αβ ( ω ) D Z st − Z pt D − D Z sb (cid:19) V b ( x, ω ) = 0 (22)Multiplication by the relevant denominators of the elimination matrices Z sb ( ω ) (cid:18) Z st ( ω ) − Z pt ( ω ) D (cid:19) ,yields the higher order governing equation (Equation 14), D − (cid:18) Z st ( ω ) Z pt ( ω ) + αβ ( ω ) Z sb ( ω ) Z pt ( ω ) + Z sb ( ω ) Y bp ( ω ) (cid:19) D + Z st ( ω ) Z sb ( ω ) Y bp ( ω ) Z pt ( ω ) ! V b ( x, ω ) = 0 (23)which is conducive to studying the system. Note that we may have arrived at this equation directly fromthe sys-ccLODAEs using P ( D ) =det M ( D ) as described in Section 3.1 because the system is homogeneous,but sought to exemplify the method for the more general case.To return to the block diagonal case: in the case of α = β = 0, the P ( D ) of the above equation may befactored into two second order differential operator polynomials - or two wave equations with wavenumbers q Z st ( ω ) Z pt ( ω ) , p Z sb ( ω ) Y sb ( ω ). Each of which provides two solutions. The boundary conditions impose that oneset of solutions describes behavior on one transmission line and the other set describes behavior on the othertransmission line. In this section, we illustrate that our method for deriving the governing equation applies not only tosystems of sys-ccLODEs but also to sys-ccLODAEs which also occur naturally in systems and phenomenon- though they are less-often the focus of development of even numerical methods. The sys-ccLODAE isthe constant coefficient version of a classical model of auditory filters - e.g. [14], that may be appropriatefor design purposes (such as machine hearing) rather than a physiological model of the cochlea. We define D ddx and choose the variable of interest for this sys-ccLODAE to be V bm ( x, ω ). Expressed as Equation3, the system is, D Z ( ω ) D b bZ fbal ( ω ) P ( x, ω ) U ( x, ω ) V bm ( x, ω ) = (24)From Equation 19 (and division by bZ fbal ), we obtain the governing equation, (cid:18) D + ZZ fbal (cid:19) V bm ( x, ω ) = 0 (25) In this case, the governing equation is the same for all dependent variables because the system is homogeneous and constantcoefficient. If the coupling parameters are set to zero, the lowest-order governing equation is different for each of the lines. We also divide by Z pt ( ω ) to put the governing equation in the standard ODE form in which the coefficient of the highestorder differential operator is unity. . Methods for Solving sys-ccLODAEs In this section, we develop methods for objective III - deriving analytic expressions for the full solutionsfor a system of multivariate constitutive equations that form a sys-ccLODAE. The problem is solved in thenative time or space domain of the problem. In Appendix B, we review differential operator methods foranalytically solving uni-ccLODEs and make some modifications. Those methods are utilized as part of ourapproach for solving sys-ccLODAEs.
For objective III, we again use concepts from linear algebra. Recall that to solve a sys-LAE as in Equation2, we use inversion of the nonsingular matrix A , which results in solutions for all x i , x i = n X j =1 h j ; i ( l ) b j (26)The Gauss-Jordan method for inversion takes an algebraic system of equations, applies row operationsto both sides, up to arriving at pivot columns of A . For the case of the full-rank matrix A , after row-normalization, this leads to the identity matrix, which is then used to determine the solution of each ofthe variables in the column vector x . Inversion, and particularly the Gauss-Jordan method is relevant fordetermining the full solutions for each of the dependent variables in a system of sys-LDAEs. The resultantexpressions are analytic (but not necessarily closed-form).The key to developing our methods for solving sys-ccLODAEs is to formulate the sys-ccLODAEs as inEquation 3. If the only interest is in the complementary solutions (which are the same for all dependentvariables for sys-ccLODAEs) or the system is homogeneous, then the solution can be obtained by determiningthe roots of the characteristic equation, det M ( a ) = 0 (27)If we are interested in the full solution (including complementary solutions and a particular) to a non-homogeneous sys-ccLODAE, then we invert the matrix of operators. Due to the constant nature of thecoefficients, we may define the inverse of the matrix of operators M ( D ) by extrapolating standard matrixinversion methods.If the inverse is obtained using Gauss-Jordan elimination, then these operations result in an equationparallel in form to Equation 26, x i ( t ) = n X j =1 h j ; i ( D ) f j ( t ) (28)where h j ; i ( D ) has a form due to the row operations for inversion described by Gauss-Jordan elimination.We then multiply each h j ; i by fractions of equal numerator and denominator (unity) repeatedly to put x i ( t )in the form, x i ( t ) = n X j =1 P j ; i ( D ) Q j ; i ( D ) f j ( t ) = n X j =1 Q j ; i ( D ) φ j ; i ( t ) (29)which can then be solved for each of the dependent variables using differential operator methods foruni-ccLODEs described in Appendix B (Equation B.11 or B.12) with Definition B.8 for the full solution orDefinition B.7 for the particular solutions only.The methods are appropriate for sys-ccLODAEs of any size, and are particularly simple and intuitivedue to their relation to standard matrix inversion (which is possible due to the constant coefficient natureof the equations). Because the methods we develop in this paper build on fundamental concepts of linearalgebra, existing algorithms may be utilized for symbolic code implementations. However, we note thatmethods for objective II only requires elimination and are therefore more computationally efficient thanthose for objective III which involves matrix inversion.10he methods for objective III provide analytic solutions which are useful in studying the system in waysthat cannot be achieved if we are limited to numerical simulation of systems of differential equations whichsuffer from issues such as accuracy, stability, specificity to parameters, boundary and initial conditions, anddiscretization schemes, and limitations to types of equations. Note that the solution method does not requiretransforms to move to a different domain but rather solves the problem in its native domains. The analyticmethod is not preferential to any particular type of problem (e.g. IVPs or BVPs). To exemplify our method for solving sys-ccLODAEs, consider a typical mass-spring system with twomasses m , m , connected via a spring with constant k . The second mass is connected to a fixed wallwith a spring with constant k . An external upward force f ( t ) is applied to the m , and the upwarddisplacements of the two bodies are x ( t ) , x ( t ).The system of constitutive equations is, ( m ¨ x ( t ) + k x ( t ) − k x ( t ) = f ( t ) m ¨ x ( t ) − k x ( t ) + ( k + k ) x ( t ) = 0 (30)Our proposed method for solving sys-ccLODAEs is applicable to vibration systems of any configurationand admits forcing terms. It is not restricted to pure mass-spring systems but can include dashpots. Tosolve for x ( t ) , x ( t ) using our method, we first formulate the system as in Equation 4. With D ddt , weexpress this system as, (cid:20) m D + k − k − k m D + k + k (cid:21) (cid:20) x ( t ) x ( t ) (cid:21) = (cid:20) f ( t )0 (cid:21) (31)Inverting for x ( t ) as described in Section 4, we obtain, x ( t ) = 1( m D + k )( m D + k + k ) − k (cid:20) m D + k + k k k m D + k (cid:21) (cid:20) f ( t )0 (cid:21) (32)We can solve this for each dependent variable by using the methods described in more detail in Appendix B.Factoring the polynomial in the denominator as1( m D + k )( m D + k + k ) − k = 1 m m D − a )( D + a )( D − b )( D + b ) (33)with a, b = i √ r ( k m + k m + k m ) ± q k m + k m + k m + k k m + k m m − k k m m , we then obtain analyticexpressions for the full solution for each of the dependent variables x ( t ) , x ( t ) . As there are no repeatedfirst order factors of the polynomial in the denominator, we may use either Equation B.11 or Equation B.12to obtain these expressions (only Equation B.12 is appropriate in the case of repeated first order factors -see Appendix B). For example, using Equation B.11 to solve for x ( t ), we obtain, x ( t ) = 1( m D + k )( m D + k + k ) − k k f ( t ) (34)= k m m a − b a e at Z e − at f ( t ) dt − a e − at Z e at f ( t ) dt − b e bt Z e − bt f ( t ) dt + 1 b e − bt Z e bt f ( t ) dt ! (35) Note that this does not provide a relationship between the constants of the complementary components of the two dependentvariables. x ( t ) = 1( m D + k )( m D + k + k ) − k k f ( t ) (36)= k m m e at Z e − at Z e ( a + b ) t Z e − bt Z e bt f ( t ) dtdtdtdt (37)For instance, for f ( t ) = f , we get the full solution, using either expression, x ( t ) = k m m a b f + C e at + C e − at + C e bt + C e − bt (38)which clearly includes the complementary solutions and a particular solution. The method yields bothparticular and complementary parts of the solution together rather than solving for each separately.
5. Methods and Limitations for sys-ccLPDAEs
In this section, we present our methods for objectives I-III (formulation, governing equation derivation,and system solutions) starting from systems of constitutive equations that are not sys-ccLODAEs but rathersys-ccLPDAEs where we have multiple independent variables and the corresponding governing equation is auni-ccLPDE as opposed to a uni-ccLODE. In order to extend our methods and determine their applicabilityto, and limitation for, sys-ccLPDAEs, we must first extend differential operator theory to the multi-operatorcase.
We may define each differential operator D i to be associated with an independent variable t i , such that, D i x ( t , t , . . . t m ) , ∂∂t i x ( t , t , . . . t m ) (39)Consequently, we may generally express a uni-ccLPDE as, P ( D , D , . . . D m ) x ( t , t , . . . t m ) = φ ( t , . . . t m ) (40)with P ( D , . . . D m ) being a polynomial in the differential operators.The multi-operator counterparts to the properties of Appendix B apply for any D i .Accordingly, for instance, P ( D , D , . . . D m ) e a t + a t + ...a m t m = P ( a , a , . . . a m ) e a t + a t + ...a m t m (41)Note that this property can be used to simply find the relationship between the different components ofthe complementary solution to a uni-ccLPDE by setting P ( a , a , . . . a m ) = 0. This is the counterpart tothe characteristic equation obtained from a uni-ccLODE (see Equation B.4).Similar to Equation B.13, the order of integration and differentiation does not matter when applied tosmooth functions regardless of the independent variables and whether they are the same or mixed indepen-dent variables. This is due to the constant coefficient nature of the system. To formulate the sys-ccLPDAEs (objective I), we express the system of equations as in Equation 3,but with t a vector of n independent variables t , t , · · · t n (e.g. time and spatial coordinates), and D avector of differential operators. Accordingly, the uni-operator polynomials P ij ( D ) of sys-ccLODAEs arereplaced by multi-operator polynomials in multiple differential operators P ij ( D ) = P ij ( D t , D t , · · · D t n ) forsys-ccLPDAEs. We therefore express the counterpart of Equation 3 for sys-ccLPDAEs as,12 ( D , D , · · · D t n ) x ( t , t , · · · t n ) = f ( t , t , · · · t n ) (42)or, compactly, M ( D ) x ( t ) = f ( t ) (43)with M ( D ) the matrix of operators, x ( t ) the vector of dependent variables, and f ( t ) the vector of forcingterms. The method for deriving higher order governing equations (objective II) for constitutive sys-ccLPDAEs,transfers directly from sys-ccLODAEs (Section 3) with no limitations or approximations but with the afore-mentioned changes to the independent variables and operators.
In this section, we exemplify the method for deriving the higher order governing equation from a systemof constitutive sys-ccLPDAEs. For our example, we consider a transmission line (which is of interest inmany propagation systems [15, 16, 17]) in the time domain. A commonly encountered transmission line hasan inductor with inductance L δx for the series component, and a simple resonant harmonic oscillator orRLC in series for the parallel component . The constitutive equations arise from Kirchhoff’s current law(continuity) and Kirchoff’s voltage law (force balance), and the forcing terms are chosen to be zero in thisexample. The system responds to a set of initial and boundary conditions. The differential equations aredefined by taking the limit of differential equations as δx −→
0. The set of constitutive equation for thissystem is, ( L ∂ ∂t∂x i ( x, t ) + R ∂∂x i ( x, t ) + C R t −∞ ∂∂x i ( x, τ ) dτ + v ( x, t ) = 0 L ∂∂t i ( x, t ) + ∂∂x v ( x, t ) = 0 (44)We use D t and D x to denote time and space differential operators, D t ∂∂t , D x ∂∂x . We express theconstitutive equations as a linear system. This is a system in which the integro-differential equation can bereformulated as a differential equation simply by multiplication of the first row by D t to use the methodsdeveloped in this paper. This yields for objective I (Equation 43), (cid:20) ( LD t + RD t + C ) D x L D t D x (cid:21) (cid:20) i ( x, t ) v ( x, t ) (cid:21) = (cid:20) (cid:21) (45)We then apply our method to derive the governing equation (objective II) for v ( x, t ). Using Equation16, we obtain, (cid:18)(cid:0) LD t + RD t + 1 C (cid:1) D x − L D t (cid:19) v ( x, t ) = 0 (46)which can then be characterized according to its parameter values which may be chosen or controlledfor design problems.As mentioned previously, derivation of the governing equation allows us to characterize the system anddetermine parametric conditions for various behaviors without the need for solving the constitutive equations- e.g. conditions under which wave propagation is supported and form of exchange of energy between thecoupled two transmission line system of Section 3.3. Trivially, setting R = L = 0 leads to the simple 1D wave equation as becomes apparent in Equation 46. Addition ofnonzero R, L adds dissipative (or, depending on the sign, growing) and dispersive behavior. .4. Solving for Dependent Variables: Methods, Approximations, and Limitations In this section, we discuss methods and limitations for objective III - obtaining exact or approximateanalytic solution, for sys-ccLPDAEs. We first formulate the sys-ccLPDAEs as in Equation 43 and thencarry out the methods of Section 4 to arrive at the multi-operator counterpart to Equation 29. This resultsin, for each dependent variable x i ( t ), x i ( t ) = n X j =1 Q j ; i ( D ) φ j ; i ( t ) (47)This equation cannot generally be solved using methods of differential operator theory. The issue is notdue to dealing with a system, but is rather inherent in partial differential equations. To our knowledge,there are no general analytic solutions for nonhomogeneous uni-ccLPDEs (Equation 40) . Certain classicaluni-ccLPDEs (predominantly homogeneous and without mixed derivatives) can be solved using the methodof separation of variables.We may solve a uni-ccLPDE if we can, exactly or approximately, factor the multi-operator polynomialinto uni-operator polynomials (which may require transformations) - i.e. we can express P ( D ,D , ··· ,D n ) ≈ P ( ˜ D ) 1 P ( ˜ D ) · · · P m ( ˜ D m ) . This allows us to then apply equations B.11 or B.12 to arrive at the analyticsolution for such a uni-ccLPDE.Consequently, we may express the analytic solutions for sys-ccLPDAEs for cases of naturally weakcoupling in the native independent variables as, x i ( t ) ≈ n X j =1 m Y k =1 Q j ; i ; k ( D k ) φ j ; i ( t ) (48)which can then be solved analytically using equations B.11 or B.12. This requires factoring the multi-operator polynomials Q j ; i ( D ) of Equation 47 such that the absolutely irreducible factors are uni-operatorpolynomials (i.e. each factor only contains D k for a single k ).Only a subset of sys-ccLPDAEs naturally have weak coupling in their independent variables. Therefore,in order to factor Q j ; i ( D ) into uni-operator polynomials, we may additionally require manipulating theequation using transformations of the independent and/or dependent variables.˜ x i ( ˜t ) ≈ N X j =1 M Y k =1 Q j ; i ; k ( ˜ D k ) ˜ φ j ; i ( ˜t ) (49)
6. On sys-vcLODAEs and sys-vcLPDAEs
In this section, we present a formulation for sys-vcLODAEs which is also extendable to sys-vcLPDAEs(objective I). We then consider special cases of sys-vcLODAEs for which we can derive the governing equation(objective II) using the same methods as sys-ccLODAEs. This also applies to the corresponding special casesof sys-vcLPDAEs. In this paper, we do not discuss methods for determining solutions to sys-vcLODAEsand sys-vcLPDAEs, and it is generally not possible to provide exact analytic solutions for these categories.
We may formulate the sys-vcLODAEs as follows, M ( D t , t ) x ( t ) = f ( t ) (50) This can be seen by the fact that unlike P ( D ) of Appendix B, P ( D ,D , ··· ,D n ) is not defined as an operator. For instance,we cannot solve an equation ( D x + D t ) u ( x, t ) = f ( x, t ), by simply expressing it as, u ( x, t ) = D x + D t f ( x, t ) D t ddt . The elements of the matrix M are H ij ( D t , t ), which take the form, H ( D t , t ) = n X k =0 a k ( t ) D kt (51)Note that sys-ccLODAEs are special cases of sys-vcLODAEs when the matrix elements H ij ( D t , t ) = P ij ( D t ).As a particularly simple example, consider a continuous transmission line with incremental sources andin which the series components are constant along x but the parallel components vary with x . The systemof equations (sys-vcLODAEs) is formulated as, (cid:20) Z ( ω ) DD Y ( x, ω ) (cid:21) (cid:20) I ( x, ω ) V ( x, ω ) (cid:21) = (cid:20) f ( x, ω ) f ( x, ω ) (cid:21) (52) In this section, we discuss methods for deriving the governing equation (objective II) for special casesof sys-vcLODAEs and sys-vcLPDAEs. We cannot generally extrapolate the methods we developed forderiving the governing equation for sys-ccLODAEs to the case of sys-vcLODAEs. Two steps of the methodsfor sys-ccLODAEs account for this:1. The order of operations of Gaussian elimination to arrive at Equation 12 does not matter in the caseof constant coefficients, but does matter for variable coefficients.2. Multiplication of Equation 12 by denominators results in the governing uni-ccLODE (Equation 14) forsys-ccLODAEs because multiplication and division are interchangable for the sys-ccLODAEs. Howeverthis is not the case for general sys-vcLODAEs.There is a subset of sys-vcLODAEs (Equation 50) for which these two issues do not arise. For these cases,we simply use the methods we constructed for sys-ccLODAEs (Section 3) for deriving the correspondinggoverning equation. The criteria is that the sys-vcLODAEs are of the form, P , ( D t ) · · · P ,n − ( D t ) H ,n ( D t , t )... . . . ... ... P n, ( D t ) · · · P n,n − ( D t ) H n,n ( D t , t ) x ( t )... x n ( t ) = f ( t )... f n ( t ) (53)In other words, only the final column - operating on x n ( t ), contains the variable coefficients. In otherwords, all variable coefficients in the system of equations must correspond to the same dependent variable.For such systems, issues (1) and (2) do not arise because the corresponding elimination matrices used forGaussian elimination contain only constant coefficient terms (the elements in the final column of M donot appear as part of the elimination matrices).As for sys-vcLPDAEs, the governing equation is derived in a similar manner and under the same conditionas sys-vcLODAEs. Future directions may include methods for deriving the governing equation for moregeneral cases of sys-vcLODAEs and sys-vcLPDAEs in which more than one dependent variable has variablecoefficients. Future directions also include developing methods for deriving approximate expressions for thesolution for each dependent variable in a sys-vcLODAE. Equivalently, the L in LU matrix decomposition contains only constant coefficient terms .2.2. Example: Continuous Transmission Line with Varying Parallel Admittance To give an example, we use the system in Equation 52 and choose V ( x, ω ) to be the dependent variableof the governing equation. Because it satisfies the condition in Equation 53, we may express the governingequation using Equation 16 as, (cid:18) D x − Z ( ω ) Y ( x, ω ) (cid:19) V ( x, ω ) = D x f ( x, ω ) − Z ( ω ) f ( x, ω ) (54)Note that the condition for extrapolating sys-ccLODAE methods for the governing equation is not metfor the same system if the dependent variable of interest is I ( x, ω ). Consequently, if given flexibility regardingthe choice of dependent variable with regards to governing equations, it is advantageous to choose dependentvariables for which the condition (Equation 53) holds.
7. Conclusion
Analytic solutions to the system and analytic characterization through governing equations provideinsights into the behaviors of DAEs as well as the parametric regions of operation for each potential behavior.This is relevant for interpretable modeling, analytic characterization and estimation, and engineering designin which the objective is to tune parameter values to achieve specific behavior. Such insights cannot directlybe obtained using numerical simulations. In this paper, we introduced methods for deriving analytic solutionsfrom DAEs as well as methods for deriving governing equations for analytic characterization.We first introduced a direct formulation for representing systems of linear differential-algebraic scalarconstitutive equations that are homogeneous or nonhomogeneous with constant or nonconstant coefficients.Using this formulation for the system of constitutive equations (DAEs), we (a) developed simple methods forderiving analytic expressions for the full solution (complementary and particular) for all dependent variablesof systems that consist of constant coefficient ordinary-DAEs and special cases of partial-DAEs. We also(b) developed simple methods for deriving the governing equation for a chosen dependent variable for theconstant coefficient ordinary-DAEs and partial-DAEs as well as special cases of variable coefficient DAEs.The proposed formulation to represent systems of DAEs uses differential operators. Unlike state spaceformulations, our formulation follows very directly from the system of constitutive equations, directly han-dles algebraic equations, handles higher order derivatives without requiring defining new variables, and isappropriate for systems containing PDEs. To develop the methods for deriving the governing equationand the analytic solutions from the system of constitutive equations, we utilized the concepts of Gaussianelimination and inversion from linear algebra and (somewhat modified) differential operator theory.The methods are simple, build on fundamental linear algebra and hence existing algorithms may beused, can be automated with symbolic coding environments thereby allowing for dealing with systems ofany size while retaining analytic nature, and are intended for use with a variety of engineering, physical, andfinancial systems. We exemplified our methods using dynamic and propagation systems such as vibratingmasses and coupled continuous transmission lines.
Future directions include addressing limitations of the methods presented here. This includes developingmethods for deriving the exact or approximate governing equation for the general case of sys-vcLDAEs asopposed to only those in which all variable coefficients correspond to a single dependent variable. Thisalso includes methods for deriving approximate analytic expressions for the full solutions for dependentvariables in sys-vcLODAEs. Another limitation is regarding solving the system of sys-ccLPDAEs, andfuture work includes developing approaches for determining transformations and approximate uni-operatorfactorization methods in order to derive analytic solutions for the sys-ccLPDAE. Methods for exact andapproximate factorization of multivariate polynomials exist in the literature but to our knowledge these donot generally result in univariate factors. Instead, they represent or approximate a multivariate polynomialas a multiplication of multivariate (rather than univariate) polynomial factors [18, 19, 20, 21]. Hence,16uture directions include developing systematic methods for constructing the transformations of t , x i , φ j ; i into ˜t , ˜ x i , ˜ φ j ; i as well as constructing approximate factorization methods. Doing so would enable us to obtainapproximate analytic solutions for sys-ccLPDAEs other than those that are weakly coupled in their nativeindependent variables and are for which the factorization is apparent.Future directions also include extending the framework to develop methods for deriving the governingequation and analytic expressions for full solutions for related classes of problems - specifically linear systemswith other operators which share certain properties with differential operators (e.g. linear difference-algebraicequations). As our formulation is in a simple matrix form M x = f , we may, in the future, be able to extendlinear algebra concepts such as factorization and pseudoinverse to M of operator polynomials. Finally, theformulation may also be used for developing numerical solution methods for boundary value problems.17 ppendix A. Classical Formulation and Solution Methods In this appendix, we discuss the main analytic methods for solving some classes of sys-LDAEs. Asmentioned previously, we are only concerned with smooth linear systems.The class of sys-LODEs are the most studied. They are represented using the state-space formulation,which is of the form, ˙ x ( t ) = A ( t ) x ( t ) + B ( t ) u ( t ) (A.1)Constitutive equations are not necessarily native to this form. Instead, state variables are defined sothat constitutive equations of a system can be formulated as in Equation A.1. Note that the formulationis a set of first order ODEs where the derivative is in a single dependent variable. Second order undampedmass-spring systems M ¨ x ( t ) + Kx ( t ) = 0 effectively fall into this category. Current methods for solving oranalyzing these sys-LODEs are built on the state-space formulation.There are several numerical algorithms for solving sys-LODE initial value problems (IVPs) of this for-mulation. The constant coefficient problems, sys-ccLODE, in the state-space form have simple analyticsolutions. In the time-domain, the solution is achieved through matrix exponentials and integrating fac-tors x ( t ) = e At x (0) + R t e A ( t − τ ) Bu ( τ ) dτ . Alternatively, the analytic solution can be achieved by Laplacetransform sX ( s ) − x (0) = AX ( s ) + BU ( s ) which is then solved as an eigenvalue problem or converted intotransfer functions for a single variable and solved accordingly . Note that these are solved in a domainthan the independent variables. Laplace transform methods are used for LTI systems and cannot be usedfor systems with variable coefficients. Beyond analytic solutions for sys-ccLODEs, there is literature on theanalysis of such systems from a controls perspective (controllability, stability, solvability, observability). Thestudy of control properties also builds on the state space formulation and is useful for designing observersand controllers - which are not the objectives of this paper.Variable coefficient sys-LODE IVPs that can be expressed in the form ˙ x ( t ) = A ( t ) x ( t ) + B ( t ) u ( t ); x (0) = x o have a solution x ( t ) = φ ( t, x (0)+ R t φ ( t, τ ) B ( τ ) u ( τ ) dτ . However, φ ( t, τ ) is a problem-specific transitionmatrix that is not always easy to determine. For the case of sys-ccLODEs in the previous paragraph, φ ( t, τ ) = e A ( t − τ ) . Special cases of sys-vcLODEs admit closed-form expressions for φ ( t, τ ), but this isgenerally not the case. The solution to certain homogeneous sys-vcLODE IVPs can be expressed usingMagnus expansion which is useful for analytic study when the series converges [22]. The higher terms inthis series are increasingly complex, and hence it is desirable to approximate by truncation.Sys-LODAEs are also formulated using a more general state-space representation, E ( t ) ˙ x ( t ) = A ( t ) x ( t ) + B ( t ) u ( t ) (A.2)where E ( t ) is a singular matrix due to rows with only zeros corresponding to the algebraic equations fortrue sys-LODAEs . This representation is also referred to as semi-state equations, singular systems, orcontinuous descriptor systems. As is the case with sys-LODEs formulated using state-space, the constitutivelaws are not native to this form, and it is not always easy (or direct) to covert them into the state-spaceformulation. Notice that expressing sys-LODAEs in this particular form leads to dealing with a singularmatrix, E ( t ). This is distinct from, and does not leverage, the fact that the constitutive equations arise frominherently independent physical laws. Consequently, sys-LODAEs formulated using the state-space repre-sentation are often referred to as singular systems (especially in the systems and control theory literature)[23].Certain constant coefficient (sys-ccLODAE) are solvable in the time domain or using Laplace transform.The condition for these systems are those for justifying inverting a resultant matrix pencil sE − A . For Conversion to a transfer function is in some ways similar to deriving a governing equation from the system of constitutiveequations. This is because they are in terms of a single dependent variable, and it is possible to convert some governingequations into transfer functions by taking the Laplace transform after deriving the governing equation from the system ofconstitutive equations sys-LODEs that contain equations with derivatives in more than one dependent variable may also be written implicitly inthis form. In this case, E ( t ) is not singular. .For solution in the time domain, the systems in Equation A.2 are generally transformed into a canonicalform. This reduces the system into two decoupled subsystems in transformed variables. ( ˙ y ( t ) = A y ( t ) + B u ( t ) N ˙ z ( t ) = z ( t ) + B u ( t ); N m = 0 , N m − = 0 (A.3)The above can then be solved analytically for sys-ccLODAEs with sufficiently differentiable u ( t ). Thisis due to the fact that the first equation is a sys-ccLODE, and the second involves a matrix N that has aspecial structure and is nilpotent [25, 26]. The nilpotency index is an example of a differential index whichis roughly associated with its distance from sys-LODEs and is used to describe the difficulty associated withnumerically solving sys-ccLODAEs [10]. Specified initial values must be consistent.In the s -domain, sys-ccLODAE are expressed and solved as generalized eigenvalue problems. Like eigen-value problems associated with sys-ccLODE, there are existing efforts for developing accurate and efficientalgorithms for solving generalized eigenvalue problems which are computationally expensive for large sys-tems [27]. Sys-LODAEs are also studied from a controls perspective for controllability and observability andto design optimal controls [25, 28, 29].Note that existing methods for studying sys-LODEs and sys-LODAEs generally require shifting awayfrom the natural structure of the system. The formulation requires that state variables be defined such thatthe system can be represented as in Equation A.1 or A.2. Appendix B. Differential Operator Theory for uni-ccLODEs
In this section, we review differential operator theory [30] (and introduce modifications) as it is essentialin developing our methods for objectives I-III. Classically, differential operator theory methods are extremelypowerful and are applied to solve uni-ccLODEs - these are of the form, P ( D ) x ( t ) = f ( t ) (B.1)as well as other special uni-ODEs. The complementary solution and particular solution are determinedseparately.Building on its algebraic abstraction and application to solving uni-ccLODEs, we utilize differentialoperator theory to formulate constitutive equations in Section 2, derive governing equations in Section 3and obtain solutions for a system of sys-ccLODAEs in Section 4.The following operator definitions and rules regarding their manipulation to solve a uni-ccLODE (Equa-tion B.1) are described in the literature [30]. We also make a modification to operator definitions to obtainthe full solution more directly, and to more directly handle repeated root cases. • Dx ( t ) , ddt x ( t ) (B.2) • Du ( t ) v ( t ) = ddt (cid:0) u ( t ) v ( t ) (cid:1) = u ( t ) Dv ( t ) + v ( t ) Du ( t ) (B.3) • P ( D ) e at = P ( a ) e at (B.4)where P ( . ) signifies a polynomial in the argument. Classically, this is used to determine the comple-mentary solution of a uni-ccLODE via finding the roots of the characteristic equation P ( a ) = 0. • P ( D ) u ( t ) v ( t ) = u ( t ) P (cid:18) D + u ′ ( t ) u ( t ) (cid:19) v ( t ) (B.5).Of particular interest is, Generalized matrix inverses (e.g. Drazin inverse) are used for some other resultant matrix pencil. De at x ( t ) = e at (cid:0) D + a (cid:1) x ( t ) (B.6) • Definition of D - classically and using our modification: – D f ( t ) , Z tt o f ( τ ) dτ (B.7).Is the classical definition in the literature. The lower limit of integration t o can take on any valueand is chosen based on convenience (e.g. −∞ ). This is because, classically, the constant ofintegration is neglected. With this, only the particular solution is obtained from Equation B.10. – D f ( t ) , Z f ( t ) dt ; with 1 D c (B.8)We redefine Equation B.7 by mapping D onto an indefinite integral and retain constants. Usingthis modification, the full solution is obtained from Equation B.10 rather than just the particularsolution . The solution to the homogeneous problem can also be obtained by x ( t ) = P ( D ) f ( t ) = 0 . • D + a f ( t ) = e − at D e at f ( t ) (B.9)which, in turn, is solvable using Definition B.7 / B.8. • x ( t ) = 1 P ( D ) f ( t ) (B.10)is used to solve for the solution of the uni-ccLODE P ( D ) x ( t ) = f ( t ) as follows: – P ( D ) f ( t ) = n X i =1 γ i D + α i f ( t ) (B.11)Partial fraction expansion then invoking Equation B.9 to obtain the solution. This expansion isappropriate unless there are repeated roots α i . – P ( D ) f ( t ) = n Y i =1 D + α i f ( t ) (B.12)Factorization and sequentially invoking Equation B.9 to obtain the solution. This is alwaysappropriate. • P ( D ) 1 Q ( D ) u ( t ) = 1 Q ( D ) P ( D ) u ( t ); P ( D ) Q ( D ) u ( t ) = Q ( D ) P ( D ) u ( t ); 1 P ( D ) 1 Q ( D ) u ( t ) = 1 Q ( D ) 1 P ( D ) u ( t )(B.13)where u ( t ) is a smooth function and P ( D ) , Q ( D ) are polynomials in the differential operator D .When applied to smooth functions the order of operations does not matter. This is the case due tothe constant nature of polynomial coefficients.Note that the factorization of P ( D ) is over complex numbers and hence the irreducible factors are alwaysfirst order polynomials. To summarize, Equations B.11 / B.12 are used in conjunction with Equation B.9and Definition B.7 / B.8 to solve uni-ccLODEs P ( D ) x ( t ) = f ( t ). Classically, the complementary andparticular solutions are found separately: the complementary solution is obtained from the characteristic Consider, for instance, the solution to the concocted uni-ccLODE (cid:0) D − D + 2 (cid:1) x ( t ) = e it . Applying Equation B.12or B.11 using the classical definition (Equation B.7) we obtain only a particular solution, x ( t ) = − i e it and must obtainthe complementary solution separately. Whereas using the modified definition (Equation B.8), we obtain the full solution x ( t ) = c e t + c e t + − i e it For example, solving ( D − a ) x ( t ) = f ( t ), with Equation B.12 and the modified definition, very directly yields a comple-mentary solution ( c + c t ) e at P ( a ) = 0, and the particular solution is obtained x ( t ) = P ( D ) f ( t ) by applying B.11 / B.12 anddefinition B.7. In contrast, using the modified definition (Equation B.8), the full solution is obtained from x ( t ) = P ( D ) f ( t ) using either Equation B.11 / B.12.In what follows, we demonstrate that the modified Definition B.8 allows us to arrive at the full solution.We also show that Equation B.12 along with Definition B.8 allows us to handle repeated root cases directly.In the main text of this paper, we use Definition B.8 (rather than Definition B.7), and prefer Equation B.12to Equation B.11 for deriving solutions. Appendix B.1. Full Solution Using Modified Definition
Using classical methods of differential operator theory (with Definition B.7), the complementary andparticular solutions are solved for separately and then summed to give the full solution. Using our simplemodification (Definition B.8), we may directly obtain the full solution. We do so using Equation B.11 orB.12.We illustrate this by solving the following uni-ccLODE example in detail so that the reader can easilyfollow and apply our methods for objective III. (cid:18) D + ( b − a ) D − ab (cid:19) x ( t ) = f ( t ) (B.14)Solving using Equation B.11, we get, x ( t ) = 1( D − a )( D + b ) f ( t ) (B.15)= 1 a + b (cid:18) D − a − D + b (cid:19) f ( t ) (B.16)= 1 a + b (cid:18) e at D e − at − e − bt D e bt (cid:19) f ( t ) (B.17)= 1 a + b (cid:18) e at Z e − at f ( t ) dt − e − bt Z e bt f ( t ) dt (cid:19) (B.18)For instance, for the case of f ( t ) = f o , this gives: x ( t ) = f o a + b (cid:18) e at Z e − at dt − e − bt Z e bt dt (cid:19) (B.19)= f o a + b e at (cid:18) e − at − a + c (cid:19) − e − bt (cid:18) e bt b + c (cid:19)! (B.20)= C e at + C e − bt − f o ab (B.21)which is the full solution. Appendix B.2. Handling Repeated Factors
In some cases, a uni-ccLODE P ( D ) x ( t ) = f ( t ) has a P ( D ) that, when factorized into first degreepolynomial factors, has repeated factors - this corresponds to repeated roots in the characteristic equationof the associated homogeneous equation.In this situation, the complementary solution (which we can obtain as part of the full solution due toDefinition B.8) has a special form. In addition, the uni-ccLODE must be solved using Equation B.12 andcannot be solved using B.11 (as is obvious by division by zero that is due to the coefficients of partial fractionexpansion). 21ue to repeated factors, the method of Equation B.12 applies more generally to solving uni-ccLODEsthan the method of Equation B.11 which cannot handle repeated first order factors. Therefore EquationB.12 is our method of choice throughout.To give show that Definition B.8 and Equation B.12 handle repeated factors, consider the followingexample, (cid:18) D + 2 D + D (cid:19) x ( t ) = f ( t ) (B.22)Again, we solve the uni-ccLODE in detailed steps so that the reader can easily follow and apply ourmethods for objective III. We get, x ( t ) = 1 D ( D + 1)( D + 1) f ( t ) (B.23)= 1 D D + 1 1 D + 1 f ( t ) (B.24)= 1 D D + 1 e − t D e t f ( t ) (B.25)= 1 D e − t D D e t f ( t ) (B.26)= Z e − t Z Z e t f ( t ) dtdtdt (B.27)For instance, if f ( t ) = t , this gives, x ( t ) = Z e − t Z Z te t dtdtdt (B.28)= Z e − t Z ( te t − e t + c ) dtdt (B.29)= Z e − t ( c t + te t − e t + c ) dt (B.30)= C e − t + C te − t + C + t − t (B.31)where the t in the second term of the complementary solution is due to repeated factors.22 eferences [1] A. Schutte, F. Udwadia, New approach to the modeling of complex multibody dynamical systems, Journal of AppliedMechanics 78 (2) (2011).[2] C. M. Pappalardo, D. Guida, On the computational methods for solving the differential-algebraic equations of motion ofmultibody systems, Machines 6 (2) (2018) 20.[3] A. M. Karam, T. M. Laleg-Kirati, Electrical equivalent thermal network for direct contact membrane distillation modelingand analysis, Journal of Process Control 47 (2016) 87–97.[4] A. M. Karam, T. M. Laleg-Kirati, Membrane fouling modeling and detection in direct contact membrane distillation,Journal of Process Control 81 (2019) 190–196.[5] L. Biegler, J. Damiano, G. Blau, Nonlinear parameter estimation: a case study comparison, AIChE Journal 32 (1) (1986)29–45.[6] F. Albalawi, Lyapunov-based economic model predictive control for nonlinear descriptor systems, Chemical EngineeringResearch and Design 163 (2020) 263–272.[7] C.-W. Ho, A. Ruehli, P. Brennan, The modified nodal approach to network analysis, IEEE Transactions on circuits andsystems 22 (6) (1975) 504–509.[8] Q. Chen, S.-H. Weng, C.-K. Cheng, A practical regularization technique for modified nodal analysis in large-scale time-domain circuit simulation, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 31 (7) (2012)1031–1040.[9] R. M¨arz, C. Tischendorf, Recent results in solving index-2 differential-algebraic equations in circuit simulation, SIAMJournal on Scientific Computing 18 (1) (1997) 139–159.[10] M. G¨unther, U. Feldmann, The dae-index in electric circuit simulation, Mathematics and Computers in Simulation 39 (5-6)(1995) 573–582.[11] F. Vicario, S. Alkhairy, R. Buizza, W. A. Truschel, Two-parameter leak estimation in non-invasive ventilation, in: 201739th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE, 2017,pp. 153–156.[12] A. Hubbard, A traveling-wave amplifier model of the cochlea, Science 259 (5091) (1993) 68–71.[13] E. L. Ginzton, W. R. Hewlett, J. H. Jasberg, J. D. Noe, Distributed amplification, Proceedings of the IRE 36 (8) (1948)956–969.[14] G. Zweig, R. Lipes, J. Pierce, The cochlear compromise, The Journal of the Acoustical Society of America 59 (4) (1976)975–982.[15] K. N. Stevens, Acoustic phonetics, Vol. 30, MIT press, 2000.[16] A. Alkhairy, Mathematical models of vocal tract with distributed sources, in: 2003 IEEE International Conference onAcoustics, Speech, and Signal Processing, 2003. Proceedings.(ICASSP’03)., Vol. 1, IEEE, 2003, pp. I–I.[17] S. J. Elliott, G. Ni, An elemental approach to modelling the mechanics of the cochlea, Hearing research 360 (2018) 14–24.[18] E. Kaltofen, J. P. May, Z. Yang, L. Zhi, Approximate factorization of multivariate polynomials using singular valuedecomposition, Journal of Symbolic Computation 43 (5) (2008) 359–376.[19] S. Gao, Factoring multivariate polynomials via partial differential equations, Mathematics of computation 72 (242) (2003)801–822.[20] C. Bajaj, J. Canny, T. Garrity, J. Warren, Factoring rational polynomials over the complex numbers, SIAM Journal onComputing 22 (2) (1993) 318–331.[21] T. Sasaki, M. Suzuki, M. Kol´ar, M. Sasaki, Approximate factorization of multivariate polynomials and absolute irre-ducibility testing, Japan journal of industrial and applied mathematics 8 (3) (1991) 357.[22] W. Magnus, On the exponential solution of differential equations for a linear operator, Communications on pure andapplied mathematics 7 (4) (1954) 649–673.[23] G. Zheng, D. Boutat, H. Wang, A nonlinear luenberger-like observer for nonlinear singular systems, Automatica 86 (2017)11–17.[24] F. Milano, I. Dassios, Primal and dual generalized eigenvalue problems for power systems small-signal stability analysis,IEEE Transactions on Power Systems 32 (6) (2017) 4626–4635.[25] E. Yip, R. Sincovec, Solvability, controllability, and observability of continuous descriptor systems, IEEE transactions onAutomatic Control 26 (3) (1981) 702–707.[26] L. Petzold, Differential/algebraic equations are not ode’s, SIAM Journal on Scientific and Statistical Computing 3 (3)(1982) 367–384.[27] J. Rommes, N. Martins, F. D. Freitas, Computing rightmost eigenvalues for small-signal stability assessment of large-scalepower systems, IEEE transactions on power systems 25 (2) (2009) 929–938.[28] A. Ilchmann, L. Leben, J. Witschel, K. Worthmann, Optimal control of differential-algebraic equations from an ordinarydifferential equation perspective, Optimal Control Applications and Methods 40 (2) (2019) 351–366.[29] T. Reis, M. Voigt, Linear-quadratic optimal control of differential-algebraic systems: the infinite time horizon problemwith zero terminal state, SIAM Journal on Control and Optimization 57 (3) (2019) 1567–1596.[30] H. Cheng, Advanced analytic methods in applied mathematics, science, and engineering, LuBan Press, 2007.[1] A. Schutte, F. Udwadia, New approach to the modeling of complex multibody dynamical systems, Journal of AppliedMechanics 78 (2) (2011).[2] C. M. Pappalardo, D. Guida, On the computational methods for solving the differential-algebraic equations of motion ofmultibody systems, Machines 6 (2) (2018) 20.[3] A. M. Karam, T. M. Laleg-Kirati, Electrical equivalent thermal network for direct contact membrane distillation modelingand analysis, Journal of Process Control 47 (2016) 87–97.[4] A. M. Karam, T. M. Laleg-Kirati, Membrane fouling modeling and detection in direct contact membrane distillation,Journal of Process Control 81 (2019) 190–196.[5] L. Biegler, J. Damiano, G. Blau, Nonlinear parameter estimation: a case study comparison, AIChE Journal 32 (1) (1986)29–45.[6] F. Albalawi, Lyapunov-based economic model predictive control for nonlinear descriptor systems, Chemical EngineeringResearch and Design 163 (2020) 263–272.[7] C.-W. Ho, A. Ruehli, P. Brennan, The modified nodal approach to network analysis, IEEE Transactions on circuits andsystems 22 (6) (1975) 504–509.[8] Q. Chen, S.-H. Weng, C.-K. Cheng, A practical regularization technique for modified nodal analysis in large-scale time-domain circuit simulation, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 31 (7) (2012)1031–1040.[9] R. M¨arz, C. Tischendorf, Recent results in solving index-2 differential-algebraic equations in circuit simulation, SIAMJournal on Scientific Computing 18 (1) (1997) 139–159.[10] M. G¨unther, U. Feldmann, The dae-index in electric circuit simulation, Mathematics and Computers in Simulation 39 (5-6)(1995) 573–582.[11] F. Vicario, S. Alkhairy, R. Buizza, W. A. Truschel, Two-parameter leak estimation in non-invasive ventilation, in: 201739th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE, 2017,pp. 153–156.[12] A. Hubbard, A traveling-wave amplifier model of the cochlea, Science 259 (5091) (1993) 68–71.[13] E. L. Ginzton, W. R. Hewlett, J. H. Jasberg, J. D. Noe, Distributed amplification, Proceedings of the IRE 36 (8) (1948)956–969.[14] G. Zweig, R. Lipes, J. Pierce, The cochlear compromise, The Journal of the Acoustical Society of America 59 (4) (1976)975–982.[15] K. N. Stevens, Acoustic phonetics, Vol. 30, MIT press, 2000.[16] A. Alkhairy, Mathematical models of vocal tract with distributed sources, in: 2003 IEEE International Conference onAcoustics, Speech, and Signal Processing, 2003. Proceedings.(ICASSP’03)., Vol. 1, IEEE, 2003, pp. I–I.[17] S. J. Elliott, G. Ni, An elemental approach to modelling the mechanics of the cochlea, Hearing research 360 (2018) 14–24.[18] E. Kaltofen, J. P. May, Z. Yang, L. Zhi, Approximate factorization of multivariate polynomials using singular valuedecomposition, Journal of Symbolic Computation 43 (5) (2008) 359–376.[19] S. Gao, Factoring multivariate polynomials via partial differential equations, Mathematics of computation 72 (242) (2003)801–822.[20] C. Bajaj, J. Canny, T. Garrity, J. Warren, Factoring rational polynomials over the complex numbers, SIAM Journal onComputing 22 (2) (1993) 318–331.[21] T. Sasaki, M. Suzuki, M. Kol´ar, M. Sasaki, Approximate factorization of multivariate polynomials and absolute irre-ducibility testing, Japan journal of industrial and applied mathematics 8 (3) (1991) 357.[22] W. Magnus, On the exponential solution of differential equations for a linear operator, Communications on pure andapplied mathematics 7 (4) (1954) 649–673.[23] G. Zheng, D. Boutat, H. Wang, A nonlinear luenberger-like observer for nonlinear singular systems, Automatica 86 (2017)11–17.[24] F. Milano, I. Dassios, Primal and dual generalized eigenvalue problems for power systems small-signal stability analysis,IEEE Transactions on Power Systems 32 (6) (2017) 4626–4635.[25] E. Yip, R. Sincovec, Solvability, controllability, and observability of continuous descriptor systems, IEEE transactions onAutomatic Control 26 (3) (1981) 702–707.[26] L. Petzold, Differential/algebraic equations are not ode’s, SIAM Journal on Scientific and Statistical Computing 3 (3)(1982) 367–384.[27] J. Rommes, N. Martins, F. D. Freitas, Computing rightmost eigenvalues for small-signal stability assessment of large-scalepower systems, IEEE transactions on power systems 25 (2) (2009) 929–938.[28] A. Ilchmann, L. Leben, J. Witschel, K. Worthmann, Optimal control of differential-algebraic equations from an ordinarydifferential equation perspective, Optimal Control Applications and Methods 40 (2) (2019) 351–366.[29] T. Reis, M. Voigt, Linear-quadratic optimal control of differential-algebraic systems: the infinite time horizon problemwith zero terminal state, SIAM Journal on Control and Optimization 57 (3) (2019) 1567–1596.[30] H. Cheng, Advanced analytic methods in applied mathematics, science, and engineering, LuBan Press, 2007.