Artifact 399ba2ce07f5f40fc40e8f24b55ed422c789e9d1809bd1a883ebeb5368e466b8:
- Executable file
r37/packages/crack/liepde.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 52453) [annotate] [blame] [check-ins using] [more...]
%******************************************************************** % * % The program LIEPDE for computing point-, contact- and higher * % order symmetries of individual ODEs/PDEs or systems of ODEs/PDEs * % * % Author: Thomas Wolf * % Date: 20.July 1996 * % * % For details of how to use LIEPDE see the file LIEPDE.TXT or the * % header of the procedure LIEPDE below. * % * %******************************************************************** create!-package('(liepde), nil); % FJW % FJW: Load support packages, but not when compiling: !#if (getd 'packages_to_load) packages_to_load crack; !#else % for REDUCE 3.6 apply1('load!-package, 'crack); !#endif symbolic procedure adjoin(elm, set); %% Add Element to SET if it is not already a member. %% FJW: Defined and flagged lose in PSL only. %% FJW: This implementation from the PSL manual. if member(elm, set) then set else elm . set$ symbolic fluid '(print_ logoprint_ nfct_ fname_ adjust_fnc proc_list_ prelim_ individual_)$ lisp << prelim_:=nil$ individual_:=nil >>$ %---------------------------- % Nachfolgende 2 Prozeduren aus crack1.48: symbolic procedure diffdeg(p,v)$ % liefert Ordnung der Ableitung von p nach v$ % p Liste Varible+Ordnung der Ableitung, v Variable (Atom) if null p then 0 % alle Variable bearbeitet ? else if car p=v then % v naechste Variable ? if cdr p then if numberp(cadr p) then cadr p % folgt eine Zahl ? else 1 else 1 else diffdeg(cdr p,v)$ % weitersuchen symbolic procedure ldiff1(l,v)$ % liefert Liste der Ordnungen der Ableitungen nach den Variablen aus v % l Liste (Variable + Ordnung)$ v Liste der Variablen if null v then nil % alle Variable abgearbeitet ? else cons(diffdeg(l,car v),ldiff1(l,cdr v))$ % Ordnung der Ableitung nach % erster Variable anhaengen %---------------------------- algebraic procedure equ_to_expr(a)$ % converts an equation into an expression begin scalar lde; return if a=nil then a else <<lisp(lde:=reval algebraic a); if lisp(atom lde) then a else num if lisp(car lde = 'EQUAL) then lhs a - rhs a else a >> end$ % of equ_to_expr %---------------------------- algebraic procedure maklist(ex)$ % making a list out of an expression if not already if lisp(atom algebraic ex) then {ex} else if lisp(car algebraic ex neq 'LIST) then ex:={ex} else ex$ %******************************************************************** module pdesymmetry$ %******************************************************************** % Routines for finding Symmetries of single or systems of ODEs/PDEs % Author: Thomas Wolf % July 1996 % The algebraic operator NPRIMITIVE returns the % numerically-primitive part of any expression. % It is defined as a simpfn in EZGCD. %% FJW: This procedure is already defined (differently but %% equivalently!) in module crsimpso: %% symbolic operator ncontent$ %% symbolic procedure ncontent p$ %% % Return numeric content of expression p %% % based on simpnprimitive in ezgcd. %% << p := simp!* p; %% if polyzerop(numr p) then 0 else %% mk!*sq(numeric!-content numr p ./ numeric!-content denr p) %% >>$ symbolic operator totdeg$ symbolic procedure totdeg(p,f)$ % Ordnung (total) der hoechsten Ableitung von f im Ausdruck p eval(cons('PLUS,ldiff1(car ldifftot(reval p,reval f),fctargs reval f)))$ symbolic procedure diffreltot(p,q,v)$ % liefert komplizierteren Differentialausdruck$ if diffreltotp(p,q,v) then q else p$ symbolic procedure diffreltotp(p,q,v)$ % liefert t, falls p einfacherer Differentialausdruck, sonst nil % p, q Paare (liste.power), v Liste der Variablen % liste Liste aus Var. und Ordn. der Ableit. in Diff.ausdr., % power Potenz des Differentialausdrucks begin scalar n,m$ m:=eval(cons('PLUS,ldiff1(car p,v)))$ n:=eval(cons('PLUS,ldiff1(car q,v)))$ return if m<n then t else if n<m then nil else diffrelp(p,q,v)$ end$ symbolic procedure ldifftot(p,f)$ % leading derivative total degree ordering % liefert Liste der Variablen + Ordnungen mit Potenz % p Ausdruck in LISP - Notation, f Funktion ldifftot1(p,f,fctargs f)$ symbolic procedure ldifftot1(p,f,vl)$ % liefert Liste der Variablen + Ordnungen mit Potenz % p Ausdruck in LISP - Notation, f Funktion, lv Variablenliste begin scalar a$ a:=cons(nil,0)$ if not atom p then if member(car p,list('EXPT,'PLUS,'MINUS,'TIMES, 'QUOTIENT,'DF,'EQUAL)) then % erlaubte Funktionen <<if (car p='PLUS) or (car p='TIMES) or (car p='QUOTIENT) or (car p='EQUAL) then <<p:=cdr p$ while p do <<a:=diffreltot(ldifftot1(car p,f,vl),a,vl)$ p:=cdr p>> >> else if car p='MINUS then a:=ldifftot1(cadr p,f,vl) else if car p='EXPT then % Exponent if numberp caddr p then <<a:=ldifftot1(cadr p,f,vl)$ a:=cons(car a,times(caddr p,cdr a))>> else a:=cons(nil,0) % Poetenz aus Basis wird mit % Potenz multipliziert else if car p='DF then % Ableitung if cadr p=f then a:=cons(cddr p,1) % f wird differenziert? else a:=cons(nil,0)>> % sonst Konstante bzgl. f else if p=f then a:=cons(nil,1) % Funktion selbst else a:=cons(nil,0) % alle uebrigen Fkt. werden else if p=f then a:=cons(nil,1)$ % wie Konstante behandelt return a end$ %--------------------- % Bei jedem totdf2-Aufruf pruefen, ob evtl. kuerzere dylist reicht % evtl die combidiff-Kette und combi nicht mit in dylist, sond. erst in % prolong jedesmal frisch generieren. %symbolic operator desort$ algebraic procedure nextdy(revx,xlist,dy)$ % generates all first order derivatives of lhs dy % revx = reverse xlist; xlist is the list of variables; % dy the old derivative begin scalar x,n,ldy,rdy,ldyx,sublist; x:=first revx; revx:=rest revx; sublist:={}; ldy:=lhs dy; rdy:=rhs dy; while lisp(not member(prepsq simp!* algebraic x, prepsq simp!* algebraic ldy)) and (revx neq {}) do <<x:=first revx; revx:=rest revx>>; n:=length xlist; if revx neq {} then % dy is not the function itself while first xlist neq x do xlist:=rest xlist; xlist:=reverse xlist; % New higher derivatives while xlist neq {} do <<x:=first xlist; ldyx:=df(ldy,x); sublist:=cons((lisp reval algebraic ldyx)= mkid(mkid(rdy,!`),n), sublist); n:=n-1; xlist:=rest xlist >>; return sublist end$ %--------------------- algebraic procedure subdif1(xlist,ylist,ordr)$ % A list of lists of derivatives of one order for all functions begin scalar allsub,revx,i,el,oldsub,newsub; revx:=reverse xlist; allsub:={}; oldsub:= for each y in ylist collect y=y; for i:=1:ordr do % i is the order of next substitutions <<oldsub:=for each el in oldsub join nextdy(revx,xlist,el); allsub:=cons(oldsub,allsub) >>; return allsub end$ %--------------------- symbolic procedure combidif(s)$ % extracts the list of derivatives from s: % u`1`1`2 --> (u, 1, 1, 2) begin scalar temp,ans,no,n1; s:=reval s; % to guarantee s is in true prefix form temp:=reverse explode s; while not null temp do <<n1:=<<no:=nil; while (not null temp) and (not eqcar(temp,'!`)) do <<no:=car temp . no;temp:=cdr temp>>; compress no >>; if (not fixp n1) then n1:=intern n1; ans:=n1 . ans; if eqcar(temp,'!`) then <<temp:=cdr temp; temp:=cdr temp>>; >>; return ans end$ %--------------------- symbolic procedure comparedif1(u1l,u2l)$ % u1l, u2l are lists of indicees of differentiation variables % checks whether u2l has more or at least equally many 1's, 2's, ... % contained as u1l. % returns a list of 1's, 2's, ... which are in excess in u2l % compared with u1l. The returned value is 0 if both are identical begin scalar ul; if u2l=nil then if u1l neq nil then return nil else return 0 else if u1l=nil then return u2l else % both are non-nil if car u1l < car u2l then return nil else if car u1l = car u2l then return comparedif1(cdr u1l,cdr u2l) else << ul:=comparedif1(u1l,cdr u2l); return if not ul then nil else if zerop ul then list car u2l else cons(car u2l,ul) >> end$ % of comparedif1 %--------------------- symbolic procedure comparedif2(u1,u1list,du2)$ % checks whether du2 is a derivative of u1 differentiated % wrt. u1list begin scalar u2l; u2l:=combidif(du2)$ % u2l=(u2, 1, 1, ..) if car u2l neq u1 then return nil else return comparedif1(u1list, cdr u2l) end$ % of comparedif2 %--------------------- symbolic procedure comparedif3(du1,u2,u2list)$ % checks whether u2 differentiated wrt. u2list % is a derivative of du1 begin scalar u1l; u1l:=combidif(du1)$ % u1l=(u1, 1, 1, ..) if car u1l neq u2 then return nil else return comparedif1(cdr u1l, u2list) end$ % of comparedif3 %--------------------- symbolic procedure listdifdif1(du1,deplist)$ % lists all elements of deplist which are *not* derivatives % of du1 begin scalar u1,u1list,res,h$ h:=combidif(du1); u1:=car h; u1list:=cdr h; for each h in deplist do if not comparedif2(u1,u1list,h) then res:=cons(h,res); return res end$ % of listdifdif1 %--------------------- symbolic operator dif$ symbolic procedure dif(s,n)$ % e.g.: dif(fnc!`1!`3!`3!`4, 3) --> fnc!`1!`3!`3!`3!`4 begin scalar temp,ans,no,n1,n2,done; s:=reval s; % to guarantee s is in true prefix form temp:=reverse explode s; n2:=reval n; n2:=explode n2; while (not null temp) and (not done) do <<n1:=<<no:=nil; while (not null temp) and (not eqcar(temp,'!`)) do <<no:=car temp . no;temp:=cdr temp>>; compress no >>; if (not fixp n1) or ((fixp n1) and (n1 leq n)) then <<ans:=nconc(n2,ans); ans:='!` . ans; ans:='!! . ans; done:=t>>; ans:=nconc(no,ans); if eqcar(temp,'!`) then <<ans:='!` . ans; ans:='!! . ans; temp:=cdr temp; temp:=cdr temp>>; >>; return intern compress nconc(reverse temp,ans); end$ %--------------------- %symbolic procedure orderofderiv(du)$ %if atom du then (length combidif(du))-1 % else nil$ %--------------------- symbolic procedure mergedepli(li1,li2)$ % li1,li2 are lists of lists % mergedepli merges the sublists to make one list of lists begin scalar newdep; while li1 and li2 do << newdep:=union(car li1, car li2) . newdep; li1:=cdr li1; li2:=cdr li2 >>; return if li1 then nconc(reversip newdep,li1) else if li2 then nconc(reversip newdep,li2) else reversip newdep end$ %--------------------- symbolic procedure adddepli(ex,revdylist)$ begin scalar a,b,c,d; for each a in revdylist do << c:=nil; for each b in a do if not my_freeof(ex,b) then c:=b . c; if c or d then d:=c . d; >>; return list(ex,d) end$ %--------------------- symbolic procedure add_xi_eta_depli(xilist,etalist,revdylist)$ begin scalar e1,g,h$ for e1:=1:(length xilist) do << g:=nth(xilist,e1); h:=pnth(g,4); rplaca(h,cadr adddepli(car g,revdylist)) >>; for e1:=1:(length etalist) do << g:=nth(etalist,e1); h:=pnth(g,3); rplaca(h,cadr adddepli(car g,revdylist)) >> end$ %--------------------- symbolic procedure subtest(uik,sb,xlist,ordok,subordinc)$ begin scalar el5,el6,el7,el8,el9,el10,sbc$ el5:=combidif(uik); el6:=car el5; el5:=cdr el5; % el6..function name, el5..var.list el7:=nil; el8:=100; el9:=nil; sbc:=sb; while sbc and ((caaar sbc neq el6) or (0 neq <<el7:=comparedif1(cdaar sbc,el5); if el7 and (not zerop el7) and (length(el7)<el8) then <<el8 :=length el7; el9 :=el7; el10:=car sbc>> else el7 >>) ) do sbc:=cdr sbc; return if sbc then (cadar sbc . caddar sbc) % simple substitution else if el9 then << % uik is total deriv of car el10 wrt el9 uik:=cadr el10 . caddr el10; % car uik becomes the expr. to replace the former uik while el9 do << uik:=totdf3(car uik, cdr uik, nth(xlist, car el9), car el9, sb, xlist, ordok, subordinc); el9:=cdr el9 >>; uik >> else % no substitution nil end$ %--------------------- symbolic procedure totdf3(s,depli,x,n,sb,xlist,ordok,subordinc)$ % s is the expression to be differentiated wrt. x which is the nth % variable in xlist % depli is a list of lists, each being the list of jet-variables % of order 0,1,2,..., such that s=s(xlist,depli), but % as little as possible jet-variables in depli % xlist, depli are lisp lists, i.e. without 'LIST % - totdf3 calculates total derivative of s(xlist,depli) w.r.t. x which % is the n'th variable, it returns (df(s,x), newdepli) % - totdf3 drops jet-variables on which s does not depend % - totdf3 automatically does substitutions using the list sb which % is updated if prolongations of substitutions are calculated, % i.e. sb is destructively changed!! % - structure of sb: lisp list of lisp lists: % ((to_be_replaced_jet_var_name, to_be_replaced_jet_var_deriv_1,..), % subst_expr_in_jet_space_coord, list_of_jet_vars_in_subst_expr) % - subordinc is a number by how much the order may increase due to % substitutions sb. % - ordok is the lowest order which must be accurate. If ordok>0 and % s is of lower order than ordok then from depli only derivatives % of order ordok-1-subordinc to ordok-1 are used. begin scalar tdf,el1,el2,el3,el4,el5,newdepli, newdy,dy,ddy,s; newdepli:=nil; % the new dependence list newdy:=nil; % the new dep.list due to chain rule ddy:=nil; % ddy .. derivatives of jet-variables % resulting from diff. of lower order %--- Should only terms in the result be acurate that include %--- derivatives of order>=ordok? if ordok>0 then << tdf:=simp!* 0; depli:=copy depli; el2:=length depli; if el2<(ordok-subordinc) then depli:=nil else for el1:=1:(ordok-1-subordinc) do << dy:=pnth(depli,el1); rplaca(dy,nil); >> >> else tdf:=simpdf {s,x}; %--- The differentiations wrt. u-derivatives for each el1 in depli do % for each order do <<dy:=union(ddy,el1); ddy:=nil;% dy .. occuring jet-var. of this order while el1 do <<el2:=car el1; el1:=cdr el1;% el2 is one jet-variable of this order el3:=simpdf {s,el2}; if zerop el3 then dy:=delete(el2,dy) else << el4:=dif(el2,n); % el4=df(el2,x) %----- Is el4 to be substituted by sb? if el5:=subtest(el4,sb,xlist,ordok,subordinc) then << el4:=car el5; newdepli:=mergedepli(newdepli,cdr el5) >> else ddy:=el4 . ddy; tdf:=addsq(tdf, multsq(simp!* el4, el3)) >> >>; newdy:=dy . newdy >>; if ddy then newdy:=ddy . newdy; newdepli:=mergedepli(reversip newdy,newdepli); % possibly drop at the end return (prepsq tdf . newdepli) end$ % of totdf3 %--------------------- symbolic procedure joinsublists(a)$ % It is assumed, a is either nil or a list of lists or nils which % have to be joined if null a then nil else append(car a,joinsublists(cdr a))$ %--------------------- symbolic procedure depnd(y,xlist)$ for each xx in xlist do for each x in xx do depend y,x$ %--------------------- algebraic procedure transeq(eqn,xlist,ylist,sb)$ <<for each el1 in sb do eqn:=sub(el1,eqn); for each el1 in ylist do for each el2 in xlist do nodepend el1,el2; eqn>>$ %--------------------- symbolic operator drop$ symbolic procedure drop(a,vl)$ % liefert summe aller terme aus a, die von elementen von vl abhaengen begin scalar b$ if not((pairp a) and (car a='PLUS)) then b:=a else <<vl:=cdr vl; % because vl is an algebraic list for each c in cdr a do if not freeoflist(c,vl) then b:=cons(c,b)$ if b then b:=cons('PLUS,reverse b)>>$ return b$ end$ %--------------------- symbolic procedure etamn(u,indxlist,xilist,etalist, ordok,truesub,subordinc,xlist)$ % determines etamn recursively % At the end, ulist= list of df(u,i,cdr indxlist) for all i begin scalar etam,x,h1,h2,h3,h4,ulist,el,r,cplist,depli; if (null indxlist) or ((length indxlist)=1) then <<cplist:=etalist; while u neq cadar cplist do cplist:=cdr cplist; etam:=(caar cplist . caddar cplist) . nil; >> else etam:=etamn(u,cdr indxlist,xilist,etalist,ordok,truesub, subordinc,xlist)$ return if null indxlist then etam else << ulist:=nil; x:=cdr nth(xilist,car indxlist); % e.g. x := (v3,3,dylist) r:=if zerop caar etam then simp!* <<depli:=nil;0>> else simp!* << h2:=totdf3(caar etam,cdar etam,car x,cadr x,truesub,xlist, ordok,subordinc)$ depli:=cdr h2; car h2 >>; etam:=cdr etam; % = reverse ulist cplist:=xilist; h3:=nil; while cplist do <<el:=car cplist; % e.g. el=xi_z cplist:=cdr cplist; if (length indxlist)=1 then h1:=dif(u,caddr el) else << h1:=dif(car etam,cadr indxlist); % e.g. h1:=u!`i!`n etam:=cdr etam; >>; ulist:=h1 . ulist; if not zerop car el then << %--- substitution of h1? if h4:=subtest(h1,truesub,xlist,ordok,subordinc) then h1:=car h4; r:=subtrsq(r, multsq(simp!* h1, simp!* <<h2:=totdf3(car el,cadddr el,car x, cadr x,truesub,xlist,0,0)$ if zerop car h2 then 0 else <<if h4 then depli:=mergedepli(depli,cdr h4) else h3:=h1 . h3; depli:=mergedepli(depli,cdr h2); car h2 >> >> ) ); >> >>; if h3 then << h3:=list h3; for h2:=1:(length indxlist) do h3:=nil . h3; depli:=mergedepli(depli,h3); >>; % (if not full then drop(r,'LIST . car revdylist) else r) . % (reverse ulist) (prepsq r . depli) . (reverse ulist) >> end$ % of etamn %--------------------- symbolic procedure prolong(uik,xilist,etalist,ordok,truesub,subordinc, xlist)$ begin scalar h; h:=combidif(uik); h:=car etamn(car h,cdr h,xilist,etalist,ordok,truesub, subordinc,xlist)$ return (simp!* car h) . cdr h end$ % of prolong %--------------------- symbolic procedure callcrack(!*time,cpu,gc,lietrace_,symcon, flist,vl,xilist,etalist)$ begin scalar g,h; if !*time then <<terpri()$ write "time to formulate conditions: ", time() - cpu, " ms GC time : ", gctime() - gc," ms"$ >>; if lietrace_ then algebraic << write"Symmetry conditions before CRACK: "; write lisp ('LIST . symcon); >>; h:=crack('LIST . symcon,list('LIST),'LIST . flist,'LIST . vl); if h neq list('LIST) then <<h:=cadr h; symcon:=cdadr h; for each g in cdaddr h do << xilist :=subst(caddr g, cadr g, xilist); etalist:=subst(caddr g, cadr g, etalist); %--> Erkennung von 'e, 'x siehe: % h:=intern car explode cadr e1; %write"h=",h;terpri()$ % if (h='x) or (h='X) then % xilist :=subst(caddr e1, cadr e1, xilist) else % if (h='e) or (h='E) or (h="e") or (h="E") then % etalist:=subst(caddr e1, cadr e1, etalist) else % rederr("One ansatz does not specify XI_ nor ETA_.") >>; if lietrace_ then << write"symcon nachher: ",symcon; write"xilist=", xilist; write"etalist=", etalist; >>; flist:=cdr reval cadddr h; if print_ then <<terpri()$terpri()$ write"Remaining free functions after the last CRACK-run:"; terpri()$ fctprint flist;terpri()$terpri()>>; >>; return list(symcon,xilist,etalist,flist) end$ % of callcrack %--------------------- symbolic operator liepde$ symbolic procedure liepde(problem,symtype,flist)$ comment problem: {{eq1,eq2, ...}, % equations { y1, y2, ...}, % functions { x1, x2, ...} } % variables % Equations `eqi' can be given as single differential % expressions which have to vanish or they can be given % in the form df(..,..[,..]) = .. . % If the equations are given as single differential % expressions then the program will try to bring it into % the `solved form' ..=.. automatically. % The solved forms (either from input or generated within % LIEPDE) will be used for substitutions, to find % all symmetries satisfied by solutions of the equations. % Sufficient conditions for this procedure to be correct, % i.e. to get *all* symmetries of the specified type on the % solution space are: % - There are equally many equations and functions. % - Each function is used once for a substitution and % each equation is used once for a substitution. % - All functions differentiated on the left hand sides % (lhs) depend on all variables. % - In no equation df(..,..[,..]) = .. does the right hand % side contain the derivative on the lhs nor any % derivative of it. % - No equation does contain a lhs or the derivative % of a lhs of another equation. % These conditions are checked in LIEPDE and execution % is stoped if they are not satisfied, i.e. if the input % was not correct, or if the program was not able to % write the input expressions properly the solved % form ..=.. . One then should find for each function % one derivative which occurs linearly in one equation. % The chosen derivatives should be as high as possible, % at least there must no derivative of them occur in any % equation. An easy way to get the equations in the % desired form is to use % FIRST SOLVE({eq1,eq2,...},{list of derivatives}) % NOTE that to improve efficiency it is advisable not to % express lower order derivatives on the left hand side % through higher order derivatives on the right hand side. % SEE also the implications on completeness for the % determination of generalized symmetries with % characteristic functions of a given order, as described % below and the two examples with the Burgers equation. symtype: {"point"} % for point symmetries {"contact"} % for contact symmetries, is only % applicable if only one function, % only one equation of order>1 {"general",order} % for generalized symmetries of % order `order' which is an integer>0 % NOTE: Characteristic functions of % generalized symmetries (i.e. the % eta_.. if xi_..=0) are equivalent % if they are equal on the solution % manifold. Therefore all dependencies % of characteristic functions on % the substituted derivatives and their % derivatives are dropped. This has the % consequence that if, e.g. for the heat % equation df(u,t)=df(u,x,2), df(u,t) is % substituted by df(u,x,2) then % {"general",2) would not include % characteristic functions depending % on df(u,t,x), or df(u,x,3). THEREFORE: % If you want to find all symmetries up % to a given order then % - either avoid substituting lower % order derivatives by expressions % involving higher derivatives, or, % - go up in the order specified in % symtype. % % Example: % % depend u,t,x % liepde({{df(u,t)=df(u,x,2)+df(u,x)**2}, % {u},{t,x}}, % {"general",3},{}) % % will give 10 symmetries + one infinite % family of symmetries whereas % % liepde({{df(u,x,2)=df(u,t)-df(u,x)**2}, % {u},{t,x}}, % {"general",3},{}) % % will give 28 symmetries + one infinite % family of symmetries. {xi!_x1 =..., ... eta!_y3=... } % - An ansatz must specify all xi!_.. for % all indep. variables and all eta!_.. % for all dep. variables in terms of % differential expressions which may % involve unknown functions/constants. % The dependencies of the unknown % functions have to declared earlier % using the DEPEND command. % - If the ansatz should consist of the % characteristic functions then set % all xi!_..=0 and assign the charac- % teristic functions to the eta!_.. . flist: {- all parameters and functions in the equations which are to be determined, such that there exist symmetries, - if an ansatz has been made in symtype then flist contains all unknown functions and constants in xi!_.. and eta!_..} Further comments: The syntax of input is the usual REDUCE syntax. For example, the derivative of y3 wrt x1 once and x2 twice would be df(y3,x1,x2,2). --> One exception: If in the equations or in the ansatz the dependence of a free function F on a derivative, like df(y3,x1,x2,2) shall be declared then the declaration has to have the form: DEPEND F, Y3!`1!`2!`2 - the ! has to preceede each special character, like `, - `i stands for the derivative with respect to the i'th variable in the list of variables (the third list in the problem above) If the flag individual_ is t then conditions are investigated for each equation of a system of equations at first individually before conditions resulting from other equations are added. If the flag prelim_ is t then preliminary conditions for equations of higher than 1st order are formulated and investigated before the full condition is formulated and investigated by CRACK. If the REDUCE switch TIME is set on with ON TIME then times for the individual steps in the calculation are shown. Further switches and parameters which can be changed to affect the output and performance of CRACK in solving the symmetry conditions are listed in the file CRINIT.RED. ; begin scalar cpu, gc, lietrace_, oldadj, eqlist, ylist, xlist, pointp, contactp, generalp, ansatzp, symord, e1, e2, ordr, sb, dylist, revdylist, xi, eta, eqordr, eqordrcop, no, eqcopy1, truesub, deplist, xilist, etalist, dycopy, freelist, eqlen, dylen, truesubno, minordr, n1, n2, n3, n4, n, h, jetord, allsub, subdy, lhslist, symcon, subordinc, eqn, depli, vl, occli, revdycopy, subordinclist, xicop, etacop, flcop; cpu:=time()$ gc:=gctime()$ lietrace_:=nil; oldadj:=adjust_fnc; adjust_fnc:=nil; %--------- extracting input data eqlist:= cdr maklist cadr problem; ylist := reval maklist caddr problem; xlist := reval maklist cadddr problem; eqlen:=length eqlist; if 1+eqlen neq length(ylist) then rederr( "Number of equations does not match number of unknown functions."); for each e1 in cdr ylist do for each e2 in cdr xlist do if my_freeof(e1,e2) then rederr( "Not all functions do depend on all variables."); if atom cadr symtype then % default case if cadr symtype = "point" then <<pointp :=t;symord:=0>> else if cadr symtype = "contact" then <<contactp:=t;symord:=1; if eqlen>1 then rederr( "Contact symmetries only in case of one equation for one function.") >> else if cadr symtype = "general" then <<generalp:=t;symord:=caddr symtype; if (not fixp symord) or (symord<1) then rederr( "The order of the generalized symmetry must be an integer > 0.") >> else rederr("Inconclusive symmetry type.") else << ansatzp:=t; % an ansatz has been made if length(ylist)+length(xlist) neq length(symtype)+1 then rederr("Number of assignments in the ansatz is wrong."); symord:=0; for each e1 in cdr symtype do for each e2 in ylist do <<n:=totdeg(e1,e2); if n>symord then symord:=n>>; if symtype = 0 then pointp :=t else if (symtype = 1) and (length(ylist)=2) then contactp:=t else generalp:=t >>$ if flist then flist:=cdr flist; problem:=0; %---- Are substitutions already given in the input? eqcopy1:=eqlist; while eqcopy1 and (pairp car eqcopy1) and (caar eqcopy1='EQUAL) and (pairp cadar eqcopy1) and (caadar eqcopy1='DF) do eqcopy1:=cdr eqcopy1; if null eqcopy1 then truesub:=eqlist; eqcopy1:=nil; %--------- initial printout if print_ and logoprint_ then <<terpri()$ write "-----------------------------------------------", "---------------------------"$ terpri()$terpri()$ write"This is LIEPDE - a program for calculating infinitesimal", " symmetries"; terpri()$ write "of single differential equations or systems of de's"; >>; terpri();terpri(); if length xlist=2 then write"The ODE" else write"The PDE"; if length ylist>2 then write"-system"; write " under investigation is :";terpri(); % for each e1 in eqlist do algebraic write"0 = ",lisp e1; for each e1 in eqlist do algebraic write lisp e1; terpri()$write "for the function(s) : ";terpri()$ terpri()$fctprint cdr reval ylist; terpri()$terpri(); eqlist:=for each e1 in eqlist collect algebraic equ_to_expr(lisp e1); if eqlen > 1 then eqlist:=desort eqlist; if !*time then <<terpri()$ terpri()$terpri()$ write"=============== Initializations" ; >>; %--------- initializations ordr:=0; for each e1 in eqlist do << h:=0; for each e2 in cdr ylist do <<n:=totdeg(e1,e2); if n>h then h:=n>>; eqordr:=h . eqordr; if h>ordr then ordr:=h >>; eqordr:=reversip eqordr; if ordr>symord then jetord:=ordr else jetord:=symord$ sb:=subdif1(xlist,ylist,jetord)$ eqlist:=cons('LIST,eqlist); if ansatzp then eqlist:=list('LIST,symtype,eqlist); if truesub then eqlist:=list('LIST,cons('LIST,truesub),eqlist); on evallhseqp; eqlist:=transeq(eqlist,xlist,ylist,sb); off evallhseqp; if truesub then <<truesub:=cdadr eqlist;eqlist:=caddr eqlist>>; if ansatzp then <<symtype:=cdadr eqlist;eqlist:=cdaddr eqlist>> else eqlist:=cdr eqlist; ylist:=cdr ylist; xlist:=cdr xlist; if lietrace_ and ansatzp then write"ansatz=",symtype; dylist:=ylist . reverse for each e1 in cdr sb collect for each e2 in cdr e1 collect caddr e2; revdylist:=reverse dylist; % dylist with decreasing order vl:=xlist; for each e1 in dylist do vl:=append(e1,vl); vl:='LIST . vl; if not ansatzp then deplist:=for n:=0:symord collect nth(dylist,n+1); % list of variables the xi_, eta_ depend on xi :=reval algebraic xi!_; eta:=reval algebraic eta!_; n:=0; xilist :=for each e1 in xlist collect <<n:=n+1; if pointp or ansatzp then << h:=mkid(xi,e1); if not ansatzp then << depnd(h,xlist . deplist); flist:=h . flist; depli:=deplist; >> else depli:=nil >> else <<h:=0;depli:=nil>>; {h,e1,n,depli} >>; depli:=if (not ansatzp) and (not generalp) then deplist else nil; n:=0; etalist:=for each e1 in ylist collect <<n:=n+1; h:=mkid(eta,e1); if not ansatzp then << if not generalp then depnd(h,xlist . deplist); % the generalp-case is done below when substitutions are known flist:=h . flist; >>; {h,e1,depli} >>; if ansatzp then << for each e1 in symtype do << xilist :=subst(caddr e1, cadr e1, xilist); etalist:=subst(caddr e1, cadr e1, etalist); %--> Erkennung von 'e, 'x siehe: % h:=intern car explode cadr e1; %write"h=",h;terpri()$ % if (h='x) or (h='X) then % xilist :=subst(caddr e1, cadr e1, xilist) else % if (h='e) or (h='E) or (h="e") or (h="E") then % etalist:=subst(caddr e1, cadr e1, etalist) else % rederr("One ansatz does not specify XI_ nor ETA_.") >>; add_xi_eta_depli(xilist,etalist,revdylist)$ >>; if lietrace_ then write"xilist=",xilist," etalist=",etalist; %---- Determining a substitution list for highest derivatives %---- from eqlist. Substitutions may not be optimal if starting %---- system is not in standard form comment: Counting in how many equations each highest derivative occurs. Those which do not occur allow Stephani-Trick, those which do occur once and there linearly are substituted by that equation. Because one derivative shall be assigned it must be one of the highest derivatives from each equation used, or one such that no other derivative in the equation is a derivative of it. Each equation must be used only once. Each derivative must be substituted by only one equation. At first determining the number of occurences of each highest derivative. The possiblity of substitutions is checked in each total derivative. $ if truesub then << %--- determination of freelist %and occurlist dycopy:=car revdylist; %--- the highest derivatives while dycopy do <<e1:=car dycopy; dycopy:=cdr dycopy; eqcopy1:=eqlist; while eqcopy1 and my_freeof(car eqcopy1,e1) do eqcopy1:=cdr eqcopy1; if null eqcopy1 then freelist :=e1 . freelist %else occurlist:=e1 . occurlist; >> >> else << no:=0; % counter of the following repeat-loop % freelist (and occurlist) are determined % only in the first run eqordrcop:=copy eqordr; % for bookkeeping which equation have been used repeat << no:=no+1; %--- incrementing the loop counter %--- truesubno is the number of substitutions so far found. %--- It is necessary at the end to check whether new substitutions %--- have been found. if null truesub then truesubno:=0 else truesubno:=length truesub; %--- substitutions of equations of minimal order are searched first minordr:=1000; %--- minimal order of the so far unused equations for each e1 in eqordrcop do if (e1 neq 0) and (e1<minordr) then minordr:=e1; dycopy:=copy nth(dylist,minordr+1); %-- all deriv. of order minordr dylen:=length dycopy; allsub:=nil; for n1:=1:dylen do %--- checking all deriv. of order minordr <<e1:=nth(dycopy,n1); %--- e1 is the current candidate %--- here test, whether e1 is not a derivative of a lhs of one %--- of the substitutions car e2 found so far h:=combidif(e1); n:=car h; h:=cdr h; e2:=truesub; while e2 and (null comparedif3(cadar e2,n,h)) do e2:=cdr e2; if null e2 then << n2:=0; %-- number of equations in which the derivative e1 occurs subdy:=nil; for n3:=1:eqlen do if not my_freeof(nth(eqlist,n3),e1) then % here should also be tested whether derivatives of e1 occur % and not just my_freeof %--> <<n2:=n2+1; if nth(eqordrcop,n3)=minordr then %--- equation is not used yet and of the right order <<e2:=cdr algebraic coeff(lisp nth(eqlist,n3),lisp e1); if hipow!*=1 then subdy:=list(n1,n3,list('EQUAL,e1,list('MINUS, list('QUOTIENT, car e2, cadr e2)))) . subdy >> >>; if n2=0 then if no=1 then freelist:=e1 . freelist else else <<%if no=1 then occurlist:=e1 . occurlist; if subdy then if n2=1 then <<h:=car subdy; truesub:=(caddr h) . truesub; n:=pnth(dycopy ,car h);rplaca(n,0); n:=pnth(eqordrcop,cadr h);rplaca(n,0); >> else allsub:=nconc(allsub,subdy); >> >> >>; %---- Taking the remaining known substitutions of highest deriv. h:=subdy:=0; for each h in allsub do if (nth(dycopy , car h) neq 0) and (nth(eqordrcop,cadr h) neq 0) then <<truesub:=(caddr h) . truesub; n:=pnth(dycopy ,car h);rplaca(n,0); n:=pnth(eqordrcop,cadr h);rplaca(n,0); >>; >> until (truesub and (length(truesub)=eqlen)) % complete or (truesubno=length(truesub))$ % no progress allsub:=eqordrcop:=dycopy:=nil; if (null truesub) or (eqlen neq length(truesub)) then rederr( "Unable to find all substitutions. Input equations as df(..,..)=..!"); >>; lhslist:=for each e1 in truesub collect cadr e1; %-- Bringing truesub into a specific form: lisp list of lisp lists: % ((to_be_replaced_jet_var_name, to_be_replaced_jet_var_deriv_1,..), % subst_expr_in_jet_space_coord, list_of_jet_vars_in_subst_expr) truesub:=for each e1 in truesub collect cons(combidif cadr e1, adddepli(caddr e1,revdylist))$ %--- Checking that no rhs of a substitution contains any lhs or %--- derivative of a lhs h:=t; %--- h=nil if lhs's are derivatives of each other no:=t; %--- no=nil if one lhs can be substituted in a rhs for each e1 in truesub do if h and no then << n1:=caar e1; n2:=cdar e1; dylen:=length n2; for each e2 in truesub do << %--- comparison of two lhs's if not(e1 eq e2) and (n1=caar e2) and comparedif1(n2,cdar e2) then h:=nil; %--- truesub is not ok %--- can the lhs of e1 be substituted on the rhs? dycopy:=caddr e2; for n:=1:dylen do if dycopy then dycopy:=cdr dycopy; for each e3 in dycopy do for each e4 in e3 do if comparedif2(n1,n2,e4) then no:=nil; >> >>; if null h then rederr( "One substitution can be made in the lhs of another substitution!"); if null no then rederr( "One substitution can be made in the rhs of another substitution!"); %--- Checking that a derivative of each dependent variable is %--- substituted once. This is a sufficient condition for having %--- a de-system that is a differential Groebner basis h:=nil; for each e1 in lhslist do h:=adjoin(car combidif e1,h); if length(h) neq length(lhslist) then rederr( "For at least one function there is more that one substituion!")$ %--- Determine of how much the order may increase by a substitution subordinc:=0; subordinclist:=for each h in truesub collect << n:=(length caddr h) - (length car h); if n>subordinc then subordinc:=n; n >>; if lietrace_ then <<terpri()$write"truesub=",truesub; terpri()$write"freelist=",freelist; %terpri()$write"occurlist=",occurlist >>; %--- To avoid non-uniqueness in the case of the investigation of %--- generalized symmetries let the characteristics eta_.. (xi_..=0) %--- not depend on substituted derivatives if generalp and (null ansatzp) then << deplist:=ylist . for each dycopy in cdr deplist collect << for each h in lhslist do %---- delete h and derivatives of h dycopy:=listdifdif1(h,dycopy); dycopy >>; for e1:=1:(length etalist) do << h:=nth(etalist,e1); depnd(car h,xlist . deplist); h:=pnth(h,3); rplaca(h,deplist) >> >>; % reduced set of solution techniques for preliminary conditions proc_list_:=delete('multintfac,proc_list_)$ if !*time then <<terpri()$ write "time for initializations: ", time() - cpu, " ms GC time : ", gctime() - gc," ms"$ cpu:=time()$ gc:=gctime()$ >>; %------ Determining first short determining equations and solving them symcon:=nil; n1:=0; if prelim_ then for each eqn in eqlist do <<n1:=n1+1; if !*time then <<terpri()$ terpri()$terpri()$ write"=============== Preconditions for the ",n1,". equation" ; >>; revdycopy:=revdylist; for e1:=(nth(eqordr,n1) + 1):ordr do revdycopy:=cdr revdycopy; n2:=cadr adddepli(eqn,revdycopy); % jet-variables in eqn vl:=n2; occli:=lastcar n2; freelist:=setdiff(car revdycopy,occli); if pointp and (subordinc=0) then eqn:=drop(eqn,occli) % dropp all terms without a highest deriv. else occli:=joinsublists n2$ % freelist must not contain substituted variables freelist:=setdiff(freelist,lhslist); % It must be possible to separate wrt freelist variables for each n4 in freelist do if not freeof(depl!*,n4) then freelist:=delete(n4,freelist); If freelist then << n:=nth(eqordr,n1); % order of this equation h:=simp!* 0; for each e1 in xilist do if (cadddr e1) and ((length cadddr e1) > n) then % xi (=car e1) is of order n h:=addsq(h, if car e1 = 0 then simp!* 0 else <<n3:=mergedepli(n3,cadddr e1); multsq(simp!* car e1, simpdf {eqn,cadr e1}) >> ); for each e2 in occli do h:=addsq(h, multsq(<<n4:=prolong(e2,xilist,etalist,nth(eqordr,n1), truesub,subordinc,xlist); vl:=mergedepli(vl,cdr n4); car n4 >>, simpdf {eqn,e2} ) ); for each e2 in freelist do << e1:=algebraic num lisp coeffn(prepsq h,e2,1); if not zerop e1 then symcon:=e1 . symcon>>; vl:=joinsublists(xlist . vl)$ for n2:=1:eqlen do <<n4:=nth(lhslist,n2); if not my_freeof(eqn,n4) then symcon:=subst(cadr nth(truesub,n2), n4, symcon); vl:=delete(n4,vl) >>; if symcon and (individual_ or (n1=eqlen)) then << flist:=callcrack(!*time,cpu,gc,lietrace_,symcon, flist,vl,xilist,etalist)$ symcon :=car flist; xilist:=cadr flist; etalist:=caddr flist; flist :=cadddr flist; cpu:=time()$ gc:=gctime()$ >> >> >>; %------------ Determining the full symmetry conditions n1:=0; vl:=nil; for each eqn in eqlist do <<n1:=n1+1; if !*time then <<terpri()$ terpri()$terpri()$ write"=============== Full conditions for the ",n1,". equation" ; >>; n2:=cadr adddepli(eqn,revdylist); n3:=n2; % n3 are the variables in the new condition symcon:=(reval algebraic num lisp prepsq addsq( <<h:=simp!* 0; for each e1 in xilist do h:=addsq(h, if car e1 = 0 then simp!* 0 else <<n3:=mergedepli(n3,cadddr e1); multsq(simp!* car e1, simpdf {eqn,cadr e1}) >> ); h >>, <<h:=simp!* 0; for each e1 in n2 do for each e2 in e1 do h:=addsq(h, multsq(<<n4:=prolong(e2,xilist,etalist,0,truesub, 0,xlist ); n3:=mergedepli(n3,cdr n4); car n4 >>, simpdf {eqn,e2} ) ); h >> )) . symcon; n3:=joinsublists(xlist . n3)$ for n2:=1:eqlen do <<n4:=nth(lhslist,n2); if not my_freeof(eqn,n4) then symcon:=subst(cadr nth(truesub,n2), n4, symcon); n3:=delete(n4,n3) >>; vl:=union(vl,n3); if individual_ or (n1=eqlen) then << flist:=callcrack(!*time,cpu,gc,lietrace_,symcon, flist,vl,xilist,etalist)$ symcon :=car flist; xilist:=cadr flist; etalist:=caddr flist; flist :=cadddr flist; cpu:=time()$ gc:=gctime()$ >> >>; eqn:=sb:=dylist:=e1:=e2:=n:=h:=deplist:=symord:=nil;%occurlist:=nil; lisp(adjust_fnc:=oldadj); %------- Calculation finished, simplification of the result h:='LIST . append(for each el in xilist collect car el, for each el in etalist collect car el ); %------- droping redundant constants or functions if null symcon then sb:=reval dropredundant(h,'LIST . flist,vl,list('LIST)); if sb then << flist:=cdr cadddr sb; h:=caddr sb; sb:=cadr sb; e1:=nil >>; %------- absorbing numerical constants into free constants h:=reval absorbconst(h,'LIST . flist); if h then if sb then sb:=append(sb,cdr h) else sb:='LIST . cdr h; %------- doing the substitutions if sb then << if print_ then <<terpri()$ write"Free constants and/or functions have been rescaled. ">>$ for each e1 in cdr sb do << xilist :=subst(caddr e1, reval cadr e1, xilist); etalist:=subst(caddr e1, reval cadr e1, etalist); symcon :=subst(caddr e1, reval cadr e1, symcon); % h:=intern car explode cadr e1; %write"h=",h;terpri()$ % if (h='x) or (h='X) then % xilist :=subst(caddr e1, cadr e1, xilist) else % if (h='e) or (h='E) or (h="e") or (h="E") then % etalist:=subst(caddr e1, cadr e1, etalist) else % rederr("One ansatz does not specify XI_ nor ETA_.") >>; >>; %------- output if length flist>1 then n:=t else n:=nil; terpri()$terpri()$ if null flist then write"No such symmetry does exist. " else write"The symmetr", if n then "ies are:" else "y is:"; terpri()$ xilist:=for each el in xilist collect <<e1:=mkid( xi,second el); for each e2 in fctargs(e1) do nodepend(e1,e2); e1:=list('EQUAL,e1,reval car el); e1>>; etalist:=for each el in etalist collect <<e1:=mkid(eta,second el); for each e2 in fctargs(e1) do nodepend(e1,e2); e1:=list('EQUAL,e1,reval car el); e1>>; %--- Backsubstitution of e.g. u`1`1 --> df(u,x,2) for each e1 in ylist do depnd(e1,list(xlist)); on evallhseqp; sb:=subdif1(cons('LIST,xlist),cons('LIST,ylist),jetord)$ algebraic( sb:=for each e1 in sb join for each e2 in e1 collect(rhs e2 = lhs e2)); off evallhseqp; xicop :=cdr algebraic(sub(sb,cons('LIST, xilist))); etacop:=cdr algebraic(sub(sb,cons('LIST,etalist))); sb:=nil$ flcop:=flist; n1:=0; for each e1 in flcop do <<n2:=cdr fargs e1; if null n2 then n2:=not freeof(symcon,e1); if null n2 then n2:=freeof(xicop,e1) and freeof(etacop,e1)$ if null n2 then << n1:=n1+1; xi:=xicop;eta:=etacop; for each e2 in flcop do if e2 neq e1 then <<xi:=subst(0,e2,xi);eta:=subst(0,e2,eta)>> else <<xi:=subst(1,e2,xi);eta:=subst(1,e2,eta)>>; terpri()$write"-------- ",n1,". Symmetry:";terpri()$ for each e2 in xi do algebraic write reval e2; for each e2 in eta do algebraic write reval e2; xicop :=subst(0,e1,xicop ); etacop:=subst(0,e1,etacop); flcop:=delete(e1,flcop); depl!*:=delete(assoc(e1,depl!*),depl!*)$ >>; >>; if flcop neq flist then << terpri()$write"-------- ";terpri()$ >>; if flcop then << if length flist>1 then n:=t else n:=nil; terpri()$ if flcop=flist then write"S" else write"Further s"; write"ymmetr",if n then "ies:" else "y:"$ terpri()$ for each e1 in xicop do algebraic write reval e1; for each e1 in etacop do algebraic write reval e1; >>; terpri()$ if flcop then <<write"with ";fctprint cdr reval ('LIST . flcop)>>; if null symcon then if flcop then <<write" which ",if n then "are" else "is"," free. "; terpri()>> else else <<h:=print_;print_:=500$ if print_ then <<terpri()$ write"which still ",if n then "have" else "has"," to satisfy: "; terpri()$ deprint symcon; >> else <<terpri()$ write"which ",if n then "have" else "has", " to satisfy conditions. To see them set "; terpri()$ write "lisp(print_= max. number of terms of an equation to print);"; terpri()$ >>; print_:=h >>; return list('LIST,'LIST . symcon,'LIST . append(xilist,etalist), 'LIST . flist); end$ % of liepde endmodule$ end$