Artifact 728daedd9b4201601153143791223133b71167e168128a59cea5f54f06779436:
- Executable file
r38/packages/crack/crident.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: 19729) [annotate] [blame] [check-ins using] [more...]
%******************************************************************** module identities$ %******************************************************************** % Routines for dealing with differential identities % Author: Thomas Wolf % May 1999 symbolic procedure drop_idty(id)$ % recycles a name of an identity <<setprop(id,nil)$ recycle_ids:=id . recycle_ids$ idnties_:=delete(id,idnties_)>>$ symbolic procedure new_id_name$ % provides a name for a new identity begin scalar id$ if null recycle_ids then << id:=mkid(idname_,nid_)$ nid_:=add1 nid_$ >> else << id:=car recycle_ids$ recycle_ids:=cdr recycle_ids >>$ setprop(id,nil)$ return id end$ symbolic procedure replace_idty$ begin scalar p,ps,ex$ ps:=promptstring!*$ promptstring!*:=""$ terpri()$ write "If you want to replace an identity then ", "type its name, e.g. id_2 <ENTER>."$ terpri()$ write "If you want to add an identity then type `new_idty' <ENTER>. "$ p:=termread()$ if (p='NEW_IDTY) or member(p,idnties_) then <<terpri()$write "Input of a value for "$ if p='NEW_IDTY then write "the new identity." else write p,"."$ terpri()$ write "You can use names of other identities, e.g. 3*id_2 - df(id_1,x); "$ terpri()$ write "Terminate the expression with ; or $ : "$ terpri()$ ex:=termxread()$ for each a in idnties_ do ex:=subst(get(a,'val),a,ex)$ if p neq 'NEW_IDTY then drop_idty(p)$ new_idty(reval ex,nil,nil)$ terpri()$write car idnties_$ if p='NEW_IDTY then write " is added" else write " replaces ",p$ >> else <<terpri()$ write "An identity ",p," does not exist! (Back to previous menu)">>$ promptstring!*:=ps$ end$ symbolic procedure trivial_idty(pdes,idval)$ if zerop idval or (pdes and null smemberl(pdes,search_li(idval,'DF))) % identity is purely alg. then t else nil$ symbolic procedure new_idty(idty,pdes,simp)$ % idty is the value of a new differential identity between equations begin scalar id,idcp$ if simp then idty:=simplifypde(idty,pdes,t,nil)$ if not trivial_idty(pdes,idty) then << idcp:=idnties_$ while idcp and (get(car idcp,'val) neq idty) do idcp:=cdr idcp; if null idcp then << id:=new_id_name(); put(id,'val,idty)$ flag1(id,'to_subst)$ flag1(id,'to_int)$ idnties_:=cons(id,idnties_) >> >> end$ symbolic procedure show_id$ begin scalar l,n$ terpri()$ l:=length idnties_$ write if l=0 then "No" else l, if l=1 then " identity." else " identities"$ if l=0 then terpri() else << n:=1; for each l in reverse idnties_ do << terpri()$ algebraic write n,") ",l," : 0 = ",lisp(get(l,'val)); n:=add1 n; if print_all then << terpri()$write " to_int : ",flagp(l,'to_int)$ terpri()$write " to_subst : ",flap(l,'to_subst)$ >> >> >> end$ symbolic procedure del_red_id(pdes)$ begin scalar oldli,pl,s,idty,news,succ,p,l$ % ,r,newr$ if idnties_ then << oldli:=idnties_$ while oldli do if not flagp(car oldli,'to_subst) then oldli:=cdr oldli else << idty:=get(car oldli,'val)$ pl:=smemberl(pdes,idty)$ for each p in pl do l:=union(get(p,'vars),l)$ if l then l:=length l else l:=0$ pl:=setdiff(pl,search_li(idty,'DF)); % now all pdes in pl are redundand, drop the longest if null pl then remflag1(car oldli,'to_subst) else << drop_idty(car oldli); % find the longest equation s of those with the most variables s:=nil; while pl do << if (null get(car pl,'starde) ) and (get(car pl,'nvars)=l ) and (null(% flagp(s,'to_int) or % flagp(s,'to_fullint) or % flagp(s,'to_sep) or % flagp(s,'to_gensep) or % flagp(s,'to_decoup) or flagp(s,'to_eval)) ) and ((null s ) or (get(car pl,'nvars)>get(s,'nvars) ) or ((get(car pl,'nvars)=get(s,'nvars)) and (get(car pl,'terms)>get(s,'terms)) ) ) then s:=car pl; pl:=cdr pl >>; if null s then remflag1(car oldli,'to_subst) else << if print_ then << write "Equation ",s," is dropped as it is a consequence of others: "$ algebraic write "0 = ",lisp(idty)$ >>$ % assuming s occurs linearly: pl:=coeffn(idty,s,1)$ news:=reval {'QUOTIENT,{'DIFFERENCE,{'TIMES,pl,s},idty},pl}; %for each r in idnties_ do %if not freeof(get(r,'val),s) then << % newr:=reval subst(news,s,get(r,'val)); % newr:=simplifypde(newr,pdes,t,nil)$ % put(r,'val,newr)$ % flag1(r,'to_subst)$ % flag1(r,'to_int)$ %>>$ succ:=t$ pdes:=drop_pde(s,pdes,news)$ oldli:=cdr oldli >> >> >> >>; if succ then return pdes end$ symbolic procedure del_redundant_de(argset)$ begin scalar pdes; if pdes:=del_red_id(car argset) then return {pdes,cadr argset}$ end$ symbolic procedure write_id_to_file(pdes)$ begin scalar s,p,h,pl,ps$ if idnties_ then << ps:=promptstring!*$ promptstring!*:=""$ write"Please give the name of the file in double quotes"$terpri()$ write"without `;' : "$ s:=termread()$ out s; off nat$ write"load crack$"$terpri()$ write"lisp(nequ_:=",nequ_,")$"$terpri()$ write"off batch_mode$"$terpri()$ write"list_of_variables:="$ algebraic write lisp cons('LIST,vl_)$ write"list_of_functions:="$ algebraic write lisp cons('LIST,pdes)$ for each h in pdes do if pl:=assoc(h,depl!*) then for each p in cdr pl do algebraic write "depend ",lisp h,",",lisp p$ write"list_of_equations:="$ algebraic write lisp( cons('LIST,for each h in idnties_ collect get(h,'val))); terpri()$ write"solution_:=crack(list_of_equations,{},"$ terpri()$ write" list_of_functions,"$ terpri()$ write" list_of_variables)$"$ terpri()$ terpri()$ write"end$"$terpri()$ shut s; on nat; promptstring!*:=ps$ >> end$ symbolic procedure remove_idl$ <<for each h in idnties_ do setprop(h,nil); idnties_:=nil>>$ symbolic procedure start_history(pdes)$ begin scalar l,ps$ ps:=promptstring!*$ promptstring!*:=""$ write"For recording the history of equations all currently"$ terpri()$ write"recorded histories would be deleted as well as all"$ terpri()$ write"present decoupling information, i.e. `dec_with'"$ terpri()$ write"would be set to nil. Please confirm (y/n). "$ l:=termread()$ if (l='y) or (l='Y) then << record_hist:=t; for each l in pdes do put(l,'histry_,l)$ for each l in pdes do put(l,'dec_with,nil)$ >>; promptstring!*:=ps$ end$ symbolic procedure stop_history(pdes)$ <<record_hist:=nil; for each l in pdes do put(l,'histry_,l)>>$ % write"Do you want to delete all dec_with information? (y/n) "$ % l:=termread()$ % if (l='y) or (l='Y) then % for each l in pdes do put(l,'dec_with,nil)$ symbolic procedure idty_integration(argset)$ begin scalar l,pdes,idcp; pdes:=car argset; idcp:=idnties_; while idcp do if not flagp(car idcp,'to_int) then idcp:=cdr idcp else if l:=integrate_idty(car idcp,pdes,%cadr argset, ftem_,vl_) then << pdes:=l;idcp:=nil>> else << remflag1(car idcp,'to_int); idcp:=cdr idcp; >>; if l then return {pdes,cadr argset} end$ symbolic procedure integrate_idty(org_idty,allpdes,%forg, fl,vl)$ % idty is a differential identity between equations % allpdes, fl, vl are lisp lists of equation names, functions and variables % ways to optimize: use conlaw instead of the initial intcurrent2 % use more general methods to take advantage of % non-conservation laws if idnties_ then begin scalar cl,ncl,vlcp,xlist,eql,a,f,newpdes,ftem_bak, nx,dl,l,k,ps,idty,pdes,extrapdes,newidtylist$ %nclu if null org_idty then if null cdr idnties_ then org_idty:=car idnties_ else << show_id()$ ps:=promptstring!*$ promptstring!*:=""$ write"Which of the identities shall be integrated? (no) "$ k:=length(idnties_); repeat l:=termread() until (fixp l) and (0<l) and (l<=k); org_idty:=nth(idnties_,k+1-l)$ promptstring!*:=ps >>$ idty:=reval num reval get(org_idty,'val)$ if trivial_idty(allpdes,idty) then return nil$ pdes:=smemberl(allpdes,idty)$ a:=all_deriv_search(idty,pdes)$ xlist:=smemberl(vl,a)$ cl:=intcurrent3(idty,cons('LIST,pdes),cons('LIST,xlist))$ % intcurrent3 is only successful if only 2 derivatives found if (not zerop caddr cl) and inter_divint then cl:=intcurrent2(idty,cons('LIST,pdes),cons('LIST,xlist))$ if zerop caddr cl then << cl:=cdadr cl; vlcp:=xlist; xlist:=nil; while vlcp do << if not zerop car cl then << ncl:=cons(car cl,ncl); xlist:=cons(car vlcp,xlist) >>; cl:=cdr cl; vlcp:=cdr vlcp >>; % ncl:=reverse ncl; % xlist:=reverse xlist; cl:=ncl; % % Now try to get a divergence in less differentiation variables. % % Each component of the divergence is tried to be written as % % a divergence in the other (right) variables % while ncl do << % a:=intcurrent2(car ncl,cons('LIST,pdes),cons('LIST,cdr xlist))$ % if not zerop caddr a then << % cl:=cons(car ncl,cl); ncl:=cdr ncl; % vlcp:=cons(car xlist,vlcp); xlist:=cdr xlist % >> else << % % It was possible to integrate car ncl to div(cdadr a,cdr xlist). % % distribute {'DF,car a,car xlist} to the divergence of cdr ncl % ncl:=cdr ncl; % a:=cdadr a; % nclu:=nil; % while ncl do << % nclu:=cons(reval {'PLUS,car ncl,{'DF,car a,car xlist}}, nclu); % ncl:=cdr ncl; % a:=cdr a % >>; % ncl:=reverse nclu; % xlist:=cdr xlist % >> % >>$ % ncl:=cl; % xlist:=vlcp; nx:=length xlist; while pdes do << ncl:=subst(get(car pdes,'val),car pdes,ncl); pdes:=cdr pdes >>$ ftem_bak:=ftem_; eql:=int_curl(reval cons('LIST,ncl), cons('LIST,fl), cons('LIST,xlist),cons('LIST,varslist(ncl,ftem_,vl)) )$ % eql has the form {'LIST,reval cons('LIST,resu),cons('LIST,neweq)} if (null eql) or (null cdadr eql) or (zerop cadadr eql) then return nil; eql:=cdr eql; if print_ then << ncl:=for i:=1:nx collect {'DF,nth(cl,i),nth(xlist,i)}; ncl:=if cdr ncl then cons('PLUS,ncl) else car ncl; terpri()$ write"The identity "$ % mathprint idty$ mathprint reval ncl; write"can be integrated to "$terpri()$ deprint(cdar eql)$ >>$ if nx < 3 then a:='y else if (null inter_divint) or !*batch_mode then << a:='n; if print_ then << write"The integrated divergence is not used because it ", "has more than 2 terms and"$ terpri()$ if !*batch_mode then write"`inter_divint' is nil." else write"batch_mode is on."$ >>$ terpri() >> else << ps:=promptstring!*$ promptstring!*:=""$ write"Shall this integration be used? (y/n) "$ repeat a:=termread() until (a='y) or (a='n); promptstring!*:=ps >>; if a='n then << a:=setdiff(ftem_,ftem_bak); for each f in a do drop_fct(f)$ ftem_:=ftem_bak >> else << % the extra conditions from the generalized integration: extrapdes:=cdadr eql$ eql:=cdar eql; % eql are now the integrated curl conditions drop_idty(org_idty)$ while eql do << if not zerop car eql then << a:=mkeq(car eql,ftem_,vl,allflags_,nil,list(0),nil,allpdes); newpdes:=cons(a,newpdes); >>; eql:=cdr eql; >>; newpdes:=reverse newpdes; % formulate the new identities for i:=1:nx do << idty:=nth(cl,i); if nx=1 then a:=car newpdes else << % at first sum over df(q^{ji},j), j<i l:=i-1; dl:=nx-2; a:=for j:=1:(i-1) collect << k:=l; l:=l+dl; dl:=sub1 dl; {'DF,nth(newpdes,k),nth(xlist,j)} >>; a:=if null a then 0 else if cdr a then cons('PLUS,a) else car a; idty:={'PLUS,idty,a}; % then sum over -df(q^{ij},j), j>i if i=1 then l:=1 else l:=k+nx-i+1; a:=for j:=(i+1):nx collect << k:=l; l:=l+1; {'DF,nth(newpdes,k),nth(xlist,j)} >>; a:=if null a then 0 else if cdr a then cons('PLUS,a) else car a; >>$ newidtylist:=cons({'DIFFERENCE,idty,a},newidtylist); >>; eql:=nil; for each a in extrapdes do << a:=mkeq(a,ftem_,vl,allflags_,t,list(0),nil,allpdes); allpdes:=eqinsert(a,allpdes); to_do_list:=cons(list('subst_level_35,%allpdes,forg,vl_, list a), to_do_list); eql:=cons(a,eql) >>; if print_ then << write"Integration gives: "$ listprint(newpdes)$terpri()$ if eql then << write"with extra conditions: "$ listprint(eql) >>$ >>; for each a in newpdes do allpdes:=eqinsert(a,allpdes)$ % now that allpdes is updated: for each a in newidtylist do new_idty(a,allpdes,t)$ return allpdes >> >> end$ symbolic procedure sortpermuli(a)$ % a is a list of numbers to be sorted and the exchanges of neighbours % are to be counted begin scalar flp,conti,newa; repeat << newa:=nil; conti:=nil; while cdr a do if car a < cadr a then <<newa:=cons(car a,newa); a:=cdr a>> else << conti:=t; flp:=not flp; newa:=cons(cadr a,newa); a:=cons(car a,cddr a); >>$ newa:=cons(car a,newa); a:=reverse newa >> until null conti; return flp . a end$ symbolic procedure curlconst(xlist,vl)$ % generates a list q^ij=r^ijk,_k with r^ijk totally antisymmetric % in the order q^(n-1)n,... % xlist is the list of xi,xj,xk % vl is the list of all variables new functions should depend on begin scalar n,qli,i,j,k,qij,a,flp,f,resu,qlicp$ n:=length xlist$ for i:=1:(n-1) do for j:=(i+1):n do << % generation of r^ijk,k qij:=nil; for k:=1:n do if (k neq i) and (k neq j) then << a:=sortpermuli({i,j,k}); flp:=car a; a:=cdr a; qlicp:=qli; while qlicp and (caar qlicp neq a) do qlicp:=cdr qlicp; if qlicp then f:=cdar qlicp else << f:=newfct(fname_,vl,nfct_); nfct_:=add1 nfct_; ftem_:=fctinsert(f,ftem_); qli:=cons(a . f,qli) >>; f:={'DF,f,nth(xlist,k)}; if flp then f:={'MINUS,f}; qij:=cons(f,qij) >>$ if null qij then <<qij:=newfct(fname_,setdiff(vl,xlist),nfct_); nfct_:=add1 nfct_; ftem_:=fctinsert(qij,ftem_)>> else if cdr qij then qij:=reval cons('PLUS,qij) else qij:=car qij; resu:=cons(qij,resu) >>$ return resu end$ symbolic procedure updt_curl(h2,rmdr,fl,done_xlist,x,cdrxlist,n,k)$ % a subroutine of int_curl begin scalar i,h4,h5,h6,h7,rmdr,y,pint,succ$ if (not zerop reval reval {'DF,rmdr,x}) then << if print_ then <<terpri()$write"No success."$terpri()>>$ succ:=nil >> else << succ:=t; if done_xlist then << % there is one computed curl component to be updated % integration wrt done_xlist h7:=intcurrent2(rmdr,fl,cons('LIST,done_xlist)); rmdr:=caddr h7; h7:=cdadr h7; % update the already computed h2-components with the new h7-comp. h4:=nil; h5:=-1; for i:=1:(k-1) do << h5:=add1 h5; for h6:=1:(n-k) do <<h4:=cons(car h2,h4);h2:=cdr h2>>; h4:=cons({'DIFFERENCE,car h2,car h7},h4); h2:=cdr h2; h7:=cdr h7; for h6:=1:h5 do <<h4:=cons(car h2,h4);h2:=cdr h2>> >>; h2:=reverse h4; >>$ % now generalized integration of the remainder if zerop rmdr then pint:=cons(0,nil) else << y:=if cdrxlist then car cdrxlist else car done_xlist; fnew_:=nil$ pint:=partint(rmdr,fl,vl_,y,genint_); % genint is max number of new terms if null pint then succ:=nil else for each h4 in fnew_ do ftem_:=fctinsert(h4,ftem_) >> >>; return if null succ then nil else cons(h2,pint) % pint=cons(generalized integral of rmdr,list of new eqn) end$ symbolic procedure int_curl(pli,fl,xlist,vl)$ % given a vector p^i satisfying p^i,_i=0, find q^{ij}=-q^{ji} % such that p^i=q^{ij},j % car result: (q^{12}, q^{13},.., q^{1n}, q^{23},.., q^{2n},.., q^{(n-1)n}) % each q^{ij} comes with r^{ijk},k % cdr result: list of new conditions in fewer variables % works only if identically satisfied, not modulo some equations % vl is the list of all relevant variables % during computation is h2 = % (q^{kn},.., q^{k(k+1)},q^{(k-1)n},.., q^{(k-1)k},.., % q^{2n},.., q^{23}, q^{1n},.., q^{13},q^{12}} ) begin scalar h1,h2,h3,resu,newpli,xcp,done_xlist,n,k,ok,neweq,ftem_bak$ % conversion from algebraic mode lists to lisp lists: pli:=cdr pli$ xlist:=cdr xlist$ vl:=cdr vl; xcp:=xlist$ n:=length(xlist); k:=0; ok:=t; ftem_bak:=ftem_; if n=1 then return {'LIST,reval cons('LIST,pli),{'LIST}}$ while cdr pli and ok do << k:=add1 k; % the integration has to be done first wrt cdr xlist. The resulting % curl will be used to change the remining pli to be integrated h3:=intcurrent2(reval car pli,fl,cons('LIST,cdr xlist)); pli:=cdr pli; h1:=cdadr h3; h3:=reval reval caddr h3; % h3 now = the remainder of the integration wrt cdr xlist if not zerop h3 then << % here the integration wrt the done_xlist. These curl updates will % not be used to update pli, because df(h3,car xlist)=0 is assumed h3:=updt_curl(h2,h3,fl,done_xlist,car xlist,cdr xlist,n,k)$ if null h3 then ok:=nil else << % generalized integration of the remainder neweq:=append(cddr h3,neweq); h2:=car h3; h1:=cons({'PLUS,car h1,cadr h3},cdr h1); % because of cdr xlist neq nil here q^{k(k+1)} is updated >> >>$ if ok then << % In the first round h2 is naturally nil --> use ok for test % append (q^{kn},.., q^{k(k+1)}) and h2 h2:=append(reverse h1,h2); % update the remaining pli to be integrated newpli:=nil; while h1 do << newpli:=cons({'PLUS,{'DF,car h1,car xlist},car pli},newpli); h1:=cdr h1; pli:=cdr pli >>; pli:=reverse newpli >>; done_xlist:=cons(car xlist,done_xlist); xlist:=cdr xlist >>$ if ok then << pli:=reval car pli; % to get the remainder of the last component of pli integrated if pli neq 0 then << k:=k+1; h3:=updt_curl(h2,pli,fl,done_xlist,car xlist,nil,n,k)$ if null h3 then ok:=nil else << neweq:=append(cddr h3,neweq); h2:=car h3; h2:=cons({'DIFFERENCE,car h2,cadr h3},cdr h2) % because of null xlist here car h2=q^{n-1,n} is updated >>$ >> >>; if null ok then << % drop all new functions h1:=setdiff(ftem_,ftem_bak); for each h2 in h1 do drop_fct(h2)$ ftem_:=ftem_bak >> else << h1:=curlconst(xcp,vl)$ while h1 do << resu:=cons({'PLUS,car h2,car h1},resu); h1:=cdr h1; h2:=cdr h2; >>$ return {'LIST,reval cons('LIST,resu),cons('LIST,neweq)} >> end$ endmodule$ end$