The Parametric Solution of Underdetermined linear ODEs
aa r X i v : . [ c s . S C ] A ug The Parametric Solution of Underdetermined linearODEs
Thomas WolfDepartment of Mathematics, Brock University500 Glenridge Avenue, St.Catharines,Ontario, Canada L2S 3A1email: [email protected] 9, 2018
Abstract
The purpose of this paper is twofold. An immediate practical use of the presentedalgorithm is its applicability to the parametric solution of underdetermined linearordinary differential equations (ODEs) with coefficients that are arbitrary analyticfunctions in the independent variable. A second conceptual aim is to present analgorithm that is in some sense dual to the fundamental Euclids algorithm, and thusan alternative to the special case of a Gr¨obner basis algorithm as it is used for solvinglinear ODE-systems. In the paper Euclids algorithm and the new ‘dual version’ arecompared and their complementary strengths are analysed on the task of solvingunderdetermined ODEs. An implementation of the described algorithm is interactivelyaccessible under [7].
Underdetermined ordinary and partial differential equations (ODEs/PDEs) are typical ob-jects of investigation in control theory but also overdetermined systems of equations whereoriginally the number of equations is larger than the number of functions may turn intounderdetermined problems in the course of their solution. Examples are the conditions forLie-symmetries and conservation laws of linearizable differential equations. These overdeter-mined systems of conditions have solutions that involve arbitrary functions of one or morevariables, i.e. in the process of solving these systems towards its conclusion, underdeterminedproblems often occur.Algorithms for solving underdetermined linear equations and systems (ODEs, PDEs,multidimensional discrete,..) are known ([1],[2],[3]). There also exists efficient softwareby Daniel Robertz et al. [4], [5]. The purpose of this paper is to describe an alternativealgorithm that is elegant, efficient and in some sense complementary to the fundamentalEuclids algorithm which is the basis for the Gr¨obner basis algorithms that are used in otherimplementations. An open question is whether the new algorithm can be generalized topartial differential equations (PDEs) and thus not only be a complimentary algorithm forthe non-commutative differential Gr¨obner basis algorithm for solving ODEs but also forsolving PDEs.An implementation of algorithms described in this paper is accessible online under [7]and is part of the package
Crack ([8]) for solving overdetermined systems. As we will1escribe, a new feature of these programs are flags which allow to prevent denominators orthat allow to reduce the number of terms in the solution. This can be achieved by absorbingexplicit x -dependent factors into new functions that are introduced during the computation.We consider a single linear ODE0 = r X i D i f i + c ( x ) , r > , (1) D i = n i X j =0 c ij ( x ) ddx ! j , n i > f . . . f r of the independent variable x . The equation may be homogeneous( c = 0) or inhomogeneous ( c = 0). With the requirement that the ODE is underdeter-mined we only assume that at least two functions f i are involved, i.e. r >
1. The differentialorder of each f i is n i . For simplicity we assume that all coefficients c ij are sufficiently oftendifferentiable, at least P i n i times. An upper index in round brackets indicates the numberof differentiations, i.e. d n f i dx n = f ( n ) i and low order derivatives are denoted by apostrophe, e.g. d f i dx = f ′′ i .The task is to find the general solution of the ODE (1) in the form of explicit differentialexpressions F i f i = F i ( x, h ( x ) , . . . , h r − ( x )) , i = 1 , . . . , r (2)in terms of parametric functions h ( x ) , . . . , h r − ( x ) which are either all free or one of them,say h ( x ), having to satisfy an ODE and the other h j ( x ) being free.For example, the general solution of the ODE0 = f ′′ x + g ′′ x − g ′ x + f + 3 x for f = f ( x ) , g = g ( x ) can be written in the form f = xx − x + 7 x − x + 9 (cid:16) ( x − x + 3 x ) h ′′ − ( x + x + 3 x − h ′ (3)+(3 x + 3 x + 17 x ) h − x + 3 x − x + 9 x (cid:17) ,g = x x − x + 7 x − x + 9) (cid:16) ( − x + 2 x − x ) h ′′ (4)+(8 x − x ) h ′ − (14 x + 14 x + 6) h + 4 x + x + 3 x − x (cid:17) (5)where h = h ( x ) is an arbitrary function. The form of the solution is not unique. Because h ( x ) is arbitrary, replacing h by a differential expression in one or more arbitrary functionswould give a solution too, but a solution with expressions of higher differential order inmore arbitrary functions without the solution being more general. One naturally seeks ageneral solution which involves only as few as possible arbitrary functions and lowest orderderivatives of them. But even this requirement does not give a unique solution. One mightwant to minimize the highest order of all parametric functions or, for example, the sum ofall orders of all parametric functions.Furthermore, the occurring function(s) could be scaled to modify the form of the solution,for example, to make it denominator free. In the above case replacing h ( x ) by ( x − x +7 x − x + 9) p ( x ) with an arbitrary function p ( x ) would make the solution polynomial butincrease its size, i.e. the total number of terms in the coefficients of the parametric functionson the right hand sides.In the following section the algorithm is given, followed by comments on its character-istics. In section 3 and 4 we formulate the new algorithm and Euclids algorithm in vectornotation in order to compare them in detail in section (5). We show that both are essentiallyifferent but also that they compliment each other in the sense that one can give a criterionunder which circumstances which of both is better suited. Another option is to combineboth algorithms in a hybrid version.Finally, in section (7) an application is described which arose from a classification ofhyperbolic evolutionary vector PDEs. The essence of the algorithm is • to partition the ODE into a total derivative and an algebraic remainder, • to introduce a new function f r +1 such that the total derivative part is df r +1 /dx andthus to write the ODE as a system of 2 equations: one equation defining f r +1 and onere-formulating the ODE in terms of f r +1 , • to use the re-formulated ODE to eliminate and substitute another function and thusto arrive again at a single ODE for the same number of functions which in some senseis closer to be solvable algebraically than the ODE before.These steps are repeated until either the ODE involves only one function and thus is notunderdetermined anymore, or until one function occurs purely algebraically and thus allow-ing the ODE to be solved for that function. The following is a more detailed pseudo codedescription. Input • list of functions f , . . . , f r of x , • linear ODE 0 = ω in f i (like (1)) Body L := {} % L will be a list of substitutions s := r while (the ODE 0 = ω involves at least two f i ) and (differential orders n i > ∀ f i ) do • Factor out ddx from ω once as far as possibleby introducing a new function f s +1 ( x ) andby computing expressions b i :0 = ω = f ′ s +1 + s X i =1 f i b i + a (6)where f s +1 , b i are given through: f s +1 = s X i =1 n i − X j =0 f ( n i − − j ) i j X k =0 ( − ( j + k ) a ( j − k ) i,n i − k (7) b i = n i X k =0 ( − ( n i + k ) a ( n i − k ) i,n i − k (8) if b i = 0 for all i then the ODE 0 = ω is exact (apart from a ), i.e.consider new ODE 0 = ˆ ω := f s +1 + R a dx where f s +1 is only an abbreviation defined in (7) else ◦ regard f s +1 as a new unknown function, ◦ solve (6) for one function f j of f , .., f s that − has a non-vanishing coefficient b i in (6), and − is of lowest possible order in ω and thus in (7), f j = − b j f ′ s +1 + X i = j f i b i + a (9) ◦ use (9) to substitute f j in (7) to get a new ODE0 = ˆ ω := − f s +1 + P si =1 P n i − j =0 . . . for functions f , .., f j − , f j +1 , .., f s , f s +1 ◦ update ω := ˆ ω, s := s + 1 , L := { f j = .. } ∪ L endif the ODE involves a function f j purely algebraically then solve for f j and add it to L : L := { f j = .. } ∪ Lh ( x ) := f s ( x ) % to remember the last introduced function P := L % to remember the complete list of substitutions L % next substitutions each as stored in the first% element of P are performed in the rest of P while s > r do P := rest ( P ) | f s = ... % as given in first ( P ) s := s − endend % of Body Output • l % list of new functions • L % the complete list of substitutions • P % All initial functions f , .., f r are either% parametric, or are given in P in terms of h • ω ( x, h ( x )) % only if the first while loop% terminates due to ω not% being underdetermined anymore Before comparing the algorithm with Euclids algorithm a few comments are necessary. • All steps involve only algebra or differentiations with only one exception: if the ho-mogeneous part of the ODE is exact then the integral of the inhomogeneous part istaken. But this integral does not have to be evaluated, i.e. to be expressed in termsof elementary functions. It can stay in a symbolic unevaluated form for the algorithmto continue. • The coefficients a ij can be arbitrary explicit functions of x , and involve, for example,sin or log with the only condition that it must be decidable whether expressionsinvolving these functions and their derivatives are zero or not. The transformation of the ODE in each step of the first while loop is reversible, i.e.the new ODE is equivalent to the previous one. Therefore the obtained solution is thegeneral one. • The algorithm terminates.Proof:We consider how the total sum P si =1 n i of differential orders n i of all functions f i changes during execution. After substitution of f j with (9) in (7) the new differentialorders ˆ n i of functions f i that occur in the new ODE are:ˆ n i = n j for i = s + 1 ≤ max( n j − , n i −
1) = n i − i ≤ s, b i = 0= n i − i ≤ s, b i = 0 . The function f j of order n j gets replaced by a new function f s +1 which then occurswith same order n j , but the order of all other functions is lowered by at least onebecause we choose j such that n j = min( n k ) ( ∀ k with b k = 0), i.e. we have ˆ n i ≤ n i − f i with b i = 0. Because the ODE has at least two functions, the total sumof derivatives is decreasing. The algorithm is therefore finite. • The algorithm results in an ODE for a single function iff the differential operators D i in (1) have a common differential factor. The remaining ODE is of the same order asthe common factor, i.e. its order can not be higher than the order of the original ODE(1). • The algorithm naturally splits into two parts, a part A) establishing a list L of sub-stitutions (9) in the first while loop and part B) performing the substitutions in thesecond while loop to obtain an explicit solution. The first part is executed very fast(see next section for more details), the second may take longer for higher order ODEsbecause expressions typically grow with each substitution, often exponentially. Anexample is given in the appendix.For many applications the list L of substitutions may even be of higher practical valuethan the explicit formulas f i = F i ( x, h ( x ) , . . . , h r − ( x )) , i = 1 , . . . , r resulting fromB). In general, the list L is a much shorter representation of the parametric solutionof the ODE than the explicit solution itself and thus is more useful as a solution, likein the following scenario.Let us assume that we have a large algebraic expression in terms of the originalfunctions f , .., f r that has to be evaluated modulo the underdetermined ODE (1).Instead of replacing the f i directly by their large explicit expressions as given in thelist P it is usually much better to perform successively the substitutions stored in L in the order they were derived, allowing cancellations to happen after each individualsubstitution. Also, if numerical computations are to be done, it is much faster tocompute the sequence of substitutions than the explicit expressions. To perform the iteration process efficiently, i.e. to replace functions by algebraic combina-tions of other functions (9) in the ODE (7) without needing to perform any differentiations,the ODE has to be represented in a form where D := ddx is factored out as far as possible. Aswe will see further down, this is also the appropriate representation for using the Euclideanlgorithm to solve the underdetermined ODE. We write the ODE in the form0 = r X i A i f i + a ( x ) (10) A i = D ˜ A i + a i ( x ) (11)˜ A i = D n i − a in i ( x ) + . . . + Da i ( x ) + a i ( x ) (12)Replacing, for example, f = w ( x ) ˆ f would only require multiplications ˆ a k := wa k toupdate this representation and no differentiations.In this notation a single step in the first while loop of the new method consists of • introducing a new function f r +1 through f r +1 = X i ˜ A i f i (13)giving the ODE the form 0 = Df r +1 + X i a i ( x ) f i + a . (14) • regarding the defining relation (13) for the new function as the new ODE and usingthe old ODE (14) for a substitution of an ‘old’ function. For that we choose one a j ofthe non-vanishing a i in (11) for which the order n j of the corresponding f j is minimal.Thus, performing the substitution f j = − a j Df r +1 + X i = j a i f i + a (15)in the definition (13) gives the new ODE0 = − f r +1 + X i = j ˜ A i f i − X i = j ˜ A j a i a j f i ! − ˜ A j a j Df r +1 ! − ˜ A j a a j ! . (16)In vector notation the above two steps change the differential operators (assuming for sim-plicity of notation j = 1) according to: A A : A r → A r +1 A : A r = − − ˜ A a D ˜ A − ˜ A a a :˜ A r − ˜ A a r a (17) a → − ˜ A (cid:18) a a (cid:19) (18)and add (15) to the list of substitutions. Updating components 2 . . . r in (17) is done purelyby multiplications and additions, only updating the first component takes differentiationswhen factoring out D .A more conventional algorithm solving underdetermined linear ODEs tries to lower theorder of the ODE by performing Euclids algorithm until the ODE is algebraic for onefunction and thus can be solved algebraically. In the following section we formulate such anEuclidean step to be able to compare it with the new algorithm above. The Euclid version in vector notation
For applying the ‘right’ Euclid algorithm the ODE is given as well in a form with D com-pletely factored out: 0 = r X i A i f i + a ( x ) (19) A i = D n i a in i ( x ) + . . . + Da i ( x ) + a i ( x ) . (20)One iteration step is performed by • choosing two functions f i , f j (w.l.o.g. n i ≥ n j ) and introducing a new function f r +1 through f j = f r +1 − D n i − n j a in i a jn j f i ! (21) • and performing this substitution in the ODE (19).The substitution is chosen such that the differential order of f i is lowered by at least one.This iteration process stops like in the new algorithm when the first function appears purelyalgebraically or when the ODE involves only one function. Therefore, in order to minimizethe number of steps one would choose f j as one of the lowest order functions and f i asone of the lowest order functions from the remaining ones. (From all these choices it isrecommendable to choose a pair so that multiplication is minimized and A i , A j have as fewas possible terms.) For j = 1 , i = 2 the vector notation gives A A A : A r → A r +1 A A : A r = A A − A D n − n a n a n A : A r (22) a → a . (23)The update of the 2 nd component in (22) requires differentiations when D is factored outas far as possible. In this section we want to justify our claim that both algorithms differ significantly fromeach other and that they are not merely variations of one and the same procedure.Both algorithms differ conceptually in that Euclids algorithm pairs two differential op-erators to get a new differential operator of lower order whereas in the new algorithm therelation that defines a new function becomes the ODE and the old ODE is used for substi-tuting one of the ‘old’ functions.Another difference is that Euclids algorithm lowers in one step the order of only onedifferential operator whereas the new algorithm lowers the order of all but one differentialoperators. Both types of steps are performed at about the same cost as for the new methodonly the update of the first component in (17) is potentially size increasing and for Euclidsmethod it is only the update of the second component of (22). Using Euclids algorithm onean of course take one of the lowest order differential operators and decrease the order ofall others but that would take r − A i with each other, all other higher order operators keeptheir high order. Moreover, when the algorithm stops because the final equation containsone function algebraically, then back substitutions start (part B of the algorithm) whicheven increase the highest orders.Differently with the new method where the differential order of all functions but onegets lowered in one step. Thus the obtained parametric solution is typically of lower orderthan the solution obtained by Euclids algorithm if this is executed by pairing only lowestorder derivatives.For example, the size of the parametric solution of the ODE0 = x f ax + ( x − g bx + h x , for f, g, h of x (24)and the differential order of parametric functions in it depend strongly on the differentialorders a, b of f and g in the ODE and the method that is used. As a, b increase from a = b = 1 to a = b = 3 one can see the following trends: • In Euclids solution the differential order of parametric functions increases from 5 to7 whereas in the solution of the new algorithm the order decreases from 4 to 3 in theexpression for g , from 4 to 2 in the expression for f and only increases from 0 to 1 inthe expression for h . • The size of expressions in Euclids solution steadily increases and in the solution of thenew method it decreases.For a = b = 1 the solutions have comparable size:- When explicit x -dependent factors are not absorbed (for absorbing factors see sec-tion 5.3) then h is parametric and the rational expressions for f, g have the form (5terms)/(2 terms) in Euclids solution and in the solution of the new method f = (4terms)/(2 terms), g = (3 terms)/(2 terms).For a = b = 3:- When factors are not absorbed then in Euclids solution h is parametric, f = (85terms)/(15 terms), g = 13 terms and in the solution of the new algorithm we have f = 1 term, g = 7 terms, h = 4 terms.- When factors are absorbed then Euclids solution gives f = 27 terms, g = 16 termsand in the solution of the new algorithms we have f = 1 term, g = 7 terms, h = 4terms.Not only the solutions of both algorithms may differ significantly but also the numberof steps to reach them. Starting, for example, with0 = f ′ + a ( x ) h ( n ) = 0to be solved for f = f ( x ) , h = h ( x ) this takes only one iteration step with the new algorithmgiving, for example for n = 5, a solution of the form f = (27 terms)/(1 term), h = (1term)/(1 term), whereas it takes n steps for Euclids algorithm giving, for example for n = 5,a more spacious solution of the form f = (40 terms)/(1 term), h = (2 terms)/(1 term).f the ODE is inhomogeneous then another difference between both algorithms becomesapparent. With non-zero inhomogeneity a in the ODE (10) the transformation (15) be-comes inhomogeneous and changes a in (18) whereas Euclids algorithm does not change a in the homogeneous transformation (21).Two ODEs that highlight the duality between both algorithms are0 = f ′ + f + g ′ + ( ah ) (20) (25)0 = f ′ + f + g ′ + a ( h ) (20) (26)where a = a ( x ) is a given function of x and the ODE has to be solved for f, g, h . As shownin table 1 for the ODE (25) the new algorithm is to be preferred whereas for ODE (26)Euclids algorithm gives a shorter solution.algorithm equation (25) equation (26)new f =1 term g =22 terms f =2 terms g =23 termsEuclid f =22 terms g =23 terms f =2 terms g =3 termsTable 1. The size of solutions given by both algorithms (both executed without absorbingof factors).An explanation of the above behaviour comes from the fact that Euclids method is bettersuited when factors a ini a jnj are small, i.e. when the coefficients of the leading derivatives in therepresentation (10) - (12) effectively cancel each other in a quotient, and less suited whenthese factors are large prime expressions in x . Differently, the new method works best iffactors a i a j are small, i.e. when the coefficients of the algebraic terms effectively cancel eachother, and less suited when these factors are large prime expressions in x .Although both ODEs look similar, they are rather different as ( ah ) (20) has many termsif the product rule of differentiation is fully applied and a ( h ) (20) has many terms if ddx isfactored out as far as possible. Thus both algorithms differ in their suitability for bothODEs.The above computations and any other tests can easily be performed online (see [7]),where access to the computer algebra system REDUCE and to the procedures uode and print uode solution is provided. The advantage of having these two, in some sense, complementary algorithms lies in thefact that both operate on the same data structure, i.e. both their input and output consistsof an ODE in the form (10) - (12) and both add one substitution to a list of substitutionsof functions in terms of newly introduced functions. Consequently, both algorithms areinterchangeable, i.e. one step performed with one algorithm could be followed by anotherstep performed with the other algorithm, in order to minimize the size of x -dependent factorsand thus to minimize growth. A different strategy could be to perform each step with bothalgorithms in parallel, to compare the size of the ODE and/or size of derived substitutionthat both methods give and to choose which one to adopt for this step. This strategy does atmost double the amount of computation but more likely will lead to savings due to workingwith smaller expressions.Any hybrid algorithm performing a mixture of Euclid steps and new algorithm steps isstill finite because in any such step the sum of differential orders of all functions is decreasing. The following efficiency improving measures work for both methods. A representation (10)- (12) (and identically (19) - (20)) of the ODE where D is maximally factored out has thedvantage that any change of functions f i → h ( x ) ¯ f i does require only multiplications and isdone very easily in both algorithms. This freedom of efficiently multiplying functions with x -dependent factors can be used for different purposes.If in a newly generated ODE all coefficients of a function f i have a non-trivial GCDˆ c ( x ) := GCD ( a in i , . . . , a i ) = 1 then this can be absorbed into a new function f r +1 and theODE be simplified by performing the substitution f i = f r +1 / ˆ c (27)and adding it to the accumulating list of substitutions. Such non-trivial GCDs occur rel-atively frequently as the GCD is taken only from the few coefficients a ij of any single oneoperator A i that changed in an iteration step.A different purpose of introducing new functions multiplied with an explicit x -dependingfactor is to avoid a denominator (den). This would arise in the first component of (17)and it can be prevented by introducing another new function f r +2 and by performing thesubstitution f r +1 = a f r +2 . (28)Similarly, in the other components of (17) one can scale f i = den (cid:18) a i a (cid:19) f r + i . (29)In Euclids algorithm the appearance of denominators in the second component of (22)can be prevented by introducing the scaling f = den a n a n ! f r +1 . (30)Another situation where a dominator occurs is at the end of both algorithms when onefunction occurs purely algebraically0 = a i f i + X j = i A j f j + a ( x ) . (31)To avoid the denominator a i one can scale all f j depending on their differential order andon their coefficients in A j . Avoiding this denominator is especially helpful, as f i = . . . isthe last substitution and thus substituted successively in all previous substitutions leadingeasily to a substantial growth of denominators. By avoiding denominators in the above wayboth methods produce denominator free solutions if the ODE is homogeneous. The case that an underdetermined ODE factorizes, i.e. that the differential operators A i in(10) have a non-trivial GCD, or in other words, that the ODE can be written in form of twonested ODEs 0 = Ω( x, ω ( x, f i )) is discovered by both algorithms. The inner ODE 0 = ω isonly determined up to a linear change ω = α ( x )ˆ ω and both algorithms will usually find Ω’sthat differ by some α ( x ). A slight advantage of the new algorithm is that the D part ofthe ODE is used within the algorithm, thus it is automatically recognized if this vanishes,i.e. if the ODE is exact (up to an inhomogeneity). ODE-systems
The introduced methods of solving a single underdetermined ODE (or converting it to anODE for a single function) can be used to convert an ODE system into an equivalent set offully decoupled ODEs, each for a single function.After treating any one equation of the original ODE system the original functions in itare expressed in terms of fewer functions which are either all free or at most one has tosatisfy a single ODE. All functions can be replaced in the remaining ODEs. This can go onas long as we have equations of at least two functions. When the procedure stops we havesolved the system or got ODEs, each containing only one function. If more than one ODEcontain only one and the same single function, then these ODEs form an over-determinedsubsystem which can be treated by a Gr¨obner bases computation and result either in theexplicit solution for this function or a single ODE for this function of an order not higherthan the lowest order of the ODEs for this function. This whole procedure terminates witheither the explicit parametric solution of the original system or a decoupled set of ODEs,each ODE for a single function.
A class of applications where underdetermined linear ODEs occur frequently is the classi-fication of hyperbolic evolutionary PDE systems. The aim of such an investigation is tofind integrable systems of PDEs by determining those systems which have a higher ordersymmetry (see below). More information about the mathematical background is given in[6] where a classification of hyperbolic vector
PDEs is discussed.Let us look at two hyperbolic scalar PDEs for functions u ( x, t, τ ) , v ( x, t, τ ) (where x, t arethe usual independent variables and τ is a symmetry parameter). The following ansatz forthe system and symmetry is generated based on homogeneity considerations. By requiringthe same homogeneity weights as the potential nonlinear Schr¨odinger equation we get forthe system the ansatz u tx = a u x + a u x v x u + a u x u x v + a u x v x u + a u x v x + a u x v + a u x v x uv tx = a v x + a u x v x v + a v x v x u + a u x v x v + a u x v x + a v x u . (32)Assuming the same differential order for the symmetry we get the ansatz u τ = b u x + . . . , v τ = b v x + . . . , (33)with right hand sides identical to those of (32), only with coefficients b i instead of a i . Allcoefficients a i , b i are undetermined functions of the product uv . The relations (33) are considered to be a symmetry of the system (32) if the symmetryconditions ∂ τ ( u tx ) − ∂ t ∂ x ( u τ ) = 0 , ∂ τ ( v tx ) − ∂ t ∂ x ( v τ ) = 0 (34)are fulfilled identically in u, v and derivatives of u and v both modulo substitutions basedon (32) and (33). In performing the differentiations in (34), doing repeatedly substitutions(32), (33) and finally setting all coefficients of different products of powers of derivatives of u, v individually to zero gives 27 ODEs with a total of 1334 terms for 26 functions a i , b j of z := uv . The length of equations ranges from 2 to 266 terms. In the course of solving thissystem the program Crack performs integrations, substitutions, splittings (i.e. separations For hyperbolic systems the homogeneity weights include negative values, for example here weight( u )=1,weight( v )=-1, so that uv has weight zero and therefore we have no limitation on the degree of powers of uv and thus all unknown coefficients are arbitrary functions of uv . hen z occurs only explicitly) and a number of case distinctions. In one of the sub cases theresulting conditions can be integrated successively up to the underdetermined linear ODE0 = 3 b ′ z − b ′ z − b ′′ z + b ′ z − b z + 2 b . (35) Step 1:
To start, partial integration gives0 = (3 b z − b z − b ′ z + 5 b z ) ′ − b + 6 b z − b . (36)which can be written as 0 = c ′ − b + 6 b z − b (37)by introducing c ( z ) through c = 3 b z − b z − b ′ z + 5 b z. (38)From the functions that occur only algebraically in (37) (i.e. b , b , b ) the ones thathave lowest derivatives in (35) are b , b . Eliminating one of them from (37), say b andsubstituting it in (38) gives 0 = 2 b ′ z − b z − c ′ z + c (39)which is not algebraic in any function yet, but already of first order, so one more step hasto be performed. Step 2:
Partial integration of (39) and introduction of c = 2 b z − c z (40)results in 0 = c ′ − b z + 2 c (41)which allows to solve for c and replace it in (40) giving0 = 12 c ′ z − c − b z . (42)This condition is purely algebraic for one function, b , and therefore the algorithm stops. Cleanup:
We need explicit solutions of (35), so what remains to be done are back sub-stitutions: (42) provides b = 12 z c ′ − z c . (43)The second substitution expressing b in terms of the parametric function c is obtainedafter backward substituting c from (41) into (37) and b from (43) in (37) to get theexplicit solution consisting of (43) and b = 13 c ′′ − z c ′ + 2 b z + 2 z c (44)involving the free function c ( z ).Another underdetermined equation resulting in this integrability problem is0 = 3 b ′ z − b ′ z + 3 b ′ z − b ′′ z − b ′ z − b z + 6 b z − b (45)which has the solution b = ( − b ′ z + 3 c ′ z + 6 b z − b z − c ) / (3 z ) b = ( − b ′ z + c ′ z + 3 b z − b z − c ) / (3 z ) . We finally obtain as a system with higher order symmetries: u tx = a uv (cid:16) ( uv ) x v − uv ( uv ) x v x (cid:17) v tx = a uv ( uv ) x v x . Appendix
The following example illustrates the comments made in the last paragraphs of section 2.3.It shows two possible representations of a solution to an underdetermined ODE. For theequation ( x − f (5)1 + 3 f (3)1 + xf ′′ + (1 − x ) f ′ + f − ( x − x − f ′′ − xf ′ = 0the generated list L of substitutions is given through f = (7289 / f ′ x − / f ′ x + 725945 / f ′ x − / f ′ x + 2665016 / f ′ x − / f ′ x + 40582 f ′ − / f x + 393212 / f x − / f x +1404610 f x − / f x + 863648 / f ) / ( x − x x − x +12566 x − x + 111667 x − x + 332143 x − x + 625912 x − x + 117948) f = ( − / f ′ − / f x + 40655 / f x − / f x +937545 / f x − / f x + 236215 / f x − / f ) / ( x − / x + 30423 / x − / x + 46012 / x − / f = ( − / f ′ − / f ) / ( x − / x + 15683 / x − / x +18956 / x − / f = − / f ′ + 37 / f x − / f x + 347 / f x − / f x + 261 / f f = − / f ′ − / f x + 124 / f x − / f x + 11 f f = − / f ′ + 5 / f x − / f x − / f f = f ′ + 2 f x + f which is a much shorter representation of the solution than the 10 page explicit form whichresults from substituting f , f , f , f , f , f in this order into each other and takes the form f =(42 terms)/(25 terms), f =(304 terms)/(61 terms) involving up to 29-digit integers.The list of 7 substitutions is not only shorter than the list of 2 substitutions for f , f , italso is much faster to derive and more useful if a differential expression is to be simplifiedmodulo the solution of the above ODE by substituting f , f , f , f , f , f , f in this order.The difference in size of both solution representations can be arbitrarily amplified byhaving an input ODE of higher order and higher degree polynomials as coefficients. We present an algorithm for the solution of underdetermined linear ODE that is compatiblebut structurally different from the (right) Euclidean algorithm. Because both algorithmsoperate on the same data structure and because both complement each other in the sensethat each one is most efficient for ODEs of different form, the combination of both algorithmsis superior to each individual one.
Acknowledgements
I would like to thank Sergey Tsarev for many comments. Daniel Robertz is thanked forcomparative runs with the OreModules package. eferences [1] Chyzak, F., Quadrat, A., Robertz, D. (2005). ”Effective algorithms for parametrizinglinear control systems over Ore algebras”, Rapport de Recherche INRIA , ApplicableAlgebra in Engineering, Communications and Computing , no 5, 319-376.[2] Pommaret, J.-F., Quadrat, A. (2004). ”A differential operator approach to multidimen-sional optimal control”, International Journal of Control , 821-836.[3] Fr¨ohler, S., Oberst, U. (1998). ”Continuous time-varying linear systems”, Systems &Control Letters , 97-110.[4] Chyzak, F., Quadrat, A., Robertz, D. OreModules project [5] Chyzak, F., Quadrat, A., Robertz, D. ”OreModules: A symbolic package for the studyof multidimensional linear systems” in: J. Chiasson, J.-J. Loiseau, ”Applications ofTime-Delay Systems”, Springer, to appear.[6] Anco, S. and Wolf, T. (2005). ”Some symmetry classifications of hyperbolic vectorevolution equations”, JNMP, Volume 12, Supplement 1, p 13-31 (also nlin.SI/0412015).[7] Wolf, T. (2007). Online demo for solving underdetermined ODEs. http://lie.math.brocku.ca/crack/uode [8] Wolf, T.: Applications of Crack in the Classification of Integrable Systems, CRMProceedings and Lecture Notes, vol 37 (2004) pp. 283-300. (arXiv nlin.SI/0301032)and online underin the Classification of Integrable Systems, CRMProceedings and Lecture Notes, vol 37 (2004) pp. 283-300. (arXiv nlin.SI/0301032)and online under