module applysym;
load_package crack, crackapp;
create!-package('(applysym),'(contrib crack));
fluid '(print_ fname_ nfct_ time_ facint_ logoprint_);
lisp flag('(yesp),'boolean)$
lisp(logoprint_:=t)$
symbolic operator freeoflist$
%----------------------------
algebraic procedure applysym(problem,symtry);
% problem ... {{equations},{functions},{variables}}
% Symtry ... {{xi_..=..,..,eta_..=..,..},{constants in first list}}
begin
scalar genlist,con,e1,e2,h1,h2,h3,h4,h5,h6,h7,modus,u,v,xlist,ylist,n,
cop1,cop2,symanz,oldsol,oldmodus,trafoprob,altlogo$
clear sy_,sym_;
array sym_(length(second symtry));
symbolic put('sy_,'simpfn,'simpiden)$
symbolic put('ff,'simpfn,'simpiden)$
symbolic put('ffi,'simpfn,'simpiden)$
ylist := maklist second problem;
xlist := maklist third problem;
con:=second symtry;
symanz:=0;
for each e1 in con do % i.e. for all symmetries do:
if freeoflist(e1,xlist) and
freeoflist(e1,ylist) then %------ no pseudo-Lie-symm.
<<genlist:=sub(e1=1,first symtry); %------ calculate symmetry
for each el2 in con do
if el1 neq el2 then genlist:=sub(el2=0,genlist);
symanz:=symanz+1;
sym_(symanz):=genlist
>>;
repeat <<
% Application
oldmodus:=modus;
repeat <<
lisp <<terpri()$
write"Do you want to find similarity and symmetry variables ",
"(enter `1;')";terpri()$
write"or generalize a special solution with new parameters ",
"(enter `2;')";terpri()$
write"or exit the program ",
"(enter `;')";terpri()
>>;
modus:=aread(nil)
>> until (modus=1) or (modus=2) or (modus=nil);
if modus neq nil then <<
% Preparing a combination of symmetries
if symanz=1 then genlist:=sym_(1)
else <<
for n:=1:symanz do <<
write"---------------------- The ",n,". symmetry is:"$
for each e2 in sym_(n) do
if rhs e2 neq 0 then write e2
>>;
write"----------------------"$
repeat
<<lisp
<<terpri()$
write
"Which single symmetry or linear combination of symmetries"$
terpri()$write"do you want to apply? "$
terpri()$
write"Enter an expression with `sy_(i)' for the i'th ",
"symmetry. Terminate input with `$' or `;'."$
terpri()
>>$
h1:=aread(nil);
if h1 then <<
for each h2 in xlist do
if h1 then if df(h1,h2) neq 0 then h1:=nil;
for each h2 in ylist do
if h1 then if df(h1,h2) neq 0 then h1:=nil;
if h1=nil then lisp <<
terpri();
write"The coefficients of the sy_(i) must be constant, ",
"i.e. numbers or constants";terpri()
>>
>>
>> until h1;
genlist:={};
cop1:=sym_(1);
while cop1 neq {} do <<
h6:=lhs first cop1;cop1:=rest cop1;
genlist:=cons(h6 = 0, genlist)
>>;
genlist:=reverse genlist;
for h2:=1:symanz do <<
h3:=coeffn(h1,sy_(h2),1);
if h3 neq 0 then <<
cop1:=genlist;cop2:=sym_(h2);
genlist:={};
while cop1 neq {} do <<
h4:=first cop1; cop1:=rest cop1$
h5:=first cop2; cop2:=rest cop2$
h6:=lhs h4;
genlist:=cons(h6 = rhs h4 + h3*(rhs h5),genlist)
>>;
genlist:=reverse genlist
>>
>>
>>;
write"The symmetry to be applied in the following is ";
write genlist;
if modus=1 then <<
write"Enter the name of the new dependent variable(s):";
u:=aread(nil);
write"Enter the name of the new independent variable(s):";
v:=aread(nil);
altlogo:=logoprint_;
logoprint_:=nil;
trafoprob:=similarity(problem,genlist,{},u,v);
logoprint_:=altlogo
>> else <<
if length h1 < 2 then <<
lisp <<terpri()$
write"What shall the name of the new constant parameter",
" be? ";terpri()>>$
h2:=aread(nil)
>>;
repeat <<
lisp <<terpri()$
write"Enter the solution to be generalized in form of an ",
"expression, which vanishes";terpri()$
write"or in form of an equation `... = ...' ";
if oldsol=nil then write":" else <<
terpri()$
write"or enter semicolon `;' to work on the solution ",
"specified before:"
>>;terpri()
>>;
h3:=aread(nil);
if h3 neq nil then oldsol:=h3
>> until oldsol neq nil;
h3:=newparam(oldsol,genlist,h2);
if h3 neq nil then oldsol:=h3
>>
>>
>> until modus=nil;
clear sym_;
return if oldmodus=1 then trafoprob else
if oldmodus=2 then oldsol else nil
end$ % of ApplySym
%----------------------------
algebraic procedure storedepend(xlist)$
% stores the dependencies of the elements of xlist in the returned
% list and clears the dependencies
begin scalar q,e1,e2$
return
for each e1 in xlist collect
<< q:=fargs e1;
for each e2 in q do nodepend e1,e2;
cons(e1,q)>>
end$ % of storedepend
%----------------------------
algebraic procedure restoredepend(prev_depend)$
% assigns the dependencies stored in the argument
begin scalar q,s,x;
for each s in prev_depend do
<<q:=first s; s:=rest s;
for each x in s do depend q,x>>
end$ % of restoredepend
%----------------------------
algebraic procedure newparam(oldsol,genlist,u_)$
% u_ is the name of the new parameter
begin
scalar h1,h2,h20,h3,h30,h4,vari,pde,e1,printold,clist,newsol,
oldsol_ex,prev_depend;
vari:=makepde(genlist,u_)$
pde:=first vari;
vari:=rest vari;
% is oldsol an invariant?
oldsol_ex:=equ_to_expr(oldsol)$ % oldsol as vanishing expression
prev_depend:=storedepend(vari)$
h2:=sub(u_=oldsol_ex,pde);
if h2 neq 0 then <<
h1:=solve(oldsol_ex,vari);
if h1 neq {} then <<
h1=first h1;
for each h3 in h1 do
if freeof(h3,arbcomplex) then h2:=sub(h3,h2)
>>
>>;
restoredepend(prev_depend)$
if 0=h2 then return lisp
<<write"The special solution to be generalized is an invariant ",
"with respect to";terpri()$
write"this symmetry, therefore no generalization is possible.";
terpri()$
for each h1 in fargs u_ do nodepend u_,h1;
algebraic oldsol
>>;
pde:=pde-1;
h1:= quasilinpde1(pde,u_,vari)$
if h1 neq {} then <<
h1:=first h1;
% h2 is expressing the constants in terms of u_ and the xlist,ylist
clist:={};
lisp<<printold:=print_; print_:=nil>>;
h2:= for each e1 in h1 collect
<<h3:=lisp(newfct(fname_,nil,nfct_))$
lisp(nfct_:=add1 nfct_)$
clist:=cons(h3,clist);
h3 = e1>>;
lisp(print_:=printold);
h20:=sub(u_=0,h2);
h3:=solve(h2,vari);
if h3 neq {} then <<
h3:=first h3;
h30:=sub(h20,h3);
write"The substitutions to generalize the solution are: "$
for each h4 in h30 do write h4$
newsol:=sub(h30,oldsol_ex);
% newsol:=second dropredundant(newsol,
% alle Konstanten und
% Funktionen die nicht in vari sind,
% vari)$
lisp <<write"The new solution";
% if length algebraic h >2 then write"s are:" else write" is:";
% terpri()$
>>;
write"0 = ",newsol
>>
>>;
for each h1 in fargs u_ do nodepend u_,h1;
return newsol
end$ % of NewParam
%----------------------------
%symbolic operator aread$
%symbolic procedure aread(prompt)$
%begin scalar val, !*echo; % Don't re-echo tty input
% rds nil; wrs nil$ % Switch I/O to terminal
% terpri()$
% if null prompt then prompt:= "Input:"; prin2 prompt;
% val := xread(nil)$
% if ifl!* then rds cadr ifl!*$ % Resets I/O streams
% if ofl!* then wrs cdr ofl!*$
% return val
%end$
%----------------------------
symbolic operator einfachst$
symbolic procedure einfachst(a,x)$
% a is an algebraic list and the element where x appears, but
% appears simplest, is found
begin
scalar el1,el2,hp,printold;
printold:=print_;
print_:=nil;
hp:=10000;
a:=cdr a;
while a do <<
el2:=car a;a:=cdr a;
if not freeof(el2,x) and
((not el1 ) or
(el2 = x ) or
<<if not polyp(el2,cons(x,nil)) then nil
else
<<coeff1(el2,x,nil)$
if hipow!*<hp then <<hp:=hipow!*; t>>
else nil >> >> ) then el1:=el2
>>;
print_:=printold;
return el1
end$
%----------------------------
algebraic procedure transdf(y,yslist,vlist,indxlist)$
begin
scalar m,n,e1,dfy;
return
if indxlist={} then sub(yslist,y)
else <<
m:=first indxlist;
n:=0;
dfy:=transdf(y,yslist,vlist,rest indxlist);
for each e1 in vlist sum <<
n:=n+1;
df(dfy,e1)*dv!/dx(n,m)
>>
>>
end$ % of TransDf
%----------------------------
algebraic procedure transderiv(yik,yslist,vlist)$
begin
scalar indxlist,y,l1,l2;
indxlist:=lisp cons('list,combidif(yik))$
return transdf(first indxlist,yslist,vlist,rest indxlist)
end$ % of TransDeriv
%----------------------------
algebraic procedure detrafo(eqlist,yslist,xslist,ulist,vlist)$
% Transformations of all orders are performed (point-, contact-,...)
% but only x-derivatives of y's are transformed. To transform other
% any other derivatives, subdiff1 must be extended to include all other
% occuring x-derivatives.
begin
scalar avar,nvar,detpd,n,m,ordr,e1,e2,e3,sb;
m:=length(xslist); n:=length(yslist)+m;
clear dyx!/duv,dv!/dx;
matrix dyx!/duv(n,n);
matrix dv!/dx(m,m);
avar:=append(yslist,xslist);
nvar:=append(ulist,vlist);
n:=0;
for each e1 in avar 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);
%write"detpd=",detpd;
if detpd=0 then return
<<write"The proposed transformation is not regular!";{}>>;
clear dyx!/duv;
ordr:=0;
for each e1 in eqlist do
for each e2 in yslist do
<<n:=totdeg(e1,lhs e2);
if n>ordr then ordr:=n>>;
sb:=subdif1(for each e1 in xslist collect lhs e1,
for each e1 in yslist collect lhs e1,ordr);
% 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):=df(rhs e1,e2)
>>
>>;
dv!/dx:=dv!/dx**(-1);
% Substitution of all derivatives
for each e1 in sb do
for each e2 in e1 do <<
if not freeof(eqlist,lhs e2) then <<
% which function is to be differentiated wrt. which variable
eqlist:=sub(lhs e2=transderiv(rhs e2,yslist,vlist),eqlist)
>>
>>;
clear dv!/dx;
return sub(xslist,sub(yslist,eqlist))$
end$ % of DeTrafo
%----------------------------
algebraic procedure grouping(el1,el2,xlist,ylist,nx,ny)$
begin scalar h,el3,xslist,yslist$
%------- Grouping the new variables to ulist and vlist
h:={};
xslist:={}; % list of expressions to calculate new indep. var.
yslist:={}; % list of expressions to calculate new dep. var.
%---- at first the obvious allocations
for each el3 in el1 do %-- all similarity variables
if freeoflist(el3,ylist) then xslist:=cons(el3,xslist) else
if freeoflist(el3,xlist) then yslist:=cons(el3,yslist) else
h:=cons(el3,h);
%---- now the symmetry variable
if freeoflist(el2,ylist) or (length(yslist) = ny) then
xslist:=cons(el2,xslist) else
if freeoflist(el2,xlist) or (length(xslist) = nx) then
yslist:=cons(el2,yslist) else
xslist:=cons(el2,xslist);
%---- now the remaining cases
for each el3 in h do
if length(yslist) < ny then yslist:=cons(el3,yslist)
else xslist:=cons(el3,xslist);
return {xslist,yslist}
end$ % of grouping
%----------------------------
algebraic procedure rename_u_(xslist,yslist,el2,u_,u,v)$
begin scalar i,vlist,ulist,el3,h,smv$
%---- Renaming the u_ to ui in yslist and to vi in xslist
i:=0;
vlist:={};
xslist:=for each el3 in xslist collect
<<i:=i+1;
if length xslist>1 then h:=mkid(v,i)
else h:=v;
vlist:=cons(h,vlist);
if el3=el2 then smv:=h;
sub(u_=h,el3)
>>;
i:=0;
ulist:={};
yslist:=for each el3 in yslist collect
<<i:=i+1;
if length yslist>1 then h:=mkid(u,i)
else h:=u;
ulist:=cons(h,ulist);
if el3=el2 then smv:=h;
sub(u_=h,el3)
>>;
return {xslist,yslist,reverse vlist,reverse ulist,smv}
end$ % of rename_u_
%----------------------------
algebraic procedure solve_for_old_var(xslist,yslist,xlist,ylist,nx,ny)$
begin scalar h1,h2$
%---- Solve for old variables
h1:=nil;
h2:=solve(append(yslist,xslist),append(xlist,ylist));
if h2={} then h1:=t
else h2:=first h2; %--- possibly other solutions
if list neq lisp(car algebraic h2) then h1:=t else
if length(h2)<(nx+ny) then el2:=t;
if h1 then repeat lisp
<<write"The algebraic system ",append(xslist,yslist),
" could not be solved for ",append(xlist,ylist),".";
write"Please enter the solution in form of a list {",
reval algebraic first xlist,"=...,...",
reval algebraic first ylist,"=...,...} or enter a ",
"semicolon ; to end this investigation:";
algebraic(h2:=aread(nil))
>> until h2=nil or ( lisp(pairp algebraic h2) and
(list=lisp(car algebraic h2)) and
(length(h2)=(nx+ny)) )
else
<<lisp<<terpri()$
write"The suggested solution of the algebraic system which will";
terpri()$
write"do the transformation is: ";
terpri()
>>;
write h2;
if yesp "Is the solution ok?" then else lisp <<
write"Please enter the solution in form of a list {",
reval algebraic first xlist,"=...,...",
reval algebraic first ylist,"=...,...} or enter a ",
"semicolon ; to end this investigation:";
algebraic(h2:=aread(nil))
>>
>>;
return h2
end$ % of solve_for_old_var$
%----------------------------
algebraic procedure switch_r_s(h2,smv,ylist,u,v)$
begin scalar xslist,yslist,el3,h$
%---- Exchange of dependent and independent variables
xslist:={};
yslist:={};
for each el3 in h2 do if freeof(ylist,lhs el3) then
xslist:=cons(el3,xslist) else
yslist:=cons(el3,yslist);
lisp <<terpri()$
write"In the intended transformation shown above",
" the dependent ";terpri()$
if length yslist>2 then
write"variables are the ",reval algebraic u,"i and " else
write"variable is ",reval algebraic u," and ";
if length xslist>2 then
write"the independent variables are the ",
reval algebraic v,"i." else
write"the independent variable is ",reval algebraic v,".";
terpri()$
write
"The symmetry variable is ",reval algebraic smv,", i.e. the ",
"transformed expression";terpri();
write"will be free of ",reval algebraic smv,".";
>>;
h:=if yesp
"Is this selection of dependent and independent variables ok?"
then nil else <<
lisp <<write"Please enter a list of substitutions. For example, to";
terpri()$
write"make the variable, which is so far call u1, to an";
terpri()$
write"independent variable v2 and the variable, which is ";
terpri()$
write"so far called v2, to an dependent variable u1, ";
terpri()$
write"enter: `{u1=v2, v2=u1};'">>;
aread(nil)
>>;
if h and (h neq {}) then <<xslist:=sub(h,xslist);
yslist:=sub(h,yslist);
smv :=sub(h,smv)>>;
return {xslist,yslist,smv}
end$ % of switch_r_s
%----------------------------
algebraic procedure makepde(genlist,u_)$
begin scalar h,el2,el3,vari,bv;
vari:={};
return
cons(
num for each el2 in genlist sum
<<h:=lhs el2;
h:=lisp <<
el3:=explode reval algebraic h;
bv:=t;
while bv do <<
if car el3 ='!_ then bv:=nil;
el3:=cdr el3
>>;
intern compress el3
>>;
depend u_,h;
vari:=cons(h,vari);
(rhs el2) * df(u_,h)
>>,
vari)
end$ % of makepde
%----------------------------
algebraic procedure totdeglist(eqlist,ylist)$
begin scalar n,ordr,e1,e2;
ordr:=0;
for each e1 in eqlist do
for each e2 in ylist do
<<n:=totdeg(e1,e2);
if n>ordr then ordr:=n>>;
return ordr
end$ % of totdeglist
%----------------------------
algebraic procedure similarity(problem,genlist,con,u,v)$
% con ... the free constants/functions in the general symmetry
% u ... the name of the new independent variables
% v ... the name of the new dependent variables
begin scalar vari,pde,el1,el2,el3,el4,copgen,symvarfound,
trans1,trans2,i,j,h,h2,n,denew,xlist,ylist,
eqlist,ulist,vlist,nx,ny,xslist,yslist,smv,
trafoprob;
cpu:=lisp time()$ gc:=lisp gctime()$
%--------- extracting input data
eqlist:=maklist first problem;
ylist :=maklist second problem; ny:=length ylist;
xlist :=maklist third problem; nx:=length xlist;
trafoprob:={problem,nil}; % to be returned if trafo not possible
problem:=nil;
eqlist:=for each el1 in eqlist collect equ_to_expr el1;
if length eqlist > 1 then eqlist:=desort eqlist;
%--------- initial printout
lisp(
% if print_ and logoprint_ then <<
% write "-----------------------------------------------",
% "---------------------------"$ terpri()$terpri()$
% write "This is SIMILARITY - a program for calculating similarity";
% write
% "variables and performing corresponding point-transformations";
% write "for given infinitesimal symmetries.";
% terpri()
% >> else
if print_ then terpri());
%--------- initializations
ordr:=totdeglist(eqlist,ylist)$
vari:=append(ylist,xlist);
for each el1 in xlist do
for each el2 in ylist do
nodepend el2,el1;
lisp(
if print_ then <<
write "The ODE/PDE (-system) under investigation is :";terpri()$
for each el1 in cdr eqlist do algebraic write"0 = ",el1;
terpri()$write "for the function(s) : ";
fctprint( cdr reval algebraic ylist);write".";
terpri()$terpri()
>>);
lisp(
if print_ then <<
if length ylist >2 then % not >1 because of alg. list in symb. mode
write"It will be looked for new dependent variables ",u,"i "
else
write"It will be looked for a new dependent variable ",u;
terpri()$
if length xlist >2 then
write"and independent variables ",v,"i"
else
write"and an independent variable ",v;
write" such that the transformed";
terpri()$
write"de(-system) does not depend on ",u;
if length ylist >2 then write"1";
write" or ",v$
if length xlist >2 then write"1";
write".";
terpri()
>>);
% for each el1 in con do % i.e. for all symmetries do:
% if freeoflist(el1,xlist) and
% freeoflist(el1,ylist) then %------ no pseudo-Lie-symm.
% <<copgen:=sub(el1=1,genlist); %------ calculate symmetry
% for each el2 in con do
% if el1 neq el2 then copgen:=sub(el2=0,copgen);
copgen:=genlist;
% write"The symmetry now under investigation is:";
% for each el1 in copgen do write el1;
%---------- formulate the PDE for the similarity variables
pde:=first makepde(copgen,u_);
%--------- find similarity variable
trans2 :={};
lisp<<terpri()$write"1. Determination of the similarity variable";
if nx+ny>2 then write"s">>;
trans1 := quasilinpde1(pde,u_,vari); % for the similarity variable
if trans1 neq {} then
<<
%-------------- Determining the similarity variables ui_
i:=0;
trans1:=for each el1 in trans1 collect
<<i:=i+1;
h:=length(genlist)-1;
if h=1 then % one single ODE
<<
% write"In the following 1 similarity variable U_ has to be";
% write"determined through 0 = ff where ff is an"$
% write"arbitrary function of arguments given in the ",
% "following list:"$
% write el1;
el2:=num(first el1 - second el1);
if freeof(el2,u_) then
el2:=num(first el1 - 2*second el1);
lisp<<write"A suggestion for this function ff provides:"$
terpri()>>$
write"0 = ", el2$
if yesp "Do you like this choice?" then {el2}
else <<
repeat <<
lisp <<
write"Put in an alternative expression which "$terpri()$
write"- is functionally dependent only on elements of",
" ff given above and"$ terpri()$
write"- depends on U_ and if set to zero determines U_"$
terpri()>>$
h:=aread(nil)$
>> until not freeof(h,u_)$
{h}
>>
>> else % a PDE or a system of DEs
lisp
<<
% terpri()$
% write"Now the similarity variables U_i, i=1,...",h,
% " have to be"$terpri()$
% write"determined through conditions 0 = ffi, i=1,...",h,
% ", where ffi are"$terpri()$
% write"arbitrary functions ";terpri()$
% algebraic(write "ffi = ",
% lisp( cons('ffi,cdr reval algebraic el1)));
% terpri()$
% write"such that the functional determinant of these ",
% "expressions"$ terpri()$
% write"including u_ from above taken w.r.t. ";
% for each el3 in cdr algebraic ylist do write reval el3,",";
% write reval cadr algebraic xlist;
% for each el3 in cddr algebraic xlist do write ",",reval el3;
% terpri()$
% write"must not vanish."$terpri()$
algebraic <<
el2:=einfachst(el1,u_);
h2:={};
for each el3 in el1 do
if el3 neq el2 then h2:=cons(num(el2-el3),h2)
>>;
write
"A suggestion for these functions ffi in form of a list ",
"{ff1,ff2,... } is: "$terpri()$
deprint(cdr reval algebraic h2);
if yesp "Do you like this choice?" then algebraic h2
else <<
write"Put in an alternative list of expression which"$
terpri()$
write"- are functionally dependent only on the above ",
"arguments and"$terpri()$
write"- which if set to zero determine U_i, i.e."$terpri()$
write"- the functional determinant of these expressions"$
terpri()$
write" including U_ from above taken w.r.t. ",
cdr reval algebraic append(ylist,xlist)$
terpri()$
write" must not vanish."$terpri()$
algebraic(h2):=aread(nil)
>>
>>
>>$
%--------- find symmetry variable
pde:=pde-1;
lisp<<terpri()$write"2. Determination of the symmetry variable">>;
trans2 := quasilinpde1(pde,u_,vari); % for the symmetry variable
if (length xlist=1) and (trans2={}) then
for each e1 in fargs u_ do nodepend u_,e1
else
<<% If no symmetry variable is found (trans2={}) then proceed
% only if special solutions of PDEs are to be found with the
% solution being only a function of the similarity variables
if trans2={} then << % take any variable
h:=reverse vari; % reverse to have not a function as symvar
while 1+sub(u_=first h,pde)=0
do h:=rest h; % no similarity var.
for each e1 in fargs u_ do nodepend u_,e1;
h:=first h;
trans2:={u_ - h};
lisp<<write"Because the correct symmetry variable was not ",
"found, the program will";terpri()$
write"take ",reval algebraic h,
" instead with the consequence ",
"that not the whole transformed ";terpri()$
write"PDE will be free of ",
reval algebraic h," but only those ",
"terms without ",
reval algebraic h,"-derivative";terpri()$
write"which is still of use for finding special ",
reval algebraic h,"-independent solutions ";
terpri()$
write"of the PDE."
>>
>> else <<
%--------------- Determining an optimal symmetry variable
for each e1 in fargs u_ do nodepend u_,e1;
symvarfound:=t;
i:=0;
trans2:=for each el1 in trans2 collect
<<i:=i+1;
lisp
<<
% terpri()$
% write"In the following the symmetry variable U_",
% " has to be"$terpri()$
% write"determined through a condition 0 = ff",
% ", where ff is"$terpri()$
% write"an arbitrary function ";terpri()$
% algebraic(write"ff = ",
% lisp( cons('ff,cdr reval algebraic el1)));terpri()$
write
"A suggestion for this function ff(..) yields:";terpri()
>>;
h:=einfachst(el1,u_);
if lisp<<h2:=reval algebraic(num h);
(not pairp h2) or (car h2 neq 'plus)>> then
h:=num(h+1);
% if h= first el1 then
% if freeof(num(h+second el1),u_) then h:=num(h+2*second el1)
% else h:=num(h+ second el1)
% else
% if freeof(num(h+ first el1),u_) then h:=num(h+2* first el1)
% else h:=num(h+ first el1);
write"0 = ",h$
if yesp "Do you like this choice?" then h
else <<
repeat <<
lisp <<
write
"Put in an alternative expression which "$terpri()$
write
"- is functionally dependent only on arguments of",
" ff given above and"$terpri()$
write
"- depends on u_ and if set to zero determines u_"$
terpri()>>$
h:=aread(nil)$
>> until not freeof(h,u_)$
h
>>
>>
>>$
%for each el1 in trans1 do
%for each el2 in trans2 do
el1:=first trans1;
el2:=first trans2;
% <<
%------- Grouping the new variables to ulist and vlist
yslist:=grouping(el1,el2,xlist,ylist,nx,ny)$
xslist:= first yslist;
yslist:=second yslist;
%---- Renaming the u_ to ui in yslist and to vi in xslist
smv:=rename_u_(xslist,yslist,el2,u_,u,v)$
xslist:= first smv;
yslist:=second smv;
vlist := third smv;
ulist:=part(smv,4);
smv:=part(smv,5); % the symmetry variable
%---- Solve for old variables
h2:=solve_for_old_var(xslist,yslist,xlist,ylist,nx,ny);
if h2 neq nil then <<
%---- Exchange of dependent and independent variables
smv:=switch_r_s(h2,smv,ylist,u,v)$
xslist:= first smv;
yslist:=second smv;
smv :=third smv;
%---- Doing the point transformation
for each el3 in ulist do
<<for each el4 in fargs el3 do nodepend el3,el4;
for each el4 in vlist do
%if el4 neq smv then % if new DEs without symm. var. smv
depend el3,el4>>;
for each el3 in ylist do
for each el4 in xlist do depend el3,el4;
eqlist:=detrafo(eqlist,yslist,xslist,ulist,vlist);
lisp(
if print_ then <<
terpri()$
write"The transformed equation";
if length(algebraic eqlist)>2 then write"s";
if symvarfound then
write" which should be free of ",reval algebraic smv,":"
else
write" in which the terms without ",reval algebraic smv,
"-derivative are free of ",reval algebraic smv,":";
terpri()
>>);
eqlist:=for each el3 in eqlist collect <<
el3:=factorize num el3;
for each el4 in el3 product
if 0=totdeglist({el4},ulist) then 1
else el4
>>;
lisp deprint(cdr reval algebraic eqlist);
if (length(vlist)>1) and (not freeof(vlist,smv)) then <<
vlist:=cons(smv,lisp(delete(reval algebraic smv,
reval algebraic vlist)));
if yesp
"Shall the dependence on the symmetry variable be dropped?"
then
<<for each el3 in ulist do nodepend el3,smv;
vlist:=rest vlist>>;
eqlist:=for each el3 in eqlist collect <<
el3:=factorize num el3;
for each el4 in el3 product
if 0=totdeglist({el4},ulist) then 1
else el4
>>
>>;
trafoprob:={{eqlist,ulist,vlist},append(xslist,yslist)}
>>
% >>
>>
>>;
% >>;
for each el1 in xlist do
for each el2 in ylist do
depend el2,el1;
clear ff,ffi;
return trafoprob
end$ % of similarity
%----------------------------
symbolic procedure simplifiziere(q)$
begin
scalar n,one_argument_functions_;
one_argument_functions_:= '(abs acos acosd acosh acot acotd acoth acsc
acscd acsch asec asecd asech asin asind asinh
atan atand atanh cbrt cos cosd cosh cot
cotd coth csc cscd csch exp hypot ln
log log10 sec secd sech sin sind sinh sqrt
tan tand tanh minus)$
return
if not pairp q then q else
if member(car q,one_argument_functions_) or
((car q = 'expt) and (numberp caddr q)) then simplifiziere(cadr q)
else <<
n:=ncontent(q);
if n=1 then q
else simplifiziere(reval list('quotient, q, reval n))
>>
end$
%----------------------------
algebraic procedure select_indep_var(clist,xlist)$
begin scalar s,scal,notfound,cq,x,xx,xanz,sanz,xcop,ok;
% Find the simplest non-zero element of clist
notfound:=t;
xcop:=xlist;
while notfound and (xcop neq {}) do <<
x :=first xcop ; xcop :=rest xcop ;
cq:=first clist; clist:=rest clist;
if cq neq 0 then <<
xanz:=0;
for each xx in xlist do
if %df(cq,xx)
not freeof(cq,xx) then xanz:=xanz+1;
ok:=nil;
if not s then ok:=t else % to guarantee s neq nil
if xanz<sanz then ok:=t else % as few dependencies as poss.
if xanz=sanz then
if length(cq)=1 then ok:=t else
if length(scal)>1 then
if part(scal,0)=plus then % if possible not a sum
if (part(cq,0) neq plus) or
(length(cq)<length(scal)) then ok:=t;
if ok then <<s:=x; scal:=cq; sanz:=xanz>>;
if scal=1 then notfound:=nil
>>
>>;
return {s,scal}
end$ % of select_indep_var
%----------------------------
procedure charsyspsfi(xlist,cqlist,u,ode,denf);
begin
scalar h,s;
h:=cqlist;
cqlist:={};
while h neq {} do <<cqlist:=cons(first h*denf,cqlist);h:=rest h>>;
cqlist:=cons(-ode*denf,cqlist);
xlist:=cons(u,reverse xlist);
s:=lisp gensym();
for each h in xlist do depend h,s;
h:=psfi(cqlist,xlist,s);
write "soln=",h;
for each h in xlist do nodepend h,s;
return h
end$ % of charsyspsfi
%----------------------------
algebraic procedure charsyscrack(xlist,cqlist,s,scal,u,ode)$
% formulation and solution of the characteristic ODE-system
begin
scalar lcopy,x,cs,flist,soln,printold,timeold,facintold,
e1,e2,e3,e4,n,nmax,dff;
% formulating the ode/pde - system .
lcopy := cqlist ;
cs := {} ;
for each x in xlist do
<< cq := first lcopy ; lcopy := rest lcopy ;
if x neq s then
<< depend x,s;
cs := .(scal*df(x,s)-cq,cs) >>
>>;
if s neq u then <<flist:={u}; depend u,s;
cs := .(scal*df(u,s)+ode,cs) >>
else flist:={};
for each x in xlist do
if s neq x then flist:=.(x,flist);
% solving the odes using crack.
if lisp(print_ neq nil) then
lisp <<
write "The equivalent characteristic system:";terpri();
deprint cdr algebraic cs$
%terpri()$
write "for the functions: ";
fctprint( cdr reval algebraic flist);write".";
printold:=print_; print_:=nil;
timeold := time_; time_ :=nil;
facintold:=facint_;facint_:=t;
>>;
% lisp(print_:=100); %############
soln:=crack(cs,flist,flist,{});
lcopy:={};
for each x in soln do <<
e1:=first x;
if (e1 neq {}) and (length e1 + length second x = length flist) then
% all remaining conditions are algebraic (non-differential) in all
% the functions of flist?
<<e2:={};
e3:={};
for each e4 in third x do if freeof(flist,e4) then e3:=e4 . e3
else e2:=e4 . e2;
if (length cs) = (length e3) then <<
% sufficiently many integrations done
for each e4 in e2 do
for each e5 in e1 do
if lisp(not freeof(lderiv(reval algebraic e5,
reval algebraic e4,
list(reval algebraic s)),'df))
then <<e2:={};e1:={}>>;
if e2 neq {} then <<
% It may be possible that derivatives of unsolved functions
% occur in the expressions of the evaluated functions: second x
nmax:=0;
for each e4 in e2 do << % for each unsolved function
for each e5 in second x do << % for each solved expression
n:=lisp(lderiv(reval algebraic rhs e5,reval algebraic e4,
list(reval algebraic s)));
if (n neq 0) and (n neq nil) then
if lisp length algebraic n = 3 then n:=1
else n:=lisp cadddr algebraic n;
if n>nmax then nmax:=n;
>>;
if nmax>0 then << % adding further conditions
e5:=e1;
while freeof(first e5,e4) do e5:=rest e5;
e5:=first e5;
dff:=e4;
for n:=1:nmax do <<
e5 :=df(e5,s); e1:= e5 . e1;
dff:=df(dff,s); e3:=dff . e3
>>
>>
>>;
lcopy:=cons({append(second x,e1),e3},lcopy);
>>
>>
>> else
if (first x = {}) and (length cs = length third x) then
lcopy:=cons({second x,third x},lcopy)
>>;
lisp <<print_:=printold; time_:=timeold; facint_:=facintold>>;
return
if lcopy={} then <<for each x in flist do nodepend x,s; {}>>
else s . lcopy
% { {{x=..,y=..,u=..,..,0=..},{df(z,s),df(z,s,2),..,c1,c2,c3,..}},..}
% df(z,s,..) only if df(z,s,..) turns up in x, y, u. .. .
end$ % of charsyscrack
%----------------------------
algebraic procedure quasilinpde1(pde,u_,vari)$
begin scalar trans1,e1,e2,q;
trans1 := quasilinpde(pde,u_,vari); % for the similarity variable
if trans1={} then <<
write"The program was not able to find the general solution ",
"of the PDE: ",pde," for the function ",u_,".";
lisp <<
write"Please enter either only a semicolon if no solution ",
"is known or enter "$terpri()$
write"the solution of the PDE in form ",
"of a list {A1,A2,...} where ";terpri()$
write"the Ai are algebraic expressions in ",
cdr reval algebraic cons(u_,vari);terpri()$
write"such that any function ff(A1,A2,...) which is not ",
"independent of `",u_;write"'";terpri()$
write"determines a solution `",u_,"' of the PDE through 0=ff: "
>>;
trans1:={aread(nil)};
if trans1={nil} then trans1:={}
>>;
return trans1
end$ % of quasilinpde1
%----------------------------
algebraic procedure quasilinpde(f,u,xlist);
% f ... PDE, u ... unknown function, xlist ... variable list
begin scalar i, q, qlist, cq, cqlist, ode, soln, tran,
xcop, s, s1, x, xx, scal, qlin, prev_depend;
symbolic put('ff,'simpfn,'simpiden)$
% if lisp(print_ neq nil) then
write"The quasilinear PDE: 0 = ",f,".";
% changing the given pde into a quasi-linear ode .
i := 0 ;
ode := f;
qlist := {}; cqlist:={};
for each x in xlist do
<<i := i+1 ;
q := mkid(p_,i) ;
qlist := .(q,qlist) ;
f := sub(df(u,x) = q , f) ;
cq:=df(f,q);
cqlist:=.(df(f,q),cqlist) ;
ode := ode - df(u,x)*df(f,q) ;
nodepend u, x ;
>> ;
qlist := reverse qlist ;
cqlist := reverse cqlist ;
% checking for linearity.
qlin:=t;
for each cq in cqlist do
for each q in qlist do
if df(cq,q) neq 0 then qlin:=nil;
if not qlin then return {};
% soln:=charsyspsfi(xlist,cqlist,u,ode,den f);
% Determination of the independent variable for the ODEs
pcopy:=cons(-ode,cqlist);xcop:=cons(u,xlist);
scal:=select_indep_var(pcopy,xcop)$
s1:=first scal;
prev_depend:=storedepend(xlist)$
soln:=charsyscrack(xlist,cqlist,s1,second scal,u,ode);
if soln={} then << % try all other coefficients as factors
repeat <<
repeat <<s :=first xcop ;xcop :=rest xcop;
scal:=first pcopy;pcopy:=rest pcopy>>
until (pcopy={}) or ((scal neq 0) and (s neq s1));
if (s neq s1) and (scal neq 0) then << lisp
<<write"New attempt with a different independent variable";
terpri()>>$
soln:=charsyscrack(xlist,cqlist,s,scal,u,ode)
>>
>>
until (soln neq {}) or (xcop={})
>>;
% solving for the constants(eg..c1,c2 etc.) and put them in a
% linear form.
% The first element of soln is the ODE-variable
tran:={};
if soln neq {} then <<
s1:=first soln;
for each s in rest soln do
<<cq:=first solve(first s,second s);
tran:=.(<<x:={};
for each xx in cq do
if lisp atom algebraic lhs xx then
x:=.(lisp <<q:=reval algebraic rhs xx;
simplifiziere q >> ,x);
x
>>,
% intsimplify(q,cons(reval algebraic u,
% cdr reval algebraic xlist)) >>,
tran);
>>;
for each s in xlist do if s neq s1 then nodepend s,s1
>>;
for each x in xlist do depend u,x;
if lisp(print_ neq nil) then
if tran neq {} then <<
write"The general solution of the PDE is given through";
for each x in tran do write"0 = ",
lisp( cons('ff,cdr reval algebraic x));
if length tran>1 then write"with arbitrary function(s) ff(..)."
else write"with arbitrary function ff(..)."
>>;
% restoring the dependencies of the variables of the PDE
restoredepend(prev_depend)$
return tran;
end$ % of quasilinpde
endmodule;
end;