Artifact bdf7e7d01a8c04616bd269301d59689f98908d1bd83e37bb658f9025abc8ab3d:
- Executable file
r38/packages/crack/conlaw3.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: 15355) [annotate] [blame] [check-ins using] [more...]
% CONLAW version 3, to calculate conservation laws of systems of PDEs % by calculating characteristic functions and conserved currents % by Thomas Wolf, June 1999 %---------------------------------------------------------------------- symbolic fluid '(print_ logoprint_ potint_ facint_ adjust_fnc quasilin_rhs)$ %------------- algebraic procedure conlaw3(problem,runmode)$ begin scalar contrace,eqlist,ulist,xlist,dequ,cllist,pllist, sb,densord,flist,eqord,maxord,dulist,revdulist,vl,expl, deplist,e1,e2,e3,n,h1,h2,h3,h4,h5,h6,condi,soln, adjust_old,potold,adjustold,udens,gensepold,inequ0, inequ,logoold,treqlist,fl,facold,u,nodep,cpu,gc, cpustart,gcstart,found,cf0,rtnlist,solns,nontriv, extraline,cf,cfcopy,nx,nde,mindensord,mindensord0, maxdensord,absmaxord,nonconstc; backup_reduce_flags()$ 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:= reverse 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 CONLAW3 - 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 reverse eqlist do write e1; lisp<<terpri()$write "for the function(s): "$ fctprint cdr reval ulist;terpri()>>$ write"======================================================"$ %--- nodep is a list of derivatives the Q do not depend on nodep:=first lhsli(eqlist)$ %--- Here comes a test that lhs's are properly chosen chksub(eqlist,ulist)$ %--- Checking whether an ansatz for characteristic functions %--- has been made, then denominator of equations is not dropped for n:=1:nde do if not lisp(null get(mkid('q_,n),'avalue)) then cf0:=t; eqlist:=reverse for each e1 in eqlist collect if part(e1,0)=EQUAL then if cf0 then lhs e1 - rhs e1 else num(lhs e1 - rhs e1) else if cf0 then e1 else num e1; if contrace then write"ulist=",ulist," eqlist=",eqlist; %--- initializations to be done only once rtnlist:={}; %------ 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; for each e1 in eqlist do for each e2 in ulist do << % h1:=totordpot(e1,e2); % if car h1>eqord then <<eqord:=car h1;eqlddeg:=cdr h1>> else % if car h1=eqord then % if cdr h1>eqlddeg then eqlddeg:=cdr h1 h1:=totdeg(e1,e2); if h1>eqord then eqord:=h1 >>; h3:=eqord; mindensord0:=mindensord$ 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; >> >>; for n:=1:nx do << % if the index of p_ should be a number then use n instead of h4 h4:=lisp(reval algebraic(part(xlist,n))); h1:=mkid(p_,h4); if not lisp(null get(mkid('p_,h4),'avalue)) then << for each e2 in ulist do << h2:=totdeg(h1,e2); if h2>eqord then eqord:=h2; if (h2>=h3) and (mindensord<=h2) then mindensord:=h2+1 >>; cf0:=t; >> >>; if maxdensord<mindensord then maxdensord:=mindensord; if contrace then write"eqord=",eqord," cf0=",cf0; %------ all transformations of input data 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); nodep:=sub(e1,nodep) >>; if contrace then write"treqlist=",treqlist, "nodep=",nodep; 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; >> >>; if cf0 then for n:=1:nx do << % if the index of p_ should be a number then use n instead of h4 h4:=lisp(reval algebraic(part(xlist,n))); h1:=mkid(p_,h4); if not lisp(null get(mkid('p_,h4),'avalue)) then << for each e1 in sb do h1:=sub(e1,h1); lisp(mkid('p_,h4)):=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 order of a derivative on the right hand side %--- absmaxord is maximal order of a derivative occuring at least % temprorily in the condition (in the null divergence) % If q*delta is linear in highest derivatives then P is of one % order lower otherwise P is of same order but one degree less maxord:=if eqord>densord then eqord else densord; absmaxord:=if lisp(null quasilin_rhs) then maxord+1 else maxord; if contrace then << write"maxord=",maxord; write"absmaxord=",absmaxord; >>; if {}=fargs first ulist then for each e1 in ulist do depnd(e1,{xlist}); sb:=subdif1(xlist,ulist,absmaxord)$ 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, %--- the conserved current pl and the condition condi condi:=0; deplist:=lisp(cons('LIST,setdiff(cdr ulist,cdr nodep))) . for n:=1:densord collect listdifdif2(nodep,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); condi:=condi+h1*part(treqlist,n); >>; cf:=reverse cf$ deplist:=for h3:=0:(absmaxord-1) collect part(dulist,h3+1); if expl then deplist:=xlist . deplist; deplist:=reverse deplist; pl:={}; for n:=1:nx do << % if the index of p_ should be a number then use n instead of h4 h4:=lisp(reval algebraic(part(xlist,n))); h1:=mkid(p_,h4); if lisp(null get(mkid('p_,h4),'avalue)) then << nodepnd({h1}); depnd(h1, deplist); fl:=h1 . fl; >>; pl:=cons(h1,pl); condi:=condi-totdif(h1,h4,n,dulist) >>; sb:=0; if contrace then write"fl=",fl," cf=",cf," pl=",pl; if contrace then lisp (write" depl*=",depl!*); if contrace then write"condi=",condi; vl:=reverse 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 << % for the investigation to 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 intern gensym()); h1:=rest h1 >>; >>; inequ:=cons(dequ,inequ) >>; if contrace then write"inequ=",inequ; condi:={condi}; if (not lisp(null get('cl_condi,'avalue))) and (part(cl_condi,0)=LIST) then condi:=append(condi,cl_condi)$ %--- freeing some space sb:=dulist:=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()$ pllist:={}; cllist:={}; found:=nil; while solns neq {} do << % for each solution (if param. are determ.) soln:=first solns; solns:=rest solns; condi:=first soln; % filtering out conservation laws found in the previous run cfcopy:=sub(second soln,cf); % any non-trivial conservation law? h1:=0; for each h2 in cfcopy do if h2 neq 0 then h1:=1; if h1 neq 0 then << pl:=sub(second soln,pl); if contrace then write"cfcopy=",cfcopy," pl=",pl; h1:=third soln; if contrace then 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 << % for each potential conservation law % h1 is the list of lists of constants/functions % depending on each other h2:=first h1;h1:=rest h1; if contrace then write"h2=",h2; if contrace then write"cfcopy=",cfcopy; nontriv:=nil; %--- is the constant/function in the characteristic functions? h3:=for each e2 in cfcopy collect << e3:=for each e1 in h2 sum fdepterms(e2,e1); if e3 neq 0 then nontriv:=t; e3 >>; if nontriv then << for each e1 in h2 do if fargs e1 neq {} then lisp << nonconstc:=cons('LIST,cons(reval e1,cdr nonconstc)); write"The function "$ fctprint list reval e1$ write" is not constant!"; extraline:=t; terpri() >>; %--- the current h4:=reverse for each e2 in pl collect for each e1 in h2 sum fdepterms(e2,e1); if contrace then write"h3-1=",h3," h4-1=",h4; sb:=absorbconst(h3,h2)$ if (sb neq nil) and (sb neq 0) then << h3:=sub(sb,h3); h4:=sub(sb,h4) >>; if contrace then write"h3-2=",h3," h4-2=",h4; if (length(h2)=1) and (fargs first h2 = {}) then << e1:=first h2; h4:=sub(e1=1,h4); h3:=sub(e1=1,h3) >>; h5:=udens; if (densord > 0) and ((cf0=nil) or (mindensord0 neq 0)) then while (h5 neq {}) and freeof(h3,first h5) do h5:=rest h5; if h5 neq {} then << % h3 is of order densord cllist:=cons(h3,cllist); pllist:=cons(h4,pllist) >> >> >>; 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() >>; for each e1 in ulist do depnd(e1,{xlist}); if contrace then write"cllist2=",cllist," pllist2=",pllist$ on evallhseqp; sb:=subdif1(xlist,ulist,maxord)$ sb:=for each e1 in sb join for each e2 in e1 collect(rhs e2 = lhs e2); if contrace then write"sb=",sb$ off evallhseqp; cllist:=sub(sb,cllist); if contrace then write"cllist3=",cllist$ pllist:=sub(sb,pllist); if contrace then write"pllist3=",pllist$ % if nx=2 then % pllist:=simppl(pllist,ulist,first xlist,second xlist)$ if contrace then << write"cllist3=",cllist; write"pllist3=",pllist; write"eqlist=",eqlist; write"xlist=",xlist >>; while pllist neq {} do << found:=t; write"Conservation law:"; h2:=first pllist; 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 h4:=eqlist; if paralist then h4:=sub(second soln,h4); print_claw(h4,h3,h2,xlist)$ %--- factoring out diff operators? h6:={}; for each h5 in nonconstc do if not freeof(h3,h5) then h6:=cons(h5,h6); if (h6 neq {}) and (h2 neq nondiv) then partintdf(h4,h3,h2,xlist,h6,vl,sb); write"======================================================"$ pllist:=rest pllist; cllist:=rest cllist; >>$ >>; >>; 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 conlaw3: ", lisp time() - cpustart, " ms GC time : ", lisp gctime() - gcstart," ms"$ lisp <<adjust_fnc:=adjustold; logoprint_:=logoold; potint_:=potold; facint_:=facold>>; recover_reduce_flags()$ return rtnlist end$ % of conlaw3 end$