% CONLAW version 1, to calculate conservation laws of systems
% of PDEs by calculating the conserved currents
% by Thomas Wolf, June 1999
%----------------------------------------------------------------------
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,nodep,subl,
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,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:= 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;
lisp<<terpri()$write "for the function(s): "$
fctprint cdr reval ulist;terpri()>>$
write"======================================================"$
%--- nodep is a list of derivatives the P do not depend on
%--- subl is the list of lhs-derivatives to be substituted
h1:=lhsli(eqlist)$
nodep:=first h1$
subl:=second h1$
%--- Here comes a test that lhs's are properly chosen
chksub(eqlist,ulist)$
eqlist:=reverse for each e1 in eqlist collect
if part(e1,0)=EQUAL then num(lhs e1 - rhs e1)
else num e1;
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(subl,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);
nodep:=sub(e1,nodep);
subl:=sub(e1,subl);
>>$
if contrace then write"treqlist=",treqlist," nodep=",nodep,
" subl=", subl;
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:=lisp(cons('LIST,setdiff(cdr ulist,cdr nodep))) .
for n:=1:highdensord collect
listdifdif2(nodep,part(dulist,n+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:=subl;
h5:=0; % h5 is an index of the equations
while h1 neq {} do <<
h5:=h5+1;
h4:=lisp intern gensym();
depnd(h4,{xlist});
desyli :=cons(h4, desyli);
ddesyli:=cons({h4,{},h5},ddesyli);
h3:=first h2;
if h3 neq 0 then
sbreserve:=cons(h3 = h3 - (first h1)/coeffn(first h1,h3,1) + h4,
sbreserve);
h1:=rest h1;
h2:=rest h2;
>>;
sbreserve:=reverse sbreserve;
desyli:=reverse desyli;
if contrace then write"rev desyli=",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:=reverse append(xlist,vl); % 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 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:=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;
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;
%--- 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 <<
h6:=h4;
while (h6 neq {}) and
freeof(first h6,first h5) do h6:=rest h6;
if h6={} then t else nil
>> 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 <<
nonconstc:=cons('LIST,cons(reval e1,cdr nonconstc));
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,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 <<
h2:=first pllist;
h3:=first cllist;
%-- Is any P_x of order densord? To be tested only now after simppl
e1:=ulist;
if (cf0=nil) and loworderlimit then
while (e1 neq {}) and <<
h4:=h2;
while (h4 neq {}) and
(totdeg(first h4,first e1) < densord) do h4:=rest h4;
if h4={} then t else nil
>> 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 "$
h4:=xlist$
while h2 neq {} do <<
if length h2 < nx then write "+"$
write"df( ",first h2,", ",first h4," )"$
h2:=rest h2;
h4:=rest h4
>>
>> else <<
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 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 conlaw1: ", 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 conlaw1
end$