%********************************************************************
module transform$
%********************************************************************
% Routines for performing transformations
% Author: Thomas Wolf
% March 1999
%
symbolic procedure input_trafo$
begin scalar ulist,vlist,u,v,ylist,xlist,yslist,xslist,oldprompt,xl2,
notallowed,full_simplify$
oldprompt:=promptstring!*$
promptstring!*:=""$
write"Under the following conditions this program performs arbitrary"$
terpri()$
write"transformations."$
terpri()$terpri()$
write"If not only variables but also functions are transformed then it"$
terpri()$
write"is assumed that all new functions depend on the same new variables"$
terpri()$
write"and that all old functions depend on the same old variables."$
terpri()$ terpri()$
write"For these procedures to be applicable the old functions and variables"$
terpri()$
write"must be given explicitly in terms of the new ones, not involving"$
terpri()$
write"unspecified functions of the new ones. Also the differential "$
terpri()$
write"equations to be transformed must contain the old independent and"$
terpri()$
write"dependent variables and their partial derivatives explicitly."$
terpri()$
% write"Hint: Splitting a single transformation involving many"$
% terpri()$
% write"variables into many transformations involving each only few"$
% terpri()$
% write"variables speeds the whole transformation up."$
% terpri()$
terpri()$
write"Give a list of new functions, e.g. `u1,u2,u3;' in the order to"$
terpri()$
write"be used to sort dervatives. If there are no new functions enter ;"$
terpri()$
ulist := termlistread()$
terpri()$
write"Give a list of all new variables, e.g. 'v1,v2,v3;' in the order to"$
terpri()$
write"be used to sort derivatives. If there are no new variables enter ;"$
terpri()$
vlist := termlistread()$
if ulist then <<
for each u in ulist do
for each v in vlist do depend u,v$
terpri()$
write"Give a list of all substitutions of old functions in terms of"$
terpri()$
write"new functions and new variables, e.g. y1=2*u1+u2*v2, y2=u2-u1*v1;"$
terpri()$
write"If there are no substitutions of old functions enter ;"$
terpri()$
yslist := termlistread()$
% Check whether all old functions do depend on the same variables
ylist:=for each u in yslist collect cadr u$
xlist:=fctargs car ylist$
for each u in cdr ylist do <<
xl2:=fctargs u$
if not_included(xlist,xl2) or
not_included(xl2,xlist) then <<
notallowed:=t$
terpri()$
write"Functions ",car ylist,",",u," do not depend on the same variables!"$
>>
>>
>>$
if notallowed then return nil;
if vlist then <<
terpri()$
write"Give a list of all substitutions of old variables in terms of"$
terpri()$
write"new functions and new variables, e.g. x1=v1-v2*u2, x2=3*v2+v1*u1;"$
terpri()$
xslist := termlistread()$
>>$
terpri()$
write"Shall the transformed equation be fully simplified,"$
terpri()$
write"i.e. redundand non-vanishing factors be dropped y/n : "$
full_simplify:=termread()$
if (full_simplify='n) then full_simplify:=nil;
terpri()$
% Dependence of the new dependent variables on old independent variables
% which are not transformed
for each v in xslist do xlist:=setdiff(xlist,list cadr v);
for each u in ulist do
for each v in xlist do depend u,v$
% Also non-changing old variables must enter the transformation as
% partial derivatives wrt them will have a different meaning
vlist:=append(vlist,xlist)$
for each v in xlist do xslist:=cons({'EQUAL,v,v},xslist)$
% If a test is necessary that all old variables are replaced then do
% if (not not_included(ftem_,newli)) and
% (not not_included(newli,ftem_)) then
promptstring!*:=oldprompt$
if print_ then <<
write"The transformation:"$terpri()$
if vlist then <<
write"The new variables: "$
listprint(vlist)$terpri()
>>;
if ulist then <<
write"The new functions: "$
listprint(ulist)$terpri()
>>;
if xslist then <<
write"The old variables expressed:"$terpri()$
mathprint cons('LIST,xslist)
>>;
if yslist then <<
write"The old functions expressed:"$terpri()$
mathprint cons('LIST,yslist)
>>;
>>;
return {'LIST,cons('LIST,ulist),
cons('LIST,vlist),
cons('LIST,yslist),
cons('LIST,xslist),
full_simplify }
end$
%----------------------------
symbolic procedure adddep(xlist)$
% all functions depending on old variables get a dependency on
% the new variables
% xlist is a lisp list ((x1,v1,v4,v5),(x2,v2,v3,v4),...)
begin scalar newdep,xs,dp;
for each xs in xlist do <<
newdep:=nil$
while depl!* do <<
dp:=car depl!*;
depl!*:=cdr depl!*;
if not freeof(dp,car xs) then
dp:=cons(car dp,union(cdr xs,cdr dp))$
newdep:=cons(dp,newdep);
>>;
depl!*:=reverse newdep;
>>;
end$
%----------------------------
symbolic procedure dropdep(xlist)$
% xlist is a lisp list
begin scalar x,dp,newdep$
for each x in xlist do <<
newdep:=nil$
while depl!* do <<
dp:=car depl!*;
depl!*:=cdr depl!*;
if not freeof(dp,x) then
dp:=delete(x,dp)$
newdep:=cons(dp,newdep);
>>;
depl!*:=reverse newdep
>>;
end$
%----------------------------
%symbolic operator TransfoDf$
symbolic procedure TransfoDf(dy,yslist,xlist,vlist)$
% - dy is the derivative to be transformed
% - yslist is a list of already computed substitutions for the old
% functions and their derivatives
% - xlist is a list of the old variables
% - vlist is a list of the new variables
% All parameters are in prefix form.
% yslist,xlist,vlist are lisp lists
% returns cons(substitution for dy, complete list of substitutions)
begin
scalar cpy,x,dym1,m,n,newdy,v$
cpy:=yslist$
while cpy and
(dy neq cadar cpy) do cpy:=cdr cpy;
return
if not null cpy then cons(car cpy,yslist) else % found rule
if not pairp dy then cons({'EQUAL,dy,dy},yslist) else << % no dy-rule
% dym1 is one lower x-derivative than dy
if ( length dy = 3 ) or
((length dy = 4) and
(cadddr dy = 1) ) then <<x:=caddr dy;dym1:=cadr dy>>
else <<
cpy:=reverse dy;
dym1:=reverse
if not numberp car cpy then <<x:= car cpy; cdr cpy >> else
if (car cpy = 1) then <<x:=cadr cpy; cddr cpy >> else
if (car cpy = 2) then <<x:=cadr cpy; cdr cpy >> else
<<x:=cadr cpy; cons(sub1 car cpy,
cdr cpy)>>
>>;
yslist:=TransfoDf(dym1,yslist,xlist,vlist);
dym1:=car yslist; % dym1 is now a substitution rule for dym1 above
dym1:=caddr dym1; % dym1 is now the expression to be substituted
yslist:=cdr yslist; % the new substitution list
% computation of the subst. rule for dy
m:=1;
while xlist and (x neq car xlist) do <<m:=add1 m; xlist:=cdr xlist >>$
if null xlist then newdy:=reval {'DF,dym1,x}
else <<
n:=0;
for each v in vlist do <<
n:=add1 n;
if not zerop algebraic(Dv!/Dx(n,m)) then
newdy:=cons({'TIMES,{'DF,dym1,v},algebraic(Dv!/Dx(n,m))},
newdy)$
% {'DF,dym1,v} is the full total derivative as it should be
% provided all functions depend directly on v (as stored in depl!*)
% or they do not depend on v but not like f(u(v)) with an
% unspecified f(u)
>>;
newdy:=if cdr newdy then reval cons('PLUS,newdy)
else if newdy then reval car newdy
else 0
>>$
% return the new subst. rule and the new yslist
cons({'EQUAL,dy,newdy},cons({'EQUAL,dy,newdy},yslist))
>>
end$ % of TransfoDf
%----------------------------
symbolic procedure Do_Trafo(arglist,x)$
begin
scalar yslist,xslist,ulist,vlist,xlist,ylist,m,n,ovar,nvar,e1,e2,e3,
x,detpd,pdes,hval,trfo,newforg,newineq_,drvs,full_simplify$
%dyx!/duv,Dv!/Dx
algebraic <<
% input of the transformation
ulist :=first x$ x:=rest x$
vlist :=first x$ x:=rest x$
yslist:=first x$
xslist:=second x$
full_simplify:=third x$ x:=nil$
>>$
% update of depl!*
xlist:=
for each e1 in cdr xslist collect <<
x:=caddr e1$
e3:=nil;
for each e2 in cdr vlist do
if not freeof(x,e2) then e3:=cons(e2,e3);
cons(cadr e1,e3)
>>$
adddep(xlist)$
algebraic <<
% checking regularity of the transformation
m:=length(xslist); n:=length(yslist)+m;
clear dyx!/duv,Dv!/Dx;
matrix dyx!/duv(n,n);
matrix Dv!/Dx(m,m);
ovar:=append(yslist,xslist);
nvar:=append(ulist,vlist);
n:=0;
for each e1 in ovar do <<
n:=n+1;m:=0;
for each e2 in nvar do <<
m:=m+1;
dyx!/duv(m,n):=df(rhs e1,e2)
>>
>>;
detpd:=det(dyx!/duv);
if detpd=0 then <<write"The proposed transformation is not regular!"$
return nil>>;
clear dyx!/duv;
% computation of Dv/Dx:=(Dx/Dv)^(-1)
n:=0;
for each e1 in xslist do <<
n:=n+1;m:=0;
for each e2 in vlist do <<
m:=m+1;
Dv!/Dx(n,m):=total_alg_mode_deriv(rhs e1,e2)
% It is assumed that ulist does depend on vlist
>>
>>;
Dv!/Dx:=Dv!/Dx**(-1);
>>$
xslist:=cdr xslist$
yslist:=cdr yslist$
vlist :=cdr vlist$
ulist :=cdr ulist$
% update of global data ftem_, vl_
if ulist then <<
for each e1 in yslist do ftem_:=delete(cadr e1,ftem_);
for each e1 in ulist do ftem_:=fctinsert(e1,ftem_)$
>>$
xlist:=for each e1 in xslist collect cadr e1$
for each e1 in xlist do vl_:=delete(e1,vl_);
vl_:=append(vl_,vlist)$
ylist:=for each e1 in yslist collect cadr e1$
% update of the pdes
pdes:=car arglist$
for each e1 in pdes do <<
hval:=get(e1,'val)$
drvs:=append(search_li2(hval,'DF),ylist)$
for each e3 in drvs do <<
trfo:=TransfoDf(e3,yslist,xlist,vlist)$
hval:=subst(caddar trfo,cadar trfo,hval);
yslist:=cdr trfo
>>$
for each e2 in xslist do
if not freeof(hval,cadr e2) then
hval:=subst(caddr e2,cadr e2,hval);
put(e1,'val,hval);
>>$
% update of forg
for each e1 in cadr arglist do
if (pairp e1) and (car e1 = 'EQUAL) then <<
hval:=caddr e1;
drvs:=append(search_li2(hval,'DF),ylist)$
for each e3 in drvs do <<
trfo:=TransfoDf(e3,yslist,xlist,vlist)$
hval:=subst(caddar trfo,cadar trfo,hval);
yslist:=cdr trfo
>>$
for each e2 in xslist do
if not freeof(hval,cadr e2) then
hval:=subst(caddr e2,cadr e2,hval);
hval:=reval hval;
newforg:=cons({'EQUAL,cadr e1,hval},newforg)$
e2:=nil;
for each e3 in ftem_ do
if not freeof(hval,e3) then e2:=cons(e3,e2);
put(cadr e1,'fcts,e2)
>> else
if not freeof(ylist,e1) then <<
e3:=yslist;
while e3 and cadar e3 neq e1 do e3:=cdr e3$
if e3 then newforg:=cons(car e3,newforg)
else newforg:=cons(e1,newforg)
>> else
newforg:=cons(e1,newforg);
% update of ineq_
newineq_:=nil;
for each e1 in ineq_ do <<
drvs:=append(search_li2(e1,'DF),ylist)$
for each e3 in drvs do <<
trfo:=TransfoDf(e3,yslist,xlist,vlist)$
e1:=subst(caddar trfo,cadar trfo,e1);
yslist:=cdr trfo
>>$
for each e2 in xslist do
if not freeof(e1,cadr e2) then
e1:=subst(caddr e2,cadr e2,e1);
newineq_:=cons(reval e1,newineq_)
>>$
ineq_:=nil;
for each e1 in newineq_ do addineq(pdes,e1);
xlist:=nil;
for each e1 in xslist do
if cadr e1 neq caddr e1 then xlist:=cons(cadr e1,xlist);
dropdep(xlist)$
for each e1 in pdes do <<
for each e2 in allflags_ do flag1(e1,e2)$
update(e1,get(e1,'val),ftem_,vl_,full_simplify,list(0),pdes)$
drop_pde_from_idties(e1,pdes,nil);
drop_pde_from_properties(e1,pdes)
>>$
% cleanup
algebraic clear Dv!/Dx;
return {pdes,newforg,vl_}
end$ % of Do_Trafo
%----------------------------
symbolic procedure Find_Trafo(arglist)$
begin
scalar dli,avf,f,ps,sol,pde,pdes,forg,batch_bak,print_bak,vlist,
xslist,vl,h1,h2,h3,h4,trtr,eligfncs,eligpdes,epdes,remain,
maxvno;
% trtr:=t$
ps:=promptstring!*$
promptstring!*:=""$
pdes:=car arglist$
% If there are functions of fewer variables then transformations can
% make them to functions of more variables which can add solutions.
% One could first compute the transformation and then check whether
% there is an ftem_ function which has an enlarged set of dependent
% variables and in this case either drops the transformation or one
% adds extra conditions df(f,y)=0 (where d/dy is to be transformed)
% for these functions. Instead a preliminary simpler routs is taken
% in the following, ftem_ may contain only constants or functions
% of the same number of variables.
maxvno:=0;
h1:=ftem_;
while h1 and <<
h3:=fctlength car h1$
if h3=0 then t else
if maxvno=0 then <<maxvno:=h3;t>> else
if h3=maxvno then t else nil
>> do h1:=cdr h1;
if h1 then return <<
if print_ then <<
write"Non-constant functions of fewer variables prevent"$terpri()$
write"the application of this technique."$terpri()
>>$
nil
>>$
if trtr then <<write"111"$terpri()>>$
% Find eligible PDEs
while pdes do <<
pde:=car pdes;pdes:=cdr pdes;
if get(pde,'nvars)=maxvno then <<
eligfncs:=nil;
avf:=get(pde,'allvarfcts)$
if avf and null cdr avf then <<
% There must only be one function of all variables because
% the other one would be part of the inhomogeneity and
% derivatives of this function would give errors in quasilinpde
% when the differentiation variable becomes a function in the
% characteristic ODE system and substitutions are done where
% the function is substituted by an expression that has been
% computed. But also if no derivatives occur, crack is strictly
% speaking not able to deal with funtions of functions.
% Therefore only one function apart from constants is allowed.
f:=car avf;
dli:=get(pde,'derivs);
h1:=t; h2:=0; % h2 counts the first order derivatives of f
while dli and h1 do
if (not pairp caar dli) or
(caaar dli neq f) then dli:=cdr dli else
if null cdaar dli then dli:=cdr dli else % f algebraic
if null cddaar dli then <<h2:=add1 h2;dli:=cdr dli>>
else h1:=nil;
if null dli and (h2 > 1) then eligfncs:=cons(f,eligfncs)
>>$
if eligfncs then <<
eligpdes:=cons(cons(pde,eligfncs),eligpdes);
epdes:=cons(pde,epdes)
>>
>>
>>$
if trtr then <<write"222"$terpri()>>$
if null epdes then return nil;
if expert_mode then pde:=selectpdes(epdes,1)
else <<
if trtr then <<write"333"$terpri()>>$
% Find PDEs with min number of allvar functions
h2:=10000;
for each h1 in epdes do <<
h3:=length get(h1,'allvarfcts);
if h3<h2 then <<h2:=h3;remain:={h1}>> else
if h3=h2 then remain:=cons(h1,remain);
>>;
epdes:=remain;
if trtr then <<write"444"$terpri()>>$
% Find PDEs with max number of variables
h2:=0;
for each h1 in epdes do <<
h3:=get(h1,'nvars);
if h3>h2 then <<h2:=h3;remain:={h1}>> else
if h3=h2 then remain:=cons(h1,remain);
>>;
epdes:=remain;
if trtr then <<write"555"$terpri()>>$
% Find shortest of these PDEs
h2:=10000;
for each h1 in epdes do <<
h3:=get(h1,'terms);
if h3<h2 then <<h2:=h3;remain:={h1}>> else
if h3=h2 then remain:=cons(h1,remain);
>>;
epdes:=remain;
if trtr then <<write"666"$terpri()>>$
pde:=car epdes$ % One could select further the one with the
% fewest variables involved in the transformation
while eligpdes and caar eligpdes neq pde do eligpdes:=cdr eligpdes;
f:=cadar eligpdes;
>>$
if trtr then <<write"777"$terpri()>>$
if null pde then return nil;
if trtr then <<write"888"$terpri()>>$
if print_ then <<
write"Finding a transformation to integrate the 1st order PDE ",pde,":"$
terpri()$
>>$
print_bak:=print_; print_:=nil$
batch_bak:=!*batch_mode; !*batch_mode:=t$
pdes:=car arglist$
forg:=cadr arglist$
h1:=level_string(session_)$
h1:=bldmsg("%s%s.",h1,"qlp")$
backup_to_file(pdes,forg,h1)$ % moved before again:, should be ok
if trtr then <<write"999"$terpri()>>$
sol:=reval algebraic(quasilinpde(lisp(get(pde,'val)),f,
lisp(cons('LIST,get(pde,'vars)))))$
restore_backup_from_file(pdes,forg,h1)$
delete!-file h1;
if trtr then <<write"000"$terpri()>>$
if trtr then <<write"sol0="$mathprint sol$terpri()>>$
!*batch_mode:=batch_bak$ print_:=print_bak$
if null sol or null cdr sol or null cdadr sol then return nil$
sol:=cadr sol;
h1:=cdr sol;
for each h2 in h1 do
if not freeof(h2,f) then sol:=delete(h2,sol);
% One could use lin_check(h2,{f}) to test linearity if needed
h1:=cdr sol;
if trtr then <<write"f=",f$terpri()$
write"h1=",h1$terpri()$
write"sol0="$mathprint sol
>>$
% make a list of all variables occuring in these expressions
% these are all the variables to occur in the transformation
h2:=get(pde,'vars)$
for each f in h2 do if member(f,sol) then h2:=delete(f,h2);
% Keep only the algebraic expressions, drop the single variables
for each f in h1 do if atom f then sol:=delete(f,sol);
if trtr then <<write"sol1="$mathprint sol>>$
% find the variable for which the algebraic expressions are
% most easily solved
if trtr then <<write"h2=",h2$terpri()>>$
if trtr then <<write"xslist=",xslist$terpri()>>$
xslist:=err_catch_solve(sol,cons('LIST,h2))$
if null xslist then return <<
write"REDUCE was not able to solve"$mathprint sol$
write"for one of "$listprint(h2);
nil
>> else xslist:=cdr reval car xslist$
if trtr then <<write"xslist=",xslist$terpri()>>$
h3:=nil;
while xslist do <<
f:=car xslist; xslist:=cdr xslist;
if (car f='EQUAL) and
((pairp caddr f) and
(caaddr f = 'ARBCOMPLEX)) then <<
h2:=delete(cadr f,h2);
h3:=cons(cadr f,h3);
xslist:=subst(1,caddr f,xslist)
>>
>>$
if trtr then <<write"h3new=",h3$terpri()$
write"h2new=",h2$terpri()>>$
for each f in get(pde,'vars) do
if not freeof(sol,f) and
not member(f,h3) and
not member(f,h2) then h3:=cons(f,h3);
if trtr then <<write"sol2=",sol$terpri()>>$
sol:=append(sol, h3);
if trtr then <<write"sol3=",sol$terpri()>>$
%terpri()$
%write"Give a list of ",length cdr sol," new variable names, like: u,v,w;"$
%terpri()$
%write"one for each of the expressions in the following list:"$
%mathprint sol;
%promptstring!*:="list of new variable names: "$
%vlist:=termlistread()$
%promptstring!*:=ps$
h3:=append(h2,h3);
h4:=h3;
vlist:=for each f in h3 collect mkid(f,'!%);
if trtr then <<write"vlist=",vlist$terpri()>>$
h1:=vlist$
if trtr then <<write"sol-2="$mathprint sol$terpri()>>$
sol:=cdr sol;
h2:=sol;
while sol do <<
vl:=cons({'DIFFERENCE,car sol,car h1},vl);
sol:=cdr sol;
h1:=cdr h1;
h4:=cdr h4;
>>$
while h4 do <<
vl:=cons({'DIFFERENCE,car h4,car h1},vl);
h1:=cdr h1;
h4:=cdr h4
>>$
if trtr then <<write"2nd SOLVE: vl="$mathprint cons('LIST,vl)$terpri()$
write"h3=",h3>>$
xslist:=err_catch_solve(cons('LIST,vl),cons('LIST,h3))$
if null xslist then return <<
write"REDUCE was not able to solve"$mathprint cons('LIST,vl)$
write"for the variables "$listprint(h3);
nil
>> else xslist:=car xslist$
if trtr or print_ then <<
write"The following variable transformation expresses variables"$
terpri()$
listprint(h3);
write" through variables "$
listprint(vlist); write" :"$terpri()$
for each f in cdr xslist do mathprint f
>>$
h3:=for each h1 in vl collect {'EQUAL,caddr h1,cadr h1};
done_trafo:=cons('LIST,cons(cons('LIST,h3),cdr done_trafo));
return Do_Trafo(arglist,{'LIST,{'LIST},cons('LIST,vlist),
{'LIST},xslist,t %full_simplify
});
end$ % of Find_Trafo
%----------------------------
symbolic procedure General_Trafo(arglist)$
% Doing a transformation for all data relevant in CRACK
% Tramsformation rule for partial derivatives d using total
% derivatives D:
%
% / p \ -1
% d | Dx | D
% --- = | --- | * ---
% p | i | i
% dx \ Dv / Dv
%
begin
scalar x;
x:=input_trafo()$
if null x then return
<<terpri()$write"No proper input --> no transformation"$nil>>$
return Do_Trafo(arglist,x)
end$
%----------------------------
endmodule$
end$