Artifact d8aadfac8ed3618ba56780ca6c4f9b6124fc90a567635c6a2f25a5e3f1872907:
- Executable file
r37/packages/crack/conlaw2.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: 25697) [annotate] [blame] [check-ins using] [more...]
% CONLAW version 2, to calculate conservation laws of systems % of PDEs by calculating characteristic functions % by Thomas Wolf, September 1997 %---------------------------------------------------------------------- symbolic fluid '(print_ logoprint_ potint_ facint_ adjust_fnc)$ %symbolic fluid '(div_p)$ %------------- symbolic procedure newil(il,mo,nx)$ if (null il) or (length il<mo) then cons(1,il) else if car il<nx then cons(add1 car il,cdr il) else <<while il and (car il = nx) do il:=cdr il; if null il then nil else cons(add1 car il,cdr il)>>$ %------------- symbolic procedure sortli(l)$ % sort a list of numbers begin scalar l1,l2,l3,m,n$ return if null l then nil else << n:=car l$ l2:=list car l$ l:=cdr l$ while l do << m:=car l$ if m>n then l1:=cons(car l,l1) else if m<n then l3:=cons(car l,l3) else l2:=cons(car l,l2)$ l:=cdr l >>$ append(sortli(l1),append(l2,sortli(l3))) >> end$ %------------- %symbolic operator combi$ symbolic procedure combi(ilist)$ % ilist is a list of indexes (of variables of a partial derivative) % and returns length!/k1!/k2!../ki! where kj! is the multiplicity of j. begin integer n0,n1,n2,n3; n1:=1; % ilist:=cdr ilist; while ilist do <<n0:=n0+1;n1:=n1*n0; if car ilist = n2 then <<n3:=n3+1; n1:=n1/n3>> else <<n2:=car ilist; n3:=1>>; ilist:=cdr ilist>>; return n1 end$ %------------- symbolic procedure derili(il)$ % make a derivative index list from a list of numbers if null il then nil else begin scalar h1,h2,h3$ h1:=sortli(il); while h1 do << h2:=reval algebraic mkid(!`,lisp car h1); h3:=if h3 then mkid(h2,h3) else h2; h1:=cdr h1 >>; return h3 end$ %------------- symbolic operator intcurrent1$ symbolic procedure intcurrent1(divp,ulist,xlist,dulist, nx,nde,eqord,densord)$ % computes a list in reverse order from which the conserved % current is computed through integration begin scalar ii,il,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11,contrace,u, nega,pii,mo,pl$ %contrace:=t; ulist:=cdr ulist; xlist:=cdr xlist; mo:=if eqord>densord then eqord-1 else densord-1$ pl:=for ii:=1:nx collect << % the components W^i il:=nil; pii:=nil; repeat << h11:=cons(ii,il); h1:=derili(h11); h11:=combi(sortli(h11)); if contrace then <<write"==== ii=",ii," il=",il," h1=",h1," h11=",h11;terpri()>>; h2:=il; h3:=nil; if null h2 then h4:=list(nil . nil) else << h4:=list(car h2 . nil); while h2 do << h3:=cons(car h2,h3);h2:=cdr h2; h4:=cons((if h2 then car h2 else nil ) . derili(h3),h4)$ >>; >>; if contrace then <<write"h3=",h3," h4=",h4;terpri()>>; for e1:=1:nde do << % for each function u u:=nth(ulist,e1); h5:=mkid(u,h1); if contrace then <<write"h5=",h5;terpri()>>; h6:=nil; if contrace then <<write"h6-1=", list('DF,divp,h5); terpri()>>; h6:=cons(reval list('QUOTIENT,list('DF,divp,h5),h11), h6); if nde=1 then h6:=car h6 else h6:=cons('PLUS,h6); if contrace then <<write"h6=",h6;terpri()>>; nega:=nil; h7:=h4; % h7 is in such an order that the first term is always positive while h7 and not zerop h6 do << h8:=car h7; h7:=cdr h7; h9:=car h8; h8:=cdr h8; if contrace then <<write if nega then "--" else "++", " h8=",h8," h9=",h9;terpri()>>; if contrace then <<write"h9-2=",h9;terpri()>>; if h9 then h6:=totdif(h6,nth(xlist,h9),h9,dulist); if contrace then <<write"h6-3=",h6;terpri()>>; h10:=if h8 then mkid(u,h8) else u; if contrace then <<write"h10=",h10;terpri()>>; h10:=list('TIMES,h6,h10); if nega then h10:=list('MINUS,h10); if contrace then algebraic write">>>> h10=",h10; pii:=cons(h10,pii)$ nega:=not nega; >> >>; % for each function u il:=newil(il,mo,nx); >> until null il; pii:=reval if length pii=1 then car pii else cons('PLUS,pii); if contrace then algebraic write"pii-end=",pii; pii >>; % for all ii return cons('LIST,pl) end$ % of intcurrent1 %------------- symbolic operator intcurrent2$ symbolic procedure intcurrent2(divp,ulist,xlist,nondiv)$ % computes the conserved current P_i from the divergence through % partial integration begin scalar h2,h3,h4,h5,h6,h7,h8,e2,e3; % repeated partial integration to compute P_i ulist :=cdr reval ulist; xlist :=cdr reval xlist; h4:=list xlist; % dequ is here a list containing only the conserved density % and flux to drop commen factors repeat << e3:=divp; h3:=car h4; % h3 is the list of variables is a spec. order h4:=cdr h4; h5:=for e2:=1:length h3 collect 0; % h5 is old list of the conserved quantities h8:=0; % h8 counts integrations wrt. all variables repeat << % integrate several times wrt. all variables h8:=h8+1; h2:=h3; % h2 is a copy of h3 h6:=nil; % h6 is new list of the conserved quantities h7:=nil; % h7 is true if any integration was possible while h2 neq nil do << % integrating wrt. each variable e2:=intpde(e3,ulist,h3,car h2,t); e3:=cadr e2; if not zerop e2 then h7:=t; h6:=cons(list('PLUS,car e2,car h5),h6); h5:=cdr h5; h2:=cdr h2 >>; h5:=reverse h6; >> until (h7=nil) % no integration wrt. no variable was possible or (e3=0) % complete integration or (h8=10); % max. 10 integrations wrt. all variables >> until (e3=0) or (h4=nil); return if e3=0 then cons('LIST, h5) else nondiv % end of the computation of the conserved current P % result is a list of P_i end$ % of intcurrent2 %------------- algebraic procedure conlaw2(problem,runmode)$ begin scalar contrace,eqlist,ulist,xlist,dequ,cllist,divlist, sb,densord,flist,eqord,maxord,dulist,revdulist,vl,expl, deplist,e1,e2,e3,n,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11, condi,soln,adjust_old,potold,adjustold,udens,gensepold, inequ0,inequ,logoold,treqlist,fl,facold,u,lhslist,cpu, gc,cpustart,gcstart,nontriv,cf0,rtnlist,paralist,solns, found,clcopy,extraline,nondiv,nx,nde,nonconstc, mindensord,mindensord0,maxdensord,rules$ lisp <<adjustold:=adjust_fnc; adjust_fnc:=t; logoold:=logoprint_; logoprint_:=t; potold:=potint_; potint_:=t; facold:=facint_; facint_:=1000>>; cpustart:=lisp time()$ gcstart:=lisp gctime()$ % contrace:=t; %--- extracting input data eqlist:= maklist first problem; ulist := maklist second problem; xlist := maklist third problem; nx:=length xlist; nde:=length eqlist; if contrace then write"eqlist=",eqlist, " ulist=",ulist," xlist=",xlist; mindensord:=part(runmode,1)$ maxdensord:=part(runmode,2)$ expl :=part(runmode,3)$ flist :=part(runmode,4)$ inequ0 :=part(runmode,5)$ problem:=runmode:=0; %--- initial printout lisp(if logoprint_ then <<terpri()$ write "--------------------------------------------------", "------------------------"$ terpri()$terpri()$ write "This is CONLAW2 - a program for calculating conservation", " laws of DEs"; terpri() >> else terpri()); if nde = 1 then write "The DE under investigation is :" else write "The DEs under investigation are :"; for each e1 in eqlist do write e1; write "for the function(s): ",ulist; write"======================================================"$ %--- lhslist is the list of derivatives to substitute %--- is needed in order not to let the Q depend on them and their deriv. eqlist:=reverse eqlist; lhslist:={}; for each h1 in eqlist do lhslist:=cons(lhs h1,lhslist); %--- bringing eqlist and ulist in line h1:={}; h2:={}; if nde neq length ulist then n:=t else for each e1 in eqlist do << e2:=part(lhs e1,1); if freeof(ulist,e2) then n:=t else h1:=cons(e2,h1); h2:=cons(lhs e1 - rhs e1, h2) >>; if n then rederr("The lists of equations and functions are not consistent!") else <<ulist:=h1; eqlist:=h2>>; if contrace then write"ulist=",ulist," eqlist=",eqlist; %--- initializations to be done only once rtnlist:={}; nondiv:=lisp gensym(); % as a marker if p-computation was not succ. %------ the list of parameters of the equation to be determined paralist:={}; for each e1 in flist do if not freeof(eqlist,e1) then paralist:=cons(e1,paralist); %------ determination of the order of the input equations eqord:=0; mindensord0:=mindensord; for each e1 in eqlist do for each e2 in ulist do << h1:=totdeg(e1,e2); if h1>eqord then eqord:=h1 >>; for n:=1:nde do << h1:=mkid(q_,n); if not lisp(null get(mkid('q_,n),'avalue)) then << for each e2 in ulist do << h2:=totdeg(h1,e2); if h2>eqord then eqord:=h2; if h2>mindensord then mindensord:=h2 >>; cf0:=t; >> >>; if contrace then write"eqord=",eqord; if maxdensord<mindensord then maxdensord:=mindensord; %------ all transformations into jet-space sb:=subdif1(xlist,ulist,eqord)$ if contrace then write"sb=",sb; treqlist:=eqlist; for each e1 in sb do treqlist:=sub(e1,treqlist); for each e1 in sb do lhslist:=sub(e1, lhslist); if contrace then write"treqlist=",treqlist, " lhslist=",lhslist; if cf0 then for n:=1:nde do << h1:=mkid(q_,n); if not lisp(null get(mkid('q_,n),'avalue)) then << for each e1 in sb do h1:=sub(e1,h1); lisp(mkid('q_,n)):=h1; >> >>; for each e1 in sb do inequ0:=sub(e1,inequ0); %--- investigate conservation laws of increasing order for densord:=mindensord:maxdensord do << nodepnd(ulist); cpu:=lisp time()$ gc:=lisp gctime()$ if cf0 then lisp<<write"A special ansatz of order ",densord, " for the characteristic"$terpri()$ write"function(s) is investigated.";terpri() >> else lisp<< write"Currently conservation laws with characteristic"; terpri(); write"function(s) of order ",densord," are determined"; terpri(); write"======================================================"$ >>; %--- repeated initializations %--- maxord is maximal derivative in condition maxord:=eqord % from the total derivatives + 1 % for safety + if eqord>densord then eqord else densord$ %######## possibly to be increased due to substitutions if contrace then write"maxord=",maxord; if {}=fargs first ulist then for each e1 in ulist do depnd(e1,{xlist}); sb:=subdif1(xlist,ulist,maxord)$ nodepnd ulist; if contrace then write"sb=",sb; dulist:=ulist . reverse for each e1 in sb collect for each e2 in e1 collect rhs e2; sb:=0; revdulist:=reverse dulist; % dulist with decreasing order udens:=part(dulist,densord+1); % derivatives of order densord vl:=for each e1 in dulist join e1; if contrace then write"vl=",vl," udens=",udens; if not flist then fl:={} else fl:=flist; %--- initializing characteristic functions cf, the list of functions fl deplist:=ulist . for n:=1:densord collect listdifdif2(lhslist,part(dulist,n+1)); if expl then deplist:=xlist . deplist; deplist:=reverse deplist; cf:={}; for n:=1:nde do << h1:=mkid(q_,n); if lisp(null get(mkid('q_,n),'avalue)) then << nodepnd({h1}); depnd(h1, deplist); fl:=cons(h1,fl); >>; cf:=cons(h1,cf); >>; cf:=reverse cf; if contrace then write"fl=",fl; if contrace then lisp (write" depl*=",depl!*); %--- generation of the conditions condi:={}; for each u in ulist do << if contrace then write"function=",u; h1:=treqlist; h2:=cf; h3:=0; while h1 neq {} do << % sum over all equations if contrace then write"equation :",first h1; for each e1 in vl do % sum over u and all its derivatives if lisp(reval algebraic(u) = car combidif algebraic(e1)) then << % for u and all its derivatives e2:=df(first h1, e1); if e2 neq 0 then << if contrace then write"e1=",e1; dequ:=first h2 * e2; e2:=1; for each e3 in lisp cons('LIST,cdr combidif(algebraic e1)) do <<dequ:=totdif(dequ,part(xlist,e3),e3,dulist)$ e2:=-e2; if contrace then write"dequ=",dequ," e3=",e3>>; h3:=h3+e2*dequ; if contrace then write"h3=",h3; >>; >>; h1:=rest h1;h2:=rest h2 >>; condi:=cons(h3,condi) >>; if contrace then write"condi=",condi; %--- generating a substitution list %--- at first using the equations themselves sb:={}; rules:={}; h1:=treqlist; h2:=lhslist; h4:=nil; % h4 is list of undifferentiated substitutions while h1 neq {} do << h3:=first h2; h5:=h3-first h1; rules:=cons(h3 => h5,rules)$ lisp(e3:=combidif h3); % extracts the list of derivatives: u`1`1`2 --> (u, 1, 1, 2) lisp(h4:=cons(list(car e3, cdr e3, h5), h4))$ % function name of h3, derivative list of h3, value of h3 h1:=rest h1; h2:=rest h2; >>; %--- then their derivatives for each e1 in vl do lisp << e1:=reval e1; % is e1 a derivative of any of the undifferentiated substitutions? h1:=h4; while h1 neq nil do << h3:=comparedif2(caar h1, cadar h1, e1); if (h3 neq nil) and (h3 neq 0) then algebraic << h3:=lisp(cons('LIST,h3)); dequ:=lisp caddar h1; % rhs which is to be differentiated for each n in h3 do dequ:=totdif(dequ,part(xlist,n),n,dulist); % new highest derivatives should be added to vl afterwards % if lower derivatives are substituted by higher derivatives rules:=cons(e1 => dequ,rules)$ lisp(h1:=nil) >> else lisp(h1:=cdr h1); >> >>; if contrace then write"rules=",rules; let rules; condi:=condi; clearrules rules$ if contrace then write"condi=",condi; vl:=append(xlist,vl); % now the full list inequ:=inequ0; %--- inequ is to stop crack if order of cf is too low if (densord neq 0) and ((cf0=nil) or (mindensord0 neq 0)) then << % investigation should stop if % cf is independent of highest order derivatives dequ:=0; for each e1 in cf do << h1:=udens; while h1 neq {} do << dequ:=dequ+df(e1,first h1)*(lisp gensym()); h1:=rest h1 >>; >>; inequ:=cons(dequ,inequ) >>; let rules; inequ:=inequ; clearrules rules$ rules:=0; if contrace then write"inequ=",inequ; %--- freeing some space sb:=revdulist:=deplist:=e1:=e2:=e3:= n:=h1:=h2:=h3:=soln:=u:=dequ:=0; %--- the real calculation if lisp(!*time) then write "time to formulate condition: ", lisp time() - cpu, " ms GC time : ", lisp gctime() - gc," ms"$ solns:=crack(condi,inequ,fl,vl); %--- postprocessing lisp terpri()$ found:=nil; while solns neq {} do << divlist:={}; cllist:={}; soln:=first solns; solns:=rest solns; condi:=first soln; cfcopy:=sub(second soln,cf); h1:=third soln; if contrace then << write"cfcopy=",cfcopy; write"soln=",soln; write"third soln=",h1; >>; fl:={}; h2:={}; for each e1 in h1 do << if not freeof(condi,e1) then fl:=cons(e1,fl); % fl to output remaining conditions later if freeof(paralist,e1) then h2:=cons(e1,h2) >>; h1:=parti_fn(h2,condi)$ if contrace then write"h1(partitioned)=",h1; extraline:=nil; nonconstc:={}; while h1 neq {} do << e1:=first h1;h1:=rest h1; for each h4 in e1 do if fargs h4 neq {} then << nonconstc:=cons(h4,nonconstc); lisp << write"The function "$ fctprint list reval h4$ write" is not constant!"; extraline:=t; terpri() >> >>; dequ:=0; % to compute rhs h2:=treqlist; % " if paralist then h2:=sub(second soln,h2); % " if contrace then write"h2=",h2; % " nontriv:=nil; h3:=for each e2 in cfcopy collect << e3:=for each h4 in e1 sum fdepterms(e2,h4); dequ:=dequ+e3*first h2; h2:=rest h2; % computes rhs if e3 neq 0 then nontriv:=t; e3 >>; if nontriv then << found:=t; cllist:=cons(<<if contrace then write"h3-1=",h3," dequ=",dequ; sb:=absorbconst(h3,e1)$ if (sb neq nil) and (sb neq 0) then << h3:=sub(sb,h3); dequ:=sub(sb,dequ) >>; if contrace then write"h3-2=",h3," dequ=",dequ; if (length(e1)=1) and (fargs first e1 = {}) then << h4:=first e1; dequ:=sub(h4=1,dequ); sub(h4=1,h3) >> else h3 >>, cllist); divlist:=cons(dequ,divlist) >> >>; if contrace then << write"characteristic functions found so far:"; write cllist; >>$ if condi neq {} then << write"There are remaining conditions: ", condi; lisp << write"for the functions: "; fctprint cdr reval algebraic fl;terpri(); write"Corresponding CLs might not be shown below as they"; terpri()$write"could be of too low order.";terpri()>>; extraline:=t; >>; if extraline then lisp << write"======================================================"$ terpri() >>; %--- Dropping conservation laws of too low order if (densord > 0) and ((cf0=nil) or (mindensord0 neq 0)) then << h1:={}; h2:={}; for each e1 in cllist do << h5:=udens; while (h5 neq {}) and freeof(e1,first h5) do h5:=rest h5; if h5 neq {} then << h1:=cons(e1,h1); h2:=cons(first divlist,h2) >>; divlist:=rest divlist; >>; cllist:=h1; divlist:=h2 >>; if contrace then write"cllist=",cllist; if cllist neq {} then << %--- Below h1 is the list of W^i in the Anco/Bluman formula h1:=for e1:=1:(length cllist) collect intcurrent1(part(divlist,e1),ulist,xlist,dulist,nx,nde, eqord,densord); %--- Backsubstitution of e.g. u`1`1 --> df(u,x,2) for each e1 in ulist do depnd(e1,{xlist}); on evallhseqp; sb:=subdif1(xlist,ulist,maxord)$ sb:=for each e1 in sb join for each e2 in e1 collect(rhs e2 = lhs e2); off evallhseqp; cllist:=sub(sb,cllist); h1:=sub(sb,h1); %--- lambda integration of h1 to compute P_i h2:=lisp gensym()$ for each e1 in ulist do h1:=sub(e1=h2*e1,h1); if not lisp(freeof(h1,'SUB)) then h1:={} else h1:=for each e1 in h1 collect << % i.e. for each cl h10:=sub(sb,first divlist); divlist:=rest divlist; % at first try direct integration to compute p h9:=intcurrent2(h10,append(nonconstc,ulist),xlist,nondiv); % if no success use lambda-integration if h9=nondiv then << h8:=t; % whether intcurrent1 is still ok %--- at first the term h10 = T^i/x^i in conca.tex for each e2 in ulist do << if h8 then h10:=err_catch_sub(e2,0,h10); if h10 eq nil then h8:=nil >>$ if contrace then write"h10-1=",h10$ if h8 and (h10 neq 0) then << for each e2 in xlist do << if h8 then h10:=err_catch_sub(e2,h2*e2,h10); if h10 eq nil then h8:=nil >>$ if h8 then << if contrace then write"h10-2=",h10$ h10:=int(h10*h2**(nx-1),h2); if contrace then write"h10-3=",h10$ %--- the following is to catch errors in: %--- sub(h2=1,h10)-sub(h2=0,h10) h6:=err_catch_sub(h2,1,h10); if contrace then write"h6=",h6$ if h6 eq nil then h7:=nil else h7:=err_catch_sub(h2,0,h10); if contrace then write"h7=",h7$ if h7 eq nil then h8:=nil else h10:=h6-h7 >> >>$ if contrace then write"h10-4=",h10$ h4:={}; % h4 becomes the inverse list of P^i h11:=0; while h8 and (e1 neq {}) do << h11:=h11+1; e2:=first e1; e1:=rest e1; if contrace then write"e2=",e2$ h3:=int(e2/h2,h2); if contrace then write"h3-1=",h3$ %--- the following is to catch errors in: %--- sub(h2=1,h3)-sub(h2=0,h3) h6:=err_catch_sub(h2,1,h3); if h6 eq nil then h7:=nil else h7:=err_catch_sub(h2,0,h3); if h7 eq nil then h8:=nil else h4:=cons(h6-h7+h10*part(xlist,h11),h4) >>; if h8 then h9:=reverse h4 >>; h9 >>; if contrace then write"h1-1=",h1$ if h1={} then << lisp << write"The conserved quantities could not be found."$ terpri() >>$ if condi neq {} then lisp << write"For that the remaining conditions should be solved."; terpri() >>; lisp << write"The adjoined symmetries are:"$terpri() >>$ for each e1 in cllist do write e1$ >>$ if contrace then << write"h1=",h1;write"cllist=",cllist;write"eqlist=",eqlist >>; while h1 neq {} do << h2:=first h1; h3:=first cllist; rtnlist:=cons({h3,h2},rtnlist); %--- conditions on parameters if paralist neq {} then for each e2 in second soln do if not freeof(paralist,lhs e2) then <<write e2,",";lisp(terpri())>>$ %--- the conservation laws %--- Test whether actually only an adjoint symmetry has been %--- computed and not a conservation law h4:=eqlist; if paralist neq {} then h4:=sub(second soln,h4); h8:=0; if h2 neq nondiv then << h5:=h4; for each e1 in h3 do << h8:=h8 + e1*(first h5)$ h5:=rest h5 >>$ for e1:=1:nx do << h8:=h8-df(part(h2,e1),part(xlist,e1))$ % for test purposes >>; if h8 neq 0 then h2:=nondiv >>; if h2=nondiv then write"Adjoint symmetry:" else write"Conservation law:"; while h3 neq {} do << if length h3 < length first cllist then write "+"$ write"(",first h3,") * (",first h4,")"$ h3:=rest h3; h4:=rest h4 >>$ if h2=nondiv then << lisp << write"could not be written as a divergence but solves the"$ terpri()$ write"adjoint symmetry condition and therefore represents"$ terpri()$ write"an adjoint symmetry."$ terpri()$ >>$ if (h8 neq 0) and (condi neq {}) then << write"Please check: if the remaining conditions guarantee "$ write" 0 = ",h8$ write"then the conservation law is"$ write" 0 = "$ for e1:=1:nx do << write"df( ",part(h2,e1),", ",part(xlist,e1)," )"$ if e1<nx then write "+"$ >> >>$ >> else << write"= "$ for e1:=1:nx do << write"df( ",part(h2,e1),", ",part(xlist,e1)," )"$ if e1<nx then write "+"$ >> >>; h1:=rest h1; cllist:=rest cllist; write"======================================================"$ >>$ >>; % if cllist neq {} then << nodepnd(ulist); >>; % while solns neq {} do << if found=nil then << write"There is no conservation law of this order."; write"======================================================"$ >> >>; % for densord:=mindensord:maxdensord if fargs first ulist = {} then for each e1 in ulist do depnd(e1,{xlist}); if lisp(!*time) then write "time to run conlaw2: ", lisp time() - cpustart, " ms GC time : ", lisp gctime() - gcstart," ms"$ lisp <<adjust_fnc:=adjustold; logoprint_:=logoold; %gensep_:=gensepold; potint_:=potold; facint_:=facold>>; return rtnlist end$ % of conlaw2 end$