Artifact 506c03b81e4c4e6d48b22611d90d408387b320257a38a2520e2dcb628bea7a23:
- Executable file
r37/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: 14481) [annotate] [blame] [check-ins using] [more...]
comment consider many solutions of crack ; % CONLAW version 3, to calculate conservation laws of systems of PDEs % by calculating characteristic functions and conserved currents % by Thomas Wolf, September 1997 %---------------------------------------------------------------------- symbolic fluid '(print_ logoprint_ potint_ facint_ adjust_fnc)$ %------------- 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,h7,h8,condi,soln, adjust_old,potold,adjustold,udens,gensepold,inequ0, inequ,logoold,treqlist,fl,facold,u,lhslist,cpu,gc, cpustart,gcstart,found,cf0,rtnlist,solns,nontriv, extraline,cf,cfcopy,nx,nde,mindensord,mindensord0, maxdensord; 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 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 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:={}; %------ 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:=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 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; >> >>; 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 derivative in condition maxord:=if eqord>densord then eqord else densord; 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 %--- the conserved current pl and the condition condi condi:=0; deplist:=ulist . for h3:=1:densord collect listdifdif2(lhslist,part(dulist,h3+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:(maxord-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:=append(vl,xlist); % 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 gensym()); h1:=rest h1 >>; >>; inequ:=cons(dequ,inequ) >>; if contrace then write"inequ=",inequ; %--- 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; 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 << write"The function "$ fctprint list reval e1$ write" is not constant!"; extraline:=t; terpri() >>; %--- the currant 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); n:=length h4$ while h3 neq {} do << if length h3 < n then write "+"$ write"( ",first h3," ) * ( ",first h4," )"$ h3:=rest h3; h4:=rest h4 >>$ write" = "$ h4:=xlist$ n:=length h4$ while h2 neq {} do << if length h2 < n then write "+"$ write"df( ",first h2,", ",first h4," )"$ h2:=rest h2; h4:=rest h4 >>; 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>>; return rtnlist end$ % of conlaw3 end$