reduce-historical/r38/packages/crack/conlaw2.red

%       CONLAW version 2, to calculate conservation laws of systems
%            of PDEs by calculating characteristic functions

%                   by Thomas Wolf, June 1999

%----------------------------------------------------------------------

symbolic fluid '(print_ logoprint_ potint_ facint_ adjust_fnc)$ 

%-------------

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$

%-------------

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,potold,adjustold,udens,gensepold,
  inequ0,inequ,logoold,treqlist,fl,facold,u,nodep,subl,cpu,
  gc,cpustart,gcstart,nontriv,cf0,rtnlist,paralist,solns,
  found,clcopy,extraline,nondiv,nx,nde,nonconstc,
  mindensord,mindensord0,maxdensord,rules$

  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 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 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 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)$

  %--- 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:={};
  nondiv:=lisp intern 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);
    nodep:=sub(e1,nodep);
    subl:=sub(e1,subl);
  >>$
  if contrace then write"treqlist=",treqlist," nodep=",nodep,
                        " subl=", subl;
  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:=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);
    >>;
    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:=subl;
    h4:=nil;       % h4 is list of undifferentiated substitutions
    while h1 neq {} do <<
      h3:=first h2;
      h5:=h3-(first h1)/coeffn(first h1,h3,1);
      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:=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 <<
      % 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 intern gensym());
          h1:=rest h1
        >>;
      >>;
      inequ:=cons(dequ,inequ)
    >>;
    let rules;
    inequ:=inequ;
    clearrules rules$
    rules:=0;
    if contrace then write"inequ=",inequ;

    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:=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,
                    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);

        if not lisp(freeof(h1,'SUB)) then h1:={}
                                     else <<
	  %--- lambda integration of h1 to compute P_i
	  h2:=lisp intern gensym()$
          h10:=ulist;
          while h10 neq {} do
          if not lisp(freeof(h1,'SUB)) then h10:={} 
                                       else <<
            e1:=first h10; h10:=rest h10;
            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);
	    if second h9 = 0 then h9:=first h9
			     else << % no success --> use lambda-integration
	      h9:=nondiv;
	      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$ 
		  %--- the following is to catch errors in: 
		  %--- int(h10*h2**(nx-1),h2)
		  h10:=if not lisp freeof(h10,'SUB) then nil else 
		   err_catch_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:=err_catch_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 neq nondiv then <<
            if nx=2 then
            h2:=first simppl({h2},ulist,first xlist,second xlist)$
            write"Conservation law:";
            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 {} then partintdf(h4,h3,h2,xlist,h6,vl,sb) 
          >>             else <<
            write"Adjoint symmetry:";
 	    while h3 neq {} do <<
	      if length h3 < length first cllist then write "+"$
	      write"(",first h3,") * (",first h4,")"$
	      h3:=rest h3; h4:=rest h4
	    >>$

            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 found characteristic function represents ",
                   "a conservation law"$
            >>$
	  >>;
	  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>>;
  recover_reduce_flags()$

  return rtnlist

end$ % of conlaw2

end$