%*********************************************************************
module integration$
%*********************************************************************
% Routines for integration of pde's
% Author: Thomas Wolf, Andreas Brand
% 1993 1995
%
% $Id: crint.red,v 1.13 1998/06/25 18:21:25 tw Exp tw $
%
symbolic procedure ldlist(p,f,vl)$
% provides a reverse list of leading derivatives of f in p, vl is list
% of variables
begin scalar a$
if not atom p then
if member(car p,list('EXPT,'PLUS,'MINUS,'TIMES,'QUOTIENT,'DF,'EQUAL))
then <<
if (car p='PLUS) or (car p='TIMES) or
(car p='QUOTIENT) or (car p='EQUAL) then
<<p:=cdr p$
while p do
<<a:=diffincl(a,ldlist(car p,f,vl))$
p:=cdr p>>
>> else
if car p='MINUS then a:=ldlist(cadr p,f,vl) else
if car p='EXPT then % if numberp caddr p then
a:=ldlist(cadr p,f,vl) else % fuehrende Abl. aus der Basis
% else a:=nil
if car p='DF then if cadr p=f then
<<p:=cddr p;
while vl do
<<a:=cons(dfdeg(p,car vl),a);
vl:=cdr vl>>;
a:=list a
>>
>>$
return a
end$
symbolic procedure diffincl(a,b)$
% a,b are lists of leading derivatives which are to be united
begin
scalar p;
while a and b do
<<a:=ddroplow(a,car b);
if car a then p:=cons(car a,p);
a:=cdr a;
b:=cdr b>>;
return
if null a then if p then nconc(p,b)
else b
else if p then a:=nconc(p,a)
else a
end$
symbolic procedure ddroplow(a,cb)$
% loescht Elemente von a, von denen cb eine Ableitung ist, loescht cb,
% wenn ein Element von a eine Ableitung von cb ist
begin
scalar h;
return
if null a then list(cb)
else
<<h:=compdiffer(car a,cb);
if numberp(h) then if h>0 then cons(nil,a) % car a=df(cb,..
else ddroplow(cdr a,cb) % cb=df(car a,..
else <<h:=ddroplow(cdr a,cb); % neither
cons(car h,cons(car a,cdr h))>>
>>
end$
symbolic procedure compdiffer(p,q)$
% finds whether p is a derivative of q or q of p or neither
begin
scalar p!>q,q!>p;
while p and ((null p!>q) or (null q!>p)) do
<<if car p>car q then p!>q:=t else % compare orders of derivatives
if car p<car q then q!>p:=t;
p:=cdr p;q:=cdr q
>>;
return
if q!>p then if p!>q then nil % neither
else 0 % q is derivative of p
else if p!>q then 2 % p is derivative of q
else 1 % p equal q
end$
symbolic procedure ldintersec(a)$
% determines the intersection of derivatives in the list a
begin
scalar b;
return
if null a then nil else
<<b:=car a;a:=cdr a;
while a do
<<b:=isec(b,car a)$
a:=cdr a
>>;
b
>>
end$
symbolic procedure isec(ca,b)$
% determines the minimum derivatives between ca and b
begin
scalar newb;
while ca do
<<newb:=cons(if car b<car ca then car b else car ca, newb);
ca:=cdr ca;b:=cdr b
>>;
return reverse newb
end$
symbolic procedure disjun(a,b)$
<<while a do
if (car a neq 0) and (car b neq 0) then a:=nil
else <<a:=cdr a;b:=cdr b>>;
if b then nil else t
>>$
symbolic procedure ddrophigh(a,cb)$
% loescht Elemente von a, die Ableitung von cb sind,
% loescht cb, wenn cb Ableitung eines Elements von a ist oder =a ist,
% haengt cb ansonsten an
begin
scalar h;
return
if null a then list(cb)
else
<<h:=compdiffer(car a,cb);
if numberp(h) then if h<2 then a % cb is deriv or = car a
else ddrophigh(cdr a,cb) % car a=df(cb,..
else cons(car a,ddrophigh(cdr a,cb)) % neither
>>
end$
symbolic procedure elibar(l1,l2,lds)$
begin
scalar found1,found2,h;
% found1=t if an LD=l1 is found, found2=t if contradiction found
while lds and (not found2) do
<<if car lds=l1 then found1:=t else
<<h:=compdiffer(l2,car lds);
if (null h) or (h eq 2) then found2:=t
>>;
lds:=cdr lds
>>;
return found1 and (not found2)
end$
symbolic procedure intminderiv(p,ftem,vlrev,maxvanz,nfsub)$
% yields a list {nv1,nv2, ... } such that nvi is the minimum
% of the highest derivatives w.r.t. vi of all the leading derivatives
% of all functions of ftem which are functions of all maxvanz variables.
% Only those are kept for which nvi>0.
% further a list ((f1 ld's of f1) (f2 ld's of f2)...),
begin scalar l,a,listoflds$
while ftem do
<<if (maxvanz=(fctlength car ftem)) or (nfsub=0) then
<<l:=ldlist(p,car ftem,vlrev);
listoflds:=cons(cons(car ftem,l),listoflds)$
a:=if a then ldintersec(cons(a,l))
else ldintersec(l)
>>$
ftem:=cdr ftem
>>$
return list(a,listoflds)
end$
symbolic procedure potintegrable(listoflds)$
begin
scalar h,l1,l2,f,n,lds,f1,f2$
if tr_genint then write "Does a potential exist?"$
%------- Determining 2 minimal sets of integration variables
% There must be two disjunct LDs such that all others are in their
% ideal. This is necessary if we want to eliminate 2 functions.
h:=listoflds;
l1:=nil;
while h do
<<l2:=cdar h; % the list of LDs for the function caar h
while l2 do
<<l1:=ddrophigh(l1,car l2)$
l2:=cdr l2>>;
h:=cdr h
>>;
return
if length l1 neq 2 then nil else
if not disjun(car l1,cadr l1) then nil else
% if there would be an overlap between l1 and l2 then it would have
% to be integrated for elimination but it cannot --> no overlap
% possible
<<% selecting interesting functions for which one LD is = l1 and all
% others are derivatives of l2 or equal l2 and vice versa. Two
% necessary (one with an LD=l1 and one with an LD=l2) from which at
% least one function (f) has no further LD.
% Exception: 2 functions with each 2 LDs equal to (l1,l2) (these
% functions are counted by n)
l2:=cadr l1;l1:=car l1;
f:=nil;f1:=nil;f2:=nil;n:=0;
h:=listoflds;
while h and ((not f1) or (not f2) or ((not f) and (n neq 2))) do
<<lds:=cdar h;
if (not f1) or (not f) then
if elibar(l1,l2,lds) then
<<f1:=cons(caar h,f1);
if length lds eq 1 then f:=caar h else
if length lds eq 2 then
if (car lds = l2) or (cadr lds = l2) then n:=n+1
>>;
if (not f2) or (not f) then
if elibar(l2,l1,lds) then
<<f2:=cons(caar h,f2);
if length lds eq 1 then f:=caar h
>>;
h:=cdr h
>>;
if f1 and ((n>1) or (f2 and f)) then list(l1,l2)
else nil
>>
end$ % of potintegrable
symbolic procedure vlofintlist(vl,intlist)$
% provides a list of elements of vl for which the corresponding
% elements of intlist are non-zero
begin scalar a;
while intlist do
<<if (car intlist) and (not zerop car intlist) then a:=cons(car vl,a);
vl:=cdr vl;
intlist:=cdr intlist
>>;
return a
end$
symbolic procedure vlofintfaclist(vl,intfaclist)$
% determining the list of all variables of vl in intfaclist
begin scalar e1,a;
for each e1 in vl do
if not my_freeof(intfaclist,e1) then a:=cons(e1,a);
return a
end$
symbolic procedure multipleint(intvar,ftem,q,vari,vl,genflag,
potflag,partial,doneintvar)$
% multiple integration of q wrt. variables in vl, max. number of
% integrations specified by intvar
% integrating factors must not depend on doneintvar, doneintvar is
% a list of integration variables so far
% partial=t then as much as possible of an expression is integrated
% even if there is a remainder
begin
scalar pri,vlcop,v,nmax,geni,intlist,iflag,n,nges,newcond,
intfaclist,ph,pih,qh,qih,intfacdepnew,intfacdep$
% intfacdep is a list of variables on which factors of integration
% depend so far, other than the integration variable in their
% integration --> no integration wrt. these variables by potint
% because there the diff. operators wrt. to different variables
% need not commute because the integrations are not done
% pri:=t$
if (not vari) and (zerop q) then return nil;
nges:=0;
vlcop:=vl;
pih:=t;
% Splitting of the equation into the homogeneous and inhomogeneous
% part which is of advantage for finding integrating factors
q:=splitinhom(q,ftem,vl)$
qh:=car q; qih:=cdr q; q:=nil;
while (vari or vlcop) and (pih or (not potflag)) do
%------- if for potflag=t one variable can not be integrated the
%------- maximal number of times (nmax) then immediate stop because
%------- then no elimination of two functions will be possible
<< %-- The next integration variable: x, no of integrations: nmax
if vari then <<v:=vari;nmax:=1>>
else <<v:=car vlcop; vlcop:=cdr vlcop;
nmax:=car intvar; intvar:=cdr intvar>>;
if zerop nmax then intlist:=cons(nil,intlist)
else
<<if pri then write"anf: intvar=",intvar," vari=",vari," q=",q$
if vari and (not member(v,vl)) then
<<qh :=reval list('INT,qh ,v)$
if freeof(qh,'INT) then <<
qih:=reval list('INT,qih,v)$
iflag:=if freeint_ and
(null freeof(qih,'INT)) then nil
else <<
intlist:=cons(list(1),intlist)$
'success>>$
if pri then <<write"232323 qh=",qh;terpri();
write"qih=",qih;terpri()>>
>>
>> else
<<n:=0$
if pri then write"333"$
intfaclist:=nil; %-- the list of integr. factors in v-integr.
if potflag or my_freeof(intfacdep,v) then
% otherwise v-integration not allowed because one integrating
% factor already depends on v
% for potflag=t this `commutativity'-demand plays no role
repeat << %--- max nmax integrations of qh and qih wrt. v
if pri then <<write"444 vor intpde:"$eqprint q$terpri()$
write"potflag=",potflag," v=",v,
" ftem=",ftem>>$
% At first trying a direct integration of the homog. part qh
ph:=intpde(qh,ftem,vl,v,partial)$ % faster if potflag=nil
if pri then <<write"nach intpde(qh):"$deprint ph>>$
%------ At first the integration of the homogeneous part
intfacdepnew:=intfacdep;
if ph and (partial or (zerop cadr ph)) then <<
%---- For the homogen. part cadr ph must be zero
intfaclist:=cons(1,intfaclist);
ph:=car ph;
if pri then <<write"565656 ph=",ph;terpri()>>;
>> else
if vari then ph:=nil
else
if facint_ then <<
ph:=findintfac(list(qh),ftem,vl,v,doneintvar,intfacdep,
not zerop n,not potflag);
% factorize before ivestig., no report of int. factors
if ph then << %--- Complete integr. of qh was possible
if print_ then write"of the homogeneous part"$terpri()$
%--- update the list of variables on which all integr.
%--- factors depend apart from the integration variable
intfacdepnew:=caddr ph;
%--- extend the list of integrating factors, cadr ph
%--- is a list of integr. factors, here only one
intfaclist:=cons(caadr ph,intfaclist);
%--- multiply the inhomogeneous part with integ. factor
qih:=reval reval reval list('TIMES,car intfaclist,qih);
if pri then <<write"454545 qih=",qih;terpri()>>;
ph:=car ph % the integral of qh
>>
>>;
%------ Now the integration of the inhomogeneous part
if not ph then pih:=nil %--- no integration possible
else <<
if zerop qih then pih:=list(0,0) else
pih:=intpde(qih,ftem,vl,v,partial)$
if null pih then
write"error2: inhom. expr. can not be integrated"$
if pri then <<write"nach intpde(qih):",pih$terpri()$
write"genflag=",genflag$terpri()>>$
if pih then
if zerop cadr pih then
<<qih:=car pih$n:=add1 n$iflag:='success$
if pri then write" success "$
>>
else if not genflag then pih:=nil
else
<<if pri then write"555"$
geni:=partint(cadr pih,smemberl(ftem,cadr pih),
vl,v,genflag)$
if geni then
<<qih:=reval list('PLUS,car pih,car geni)$
n:=add1 n$
ftem:=union(fnew_,ftem)$
newcond:=append(cdr geni,newcond)$ % additional de's
if pri then
<<terpri()$write"nach gen newcond:",newcond>>$
iflag:='genint
>> else
if partial then qih:=car pih
else pih:=nil
>>;
if pih then <<
if pri then write"AAA"$
qh:=ph;
if (not potflag) and (n neq 1) then
intfacdep:=intfacdepnew
%-The first integr. factor of all v-integrations does not
% depend on any earlier integration variables and can
% therefore be taken out of all integrals --> no incease
% of intfacdep for n=1.
%-For potential integration the integration variables and
% extra-dependencies-variables of integr. factors need not
% be disjunct therefore the intfacdep-update only for
% not-potflag
>> else <<
if pri then write"BBB"$
% inhomogeneous part can not be integrated therefore
% reversing the succesful integration of the hom. part
if car intfaclist neq 1 then
qih:=reval list('QUOTIENT,qih,car intfaclist);
intfaclist:=cdr intfaclist
>>;
>>; %-- end of the integration of the inhomog. part
if pri then write"n=",n," nmax=",nmax," not pih=",not pih$
>> until (n=nmax) or (not pih)$ %---- end of v-integration
%------- variables of done integrations are collected for
%------- determining integrating factors in later integr.
if not zerop n then doneintvar:=cons(v,doneintvar)$
nges:=nges+n;
intlist:=cons(intfaclist,intlist)
>>$ % of not ( vari and (not member(v,vl)))
if vari then <<vari:=nil;vlcop:=nil>>;
if pri then write"ende: intvar=",intvar," vari=",vari,
" qh=",qh," qih=",qih$
>> % of (nmax neq zero)
>>$ % of ( while (vari or vlcop) and (pih or (not potflag)) )
% intlist and its elements intfaclist are in the inverse order
% to vl and the sequence of integrations done
q:=reval list('PLUS,qh,qih)$ %--- adding homog. and inhomog. part
if pri then <<terpri()$write"\\\ newcond:"$deprint newcond;
write"multiple result:",if null iflag then nil
else list(q,intlist,newcond,nges)
>>;
return if null iflag then nil
else list(q,intlist,newcond,nges)
end$ % of multipleint
symbolic procedure uplistoflds(intlist,listoflds)$
begin
scalar f,h1,h2,h3,h4,lds,itl;
while listoflds do
<<f:=caar listoflds;
lds:=cdar listoflds;
listoflds:=cdr listoflds;
h2:=nil; % h2 becomes the new list of lds of f
while lds do
<<h3:=car lds; lds:=cdr lds;
itl:=intlist;
h4:=nil; % h4 becomes one new ld of f
while h3 do
<<h4:=cons(car h3 - if null car itl then 0
else length car itl, h4);
h3:=cdr h3;itl:=cdr itl
>>;
h2:=cons(reverse h4,h2)
>>;
h1:=cons(cons(f,h2),h1)
>>;
return h1 % updated listoflds
end$ % of uplistoflds
symbolic procedure addintco(q, ftem, ifac, vl, vari)$
begin scalar v,f,l,vl1;
% multi.ing factors to the constants/functions of integration
if zerop q then l:=1
else
<<ftem:=fctsort ftem$
while ftem do
if fctlength car ftem<length vl then ftem:=nil
else if fctlinear(q,f) then
<<f:=car ftem$ ftem:=nil>> else
ftem:=cdr ftem$
if f then
<<l:=lderiv(q,f,fctargs f)$
l:=reval coeffn(q,reval car l,cdr l)
>> else l:=1
>>;
% the constants and functions of integration
if vari then q:=list('PLUS,q,intconst(l,vl,vari,list(1)))
else
<<vl1:=vl;
while vl1 do
<<v:=car vl1; vl1:=cdr vl1;
if car ifac then
q:=list('PLUS,q,intconst(l,vl,v,car ifac))$
% l..product of factors in the coefficient of the function to be
% eliminated, car ifac .. list of integrating factors
ifac:=cdr ifac;
>>
>>$
return reval q
end$ % of addintco
symbolic procedure integratepde(p,ftem,vari,genflag,potflag)$
% Generalized integration of the expression p
% if not genflag then "normal integration"
% Equation p must not be directly separable, otherwise the depen-
% dencies of functions of integration on their variables is wrong,
% i.e. no dependence on explicit variables
% ftem are all functions from the `global' ftem which occur in p, i.e.
% ftem:=smemberl(ftem,p)$
% if vari=nil then integration w.r.t. all possible variables within
% the equation
% else only w.r.t. vari one time
begin
scalar vl,vlrev,v,intlist,
ili1a,ili2a,maxvanz,fsub,h,hh,nfsub,iflag,newcond,
n1,n2,pot1,pot2,p1,p2,listoflds,secnd,ifac0,
ifac1a,ifac1b,ifac2a,ifac2b,cop,v1a,v2a,pri$
% pri:=t;
if pri then <<terpri()$write"Start Integratepde">>$
vl:=argset ftem$
vlrev:=reverse vl;
if vari then <<potflag:=nil;
if zerop p then iflag:='success>>
else
<<%------- determining fsub=list of functions of all variables
maxvanz:=length vl$
fsub:=nil;
h:=ftem;
while h do
<<if fctlength car h=maxvanz then
fsub:=cons(car h,fsub)$
h:=cdr h
>>$
nfsub:=length fsub$ % must be >1 for potential-integration
h:=intminderiv(p,ftem,vlrev,maxvanz,nfsub)$ % fsub is also for below
intlist:=car h$
%-- list of necessary integrations of the whole equation to solve
%-- for a function of all variables
listoflds:=cadr h$ %-- list of leading derivatives
>>$
if pri then <<terpri()$
write"complete integrations:",intlist," for:",vl>>;
%-- n1 is the number of integrations which must be done to try
%-- potential integration which must enable to eliminate 2 functions
%-- n2 is the number of integrations actually done
n1:=for each h in intlist sum h;
if (not vari) and (zerop n1) then
<<n2:=0;
if potflag then % else not necessary
for h:=1:(length vl) do ifac0:=cons(nil,ifac0)
>> else
<<if tr_genint then
<<terpri()$write "integration of the expression : "$
eqprint p>>$
if pri then
<<terpri()$write"at first all multiple complete integration">>;
%-- At first if possible n2 integrations of the whole equation
h:=multipleint(intlist,ftem,p,vari,vl,genflag,nil,nil,nil)$
% potflag=nil, partial=nil, doneintvar=nil
if h then
<<p:=car h;
ifac0:=cadr h; % ifac0 is the list of lists of integr. factors
newcond:=caddr h;
n2:=cadddr h; % number of done integrations
% doneintvar & intfacdep for the two halfs of integrations
% from the two parts of ifac0
h:=nil;
iflag:='success;
>> else n2:=0;
ftem:=union(fnew_,ftem)$
>>;
%------------ Existence of a potential ?
if (n1=n2) and potflag and (nfsub>1) then
%---- at least 2 functions to solve for
<<if not zerop n2 then %---- update listoflds
listoflds:=uplistoflds(reverse ifac0,listoflds)$
if pri then <<terpri()$write"uplistoflds:",listoflds>>$
if h:=potintegrable(listoflds) then
<<ili1a:=car h; ili2a:=cadr h;
% The necess. differentiations of the potential
if pri then
<<terpri()$write"potintegrable:",ili1a," ",ili2a>>$
if pri then <<write"+++ intlist=",intlist,
" ili1a=",ili1a,
" ili2a=",ili2a>>$
%-- distributing the integrating factors of ifac0 among
%-- the two lists ifac1b and ifac2b which are so far nil
%-- such that (ifac1b and ili1a are disjunct) and
%-- (ifac2b and ili2a are disjunct)
v1a:=vlofintlist(vl,ili1a);
v2a:=vlofintlist(vl,ili2a);
hh:=t;
cop:=reverse ifac0;
ifac1a:=ili1a;ifac2a:=ili2a;
while hh and cop do <<
% cop is a list of lists of integr. factors
if car cop then h:=vlofintfaclist(vl,cdar cop)
else h:=nil;
if freeoflist(h,v2a) and (car ifac2a=0) then <<
ifac1b:=cons( nil, ifac1b);
ifac2b:=cons( reverse car cop, ifac2b)
>> else
if freeoflist(h,v1a) and (car ifac1a=0) then <<
ifac2b:=cons( nil, ifac2b);
ifac1b:=cons( reverse car cop, ifac1b)
>> else
if car cop then hh:=nil;
ifac1a:=cdr ifac1a;
ifac2a:=cdr ifac2a;
cop:=cdr cop;
>>;
% the elements of ifac1b,ifac2b are in reverse order to
% ifac1a,ifac2a and are in the same order as vl, also
% the elements in the infac-lists are in inverse order,
% i.e. in the order the integrations have been done
if pri then <<terpri()$
write "ifac1a=",ifac1a," ifac1b=",ifac1b,
" ifac2a=",ifac2a," ifac2b=",ifac2b >>$
%-- lists of integrations to be done to both parts
if hh then
repeat % possibly a second try with part2 integrated first
<<n1:=for each n1 in ili1a sum n1;
% n1 .. number of remaining integrations of the first half
p1:=multipleint(ili1a,ftem,p,nil,vl,genflag,t,t,
% potflag=t, partial=t
union(vlofintlist(vl,ili2a),
vlofintlist(vl,ifac1b)))$
% that the variables of integration are not in ifac1b
% was already checked. Only restriction: the integrating
% factors must not depend on the last argument.
ftem:=union(fnew_,ftem)$
if p1 then <<
ifac1a:=cadr p1;
% ifac1a is now the list of integrating factors
if newcond then newcond:=nconc(newcond,caddr p1)
else newcond:=caddr p1;
if pri then <<terpri()$write"mul2: newcond=",newcond>>$
n2:=cadddr p1;
p1:=car p1
>>;
if p1 and (n1=n2) then
%--- if the first half has been integrated suff. often
<<%--- integrating the second half sufficiently often
n1:=for each n1 in ili2a sum n1;
% calculation of the 2. part which is not contained in p1
p2:=p1;
cop:=ifac1a; hh:=vlrev; % because ifac1a is reversed
while cop do <<
h:=car cop;cop:=cdr cop;
v:=car hh;hh:=cdr hh;
% h is the list of integrating factors of the v-integr.
while h do <<
p2:=reval list('QUOTIENT,list('DF,p2,v),car h);
h:=cdr h
>>
>>;
p2:=reval reval list('PLUS,p,list('MINUS,p2));
p2:=multipleint(ili2a,ftem,p2,nil,vl,genflag,t,nil,
% potflag=t, partial=nil
union(vlofintlist(vl,ili1a),
vlofintlist(vl,ifac2b)))$
ftem:=union(fnew_,ftem)$
if p2 then <<
ifac2a:=cadr p2;
% ifac2a is now list of integrating factors
if newcond then newcond:=nconc(newcond,caddr p2)
else newcond:=caddr p2;
if pri then <<terpri()$write"mul3: newcond=",newcond>>$
n2:=cadddr p2;
p2:=car p2
>>;
if p2 and (n1=n2) then
% if the second half has been integrated sufficiently often
<<% can both halfes be solved for different functions
% i.e. are they algebraic and linear in different functions?
pot1:=nil;
pot2:=nil;
for each h in fsub do <<
if ld_deriv_search(p1,h,vl) = (nil . 1) then
pot1:=cons(h,pot1);
if ld_deriv_search(p2,h,vl) = (nil . 1) then
pot2:=cons(h,pot2);
>>$
if (null not_included(pot1,pot2)) or
(null not_included(pot2,pot1)) then p2:=nil
>>$
if p2 and (n1=n2) then
<<iflag:='potint;
% ifac1b,ifac2b are in reverse order to ifac1a,ifac2a!
pot1:=newfct(fname_,vl,nfct_)$ % the new potential fct.
pot2:=pot1;
nfct_:=add1 nfct_$
fnew_:=cons(pot1,fnew_);
v:=vlrev;
while v do
<<cop:=car ifac1a; ifac1a:=cdr ifac1a;
while cop do <<
pot1:=reval list('QUOTIENT,list('DF,pot1,car v),car cop);
cop:=cdr cop
>>;
cop:=car ifac2a; ifac2a:=cdr ifac2a;
while cop do <<
pot2:=reval list('QUOTIENT,list('DF,pot2,car v),car cop);
cop:=cdr cop
>>;
v:=cdr v;
>>;
p:=addintco(list('PLUS,p1,reval pot2),
ftem,ifac1b,vlrev,nil)$
newcond:=cons(addintco(list('PLUS,p2,
list('MINUS,reval pot1)),
ftem,ifac2b,vlrev,nil),
newcond) % vari=nil
>>
;if pri then write":::"$
>>;
secnd:=not secnd;
% retry in different order of integration, p is still the same
if (iflag neq 'potint) and secnd then
<<cop:=ili1a;ili1a:=ili2a;ili2a:=cop>>
>> until (iflag eq 'potint) or (not secnd)
>>;
>>$
%--------- returning the result
return if not iflag then nil
else
<<if iflag neq 'potint then % constants of integration
p:=addintco(p, ftem, % the completely reversed ifac0
<<h:=nil;
while ifac0 do <<h:=cons(reverse car ifac0,h);ifac0:=cdr ifac0>>;
h
>>, vl, vari)$
if pri then
<<terpri()$write"ENDE INTEGRATEPDE"$deprint(cons(p,newcond))>>$
cons(p,newcond)
>>
end$ % of integratepde
symbolic procedure intpde(p,ftem,vl,x,potint)$
% integration of an polynomial expr. p w.r.t. x
begin scalar f,ft,l,l1,l2,vl,k,s,a,iflag,flag$
ft:=ftem$
vl:=cons(x,delete(x,vl))$
while ftem and not flag do
<<f:=car ftem$ % integrating all terms including car ftem
if member(x,fctargs f) then
<<l1:=l:=lderiv(p,f,vl)$
while not (flag or (iflag:=intlintest(l,x))) do
<<k:=reval coeffn(p,car l,cdr l)$
if intcoefftest(car lderiv(k,f,vl),car l,vl) then
<<a:=decderiv(car l,x)$
k:=reval list('INT,subst('v_a_r_,a,k),'v_a_r_)$
k:=reval subst(a,'v_a_r_,k)$
s:=cons(k,s)$
p:=reval aeval list('DIFFERENCE,p,list('DF,k,x))$
if (l:=lderiv(p,f,vl))=l1 then flag:='neverending
else l1:=l
>> else
flag:='coeffld
>>$
% iflag='nofct is the so far integrable case
% the non-integrable cases are: flag neq nil,
% (iflag neq nil) and (iflag neq 'nofct), an exception to the
% second case is potint where non-integrable rests are allowed
if iflag='nofct then ftem:=smemberl(ftem,p)
else if potint or (fctlength f<length vl) then
<<ftem:=cdr ftem$flag:=nil>> else
flag:=(iflag or flag)
>> else
ftem:=cdr ftem
>>$
return
if not flag then
<<l:=explicitpart(p,ft,x)$
l1:=list('INT,l,x)$
s:=reval aeval cons('PLUS,cons(l1,s))$
if freeint_ and
(null freeof(s,'INT)) then nil
else <<
k:=start_let_rules()$
l2 := reval {'DF,l1,x} where !*precise=nil;
if 0 neq (reval {'DIFFERENCE,l,l2} where !*precise=nil) then <<
write"REDUCE integrator error:"$terpri()$
algebraic write "int(",l,",",x,") neq ",l1;terpri()$
write"Result ignored.";terpri()$
stop_let_rules(k)$
nil
>> else <<
p:=reval reval aeval list('DIFFERENCE,p,l2)$
stop_let_rules(k)$
if poly_only then if ratexp(s,ft) then list(s,p)
else nil
else list(s,p)
>>
>>
>> else nil$
end$ % of intpde
symbolic procedure explicitpart(p,ft,x)$
% part of a sum p which only explicitly depends on a variable x
begin scalar l$
if not member(x,argset smemberl(ft,p)) then l:=p
else if pairp p then
<<if car p='MINUS then
l:=reval list('MINUS,explicitpart(cadr p,ft,x))$
if (car p='QUOTIENT) and not member(x,argset smemberl(ft,caddr p))
then
l:=reval list('QUOTIENT,explicitpart(cadr p,ft,x),caddr p)$
if car p='PLUS then
<<for each a in cdr p do
if not member(x,argset smemberl(ft,a)) then l:=cons(a,l)$
if pairp l then l:=reval cons('PLUS,l)>> >>$
if not l then l:=0$
return l$
end$
symbolic procedure intconst(co,vl,x,ifalist)$
% The factors in ifalist must be in the order of integrations done
if null ifalist then 0 else
begin scalar l,l2,f,coli,cotmp$
while ifalist do <<
cotmp:=coli;
coli:=nil;
while cotmp do <<
coli:=cons(list('INT,list('TIMES,car ifalist,car cotmp),x),coli);
cotmp:=cdr cotmp
>>;
coli:=cons(1,coli);
ifalist:=cdr ifalist
>>;
while coli do
<<f:=newfct(fname_,delete(x,vl),nfct_)$
nfct_:=add1 nfct_$
fnew_:=cons(f,fnew_)$
l:=cons(list('TIMES,f,car coli),l)$
coli:=cdr coli
>>$
if length l>1 then l:=cons('PLUS,l)
else if pairp l then l:=car l
else l:=0$
if co and co neq 1 then
if pairp co then
<<if car co='TIMES then co:=cdr co
else co:=list co$
while co do <<if my_freeof(car co,x) then l2:=cons(car co,l2)$
co:=cdr co>>
>>
else if co neq x then l2:=list co$
return reval if l2 then cons('TIMES,cons(l,l2))
else l
end$
symbolic procedure intcoefftest(l,l1,vl)$
begin scalar s$
return if not pairp l then t
else if car l='DF then
<<s:=decderiv(l1,car vl)$
if pairp s and pairp cdr s then s:=cddr s
else s:=nil$
if diffrelp(cons(cddr l,1),cons(s,1),vl) then t
else nil>>
else t$
end$
symbolic procedure fctlinear(p,f)$
<<p:=ldiffp(p,f)$
(null car p) and (cdr p=1)>>$
symbolic procedure intlintest(l,x)$
% Test,ob in einem Paar (fuehrende Ableitung.Potenz)
% eine x-Ableitung linear auftritt
if not car l then
if zerop cdr l then 'nofct
else 'nonlin
else % Fkt. tritt auf
if (car l) and (cdr l=1) then % fuehr. Abl. tritt linear auf
if pairp car l then % Abl. der Fkt. tritt auf
if caar l='DF then
if member(x,cddar l) then nil
% x-Abl. tritt auf
else if member(x,fctargs cadar l) then 'noxdrv
else 'noxdep
else 'nodrv
else 'nodrv
else 'nonlin$
symbolic procedure decderiv(l,x)$
% in Diff.ausdr. l wird Ordn. d. Abl. nach x um 1 gesenkt
begin scalar l1$
return if l then if car l='DF then
<<l1:=decderiv1(cddr l,x)$
if l1 then cons('DF,cons(cadr l,l1))
else cadr l>>
else l
else nil$
end$
symbolic procedure decderiv1(l,x)$
if null l then nil
else
if car l=x then
if cdr l then
if numberp cadr l then
if cadr l>2 then cons(car l,cons(sub1 cadr l,cddr l))
else cons(car l,cddr l)
else cdr l
else nil
else cons(car l,decderiv1(cdr l,x))$
symbolic procedure integratede(q,ftem,genflag)$
% Integration of a de
% result: newde if successfull, nil otherwise
begin scalar l,l1,l2,fl$
ftem:=smemberl(ftem,q)$
again:
if l1:=integrableode(q,ftem) then % looking for an integrable ode
if l1:=integrateode(q,car l1,cadr l1,caddr l1,ftem) then
% trying to integrate it
<<l:=append(cdr l1,l);
q:=simplifypde(car l1,ftem,nil)$
ftem:=smemberl(union(fnew_,ftem),q)$
fl:=t
>>$
if l1:=integratepde(q,ftem,nil,genflag,potint_) then
% trying to integrate a pde
<<q:=car l1$
for each a in cdr l1 do
<<ftem:=union(fnew_,ftem)$
if (l2:=integratede(a,ftem,nil)) then l:=append(l2,l)
else l:=cons(a,l)
>>$
fl:=t$
if null genflag then l1:=nil$
ftem:=smemberl(union(fnew_,ftem),q);
goto again
>>$
if fl then
<<l:=cons(q,l)$
l:=for each a in l collect reval aeval a$
l:=for each a in l collect
if pairp a and (car a='QUOTIENT) then cadr a
else a>>$
return l$
end$
symbolic procedure intflagtest(q,fullint)$
if flagp(q,'to_int) then
if fullint then
if get(q,'full_int) then t else
if get(q,'starde) then nil else
begin scalar fl,vl,dl,l,n,mi$
n:=get(q,'nvars)$
for each f in get(q,'rational) do % only rational fcts
if fctlength f=n then fl:=cons(f,fl)$
if null fl then return nil$
vl:=get(q,'vars)$
dl:=get(q,'derivs)$
for each d in dl do
if member(caar d,fl) then
put(caar d,'maxderivs,maxderivs(get(caar d,'maxderivs),cdar d,vl))$
dl:=fl$
while vl do
<<mi:=car get(car fl,'maxderivs)$
l:=list car fl$
put(car fl,'maxderivs,cdr get(car fl,'maxderivs))$
for each f in cdr fl do
<<if (n:=car get(f,'maxderivs))=mi then l:=cons(f,l)
else if n<mi then
<<l:=list f$
mi:=n>>$
put(f,'maxderivs,cdr get(f,'maxderivs))
>>$
dl:=intersection(l,dl)$
if dl then vl:=cdr vl
else vl:=nil>>$
for each f in fl do remprop(f,'maxderivs)$
return dl
end
else t$
symbolic procedure integrate(q,genintflag,fullint)$
% integrate pde q; if genintflag is not nil then indirect
% integration is allowed
% if fullint is not nil then only full integration is allowed
% Es wird noch nicht ausgenutzt:
% 1) Fcts, die rational auftreten
% 2) starde
begin scalar l$
if intflagtest(q,fullint) then
<<if (l:=integratede(get(q,'val),get(q,'fcts),genintflag)) then
<<ftem_:=reverse ftem_$
for each f in reverse fnew_ do
ftem_:=fctinsert(f,ftem_)$
ftem_:=reverse ftem_$
fnew_:=nil$
flag1(q,'to_eval)$
update(q,car l,ftem_,get(q,'vars),t,list(0))$
l:=cons(q,mkeqlist(cdr l,ftem_,get(q,'vars),
allflags_,t,get(q,'orderings)))$
remflag1(q,'to_int)$
if print_ then
<<terpri()$
if cdr l
then write "Indirect integration of ",q," yields ",l
else write "Integration of ",q>>$
>>$
remflag1(q,'to_int)>>$
return l$
end$
symbolic procedure quick_integrate_one_pde(pdes)$
begin scalar m,q,p$
% find the lowest order derivative wrt. only one variable
while pdes and
(get(car pdes,'length) = 1) do << % only 1 term
q:=caar get(car pdes,'derivs)$
if ( length q = 2 ) or
((length q = 3) and
(fixp caddr q) and
((null p) or
(caddr q<m) ) ) then
<<p:=car pdes;
m:=if (length q = 2) then 1
else caddr q>>;
pdes:=cdr pdes
>>$
if p then p:=integrate(p,nil,t)$
return p
end$
symbolic procedure integrate_one_pde(pdes,genintflag,fullint)$
% trying to integrate one pde
begin scalar l,l1,m,p$
% at first selecting all eligible de's
m:=-1$
while pdes do <<
if flagp(car pdes,'to_int) and not(get(car pdes,'starde)) then <<
l:=cons(car pdes,l);
if get(car pdes,'nvars)>m then m:=get(car pdes,'nvars)$
>>;
pdes:=cdr pdes
>>;
l:=reverse l;
% find an equation with as many as possible variables for integration
while m>=0 do <<
l1:=l$
while l1 do
if (get(car l1,'nvars)=m) and
(p:=integrate(car l1,genintflag,fullint)) then <<
m:=-1$
l1:=nil
>> else l1:=cdr l1$
m:=sub1 m
>>$
return p$
end$
endmodule$
%********************************************************************
module generalized_integration$
%********************************************************************
% Routines for generalized integration of pde's containing unknown
% functions
% Author: Andreas Brand
% December 1991
symbolic procedure gintorder(p,ftem,vl,x)$
% reorder equation p
begin scalar l,l1,q,m,b,c,q1,q2$
if pairp(p) and (car p='QUOTIENT) then <<
q:=caddr p$
p:=cadr p$
l1:=gintorder1(q,ftem,x,t)$
if DepOnAllVars(car l1,x,vl) then return nil;
q1:=car l1;
q2:=cadr l1;
>>$
if pairp(p) and (car p='PLUS) then p:=cdr p % list of summands
else p:=list p$
while p do
<<l1:=gintorder1(car p,ftem,x,nil)$
if DepOnAllVars(car l1,x,vl) then l:=p:=nil
else <<l:=termsort(l,l1)$p:=cdr p>> >>$
if l then
<<l:=for each a in l collect if cddr a then
<<b:=car a$
c:=cdr reval coeff1(cons('PLUS,cdr a),x,nil)$
m:=0$
while c and (car c=0) do <<c:=cdr c$m:=add1 m>>$
if m>0 then b:=list('TIMES,list('EXPT,x,m),b)$
cons(reval b,c)>>
else
cons(reval car a,cdr reval coeff1(cadr a,x,nil))$
if q then <<
l:=for each a in l collect
cons(car a,for each s in cdr a collect
reval list('QUOTIENT,s,q2))$
l:=for each a in l collect
cons(reval list('QUOTIENT,car a,q1),cdr a)
>>$
>>$
return l$
end$
symbolic procedure DepOnAllVars(c,x,vl)$
% tests for occurence of vars from vl in factors of c depending on x
begin scalar l$
if pairp c and (car c='TIMES) then c:=cdr c
else c:=list c$
while c and vl do
<<if not my_freeof(car c,x) then
for each v in vl do if not my_freeof(car c,v) then l:=cons(v,l)$
vl:=setdiff(vl,l)$
c:=cdr c
>>$
return (null vl)$
end$
symbolic procedure gintorder1(p,ftem,x,mode2)$
% reorder a term p
begin scalar l1,l2,sig$
% mode2 = nil then
% l2:list of factors of p not depending
% on x or beeing a power of x
% l1:all other factors
% mode2 = t then
% l2:list of factors of p not depending on x
% l1:all other factors
if pairp p and (car p='MINUS) then <<sig:=t$p:=cadr p>>$
if pairp p and (car p='TIMES) then p:=cdr p
else p:=list p$
for each a in p do
<<if my_freeof(a,x) and freeoflist(a,ftem) then l2:=cons(a,l2)
% freeoflist(a,ftem) to preserve linearity
else if mode2 then l1:=cons(a,l1)
else if a=x then l2:=cons(a,l2)
else if pairp a and (car a='EXPT) and (cadr a=x) and fixp caddr a
then l2:=cons(a,l2)
else l1:=cons(a,l1)>>$
if pairp l1 then
if cdr l1 then l1:=cons('TIMES,l1)
else l1:=car l1$
if pairp l2 then
if cdr l2 then l2:=cons('TIMES,l2)
else l2:=car l2$
if sig then if l2 then l2:=list('MINUS,l2)
else l2:=list('MINUS,1)$
return list(if l1 then l1 else 1,if l2 then l2 else 1)$
end$
symbolic procedure partint(p,ftem,vl,x,genint)$
begin scalar f,neg,l1,l2,n,k,l,h$
if tr_genint then <<
terpri()$
write "generalized integration of the unintegrated rest : "$
eqprint p
>>$
l:=gintorder(p,ftem,vl,x)$
% would too many new equations and functions be necessary?
if pairp(l) and (length(l)>genint) then return nil;
l:=for each s in l collect <<
h:=varslist(car s,ftem,vl)$
if h=nil then <<
list('TIMES,x,car s,cons('PLUS,cdr s))
>> else <<
f:=newfct(fname_,h,nfct_)$
nfct_:=add1 nfct_$
fnew_:=cons(f,fnew_)$
neg:=t$
n:=sub1 length cdr s$
k:=-1$
if (pairp car s) and (caar s='DF) then
<<h:=reval reval list('DIFFERENCE,cadar s,list('DF,f,x,add1 n))$
if not zerop h then l1:=cons(h,l1)$
l2:=cddar s>>
else
<<h:=signchange reval reval list('DIFFERENCE,car s,
list('DF,f,x,add1 n))$
if not zerop h then l1:=cons(h,l1)$
l2:=nil>>$
reval cons('PLUS, for each sk on cdr s collect
<<neg:=not neg$
k:=add1 k$
reval list('TIMES,if neg then -1 else 1,
append(list('DF,f,x,n-k),l2),
tailsum(sk,k,x) )
>>
)
>>
>>$
if l then l:=cons(reval cons('PLUS,l),l1)$
if tr_genint then
<<terpri()$
write "result (without constant or function of integration): "$
if l then <<
eqprint car l$
write "additional equations : "$
deprint cdr l
>> else write " nil "$
>>$
return l$
end$
symbolic procedure tailsum(sk,k,x)$
begin scalar j$
j:=-1$
return reval cons('PLUS,
for each a in sk collect
<<j:=j+1$
list('TIMES,a,prod(j+1,j+k),list('EXPT,x,j)) >> )
end$
symbolic procedure prod(m,n)$
if m>n then 1
else for i:=m:n product i$
endmodule$
%********************************************************************
module intfactor$
%********************************************************************
% Routines for finding integrating factors of PDEs
% Author: Thomas Wolf
% July 1994
% The following without factorization --> faster but less general
%symbolic procedure fctrs(p,vl,v)$
%begin scalar fl1,fl2,neg;
%
%write"p=",p;
%
% if car p='MINUS then <<neg:=t;p:=cdr p>>$
% return
% if not pairp p then if my_freeof(p,v) and (not freeoflist(p,vl)) then
% list(p,1,neg) else
% list(1,p,neg)
% else if car p='PLUS then list(1,p,neg)
% else
% if car p='TIMES then
% <<for each el in cdr p do if my_freeof(el,v) and (not freeoflist(p,vl)) then
% fl1:=cons(el,fl1) else
% fl2:=cons(el,fl2);
% if pairp fl1 then fl1:=cons('TIMES,fl1);
% if pairp fl2 then fl2:=cons('TIMES,fl2);
% if not fl1 then fl1:=1;
% if not fl2 then fl2:=1;
% list(fl1,fl2,neg)
% >> else if my_freeof(p,v) and (not freeoflist(p,vl)) then
% list(p,1,neg) else
% list(1,p,neg)
%end$ % of fctrs
%
symbolic procedure fctrs(p,indep,v)$
begin scalar fl1,fl2;
p:=cdr reval factorize p;
for each el in p do
if freeoflist(el,indep) and
((v=nil) or (not my_freeof(el,v))) then fl1:=cons(el,fl1)
else fl2:=cons(el,fl2);
if null fl1 then fl1:=1;
if null fl2 then fl2:=1;
if pairp fl1 then if length fl1 = 1 then fl1:=car fl1
else fl1:=cons('TIMES,fl1);
if pairp fl2 then if length fl2 = 1 then fl2:=car fl2
else fl2:=cons('TIMES,fl2);
return list(fl1,fl2)
end$ % of fctrs
symbolic procedure extractfac(p,indep,v)$
% looks for factors of p dependent of v and independent of indep
% and returns a list of the numerator factors and a list of the
% denominator factors
begin scalar nu,de$
return
if (pairp p) and (car p='QUOTIENT) then
<<nu:=fctrs( cadr p,indep,v);
de:=fctrs(caddr p,indep,v);
list( reval if car de neq 1 then list('QUOTIENT, car nu, car de)
else car nu,
if cadr de neq 1 then list('QUOTIENT, cadr nu, cadr de)
else cadr nu
)
>> else fctrs(p,indep,v)
end$ % of extractfac
%----------------------------
symbolic procedure get_kernels(ex)$
% gets the list of all kernels in ex
begin scalar res,pri$
% pri:=t;
ex:=reval ex$
if pri then <<terpri()$write"ex=",ex>>;
if pairp ex then
if (car ex='QUOTIENT) or (car ex='PLUS) or (car ex='TIMES) then
for each s in cdr ex do res:=union(get_kernels s,res) else
if (car ex='MINUS) or
((car ex='EXPT) and
% (numberp caddr ex)) then % not for e.g. (quotient,2,3)
(cadr ex neq 'E) and
(cadr ex neq 'e) and
(not fixp cadr ex) ) then res:=get_kernels cadr ex
else res:=list ex
else if idp ex then res:=list ex$
if pri then <<terpri()$write"res=",res>>;
return res$
end$
%------------------
symbolic procedure specialsol(p,vl,fl,x,indep,gk)$
% tries a power ansatz for the functions in fl in the kernels
% of p to make p to zero
% indep is a list of kernels, on which the special solution should
% not depend. Is useful, to reduce the search-space, e.g. when an
% integrating factor for a linear equation for, say f, is to be
% found then f itself can not turn up in the integrating factor fl
% gk are kernels which occur in p and possibly extra ones which
% e.g. are not longer in p because of factorizing p but which are
% likely to play a role, if nil then determined below
% x is a variable on which each factor in the special solution has
% to depend on.
begin
scalar e1,e2,n,nl,h,hh,ai,sublist,eqs,startval,pri,printold,pcopy;
%pri:=t;
p:=num p;
pcopy:=p;
if pri then <<
terpri()$write"The equation for the integrating factor:";
terpri()$eqprint p;
>>;
if null gk then gk:=get_kernels(p);
for each e1 in fl do <<
h:=nil; %---- h is the power ansatz
if pri then
for each e2 in gk do <<
terpri()$write"e2=",e2;
if my_freeof(e2,x) then write" freeof1";
if not freeoflist(e2,fl) then write" not freeoflist"$
if not freeoflist(e2,indep) then write" dependent on indep"
>>;
%----- nl is the list of constants to be found
for each e2 in gk do
if (not my_freeof(e2,x)) and % integ. fac. should depend on x
freeoflist(e2,fl) and % the ansatz for the functions to be
% solved should not include these functions
freeoflist(e2,indep) then <<
n:=gensym();nl:=cons(n,nl);
h:=cons(list('EXPT,e2,n),h);
>>;
if h then <<
if length h > 1 then h:=cons('TIMES,h)
else h:=car h;
%-- the list of substitutions for the special ansatz
sublist:=cons((e1 . h),sublist);
if pri then <<terpri()$write"Ansatz: ",e1," = ",h>>;
p:= reval num reval subst(h,e1,p);
if pri then <<terpri()$write"p=";eqprint p>>
>>
>>;
if null h then return nil;
%------- An numerical approach to solve for the constants
if nil then << % numerical approach
%--- Substituting all kernels in p by numbers
on rounded;
precision 20;
terpri()$terpri()$write"Before substituting numerical values:";
eqprint p;
terpri()$terpri()$write"Constants to be calculated: ";
for each n in nl do write n," ";
for each e1 in nl do <<
h:=p;
for each e2 in gk do
if freeoflist(e2,fl) then
if pairp e2 and ((car e2 = 'DF) or (car e2 = 'INT)) then <<
n:=list('QUOTIENT,1+random 30028,30029);
terpri();write"substitution done: ",e2," = ",n;
h:=subst(list('QUOTIENT,1+random 30028,30029),e2,h)
>>;
for each e2 in gk do
if freeoflist(e2,fl) then
if not(pairp e2 and ((car e2 = 'DF) or (car e2 = 'INT))) then <<
n:=list('QUOTIENT,1+random 30028,30029);
terpri();write"substitution done: ",e2," = ",n;
h:=subst(list('QUOTIENT,1+random 30028,30029),e2,h)
>>;
terpri()$terpri()$write"The equation after all substitutions: ";
terpri()$
eqprint h;
eqs:=cons(reval h,eqs);
startval:=cons(list('EQUAL,e1,1),startval)
>>;
if length eqs = 1 then eqs:=cdr eqs
else eqs:=cons('LIST,eqs);
if length startval = 1 then startval:=cdr startval
else startval:=cons('LIST,startval);
terpri()$write"start rdsolveeval!";terpri()$terpri()$
h:=rdsolveeval list(eqs,startval);
eqs:=nil;
off rounded;
>>;
%----- An analytical approach to solve for the constants
if null pri then <<printold:=print_;print_:=nil>>;
if p and not zerop p then % uebernommen aus SEPAR
if not (pairp p and (car p='QUOTIENT) and % " " "
intersection(argset smemberl(fl,cadr p),vl)) then
p:=separ2(p,fl,vl) else
p:=nil;
if null pri then print_:=printold;
if p then <<
% possibly investigating linear dependencies of different caar p
% solve(lasse x-abhaengige und nl-unabhaengige faktoren fallen von
% factorize df(reval list('QUOTIENT, caar p1, caar p2),x),nl).
while p do
if freeoflist(cdar p,nl) then <<eqs:=nil;p:=nil>>
% singular system --> no solution
else <<
eqs:=cons(cdar p,eqs);
p:=cdr p
>>;
>>;
if pri then <<terpri()$write"eqs1=",eqs>>;
if (null eqs) or (length eqs > maxalgsys_) then return nil
else <<
if pri then <<
terpri()$write"The algebraic system to solve for ",nl," is:";
if length eqs > 1 then deprint eqs
else eqprint car eqs
>>;
if length eqs > 1 then eqs:=cons('LIST,eqs)
else eqs:=car eqs;
if pri then <<terpri()$write"eqs2=",eqs;terpri();write"nl=",nl>>$
% for catching the error message `singular equations'
hh:=cons('LIST,nl);
eqs:=<<
ai:=!!arbint;
h:=errorset({'solveeval,mkquote{eqs, hh}},nil,nil)
where !*protfg=t;
erfg!*:=nil;
if errorp h then nil else cdar h % cdr for deleting 'LIST
>>;
if pri then <<terpri()$write"eqs3a=",eqs," ai=",ai," !!arbint=",
!!arbint;terpri()>>$
if not freeof(eqs,'ARBCOMPLEX) then <<
eqs:=reval car eqs;
for h:=(ai+1):!!arbint do
eqs:=subst(0,list('ARBCOMPLEX,h),eqs);
if pri then <<terpri()$write"eqs3b=",eqs;terpri()>>$
eqs:=<<
h:=errorset({'solveeval,mkquote{eqs, hh}},nil,nil)
where !*protfg=t;
erfg!*:=nil;
if errorp h then nil else cdar h % cdr for deleting 'LIST
>>;
>>;
if pri then <<terpri()$write"eqs3c=",eqs;terpri()>>$
%--- eqs is the list of solutions for the constant exponents of the
%--- integrating factor
if null eqs then return nil;
if length nl=1 then eqs:=list eqs;
if pri then <<write"nl=",nl," eqs4=",eqs;terpri()>>;
for each e1 in eqs do << % each e1 is a list of substitutions
if pri then <<terpri()$write"e2=",e1;terpri()>>$
if car e1='LIST then e1:=cdr e1;
if pri then <<terpri()$write"e3=",e1;terpri()>>$
for each e2 in e1 do <<
if pri then algebraic write"solution:",symbolic e2;
sublist:=subst(caddr e2,cadr e2,sublist)
>>;
if pri then <<terpri()$write"The sublist is:",sublist>>
>>;
>>;
if pri then <<terpri()$write"pcopy1=",pcopy;terpri()>>$
for each e1 in sublist do <<
pcopy:=subst(cdr e1,car e1,pcopy);
if pri then <<terpri()$write"e1=",e1;terpri();
write"pcopy2=",pcopy;terpri()>>$
>>$
if pri then <<terpri()$write"pcopy3=",reval pcopy;terpri()>>$
if pri then <<terpri()$write"pcopy4=",reval reval pcopy;terpri()>>$
if not zerop reval reval pcopy then return nil else
return for each e1 in sublist collect (car e1 . reval cdr e1)
end$ % of specialsol
%------------------
symbolic operator findintfac$
symbolic procedure findintfac(pl,ftem,vl,x,doneintvar,intfacdep,
factr,verbse)$
% - pl is a list of equations from which the *-part (inhomogeneous
% terms) have been dropped.
% - each equation of pl gets an integrating factor h
% - doneintvar is a list of variables, on which the integrating factor
% should not depend. The chances to find an integrating factor
% increase if the inhomogeneous part of pl is dropped and
% separately integrated with generalized integration later.
% - if factr is not nil then the equation(s) pl is(are) at first
% factorized, e.g. if integration(s) have already been done
% and there is a chance that the equation can factorize, thereby
% simplify and giving a higher chance for integrability.
begin
scalar h,newequ,tozero,fl,e1,pri,factr,exfactors,ftr,gk;
% exfactors is the list of factors extracted at the beginning
% pri:=t;
factr:=t; % whether tozero should be factorized
if pri then <<terpri()$write"START VON FINDINTFAC">>;
%--- Generation of the condition for the integrating factor(s) in fl
for each e1 in pl do <<
%--- extracting factors dependend on x and independent of
%--- doneintvar but only if integrations have already been done,
%--- i.e. (doneintvar neq nil)
gk:=union(gk,get_kernels(e1));
if factr then <<ftr:=extractfac(e1,append(doneintvar,ftem),x);
if not evalnumberp car ftr then
gk:=union(gk,list car ftr);
>>
else ftr:=list(1,nil);
exfactors:=cons(car ftr,exfactors);
if car ftr neq 1 then <<
e1:=cadr ftr;
if pri then <<terpri()$write"extracted factor:",eqprint car ftr>>;
>>;
%--- fl is to become the list of integrating factors h
h:=gensym();
depl!*:=cons(list(h,x),depl!*)$
depend h,x;
fl:=h . fl;
e1:=intpde(reval list('TIMES,h,e1),ftem,vl,x,t);
if e1 and car e1 then <<
newequ:=car e1 . newequ;
tozero:=cadr e1 . tozero;
if pri then <<
terpri()$write" the main part of integration:"$ eqprint(car e1);
terpri()$write"car e1=",car e1;
terpri()$write" the remainder of integration:"$ eqprint(cadr e1);
terpri()$write"cadr e1=",cadr e1;
>>
>>;
>>;
if null tozero then return nil;
%-------- newequ is the integral
newequ:=if length pl > 1 then cons('PLUS,newequ)
else car newequ;
%-------- tozero is the PDE for the integrating factor
tozero:=reval if length pl > 1 then cons('PLUS,tozero)
else car tozero;
if pairp tozero and (car tozero='QUOTIENT) then tozero:=cadr tozero$
if factr then <<
h:=cdr reval list('FACTORIZE,tozero)$
if pri then <<terpri()$write"The factors of tozero:",h>>;
tozero:=nil;
for each e1 in h do
if smemberl(fl,e1) then tozero:=cons(e1,tozero)$
tozero:= reval if length tozero > 1 then cons('TIMES,tozero)
else car tozero;
>>;
if nil and pri then <<write"tozero =";eqprint tozero >>;
h:=nil;
% actually only those f in ftem, in which pl is nonlinear, but also
% then only integrating factors with a leading derivative low enough
h:=specialsol(tozero,vl,fl,x,append(ftem,doneintvar),gk);
% h:=specialsol(tozero,vl,fl,x,doneintvar,gk);
if pri then <<write"h=",h;terpri()>>;
if h then <<
for each e1 in h do << % each e1 is one integrating factor determined
if pri then <<terpri()$write"e1=",e1;
terpri()$write"newequ=",newequ;terpri()>>;
newequ:=reval subst(cdr e1,car e1,newequ);
if pri then <<terpri()$write"newequ=",newequ>>;
>>
>> else if pri then write"no integrating factor";
%--- delete all dependencies of the functions in fl
%--- must come before the following update
for each e1 in fl do <<
depl!*:=delete(assoc(e1,depl!*),depl!*)$
depl!*:=delete(assoc(mkid(e1,'_),depl!*),depl!*)$
>>;
%--- update intfacdep
for each e1 in vl do
if (e1 neq x) and my_freeof(intfacdep,e1) and
((not my_freeof(h,e1)) or (not my_freeof(exfactors,e1)))
then intfacdep:=cons(e1,intfacdep);
%--- returns nil if no integrating factor else a list of the
%--- factors and the integral
if h and print_ and verbse then <<
terpri()$write"The integrated equation:";
eqprint newequ;
terpri()$
if length pl = 1 then write"An integrating factor has been found:"
else write"Integrating factors have been found: "$
>>;
return if (null h) or (zerop newequ) then nil else
list(newequ,
for each e1 in h collect <<
ftr:=car exfactors;
exfactors:=cdr exfactors;
gk:=if ftr=1 then cdr e1
else reval list('QUOTIENT,cdr e1,ftr);
if print_ and verbse then mathprint gk;
gk
>>,
intfacdep)
end$
endmodule$
%********************************************************************
module odeintegration$
%********************************************************************
% Routines for integration of ode's containing unnown functions
% Author: Thomas Wolf
% August 1991
symbolic procedure integrateode(de,fold,xnew,ordr,ftem)$
% once equations are factorized before integration the % * lines
% can be droped or reduced
begin scalar newde,newnewde,l,l1,h,factrs,fc,changd,newcond,facnum$
if pairp de and (car de='QUOTIENT) then de:=cadr de$ % *
factrs:=cdr reval list('FACTORIZE,de); % *
facnum:=length factrs; % *
l:=for each fc in factrs collect % *
if not smember(fold,fc) then <<facnum:=facnum-1;fc>> % *
else % *
<<if facnum>1 then <<l1:=integrableode(fc,ftem); % *
if l1 then <<fold:=car l1; % *
xnew:=cadr l1; % *
ordr:=caddr l1 % *
>> % *
>> % *
else l1:=t;
h:= % the integrated equation
if not l1 then nil else % *
if not xnew then << % Integr. einer alg. Gl. fuer eine Abl.
newde:=cadr solveeval list(fc,fold)$
if not freeof(newde,'ROOT_OF) then nil
else <<
newde:=reval list('PLUS,cadr newde,list('MINUS,caddr newde))$
if (l:=integratepde(newde,ftem,nil,genint_,nil)) then
<<newcond:=append(newcond,cdr l);car l>>
%genflag=t,potflag=nil
else nil
>>
>> else % eine ode fuer ein f?
if not pairp fold then % i.e. not df(...,...), i.e. fold=f
odeconvert(fc,fold,xnew,ordr,ftem)
% --> ode fuer eine Abl. von f
else <<
newde:=odeconvert(fc,fold,xnew,ordr,ftem)$
if not newde then nil
else <<
newnewde:=cadr solveeval list(newde,fold)$
newnewde:=reval list('PLUS,cadr newnewde,list('MINUS,
caddr newnewde))$
ftem:=union(fnew_,ftem)$
newnewde:=integratede(newnewde,ftem,nil)$
if newnewde then <<newcond:=append(newcond,cdr newnewde);
car newnewde>>
else newde
>>
>>;
if h then <<changd:=t;h>> % factors to be collected % *
else fc % *
>>;
return if not changd then nil
else
cons(if length l > 1 then cons('TIMES,l) % *
else car l ,newcond) % *
end$ % of integrateode
symbolic procedure odecheck(ex,fint,ftem)$
% assumes an revaled expression ex
% Does wrong if car ex is a list!
begin scalar a,b,op,ex1$
%***** ex is a ftem-function *****
if ex=fint then % list(ex,0,0,..)
<<a:=list ex$
ex:=fctargs ex$
while ex do
<<a:=append(list(0,0),a)$
ex:=cdr ex>>$
% not checked if it is a function of an expression of x
op:=reverse a>>
else if pairp ex then
%***** car ex is 'df *****
if (car ex)='df then
<<a:=odecheck(cadr ex,fint,ftem)$
if not pairp a then op:=a
else % a is list(fctname,0,0,..,0,0)
<<op:=list(car a)$
a:=fctargs car a$ % a is list(variables), not checked
ex:=cddr ex$ % ex is list(derivatives)
while a do
<<ex1:=ex$
while ex1 and ((car a) neq (car ex1)) do ex1:=cdr ex1$
if null ex1 then op:=cons(0,cons(0,op))
else
<<if not cdr ex1 then b:=1 % b is number of deriv. of that var.
else
<<b:=cadr ex1$
if not numberp b then b:=1>>$
op:=cons(b,cons(b,op))>>$
a:=cdr a>>$
op:=reverse op>> >>
else
%***** car ex is a standard or other function *****
<<a:=car ex$ % for linearity check
ex:=cdr ex$
if a='INT then ex:=list reval car ex$
while (op neq '!_abb) and ex do
<<b:=odecheck(car ex,fint,ftem)$
if b then % function found
if b eq '!_abb then op:='!_abb % occures properly
else op:=odetest(op,b)$
ex:=cdr ex>> >>$
return op
end$
symbolic procedure integrableode(p,ftem)$
if delength p>(if odesolve_ then odesolve_ else 0) then
(if cont_ then
if yesp("expression to be integrated ? ") then
integrableode1(p,ftem))
else integrableode1(p,ftem)$
symbolic procedure integrableode1(p,ftem)$
begin scalar a,b,u,vl,le,uvar,
fint,fivar,% the function to be integrated and its variables
fold, % the new function of the ode
xnew, % the independ. variable of the ode
ordr1, % order of the ode
ordr2, % order of the derivative for which it is an ode
intlist$ % list of ode's
ftem:=smemberl(ftem,p)$
vl:=argset ftem$
% p muss genau eine Funktion aus ftem von allen Variablen enthalten.
% Die Integrationsvariable darf nicht Argument anderer in p enthaltener
% ftem-Funktionen sein.
a:=ftem$
b:=nil$
le:=length vl$
while a and vl do
<<u:=car a$
uvar:=fctargs u$
if (length uvar) = le then
if b then
<<vl:=nil$a:=list(nil)>>
else
<<b:=t$
fint:=u$
fivar:=uvar>>
else vl:=setdiff(vl,uvar)$
a:=cdr a>>$
if not b then vl:=nil$
le:=length p$
if ((1<le) and vl) then
<<a:=odecheck(p,fint,ftem)$
if not atom a then % The equation is an ode
<<ordr1:=0$
ordr2:=0$
xnew:=nil$
a:=cdr a$
b:=fivar$
while b do
<<if (car a) neq 0 then
<<fold:=cons(car b,fold)$
if (car a) > 1 then fold:=cons(car a,fold)>>$
ordr2:=ordr2+car a$
if (car a) neq (cadr a) then
<<xnew:=car b$
if not member(xnew,vl) then
<<b:=list(nil)$vl:=nil>>$
ordr1:=(cadr a) - (car a)>>$
b:=cdr b$
a:=cddr a>>$
fold:=reverse fold$
%fold is the list of diff.variables + number of diff.
if fold then fold:=cons('df,cons(fint,fold))
else fold:=fint$
if vl and ((ordr1 neq 0) or (ordr2 neq 0)) then
intlist:=list(fold,xnew,ordr1,ordr2)
>> % of variable found
>>$ % of if
return intlist
end$ % of integrableode1
symbolic procedure odetest(op,b)$
if not op then b
else % op=nil --> first function found
if (car op) neq (car b) then '!_abb else % f occurs in differ. fct.s
begin scalar dif,a$
dif:=nil$ % dif=t --> different derivatives
a:=list(car op)$ % in one variable already found
op:=cdr op$
b:=cdr b$
while op do
<<a:=cons(max(cadr op,cadr b),cons(min(car op,car b),a))$
if (car a) neq ( cadr a) then
if dif then
<<a:='!_abb$
op:=list(1,1)>>
else dif:=t$
op:=cddr op$
b:=cddr b>>$
if not atom a then a:=reverse a$
return a % i.e. '!_abb or (fctname,min x1-der.,max x1-der.,...)
end$
symbolic procedure odeconvert(de,ford,xnew,ordr,ftem)$
begin scalar j,ford_,newco,oldde,newde,newvl,null_,ruli,zd$
% trig1,trig2,trig3,trig4,trig5,trig6$
ford_:=gensym()$
depl!*:=delete(assoc(ford_,depl!*),depl!*)$
depend1(ford_,xnew,t)$
oldde:=reval subst(ford_,reval ford,de)$
if pairp ford then % i.e. (car ford) eq 'DF then
<< for j:=1:ordr do
oldde:= subst( reval list('DF,ford_,xnew,j),
reval list('DF,ford,xnew,j), oldde)>>$
algebraic !!arbconst:=0$
newde:=algebraic first
odesolve(symbolic oldde,symbolic ford_,symbolic xnew)$
ruli:= start_let_rules()$
newde:=reval newde;
% Instead of the following test one should return several cases
zd:=zero_den(newde,cons(ford_,ftem),union(list xnew,argset ftem));
% if safeint_ and zero_den(newde,ftem,argset ftem) then newde:=nil;
if freeint_ and null freeof(newde,'INT) then newde:=nil;
if newde and (cadr newde neq oldde) then << % solution found
% Test der Loesung klappt nur, wenn Loesung explizit gegeben
if cadr newde neq ford_ then <<
if print_ then
<<write "Solution of the ode is not explicitly given."$
algebraic write "Equation is: ",algebraic symbolic oldde$
algebraic write "Solution is: ",algebraic symbolic newde
>>;
if poly_only then % The solution must be rational in the
% function and constants of integration
if not rationalp(newde,ford_) then newde:=nil else <<
j:=1;
while (j leq ordr) and
rationalp(subst(ford_,list('arbconst,j),newde),ford_) do j:=j+1;
if j leq ordr then newde:=nil
>>;
if newde then
if (caadr newde = 'QUOTIENT) and (zerop caddr newde) then
newde:={'EQUAL,cadadr newde,0} else
if (caaddr newde = 'QUOTIENT) and (zerop cadr newde) then
newde:={'EQUAL,0,cadr caddr newde}
>> else <<
null_:=reval reval aeval subst(caddr newde, ford_, oldde)$
% reval reval because of a REDUCE bug for special data,
% to be dropped as soon as possible
if (null_ neq 0) then <<
% newde:=nil$
if print_ then <<
write "odesolve solves : "$
deprint list oldde$
write "to"$
deprint list newde$
Write "which inserted in the equation yields"$
deprint list null_$
>>
>>
>>
>>$
if newde then
<<newde:=list('PLUS,cadr newde,list('MINUS,caddr newde))$
if zerop reval list('PLUS,newde,list('MINUS,oldde)) then newde:=nil$
if newde and (zd neq nil) then
newde:=cons('TIMES,append(zd,list newde))$
>>$
depl!*:=delete(assoc(ford_,depl!*),depl!*)$
stop_let_rules(ruli)$
return
if null newde then nil
else
<<newde:=subst(ford,ford_,newde)$
newvl:=delete(xnew,if (pairp ford and (car ford='DF))
then fctargs cadr ford
else fctargs ford)$
% if pairp ford then newvl:=delete(xnew,cdr assoc(cadr ford,depl!*))
% else newvl:=delete(xnew,cdr assoc(ford,depl!*))$
for j:=1:ordr do <<
newco:=newfct(fname_,newvl,nfct_)$
nfct_:=add1 nfct_$
fnew_:=cons(newco,fnew_)$
newde:=subst(newco,list('arbconst,j),newde)
% newde:=subst(newco, prepf !*kk2f list('arbconst,j),newde)
% newde:=reval subst(newco,list('arbconst,j),newde)
% newde:=reval subst(newco, prepf !*kk2f list('arbconst,j),newde)
>>$
newde>>
end$
endmodule$
end$