Artifact 380114db79a5d9eda5cd21dda1f6e8a408ea2baf3e3856ed6c4194e416157943:
- Executable file
r37/packages/crack/conlaw1.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: 19451) [annotate] [blame] [check-ins using] [more...]
comment paralist, found, solns, plcopy, parti_fn ; % CONLAW version 1, to calculate conservation laws of systems % of PDEs by calculating the conserved currents % by Thomas Wolf, September 1997 %---------------------------------------------------------------------- symbolic fluid '(print_ logoprint_ potint_ facint_ adjust_fnc)$ %------------- algebraic procedure conlaw1(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,h9,h10,condi, soln,soln2,adjust_old,potold,adjustold,udens,gensepold, inequ0,inequ,logoold,treqlist,fl,fl2,facold,u,lhslist, cpu,gc,cpustart,gcstart,found,cf0,rtnlist,deglist,maxdf, qlist,extraorder,nx,nde,highdensord,print_old,paralist, plcopy,desyli,ddesyli,vl1,lhsord,loworderlimit,extraline, 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 CONLAW1 - 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 P depend on them and their deriv.s 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:={}; deglist:={}; %------ 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; % the max. degree of any equation extraorder:=0; % increase of order through substitution lhsord:=0; % max. order of an lhs for e1:=1:nde do << e3:=part(eqlist,e1); h2:=0; % degree of this equation h3:=totdeg(part(lhslist,e1),part(ulist,e1)); if h3>lhsord then lhsord:=h3; for each e2 in ulist do << h1:=totdeg(e3,e2); if h1>h2 then h2:=h1; if h1>eqord then eqord:=h1; >>; deglist:=cons(h2,deglist); if h2>h3 then extraorder:=extraorder+h2-h3; >>; deglist:=reverse deglist; %------ determination of the order of the ansatz 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>mindensord then mindensord:=h2 >>; cf0:=t; >> >>; if contrace then write"eqord=",eqord," mindensord=",mindensord, " extraorder=",extraorder," lhsord=",lhsord$; 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: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 << cpu:=lisp time()$ gc:=lisp gctime()$ if cf0 then lisp<<write"A special ansatz of order ",densord," for the", " conserved current is investigated."; terpri()>> else lisp<<write"Currently conservation laws with a conserved", " density";terpri()$ write"of order ",densord," are determined";terpri(); write"======================================================"$ >>; if densord+1<lhsord then lisp << write "The order of the ansatz is too low for substitutions of equations"$ terpri()$write"to be done --> no investigation"$terpri()$terpri() >> else << % If densord+1=lhsord then conservation laws of lower order % than the ansatz are to be considered because this is the % first time they can appear because for a lower order ansatz % substitutions can not be made and no non-trivial CLs can be detected if densord+1=lhsord then loworderlimit:=nil else loworderlimit:=t; if contrace then write"loworderlimit=",loworderlimit$ %--- repeated initializations highdensord:=densord+extraorder; %--- the max. order of P_xi % which is the the order of P_x1 (densord) plus extraorder % to get the max order P_x2,P_x3,... maxord:=highdensord+1+extraorder; %--- the maximal order of % derivatives in condition, = highdensord + extraorder due % to substitutions + 1 due to Div if contrace then write"densord=",densord," highdensord=",highdensord, " 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; 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; vl1:=for e1:=1:(highdensord+2) join part(dulist,e1); % vl1 is to generate subst. of u-deriv. up to order highdensord+1 if contrace then write"vl1=",vl1; %--- initializing the list of unknown functions fl, %--- the conserved current pl and the condition condi if not flist then fl:={} else fl:=flist; deplist:=ulist . for h3:=1:highdensord collect listdifdif2(lhslist,part(dulist,h3+1)); deplist1:=for h3:=1:(densord+1) collect part(deplist,h3); if expl then << deplist :=xlist . deplist; deplist1:=xlist . deplist1>>; deplist :=reverse deplist; deplist1:=reverse deplist1; if contrace then lisp (write"1. depl*=",depl!*); pl:={}; condi:=0; 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}); if n=1 then depnd(h1, deplist1) else depnd(h1, deplist); fl:=h1 . fl; >>;% else h1:=sub(sb,h1); pl:=cons(h1,pl); condi:=condi+totdif(h1,h4,n,dulist) >>; pl:=reverse pl; if contrace then write"fl=",fl," cf=",cf," pl=",pl; if contrace then lisp (write"2. depl*=",depl!*); if contrace then write"condi=",condi; if contrace then write"udens=",udens; %--- generating a substitution list with equations represented % by a symbol and derivatives of equations represented by % derivatives of that symbol %--- at first using the equations themselves sbreserve:={}; desyli:={}; % list of symbols each representing an equation ddesyli:={}; % list of all these symbols + their derivatives % with the following structure: each element of ddesyli has % the form {derivative of the symbol, % {numbers of the differentiation variables}, % number of the symbol} h1:=treqlist; h2:=lhslist; h5:=0; % h5 is an index of the equations while h1 neq {} do << h5:=h5+1; h4:=lisp gensym(); depnd(h4,{xlist}); desyli :=cons(h4, desyli); ddesyli:=cons({h4,{},h5},ddesyli); h3:=first h2; sbreserve:=cons(h3 = h3 - first h1 + h4, sbreserve); h1:=rest h1; h2:=rest h2; >>; sbreserve:=reverse sbreserve; desyli:=reverse desyli; %--- then their derivatives h1:=sbreserve; lisp(h2:=nil); % h2 is list of underived substitutions for each e1 in h1 do lisp << e3:=combidif(algebraic lhs e1); h2:=cons(list(car e3, cdr e3, algebraic rhs e1), h2) % function name and derivative list of e1 >>; for each e1 in vl1 do lisp << % e1 occurs in condi h1:=h2; h5:=0; % counter of the equation while h1 neq nil do << h5:=h5+1; h3:=comparedif2(caar h1, cadar h1, reval algebraic 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 h6:=part(desyli,h5); for each n in h3 do h6:=df(h6,part(xlist,n)); ddesyli:=cons({h6,h3,h5}, ddesyli); sbreserve:=cons(e1 = dequ, sbreserve); lisp(h1:=nil) >> else lisp(h1:=cdr h1); >> >>; if contrace then write"sbreserve=",sbreserve; sb:=sub(for each e1 in desyli collect e1=0,sbreserve); if contrace then write"sb=",sb; rules:={}; for each e1 in sb do << h5:=lhs e1; h6:=rhs e1; rules:=cons(h5 => h6,rules) >>$ if contrace then write"rules=",rules; let rules; condi:=condi; clearrules rules$ if contrace then write"condi=",condi; vl:=append(vl,xlist); % now the full list inequ:=inequ0; %--- inequ is to stop crack if order of pl is too low if loworderlimit then << % for the investigation to stop if % P_x1 is independent of derivatives of order densord dequ:=0; e1:=first pl; 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:=revdulist:=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; plcopy:=sub(second soln,pl); h1:=third soln; % list of free/undeterm. const./functions if contrace then << write"plcopy=",plcopy; 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; 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; %--- h4 is the currant for a single conservation law h4:=for each e2 in plcopy collect for each e1 in h2 sum fdepterms(e2,e1); if contrace then write"h4-1=",h4; sb:=absorbconst(h4,h2)$ if (sb neq nil) and (sb neq 0) then h4:=sub(sb,h4); if contrace then write"h4-2=",h4; if (length(h2)=1) and (fargs first h2 = {}) then << e1:=first h2; h4:=sub(e1=1,h4); >>; if contrace then write"udens=",udens; h5:=udens; if (cf0=nil) and loworderlimit then while (h5 neq {}) and freeof(first h4,first h5) do h5:=rest h5; if contrace then write"h5=",h5; if h5 neq {} then << % P_x1 in h4 is of order densord or no loworderlimit % h3 is the lhs of the conservation law h3:=for e1:=1:nx sum totdif(part(h4,e1),part(xlist,e1),e1,dulist); if contrace then write"h3-1=",h3; if h3 neq 0 then << % non-trivial conservation law %--- Compute the characteristic functions %--- We have already h3 = Div P h3:=sub(sbreserve,h3); if contrace then write"h3-2=",h3; if contrace then write"ddesyli=",ddesyli; divlist:={}; for each e1 in ddesyli do << h6:=coeffn(h3,first e1,1); if h6 neq 0 then << h3:=h3-h6*first e1; divlist:=cons({h6,second e1,third e1},divlist) >> >>; if contrace then write"h3-3=",h3; if contrace then write"divlist=",divlist; qlist:=for e1:=1:nde collect 0; for each e1 in divlist do << % for each derivative of an equ. h9:=first e1; % the coeff of the equn. derivative e2:=second e1; h10:=third e1; if h9 neq 0 then if e2={} then qlist:=part(qlist,h10):= part(qlist,h10)+h9 else << h6:=-1; % the alternating sign if length e2>1 then << h7:=part(treqlist,h10); if paralist neq {} then h7:=sub(second soln,h7); h8:=for each e2 in rest second e1 collect h7:=totdif(h7,part(xlist,e2),e2,dulist)$ >> else h8:={}; h8:=append(h8,{part(treqlist,h10)}); while e2 neq {} do << e3:=first e2; e2:=rest e2; h4:=part(h4,e3):= part(h4,e3)+h6*h9*(first h8); h9:=totdif(h9,part(xlist,e3),e3,dulist)$ if e2 neq {} then << h8:=rest h8; h6:=-h6; >> else qlist:=part(qlist,h10):= part(qlist,h10)+h6*h9 >>; >>; >>; %--- Is the CL trivial, i.e. all Q's are 0 ? e2:=t; for each e1 in qlist do if e1 neq 0 then e2:=nil; if e2 then h4:=nil; if h4 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() >>; cllist:=cons(qlist,cllist); pllist:=cons(h4,pllist); >>; if contrace then write"cllist=",cllist; if contrace then write"pllist=",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,highdensord)$ 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 << h2:=first pllist; h3:=first cllist; %-- Is P_x1 of order densord? To be tested only now after simppl e1:=ulist; if (cf0=nil) and loworderlimit then while (e1 neq {}) and (totdeg(first h2,first e1) < densord) do e1:=rest e1; if e1 neq {} then << found:=t; write"Conservation law:"; 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 if 0 = (for each e1 in h3 sum e1) then write" 0 = " else << 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$ while h2 neq {} do << if length h2 < nx then write "+"$ write"df( ",first h2,", ",first h4," )"$ h2:=rest h2; h4:=rest h4 >>; write"======================================================"$ >>; pllist:=rest pllist; cllist:=rest cllist; >>$ if solns neq {} then nodepnd(ulist); >>; sbreserve:=0; nodepnd(desyli); 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 conlaw_1: ", lisp time() - cpustart, " ms GC time : ", lisp gctime() - gcstart," ms"$ lisp <<adjust_fnc:=adjustold; logoprint_:=logoold; potint_:=potold; facint_:=facold>>; return rtnlist end$ % of conlaw1 end$