% This file contains all the necessary statements for computing the
% symmetries of the KdV equations up to a certain order. By changing
% some of the statements in the first part of the file, it may be easily
% adapted to compute symmetries of other (systems of) equations.

load tools,integrator,supervf;

% Give the set of dependent variables u,v,... etc.
dependent_variables := {u}$

% Since for the odd variable numbers are given instead of variables we introduce
% odd_offset, the number after which the dependent odd variables start (and continue 
% consecutively) and the number of odd variables:
odd_offset:=0$
nr_odd_variables:=0$

% Give the order of the system of pde's 
order_pde := 3$

% Give the order of the symmetries one wants to consider
order_sym := 5$

% Give the expressions for the t derivatives ut,vt,... of the dependent variables
ut:=u3+u*u1$

% Give the set of nonlocal variables to be considered as well.
nonlocal_variables := {}$
nonlocal_odd_offset := 0$
nr_odd_nonlocal := 0$

% Give the x and t derivatives px and pt of all nonlocal variables p.

% Specify the functions f(1),...,f(n), f(-1),...,f(-m) here.
% If not specified, the functions will be made dependent on 
% the proper variables further on.

algebraic operator f,c;

nr_odd_f:=nr_odd_variables$
nr_even_f:=nr_variables$
nr_odd_c:=0$
nr_even_c:=0$

% Give or give not information during construction of equations
% Useful when computing large examples.
write_mke:=nil$

%--------------------------------------------------------------------------------
%  Make no changes behind this line. We will now compute the symmetry conditions 
%  for the vectorfield f(1)*d/du + f(2)*d/dv + ... + D_x(f(1))*d/du_1 + ...,
%  where f(1),f(2),... depend on the variables 
%  x,t, u,v,... ,u1,v1,..., un,vn,..., p1,...,pm, if n is the order of the 
%  symmetry and p1,...,pm are all the nonlocal variables considered;

nr_variables := length dependent_variables$
nr_nonlocal  := length nonlocal_variables$

dim_vars := 2 + nr_nonlocal + nr_variables*(order_pde + order_sym + 1)$
dim_odd_vars := max(odd_offset + nr_odd_variables*(order_pde + order_sym + 1),
                    nonlocal_odd_offset + nr_odd_nonlocal)$

vars := for i:=1:order_pde + order_sym + 1 join 
          for j:=1:nr_variables collect mkid(part(dependent_variables,j),i)$
vars := x . t . append(nonlocal_variables, append(dependent_variables,vars))$

algebraic operator equ,var_x;

initialize_equations(equ,nr_variables+nr_odd_variables,vars,
                     {c,nr_even_c,nr_odd_c},{f,nr_even_f,nr_odd_f});

vectorfield(ddx,vars);
vectorfield(ddt,vars);

% The following procedure gives the number of D_x^n(ui) if ui is the i-th 
% local variable
algebraic procedure var_nr(i,n);
  2+nr_nonlocal+n*nr_variables+i$

algebraic procedure odd_var_nr(i,n);
  odd_offset+n*nr_odd_variables+i$

for i:=1:dim_vars do var_x i:=part(vars,i);

% We construct the components of the total derivatives D_x and D_t
ddx(0,1) := 1$ 
ddx(0,2) := 0$
for i:=1:nr_nonlocal do 
  ddx(0,2+i) := mkid(part(nonlocal_variables,i),x);
for i:=1:nr_variables do 
  for n:=0:order_pde + order_sym do 
    ddx(0,var_nr(i,n)) := var_x(var_nr(i,n+1));
for i:=1:nr_odd_nonlocal do
  ddx(1,nonlocal_odd_offset+i):=mkid(mkid(ext,nonlocal_odd_offset+i),x);
for i:=1:nr_odd_variables do 
  for n:=0:order_pde + order_sym do
    ddx(1,odd_var_nr(i,n)) := ext(odd_var_nr(i,n+1));

procedure mk_ddt;
begin
  ddt(0,1) := 0$
  ddt(0,2) := 1$
  for i:=1:nr_nonlocal do 
    ddt(0,2+i) := mkid(part(nonlocal_variables,i),t);
  for i:=1:nr_variables do 
    ddt(0,var_nr(i,0)) := mkid(part(dependent_variables,i),t);
  for i:=1:nr_variables do
    for n:=1:order_sym do 
      ddt(0,var_nr(i,n)) := ddx ddt(0,var_nr(i,n-1));
  for i:=1:nr_odd_nonlocal do
    ddt(1,nonlocal_odd_offset+i):=mkid(mkid(ext,nonlocal_odd_offset+i),t); 
  for i:=1:nr_odd_variables do 
    ddt(1,odd_var_nr(i,0)) := mkid(mkid(ext,odd_offset+i),t);
  for i:=1:nr_odd_variables do
    for n:=1:order_sym do 
      ddt(1,odd_var_nr(i,n)) := ddx ddt(1,odd_var_nr(i,n-1));
end$

% For the construction of the symmetry condition we need to compute
% the action of the linearisator of the system of pde's on 
% f(1),...,f(n),f(-1),...,f(-m)
% We will save these as equ(1),...,equ(n+m)

vectorfield(symmetry,vars);

procedure make_prolongation;
begin
  for i:=1:nr_variables do
  begin 
    symmetry(0,var_nr(i,0)) := f(i);
    for n:=1:order_pde do <<
      if write_mke then write "Prolongation of f(",i,") up to order ",n;
      symmetry(0,var_nr(i,n)):=ddx symmetry(0,var_nr(i,n-1))>>;
  end;
  for i:=1:nr_odd_variables do
  begin 
    symmetry(1,odd_var_nr(i,0)) := f(-i);
    for n:=1:order_pde do <<
      if write_mke then write "Prolongation of f(",-i,") up to order ",n;
      symmetry(1,odd_var_nr(i,n)):=ddx symmetry(1,odd_var_nr(i,n-1))>>;
  end;
end$

procedure make_equations;
begin
  for i:=1:nr_variables do 
  begin scalar evolution; 
    evolution:=mkid(part(dependent_variables,i),t);
    if write_mke then write "Computing equation for f(",i,")";
    equ(i):=ddt f(i) - symmetry evolution;
  end;
  for i:=1:nr_odd_variables do 
  begin scalar evolution;
    evolution:=mkid(mkid(ext,odd_offset+i),t);
    if write_mke then write "Computing equation for f(",-i,")";
    equ(nr_variables+i):=ddt f(-i) - symmetry evolution;
  end;
  if not write_mke then 
    if (nr_variables+nr_odd_variables)=1 then write "Introduced equation 1"
    else write "Introduced equations ",1,",...,",nr_variables+nr_odd_variables;
end$

% Check if f(-m),...,f(n) are already defined, if not make them
% dependent on the proper variables.

lisp operator has_no_definition;
lisp procedure has_no_definition(opr,i);
if assoc(list(opr,i),get(opr,'kvalue)) then nil else t$

for i:=-nr_odd_variables:nr_variables do
  if i neq 0 and has_no_definition(f,i) then
    <<for k:=1:nr_variables do
        for l:=0:order_sym do depend f(i),var_x(var_nr(k,l));
      depend f(i),x,t
    >>;

% Define some handy abbreviations
define es=integrate_equation,
       seq=integrate_equations,
       xes=integrate_exceptional_equation,
       pr=show_equation,
       preq=show_equations,
       te=equations_used(),
       pte=put_equations_used,
       fu=functions_used,
       pfu=put_functions_used;

% Compute the prolongation of the vectorfield, the components of D_t
% and finally all the equations.

make_prolongation();
mk_ddt();
make_equations();

% For the KdV, cracking the problem is utterly simple:
% (other systems require more skill !!)

auto_solve 1;

end;