Artifact c656aec7fefc579c690e24e17b9ba9aae239e80f6f0995a0267c6f47fe77d514:
- Executable file
r36/src/spde.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: 35320) [annotate] [blame] [check-ins using] [more...]
module spde; % Determine Lie symmetries of partial differential eqns. % Author: Fritz Schwarz. %*******************************************************************$ % $ % This is the REDUCE package SPDE for determining $ % Lie symmetries of partial differential equations $ % Version of November 1986 $ % $ % $ % Fritz Schwarz $ % GMD Institut F1 $ % Postfach 1240 $ % 5205 St. Augustin $ % West Germany $ % $ % Tel. 02241-142782 $ % EARN Id. DBNGMD21.GF1002 $ %*******************************************************************$ create!-package('(spde),'(contrib spde)); algebraic operator x,u,xi,eta,c,xi!*,eta!*$ algebraic operator deq,dx,du,gl,gen,sder,rule$ fluid '(depl!*); global'(pclass mm nn num!-cgen num!-dgen)$ share pclass,mm,nn$ lisp(pclass:=mm:=nn:=num!-cgen:=num!-dgen:=0)$ lisp(operator simpsys,result,prsys,prsys!*)$ fluid '(!*list kord!*)$ fluid'(uhf dfsub csub czero rdep !*rational)$ fluid'(list!-m list!-deq list!-pq)$ %symbolic procedure prload$ % begin % if not getd 'solve1 then load solve1,solvetab,quartic; % if not getd 'depend1 then load depend; % if not getd 'ratfunpri then load ratprin; % end$ symbolic procedure prload; nil; %*******************************************************************$ % Auxiliary RLISP procedures $ %*******************************************************************$ remflag('(ordp),'lose); % We must use this definition. symbolic procedure ordp(u,v)$ %Modified ordering function which orders kernels with CAR parts; %DF, ETA, XI and C ahead of anything else; if null u then null v else if null v then t else if eq(u,'df) or eq(u,'eta) and not eq(v,'df) or eq(u,'xi) and not(eq(v,'df) or eq(v,'eta)) or eq(u,'c) and not(eq(v,'df) or eq(v,'eta) or eq(v,'xi)) then t else if eq(u,'eta) and eq(v,'df) or eq(u,'xi) and (eq(v,'df) or eq(v,'eta)) or eq(u,'c) and (eq(v,'df) or eq(v,'eta) or eq(v,'xi)) or eq(v,'df) or eq(v,'eta) or eq(v,'xi) or eq(v,'c) then nil else if atom u then if atom v then if numberp u then numberp v and not(u<v) else if numberp v then t else orderp(u,v) else nil else if atom v then t else if car u=car v then ordp(cdr u,cdr v) else ordp(car u,car v)$ symbolic procedure makeset u$ if not u then nil else if member(car u,cdr u) then makeset cdr u else car u . makeset cdr u$ symbolic procedure lastmem u$ if cdr u then lastmem cdr u else car u$ symbolic procedure xmember(u,v)$ reverse member(u,reverse v)$ symbolic procedure sacar(a,u)$ if atom u then nil else if eq(a,car u) and cdr u then list u else append(sacar(a,car u),sacar(a,cdr u))$ symbolic procedure scar(a,u)$ if atom u then nil else if a=car u then u else scar(a,car u) or scar(a,cdr u)$ symbolic procedure inter(u,v); if not u then nil else if member(car u,v) then (car u) . inter(cdr u,v) else inter(cdr u,v)$ symbolic procedure compl(u,v)$ if not u then nil else if member(car u,v) then compl(cdr u,v) else car u . compl(cdr u,v)$ symbolic procedure vlist u$ %U is list of items, returns U with all integers omitted; if not u then nil else if numberp car u then vlist cdr u else (car u) . vlist cdr u$ symbolic procedure delnil u$ %U is list, returns U with all occurences of nil deleted; if not u then nil else if car u then (car u) . delnil cdr u else delnil cdr u$ symbolic procedure prlist u$ %U is list of items, returns list of all pairs in U; if not u then nil else if pairp car u then (car u) . prlist cdr u else prlist cdr u$ symbolic procedure appends(u,v,w)$ append(u,append(v,w))$ symbolic procedure propa(fn,u)$ %FN is predicate of a single argument, U a list; %Returns T if predicate is true for all elements of U; begin scalar ind; ind:=t; while ind and u do <<ind:=apply1(fn,car u); u:=cdr u>>; return ind; end$ symbolic procedure sortx(fn,u)$ begin scalar v,w; while u do<<v:=maxmem(fn,u); u:=delete(v,u); w:=v . w>>; return w; end$ symbolic procedure maxmem(fn,u)$ %FN is function of a single argument, U a list; %Returns element of U for which FN is maximal; begin scalar v; v:=car u; foreach x in cdr u do if greaterp(apply1(fn,x),apply1(fn,v)) then v:=x; return v; end$ symbolic procedure maxl u$ %U is list of integers, returns largest element of U; if not u then -10000 else max(car u,maxl cdr u)$ symbolic procedure suml u$ %U is list of integers, returns sum of all elements; if not u then 0 else plus2(car u,suml cdr u)$ symbolic procedure spde!-subsetp(u,v)$ %U and V are list representing sets; %Returns T if set U is subset of V; if not u then t else member(car u,v) and spde!-subsetp(cdr u,v)$ symbolic procedure product!-set2(u,v)$ %U and V are lists representing sets, returns list representing; %product set of sets represented by U and V; begin scalar w; foreach x in u do foreach y in v do w:=list(x,y) . w; return w; end$ symbolic procedure leqgrt(l,i,j)$ i leq j and eqn(l,i) or i geq add1 j$ symbolic procedure fidep u$ assoc(u,depl!*) and cdr assoc(u,depl!*)$ symbolic procedure mkdep u$ foreach x in cdr u do depend1(car u,x,t)$ symbolic procedure rmdep u$ <<rmsubs(); foreach x in cdr u do depend1(car u,x,nil)>>$ symbolic procedure blanks l; begin scalar u; u := '(!"); for k:=1:l do u:='! . u; return compress('!" . u) end$ symbolic procedure terpri2$ <<terpri(); terpri()>>$ %*******************************************************************$ % Auxiliary procedures for manipulating standard forms $ %*******************************************************************$ symbolic procedure lcf u$ not domainp u and lc u$ symbolic procedure minus!-f u$ %U is s.f., returns T if lnc U is negative; minusf numr simp reval u$ lisp operator minus!-f$ symbolic procedure lengthn u$ if not u then 0 else if numberp car u then plus(sub1 car u,lengthn cdr u) else plus(1,lengthn cdr u)$ symbolic procedure lengthf u$ %U is prefix s.f., returns printlength for U; if not u then 0 else if atom u then flatsizec u else if eqcar(u,'plus) then plus(times(3,sub1 length cdr u),lengthf cdr u) else if eqcar(u,'times) or eqcar(u,'minus) then plus(sub1 length cdr u,lengthf cdr u) else if eqcar(u,'quotient) then if !*rational then add1 add1 max(flatsizec cadr u,flatsizec caddr u) else add1 plus(flatsizec cadr u,flatsizec caddr u) else if eqcar(u,'expt) then add1 flatsizec cadr u else if eqcar(u,'dx) or eqcar(u,'du) then plus(flatsizec cadr u,4) else if eqcar(u,'xi) or eqcar(u,'eta) or eqcar(u,'c) or eqcar(u,'x) or eqcar(u,'u) then times(2,length u) else if eqcar(u,'df) then plus(4,lengthf cadr u,lengthf cddr u) else plus(lengthf car u,lengthf cdr u)$ lisp operator lengthf$ symbolic procedure diford u$ lengthn cddr u$ symbolic procedure adiff(u,v)$ %U is kernel with CAR part DF, V is kernel; %Returns U integrated with respect to V; if not member(v,u) then u else if length u=3 and member(v,u) then cadr u else if not cdr member(v,u) or not numberp cadr member(v,u) then delete(v,u) else if cadr member(v,u)=2 then append(xmember(v,u),cddr member(v,u)) else append(xmember(v,u),(sub1 cadr member(v,u)) . cddr member(v,u))$ symbolic procedure sub!-int!-df u$ %U is kernel with CAR part INT, returns integrated kernel if CADR; %part of U is DF and integration variable occurs as argument of DF; if eqcar(cadr u,'df) and member(lastmem u,cadr u) then adiff(cadr u,lastmem u) else u$ symbolic procedure subintf u$ %U is s.f., performs all integrations which may be done; %by cancellation of corresponding differentiation; begin foreach x in makeset sacar('int,u) do u:=subst(sub!-int!-df x,x,u); return numr simp prepf u; end$ symbolic procedure monop u$ %Returns T if u is monomial; domainp u or not red u and monop lc u$ symbolic procedure solvef(u,v)$ car solve0(prepf u,v)$ symbolic procedure comfacn u$ lnc ckrn u$ symbolic procedure remfacn u$ quotf(u,lnc ckrn u)$ %*******************************************************************$ % Procedures for manipulating l.d.f.'s, U is always l.d.f. $ % in this section $ %*******************************************************************$ symbolic procedure ldf!-mvar u$ %Returns function argument of mvar U; (if eqcar(x,'df) then cadr x else x) where x=mvar u; symbolic procedure ldf!-fvar u$ %Returns all function arguments occuring in U; makeset foreach x in u collect ldt!-tvar x$ symbolic procedure ldf!-fvar!-part(u,v)$ %V is function xi(i), eta(alpha) or c(k), returns l.d.f. of those; %terms in U with ldt-tvar x equal to V, overall factors not removed; begin scalar w; foreach x in u do if eq(ldt!-tvar x,v) then w:=x . w; return reverse w; end$ symbolic procedure ldf!-dep!-var u$ %Returns all variables x(i) or u(alpha) which occur as; %arguments of XI, ETA or C; begin scalar v; foreach x in u do if assoc(ldt!-tvar x,depl!*) then v:=append(cdr assoc(ldt!-tvar x,depl!*),v); return makeset v; end$ symbolic procedure ldf!-pow!-var u$ %Returns all variables x(i) or u(alpha) which occur as powers; begin scalar v,z; foreach x in u do v:=append(v,kernels tc x); foreach y in prlist makeset v do if eqcar(y,'x) or eqcar(y,'u) then z:=y . z; return makeset z; end$ symbolic procedure ldf!-deg(u,v)$ %V is kernel x(i) or u(alpha), returns degree of U in V; maxl foreach x in u collect degreef(tc x,v)$ symbolic procedure ldf!-spf!-var u$ %Returns all variables x(i) or u(alpha) which occur as; %arguments of any other kernel than xi, eta or c; begin scalar v,z; foreach x in u do v:=append(v,kernels tc x); foreach y in prlist makeset v do if not eqcar(y,'x) and not eqcar(y,'u) then z:=appends(sacar('x,cdr y),sacar('u,cdr y),z); return makeset z; end$ symbolic procedure ldf!-all!-var u$ %Returns all variables x(i) or u(alpha) which occur in U; makeset appends(ldf!-dep!-var u,ldf!-pow!-var u,ldf!-spf!-var u)$ symbolic procedure ldf!-sep!-var u$ %Returns all variables w.r.t. which U may be separated; compl(compl(ldf!-pow!-var u,ldf!-dep!-var u),ldf!-spf!-var u)$ symbolic procedure ldf!-int!-var u$ %Returns all variables w.r.t. which U may be integrated; if eqcar(mvar u,'df) then begin scalar v; v:=ldf!-all!-var u; while v and u do <<v:=compl(v,compl(ldt!-dep car u,ldt!-dfvar car u)); u:=cdr u>>; return v; end$ symbolic procedure ldf!-int u$ %U is l.d.f, returns U with all possible integrations performed; %or unchanged if integration is not possible; begin scalar v,w,z,test; integer nfun; a: test:=nil; w:=ldf!-int!-var u; nfun:=find!-nfun(); foreach x in w do if not smember('int,z:=caadr algebraic int(lisp prepf u,x)) or not smember('int,z:=subintf z) then <<v:=!*a2k list('c,nfun:=nfun+1); test:=t; mkdep(v . delete(x,ldf!-all!-var u)); u:=addf(z,!*k2f v)>>; if test then go to a; return u; end$ symbolic procedure ldf!-df!-diff u$ %Returns list of all df-kernels which may be obtained; %from U by differentiation or nil; begin scalar dfvar,dfsub,v,w,z0,z; integer n0,nmax; v:=compl(ldf!-dep!-var u,ldf!-spf!-var u); if not v then return; w:=foreach x in v collect list(x,add1 ldf!-deg(u,x)); nmax:=maxl foreach x in w collect cadr x; while (n0:=n0+1) leq nmax and not(z0:=nil) do <<foreach x in w do if cadr x geq n0 then z0:=(car x) . z0; z:=z0 . z>>; z:=reverse z; dfvar:=foreach x in car z collect list x; foreach x in cdr z do dfvar:= append(dfvar,foreach y in dfvar collect car product!-set2(x,y)); foreach x in dfvar do begin scalar p,q; p:=x; q:=u; while p and q and red q do <<q:=ldf!-simp numr difff(q,car p); p:=cdr p>>; if pairp q and not red q and eqcar(mvar q,'df) then dfsub:=(mvar q) . dfsub; end; return makeset dfsub; end$ symbolic procedure ldf!-sub!-var u$ %Returns function w.r.t. which U may be resolved; begin scalar v,w,z; w:=ldf!-all!-var u; foreach x in u do if not v and not eqcar(z:=tvar x,'df) and monop tc x and spde!-subsetp(w,ldt!-dep x) and not smember(z,delete(x,u)) then v:=z; return v; end$ symbolic procedure ldf!-simp u$ %Returns l.d.f. form of U; if not u then nil else if not red u then numr simp prepf !*k2f mvar u else begin scalar v; v:=numr simp prepf u; if not domainp v then v := quotf(v,cdr comfac v); return absf v end$ symbolic procedure ldf!-sep u$ %Returns list of l.d.f. into which U has been separated; begin scalar v; integer k; if not(v:=ldf!-sep!-var u) then return list u; foreach x in v do u:=subst(list('ux,1,k:=k+1),x,u); return foreach x in coeff!-all(u,'ux) collect ldf!-simp numr simp prepf x; end$ symbolic procedure ldf!-subf0 u$ %Returns U with CZERO substituted; ldf!-simp delnil foreach x in u collect ldt!-subt0 x$ %*******************************************************************$ % Procedures for manipulating l.d.t.'s, U is always l.d.t. $ % in this section $ %*******************************************************************$ symbolic procedure ldt!-tvar u$ %U is l.d.t., returns function argument of tvar U; (if eqcar(x,'df) then cadr x else x) where x=tvar u$ symbolic procedure ldt!-dfvar u$ %U is l.d.t., returns variables w.r.t. which tvar u is derived or nil; (if eqcar(x,'df) then vlist cddr x else nil) where x=tvar u$ symbolic procedure ldt!-dep u$ %U is l.d.t., returns list of variables x or y which occur as; %arguments LDT-tvar u; (if x then cdr x else nil) where x=assoc(ldt!-tvar u,depl!*)$ symbolic procedure ldt!-subt0 u$ %U is l.d.t., returns U if LDT-tvar u is not on czero; if not member(ldt!-tvar u,czero) then u else nil$ %*******************************************************************$ % Procedures for constructing the determining system $ %*******************************************************************$ symbolic procedure cresys u$ begin scalar r,v,w,lgl,lsub,depl!*!*,list!-sder; remprop('df,'kvalue); remprop('df,'klist); remprop('c,'kvalue); remprop('c,'klist); prload(); rmsubs(); depl!*:=nil; if car u then list!-deq:=foreach x in u collect assoc(x,get(car x,'kvalue)) else list!-deq:=get('deq,'kvalue); if eqn(length list!-deq,1) then begin scalar p; p:=maxmem(function length,makeset sacar('u,list!-deq)); p:=mk!*sq !*k2q p; list!-sder:=list list(list('sder, cadaar list!-deq),p); end else if car u then list!-sder:=foreach x in list!-deq collect assoc(list('sder,cadar x),get('sder,'kvalue)) else list!-sder:=get('sder,'kvalue); if not list!-deq then rerror(spde,1, "Differential equations not defined"); if not list!-sder then rerror(spde,2, "Substitutions for derivatives not defined"); mm:=find!-m list!-deq; nn:=find!-n list!-deq; list!-m:= makeset foreach x in sacar('u,list!-deq) collect cadr x; for k:=1:nn do<<w:=!*a2k list('xi,k) . w; v:=!*a2k list('x,k) . v>>; for k:=1:mm do if member(k,list!-m) then <<w:=!*a2k list('eta,k) . w; v:=!*a2k list('u,k) . v>>; for k:=1:nn do r:=(!*a2k list('dx,k)) . r; for k:=1:mm do r:=(!*a2k list('du,k)) . r; for k:=1:mm do depl!*!*:=(!*a2k list('eta,k) . v) . depl!*!*; for k:=1:nn do depl!*!*:=(!*a2k list('xi,k) . v) . depl!*!*; depl!*:=depl!*!*; kord!*:=reverse r; foreach x in list!-sder do lsub:=((mvar caadr cadr x) . prepsq caar solvef(caadr cadr assoc (list('deq,cadar x),list!-deq),mvar caadr cadr x)) . lsub; foreach x in list!-deq do begin scalar s,z,lx,lu; z:=caadr cadr x; lx:=makeset sacar('x,z); lu:=makeset sacar('u,z); foreach y in lx do s:=addf(s, multf(!*k2f !*a2k list('xi,cadr y),numr simp prepsq difff(z,y))); foreach y in lu do if length y=2 then s:=addf(s,multf (!*k2f !*a2k list('eta,cadr y),numr simp prepsq difff(z,y))) else s:=addf(s, multf(numr zeta!* cdr y,numr simp prepsq difff(z,y))); s:=numr subf(s,lsub); s:=numr subf(s,lsub); lgl:=append(coeff!-all(s,'u),lgl); end; uhf:=list(makeset lgl,foreach x in reverse w collect !*k2q x); end$ lisp rlistat'(cresys)$ symbolic procedure totder(u,i)$ begin scalar z,v,w; v:=car difff(u,!*a2k list('x,i)); z:=makeset sacar('u,u); for k:=1:mm do if member(k,list!-m) then z:=(!*a2k list('u,k)) . z; foreach x in makeset z do w:=addf(w, multf(!*k2f !*a2k append(x,list i),car difff(u,x))); return numr simp prepf addf(v,w); end$ symbolic procedure zeta!* u$ if not get('deq,'kvalue) and (eqn(mm,0) or eqn(nn,0)) then rerror(spde,3,"Number of variables not defined") else if length u geq 3 then begin scalar v,w; prload(); if eqn(nn,0) then nn:=find!-n list!-deq; v:=totder(numr zeta!* reverse cdr reverse u,car reverse u); for s:=1:nn do w:=addf(w, multf(!*k2f !*a2k('u . append(reverse cdr reverse u,list s)), totder(!*k2f !*a2k list('xi,s),car reverse u))); return simp prepsq(addf(v,negf w) ./ 1); end else begin scalar v,w; prload(); if eqn(nn,0) then <<nn :=find!-n list!-deq; mm:=find!-m list!-deq>> else begin scalar p,z; for k:=1:mm do z:=cons(k,z); for k:=1:nn do p:=(!*a2k list('x,k)) . p; for k:=1:mm do p:=(!*a2k list('u,k)) . p; for k:=1:nn do mkdep((!*a2k list('xi,k)) . p); for k:=1:mm do mkdep((!*a2k list('eta,k)) . p); list!-m:=z; end; v:=totder(!*k2f !*a2k list('eta,car u),cadr u); for s:=1:nn do w:=addf(w, multf(!*k2f !*a2k list('u,car u,s), totder(!*k2f !*a2k list('xi,s),car reverse u))); return simp prepsq(addf(v,negf w) ./ 1); end$ symbolic procedure simpu u$ !*p2q mksp(('u . (car u . reverse ordn cdr u)),1)$ put('u,'simpfn,'simpu)$ put('zeta,'simpfn,'zeta!*)$ symbolic procedure coeff!-all(u,v)$ begin scalar z; list!-pq:=nil; splitrec(u,v,1,nil); foreach x in list!-pq do z:=(ldf!-simp numr simp prepf cdr x) . z; return makeset z; end$ symbolic procedure splitrec(u,v,p,q)$ if domainp u then begin scalar y; p:=multf(u,p); if y:=assoc(q,list!-pq) then rplacd(y,addf(cdr y,p)) else list!-pq:=(q . p) . list!-pq; end else begin if eqcar(mvar u,v) and length mvar u greaterp 2 then splitrec(lc u,v,p,(lpow u) . q) else splitrec(lc u,v,!*t2f(lpow(u) .* p),q); if red u then splitrec(red u,v,p,q); end$ symbolic procedure find!-m u$ maxl makeset foreach x in sacar('u,u) collect cadr x$ symbolic procedure find!-n u$ begin scalar vx,vu,wx,wu; vx:=makeset sacar('x,u); vu:=makeset sacar('u,u); foreach x in vx do wx:=(cadr x) . wx; foreach x in vu do if length x geq 3 then wu:=append(cddr x,wu); return max(maxl wx,maxl wu); end$ %*******************************************************************$ % Procedures for solving the determining system $ %*******************************************************************$ symbolic procedure rule0$ %Searches for equations of the form C(I)=0 and stores them on CZERO; if uhf then foreach x in car uhf do if not red x and not eqcar(mvar x,'df) then czero:=(mvar x) . czero$ symbolic procedure rule1$ %Searches for equations of the form DF(function,variable)=0; %and stores it on the list RDEP; if uhf and car uhf then begin scalar dfsub; foreach x in car uhf do if not red x and eqcar(mvar x,'df) and eqn(diford mvar x,1) then rdep:=(mvar x) . rdep; if rdep then return t; end$ symbolic procedure rule1!-diff$ %Searches for equations of the form DF(function,variable)=0; %which may be obtained by a single differentiation and stores it on; %the list RDEP; if uhf and car uhf then begin scalar u,v,z; foreach x in car uhf do if(z:=ldf!-df!-diff x) then u:=append(z,u); foreach x in u do if eqn(diford x,1) then v:=x . v; rdep:=makeset v; if rdep then return t; end$ symbolic procedure rulec l$ %Searches for equations of length L which may be solved for a; %function and stores the corresponding rules on CSUB; if uhf and car uhf then begin scalar v; foreach u in car uhf do if leqgrt(length u,l,4) and (v:=ldf!-sub!-var u) and not smember(v,csub) and not inter(foreach x in csub collect car x,ldf!-fvar u) then csub:=(v . prepsq caar solvef(u,v)) . csub; if csub then return t; end$ symbolic procedure ruledf l$ %Searches for equations of the form DF(function,derivative list)=0; %the derivative beeing of order L and stores the resulting; %substitution polynomial on CSUB; if uhf and car uhf then begin scalar dfsub; foreach x in car uhf do if not red x and eqcar(mvar x,'df) and eqn(diford mvar x,l) and not smember(ldf!-mvar x,dfsub) then dfsub:=(mvar x) . dfsub; csub:=foreach x in dfsub collect(cadr x) . crepol x; if csub then return t; end$ symbolic procedure ruledf!-diff l$ %Searches for all equations of the form; %DF(function,derivative list)=0 which may be obtained by; %differentiation, picks out those of order L and stores; %the corresponding substitution polynomial on CSUB; if uhf and car uhf then begin scalar v,dfsub; foreach u in car uhf do v:=append(v,ldf!-df!-diff u); if not(v:=makeset v) then return; foreach x in v do if eqn(diford x,l) then dfsub:=x . dfsub; if not dfsub then return; csub:=((cadar dfsub) . crepol car dfsub) . csub; if csub then return t; end$ symbolic procedure rule!-int l$ %Searches for an equation of length L which may be solved for a; %function after beeing integrated and stores the corresponding; %rule on CSUB; if uhf and car uhf then begin scalar v,w; foreach u in car uhf do if not csub and leqgrt(length u,l,4) and (v:=ldf!-sub!-var(w:=ldf!-int u)) then csub:=list(v . prepsq caar solvef(w,v)); if csub then return t; end$ symbolic procedure simpsys0$ %Removes variable which are stored on list CZERO; begin scalar u,v; if pclass=2 then<<write"Entering SIMPSYS0"; terpri2()>>; u:=delnil foreach x in car uhf collect ldf!-subf0 x; v:=foreach x in cadr uhf collect ldf!-subf0 numr x ./ denr x; uhf:=list(makeset u,v); if pclass=1 then begin terpri2(); if eqn(length czero,1) then write"Substitution" else write"Substitutions"; terpri(); foreach x in czero do algebraic write (lisp aeval x),":=0"; terpri(); end; if pclass=2 then<<write"CZERO:="; prettyprint czero; terpri()>>; czero:=nil; if pclass=2 then<<write"Leaving SIMPSYS0"; terpri2()>>; end$ symbolic procedure simpsys!-rdep$ %Removes dependencies which are stored on list RDEP; begin scalar u,v; if pclass=2 then<<write"Entering SIMPSYS!-RDEP"; terpri2()>>; foreach x in rdep do rmdep cdr x; u:=makeset delnil foreach x in car uhf collect ldf!-simp x; v:=foreach x in cadr uhf collect simp prepsq x; uhf:=list(u,v); if pclass=1 then begin terpri(); write"Dependencies removed"; terpri2(); foreach x in rdep do <<maprin cadr x; prin2!*" independent of "; maprin caddr x; terpri!* t;>>; terpri(); end; if pclass=2 then<<write"RDEP:='"; prettyprint rdep; terpri()>>; if pclass=2 then<<write"Leaving SIMPSYS!-RDEP"; terpri2()>>; end$ symbolic procedure simpsys!-sep$ %Performs all possible separations; if uhf and car uhf then begin scalar u,v,test; if pclass=2 then<<write"Entering SIMPSYS!-SEP"; terpri2()>>; foreach x in car uhf do if eqn(length(v:=ldf!-sep x),1) then u:=x . u else begin u:=append(v,u); if pclass=1 or pclass=2 then begin scalar z; integer l; terpri(); l:=length car uhf-length member(x,car uhf)+1; write"Equation ",l," separated into the terms"; terpri(); if pclass=1 then for k:=1:length v do begin z:=prepf nth(v,k); !*list := lengthf z geq 50; algebraic write"Term ",k," ",z; end; if pclass=2 then foreach y in v do prettyprint y; end; test:=t; end; !*list := nil; if test then uhf:=list(reverse makeset u,cadr uhf); if pclass=2 then<<write"Leaving SIMPSYS!-SEP"; terpri2()>>; end$ symbolic procedure simpsys!-sub$ %Performs all substitutions which are stored on CSUB; if uhf and car uhf then begin scalar u,v; if pclass=2 then<<write"Entering SIMPSYS!-SUB"; terpri2()>>; if pclass=1 then prrule csub; if pclass=2 then<<write"CSUB:='"; prettyprint csub; terpri()>>; u:=makeset delnil foreach x in car uhf collect ldf!-simp numr subf(x,csub); v:=foreach x in cadr uhf collect subsq(x,csub); uhf:=list(u,v); csub:=nil; if pclass=2 then<<write"Leaving SIMPSYS!-SUB"; terpri2()>>; end$ symbolic procedure simpsys$ if not uhf then rerror(spde,4,"The determining system is not defined") else if not car uhf then rerror(spde,5,"The determining system completely solved") else begin scalar u,v; integer nfun; prload(); u:=makeset delnil foreach x in car uhf collect ldf!-simp x; v:=foreach x in cadr uhf collect simp prepsq x; uhf:=list(u,v); mark0: if pclass=1 then<<prsys!*"Entering main loop">> else if pclass=2 then prtlist"Entering main loop"; czero:=csub:=rdep:=nil; simpsys!-sep(); rule0(); if czero then<<simpsys0(); go to mark0>>; if rule1() or rule1!-diff() then<<simpsys!-rdep(); go to mark0>>; if ruledf 2 or rulec 2 or rule!-int 2 or ruledf!-diff 2 or ruledf 3 or rulec 3 or rule!-int 3 or ruledf!-diff 3 or ruledf 4 or rulec 4 or rule!-int 4 or ruledf!-diff 4 or ruledf 5 or rulec 5 or rule!-int 5 or ruledf!-diff 5 then <<simpsys!-sub(); go to mark0>>; if car uhf then <<write"Determining system is not completely solved"; terpri2(); prsys!*"The remaining equations are"; if not zerop(nfun:=find!-nfun()) then write"Number of functions is ",nfun>>; end$ symbolic procedure crepol u$ begin scalar l1,f; integer pow,nfun; nfun:=find!-nfun(); l1:=cdr assoc(car(u:=cdr u),depl!*); while (u:=cdr u) do begin scalar v; v:=car u; if length u=1 or not numberp cadr u then pow:=1 else <<pow:=cadr u; u:=delete(pow,u);>>; for k:=1:pow do begin scalar w; w:=!*a2k list('c,nfun:=nfun+1); mkdep(w . delete(v,l1)); if k=1 then f:=w . f; if k=2 then f:=list('times,w,v) . f; if k geq 3 then f:=list('times,w,list('expt,v,k-1)) . f; end; end; return append('(plus),f); end$ %*************************************************************$ % Procedures for analysing the result $ %*************************************************************$ symbolic procedure cpar u$ begin scalar v; v:=makeset appends(sacar('xi,u),sacar('eta,u),sacar('c,u)); foreach x in v do if not assoc(x,depl!*) then v:=delete(x,v); return v; end$ symbolic procedure makeset!-c!-x u$ if not u then nil else if member!-c!-x(car u,cdr u) then makeset!-c!-x cdr u else car u . makeset!-c!-x cdr u$ symbolic procedure member!-c!-x(u,v)$ if not v then nil else if equal!-c!-x(u,car v) then v else member!-c!-x(u,cdr v)$ symbolic procedure equal!-c!-x(u,v)$ begin scalar p,q; p:=scar('c,u) or scar('xi,u) or scar('eta,u); q:=scar('c,v) or scar('xi,v) or scar('eta,v); return equal(subst('cxx,p,u),subst('cxx,q,v)); end$ symbolic procedure numgen$ length get('gen,'kvalue)$ symbolic operator numgen$ symbolic procedure gengen$ begin scalar u,z,cgen,dgen; integer ngen; remprop('gen,'kvalue); remprop('gen,'klist); foreach x in cadr uhf do u:=append(ldf!-fvar numr x,u); foreach x in makeset u do begin scalar v,w; w:=nil ./ 1; if assoc(x,depl!*) then v:=foreach y in cadr uhf collect simp prepsq(ldf!-fvar!-part(numr y,x) ./denr y) else v:=foreach y in cadr uhf collect simp prepsq((lcf ldf!-fvar!-part(numr y,x)) ./denr y); for k:=1:nn do if numr nth(v,k) then w:=addsq(multsq(nth(v,k),!*k2q !*a2k list('dx,k)),w); for k:=1:mm do if numr nth(v,nn+k) then w:=addsq(multsq(nth(v,nn+k),!*k2q !*a2k list('du,k)),w); if assoc(x,depl!*) then cgen:=(absf remfacn numr simp prepf numr w) . cgen else dgen:=(absf remfacn numr simp prepf numr w) . dgen; end; dgen:=makeset dgen; cgen:=makeset!-c!-x cgen; num!-dgen:=length dgen; num!-cgen:=length cgen; for k:=1:nn do if member(z:=!*k2f !*a2k list('dx,k),dgen) then <<setk(list('gen,ngen:=add1 ngen),prepf z); dgen:=delete(z,dgen)>>; for k:=1:mm do if member(z:=!*k2f !*a2k list('du,k),dgen) then <<setk(list('gen,ngen:=add1 ngen),prepf z); dgen:=delete(z,dgen)>>; dgen:=sortx(function length,dgen); foreach x in dgen do setk(list('gen,ngen:=add1 ngen),prepf x); cgen:=sortx(function length,cgen); foreach x in cgen do setk(list('gen,ngen:=add1 ngen),prepf x); end$ symbolic operator gengen$ algebraic procedure comm(a,b)$ begin scalar z; if (lisp length list!-deq)=0 then <<write"Differential equations not defined"; return nil>>; z:= (for k:=1:nn sum df(a,dx k)*df(b,x k)-df(b,dx k)*df(a,x k)) +(for k:=1:mm sum df(a,du k)*df(b,u k)-df(b,du k)*df(a,u k))$ return z; end$ algebraic procedure result$ begin integer l; if (l:=lisp length list!-deq)=1 then write"The differential equation" else write"The differential equations"; for j:=1:l do begin scalar z; integer i,k; lisp(z:=car cadadr nth(list!-deq,j)); i:=lisp cadar nth(list!-deq,j); k:=lisp lengthf prepf z; symbolic(!*list := k>40); write"DEQ(",i,"):=",lisp prepf z; end; !*list := nil; if (lisp length car uhf) neq 0 then prsys!*"The determining system is not completely solved" else <<lisp gengen(); prgen(); comm!-tab()>>; end$ %*************************************************************$ % Procedures for displaying the output $ %*************************************************************$ symbolic procedure prsys!* u$ if uhf and car uhf then <<terpri(); write u; terpri(); prsys(); terpri()>>$ symbolic procedure prsys$ begin scalar v; terpri(); remprop('gl,'kvalue); remprop('gl,'klist); for k:=1:length car uhf do begin scalar z; integer l; z:=prepf nth(car uhf,k); l:=lengthf prepf nth(car uhf,k); !*list := l>50; algebraic write"GL(",k,"):=",z; setk(list('gl,k),z); end; terpri2(); write"The remaining dependencies"; terpri2(); v:=makeset appends(sacar('xi,car uhf),sacar('eta,car uhf),sacar('c,car uhf)); foreach x in v do write!-dep x; !*list := nil; end$ symbolic procedure prrule u$ begin terpri2(); if eqn(length u,1) then write"Substitution" else write"Substitutions"; terpri2(); foreach x in u do <<maprin car x; prin2!*" = "; maprin cdr x; terpri!* t;>>; terpri(); foreach x in u do foreach y in sacar('c,cdr x) do write!-dep y; end$ symbolic procedure prtlist u$ <<write u; terpri2(); write"DEPL!*:='"; prettyprint depl!*; write"UHF:='"; prettyprint uhf>>$ symbolic procedure write!-df!-sub$ if get('df,'kvalue) then begin scalar w; w:=get('df,'kvalue); remprop('df,'kvalue); terpri(); if length w=1 then write"Constraint" else write"Constraints"; terpri2(); foreach x in w do begin scalar u,v; u:=car x; v:=cadadr x; algebraic write lisp u,":=",lisp prepsq v; terpri2(); end; put('df,'kvalue,w); end$ algebraic procedure prgen$ begin scalar lcpar; for k:=1:nn do <<order dx k; factor dx k>>; for k:=1:mm do factor du k$ lisp(lcpar:=cpar get('gen,'kvalue)); write"The symmetry generators are"; for k:=1:numgen() do if (lisp lengthf reval list('gen,k)) leq 60 then <<symbolic(!*list := nil); write"GEN(",k,"):=",gen k>> else begin scalar z; integer r,s,nt; operator gen!*; nt:=lisp length(z:=numr simp reval list('gen,k)); r:=lisp maxl foreach x in z collect abs comfacn list x; if r=1 then r:=0 else r:=lisp flatsizec r; for l:=1:nt do gen!* l:=lisp prepf list nth(z,l); for l:=1:nt do begin symbolic(!*list := lengthf prepf tc nth(z,l) geq 56); s:=lisp abs comfacn list nth(z,l); if r=0 then s:=0 else if s=1 then s:=-1 else s:=lisp flatsizec s; if l=1 then write"GEN(",k,"):=",lisp blanks(r-s+1),gen!* 1 else if minus!-f gen!* l then write lisp blanks(r-s+6),gen!* l else write lisp blanks(r-s+6)," + ",gen!* l; end; clear gen!*; symbolic(!*list := nil); end; if (lisp length lcpar) neq 0 then <<write"The remaining dependencies"; lisp terpri()>>; for k:=1:(lisp length lcpar) do <<lisp write!-dep nth(lcpar,k);>>; if (lisp length lcpar) neq 0 then lisp terpri(); lisp write!-df!-sub(); end$ algebraic procedure comm!-tab$ if (lisp num!-dgen) geq 2 then begin integer nd; scalar v; nd:=lisp num!-dgen; write"The non-vanishing commutators of the finite subgroup"; for i:=1:nd-1 do for j:=(i+1):nd do if(v:=comm(gen i,gen j)) neq 0 then if (lisp lengthf reval v) leq 60 then <<symbolic(!*list := nil); write"COMM(",i,",",j,"):= ",v>> else begin integer r,s,nt; scalar z; operator gen!*; nt:=lisp length(z:=numr simp reval v); r:=lisp maxl foreach x in z collect abs comfacn list x; if r=1 then r:=0 else r:=lisp flatsizec r; for i:=1:nt do gen!* i:=lisp prepf list nth(z,i); for l:=1:nt do begin symbolic(!*list := lengthf reval list('gen!*,l) geq 63); s:=lisp abs comfacn list nth(z,l); if r=0 then s:=0 else if s=1 then s:=-1 else s:=lisp flatsizec s; if l=1 then write"COMM(",i,",",j,"):=",lisp blanks(r-s+1),gen!* 1 else if minus!-f gen!* l then write lisp blanks(r-s+9),gen!* l else write lisp blanks(r-s+9)," + ",gen!* l; end; clear gen!*; end; symbolic(!*list := nil); end$ symbolic procedure write!-dep u$ if assoc(reval u,depl!*) then begin scalar v; v:=cdr assoc(u,depl!*); write car u,"(",cadr u,") depends on "; write caar v,"(",cadar v,")"; foreach x in cdr v do write",",car x,"(",cadr x,")"; terpri2(); end$ symbolic operator write!-dep$ symbolic procedure find!-nfun$ if not get('c,'klist) then 0 else maxl makeset foreach x in get('c,'klist) collect cadar x$ endmodule; end;