LONG GUIDE TO THE STANDARD FORM PACKAGESpecial Purpose SystemsRELEASE 2: AUGUST 27, 1993 Gregory J. Reid and Allan D. Wittkopf REQUIREMENTS: Access to Maple: Version 5 (releases 1 or 2) or later version. The programs in the Standard Form Package which is available by anonymous ftp to: math.ubc.ca (i.e. after user type anonymous). Then change directory to pub/reid/standard and read the file README.1ST in that directory for instructions. ACKNOWLEDGEMENTS: We wish to especially thank: Cameron Connell who substantially contributed to the writing of this manual (especially to sections 3.6 and 6); David Weih (who provided the preversions of the coordinate transformation and interface routines); Alan Boulton who was an indispensable participant in many lively discussions about the package and who contributed the initial interface to the Liesymm package; and Ian Lisle whose beautiful moving frame algorithm is an inspiration for future developments. GJR gratefully acknowledges the generous support of a grant from the Natural Sciences and Engineering Research Council of Canada. We also wish to thank: Satish Reddy who performed the first computer algebra experiments back in 1987; Carl Wulfman for asking interesting questions which first turned GJR's mind to the direction of applying formal operations to differential equations to expose their properties; Sukeyuki Kumei for initial discussions; Beth Filliman for a test translation of the Macsyma prototype program into Maple; George Bluman for support and discussions, and providing systems of PDEs which were intractable enough that the development of such a formal package was kept on practical lines; Patrick Doran-Wu and Alex Ma for helpful suggestions; Rob Entwistle, Alistair Blachford, Tom Nicol and Bob Bajwa for expert computer support and Diana Kwan for her encouragement and support. CONTENTS 1. INTRODUCTION 2. GETTING STARTED WITH THE PACKAGE 2.0 HOW TO USE THIS MANUAL 2.1 SIMPLE EXAMPLES 2.2 SYSTEMS OF ALGEBRAIC EQUATIONS 2.3 NOW YOU'RE READY FOR SYSTEMS OF DIFFERENTIAL EQUATIONS 2.4 IMPORTANT NOTE ABOUT OPTIONS FOR standard_form 3. MORE ADVANCED ASPECTS OF THE PACKAGE 3.1 HOW TO CHANGE THE ORDERING - the option janet_order 3.2 CALCULATION OF INITIAL DATA AND TAYLOR SERIES SOLUTIONS OF SYSTEMS OF PDE 3.3 THE STRATEGY OPTION 3.3.1 DISPLAYING INTERMEDIATE INFO USING STRATEGY 3.3.2 CONTROLLING EQUATION FLOW USING STRATEGY'S KNOBS 3.3.3 THE STRATEGY OPTIONS INCONSOPT AND SIMPOPT 3.4 STORING INTERMEDIATE RESULTS AND THE USEREASY OPTION 3.4.1 STORING INTERMEDIATE RESULTS 3.4.2 SELECTING NICE EQUATION SUBSETS WITH USEREASY 3.5 CLASSIFICATION PROBLEMS - SYSTEMS OF PDE WITH UNSPECIFIED FUNCTIONS OR PARAMETERS 3.5.1 LOW-TECH APPROACH 3.5.2 HI-TECH APPROACH - the class_eqns option 3.6 INTERFACE TO MAPLE'S DERIVATIVE NOTATION 3.7 INTERFACE TO MANSFIELD'S NOTATION FOR HER DIFFGROB PACKAGE 3.8 THE DIRECTED GRAPH OF A PDE SYSTEM 4. APPLICATION OF THE PACKAGE TO SYMMETRY ANALYSIS OF PDES 4.1 CALCULATION OF THE STRUCTURE OF LIE SYMMETRY ALGEBRAS OF PDES FROM THEIR INFINITESIMAL DETERMINING SYSTEMS 4.1.1 SYMMETRIES OF THE HEAT EQUATION: Ut = Uxx 4.1.2 A GROUP CLASSIFICATION PROBLEM: Uxx = Utt + H(U) Ut 4.2 INTERFACE TO THE MACSYMA PACKAGE SYMGRP.MAX 4.2.1 PRODUCING THE DETERMINING SYSTEM OF Ut = Uxx 4.2.2 PRODUCING THE DETERMINING SYSTEM OF Uxx = Utt + H(U) Ut 4.3 INTERFACE TO THE MAPLE PACKAGE SYMMETRY 4.4 INTERFACE TO THE MAPLE PACKAGE LIESYMM 5. INTEGRATION OF THE EQUATIONS 5.1 USE OF taylor_exp 5.2 INTERACTIVE INTEGRATION 6. CHANGING VARIABLES IN SYSTEMS OF PDE 7. TOOLBOX OF INTERACTIVE FUNCTIONS 7.1 Clearderivdep 7.2 polyd 7.3 sfdiff 8. STRATEGIES TO DEAL WITH MEMORY EXPLOSION 8.1 ADJUSTING THE EQUATION-FLOW SETTINGS 8.2 VARYING THE ORDERING 8.3 CLASSIFICATION VARIABLES AND DIFFERENTIAL EXTENSIONS 8.4 THE BIG GUNS - STORAGE AND USEREASY 9. NONTRIVIAL EXAMPLES 9.1 A BIG EXAMPLE: THE MHD EQUATIONS 9.2 A SMALL NASTY EXAMPLE: THE NONLINEAR TELEGRAPH SYSTEM REFERENCES INDEX OF FUNCTIONS AND OPTIONS 1. INTRODUCTION The Standard Form Package consists of procedures for obtaining information about analytic systems of Partial and Ordinary Differential Equations (PDEs and ODEs) using completely algorithmic differentiation and algebraic processes. These procedures do not depend on the heuristic process of explicitly integrating such systems. The package can yield information about the number of solutions possessed by a system (if any), canonical forms for such systems and structural properties of symmetry groups of differential equations. Paradoxically such information sometimes leads to new ways of solving PDE systems [Reid91c], [Reid92b]. The algorithms of the package may also be used as precursors for the application of more heurisitic methods (e.g. canonically simplifying systems of PDEs so that existing heuristic or numerical methods have a better chance of solving such systems). We have applied the algorithms of this package for the most part to problems concerning symmetry analysis of differential equations. We stress however that the algorithms can be used to obtain information about analytic systems of PDEs regardless of their origin. The theoretical background of this package lies in the formal theory of differential equations, the algorithms of this theory and their application to symmetries of PDEs developed in [Reid90, Reid91a, Reid91b, Reid91c, Reid92a]. The general area of formal analysis of differential equations is strewn with mathematics of every variety. Amongst these: the abstract cohomological approaches of Spencer and his co-workers [Spenc69]; the algebraic work of Ritt and his descendents [Ritt50, Kol73]; Pommaret's geometric formal methods [Pomm78, Pomm90, Schu92]; the differential Buchberger algorithms of Carra-Ferro [Carra87], Ollivier [Olliv91], and Mansfield [Mans91]; approaches through the use of exterior differential systems [Cart45, BCGGG91, HaE71, HartTker91]; applications to Lie symmetry groups of differential equations [Reid90, Reid91b, Reid92a], [Sch92a, Sch92b]; Wolf's resultant algorithm [Wolf89]. This area lies at the intersection of geometrical, algebraic, analytical and algorithmic techiques. Formal theory concerns the existence of formal series solutions of PDE systems. By a formal series solution of a system we mean a series that satisfies the system assuming that all operations - differentiation and summation of series etc - can be carried out formally term by term. This theory has its roots in the works of Tresse [Tress94], Riquier [Riq10] and Janet [Jan20] and has been developed and applied by the workers above and many others. The heart of the Standard Form Package is the function standard_form, which is a function for canonically simplifying systems of PDEs and ODEs. The simplified system has the same formal solutions as the original system away from a certain degenerate set of points which are described by a set of pivots returned by the algorithm. One of the major ideas of the algorithm [Jan20], [Reid91a] is to split the set of all derivatives of a solution at a point in the space of independent variables of a given PDE system into two disjoint classes: parametric derivatives (those derivatives whose values at the point can be freely ascribed), and principal derivatives (those derivatives whose values at a point are uniquely determined through the system once the values of the parametric derivatives have been given). The dimension of the space of formal solutions at the point is equal the number of parametric derivatives and can be determined by an algorithm which does not depend on explicitly solving the system, once the system is in standard form. The collection of parametric derivatives at a point is called the parametric initial data for the system. If the system is an analytic function of its independent, dependent variables and derivatives at a point, the parametric initial data is analytic and the ordering used in the calculation is of certain type (e.g. dominated by total degree), then Riquier [Riq10, Thom29] showed that the corresponding formal solution about the given point is the unique analytic solution of the system. One of the major obstacles in implementing formal methods is that the size of the calculations grows explosively with the order of the derivatives in the systems of PDEs to which it is applied. While developing this package, we have been primarily motivated by our efforts in solving real problems, provided by our own interests, and those of other applied researchers we have worked with. Consequently, the program has been developed, not in a vacuum, but in a practical environment. We have attempted to make the existing code have enough flexibility to allow users to deal with the difficulties that encountered in applied problems (see, in particular, section 8, for ways of dealing with memory explosion). The function standard_form is most useful for simplifying systems of large numbers of linear PDEs, which are overdetermined or have many redundant equations. The strong point of the package is that it is algorithmic for linear systems of PDEs - no hit or miss heuristics are involved. The canonical form is guaranteed after a finite number of steps. For nonlinear systems the package is not fully algorithmic but may partially simplify such systems. Buchberger's algorithm has provided a method for canonically simplifying polynomially nonlinear systems of algebraic equations. Although valuable progress has been made for polynomially nonlinear systems of PDEs [Carra97, Mans91, Olliv91], as yet no fully effective differential analogue of Buchberger's algorithm is known. The procedure standard_form can be regarded as a generalization of the Gaussian Reduction algorithm for systems of linear algebraic equations to systems of linear PDEs (to which it reduces if only linear algebraic equations are involved). If you are dealing with a system of linear PDEs then the procedure standard_form can in principal guarantee you: (1) determination of whether the system is consistent or inconsistent (i.e. has no formal or analytical solutions). (2) reduction of the system to an ordered canonical form using the procedure standard_form; (3) determination of the the dimension of the solution space of the system by using the procedure initial_data; (4) determination of the Taylor series solution of the system to a given finite degree using the procedure taylor_exp; at points which do not belong to a certain degenerate set of points described by a set of equations (the "pivots") which is also returned by standard_form. The second part of the package, involves powerful algorithms for the symmetry analysis of PDEs. In particular our packages enable the dimension and structure of Lie symmetry algebras of differential equations to be directly calculated from their determining systems, without needing to integrate such systems [Reid91b, Reid92a]. The conventional approach to infinitesimal symmetries of differential equations has been to try to integrate the determining equations, explicitly determining the basis generators of the algebra, and then substituting the results into the commutation relations above to determine the structure constants [BlK89, Olv86, Ovs82]. Since this process is based on the heuristic process of integration it is not always guaranteed to succeed. Our approach to finding the structure of Lie symmetry algebras of PDEs is algorithmic: (1) The determining system for the infinitesimals by one of the computer packages designed for this purpose (see [Here93, Champ91]). (2) Algorithm standard_form is used to reduce the determining system to standard form. (3) Application of initial_data to the standard form produced in Step (2) algorithmically determines the dimension of the algebra of infinitesimal symmetries, without having to determine the explicit form of these generators. (4) Algorithm fin_com_table uses the standard form from Step (2) and its initial data from Step (3) to determine the structure of the of the finite-dimensional subalgebra of infinitesimal symmetries corresponding to the parameters in the finite initial data. This algorithm depends only on differentiation and substitution. (5) taylor series approximations of the infinitesimals or basis generators can be found to any given finite degree using algorithms taylor_exp, or finite_generators respectively. All of Section 2 is essential reading. To get started, those interested in symmetry analysis of PDEs should carefully cover all of sections 2 and 4.1. 2. GETTING STARTED WITH THE PACKAGE 2.0 HOW TO USE THIS MANUAL The main intention of this manual is to provide a low-level, example-based guide to our package, which assumes a minimum of mathematical and computer-algebra background, but tries to sneak in examples which illustrate some of the ideas behind the package. Consequently the flavour is loose and sacrifices precision in favour of analogy. A computer-oriented Maple-styled short guide is planned together with works that concentrate on the mathematics involved. A helpful way to learn about the programs is to type (or cut and paste) the examples we present (or a similar example of your choice) into a Maple session. We encourage experimentation. Many examples are simple enough to work by hand. Nontrivial examples are given in Section 9. We welcome your suggestions, criticisms and additions. The aim of the manual is to provide easy continuous reading. It is written in askey so that the find option of your editor can be used to locate all occurrences of that topic (you may find it helpful to first look at the table of contents, or at the index of functions and options at the end of the manual). 2.1 SIMPLE EXAMPLES EXAMPLE A Consider the system of nonlinear PDEs: y Zxxx - x Zxyy = Zyy - y Zy 2 2 2 2 y x Zxxx Zxyy + x Zxxx + x y Zxyy = 0 2 2 y Zxyy - x W + 2 x y Z = 0 2 2 Zyy - y Zy + 2 x y W = x W where the dependent variables W and Z are functions of the independent variables x and y, and Zxxx denotes the third partial derivative of Z with respect to x etc. We start a Maple session (Maple 5.0 or higher version) and read in the Standard Form Package: > read sfprogs; (Maple command lines begin with > and are terminated with ; or :). To use the programs we input the system as a list of equations with independent variables x1,x2,x3, etc. and dependent variables V(1), V(2), ... etc. So denoting x by x1, y by x2, W(x,y) by V(1) and Z(x,y) by V(2) the equations above are written as a list: > sys:= [x2*V(2,1,1,1)-x1*V(2,1,2,2) = V(2,2,2)-x2*V(2,2), 2*x2*x1*V(2,1,1,1)*V(2,1,2,2) + x1*V(2,1,1,1)^2 + x1*x2^2*V(2,1,2,2)^2 = 0, x2*V(2,1,2,2) - x1*V(1) + 2*x1^2*x2*V(2)^2 = 0, V(2,2,2) - x2*V(2,2) + 2*x1^2*x2*V(2)^2 = x1*V(1)]: where we have terminated the command line with : instead of ; to suppress the list sys from being displayed. Derivatives are denoted by using a repeated index notation specific to the program. For example, V(2,2) denotes the x2 partial of V(2), and, V(2,1,2,2) denotes the partial derivative with respect to x1 of the x2 partial of the x2 partial derivative of V(2) (unlike Maple's usual notation for differentiation). Notice that Maple is by default, case sensitive, so that x1 is distinct from X1 etc. The program depends on a complete ordering which is defined on derivatives. By default the ordering is the following standard one, denoted >s: (1) First ordered by total order of the derivative. (2) Then lexiographically by derivative: (d/dx1 >s d/dx2 >s ... >s d/dxm ) (3) Then by function number (V(1) >s V(2) >s ... >s V(n)) Other orderings are possible and assigning different orderings can be useful for reasons of efficiency or for forcing the system or its equations to take certain pleasant forms [Reid92b]. Hence, in the standard ordering, the derivatives in our example are ordered V(2,1,1,1) >s V(2,1,2,2) >s V(2,2,2) >s V(2,2) >s V(1) >s V(2) There is a unique derivative (or dependent variable) which is largest in each equation of a system with respect to the above ordering. To reduce this system to standard form we use the procedure standard_form: > sf:= standard_form(sys); and the resulting output is: The system parameters are as follows: There are 2 independent variables: [x1, x2] There are 2 dependent variables: [V(1), V(2)] There are no classification variables: sf := [V(2, 1, 1, 1) = 0, V(2, 1, 2) = 0, V(2, 2, 2) = x2 V(2, 2), 2 V(1) = 2 x1 x2 V(2) ] The standard form above will have the same formal solutions as the original system sys at points where none of the pivots vanish. The pivots for the calculation are stored in a global variable _pivs and are: > _pivs; 2 2 x2 x1 (x2 + x1 + x2) [x1 = 0, - x2 - x1 = 0, - -------- = 0, - x2 = 0, - ------------------ = 0] 2 x2 x2 + x1 Consequently the calculation is valid at all points (x1,x2) with x1#0, x2#0, x2^2+x1#0, x2^2+x1+x2#0. The pivots e.g. x1, -x2^2 - x1 etc above are quantities which appear in divisions during the computation (for nonlinear systems they may involve dependent variables and their derivatives). It may also yield valid results at points other than these, but this conclusion can not be deduced from the output standard form, sf, above. The sfprogs command shortenpiv will remove repeated factors from the list _pivs (in the order that they occurred): > spivs:= shortenpiv(_pivs); 2 2 spivs := [x1, x2 + x1, x2, x2 + x1 + x2] In our original notation the (considerably) simplified system is: 2 Zxxx = 0, Zxy = 0, Zyy = y Zy, W = 2 x y Z Once a system of PDEs sys is in standard form then its initial data can be calculated using algorithm initial_data. As mentioned in the Introduction the initial data consists of a set of assignments of values to parametric derivatives at an initial point which is denoted by X=(X1,X2,...,Xm). In particular if V(r) is one of these parametric derivatives, with r a multiindex, r=(i,j,...k), and b is the value of this derivative at the initial point x=X (i.e. x1=X1, x2=X2, ... , xm=Xm) then the corresponding initial condition has the form: V(r)| = b x=X For simplicitly of display we will write the above condition simply V(r) = b (with evaluation at x=X being understood) and denote the different arbitrary constants (b's) by a1, a2, a3, ... etc. For example applying the procedure initial_data to the standard form of the system above gives: > id:= initial_data(sf); yielding id := [[V(2) = a1, V(2, 2) = a2, V(2, 1) = a3, V(2, 1, 1) = a4], []] The initial data list id above is equivalent to the initial conditions: V(2)| = a1, V(2,2)| = a2, V(2,1)| = a3, V(2,1,1)| = a4 x=X x=X x=X x=X posed at x=X (i.e. x1=X1, x2=X2) where a1, a2, a3 and a4 are arbitrary. Or in terms of our original notation (with (X1,X2)=(x0,y0)) Z(x0,y0) = a1, Zy(x0,y0) = a2, Zx(x0,y0) = a3, Zxx(x0,y0) = a4 Since the initial data contains 4 arbitrary constants the general solution of the system depends on 4 parameters and the solution space is 4-dimensional. If the solution space is infinite dimensional then there would be an infinite number of initial conditions of the form above. Writing these down in a finite manner requires the clever packing of these infinities into a finite number of analytic functions (see the wave equation example B later in this section and section 3.2). Once we have computed the standard form of the original system, and computed the initial data, the taylor_series option will compute the formal taylor series solution of the system to any finite degree N, with the command > taylor_exp(sf,id,N) where sf and id are the standard form and initial data of the system, respectively. The output will be the Nth degree Taylor polynomial for the dependent variables about an initial point X = (X1,...,Xm). In these polynomials hi = xi - Xi, i=1...,m and their coefficients are functions of the the initial point X and the initial data. To form the taylor series to degree 2 about a general initial point (X1,X2) for the above example we use the command: > ts:= taylor_exp(sf,id, 2); which yields 2 2 2 2 ts := [V(1) = 2 X1 X2 a1 + (2 X1 X2 a1 a2 + 4 X1 a1 a2 + 2 X1 X2 a2 ) h2 2 2 + (4 X1 X2 a1 a3 + 2 X2 a1 ) h1 + (2 X1 a1 + 4 X1 X2 a1 a2) h2 2 2 + (2 X1 X2 a1 a4 + 4 X2 a1 a3 + 2 X1 X2 a3 ) h1 2 + (4 X1 X2 a2 a3 + 2 a1 + 4 X2 a1 a2 + 4 X1 a1 a3) h1 h2, 2 2 V(2) = a1 + a3 h1 + a2 h2 + 1/2 a4 h1 + 1/2 X2 a2 h2 ] where h1 = (x1 - X1) and h2 = (x2 - X2) (note that the initial data must be in the format presented by the procedure initial_data, so that a1, a2, ... are not given particular values). So the taylor series solution about x0=0, y0=0 to order 2 for our original system is 2 W(x,y) = 2 a1 x y + 0(3) 2 Z(x,y) = a1 + a3 x + a2 y + a4 x /2 + 0(3) Note that the coefficients are expressed entirely in terms of initially specified information: the initial point, X and the initial conditions, a1,...a4. Hence, given the initial data, we can always compute the taylor series solution to any finite degree in terms of these values. In general the initial data specified by algorithm initial_data is sufficient to uniquely determine formal solutions of systems of PDEs. For those familiar with the underlying theory, this is equivalent to the theorem by Riquier and Janet [Riq10, Jan20] that if a system is in passive form (i.e, all integrability conditions are satisfied) then specifying values for the parametric derivatives uniquely specifies values for the remaining principal ones. If the standard form sf is an analytic function of its arguments and the dimension of the solution space is finite then this series will converge to a unique analytic solution of the system in some neighbourhood of the initial point. EXAMPLE B Consider the wave equation, Uxx - Utt = 0 written in Maple's usual derivative notation: > msysb:= [diff(U(t,x),t,t) - diff(U(t,x),x,x) = 0]; It can be converted to the repeated index notation used by the program by hand or by using the program described in Section 6. Applying the program of Section 6 (which works in Maple 5 release 2 or later): > read convprogs; > sysB:=map2rep(msysB,[u],[t,x]); sysB := [V(1, 1, 1) - V(1, 2, 2) = 0] Taking standard_form of this equation gives: > sfB:= standard_form(sysB); sfB := [V(1, 1, 1) = V(1, 2, 2)] and initial data > idB:= initial_data(sfB); idB := [[], [V(1) = f1(x2), V(1, 1) = f2(x2)]] This corresponds to the usual initial data for the wave equation U(t0,x) = f1(x), Ut(t0,x) = f2(x) Notice how this really corresponds to an infinite sequence of initial conditions of the form V(r)| = b: x=X U(t0,x0) = f1(x0), Ux(t0,x0) = f1'(x0), Uxx(t0,x0) = f1''(x0), ... Ut(t0,x0) = f2(x0), Utx(t0,x0) = f2'(x0), Utxx(t0,x0) = f2''(x0), ... so the program has a clever of packing the infinite number of initial conditions into the functions f1(x) amd f2(x). In general the initial data id is presented as a list which has 2 sublists. The first entry of id (i.e. the list id[1]) contains a list of a finite number (J say) of initial conditions each of form V(r) = aj, j=1,...,J, The second list contains K (infinite) initial conditions each of form: V(s) = fk(y) where y is a sublist of [x1,x2,...,xm] If the second sublist is empty (i.e. id[2]=[]) then the solution only depends on the finite number of parameters a1, a2, ... ,aJ. A more detailed description (especially of the infinite data case) is given in section 3.2. 2.2 SYSTEMS OF ALGEBRAIC EQUATIONS In this section we apply the programs to some simple systems of algebraic equations. Maple's solve command and Grobner package are far more suitable for treating such systems of algebraic (non-differential) equations than standard_form. We consider this case first to offer the user some informal insight into the workings of our package, and how it deals with nonlinear equations. We also explain the options divpiv and the global list _pivs, which records pivoting operations of the package. This discussion should be helpful for understanding the differential equations case. Example A We enter a system of linear algebraic equations. As usual the equations must be written in terms of dependent variables V(1),V(2),... which are assumed to be functions of the (independent) variables x1,x2, ... that appear in the equations. So the system below is regarded as a linear system for V(1),V(2),V(3),V(4) which are functions of x1, with k being a parameter. > read sfprogs; > sysA:= [ 5*k*V(1) - 2*V(2) = 3*V(3)-V(4), 6*k*V(1) + 5*V(2) - x1*V(4) = 0]: We now find the standard form: > sfA:= standard_form(sysA); (- 5 + 2 x1) V(4) 15 V(3) sfA := [V(1) = 1/37 ----------------- + ---- ----, k 37 k 18 V(2) = (6/37 + 5/37 x1) V(4) - ---- V(3)] 37 The standard form algorithm solved for those equations that are highest in the standard ordering, as described in 2.1. Since V(1) >s V(2) >s V(3) >s V(4), it solved the two equations for V(1) and V(2). The initial_data algorithm yields: > idA:= initial_data(sfA); idA := [[], [V(3) = f1(x1), V(4) = f2(x1)]] so that V(3) and V(4) are arbitrary functions of x1, which makes sense since there are no restrictions on V(3) and V(4) in the above standard form. The input system contained a parameter k so we are dealing with a class of equations and the standard form obtained may not be valid for exceptional values of k. The way the standard form was achieved was through the standard gaussian elimination algorithm. Thus the program first solved the first equation for V(1) (so that k arose as a pivot), subsituted into the second equation, and then solved the second equation for V(2). However, if a pivot is equal to zero, then the algorithm has commited an error, since it has divided by it. Such cases require separate treatment. To detect such exceptional classification cases the pivots are stored in a global list called _pivs. > _pivs; [5 k = 0, 37/5 = 0] So we see that exceptional cases occur only if any of the pivot equations are satisfied. Standard Forming the only exceptional case k=0 > sfA1:= standard_form(subs(k=0, sysA)); yields sfA1 := [V(2) = 1/5 x1 V(4), V(3) = (1/3 - 2/15 x1) V(4)] with initial data > idA1:= initial_data(sfA1); idA1 := [[], [V(1) = f1(x1), V(4) = f2(x1)]] and pivots > _pivs; _pivs := [-2 = 0, -15/2 = 0] so that no further branching is possible. In the PDE case the program has sophisicated classification capabilities and these are discussed in Section 3.5. Example B We consider the system of nonlinear algebraic equations: > sysB := [ x1*V(1) - log(x2)*V(2)^2 = 0, V(3)^2 - 2*V(3)*V(4) + V(4)^2 = 0, V(3)^4*V(2)^2 - 2*x4*log(V(4))*V(3)^2*V(2) + x4^2*log(V(4))^2 =0]: > sfB:= standard_form(sysB); yields: ****************************************************** ********** The system is not fully reduced *********** ** The non-linear equations are in the second list: ** ****************************************************** 2 ln(x2) V(2) sfB := [[V(1) = ------------, V(3) = V(4)], x1 4 2 2 2 2 [0 = - V(4) V(2) + 2 x4 ln(V(4)) V(4) V(2) - x4 ln(V(4)) ]] with pivots: > _pivs; [x1 = 0, 1 = 0] The standard ordering >s on the dependent variables, as given in 2.1, is: V(1) >s V(2) >s V(3) >s ... which allows us to identify the largest variable (wrt >s) in any equation. If an equation is of form N p (V(i) - r) = 0 (*) with largest variable V(i) where p is a function only of the independent variables and r is not a function of V(i) (N a nonzero positive integer), then procedure standard_form will solve the equation for V(i) to obtain V(i) = r (**) Thus in Example B it solved the first equation for V(1) and since the second equation was (V(3) - V(4))^2 = 0, it solved this equation to obtain V(3) = V(4). It did not solve the last equation, even though this equation has form (*) 4 / x4 ln(V(4)) \ 2 V(3) | V(2) - ----------- | = 0 | 2 | \ V(3) / since p (= V(3)^4) was a function of a dependent variable. Consequently the system was not fully reduced. ( Algorithm initial_data can not be applied here since the system is not in standard form.) If the option divpiv is added to the argument list of standard_form then the procedure will solve equations of form (*) to obtain (**) even if p depends on some of the dependent variables and their derivatives, provided they do not depend on V(i) or derivatives of V(i). Thus for our example: > sfB1:= standard_form(sysB, divpiv); fully reduces the system 2 2 ln(x2) x4 ln(V(4)) x4 ln(V(4)) sfB1 := [V(1) = --------------------, V(2) = -----------, V(3) = V(4)] 4 2 V(4) x1 V(4) with pivots: > _pivs; 4 [x1 = 0, V(3) = 0, 1 = 0] alerting us that V(3) = 0 is a special case that needs separate treatment. Since the system is now in standard form, algorithm initial_data can be applied: > id:= initial_data(sfB1); which yields: id := [ V(4) = f2(x1,x2,x3,x4) ] It's easy to see that allowing division by such pivots is a dangerous operation and so we've made divpiv an option. By returning a system of unsimplified leading nonlinear equations the opportunity is open to then use other packages to simplify these nonlinear equations (e.g. Maple's algebraic Grobner package or Mansfield's package [Mans91]). Example C. Consider > sysC := [2*V(1) - a*x1 = 0, V(1) - x1 = 0]; > standard_form(sysC); leads to the error message: Error, (in maxder) Inconsistent equation in maxder:, 0 = 1/2*x1*(a-2) which indicates that if a#2 then the system is inconsistent and has no solution (no relations are allowed amongst the independent variables). The case a=2 has to be considered separately. Example D. Consider > sysD:= [ V(1) - 4*x1 = 0, V(2)/(V(1)-4*x1) = 0]: > sfD:= standard_form(sysD); leads to the error message: Error, (in clearderivdep) division by zero which indicates the system has no solutions for obvious reasons. 2.3 NOW YOU'RE READY FOR SYSTEMS OF DIFFERENTIAL EQUATIONS The algorithm standard_form can be applied to systems of PDEs with independent variables x1,x2, ... , xm and dependent variables V(1),V(2), ... , V(n). We give an informal description through examples of some of the major aspects of the program. Example A. Consider the system of PDEs Z - Wxy + y Wyy - xW = 0 (1) Zx - Wxxy + y Wxyy - W - x Wx = 0 (2) for W(x,y) and Z(x,y). The largest derivative in equation (1) is Wxy and the largest derivative in equation (2) is Wxxy. The first stage of the standard_form algorithm reduces a system as much as possible by a gaussian-elimination type process to a solved-form consisting of equations of the form derivative = rhs where the rhs derivatives are lower in the ordering than the derivative on the lhs. At the same time the program records the pivots in the global list _pivs. After this stage no derivative can appear explicitly on both the LHS and RHS of the system, and derivatives only appear once on the LHS. For the above system this triangular form is accomplished by solving the first equation for Wxy and the second equation for Wxxy giving: Wxy = Z + y Wyy - x W (1)' Wxxy = Zx + y Wxyy - W - x Wx (2)' This is the maximum elimination possible if we regard Wxy and Wxyy as completely different algebraic variables. However they aren't independent since Wxxy = (Wxy)x and the system above is equivalent to (has the same solutions as) the system: (1), (2) - d/dx1 (1). The second equation is an identity, so the system reduces to just one equation: Wxy = Z + y Wyy - x W And this is the standard form of the system. To carry out the calculation in Maple we read the package sfprogs > read sfprogs; and let W=V(1), Z=V(2), x=x1, y=x2, so that the system for example A is: > sysA:= [ V(2) - V(1,1,2) + x2*V(1,2,2) - x1*V(1) = 0, V(2,1) - V(1,1,1,2) + x2*V(1,1,2,2) - V(1) - x1*V(1,1) = 0]; and then > sfA := standard_form(sysA); yields sfA:= [ V(1,1,2) = V(2) + x2 V(1,2,2) - x1 V(1) ] agreeing with the result obtained above. Example B. In general such differential substitutions, as in the previous example will lead to some nontrivial relations amongst the derivatives of the system. These new equations can be added to the original equations, and the resulting system can again be triangularized. After a finite number of steps the system will be obtained in a form where no derivative on the LHS of the system is a derivative of any other derivative on the LHS of the system. For example: Px = P Py = x P where P=P(x,y). However this system may still not be fully reduced. We see that for consistency we need the integrability condition: (Px)y - (Py)x = 0 (P)y - (x P)x = 0 (modulo system) Py - P - x Px = 0 We substitute for Px and Py from the first two equations and the new equation reduces to P = 0. We add this to the system. But with this equation, the whole system reduces to P = 0 and this is the standard form of the system. In maple we have > sysB:= [V(1,1)=V(1), V(1,2)=x1*V(1)]; for the above system (with x=x1, P=V(1)) and applying standard_form > sfB:= standard_form(sysB); yields sfB:= [ V(1) = 0 ] as above. We encourage the reader to experiment with their own examples to get a feel for the algorithms. Example C. Consider the linear system of PDEs: Uxt - r x t = 0, Uxx - Ux = 0 Where U = U(x,t). Putting U = V(1), x=x1, t=x2 then > read sfprogs; and the system is entered as: > sysC:= [ V(1,1,2) - r*x1*x2 = 0, V(1,1,1) - V(1,1) = 0]; Applying standard_form: > sfC1:= standard_form(sysC); gives the message: Error, (in bigstandard) System is inconsistent:, r*x2*(-1+x1) = 0 Written in the original variables, the condition rt(-1+x) = 0 has resulted from taking the integrability condition (Uxt)x=(Uxx)t, after the algorithm has reduced it to the solved form Uxt = r x t , Uxx = Ux So the system is inconsistent (i.e. has no solutions) unless r=0. When r=0 > sfC2:= standard_form(subs(r=0, sysC)); sfC2 := [V(1, 1, 1) = V(1, 1), V(1, 1, 2) = 0] and the initial data is > idC2:= initial_data(sfC2); idC2 := [[V(1, 1) = a1], [V(1) = f1(x2)]] so that the general solution depends on one arbitrary function of x2 and one arbitrary constant (i.e. the solution space is infinite-dimensional). In particular the initial data is V(1,1)| = a1 and V(1)| = f1(x2) x1=X1,x2=X2 x1=X1 where a1 is an arbitrary constant and f1 is an arbitrary function of x2. Written in our original notation these initial conditions are: Ux(x0,t0) = a1 and U(x0,t) = f1(t) where (x0,t0) is a fixed point. Also > _pivs; [-1 = 0, 1 = 0] so that there are no nontrivial pivots. The taylor series expansion of the solution to third order is > ts:= taylor_exp(sfC2, idC2, 3); / d \ 2 ts := [V(1) = f1(X2) + a1 h1 + |----- f1(X2)| h2 + 1/2 a1 h1 \ dX2 / / 2 \ / 3 \ | d | 2 3 | d | 3 + 1/2 |------ f1(X2)| h2 + 1/6 a1 h1 + 1/6 |------ f1(X2)| h2 ] | 2 | | 3 | \ dX2 / \ dX2 / which is the taylor solution about a fixed initial point (X1,X2) giving the values of V(1) at X1 + h1, X2 + h2. There are some nonzero third order terms in the series so the series may not terminate at this order (it would have terminated if there had been been no third order terms - see Section 5.1). Example D. Consider the nonlinear system of PDEs for U=U(x,t) 2 Uxt = 0, Ux Uxx - Ux - Ut Uxx + Ut Ux = 0 that is putting U=V(1), x=x1, t=x2 > sysD:= [ V(1,1,2) = 0, V(1,1)*V(1,1,1)-V(1,1)^2 - V(1,2)*V(1,1,1)+V(1,2)*V(1,1) = 0]; > sfD:= standard_form(sysD); ****************************************************** ********** The system is not fully reduced *********** ** The non-linear equations are in the second list: ** ****************************************************** sfD := [[V(1, 1, 2) = 0], 2 [0 = - V(1, 1)V(1, 1, 1) + V(1, 1) + V(1, 2)V(1, 1, 1) - V(1, 2) V(1, 1)]] does not reduce the system completely to standard form. Adding the option divpiv forces standard_form to solve for V(1,1,1) in the second equation and the system is completely reduced: > sfD1:= standard_form(sysD, divpiv); sfD1 := [V(1, 1, 1) = V(1, 1), V(1, 1, 2) = 0] Great care must be taken when dealing with nonlinear systems of PDE. In particular the pivots should always be examined to identify special cases. Here > _pivs; [-1 = 0, V(1, 1) - V(1, 2) = 0] So a special case occurs if V(1,1) - V(1,2) = 0 (reflecting the fact that the sfD1 was obtained by the dangerous operation of dividing by the pivot V(1,1) - V(1,2) in the second equation to solve for V(1,1,1)). To treat this case we append this equation to sysD > sysD2:= [op(sysD), _pivs[2]]; sysD2 := [V(1, 1, 2) = 0, 2 V(1, 1) V(1, 1, 1) - V(1, 1) - V(1, 2) V(1, 1, 1) + V(1, 2) V(1, 1) = 0, V(1, 1) - V(1, 2) = 0] and reduce the resulting system to standard form: > sfD2:= standard_form(sysD2); sfD2 := [V(1, 2, 2) = 0, V(1, 1) = V(1, 2)] > _pivs; [-1 = 0, 1 = 0, 1 = 0] so there are no further special cases. 2.4 IMPORTANT NOTE ABOUT OPTIONS FOR standard_form The function standard_form can be used with a list of options: > standard_form(sys, optionlist); You've already met one of those options (divpiv) and many more will be described throughout the manual. Any number of the options can be used in any order. Effective use of the program will often entail using quite a number of them. For example > standard_form(sys, class_eqns=[[V(4),3], V(4,3,3) = V(4)], store, strategy = [1.75], divpiv); where the 3 options above can be given in any order. See the index of functions and options for a list of these options. 3. MORE ADVANCED ASPECTS OF THE PACKAGE 3.1 HOW TO CHANGE THE ORDERING - the option janet_order The standard form of a system depends on the ordering chosen. For purposes of efficiency, or to help expose simple equations in the system it is useful to be able to change the ordering. The ordering on derivatives must be left invariant under differentiation. A wide class of orderings having this property was given by Riquier [Riq10] (also see Janet [Jan20]). Such orderings can be characterised by a matrix R with m+n columns (m=#indep vars and n=#dep vars), with nonnegative integer entries. There are orderings beyond those of Janet - see [CarraSitt93], [Rust93], [Weisp93]. To define such janet orderings we represent derivatives of V(p): a1 a2 am (d/dx1) (d/dx2) ... (d/dxm) V(p) by the corresponding vectors [a1, a2, ... , am, 0,...,0,1,0,...,0] (m+p)-th (e.g. V(2,1,1,2) where m=2, n=2 is represented by [2,1,0,1]). The usual lexicographic ordering >lex on vectors is defined by: p=[p1,p2, ... , pk] >lex q=[q1,q2, ..., qk] iff the first nonzero pj - qj is positive, that is, the first nonzero element of the vector p-q is positive. The ordering on derivatives > defined by a matrix R is given by: R a1 a2 am b1 b2 bm (d/dx1) (d/dx2) ... (d/dxm) V(p) > (d/dx1) (d/dx2) ... (d/dxm) V(q) R iff T T R v >lex R w where T denotes transpose and v=[a1,...,am,0,...,1,...,0], w=[b1,...,bm,0,...,1,...,0] m+p m+q The standard ordering given in Section 2.3 has ordering matrix x1 x2 x3 ... xm V(1) V(2) V(3) .... V(n) 1 1 1 1 0 0 0 0 1 0 0 0 0 0 R = 0 1 0 ...... 0 0 0 1 ...... 0 0 0 0 ...... 0 0 0 0 1 ...... 0 0 0 0 0 n (n-1) (n-2) 1 Example A. New ordering matrices are entered as list of lists option to standard_form (i.e. not in maple's matrix or array format). For the heat equation: > sf:= standard_form([V(1,2)-V(1,1,1) = 0]); yields sf:= [V(1,1,1)=V(1,2)] in the standard ordering. To execute standard_form with nonstandard orderings these are entered as an option to the standard_form procedure, in the form: janet_order = list of lists where each list has length m+n (m=#indeps n=#deps) For example to reduce the heat equation in in the ordering defined by the ordering matrix x1 x2 V(1) 0 1 0 R = 1 0 0 0 0 1 (i.e. ordering first by the independent variable x2) the command is: > sf:= standard_form([V(1,2)-V(1,1,1)=0], janet_order=[[0,1,0], [1,0,0], [0,0,1]]); yields sf:= [V(1,2) = V(1,1,1)] Entering at least one row is enough to define a complete ordering. If your entered ordering fails to distinguish between 2 derivatives the program will use the default ordering to break the tie. For example, above, we could have obtained the result above by: > sf:= standard_form([V(1,2) - V(1,1,1) = 0], janet_order=[[0,1,0]]); that is sf:= [V(1,2) = V(1,1,1)] The effect on systems of PDEs by using new orderings can be dramatic. The most useful (but computationally most expensive orderings) are the lexicographical orderings which order first by function number rather than order of derivative: The first row of R is then [0,0,...,0,perm] where perm is some permutation of the integers {1,2,...,n}. The lexicographic orderings have the effect of sequentially decoupling a system of PDEs. An interesting extension of the work in [Reid91c] would be to investigate how orderings affect the directed graph of a PDE system. 3.2 CALCULATION OF INITIAL DATA AND TAYLOR SERIES SOLUTIONS OF SYSTEMS OF PDE Once a system of PDEs sys is in standard form > sf:= standard_form(sys); then its initial data can be calculated using algorithm initial_data: > id:= initial_data(sf); The first entry of id (i.e. the list id[1]) contains a list of J (finite) initial conditions each of form V(r) = aj which can be interpreted as: V(r)| = aj x=X i.e. the derivative V(r) evaluated at the fixed point X=(X1,X2,...,Xm) has the value aj. The second list contains K (infinite) initial conditions each of form: V(s) = fk(y) where y is a sublist of [x1,x2,...,xm] which can be interpreted as: V(s)| = fk(y) H i.e. the derivative V(s) evaluated on the hyperplane H takes the values fk(y) where fk is an arbitrary function of y and H:= {(x1,x2,...,xm)| xj = Xj, for xj not in y, xj arbitrary for xj in y}. Example A. For the heat equation Uxx - Ut = 0: > read sfprogs; > sf1:= standard_form([V(1,1,1) - V(1,2) = 0]): yields sf1:= [V(1,1,1) = V(1,2)] and > id1:= initial_data(sf1); yields the initial data id1:= [[],[V(1) = f1(x2), V(1,1) = f2(x2)]] that is U(x0,t) = f1(t), Ux(x0,t) = f2(t) and there are no finite initial conditions since the first list is empty. To obtain the usual initial data for the heat equation we have to use a nonstandard ordering (see 3.1) > sf2:= standard_form([V(1,1,1)-V(1,2)=0], janet_order=[[0,1,0]]); yields sf2:= [V(1,2) = V(1,1,1)] and > id2:= initial_data(sf2); gives id2:= [[], [V(1) = f1(x1)]] Which is the usual initial data for the heat equation (i.e. U(x,t0) = f1(x)). See 2.2 for an example where only finite initial conditions occur and Section 2.3 Example C, for a case where both finite and infinite initial conditions occur. Once the standard form of a system of PDEs sys has been determined and the initial data calculated from sf then its taylor series expansion about a point (X1,X2,...,Xm) in powers of h1(=x1-X1),...,hm(=xm-Xm) to an order N can be calculated by algorithm taylor_exp from sf and id: > read sfprogs; > sf:= standard_form(sys); > id:= initial_data(sf); > ts:= taylor_exp(sf, id, N); Example B. The taylor series expansion for the heat equation to degree 4 is given by: > ts:= taylor_exp(sf2, id2, 4); where sf2 and id2 have been calculated previously (above). Note: at present taylor_exp though simple is very inefficient and it would be worthwhile to develop strategies to increase its efficiency. 3.3 THE STRATEGY OPTION The strategy option allows the user to control many details of the operation of standard_form. Invoking the strategy option involves entering a list, at the same time as a call to the standard form program, whose entries are new values for certain parameters that control the operation of standard form. These parameters include one that control the amount of intermediate information that standard form displays as it runs (infolevel below), and also basic decisions regarding how many new, unclassified equations it should attempt to process in each iteration (controlled by the values of p1 and p2 below). To change any of the parameters in strategy from their default settings, you enter a new strategy list as an argument to standard_form, in the following way > sf:=standard_form(sys,strategy=list) where list is one of the 3 possibilities: [infolevel] [infolevel, p1, p2] [infolevel, p1, p2, inconsopt, simpopt] i.e. list has 1, 3 or 5 entries. Also infolevel, p1 and p2 have numerical values and inconsopt, simpopt are boolean-valued. default setting: atrategy = [1.0, 1.0, 2.0, true, true] 3.3.1 DISPLAYING INTERMEDIATE INFO USING STRATEGY The first entry of the strategy list, regardless of whether it contains 1,3 or 5 entries, controls the displaying of intermediate information about the current run. At any time during the running of standard form, the equations in the system will be in one of 4 classes: one-term eqns Easy eqns Leading Nonlinear eqns Unclassified eqns OT Ez NL UnC The leading nonlinear equations are those equations which are nonlinear in their highest derivative. The program tries to take advantage of one term equations as much as possible (since they always lead to simplification), and then uses a criteria (the Easy-Unclassified choice) for choosing an easy subset of the full set of equations (the easy equations). Nasty (because of their length or nonlinearity) equations are stored in an unclassified equations list. Usually the algorithm has terminated when there are no unclassified equations remaining. The program has sophisticated strategies for avoiding doing and redoing integrability conditions. Integrability conditions are executed once at each iteration. strategy = [infolevel] (or [infolevel,p1,p2] or [infolevel,p1,p2,inconsopt,simpopt]) where: infolevel < 1 no intermediate information will be displayed, except the usual Maple memory and time reports. 1 <= infolevel < 1.25 (default setting) standard_form will report the number of equations in each class at the end of each basic iteration 1.25 <= infolevel < 1.5 standard_form will also report the maximum length of equations in the Ez and UnC categories, and the number of equations in the UnC category that have been classified and put into the OT, Ez and NL categories. 1.5 <= infolevel < 1.75 in addition to the above info standard_form gives a glossary of the symbols in the output and gives a more detailed reporting of the steps of the algorithm 1.75 <= infolevel in addition to the above info standard_form gives detailed info about all its major steps (this is GR's favourite level and a good way to undestand the algorithm). Of course you can also use Maple's printlevel command. 3.3.2 CONTROLLING EQUATION FLOW USING STRATEGY'S KNOBS The behaviour of the program can be changed by twisting two knobs; the p1 knob and the p2 knob (that is by giving the parameters p1 and p2 different values) where: strategy=[infolevel, p1, p2] These (`equation flow') parameters control the the quantity and relative size of the equations being chosen from UnC and placed into the easy list Ez. Special Settings default setting: strategy = [1, 1, 2] (i.e. p1=1, p2=2) Ultra-Conservative setting: strategy=[infolevel, 1, 0.9] (i.e. p1=1, p2=0.9) Good for large or complicated systems (transfers 1 equation at a time from UNC to Ez). Ultra-HOG setting: strategy=[infolevel, 100000, 100000] (i.e. p1=100000, p2=100000) Only use Ultra-HOG on very simple systems, when you're greedy for a quick run (transfers everything from UnC to Ez). The p1 knob ----------- Our first parameter p1 controls an upper bound for the number of equations being placed into Ez relative to the number already present. This bound is defined as a function of p1 as follows: # equations currently Max # equations taken from UnC & put in Ez in Ez (j) bound(j) ----------------- ------------------------------------------ 0 <=j <=15 10 x p1 16 <=j <=25 8 x p1 26 <=j <=35 6 x p1 36 <=j <=40 5 x p1 41 <=j <=45 4 x p1 45 <=j 3 x p1 In our experience the above relationship has worked fairly well for large systems, as the memory and time required for simplification of larger systems behaves roughly as the product of the number of equations in Ez and the number of new equations. Setting p1=10000 will mean that the number of equations currently in Ez will be ignored, and the maximum number of equations allowed by the p2 knob will be transferred from UnC to Ez at each iteration. Setting p1=0 means that only 1 equation will be transferred from UnC to Ez at each iteration. As the p1 knob is turned up the number of equations transferred from UnC to Ez becomes less sensitive to the currrent number of equations in Ez. Changing the way Ez equations are chosen from UnC by altering the numbers (15,25, ...etc) in the above table - can be accomplished by changing the these numbers in the following piece of sfprogs code (par[7] is p1): $$$$$$ start of section of sfprogs code $$$$$$ # # Here we select which Hard equations are to be put in Easy # t2:=array[1..4]; if nops(Hard)>0 then if nops(Easy)<=15 then numeqn:=10*par[7]; elif nops(Easy)>15 and nops(Easy)<=25 then numeqn:=8*par[7]; elif nops(Easy)>25 and nops(Easy)<=35 then numeqn:=6*par[7]; elif nops(Easy)>35 and nops(Easy)<=40 then numeqn:=5*par[7]; elif nops(Easy)>40 and nops(Easy)<=45 then numeqn:=4*par[7]; elif nops(Easy)>45 then numeqn:=3*par[7]; fi; $$$$$$ end of section of sfprogs code $$$$$$ More extensive changes to the way equations are moved between lists require great care (note: par[8] in sfprogs corresponds to p2); but is an activity which warrants investigation and since it is central to efficiency of the program. The p2 knob ----------- Our second parameter controls the flow of equations in two ways, and is sensitive to the length of equations in UnC. The length of an equation is given by Maple's length function and is a simple measure of the complexity of an equation by how much memory is required to store it. The equations in UnC are placed in order according to their increasing length. If the length of the first equation in UnC is more than p2 times longer than the length of the longest equation in Ez, then that equation is the only equation added to Ez in that iteration. If not, then equations of UnC are considered in order of increasing length until either we have the maximum number of equations (as dictated by the p1 knob) or the length of the next equation is more than p2 times longer than the length of the first equation added to Ez. Setting p2=10000 will make standard form ignore the lengths of equations in UnC, and transfer the maximum number allowed by the p1 knob. This option is very useful for small complicated systems as the simplification of complicated equations is put off until after the simplification of the easier equations. This gives the easier equations an opportunity to simplify the complicated (UnC) equations before they are put in Ez. A note of clarification. In the process of moving equations from the UnC list to the Ez list, standard_form will automatically move equations that are not polynomials in their top derivatives into the NL list, (example, sin(V(1,1,1))+V(1,1,1)=0), and these will not be counted towards the maximum number of equations that can be transferred from the UnC list, as determined by the p1 knob. However, equations which are polynomial in their top derivative will be transferred to the Ez list, even if their degree is higher than 1, and all such equations will be counted in the total number of such equations to be allowed. However, after reducing these equations, with respect to the other equations in Ez, all those which are still polynomially nonlinear will be moved from the Ez to the NL list, and standard_form will not transfer any more equations from UnC to Ez, even if their are still some which are short enough given the settings of the p2 knob. Example 1: strategy=[1.75, 1, 2] p1=1, p2=2 (Default knob settings and infolevel set to its highest value) Current Max length in Ez: 117 Lengths in UnC (min to max): 110, 140, 151, 152, 240, ...... Current Number of equations in Ez: 7 In this case there are 7 equations in the Ez list at the beginning of the iteration, so, by the above table describing the operation of the p1 knob, the maximum number of equations that will be added to the Ez list is 10 x p1 = 10 x = 10. But the fifth equation in UnC has length 240 > 2 x 110 = p2 x 110, and 110 is the length of the first equation in UnC, therefore the first equation added to Ez, so due to the second part of the p2 test only the first 4 equations, with lengths 110, 140, 151, 152, which are all less than 2 x 110, are moved. Example 2: strategy=[1.75, 0.5, 2] p1=0.5, p2=2 Current max length in Ez: 164 Lengths in UnC(min to max): 1310, 1476, .... Current number of equations in Ez: 20 In this case there are 20 equations in Ez initially, so the maximum number of equations that will be added is 8 x p1 = 8 x 0.5 = 4. But the first equation in UnC has length 1310 > 2 x 164 > p2 x 164, so since 164 is the longest equation initially in Ez, only the first equation is moved, by the first part of the p2 test. Example 3: strategy=[1.75, 0.5, 2] p1=1, p2=2 Current max length in Ez: 54 Lengths in UnC(min to max): 58, 61, 62, 66, 75, 79, 80, 81, ...... Current number of equations in Ez: 43 In this case the maximum number of equations is 4 x p1 = 4 x 1 = 4, since their are 43 equations in Ez initially. Since all of the first four equations in Unc have lengths less than 2 x 54 = p2 x 54, all of them are chosen. Special Settings: Default Knob Settings: p1=1, p2=2 Ultra-Conservative (one equation at a time): p1=1, p2=0.9 With this setting the program will only move one equation at a time into the Ez list, since the second equation is always greater than 0.9 times the length of the first. For example including strategy=[1.75, 1, 0.9] in the option list for standard_form will set the infolevel at 1.75, and the knob settings to the ultra-conservative settings of p1=1, p2=0.9. Ultra-HOG (everything at once): p1=100000, p2=100000 Note: only use Ultra-HOG on very simple systems. 3.3.3 THE STRATEGY OPTIONS INCONSOPT AND SIMPOPT The 4th and 5th elements of the strategy list must be boolean valued, and their use is outlined in the table below; they have rarely if ever been used. strategy list element value result 4th (inconsopt) true standard_form will terminate if any inconsistent equations are found. false standard_form will save the inconsistent equation in the global variable _incons, and continue 5th (simpopt) true standard_form will collect coefficents of dependent variables prior to any call of maple`s simplify procedure. false standard_form will apply maple`s simplify procedure directly on the equations, without first collecting coefficients of the dependent variables The default value for both of these parameters is true. INCONSOPT: The 4th element of strategy is the inconsopt parameter. When its value is true, the standard_form will terminate as soon as an equation relating only the independent variables, as illustrated in example C of section 2.2.3. When the inconsopt parameter has the value false, then any such equations that standard_form finds will be stored in the global variable _incons. Standard_form will then ignore the equation, and continue to perform the algorithm on the remaining equations. This allows the standard form to be calculated when relations amongst the independent variables are found that actually reduce to zero, but this fact is not recognized by maple. For instance, the equation sin(x1)^2+cos(x1)^2-1=0 reduces to zero, but maple will not recognize this. The standard form obtained with inconsopt set to false is only valid if all the equations in incons reduce to zero. SIMPOPT: The 5th element of strategy is the simpopt parameter. If it has its default value of true, then whenever maple simplifies equations, it will first collect the coefficients of any expressions involving the dependent variables and their derivatives, and then will apply maple`s simplify procedure. If the simpopt parameter is set to false, then standard_form will apply maple's simplify procedure on the equations directly, without collecting coefficients first. 3.4 STORING INTERMEDIATE RESULTS AND THE USEREASY OPTION 3.4.1 STORING INTERMEDIATE RESULTS If the option store is included in the argument list for standard_form, e.g. > sf:= standard_form(sys, store); it will store the status of the calculations at the end of each iteration to a file called StandardStorage.m in your current directory. This file will be in stored in maple's internal format (as indicated by the .m suffix), so the file cannot be viewed as stored, but must be read into maple. To access or view the information in StandardStorage.m, go through the following steps: 1) Start a new maple session. If you want to view the information in StandardStorage.m as it is updated, that is, while standard_form is running, this will mean starting a new maple session independent of the one that standard_form is running on, by, for example, starting a new system shell. Once this new session is started, the Standardform package will be needed to read the file, so type > read sfprogs; 2) Now call the procedure readstored(). This returns the state of the calculations at the end of the last iteration (integrability condition loop) in the form of a list. For example, if we call the list L: > L:= readstored(); then the list L will have 4 elements as follows: L[1] is the list of one term equations L[2] is the list of easy equations L[3] is the list of equations nonlinear in their highest derivatives L[4] is the list of unclassified equations L[5] is the list of pivots as they were at the end of the last iteration. The file that the intermediate results are stored in can be given a name different from StandardStorage, if so desired, if, for example, the results from several runs are needed. To store the results in a file named filename, type, when standard_form is invoked > sf:=standard_form(sys,store=`filename`); The results will actually be stored in a file filename.m, in maple's internal format, and must be be read into a maple session to view the most updated results: > L:=readstored(`filename`); (note that the .m is not be included). The command > standard_form(sys,storeall); will save the results from each iteration in a different file, and thus, to have a permanent record of all of them (if you're not concerned with flooding you file-space). In particular the results at the end of the nth iteration will be stored in the file StandardStorage_n.m. To view the results of the nth iteration type > L:=readstored(`StandardStorage_n`); It is also possible to store results to files with user-assigned names: > standard_form(sys,storeall=`filename`); The results at the end of the nth iteration will then be stored in a file filename_n.m. To view these files, start a maple session, read in the standard form package, and type > L:=readstored(`filename_n`); This is useful for: recovering intermediate results if there is a memory explosion, and restarting them, (perhaps with the original system and without some of the unclassified equations), and also peeking at the current status of the system. 3.4.2 SELECTING NICE EQUATION SUBSETS WITH USEREASY Including: usereasy = N (N a positive integer) in the option list for standard_form, will force standard_form to take the first N equations in the input system to standard_form into the easy list Ez at standard_form's first iteration. This over rides standard_form's Ez/UnC choice strategy at its first iteration. It is most useful when the user has their own premonition of a set of equations that would help the program out (and places these as the first N equations in the input system list). One situaiton in which this can be helpful is when memory explosion has occurred, and the results of previous run before the explosion have been saved using one of the storage options. The user chooses N easy equations, throws out the ugly ones, and forms a new system, including the original system: > sf:= standard_form(sys, store); ..... ..... !!!! memory explosion !!!! pick N easy eqns from stored eqns, throwing out the ugly ones Advisable to start a new session if the program went down in flames > sysnew:= [N easy, op(sys)]; > sf:= standard_form(sysnew, usereasy = N); 3.5 CLASSIFICATION PROBLEMS - SYSTEMS OF PDE WITH UNSPECIFIED FUNCTIONS OR PARAMETERS Consider a system of PDEs with dependent variables V(1),...,V(n) which are functions of independent variables x1,...,xm. In addition suppose that the system involves certain parameters or functions of unspecified form ("classification parameters or functions"). In general the properties and number of solutions of the system will depend on the values and forms of these parameters and functions. We will call the "classification problem for the system" the problem of determining the standard form of the system for every functional form of the classification functions. In general a "classification tree" results, each branch of which is characterised by a certain number of differential equalities and inequalities amongst the classification functions and parameters. 3.5.1 LOW-TECH APPROACH TO CLASSIFICATION Problems where there appear only parameters or say only one unspecified function can often be handled by the standard form procedure without embellishment. Already we have classified the system of algebraic equations sysA, given in Example A of section 2.2. The resulting classification tree was: sysA | k=0 | k#0 ----------------------------------------------- | | | | sfA sfA1 idA idA1 where the classification parameter was k and separate branches were picked out by examining the pivots list _pivs. We also classified the system of PDEs sysC given in Example C of Section 2.3 with classification parameter r which yielded the classification tree: sysC | r=0 | r#0 ----------------------------------------------- | | | | sfC2 inconsistent idC2 no solutions Example A. We will classify the system of PDEs: g(t) Ux - g'(t) U = 0 , Ut = 0 where U = U(x,t), and g(t) is the classification function. Then in our notation: > read sfprogs; > sysA:= [ g(x2)*V(1,1) - diff(g(x2),x2)*V(1) = 0, V(1,2) = 0]; Case 1: No constraints on the classification function g(t) - the generic case. > sfA1:= standard_form(sysA); yields sfA1 := [ V(1) = 0] and initial data > idA1:= initial_data(sfA1); idA1:= [[],[]] The pivots are > _pivs; / 2 \ | d | / d \2 |------ g(x2)| g(x2) - |----- g(x2)| | 2 | \ dx2 / \ dx2 / [-1 = 0, g(x2) = 0, ------------------------------------- = 0] 2 g(x2) and these can be given a more transparent form by using the sfprogs procedure shortenpiv which removes repeated factors from the pivots: > spivsA1:= shortenpiv(_pivs); / 2 \ | d | / d \2 spivsA1 := [g(x2), |------ g(x2)| g(x2) - |----- g(x2)| ] | 2 | \ dx2 / \ dx2 / So g(x2) = 0 and 2 g''(x2) = g'(x2) /g(x2) are special cases. Case 2: g(x2) = 0 > g:= proc(x2) 0 end; > sfA2:= standard_form(sysA); yields sfA2 := [V(1,2) = 0] with initial data > idA2 := initial_data(sfA2); idA2 := [[], [V(1) = f1(x1)]] Case 3: g''(x2) = g'(x2)^2/g(x2) To invoke the above simplification rule we can use Maple's inbuilt capability for defining how diff acts on a function. We employ Maple's diff/ capability in two steps: > `diff/g`:= proc(x2,x2) h(x2) end; # i.e. g'(x2) = h(x2) > `diff/h`:= proc(x2,x2) h(x2)^2/g(x2) end; # i.e. h'(x2) = h(x2)^2 /g(x2) > sfA3:= standard_form(sysA); yields sfA3 := [ V(1,1) = h(x2) V(1)/g(x2), V(1,2) = 0 ] with initial data > idA3:= initial_data(sfA3); idA3 := [[V(1)=a1], []] and (shortened) pivots > spivsA3:= shortenpiv(_pivs); spivsA3:= [g(x2)] This is just Case 2 again, so there are no further branches. 3.5.2 HI-TECH APPROACH TO CLASSIFICATION - the class_eqns option General classification problems involve the classification of systems of PDEs in dependent variables V(1), ... , V(n) and independent variables x1, ... , xm and classification functions g1, ... , gs. In the hi-tech classification approach the classification functions g1, ... , gs are regarded as additional dependent variables ("classification variables") V(n+1), ... , V(n+s) which satisfy an consistent system of PDEs called the classifying equations. It is important to note that: THE CLASSIFYING EQUATIONS MUST BE IN STANDARD FORM. THE DEPENDENCIES OF THE CLASSIFYING VARIABLES MUST BE DECLARED. THE CLASSIFICATION VARIABLES MUST BE NUMBERED HIGHER THAN THE OTHER DEPENDENT VARIABLES In particular the classifying equations are entered as an option: > read sfprogs; > standard_form(sys, class_eqns = class) Here class is the list of the classifying equations in standard form and the dependencies of each classifying variable V(r) declared in the form of the list: [ V(r), k,l,...] meaning that V(r) depends on xk, xl, etc. Example A. Consider again the system of PDEs dealt with by the low-tech approach in Example A of Section 3.5.1: g(t) Ux - g'(t) U = 0 , Ut = 0 where U = U(x,t), and g(t) is the classification function. In the hi-tech approach we set U=V(1) and define the classification variable V(2)=g(t) (the classification variables must be numbered higher than the nonclassification variables). So in Maple > read sfprogs; > sysA:= [ V(2)*V(1,1) - V(2,2)*V(1) = 0, V(1,2) = 0 ]: Case 1: No constraints on the classification variable V(2) - the generic case Then the class equations are entered as a list: > classA1:= [[V(2),2]]; At the moment there are no differential constraints on V(2), so the only information we enter is the dependency info that V(2) only depends on x2. A valid alternative here would have been classA2:=[[V(2),2],V(2,1)=0], but not classA2:=[V(2,1)=0] since standard_form uses the dependency declaration to determine which are the classification variables. Reduction of the system to standard form gives: > sfA1:= standard_form(sysA, class_eqns = classA1); The system parameters are as follows: There are 2 independent variables: [x1, x2] There are 1 dependent variables: [V(1)] There are 1 classification variables: [V(2)] sfA1 := [V(1) = 0] which agrees with Case 1 of Example A in 3.5.1. Notice that, standard_form has reassured us by alerting us that V(2) has been declared as a classification variable. To calculate the initial data, we apply initial_data with the option class_eqns = classA1, alerting initial_data that V(2) is a classification variable: > idA1:= initial_data(sfA1, class_eqns=classA1); idA1 := [[], [V(2) = f1(x2)]] Note that the idA1 includes the initial data for the classification variable. The pivots are > _pivs; 2 V(2, 2, 2) V(2, 2) [-1 = 0, V(2) = 0, ---------- - -------- = 0] V(2) 2 V(2) and in their shortened form: > spivsA1:= shortenpiv(_pivs); 2 spivsA1 := [V(2, 2, 2) V(2) - V(2, 2) , V(2)] So V(2) = 0 and 2 V(2, 2, 2) V(2) - V(2, 2) = 0 are special cases. Case 2: V(2) = 0 We need to give both the dependency info and the class equations in standard form in the list class_eqns. V(2) as a constant, does not depend on anything, and to reflect this [V(2)] is included in class_eqns. Obviously the standard form of the class equations is V(2) = 0. Reduction of this case to standard form gives: > sfA2:= standard_form(sysA, class_eqns = [[V(2)],V(2)=0]); sfA2 := [V(1, 2) = 0] and > spivsA2:= shortenpiv(_pivs); spivsA2 := [] so there are no nontrivial shortened pivots. The initial data is: > idA2:= initial_data(sfA2, class_eqns= [[V(2)], V(2)=0]); idA2 := [[], [V(1) = f1(x1)]] To alert initial_data is a classification variable, we have inserted the class_eqns option in the same way as for standard_form. 2 Case 3: V(2, 2, 2) V(2) - V(2, 2) = 0 This case corresponds to the vanishing of the first pivot in spivsA1: > class_sys:= [spivsA1[1]=0]; 2 class_sys := [V(2, 2, 2) V(2) - V(2, 2) = 0] To reduce the classifying system to standard form we need the option divpiv so that standard_form will divide by the pivot V(2): > sfclass_sys:= standard_form(class_sys, divpiv); 2 V(2, 2) sfclass_sys := [V(2, 2, 2) = --------] V(2) To give the classifying equations in proper format we include the dependency info for V(2): > classA3:= [[V(2),2], op(sfclass_sys)]; 2 V(2, 2) classA3 := [[V(2), 2], V(2, 2, 2) = --------] V(2) Reduction of the original system to standard form for this case gives: > sfA3:= standard_form(sysA, class_eqns = classA3); V(2, 2) V(1) sfA3 := [V(1, 1) = ------------, V(1, 2) = 0] V(2) > spivsA3:= shortenpiv(_pivs); spivsA3 := [V(2)] so that no further branching occurs. > idA3:= initial_data(sfA3, class_eqns = classA3); idA3 := [[V(1) = a1, V(2) = a2, V(2, 2) = a3], []] Note that the initial data includes the initial data for the classification equation for V(2). In summary we have the classification tree: sysA g(t) Ux - g'(t) U = 0 Ut = 0 | | g(t)=0 | g(t)#0 ---------------------------------------- | | | | | Ut = 0 | | U(x,t0) = f1(x) | | Case 2 | | | 2 | 2 g''(t)=g'(t) /g(t) | g''(t)#g'(t) /g(t) ------------------------------------------- | | | | Ut = 0 U = 0 Ux = g'(t)U/g(t) Case 1 U(x0,t0) = a1 Case 3 Figure: Classification tree for the system of PDEs g(t) Ux - g'(t) U = 0, Ut = 0 for U=U(x,t) with classification function g(t). Each branch gives the pivot conditions, initial data and standard form of the system. 3.6 INTERFACE TO MAPLE'S DERIVATIVE NOTATION To convert a system of PDEs written in the repeated index notation specific to our programs we use the convprogs routine rep2map: Example A To change the following system which is in our repeated index notation to Maple's usual notation for derivative: > sysA := [V(2,3) = 0, V(2,1) = 0, V(1,3) = 0, x1*x2*V(3,3,3) = 0, V(2,2)-2*V(1,1) = 0, 2*V(3,1,3)+V(1,2)-V(1,1,1) = 0, V(3,2)-V(3,1,1) = 0]; sysA := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, x1 x2 V(3, 3, 3) = 0, V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 1, 3) + V(1, 2) - V(1, 1, 1) = 0, V(3, 2) - V(3, 1, 1) = 0] we read the conversion package (sfprogs is not needed): > read convprogs; To use the conversion function rep2map (repeated index to maple) which takes as its first argument the system of PDEs written as a list of equations in our repeated-index notation, and as second argument the number of independent variables (=3 here): > msys:= rep2map(sysA, 3); which yields: msys := [D[3](g)(x1, x2, x3) = 0, D[1](g)(x1, x2, x3) = 0, D[3](f)(x1, x2, x3) = 0, x1 x2 D[3, 3](h)(x1, x2, x3) = 0, D[2](g)(x1, x2, x3) - 2 D[1](f)(x1, x2, x3) = 0, 2 D[1, 3](h)(x1, x2, x3) + D[2](f)(x1, x2, x3) - D[1, 1](f)(x1, x2, x3) = 0, D[2](h)(x1, x2, x3) - D[1, 1](h)(x1, x2, x3) = 0] This function is useful for interfacing our package sfprogs with other packages and procedures e.g. Maple's dsolve command. A procedure, map2rep is also available which does the opposite conversion. If you have a system written in maple's standard derivative notation, in some dependent and independent variables, with arbitrary names, map2rep will convert it to a system in repeated index notation. The syntax for the map2rep procedure is > map2rep(expr, deps, indeps) where expr is the expression to be converted, and deps and indeps are lists of the names used in expr for the dependent and independent variables, in the order that they are to be labelled in the repeated index notation. In the situation where the independent variables are already written as xi, i=1..n, just include a list of independents consisting of the xi. The independent variables appearing explicitly in the expression will then be unaffected by the conversion procedure. As an example, we do the reverse of the above example. To show the most general situation, we enter the system with different names for the independents: we make the substitutions x1=x,x2=y,x3=z. The system, in maple notation, is then sysC := [D[3](g)(x, y, z) = 0, D[1](g)(x, y, z) = 0, D[3](f)(x, y, z) = 0, x y D[3, 3](h)(x, y, z) = 0, D[2](g)(x, y, z) - 2 D[1](f)(x, y, z) = 0, 2 D[1, 3](h)(x, y, z) + D[2](f)(x, y, z) - D[1, 1](f)(x, y, z) = 0, D[2](h)(x, y, z) - D[1, 1](h)(x, y, z) = 0] To convert it to repeated index notation, we use the program map2rep as follows: > map2rep(sysC,[f,g,h],[x,y,z]); [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, x1 x2 V(3, 3, 3) = 0, V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 1, 3) + V(1, 2) - V(1, 1, 1) = 0, V(3, 2) - V(3, 1, 1) = 0] and this is the system we started with originally. 3.7 INTERFACE TO MANSFIELD'S NOTATION FOR HER DIFFGROB PACKAGE As we have already mentioned, the program standard_form is basically a linear solver, but with some capabilities for nonlinear systems of PDEs. An approach and a Maple program for simplifying polynomially nonlinear systems of PDEs has been given by Mansfield [Mans91]. We illustrate by example how to convert back and forth between Mansfield's notation (recently she changed her notation, but we have not yet changed our conversion program to work for this changed notation) and our own: Example A To change the following system which is in our repeated index notation to the notation used by Mansfield's program: > sysA := [V(2,3) = 0, V(2,1) = 0, V(1,3) = 0, x1*x2*V(3,3,3) = 0, V(2,2)-2*V(1,1) = 0, 2*V(3,3,1)+V(1,2)-V(1,1,1) = 0, V(3,2)-V(3,1,1) = 0]; sysA := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, x1 x2 V(3, 3, 3) = 0, V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 3, 1) + V(1, 2) - V(1, 1, 1) = 0, V(3, 2) - V(3, 1, 1) = 0] we read the conversion package: > read convprogs; To use the conversion function rep2man (repeated index to mansfield) we first have to declare the global variables vars: (the independent variables) and ukns (the names of the dependent variables). > vars:= [x1,x2,x3]: > ukns:= [v1,v2,v3]: The system in Mansfield's notation is then: > msysA:= rep2man(sysA); yielding msysA := [p(v2, [x1, x2, x3], [0, 0, 1]), p(v2, [x1, x2, x3], [1, 0, 0]), p(v1, [x1, x2, x3], [0, 0, 1]), x1 x2 p(v3, [x1, x2, x3], [0, 0, 2]), p(v2, [x1, x2, x3], [0, 1, 0]) - 2 p(v1, [x1, x2, x3], [1, 0, 0]), 2 p(v3, [x1, x2, x3], [1, 0, 1]) + p(v1, [x1, x2, x3], [0, 1, 0]) - p(v1, [x1, x2, x3], [2, 0, 0]), p(v3, [x1, x2, x3], [0, 1, 0]) - p(v3, [x1, x2, x3], [2, 0, 0])] To convert a system written in Mansfield's notation to repeated index notation we use the function man2rep (again both vars - the names of the independent variables in the Mansfield system (these may be different to x1,x2,...) and ukns (the names of the dependent variables): > vars:= [x1,x2,x3]: > ukns:= [v1,v2,v3]: > rsysA := man2rep(msysA); yielding rsysA := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, x1 x2 V(3, 3, 3) = 0, V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 1, 3) + V(1, 2) - V(1, 1, 1) = 0, V(3, 2) - V(3, 1, 1) = 0] 3.8 THE DIRECTED GRAPH OF A PDE SYSTEM This procedure creates the directed graph of a PDE system as described in [Reid91c]. This procedure is still preliminary and has not been thoroughly checked. But have fun anyway! Example: xmaple (Maple V release 2) is needed for this application > read sfprogs; > read graphDE; # read the graph procedure which employs subroutines # from sfprogs > with(networks); > read `heat.sys`: # the determining system for the Heat equation > sf:= standard_form(sys, janet_order=[[0,0,1,0,0,0]]); # the standard form and the graph will depend on # the ordering which is used > msf:= make_vset(sf); #returns the vertices and edges of the graph > G:= graph(msf): # makes the data structure for the graph > show(G); # shows the graph as a data structure > draw(G); # attempts to draw the graph is directed - # but xmaple (Maple V release 2) does # not yet draw arrows on the edges 4. APPLICATION OF THE PACKAGE TO SYMMETRY ANALYSIS OF PDES A symmetry of a system of PDEs is a transformation which leaves invariant the family of solutions of the system. Symmetries have proved to be a useful tool for deriving analytical and geometrical properties of systems of PDE: invariant solutions, mappings between equations, connections with conservation laws (Noether's Theorem), and the relationship of higher order symmetries (Lie-Backlund symmetries) with the inverse-scattering method. Lie point symmetries are the simplest type of symmetries of systems of PDEs. Given a system of PDEs with independent variables x=(x1,x2,...,xp) and dependent variables u=(u1, u2, ... , uq) a 1 parameter (a)-Lie symmetry group is a group of transformations of the form: x' = X(x,u; a) u' = U(x,u; a) i.e. the transformation parameterised by a, maps a point (x,u) to a point (x', u') in the space of independent and dependent variables of the system of PDEs. Instead of dealing with the group transformations in their (nonlinear) form above Lie dealt with the transformations in their (linearized) infinitesimal form: 2 x' = X(x,u; a) = x + V(x,u) a + O(a ) 2 u' = U(x,u; a) = u + W(x,u) a + O(a ) with associated infinitesimal generator L = V(1)d/dx1 + V(2)d/dx2 + ... + V(p)d/dxp + W(1)d/du1 + ... + W(q)d/duq in the infinitesimals V(1),...,V(p),W(1),...,W(q). In order to input such infinitesimals into our programs we'll replace uk, W(k) by x.(p+k), V(p+k) (x.1 = x1, and u1, W(1) is replaced by x.(p+1), V(p+1) etc. ). In this notation the infinitesimal operator is L = V(1) d/dx1 + ... + V(r) d/dxr where r = p+q. According to a fundamental theorem of Lie the infinitesimals V(1),...,V(r) uniquely characterise the group. Lie provided a simple algorithm [BlK89], [Olv86], [Ovs82] for deriving the equations which must be satisfied by the infinitesimals. Today there are many effective symbolic programs for deriving these so-called "defining" or "determining equations" [Here93]. The determining equations of a given system of PDEs form a linear homogeneous system of PDEs (regardless of whether the original system was linear or not). This linearization of the group problem was a key to the success of Lie's methods. In addition such determining systems are usually very overdetermined. Another fundamental result of Lie is that if M and N are two infinitesimal symmetry generators of a PDE system then so too is the commutator of M and N: [ M, N ] = M N - N M So these infinitesimal generators form a Lie algebra with respect to [,]. If this Lie algebra is finite-dimensional then according to Lie it is characterised by its "structure constants" [Ovs82, Olv86, BlK89] k [L(i), L(j)] = C L(k) (sum on k) ij where L(i) is the infinitesimal basis generator corresponding to the i-th basis solution of the determining equations. (i.e. L = a1 L(1) + a2 L(2) + ... + aK L(K) where a1, a2, ... ,aK are the parameters in the general solution of the determining system). The structure of this algebra provides important coordinate-independent (geometric) information about the original system of PDEs: information about reductions to ODEs, mappings to other more tractable PDEs etc. 4.1 CALCULATION OF THE STRUCTURE OF LIE SYMMETRY ALGEBRAS OF PDES FROM THEIR INFINITESIMAL DETERMINING SYSTEMS The conventional approach to infinitesimal symmetries has been to try to integrate the determining equations, explicitly determining the basis generators L(i), and then substituting the results into the commutation relations above to determine the structure constants. Since this process is based on the heuristic process of integration it is not always guaranteed to succeed. Our approach to finding the structure of Lie symmetry algebras of PDEs is completely algorithmic: (1) The determining system for the infinitesimals is either produced by hand or by one of the many computer packages automating this process [Here93, Champ91]. This system is converted into the xi,V(j) notation used by our sfprogs package (see 4.2,4.3,4.4). (2) Algorithm standard_form is used to reduce the determining system to standard form. (3) Application of initial_data to the standard form produced in Step (2) algorithmically determines the dimension of the algebra of infinitesimal symmetries. (4) Algorithm fin_com_table uses the standard form from Step (2) and its initial data from Step (3) to determine the structure of the of the finite-dimensional subalgebra of infinitesimal symmetries corresponding to the parameters a1,a2,... in the finite initial data. This algorithm depends only on differentiation and substitution. (5) taylor series approximations of the infinitesimals or basis generators are found to any chosen finite degree using algorithms taylor_exp, or finite_generators respectively. Notes 1. Commutators of generators belonging to infinite-dimensional Lie symmetry algebras can also be computed using the procedure lie_symm (instead of fin_com_table). However the output of lie_symm is complicated and currently does not have a good structural interpretation. 2. Classes of PDEs with variable coefficients can be treated using the class_eqns option (see Section 4.4.2). 3. In many applications the explicit forms of the infinitesimals are required (e.g. in the calculation of invariant solutions) so our algorithms do not make redundant the packages which attempt to explicitly solve determining equations. The role of our algorithms is instead a complementary one and we are interested in developing interfaces with such packages. 4.1.1 SYMMETRIES OF THE HEAT EQUATION: Ut = Uxx We will not violate a longstanding tradition so we first consider the heat equation: Ut - Uxx = 0 where U=U(x,t). Step (1): In Section 4.2.1 it is shown how to produce the determining system for the heat equation using the macsyma package symgrp.max. The converted system has been stored in the file called heat.sys. > read sfprogs; We next read the file heat.sys containing the infinitesimal determining system for the Heat equation > read `heat.sys`; sys := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, V(3, 3, 3) = 0, V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 3, 1) + V(1, 2) - V(1, 1, 1) = 0, V(3, 2) - V(3, 1, 1) = 0] The infinitesimal generator of the algebra is given by L = V(1) d/dx1 + V(2) d/dx2 + V(3) d/dx3, where x1 = x (space), x2 = t (time) and x3 = U (temperature), with corresponding infinitesimals V(1), V(2) and V(3). Step (2): We reduce the system to standard form: > sf:= standard_form(sys); sf := [V(3, 2, 2, 3) = 0, V(3, 1, 1) = V(3, 2), V(3, 1, 3) = - 1/2 V(1, 2), V(1, 2, 2) = 0, V(2, 2, 2) = - 4 V(3, 2, 3), V(3, 3, 3) = 0, V(1, 1) = 1/2 V(2, 2), V(2, 1) = 0, V(1, 3) = 0, V(2, 3) = 0] (For fun also try: > sffun:= standard_form(sys, janet_order=[[0,1,0,0,0,0]]); > idfun:= initial_data(sffun); ) Step (3): We calculate the initial data of the system: > id:= initial_data(sf); id := [ [ V(1)=a1, V(1, 2)=a2, V(2)=a3, V(2, 2)=a4, V(3, 3)=a5, V(3, 2, 3)=a6 ], [V(3)=f1(x2), V(3, 1)=f2(x2)]] See Section 3.2 for interpretation of the initial data. So the dimension of the group of Lie symmetries of the system is infinite-dimensional. It is determined by 6 parameters and 2 arbitrary functions of 1 variable. Note that we did not need to integrate the determining system to uncover this information. Step (4): To determine the commutation relations amongst the generators of the finite (6-dimensional) component of the algebra we use the procedure fin_com_table: > ct:= fin_com_table(sf, id); [ 0 - 1/2 L(5) 0 1/2 L(1) 0 - 2 L(2) ] [ ] [ * 0 - L(1) - 1/2 L(2) 0 0 ] [ ] [ * * 0 L(3) 0 - 4 L(4) + L(5) ] ct := [ ] [ * * * 0 0 L(6) ] [ ] [ * * * * 0 0 ] [ ] [ * * * * * 0 ] The i-j entry of this table corresponds to [L(i), L(j)] (e.g. [L(1), L(2)] = -1/2 L(5)). To look at ct later use op(ct) or print(ct). Step (5): To obtain the taylor series approximations of the infinitesimals to degree 3 about (X1,X2,X3) in powers of h1,h2,h3: > ts:= taylor_exp(sf, id, 3); ts := [V(1) = a1 + 1/2 a4 h1 + a2 h2 - 2 a6 h1 h2, 2 V(2) = a3 + a4 h2 - 2 a6 h2 , V(3) = / d \ / d \ 2 f1(X2) + f2(X2) h1 + |----- f1(X2)| h2 + a5 h3 + 1/2 |----- f1(X2)| h1 \ dX2 / \ dX2 / / 2 \ / d \ | d | 2 + |----- f2(X2)| h1 h2 - 1/2 a2 h1 h3 + 1/2 |------ f1(X2)| h2 + a6 h2 h3 \ dX2 / | 2 | \ dX2 / / 2 \ / d \ 3 | d | 2 2 + 1/6 |----- f2(X2)| h1 + 1/2 |------ f1(X2)| h1 h2 + 1/2 a6 h1 h3 \ dX2 / | 2 | \ dX2 / / 2 \ / 3 \ | d | 2 | d | 3 + 1/2 |------ f2(X2)| h1 h2 + 1/6 |------ f1(X2)| h2 | 2 | | 3 | \ dX2 / \ dX2 / ] By linearity we can express the general infinitesimal generator in terms of basis generators L(1),...,L(6) corresponding to the initial data a1,a2,...,a6 respectively and two basis generators M1(f1(x2)), M2(f2(x2)) corresponding to the arbitrary functions f1,f2 in the initial data. That is: L = a1 L(1) + a1 L(2) + ... + a6 L(6) + M1(f1(x2)) + M2(f2(x2)) In particular we can get taylor approximations of L(1),...,L(6) by making substitutions in the taylor series above: For example: L(1) = T[1]*Dx(1) + T[2]*Dx(2) + T[3]*Dx(3) where T[i] = coefficient a.i in ts[i] and Dx(1) = d/dx1 etc. We've provided a procedure finite_generators, which does this automatically. For example to obtain the taylor approximations to degree 4 of the finite (subgroup) basis generators: > fgs:= finite_generators(sf, id, 4); fgs := [L(1) = Dx(1), L(2) = h2 Dx(1) - 1/2 h3 h1 Dx(3), L(3) = Dx(2), L(4) = 1/2 h1 Dx(1) + h2 Dx(2), L(5) = h3 Dx(3), 2 2 L(6) = - 2 h1 h2 Dx(1) - 2 h2 Dx(2) + (h3 h2 + 1/2 h3 h1 ) Dx(3)] where Dx(1) represents d/dx1 etc. In our original notation expanding about (X1,X2,X3)=(0,0,0) these generators are: L(1) = d/dx + O(5), L(2) = t d/dx - 1/2 x U d/dU + O(5), L(3) = d/dt + O(5), L(4) = 1/2 x d/dx + t d/dt + O(5) L(5) = U d/dU + O(5), L(6) = -2xt d/dx - 2 t^2 d/dt +(tU + x^2 U/2)d/dU + O(5) Similarly, although we have not automated it yet, we can obtain taylor approximations of the generators M(f1(x2)) and M2(f2(x2)) corresponding to the infinite initial data. For the moment you can do this yourself for M1(f1(x2)) by assigning f2:= proc(z) 0 end: and substituting a1=0,a2=0,...,a6=0 into ts. And similarly for M2(f2(x2)). As is well known the taylor approximations of the 6-dimensional finite symmetry subalgebra of the Heat equation are at most cubic in x,t,U, so the remainder terms O(5) in the above generators are exact. There are no 4-th order terms in the above expansions and this motivates the general question: If there are no N-th order terms in the taylor series approximations of the generators does this mean that they are exact to order N-1 ? (Answer - Yes: see Section 5.1). We emphasise, however, that the termination of the Taylor series had nothing to do with the effectiveness of the algorithm fin_com_table for calculating the structure of the finite dimensional subalgebra above. In particular, fin_com_table does not even use Taylor series in the construction of structure constants (instead it uses the approach of [Reid92a] which makes the Taylor series approach of [Reid91b] obsolete). Finally the command > lie_symm(sf, id); will display the entire 8 X 8 commutator table, between the operators L(1),L(2),L(3),L(4),L(5),L(6),M1(f1(x2)),M2(f2(x2)), in a somewhat obscure vector-format (how about that for great documentation!). In particular the algorithms can calculate the commutators for infinite-dimensional Lie symmetry groups in terms of their initial data [Reid90, Reid91b, Reid92b] (again no knowledge of the solutions is required). However we have not as yet developed a good structural (coordinate-independent) interpretation of this table. Despite its name lie_symm is not a part of or an application of the differential forms package liesymm. 4.1.2 A GROUP CLASSIFICATION PROBLEM: Uxx = Utt + H(U) Ut The group classification problem for the family of nonlinear telegraph equations Uxx = Utt + H(U) Ut where U=U(x,t) (*) is to determine the size and structure of the group of Lie symmetries of (*) for every functional form of H(U) (see [Reid91b] for complete results). To illustrate the method We will follow a couple of branches of the classification. Step (1): In Section 4.2.2 it is shown how to produce the determining system for (*) by using the macsyma package symgrp.max. The converted system has been stored in the file called telegraph.sys. > read sfprogs; We next read the file telegraph.sys containing the infinitesimal determining system of (*): > read `telegraph.sys`; sys := [V(1, 3) = 0, V(2, 3) = 0, V(3, 3, 3) = 0, V(2, 1) - V(1, 2) = 0, V(2, 2) - V(1, 1) = 0, V(3, 2) V(4) - V(3, 2, 2) + V(3, 1, 1) = 0, V(1, 2) V(4) - 2 V(3, 3, 1) - V(1, 2, 2) + V(1, 1, 1) = 0, V(3) V(4, 3) - V(2, 2) V(4) + 2 V(1, 1) V(4) + V(2, 2, 2) - V(2, 1, 1) - 2 V(3, 3, 2) = 0 ] The infinitesimal generator of the algebra is given by L = V(1) d/dx1 + V(2) d/dx2 + V(3) d/dx3, where x1 = x (space), x2 = t (time) and x3 = U, with corresponding infinitesimals V(1), V(2) and V(3). The classification variable V(4) depends only on x3. Case 1: No constraints on V(4). Step (2): Since we are treating V(4) as a classification variable (see 3.5.2) we have to enter its dependency info in the list class_eqns. Reducing the system to standard form > sf1:= standard_form(sys, class_eqns=[[V(4),3]]); sf1 := [V(1, 1) = 0, V(2, 1) = 0, V(1, 2) = 0, V(2, 2) = 0, V(1, 3) = 0, V(2, 3) = 0, V(3) = 0] Step (3): Calculation of the initial data of the system gives > id1:= initial_data(sf1, class_eqns=[[V(4),3]]); id1 := [[V(1) = a1, V(2) = a2], [V(4) = f1(x3)]] where we have entered the class_eqns option to initial_data since there is a classification variable (see 3.5.2). The group of Lie symmetries of the system is 2-dimensional. Since we are dealing with a classification problem we look at the (shortened) pivots to identify special cases: > spivs1:= shortenpiv(_pivs); 2 spivs1 := [- 5/2 V(4, 3) + V(4) V(4, 3, 3), 2 2 - 2 V(4) V(4, 3, 3) + V(4) V(4, 3, 3, 3) V(4, 3) + V(4, 3, 3) V(4, 3) , V(4), V(4, 3)] So the above standard form and initial data are valid if these pivots are not identically zero. Step (4): The commutator table is given by: > ct1:= fin_com_table(sf1, id1, class_eqns=[[V(4),3]]); [ 0 0 ] ct1 := [ ] [ * 0 ] Note that again we had to enter the class_eqns option to a procedure since we are dealing with a classification problem. Step (5): Neither taylor_exp nor finite_generators have a class_eqns option yet. So for the time being you have to manipulate your problems into a form with a non-classification form to use these functions and the right initial data. From the pivots spivs1 we see that a special case is: 2 2 Case 2: H(U) H'(U) H'''(U) - 2 H(U) H''(U) + H''(U) H(U) = 0 that is spivs1[2] = 0. The class equations must be in solved form > sfclass2:= standard_form([spivs1[2]=0], divpiv); 2 V(4, 3) V(4, 3, 3) V(4, 3, 3) sfclass2 := [V(4, 3, 3, 3) = - ------------------ + 2 -----------] V(4) V(4, 3) where the option divpiv had to be used to force standard_form to divide by the nonlinear pivot V(4)V(4,3) and solve for V(4,3,3,3). The class equations for Case 2 are: > class2:= [[V(4),3], op(sfclass2)]; 2 V(4, 3) V(4, 3, 3) V(4, 3, 3) class2 := [[V(4), 3], V(4, 3, 3, 3) = - ------------------ + 2 -----------] V(4) V(4, 3) Step (2): Standard forming the system for Case 2 gives > sf2 := standard_form(sys, class_eqns=class2); V(3) V(4, 3) sf2 := [V(1, 1) = - ------------, V(2, 1) = 0, V(3, 1) = 0, V(1, 2) = 0, V(4) V(3) V(4, 3) V(2, 2) = - ------------, V(3, 2) = 0, V(1, 3) = 0, V(2, 3) = 0, V(4) V(3) V(4, 3) V(3) V(4, 3, 3) V(3, 3) = ------------ - ---------------] V(4) V(4, 3) with shortened pivots > spivs2:= shortenpiv(_pivs); spivs2 := [V(4), V(4, 3)] Step (3): The initial data is: > id2:= initial_data(sf2, class_eqns=class2); id2 := [ [V(1)=a1, V(2)=a2, V(3)=a3, V(4)=a4, V(4, 3)=a5, V(4, 3, 3)=a6], [] ] Note that the initial data a4,a5,a6 for second order ODE satisfied by the classification variable V(4) was included in id2. The initial data corresponding to the infinitesimals V(1), V(2), V(3) is a1, a2, a3 so the symmetry algebra for Case 2 is 3-dimensional. Step (4): Calculation of the commutator table gives: > ct2 := fin_com_table(sf2, id2, class_eqns=class2); [ a5 L(1) ] [ 0 0 - ------- ] [ a4 ] [ ] ct2 := [ a5 L(2) ] [ * 0 - ------- ] [ a4 ] [ ] [ * * 0 ] Notice how initial data for the classification variable V(4) has entered into the commutator table. Step (5): As we've already mentioned taylor_exp and finite_generators don't have a class_eqns option so they can't be directly applied to this case. Cunning such as that exhibited below is required to apply taylor_exp: > ts2:= taylor_exp([op(sf2), op(sfclass2), V(4,1)=0, V(4,2)=0], id2, 2); a3 a5 h1 a3 a5 h2 ts2 := [V(1) = a1 - --------, V(2) = a2 - --------, a4 a4 /a3 a5 a3 a6\ 2 V(3) = a3 + |----- - -----| h3, V(4) = a4 + a5 h3 + 1/2 a6 h3 ] \ a4 a5 / Since there are no second order terms in the taylor expansion of the infinitesimals V(1), V(2), V(3) we can assume that these expansions have terminated at first order yielding by inspection of ts2 the exact (see 5.1) symmetries: L(1) = Dx(1), L(2) = Dx(2), L(3) = - (a5/a4) h1 Dx(1) - (a5/a4) h2 Dx(2) + (a5/a4 - a6/a5) h3 Dx(3), where a4 = H(U0), a5 = H'(U0) and a6 = H''(U0). 4.2 INTERFACE TO THE MACSYMA PACKAGE SYMGRP.MAX Requirements: Access to the symgrp.max macsyma symmetry package of Champagne-Hereman-Winternitz (available by ftp) In [Cha90, Here93] a macsyma program symgrp.max is described which produces the determining equations for the infinitesimal symmetries of systems of PDEs. A feedback mechanism is also described for the attempted solution of such systems. We illustrate through examples how the determining systems produced by symgrp.max can be converted into a form suitable for input into our maple program standard_form. For more info contact Willy Hereman at WHEREMAN@flint.mines.colorado.edu 4.2.1 PRODUCING THE DETERMINING SYSTEM OF Ut = Uxx Production of the infinitesimal determining system for the Heat Equation Ut - Uxx = 0 where U=U(x,t). To produce the determining system we batch the macsyma-syntax files heat.dat and maxconv in macsyma, both of which reside in the same directory as sfprogs: (c1) batch("heat.dat"); (c2) batch("maxconv"); heat.dat is typical of the macsyma batch files needed to produce infinitesimal determining systems of PDEs - in this case the heat equation (have a peek at it). You'll notice that heat.dat requires the macsyma program symgrp.max which is available by ftp. We have provided a macsyma program, "maxconv", which takes the determining system produced by the typical macsyma dat files and writes them to a file called lode in your current directory which is in a form more suitable for conversion into maple syntax. After executing the above commands in macsyma, a file called lode should appear in your directory which has contents: Listing of lode: [2, ['DIFF(ETA(2),U(1),1),'DIFF(ETA(2),X(1),1), 'DIFF(ETA(1),U(1),1),'DIFF(PHI(1),U(1),2), 'DIFF(ETA(2),X(2),1)-2*'DIFF(ETA(1),X(1),1), 2*'DIFF(PHI(1),U(1),1,X(1),1)+'DIFF(ETA(1),X(2),1) -'DIFF(ETA(1),X(1),2), 'DIFF(PHI(1),X(2),1)-'DIFF(PHI(1),X(1),2)]]; Before lode can be read into maple it must be stripped of ' using your editor (maybe this could be avoided by extracting lode in a smarter way from the macsyma program, but no one has figured out how to do this yet). After stripping, lode is: [2, [DIFF(ETA(2),U(1),1),DIFF(ETA(2),X(1),1), DIFF(ETA(1),U(1),1),DIFF(PHI(1),U(1),2), DIFF(ETA(2),X(2),1)-2*DIFF(ETA(1),X(1),1), 2*DIFF(PHI(1),U(1),1,X(1),1)+DIFF(ETA(1),X(2),1) -DIFF(ETA(1),X(1),2), DIFF(PHI(1),X(2),1)-DIFF(PHI(1),X(1),2)]]; ETA(1),ETA(2) correspond to the infinitesimals for the independent variables X(1),X(2) (i.e. space x and time t), and PHI(1) corresponds to the infinitesimal for the dependent variable U(1) (the temperature U(x,t)). To produce the determining equations in a form suitable for simplification by standard_form we start a maple session and read our maple package convprogs: > read convprogs; Then we read the determining system and give it a name: > read lode: > msys:= ": The conversion of msys to repeated index notation is accomplished by applying the function max2rep: > sys:= max2rep(msys); obtaining sys := [V(2, 3) = 0, V(2, 1) = 0, V(1, 3) = 0, V(3, 3, 3) = 0, V(2, 2) - 2 V(1, 1) = 0, 2 V(3, 3, 1) + V(1, 2) - V(1, 1, 1) = 0, V(3, 2) - V(3, 1, 1) = 0] Notice the correspondence: ETA(1)=V(1),ETA(2)=V(2),PHI(1)=V(3) and X(1)=x1, X(2)=x2, U(1)=x3. You can either save sys to a file e.g. > save sys, `heat.sys`; or read in the sfprogs package, and apply the functions in that package (see Section 4.1). Note: The message `Warning - DIFF had an argument which was neither ETA nor PHI' indicates that your system contains unspecified functions other than the infinitesimals ETA and PHI. These will require individual attention for conversion (i.e. you're probably dealing with a classification problem see the following example). 4.2.2 PRODUCING THE DETERMINING SYSTEM OF Uxx = Utt + H(U) Ut Production of the infinitesimal determining system for the class of nonlinear telegraph equations Uxx = Utt + H(U) Ut where U=U(x,t) and H is a classification function. We make up a dat file which is very similar to the heat.dat file of Section 4.2.1 (we have called this file telegraph.dat and its available for inspection in the same directory as heat.dat). When the PDE involves a classification function (e.g. H here) then that information is entered in the macsyma dat file by using the macsyma depends command (e.g. here the line depends([H], [U[1]]); was added to telegraph.dat). Following the same procedure as in Example A, we execute the commands: (c1) batch("telegraph.dat"); (c2) batch("maxconv"); The resulting file, lode which has been written to your current directory and which contains the determining system in macsyma format (after removing ` using your editor), is: listing of lode: [2, [DIFF(ETA(1),U(1),1),DIFF(ETA(2),U(1),1),DIFF(PHI(1),U(1),2), DIFF(ETA(2),X(1),1)-DIFF(ETA(1),X(2),1), DIFF(ETA(2),X(2),1)-DIFF(ETA(1),X(1),1), DIFF(PHI(1),X(2),1)*h-DIFF(PHI(1),X(2),2)+DIFF(PHI(1),X(1),2), DIFF(ETA(1),X(2),1)*h-2*DIFF(PHI(1),U(1),1,X(1),1)-DIFF(ETA(1),X(2),2) +DIFF(ETA(1),X(1),2), PHI(1)*DIFF(h,U(1),1)-DIFF(ETA(2),X(2),1)*h+2*DIFF(ETA(1),X(1),1)*h +DIFF(ETA(2),X(2),2)-DIFF(ETA(2),X(1),2) -2*DIFF(PHI(1),U(1),1,X(2),1)]]; To convert this into a format that maple and standard_form can handle, we start a maple session: > read convprogs: We read in the file that was created by our macsyma session > read lode: and give the determining system a name: > msys:= ": To use the hi-tech classification where h is regarded as a classification variable we substitute h with PHI(2): > msys:= subs(h=PHI(2), msys); msys := [2, [DIFF(ETA(1), U(1), 1), DIFF(ETA(2), U(1), 1), DIFF(PHI(1), U(1), 2), DIFF(ETA(2), X(1), 1) - DIFF(ETA(1), X(2), 1), DIFF(ETA(2), X(2), 1) - DIFF(ETA(1), X(1), 1), DIFF(PHI(1), X(2), 1) PHI(2) - DIFF(PHI(1), X(2), 2) + DIFF(PHI(1), X(1), 2), DIFF(ETA(1), X(2), 1) PHI(2) - 2 DIFF(PHI(1), U(1), 1, X(1), 1) - DIFF(ETA(1), X(2), 2) + DIFF(ETA(1), X(1), 2), PHI(1) DIFF(PHI(2), U(1), 1) - DIFF(ETA(2), X(2), 1) PHI(2) + 2 DIFF(ETA(1), X(1), 1) PHI(2) + DIFF(ETA(2), X(2), 2) - DIFF(ETA(2), X(1), 2) - 2 DIFF(PHI(1), U(1), 1, X(2), 1)] ] This ensures that when the conversion program is applied to msys: > sys:= max2rep(msys); sys := [V(1, 3) = 0, V(2, 3) = 0, V(3, 3, 3) = 0, V(2, 1) - V(1, 2) = 0, V(2, 2) - V(1, 1) = 0, V(3, 2) V(4) - V(3, 2, 2) + V(3, 1, 1) = 0, V(1, 2) V(4) - 2 V(3, 3, 1) - V(1, 2, 2) + V(1, 1, 1) = 0, V(3) V(4, 3) - V(2, 2) V(4) + 2 V(1, 1) V(4) + V(2, 2, 2) - V(2, 1, 1) - 2 V(3, 3, 2) = 0 ] h and its derivatives will be replaced by V(4) and the appropriate derivatives of V(4). At this stage you can either save sys to a file e.g. > save sys, `telegraph.sys`; or read in sfprogs and continue (see Section 4.1). 4.3 INTERFACE TO THE MAPLE PACKAGE SYMMETRY Mark Hickman [Hick93] has written a package for producing infinitesimal determining equations and solving them interactively. Using the map2rep procedure, explained in section 3.6, it is possible to convert the equations produced by this package into repeated index notation, interactively, to use the standard form programs on them. Currently, Hickman's package only is only available for release 2 of Maple 5 (for more info contact Mark Hickman at: M.Hickman@math.canterbury.ac.nz). His package exploits the postscript interface so that greek letters etc. can be displayed and calculations can be done in much the way one would do them by hand. As an example we use Hickman's package ``symmetry'' to generate the determining system for the class of nonlinear telegraph systems: Vt = Ux , Vx = C(U) Ut + B(U) with dependent variables U = U(x,t), V = V(x,t), and classification functions C(U) and B(U). First we read in Hickman's symmetry package. > read symmetry; To use Hickman's package, the dependencies of the arbitrary functions B and C on u must be entered by using his domain command > domain([u],B): > domain([u],C): Hickman's package has the procedure setup for setting up the maple environment to use his symmetry tools: > setup(); Welcome to Symmetry Please remember to end each input with a ; The independent variables are x,t; The dependent variables are u,v; The coefficients of the independent variables are > X,T; The coefficients of the dependent variables are > Y,Z; The DEs are (use the Diff operator) Diff(v,t)=Diff(u,x),Diff(v,x)=C*Diff(u,t)+B The independent variables are [x, t] The dependent variables are [u, v] The symmetry generator is / d \ / d \ / d \ / d \ P -> X |---- P| + T |---- P| + Y |---- P| + Z |---- P| \ dx / \ dt / \ du / \ dv / The differential equations are d d d / d \ ---- v = ---- u, ---- v = C |---- u| + B dt dx dx \ dt / Are these correct? y or n (do not use a ;) y Note that the jet form of the equations are given by the variable eqn > eqn; [vt = ux, vx = C ut + B] To find the determining system, we use the procedure geneqn: > eq:=op(geneqn(eqn)); / d \ / d \ / d \ eq := - |---- Y| + |---- X| ux + |---- T| ut \ dx / \ dx / \ dx / // d \ / d \ / d \ \ - ux ||---- Y| - |---- X| ux - |---- T| ut| \\ du / \ du / \ du / / // d \ / d \ / d \ \ / d \ - (C ut + B) ||---- Y| - |---- X| ux - |---- T| ut| + |---- Z| \\ dv / \ dv / \ dv / / \ dt / / d \ / d \ - |---- X| (C ut + B) - |---- T| ux \ dt / \ dt / // d \ / d \ / d \ \ + ut ||---- Z| - |---- X| (C ut + B) - |---- T| ux| \\ du / \ du / \ du / / // d \ / d \ / d \ \ + ux ||---- Z| - |---- X| (C ut + B) - |---- T| ux|, \\ dv / \ dv / \ dv / / / / d \ / d \\ / d \ / d \ / d \ Y |- |---- C| ut - |---- B|| + |---- Z| - |---- X| (C ut + B) - |---- T| ux \ \ du / \ du // \ dx / \ dx / \ dx / // d \ / d \ / d \ \ + ux ||---- Z| - |---- X| (C ut + B) - |---- T| ux| \\ du / \ du / \ du / / // d \ / d \ / d \ \ / d \ + (C ut + B) ||---- Z| - |---- X| (C ut + B) - |---- T| ux| - (|---- Y| \\ dv / \ dv / \ dv / / \ dt / / d \ / d \ // d \ / d \ / d \ \ - |---- X| ux - |---- T| ut + ut ||---- Y| - |---- X| ux - |---- T| ut| \ dt / \ dt / \\ du / \ du / \ du / / // d \ / d \ / d \ \ + ux ||---- Y| - |---- X| ux - |---- T| ut|) C \\ dv / \ dv / \ dv / / In Hickman's package, the names for the infinitesimals are indexed quantities in the maple environment. Such a quantity (M say) can be converted into Maple's function notation by making the substitution M=op(0,M)(op(M)). Thus to carry out this conversion on the system above we make the substitutions: > neweq:=subs(X=op(0,X)(op(X)),T=op(0,T)(op(T)),Y=op(0,Y)(op(Y)),Z=op(0,Z)(op( Z)),B=op(0,B)(op(B)),C=op(0,C)(op(C)),[eq]); neweq := [ / d \ / d \ / d \ - |---- Y(x, t, u, v)| + |---- X(x, t, u, v)| ux + |---- T(x, t, u, v)| ut - \ dx / \ dx / \ dx / // d \ / d \ / d \ \ ux ||---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut| \\ du / \ du / \ du / / - (C(u) ut + B(u)) // d \ / d \ / d \ \ ||---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut| \\ dv / \ dv / \ dv / / / d \ / d \ + |---- Z(x, t, u, v)| - |---- X(x, t, u, v)| (C(u) ut + B(u)) \ dt / \ dt / / d \ / d \ - |---- T(x, t, u, v)| ux + ut (|---- Z(x, t, u, v)| \ dt / \ du / / d \ / d \ - |---- X(x, t, u, v)| (C(u) ut + B(u)) - |---- T(x, t, u, v)| ux) + ux ( \ du / \ du / / d \ / d \ |---- Z(x, t, u, v)| - |---- X(x, t, u, v)| (C(u) ut + B(u)) \ dv / \ dv / / d \ - |---- T(x, t, u, v)| ux), \ dv / / / d \ / d \\ / d \ Y(x, t, u, v) |- |---- C(u)| ut - |---- B(u)|| + |---- Z(x, t, u, v)| \ \ du / \ du // \ dx / / d \ / d \ - |---- X(x, t, u, v)| (C(u) ut + B(u)) - |---- T(x, t, u, v)| ux + ux ( \ dx / \ dx / / d \ / d \ |---- Z(x, t, u, v)| - |---- X(x, t, u, v)| (C(u) ut + B(u)) \ du / \ du / / d \ / d \ - |---- T(x, t, u, v)| ux) + (C(u) ut + B(u)) (|---- Z(x, t, u, v)| \ du / \ dv / / d \ / d \ - |---- X(x, t, u, v)| (C(u) ut + B(u)) - |---- T(x, t, u, v)| ux) - ( \ dv / \ dv / / d \ / d \ / d \ |---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut + \ dt / \ dt / \ dt / // d \ / d \ / d \ \ ut ||---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut| \\ du / \ du / \ du / / + // d \ / d \ / d \ \ ux ||---- Y(x, t, u, v)| - |---- X(x, t, u, v)| ux - |---- T(x, t, u, v)| ut| \\ dv / \ dv / \ dv / / ) C(u) ] The infinitesimals are now expressed in maple function notation. The above command could have also been carried out using > neweq:= subs(map(proc(z) z=op(0,z)(op(z)) end, [X,T,Y,Z,B,C]), [eq]); neweq should be saved to a file now, and a new Maple session started since Hickman's array enviroment otherwise causes problems for later calculations. Once in the new session, we read in the conversion programs: >read convprogs; to apply the map2rep procedure (see 3.6) it is necessary to have B and C as functions of all the variables before using the procedure, so > sys1:= subs(C(u)=C(x,t,u,v), B(u)=B(x,t,u,v), neweq): > sys1:=map2rep(sys1,[X,T,Y,Z,B,C],[x,t,u,v,ux,ut]); sys1:= [- V(3, 1) + V(1, 1) x5 + V(2, 1) x6 - x5 (V(3, 3) - V(1, 3) x5 - V(2, 3) x6) - %1 (V(3, 4) - V(1, 4) x5 - V(2, 4) x6) + V(4, 2) - V(1, 2) %1 - V(2, 2) x5 + x6 (V(4, 3) - V(1, 3) %1 - V(2, 3) x5) + x5 (V(4, 4) - V(1, 4) %1 - V(2, 4) x5), V(3) (- V(6, 3) x6 - V(5, 3)) + V(4, 1) - V(1, 1) %1 - V(2, 1) x5 + x5 (V(4, 3) - V(1, 3) %1 - V(2, 3) x5) + %1 (V(4, 4) - V(1, 4) %1 - V(2, 4) x5) - (V(3, 2) - V(1, 2) x5 - V(2, 2) x6 + x6 (V(3, 3) - V(1, 3) x5 - V(2, 3) x6) + x5 (V(3, 4) - V(1, 4) x5 - V(2, 4) x6)) V(6) ] %1 := V(6) x6 + V(5) To apply the standard form package these need to be equations: > sys1:= map(proc(z) z=0 end, sys1): Finally, to obtain the determining equations, we polynomially decompose the above system on x5 and x6 (i.e, ux and ut) using the sfprogs procedure polyd (see section 7.2): > read sfprogs; > polyd(sys1,[x5,x6]); [- V(3, 1) - V(5) V(3, 4) - V(1, 2) V(5) + V(4, 2) = 0, V(1, 1) - V(3, 3) + V(4, 4) - V(2, 2) = 0, - V(1, 2) V(6) + V(2, 1) + V(5) V(2, 4) - V(6) V(3, 4) + V(4, 3) - V(1, 3) V(5) = 0, V(1, 3) - V(2, 4) = 0, V(6) V(2, 4) - V(1, 3) V(6) = 0, - V(3) V(5, 3) - V(3, 2) V(6) + V(4, 1) - V(1, 1) V(5) + V(5) (V(4, 4) - V(1, 4) V(5)) = 0, - V(2, 1) - (- V(1, 2) + V(3, 4)) V(6) - V(1, 3) V(5) - V(5) V(2, 4) + V(4, 3) = 0, - V(3) V(6, 3) - (- V(2, 2) + V(3, 3)) V(6) - V(1, 1) V(6) - V(5) V(1, 4) V(6) + V(6) (V(4, 4) - V(1, 4) V(5)) = 0, 2 - V(2, 3) + V(1, 4) V(6) = 0, V(2, 3) V(6) - V(6) V(1, 4) = 0, - V(6) V(2, 4) - (- V(1, 3) - V(2, 4)) V(6) - V(1, 3) V(6) = 0] These then are the determining equations for the symmetry generators of the nonlinear telegraph system. They are, in fact, quite difficult to reduce to standard form, and we go through this computation in section 9.2. We also mention that the 2 original determining equations (before polynomial decompostion) could have been entered directly into the standard_form provided that the expressions: [V(1),1,2,3], [V(2),1,2,3], [V(3),1,2,3], [V(4),1,2,3], expressing the dependency info was appended to the above equations (i.e. that each V(i) is independent of x4 and x5) are added to the system list. In addition the dependency declarations [V(5),3], [V(6),3] and are inserted in the class_eqns list since this is a classification problem (see section 9.2). Finally the Hickman package can also be used in batch mode: Batch template for the Hickman Program Symmetry: read symmetry; # read Hickman's package symmetry domain([u],B); # dependency of classification variable domain([u],C); # dependency of classification variable varlist := [x,t]; # independent variables INDEPCOEFF:= [X,T]; # their corresponding infinitesimals dependvar := [u,v]; # dependent variables DEPCOEFF := [Y,Z]; # their corresponding infinitesimals eqn:= [Diff(v,t)=Diff(u,x),Diff(v,x)=C*Diff(u,t)+B]; # differential equations (solved-form) vars:=op(varlist),op(dependvar); domain([vars],op(INDEPCOEFF),op(DEPCOEFF)); INDEPCOEFF:=[seq(op(i,INDEPCOEFF)[vars], i=1..nops(INDEPCOEFF))]; DEPCOEFF:=[seq(op(i,DEPCOEFF)[vars], i=1..nops(DEPCOEFF))]; eqn:=jetform(eqn); eqn:= geneqn(eqn); # unpolynomially decomposed determining system end of template 4.4 INTERFACE TO THE MAPLE PACKAGE LIESYMM The map2rep procedure explained in section 3.6 allows for an easy interface between standard form and the maple symmetry package liesymm. It can be used to convert the determining equations produced by the liesymm package into repeated index notation which can then be analyzed by the standard form programs. Documentation for the liesymm package can be found in [HaE71, CamDevFee92], as well as online documentation. We illustrate the procedure with the example below. Note that this is in release 2 of Maple V. As an example, we find the determining equations for the linear heat equation, using the liesymm package, and then use the map2rep procedure to put them in repeated index notation. To begin with, we start a maple session, and read in the liesymm package using with: > with(liesymm); Our equation is the linear heat equation: > sys:=[Diff(u(x,t),t)-Diff(u(x,t),x,x)=0]; / d \ sys := [|---- u(x, t)| - Diff(u(x, t), x, x) = 0] \ dt / Note that when using liesymm, it is necessary to write the DEs using maple's inert differential operator Diff. To find the determining system we use the liesymm procedure determine. > dets:=determine(sys,V,[u(x,t)],k); 2 2 d d dets := {----- V2(x, t, u) = 0, ----- V1(x, t, u) = 0, 2 2 du du / 2 \ d | d | ---- V1(x, t, u) = - |------- V2(x, t, u)|, du \ du dx / 2 d d ----- V3(x, t, u) = ---- V3(x, t, u), 2 dt dx 2 / 2 \ d | d | / d \ ----- V1(x, t, u) = 2 |------- V3(x, t, u)| + |---- V1(x, t, u)|, 2 \ du dx / \ dt / dx 2 d / d \ / d \ ----- V2(x, t, u) = |---- V2(x, t, u)| - 2 |---- V1(x, t, u)|, 2 \ dt / \ dx / dx d d ---- V2(x, t, u) = 0, ---- V2(x, t, u) = 0, du dx 2 / 2 \ d | d | ----- V3(x, t, u) = 2 |------- V1(x, t, u)|} 2 \ du dx / du The determine procedure takes 4 arguments. The first is the set or list of differential equations. The second is a name to give to the infinitesimal generators (i.e. the dependent variables in the determining equations). We have used V. The third is a list of the dependent and independent variables entered in the form of dependent variables as functions of the independent variables (as above). The fourth is a name for the extended variable names. The infinitesimals are numbered such that, if there are n independents, and m dependents, then the first n infinitesimals are for the independents, numbered in the order they appear in the argument lists of the dependents, and the remaining m are for the dependents in the order they were enetered in the list. Now that we have the determining equations, we can convert them to repeated index notation using the map2rep procedure. However, we cannot enter dets directly, because the liesymm procedure uses an inert differentiation operator, that appears in the determining equations dets. To force these operators to become active, we must use the liesymm procedure value(). So, we enter > newdets:=map2rep(value(dets),[V1,V2,V3],[x,t,u]); newdets := {V(1, 3, 3) = 0, V(2, 3, 3) = 0, V(2, 1) = 0, V(2, 3) = 0, V(2, 1, 1) = - 2 V(1, 1) + V(2, 2), V(1, 3) = - V(2, 1, 3), V(3, 2) = V(3, 1, 1), V(3, 3, 3) = 2 V(1, 1, 3), V(1, 1, 1) = V(1, 2) + 2 V(3, 1, 3)} and these are the determining equations in repeated index notation. 5. INTEGRATION OF THE EQUATIONS The focus of this package has been on canonically simplifying systems of PDEs. Other packages have gone a long way in implementing methods for integrating PDEs. It is hoped to develop interfaces to such packages in the future to take advantage of their integrators (e.g. the powerful package of Wolf [Wolf89, Wolf91, Wolf92]). In [Reid91c] we present an algorithm for integrating systems on the basis of directed graphs associated with their standard forms, and in [Reid92b] we give an algorithm for solving systems of PDEs in their coefficient fields by using standard_form to recursively decouple and solve ODEs. 5.1 USE OF taylor_exp Occasionally systems of PDEs possess polynomial solutions. In particular often the determining systems for infinitesimals possess solutions which are second order (or sometimes third order) polynomial functions. By using taylor_exp to third order (or sometimes 4-th order) one can identify whether the series has terminated at order 2 (or order 3) in view of the observation that: If all the order K terms vanish in the taylor series representation to order K of a function f(x) about an arbitrary point x=x0, (i.e. they vanish for all values of x0) then the series terminates at order K-1. e.g. in the one variable case: no K-th order terms for any x0 in (K) f(x) = f(x0) + f'(x0) (x-x0) + ... + f (x0) (x-x0)^K / K! + O(K+1) (K) (J) implies f (x0) = 0 for all x0. So f (x0) = 0 for J > K and f(x) is a polynomial of degree at most K. e.g. y''(x) = 0 has initial data y(x0) = a1, y'(x0) = a2 and taylor series to order 2: y(x) = a1 + a2 (x-x0) + O(3) This obviously isn't true if the expansion is about particular point e.g. the series for sin(x) = x + x^3/6 + O(5) doesn't terminate even though it has no 4-th order terms. A completely algorithmic method for finding all polynomial/rational solutions of linear PDE systems with polynomial/rational coefficients or more generally all solutions belonging to some extension of the rational function field is given in [Reid92b]. 5.2 INTERACTIVE INTEGRATION At any stage you can integrate simple PDEs by hand or by a package (e.g. maple's dsolve). e.g. V(3,1,1) = 0 integrates to V(3) = V(5) + V(6)*x1, V(5,1) = 0, V(6,1) = 0 Appending the above 3 equations to the system you wish to standard_form or solve may help in the process. 6. CHANGING VARIABLES IN SYSTEMS OF PDEs The procedure transform can be used to perform a general (invertible) transformation of dependent and independent coordinates on a system of PDEs. The basic formula for calculating the prolonged transformations on partial derivatives of the dependent coordinates with respect to the independent coordinates is given in [BlK89]. This highly recursive formula is the core of the transform procedure, which can express a given PDE in any coordinates given an invertible transformation from the original coordinates. It is also possible to display the little and big jacobians of the transformation. The little jacobian is the total differential (i.e. regarding the dependent variables as functions of the independents) of the transformation for the independent variables with respect to the new independent variables. The big jacobian is the jacobian of the entire transformation with respect to the independent and dependent coordinates. The transform procedure works in regular maple notation, not repeated index notation. Hence coordinate transformations can be done in a more natural way. The transform procedure can be used together with the conversion procedures described in section 5.6 to use the standard form programs on the new equations, or to change variables on the standard form of a system. The basic syntax of the transform procedure is > transform(eqns,newindeps,indeptrans,deptrans); where: eqns is an equation or system of equations to be transformed; newindeps is a list of the new independent variables; indeptrans is a list of the transformations on the independent variables; deptrans is a list of the transformations for the dependent variables The transformations in indeptrans and deptrans must be expressed as equalities, with the transformation, in terms of the new variables, on the right hand side, and the old variable to be replaced, on the left hand side. For instance, if the existing independent variable and dependent variable was x and u(x) respectively, and we wished to make a general linear change to new independent and dependent variables y and v(y), then indeptrans:=[x=a11*y+a12*v(y)], and deptrans:=[u(x)=a21*y+a22*v(y)]. The dependencies must be included for both the new and old dependent variables. Procedure transform will return the system of equations eqns, expressed in the new coordinates y, v(y). To have transform return either the little or big jacobian, a fifth argument is included whose value is either little, big, or biglittle to have either the little jacobian, big jacobian, or both, respectively, returned. If either big or biglittle is the fifth argument, then a sixth argument is included, which must be a list of the new dependent variables, written without dependencies in this case. If any of these options are invoked, transform will return a list of 2 or 3 elements whose first element is the transformed equations and whose remaining elements are the requested jacobians. These options and the general use of transform are demonstrated in the examples below. Example A: As a first example, we wish to transform the equation ut - 2 t uxx=0, under the coordinate transformation x=ay, t=K(s), and u(x,t)=v(y,s), where K(s) is an arbitrary function and a is a parameter. > read transform; # read in the transform package > sys:=[diff(u(x,t),t)-2*t*diff(u(x,t),x,x)=0]; # define the system sys := [D[2](u)(x, t) - 2 t D[1, 1](u)(x, t) = 0]. Then, we invoke transform, to create the transformed system in the new coordinates > tsys:=transform(sys,[y,s],[x=a*y,t=K(s)],[u(x,t)=w(y,s)]); D[2](w)(y, s) K(s) D[1, 1](w)(y, s) tsys := [------------- - 2 --------------------- = 0] d 2 ---- K(s) a ds Example B: Here, we will transform the equation H(Yx) Yxx = Yt where H is an arbitrary function, by the hodograph x=Y', t=t', Y=x'. Also, we will display the little and big jacobians. We enter the system > sys:=[H(diff(Y(x,t),x))*diff(Y(x,t),x,x)=diff(Y(x,t),t)]; sys := [H(D[1](Y)(x, t)) D[1, 1](Y)(x, t) = D[2](Y)(x, t)] and then invoke transform > transform(sys,[y,s],[x=Z(y,s),t=s],[Y(x,t)=y],biglittle,[Z]); table([ 1 H(-------------) D[1, 1](Z)(y, s) D[1](Z)(y, s) D[2](Z)(y, s) 1 = [- --------------------------------- = - -------------] 3 D[1](Z)(y, s) D[1](Z)(y, s) [ 1 ] [ ------------- 0 ] [ D[1](Z)(y, s) ] 2 = [ ] [ D[2](Z)(y, s) ] [ - ------------- 1 ] [ D[1](Z)(y, s) ] [ 0 0 1 ] [ ] 3 = [ 0 1 0 ] [ ] [ 1 0 0 ] ]) The first element is the transformed equation, and the second and third elements are the little and big jacobians respectively. Note that the system is linearized when H(Zy)=1/(Zy)^2. 7. TOOLBOX OF INTERACTIVE FUNCTIONS In future we'll be making some of the internal routines of the sfprogs package available for interactive work. Clearderivdep, below, is a useful example of such a routine. 7.1 Clearderivdep Given a system in standard form sf, then the command > read sfprogs; > Clearderivdep(sf, sys); clears the list of equations sys of all dependencies on sf. (sf and sys must be entered as lists of equations in repeated index notation). Example: > read sfprogs; > Clearderivdep([V(3)=V(5)+V(6)*x1,V(5,1)=0,V(6,1)=0], [V(3,1,1)=0, V(3,1) = V(3,1)^2 + V(4)]); yields: 2 [0 = V(6) + V(4) - V(6)] (note that the first equation simplified to 0=0 and was omitted from the output). Also see Section 5.2. 7.2 polyd Procedure polyd Input: list of equations which is polynomial in pvars, and the list of variables pvars. Output: list of polynomially decomposed eqns with respect to pvars. Example: > read sfprogs; > polyd([t1^2*y+t1*t2*k=h ,(t1+t2)*u=t1*y],[t1,t2]); yields [- h = 0, y = 0, k = 0, u = 0, u - y = 0] 7.3 sfdiff sfdiff(exp, string) where exp is an expression or equation in repeated index notation and string is a string of independent variables from the set of independent variables {x1,x2, ... ,xm}. Example: > read sfprogs; > sfdiff(V(1,1,1)^2, x1); 2 V(1, 1, 1) V(1, 1, 1, 1) > sfdiff(V(2)^2 + V(1,3)=0, x1,x3); 2 V(2, 3) V(2, 1) + 2 V(2) V(2, 1, 3) + V(1, 1, 3, 3) = 0 8. STRATEGIES TO DEAL WITH MEMORY EXPLOSION In general, memory explosion seems to occur when the system to which standard_form is applied is very complex (algebraically), or has a large number of equations. We discuss strategies for dealing with these problems (also see the the interactive programs of Kersten [Ker87] and the automatic programs of Wolf [Wolf89, Wolf91, Wolf92] for interesting approaches to the problem of controlling memory swell). 8.1 ADJUSTING THE EQUATION-FLOW SETTINGS Changing the equation-flow settings to their ideal can be a system dependent process. Often adjusting the knobs to their Ultra-Conservative setting (i.e. one equation at a time) will lead to termination using a reasonable quantity of memory, However, some which cases blow up in the conservative setting terminate in a reasonable time using little memory in the default setting. This occurs in systems having one or more large equations which dramatically reduce the complexity of the system once solved and substituted for. These situations are not easy to detect, and currently the best choice is an experimental approach. Another (desparate strategy!) would be to change that part of the code which controls the way in which easy equations are chosen from unclassified equations (see 3.3.2). 8.2 VARYING THE ORDERING Changing the ordering using the janet_order option can lead to drastic reductions in memory and run time. We are not aware of any general strategy for choosing the ordering to obtain such reductions. The main reason for this is the difficulty in predicting the path on which the solution will be arrived at. Again the best method currently seems to be experimental. Sometimes ordering with respect to the dependent variables last, and differentiation variable first, minimizes memory expansion and produces the best run times. Sometimes total ordering after the differentiation variable ordering can be effective. e.g. 3 independent variables , 3 dependent variables set > ordertab1:= [[1,0,0, 0,0,0], [0,1,0, 0,0,0], [0,0,1, 0,0,0]]: > standard_form(, janet_order = ordertab1 ) and if this does not work, one could permute order of the three lists inside ordertab1. 8.3 CLASSIFICATION VARIABLES AND DIFFERENTIAL EXTENSIONS (or polynomials are nicer) Maple and other computer algebra systems can always reduce polynomial or rational functions to normal form. However, when a system possesses a number of non-polynomial functions such as sin(x1), exp(x3),and x4^(1/2), there is a great deal of memory consumption involved in simplification and manipulation of such expressions (in general Maple has no canonical way of simplifying expressions involving such functions). One way to get eliminate a nonpolynomial function is to introduce a new classification variable which represents the function, and specify its behavior with a differential constraint ( i.e. define the function by a differential extension). For example consider a system of PDEs with 5 dependent variables, 5 independent variables, and a painful recurrence of the factor x4^(1/2). Substituting x4=V(6)^2 into the system and introducing the class equation class_eqns=[[V(6),4],V(6,4)=1/(2*V(6))] which defines the derivative properties of V(6) eliminates the square root and replaces it with pleasant rational occurences of V(6). Alternatively, a non-standard method (courtesy of ADW), giving a linear form for the differential constraint could be used. For example: First substitute x4^(1/2)=V(6) then set: class_eqns=[[V(6),4],V(6,4)=V(6)/(2*x4)] (i.e. V(6)=c2*x^(1/2) where c2 is an arbitrary constant). However, one must ensure that none of the pivots vanish. For example, if one of the pivots was: V(6)^2-x4=0 then an invalid operation has occurred and the calculation must be done again with a suitably different class equation. 8.4 THE BIG GUNS - STORAGE AND USEREASY Okay... So we're desperate now. Don't sweat it yet, there are still a couple of ways to obtain the answer. However, do not read further if you don't expect to have to work for it. Still Here? Okay, here it is you glutton for punishment: One can use the storage options (see 3.4) to store the system at the most recent iteration (store) or if you don't mind clogging up your filespace at every iteration (storeall). You can also set the equation flow parameters at ultra-consevative. And hack at the system a piece at a time. In order to do this, one requires a working knowledge of Maple's list handling features such as op(), susbop(), map() to manipulate the lists (also see section 7). Oh, and lest we forget that many of the options mentioned previously in this manual will come to bear on this strategy, so go read the first 7 chapters again. One strategy which GR has used with some success is to use the store option to recover the system in at some chosen stage or at the last iteration before the WOD (our favorite phrase for memory explosion). Then choose the simplest equations from the system, throw the ugly ones out and start the original system together with these simplest equations and perhaps use the usereasy option (see 3.4.2). You're Back Already?!? OK I trust you. The storeall option (3.4) allows one to retrieve the stored system before the WOD was reached. So you've now run the code in ultra-conservative setting (see (3.3.2)) with the infolevel at its maximum (i.e. strategy = [1.75, 0.9, 1.0]). You should have as many StandardStorage_i.m files in your directory as iterations the program completed (Usually 50~80 of them for large systems) DON'T PANIC! This is normal. So now you begin another Maple session (Note: it is important to start with a fresh slate if the program went down in flames) load the code, and use readstored to load the most recent version (highest numbered store-file) with: > L:=readstored(`StandardStorage_'i'`): where 'i' is the highest number available. So now you have the system, i.e. OT,Ez,UnC,_pivs and possibly NL all stored in a list with name L (you can refer to 4.3 to see which is which). A beneficial move at this point would be to issue the command: > assign(L); which will now assign the list of one term equations to the variable OT, the list of Easy equations to the variable Ez, etc. Now comes the fun part! Try: > map(length,Ez); > map(length,OT); You should now have a couple of lists displayed, and they possess the length (a crude measure of complexity) of the corresponding elements in Ez and OT respectively. You can look at each individual element (say the i'th element of Ez) with: > Ez[i]; or > op(i,Ez); Good now we're ready for the real work. Since there was an explosion, there is likely a number of elements of Ez with relatively large length. Start a new variable: > tough:=NULL; and pull undesirable or long elements out of Ez from the last to first using: > tough:=tough,Ez[i]: > Ez:=subsop(i=NULL,Ez): In this way you are moving long or complicated equations out of Ez and appending them to tough ( it is imperative not to lose equations on the way ) Okay, all the real nasties out? You can check this with another map command: > map(length,Ez) Good! Now turn Ez into an appendable variable by typing: > Ez:=op(Ez): now > map(length,UnC); Look for the simple equations in the UnClassified (it's good to manually check these also; sometimes even long equations can be simple) and a good rule is the fewer the classification variables, and the smaller the coefficients (as far as the xi polynomials go) the better. To transfer them to the Ez list do as follows: > Ez:=Ez,UnC[i]: > UnC:=subsop(i=NULL,UnC): Put it all together and give it a test drive! Okay, first step : We know that we like the Ez equations and the OT equations so let's put them together: > Ez:=op(OT),Ez: Now we want to know how many good equations there are, we use: > nops([Ez]); Remember the number you now see, we call it numeasy for reference. Now we put everything together into a variable called sys: > sys:=[Ez,op(UnC),tough,op(NL)]: > firstpivs:=_pivs: and finally: > save sys,firstpivs,temporaryfile; We got it right where we want it now! Remember that number? numeasy? You need it in the following statement: > sf:=standard_form(sys,usereasy=numeasy,storeall=`Attempt2`): This will now fire off the system and store the system at the i'th iteration in Attempt2_i.m, assuming that the first numeasy number of equations are fairly easy, and will solve those first. Hopefully the nasty equations which caused the blowup will now simplify to a reasonable, workable level. SPECIAL NOTES: If it is important to know the pivots of the system, some of those are now stored in temporaryfile in the variable firstpivs, and the others since the re-run are being built up in _pivs. Both lists are important and result from the current problem being solved. In the new code we have not yet found a linear problem of practical origin for which we've had to resort to this drastic measure, so be sure to exhaust the other possibilities first. Also if you find a problem of this magnitude, we would be interested in knowing of it. 9. NONTRIVIAL EXAMPLES 9.1 A BIG EXAMPLE: THE MHD EQUATIONS A description on how to do this type of problem was already introduced in a simple example in sections 4.2.1 and 4.1.1 of the manual. As an example of our programs working on a nontrivial (large) example we generate the determining system for symmetries of the the magneto-hydro-dynamic (mhd) equations (below). The computation of the standard form, although a large one, terminates in about 50 minutes, on a Silicon Graphics machine, and remarkably never uses more than about 2.2 Megabytes of RAM (Sparc Stations give comparable performance). MHD EQUATIONS: f + (v.Grad) f + f Div v = 0 t f ( v + (v.Grad) v ) + Grad ( p + (H.H)/2 ) - (H.Grad ) H = 0 t H + (v.Grad) H + H (Div v) - (H.Grad) v = 0 t Div H = 0 K K ( p / f ) + (v.Grad) (p / f ) = 0 t Above pressure=p=p(x,y,z,t), mass density = f = f(x,y,z,t), coefficient of viscosity = K, vector fluid velocity v=v(x,y,z,t)=(Vx,Vy,Vz) and magnetic field H=H(x,y,z,t)=(Hx,Hy,Hz). Also Grad = (d/dx,d/dy,d/dz), and . denotes dot product. For the macsyma program we use the correspondence: x=x[1], y=x[2], z=x[3], t=x[4], f = u[1], p = u[2], Vx=u[3], Vy=u[4], Vz=u[5], Hx=u[6], Hy=u[7], Hz=u[8]. with corresponding symmetry generator: L = ETA[1] d/dx + ETA[2] d/dy + ETA[3] d/dz + ETA[4] d/dt + PHI[1] d/df + PHI[2] d/dp + PHI[3] d/dVx + PHI[4] d/dVy + PHI[5] d/dVz + PHI[6] d/dHx + PHI[7] d/dHy + PHI[8] d/dHz The corresponding data file for the macsyma program is stored under mhd.dat in the same directory as sfprogs. Batching a file containing the (macsyma) commands: (c1) batch("mhd.dat"); (c2) batch("maxconv"); (c3) quit(); produces a file lode containing the mhd determining equations in macsyma syntax. Following the same path as in section 4.2.1 and 4.1.1 of the manual: Once the mhd determining equations have been generated the file lode is stripped of ' using an editor. We start a maple session by loading the conversion programs and the stripped lode: > read convprogs; > read lode; Then give the determining system a name: > msys:= ": The conversion of msys to repeated index notation is accomplished by applying the function max2rep: > sys:= max2rep(msys); Notice the correspondence: ETA(i)=V(i), i=1,2,3,4 PHI(1)=V(5),PHI(2)=V(6),...,PHI(8)=V(12), and x1=x[1], x2=x[2], x3=x[3], x4=x[4], x5=u[1], x6=u[2], x7=u[3], ... , x12=u[8] and save the converted system to a file (called mhd.sys here), which you will find in the same directory as sfprogs. > save sys, `mhd.sys`; We start a new maple session, read the standard form package > read sfprogs; and the converted system > read `mhd.sys`; Next the determining system is reduced to standard form: > sf:= standard_form(sys); The resulting standard form obtained in about 1 hour of cpu on a SGI machine, which we have not displayed has 133 equations - 102 of which are one-term equations The original system had 42 one-term equations. The shortened pivots are > spivs:= shortenpiv(_pivs); 2 2 2 [K, x11, x10, x12, x5, x8, x7, x6, x10 + x12 + x11 , - x8 x12 + x9 x11,- 2 + K, - x7 x12 + x9 x10, 4 2 2 2 2 2 2 4 x11 + 4 x11 x12 - x10 x11 - 2 x12 x10 - 2 x10 , - 1 + K, - 3/2 + K, 2 2 2 2 2 x11 + x12 + 1/4 x10 , x11 + 1/2 x12 , 4 2 2 2 2 4 2 2 x11 - 1/2 x11 x12 - x10 x11 + 5/2 x10 + x12 x10 , 2 2 x11 - 2 x10 ] So this standard form is valid if these shortened pivots do not vanish identically that is the standard form is valid for K # 0,1,3/2,2. We calculate its initial data: > id:= initial_data(sf); yielding [[V(1) = a1, V(2) = a2, V(3) = a3, V(4) = a4, V(4, 4) = a5, V(5) = a6, V(7) = a7, V(8) = a8, V(9) = a9, V(10) = a10, V(11) = a11, V(11, 12) = a12, V(12) = a13], []] So if K # 0,1,3/2,2 the Lie group is 13-dimensional (i.e. the Lie symmetry group is finite-dimensional). The commutator table for the symmetries L(1),L(2),...,L(13) corresponding to this initial data a1,a2,...,a13 respectively is given by: > ct:= fin_com_table(sf,id); The table is a 13 X 13 array, which we have not displayed because of its complexity. The quantities X1,X2,X3,...,X11,X12 appearing in the table are the coordinates of the arbitrary point through at which the initial data is posed. This table was calculated before the exact form of the generators was known and we will give it in simplified form at the end of the session. To calculate the generators we use the command finite_generators, which uses the taylor_exp routine (see 4.1.1). In particular > fg1:= finite_generators(sf,id,1); gives the taylor expansion of the generators to order 1: fg1 := [L(1) = Dx(1), L(2) = Dx(2), L(3) = Dx(3), L(4) = Dx(4), L(5) = h1 Dx(1) + h2 Dx(2) + h3 Dx(3) + h4 Dx(4), / h1 X7 h4\ / X8 h4 h2 \ L(6) = |- 1/2 ---- + 1/2 -----| Dx(1) + |1/2 ----- - 1/2 ----| Dx(2) \ X5 X5 / \ X5 X5 / / h3 X9 h4\ / h5 \ h7 Dx(7) + |- 1/2 ---- + 1/2 -----| Dx(3) + |1 + ----| Dx(5) - 1/2 -------- \ X5 X5 / \ X5 / X5 h8 Dx(8) h9 Dx(9) - 1/2 -------- - 1/2 --------, X5 X5 L(7) = h4 Dx(1) + Dx(7), L(8) = h4 Dx(2) + Dx(8), L(9) = h4 Dx(3) + Dx(9), L(10)= ..., L(11)=..., L(12)=..., L(13)=...] where h1=(x1-X1), h2=(x2-X2), etc and L(10), L(11), ... , L(13) are long expressions. These generators contain terms of order 1 in the hi's so by section 5.1 of the manual we can't be sure that they are exact (i.e. that there are no higher order terms in the expansion). Expanding the generators to order 2 > fg2:= finite_generators(sf,id,2); we find that there are no order 2 terms in the generators (by inspection or by checking that fg1=fg2 using the Maple boolean command evalb(fg1=fg2)). So the generators fg1 are exact (see section 5.1 of the manual). In terms of the independent variables x1,...,x12 and the point (X1,X2,...,X12) at which the initial data is posed the generators are given by > subslist1:= [h1=x1-X1,h2=x2-X2,h3=x3-X3,h4=x4-X4, h5=x5-X5,h6=x6-X6,h7=x7-X7,h8=x8-X8, h9=x9-X9,h10=x10-X10,h11=x11-X11,h12=x12-X12]: > fg1:= subs(subslist1, fg1): and look rather horrible. To simplify their appearance we use our freedom to choose the point (X1,X2,...,X12) at which the initial data is posed. Usually X1=0,X2=0,...,X12=0 is the best choice but the generators and standard form are singular when X5=0 and X10=0, so we put X5=1, X10=1 and the remaining Xi's to zero: > subslist2:= [X1=0,X2=0,X3=0,X4=0, X5=1,X6=0,X7=0,X8=0, X9=0,X10=1,X11=0,X12=0]: > fg1:= subs(subslist2, fg1); obtaining the generators [L(1) = Dx(1), L(2) = Dx(2), L(3) = Dx(3), L(4) = Dx(4), L(5) = x1 Dx(1) + x2 Dx(2) + x3 Dx(3) + x4 Dx(4), L(6) = - 1/2 x1 Dx(1) - 1/2 x2 Dx(2) - 1/2 x3 Dx(3) + x5 Dx(5) - 1/2 x7 Dx(7) - 1/2 x8 Dx(8) - 1/2 x9 Dx(9), L(7) = x4 Dx(1) + Dx(7), L(8) = x4 Dx(2) + Dx(8), L(9) = x4 Dx(3) + Dx(9), L(10) = x1 Dx(1) + x2 Dx(2) + x3 Dx(3) + 2 x6 Dx(6) + x7 Dx(7) + x8 Dx(8) + x9 Dx(9) + x10 Dx(10) + x11 Dx(11) + x12 Dx(12), L(11) = - x2 Dx(1) + x1 Dx(2) - x8 Dx(7) + x7 Dx(8) - x11 Dx(10) + x10 Dx(11), L(12) = x3 Dx(2) - x2 Dx(3) + x9 Dx(8) - x8 Dx(9) + x12 Dx(11) - x11 Dx(12), L(13) = - x3 Dx(1) + x1 Dx(3) - x9 Dx(7) + x7 Dx(9) - x12 Dx(10) + x10 Dx(12)] These results agree with those given by other investigators (hereman [Here93]). [in more physical notation x1=x,x2=y,x3=z,x4=t,x5=f,x6=p,x7=Vx,x8=Vy,x9=Vz, x10=Hx,x11=Hy,x12=Hz so that L(5)=x d/dx + y d/dy + z d/dz + t d/dt etc] Recall that the commutator table ct was calculated already. We simplify it by wisely choosing our initial point using the same substitution as above: > L:= proc(k) cat(L,k) end: # this procedure shortens L(1) to L1, etc. > ct:= map(simplify, subs(subslist2, op(ct))); L1 L2 L3 L4 L5 L6 L7 L8 L9 L10 L11 L12 L13 L1 [ 0 0 0 0 L1 - 1/2 L1 0 0 0 L1 L2 0 L3 ] [ ] L2 [ * 0 0 0 L2 - 1/2 L2 0 0 0 L2 - L1 - L3 0 ] [ ] L3 [ * * 0 0 L3 - 1/2 L3 0 0 0 L3 0 L2 - L1 ] [ ] L4 [ * * * 0 L4 0 L1 L2 L3 0 0 0 0 ] [ ] L5 [ * * * * 0 0 0 0 0 0 0 0 0 ] [ ] L6 [ * * * * * 0 1/2 L7 1/2 L8 1/2 L9 0 0 0 0 ] [ ] L7 [ * * * * * * 0 0 0 L7 L8 0 L9 ] [ ] L8 [ * * * * * * * 0 0 L8 - L7 - L9 0 ] [ ] L9 [ * * * * * * * * 0 L9 0 L8 - L7 ] [ ] L10[ * * * * * * * * * 0 0 0 0 ] [ ] L11[ * * * * * * * * * * 0 - L13 L12 ] [ ] L12[ * * * * * * * * * * * 0 - L11 ] [ ] L13[ * * * * * * * * * * * * 0 ] The i-j th entry in the above table corresponds to [L(i), L(j)]. In summary the exact symmetry generators for the mhd equations were obtained by applying fully automatic procedures. The above results are valid for k # 0,1,3/2, 2. We have separately run the cases K=1, 3/2, 2 and found that the standard forms were just the same as substituting the given values of K in the standard form above - and that we get the same size of group and so these cases have the same generators as above. When we ran the program on the case K=0 we found a 14-parameter invariance group, with the same generators above and an extra generator translation invariance in the pressure: d/dp. 9.2 A SMALL NASTY EXAMPLE: THE NONLINEAR TELEGRAPH SYSTEMS The nonlinear telegraph system produces a determining system which appears quite simple and easy to deal with, but leads to memory explosion in the default settings for standard_form. We believe that this system is a tough challenge for any automatic PDE simplification or integration program (and the solutions are interesting mathematically and physically! [Reid91b]). The nonlinear telegraph system is Vt = Ux , Vx = C(U) Ut + B(U) . where C and B are arbitary functions of U. The determining system for it is written in repeated index notation, as it would be entered into standard form, was derived in Section 4.3 and is given below: sys:= [V(6)*V(3,2)+V(5,3)*V(3)+V(5)^2*V(1,4)+V(5)*(V(1,1)-V(4,4))-V(4,1)=0, V(6)*(V(3,4)-V(1,2))+V(5)*(V(2,4)+V(1,3))-V(4,3)+V(2,1)=0, V(6,3)*V(3)+2*V(5)*V(6)*V(1,4)+V(6)*(V(3,3)+V(1,1)-V(2,2)-V(4,4))=0, V(2,3)-V(6)*V(1,4)=0, V(5)*(V(3,4)+V(1,2))-V(4,2)+V(3,1)=0, V(2,2)+V(3,3)-V(1,1)-V(4,4)=0, V(1,3)-V(2,4)=0, V(6)*(V(3,4)+V(1,2))+V(5)*(V(1,3)-V(2,4))-V(4,3)-V(2,1)=0]; Note that it is a classification problem. In addition to the 4 infinitesimals corresponding to the 2 independent and 2 dependent variables of the telegraph system, there are the additional dependent variables V(5) and V(6), corresponding to the arbitrary functions C and B. They are classification variables, and so are entered in the list class. class:= [[V(5),3],[V(6),3]]; Now we will find the standard form. Note that we are setting the equation control knobs (see 3.3.2) in the strategy list at their ultraconservative setting. > sf:= standard_form(sys, class_eqns=class,strategy=[1.0,1.0,0.9]); The system parameters are as follows: There are 4 independent variables: [x1, x2, x3, x4] There are 4 dependent variables: [V(1), V(2), V(3), V(4)] There are 2 classification variables: [V(5), V(6)] bytes used=589668, alloc=393144, time=6.600 Initial State of System: #OT: 0 #Ez: 0 #NL: 0 #UnC: 8 bytes used=681444, alloc=393144, time=7.666 ******End iteration #1 Time:.2.683 #OT: 0 #Ez: 1 #NL: 0 #UnC: 7 bytes used=1012948, alloc=589716, time=10.333 ******End iteration #2 Time:.2.050 #OT: 0 #Ez: 2 #NL: 0 #UnC: 6 bytes used=1174184, alloc=589716, time=12.466 ******End iteration #3 Time:.3.067 #OT: 0 #Ez: 3 #NL: 0 #UnC: 6 bytes used=1448196, alloc=589716, time=15.633 ******End iteration #4 Time:.3.400 #OT: 0 #Ez: 4 #NL: 0 #UnC: 6 bytes used=1737072, alloc=655240, time=19.083 ******End iteration #5 Time:.3.300 #OT: 0 #Ez: 5 #NL: 0 #UnC: 5 bytes used=1997928, alloc=655240, time=22.416 ******End iteration #6 Time:.5.683 #OT: 0 #Ez: 6 #NL: 0 #UnC: 6 bytes used=2512716, alloc=786288, time=28.266 ******End iteration #7 Time:.5.633 #OT: 0 #Ez: 7 #NL: 0 #UnC: 5 bytes used=5970332, alloc=1376004, time=59.416 ******End iteration #8 Time:.31.317 #OT: 0 #Ez: 8 #NL: 0 #UnC: 7 bytes used=8773096, alloc=1376004, time=83.216 ******End iteration #9 Time:.24.550 #OT: 0 #Ez: 9 #NL: 0 #UnC: 9 bytes used=11655660, alloc=1441528, time=109.000 ******End iteration #10 Time:.21.484 #OT: 0 #Ez: 10 #NL: 0 #UnC: 9 bytes used=15004008, alloc=1572576, time=138.116 ******End iteration #11 Time:.29.966 #OT: 0 #Ez: 11 #NL: 0 #UnC: 10 bytes used=18418052, alloc=1572576, time=170.883 ******End iteration #12 Time:.34.650 #OT: 0 #Ez: 12 #NL: 0 #UnC: 10 bytes used=23125076, alloc=1572576, time=215.566 ******End iteration #13 Time:.45.584 #OT: 0 #Ez: 13 #NL: 0 #UnC: 11 bytes used=39105772, alloc=2555436, time=352.183 ******End iteration #14 Time:.139.534 #OT: 0 #Ez: 14 #NL: 0 #UnC: 14 bytes used=61606716, alloc=2817532, time=538.583 ******End iteration #15 Time:.181.583 #OT: 0 #Ez: 15 #NL: 0 #UnC: 17 bytes used=69920608, alloc=2817532, time=613.500 ******End iteration #16 Time:.74.183 #OT: 1 #Ez: 15 #NL: 0 #UnC: 8 bytes used=71817492, alloc=2817532, time=630.733 ******End iteration #17 Time:.22.134 #OT: 3 #Ez: 14 #NL: 0 #UnC: 2 bytes used=86684956, alloc=2817532, time=747.966 ******End iteration #18 Time:.113.016 #OT: 10 #Ez: 7 #NL: 0 #UnC: 5 bytes used=87367372, alloc=2817532, time=755.000 ******End iteration #19 Time:.5.284 #OT: 8 #Ez: 6 #NL: 0 #UnC: 0 bytes used=87688192, alloc=2817532, time=758.750 bytes used=87689740, alloc=2817532, time=759.050 sf := [V(1, 1) = 0, V(2, 1) = 0, V(4, 1) = 0, V(1, 2) = 0, V(2, 2) = 0, V(4, 2) = 0, V(1, 3) = 0, V(2, 3) = 0, V(4, 3) = 0, V(1, 4) = 0, V(2, 4) = 0, V(4, 4) = 0, V(3) = 0] The output illuminates the process standard_form went through in completing this awesome task. The first 7 iterations are fairly routine short iterations, with standard form adding equations to the Ez list 1 equation at a time (recall that this is because we have set the knobs to their ultra-conservative setting). The next 6 iterations up to the 13th are of moderate length in terms of CPU time, and the UnC list is growing a bit faster than before, showing that as standard_form adds the longer equations from the UnC list, more difficult integrability conditions are arising. But the 14th iteration takes 139.534 seconds of CPU time, and the UnC list grows by an additional 4, as it does at the 15 iteration which takes 181.583 seconds of CPU. The equations arising at these iterations are causing far more work form standard_form than in previous iterations. At the 16th iteration however, a "miracle" happens and we get a collapse in the UnC list from 17 to 8 equations, and at the 17th, it collapses to just 2. Their is one more difficult iteration at #18, which takes 113.016 CPU seconds, but the final system has only 13 one term equations. The final standard form is very simple, but it was a great deal of work to find it. Note that while the output system is quite simple, the pivots, if we were to look at them, even in simplified form, would take up several pages. But with Ian Lisle's ultra sophisticated approach to classification problems given in [Lis92] these pivots can be simplified to a few lines. The generic result above reveals only inspectional symmetries. However by examining the pivots and following special cases where some of the pivots vanish truly nontrivial results are obtained [Reid91b]. For example a nontrivial calculation and result occurs if B=1/(1+U^2), C=1/(1+U^2)^2 > sf1:= standard_form(sys, class_eqns=[[V(5),3],[V(6),3], V(5)=1/(1+x3^2), V(6)=1/(1+x3^2)^2],strategy=[1.75,1.0,0.9]); V(3) sf1 := [V(1, 1) = ----, V(2, 1) = 0, V(3, 1) = 0, V(4, 1) = 0, V(1, 2) = 0, x3 V(3) V(2, 2) = - ----, V(3, 2) = 0, V(4, 2) = 0, V(1, 3) = 0, x3 V(3) V(3) V(2, 3) = - 2 -------------, V(3, 3) = ----, V(4, 3) = 0, 2 2 x3 (1 + x3 ) x3 V(3) V(3) V(1, 4) = - 2 ----, V(2, 4) = 0, V(3, 4) = 0, V(4, 4) = - ----] x3 x3 > id1:= initial_data(sf1, class_eqns=[[V(5),3],[V(6),3], V(5)=1/(1+x3^2), V(6)=1/(1+x3^2)^2]); id1 := [[V(1) = a1, V(2) = a2, V(3) = a3, V(4) = a4], []] Thus there is a nontrivial 4-dimensional Lie symmetry group in this case. The commutator table is > ct1:= fin_com_table(sf1, id1, class_eqns=[[V(5),3],[V(6),3], V(5)=1/(1+x3^2), V(6)=1/(1+x3^2)^2]); [ L(1) ] [ 0 0 ---- 0 ] [ X3 ] [ ] [ L(2) ] [ * 0 - ---- 0 ] ct1 := [ X3 ] [ ] [ L(1) L(4) ] [ * * 0 2 ---- + ---- ] [ X3 X3 ] [ ] [ * * * 0 ] where X3 can be taken to be any fixed nonzero value (e.g. X3 = 1). REFERENCES \documentstyle[12pt]{article} \begin{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % BIBLIOGRAPHY % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{thebibliography}{Abc99} \bibitem{BCGGG89} Bryant, R. L., Chern, S. S., Gardner R. B., Goldschmidt H. L. and Griffiths, P. A. (1991), Exterior Differential Systems, Mathematical Sciences Research Institute Publications {\bf 18} (Springer-Verlag, New York). \bibitem{BlK89} Bluman, G.W. and S. Kumei. 1989. {\em Symmetries and Differential Equations.} Springer, New York. \bibitem{Boch89} Bocharov, A. V. and Bronstein, M. L. (1989), ``Efficiently implementing two methods of the geometrical theory of differential equations'', {\it Acta. Appl. Math.} {\bf 16} 143-16. \bibitem{CamDevFee92} {\sc J. Carminati, J.S. Devitt} and {\sc G.J. Fee},\, {\em Isogroups of differential equations using algebraic computing}, J. Sym. Comp. {\bf 14} (1992) 103-120. \bibitem{Carra87} Carr\`{a}-Ferro, G. 1987, ``Gr\"{o}bner bases and differential ideals," {\it Proc. AAECC5}, Menorca, Spain, Springer-Verlag, 129-140. \bibitem{CarraSitt93} Carr\`a Ferro, G. and Sit, W. Y. (1993), ``On Term-Orderings and Rankings", to appear. \bibitem{Cart45} E. Cartan, {\it Les syst\`{e}mes diff\'{e}rentiels ext\'{e}rieurs et leurs applications g\'{e}om\'{e}triques}, (Paris, Hermann, 1945). \bibitem{Champ91} Champagne, B., Hereman, W. and Winternitz, P. (1991), ``The computer calculation of Lie point symmetries of large systems of differential equations," {\it Computer Physics Communications} {\bf 66}, 319-340. \bibitem{Cha90} Char, B.W., K.O. Geddes, G.H. Gonnet, M.B. Monagan, and S.M. Watt. 1990. (3rd ed.) {\em {\sc Maple} Reference Manual.} Watcom, Waterloo. \bibitem{Fed86} Fedorova, R.N. and V.V. Kornyak. 1986. Determination of Lie B\"acklund symmetries of differential equations using {\sc FORMAC}. {\em Comput. Phys. Comm.} {\bf 39}: 93--103. \bibitem{HaE71} Harrison, B.K. and F.B. Estabrook. 1971. Geometric approach to invariance groups and solution of partial differential equations. {\em J.~Math. Phys.} {\bf 12}: 653--666. \bibitem{HartTker91} Hartley, D. and Tucker, R. W. 1991. {\it A Constructive Implementation of the Cartan-Kahler Theory of Exterior Differential Systems}. J. Symb. Comp. {\bf 12}, 655-657. \bibitem{Here93} Hereman, W. 1993. Review of Symbolic Software for the Computation of Lie Symmetries of Differential Equations. Preprint. \bibitem{Hick93} Hickman, M. 1993. The use of MAPLE in the Search for Symmetries. Preprint. \bibitem{Jan20} Janet, M. 1920. Sur les syst\`emes d'\'equations aux d\'eriv\'ees partielles. {\em J. de math\/} {\bf 3}: 65--151. \bibitem{Ker87} Kersten, P.H.M. 1987. {\em Infinitesimal symmetries: a computational approach.} Centrum voor Wiskunde en Informatica, Amsterdam. \bibitem{Kol73} Kolchin, E. R. (1973), {\it Differential Algebra and Algebraic Groups}, Academic Press, New York. \bibitem{Lis92} Lisle, I.G. 1992. {\em Equivalence Transformations for Classes of Differential Equations\/.} Ph.D Thesis, University of British Columbia. \bibitem{Mans91} Mansfield, E. 1991. {\em Differential Gr\"{o}bner Bases}, Ph.D Thesis, University of Sydney. \bibitem{Olliv91} Ollivier, F. (1991), ``Standard bases of differential ideals'', {\it Lecture Notes in Comp. Sci.}, {\bf 508}, 304-321. \bibitem{Olv86} Olver, P.J. 1986. {\em Applications of Lie Groups to Differential Equations.} Springer, New York. \bibitem{Ovs82} Ovsiannikov, L.V. 1982. {\em Group Analysis of Differential Equations.\/} Academic Press, New York. \bibitem{Pank89} E. V. Pankrat'ev, Acta Applicandae Mathematicae {\bf 16}, 167 (1989). \bibitem{Pomm78} Pommaret, J. F. 1978. {\em Systems of Partial Differential Equations and Lie Pseudogroups} (Gordon and Breach science publishers, Inc.) \bibitem{Pomm90} Pommaret, J. F. 1990. {\em Equations aux d\'{e}riv\'{e}es partielles et th\'{e}orie des groupes} (Institute National de Recherche en Informatique et en Automatique, Roquencourt). \bibitem{Reid90} Reid, G.J. 1990. A triangularization algorithm which determines the Lie symmetry algebra of any system of PDEs. {\em J. Phys.} {\bf A23}: L853--859. \bibitem{Reid91a} Reid, G.J. 1991. Algorithms for reducing a system of PDEs to standard form, determining the dimension of its solution space and calculating its Taylor series solution. {\em Euro. J. Appl. Maths.\/} {\bf 2}: 293--318. \bibitem{Reid91b} Reid, G.J. 1991. Finding abstract Lie symmetry algebras of differential equations without integrating determining equations. {\em Euro. J. Appl. Maths.\/} {\bf 2}: 319--340. \bibitem{Reid91c} Reid, G. J. and Boulton, A. (1991c), ``Reduction of systems of differential equations to standard form and their integration using directed graphs," {\it Proceedings of ISSAC '91}, ACM Press, New York, 308-312. \bibitem{Reid92a} Reid, G.J., I.G. Lisle, A. Boulton and A. D. Wittkopf. 1992. Algorithmic determination of commutation relations for Lie symmetry algebras of PDEs. In {\em Proc. ISSAC~'92 }. ACM Press, Berkeley. \bibitem{Reid92b} Reid, G. J. and McKinnon, D. K. (1992b), {\it Solving systems of linear PDEs in their coefficient field by recursively decoupling and solving ODEs}, preprint. \bibitem{Riq10} Riquier, C. 1910. {\em Les Syst\`{e}mes d'\'{E}quations aux D\'{e}riv\'{e}es Partielles.\/} Gauthier-Villars, Paris. \bibitem{Ritt50} Ritt, J. F. (1950), {\it Differential Algebra}, (Dover). \bibitem{Rust93} Rust, C. 1993. On The Classification of Rankings of Partial Derivatives. Preprint. \bibitem{Schu92} Sch\"{u}, J., W.M. Seiler, and J. Calmet. 1992. Algorithmic Methods for Lie Pseudogroups. In {\em Proc. Modern Group Analysis.\/} Acireale. \bibitem{Sch84} Schwarz, F. 1984. The Riquier-Janet theory and its application to nonlinear evolution equations. {\em Physica\/} {\bf D11}: 243--251. \bibitem{Sch85} Schwarz, F. 1985. Automatically determining symmetries of differential equations. {\em Computing.} {\bf 34}: 91--106. \bibitem{Sch92a} Schwarz, F. 1992. An Algorithm for Determining the Size of Symmetry Groups.{\em Computing.} {\bf 49}: 95--115. \bibitem{Sch92b} Schwarz, F. 1992. Reduction and completion algorithms for Partial Differential Equations. In {\em Proc. ISSAC '92}. 49--56. ACM Press, Berkeley. \bibitem{Spenc69} Spencer, D. C. 1969. {\it Overdetermined systems of linear partial differential equations}. Bull. of the Amer. Math. Soc. {\bf 75}, 179-239. \bibitem{Thom29} Thomas, J.M. 1929. Riquier's existence theorems. {\em Annals of Maths.\/} {\bf 30}: 285--310. \bibitem{Top} Topunov, L. 1989. Reducing systems of linear differential equations to a passive form. {\em Acta Appl. Math.\/} {\bf 16}: 191--206. \bibitem{Tress94} Tresse, A. 1894. Sur les invariants diff\'{e}rentiels des groupes continus de transformations. {\em Acta Mathematica\/} {\bf 18}: 1--88. (English translation I. Lisle 1989, available from the author) \bibitem{Wolf89} Wolf, T. (1989), {\it Proc. EUROCAL 87, Lect. Notes in Comp. Sci.} {\bf 378}, 479. \bibitem{Wolf91} Wolf, T. (1991), ``The symbolic integration of exact PDEs," preprint. \bibitem{Wolf92} Wolf, T. (1992), ``CRACK Reference Manual", REDUCE Network Library. \bibitem{Weisp93} Weispfenning, V. (1993), ``Differential-Term Orders'', to appear in {\it Proc. of ISAAC 1993} (Kiev, acm press). \end{thebibliography} \end{document} INDEX OF FUNCTIONS AND OPTIONS This is a preliminary effort at making an index (it needs to be completed). To use this index, use your editor's search capability to find occurrences of the given expression. FUNCTIONS standard_form initial_data taylor_exp finite_generators lie_symm polyd Clearderivdep sfdiff make_vset shortenpiv map2rep rep2map man2rep rep2man transform readstored max2rep OPTIONS divpiv class_eqns usereasy strategy storeall store little big biglittle PACKAGES sfprogs convprogs maxconv graphDE symmetry liesymm symgrp.max OTHERS _pivs _incons