Artifact 0b4db57d7c17c5159a9ca32375f338752d87fca1d1a49d0848f91e8fcb7848c6:
- Executable file
r38/packages/crack/crgensep.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 47622) [annotate] [blame] [check-ins using] [more...]
%********************************************************************* module gensep_lin$ %********************************************************************* % Routines for generalized separation of de's % Author: Andreas Brand, Thomas Wolf 1990 1994 1997 % Thomas Wolf since 1997 symbolic procedure quick_gen_separation(arglist)$ % Indirect separation of a pde if vl_ then % otherwise not possible --> save time begin scalar p,l,l1,pdes,stp$ % pdes:=clean_up(car arglist)$ % if pdes then l:=list(pdes,cadr arglist)$ % because the bookeeping of to_drop is not complete instead: pdes:=car arglist$ % to return the new list of pdes in case gensep is not successful if expert_mode then << l1:=selectpdes(pdes,1); if get(car l1,'starde) then flag(l1,'to_gensep) >> else l1:=cadddr arglist$ if (p:=get_gen_separ_pde(l1,t,t)) then % high priority <<l:=gensep(p,pdes)$ pdes:=delete(p,pdes)$ for each a in l do << pdes:=eqinsert(a,pdes)$ if member(a,pdes) and (stp:=get(a,'starde)) then to_do_list:=cons(list(if cdr stp=0 then 'separation else 'gen_separation, %pdes,cadr arglist,caddr arglist, list a), to_do_list) >>$ l:=list(pdes,cadr arglist)>>$ return l$ end$ symbolic procedure gen_separation(arglist)$ % Indirect separation of a pde if vl_ then % otherwise not possible --> save time begin scalar p,l,l1,pdes,stp$ % pdes:=clean_up(car arglist)$ % if pdes then l:=list(pdes,cadr arglist)$ % because the bookeeping of to_drop is not complete instead: pdes:=car arglist$ % to return the new list of pdes in case gensep is not successful if expert_mode then << l1:=selectpdes(pdes,1); if get(car l1,'starde) then flag(l1,'to_gensep) >> else l1:=cadddr arglist$ if (p:=get_gen_separ_pde(l1,nil,t)) then % low priority <<l:=gensep(p,pdes)$ pdes:=delete(p,pdes)$ for each a in l do << pdes:=eqinsert(a,pdes)$ if member(a,pdes) and (stp:=get(a,'starde)) then to_do_list:=cons(list(if cdr stp=0 then 'separation else 'gen_separation, %pdes,cadr arglist,caddr arglist, list a), to_do_list) >>$ l:=list(pdes,cadr arglist)>>$ return l$ end$ symbolic procedure maxnoargs(fl,v)$ % determines the maximal number of arguments of any of the % functions of fl begin scalar f,n,m; n:=0; for each f in fl do <<m:=fctargs f; m:=if m and not_included(v,m) then length m else 0; if n<m then n:=m >>; return n end$ symbolic procedure get_gen_separ_pde(pdes,high_priority,lin)$ % looking for a pde in pdes which can be indirectly separated % if lin then only a liner PDE % p ...the next equation that will be chosen % dw...whether p was already delt with % na...number of variables in the equation % nv...maximal number of arguments of any of the functions of p % nf...min number of functions to be dropped before direct sep. % leng...length of p begin scalar p,nv,nf,dw,len,h1,h2,h3,h4,nvmax$ %na,h5 % ncmax:=nvmax$ if high_priority then << nvmax:=0; for each p in pdes do if (h1:=get(p,'nvars))>nvmax then nvmax:=h1; p:=nil >>$ while pdes do << if flagp(car pdes,'to_gensep) and (null lin or get(car pdes,'linear_)) and % not too many terms or enough terms <<h1:=get(car pdes,'length)$ (null high_priority) or (get(car pdes,'nvars)=nvmax) or ( low_gensep>h1) or (high_gensep<h1) >> and % no single function depending on all variables: (h3:=get(car pdes,'starde) ) and % not directly separable: (cdr h3 neq 0 ) and % Each time the equation is investigated and differentiated % wrt the first element of car h3, this element is dropped. % --> The equation should not already have been differentiated % wrt all variables: (not null car h3 ) and % If equations have been investigated by generalized % separation or if equations resulted from generalized % separation then they get the flag used_ to be solved % first, not to have too many unevaluated new functions % at a time ((h4:=flagp(car pdes,'used_) ) or (null dw) ) and % The variables in h3 are the ones wrt which direct separation % shall be achieved after differentiation, therefore functions % of these variables have to be thrown out. The remaining % functions shall be of as many as possible arguments to % make quick progress: ((null p ) or (len > h1 ) or % neu ((len = h1) and ( % neu (nv < (h2:=maxnoargs( get(car pdes,'fcts), car h3 )) ) or ((nv = h2) and ( % (na < (h5:=get(car pdes,'nvars)) ) or % ((na = h5) and ( ((null dw) and flagp(car pdes,'used_)) or (( nf > cdr h3 ) or ((nf = cdr h3 ) and (len > h1 ) ) ) ) )))) then <<p:=car pdes; nv:=if (null nv) or (null h2) then maxnoargs(get(p,'fcts),car h3) else h2; % na:=if (null na) or (null h5) then get(car pdes,'nvars) % else h5; if h4 then dw:=h4; nf:=cdr h3; len:=h1>>$ pdes:=cdr pdes$ >>; return p end$ %----------------- symbolic procedure gensep(p,pdes)$ % generalized separation of pde p if zerop cdr get(p,'starde) then separate(p,pdes) % be dropped? else % e.g. a=((x y z).2) % POSSIBLE REASONS FOR FAILURE: % - After the sequence of divisions and differentiations in the direct % separation, if there explicit v-dependent coefficients are taken % out which also contain later integration variables, then the integrands % are not total derivatives anymore --> no backintegration is possible. % - This corresponds to the case when all variables occur explicitly but % in a non-product expression, like sin(x*y*z) begin scalar a,pl$ if print_ then <<terpri()$write "generalized separation of ",p>>$ if tr_gensep then <<a:=get(p,'starde)$ terpri()$write "de to be separated : "$ typeeqlist list p$ terpri()$write "variable list for separation : ",car a$ terpri()$write "for each of these variables ",cdr a$ if (cdr a)=1 then write" function does" else write" functions"$ write" depend on it " >>$ %--- so far only one DE p in the pool starlist of equations pl:=partitn(get(p,'val), % expression nil, % history of divisions/diff so far get(p,'fcts), % functions get(p,'vars), % variables car get(p,'starde) % separation charac. ); if pl then << pl:=append(for each a in car pl collect cdr a,cadr pl); pl:=mkeqlist(pl,fctsort union(ftem_,get(p,'fcts)),get(p,'vars), cons('to_drop,allflags_),t,get(p,'orderings),pdes)$ drop_pde(p,nil,nil); flag(pl,'used_); if print_ then <<terpri()$ a:=length pl$ write "separation yields ",a," new equation"$ if a > 1 then write"s"$ write" : "$ if tr_gensep then typeeqlist pl else listprint pl$ terpri() >> >> else << remflag1(p,'to_gensep)$ pl:=list p >>$ return pl$ end$ %----------------- symbolic procedure partitn(q,old_histy,ftem,vl,a)$ % This procedure calls itself recursively! % q: a **-expression to be separated % old_histy: a list of elements {denom,v,{(divisor . expr_before),..}} % where a sequence of divisions through factors from the % list of divisors and differentiations wrt. v and % afterwards multiplication with denom resulted in q % ftem: functions in the expression % vl: variables in the expression % a: the variables for direct separation=car starp() % % RETURNS {{(c1.q1),(c2.q2),(c3.q3),..},{s1,s2,s3,..},{r1,r2,..},{f1,f2,..}} % where qi=0 are necessary consequences, % qi are not *-expressions, % sum_i ci*qi = q % si=0 are consistency conditions determining constants/functions % of integration % ri=0 are other cases to be checked --> case distinctions begin scalar histy,l1,l4,nv,vl1,nv1,h,x,f,ft,aa,bb,cc,y, ruli,extra_cond,par,cases,newf$ %--- ft: the list of functions to drop from q by differentiation %--- to do a direct separation wrt x: % x = any one variable in a on which a function with as % many as possible variables does not depend on % Find such a function and variable x ft:=ftem; nv:=0; while ft do << vl1:=fctargs car ft; nv1:=if vl1 then length vl1 else 0; if nv1 > nv then << h:=setdiff(a,vl1); if h then << x:=car h; % if possible find an x occuring explicitly in q while h and freeof(q,car h) do h:=cdr h; if h then x:=car h; f:=car ft; nv:=nv1 >> >>; ft:=cdr ft >>; if nv=0 then x:=car a; % no x was found if tr_gensep then <<terpri()$write "--- The aim is to separate directly w.r.t. ",x$ write " the expression : "$deprint list q >>$ % Find all functions ft in q depending on x ft:=nil$ for each f in ftem do if member(x,fctargs f) and not freeof(q,f) then ft:=cons(f,ft)$ ft:=fctsort reverse ft$ % sorting w.r.t. number of args ruli:=start_let_rules()$ %--- throwing out functions ft until ft=nil %--- or until the expression lost the *-property while ft do % for each function to get rid of % (possibly each time a different diff variable) if null (l1:=felim(q,car ft,ftem,vl)) then ft:=nil % to stop else << %prettyprint l1; for each h in cdadr l1 do % special extra cases if not freeoflist(car h,ftem) then cases:=cons(car h,cases); %write"cadr l1=",cadr l1$terpri()$ if zerop car l1 then << q:=reval {'QUOTIENT,cdr cadadr l1,car cadadr l1}; % new expression cc:=cons(car cadr l1,cddadr l1); >> else << q:=car l1$ % new expression cc:=cadr l1; >>$ if (pairp q) and (car q='QUOTIENT) then << bb:=caddr q; % we take off the denimonator q:=cadr q >> else bb:=1$ histy:=cons(cons(bb,cc),histy)$ % extended history %terpri()$write"q=",q$terpri()$ %write"histy=",histy$terpri()$ ftem:=smemberl(ftem,q)$ % functions still in q aa:=stardep(ftem,argset(ftem))$ % still *-expression? if not aa or zerop cdr aa then ft:=nil % to stop else ft:=smemberl(cdr ft,ftem) % remain. fcts of x >>$ stop_let_rules(ruli)$ if null l1 then if tr_gensep then <<terpri()$ write"felim or newgensep gave nil!!"$terpri()$ write"q=",q;terpri() >> else else RETURN << if tr_gensep then << terpri()$ write"Now ready for direct separation." >>; %--- prepare list of variables vl1 for direct separation vl1:=nil$ for each y in vl do if my_freeof(ftem,y) then vl1:=cons(y,vl1); %--- direct separation if useful (i.e. if aa(=stardep) neq nil) if vl1 and not (q=0) then <<if tr_gensep then <<terpri()$write "trying direct separation of "$ deprint list q$ write "Remaining variables: ",vl1>>$ l1:=separ(q,ftem,vl1,nil)$ % direct separation of the numerator if tr_gensep then <<terpri()$ write "The result of direct separation: "$ deprint for each bb in l1 collect cdr bb>>$ >> else l1:=list cons(1,q)$ if tr_gensep then << terpri()$ write"Separation gave ",length l1," condition(s)" >>; % Although the vaiable x does not occur anymore % (each felim call removed one function of x % and direct separation removed explicit occurences of x) % the remaining expression may still be indirectly separable % --> recursive call of partitn % l4 becomes a list of pairs (sep_coefficient . sep_remainding_factor) for each h in l1 do << ft:=smemberl(ftem,cdr h); vl1:=argset(ft)$ if null (aa:=stardep(ft,vl1)) then l4:=cons(h,l4) else << par:=partitn(cdr h, % expression append(histy, % history so far, old_histy), % needed to add new functions % of integration properly differentiated to be % able to integrate below ft, % functions vl1, % variables car aa % separation charac. ); % RETURNS {{(c1.q1),(c2.q2),(c3.q3),..},{s1,s2,s3,..}, % {r1,r2,..},{f1,f2,..} } % where qi=0 are necessary consequences, % qi are not *-expressions, % sum_i ci*qi = q % si=0 are consistency conditions determining constants/functions % of integration % ri=0 are other cases to be checked --> case distinctions l4:=append(l4,for each f in car par collect ({'TIMES,car h,car f} . cdr f)); extra_cond:=append(extra_cond,cadr par); % collecting any % appearing conditions cases:=append(cases,caddr par); newf:=cadddr par; ftem:=append(ftem,newf); >> >>$ %--- backintegration par:=backint(l4,old_histy,histy,ftem,vl)$ extra_cond:=append(extra_cond,cadr par); % collecting any conditions {car par,extra_cond,cases,append(newf,caddr par)} >> end$ %----------- symbolic procedure felim(q,f,ftem,vl)$ % returns: nil if not successful (quotient) otherwise % {the expression after all the divisions and differentiations, % {diff variable, sequence of (factor . expression before)} } begin scalar a,b,l,l1,ft1,v,prflag$ %--- getting rid of f through diff. wrt. v v:=car setdiff(vl,fctargs f)$ %--- ft1 are all v-independent functions %--- In the call to separ one has to disregard v-dep. functions ft1:=nil$ for each f in ftem do if my_freeof(f,v) then ft1:=cons(f,ft1)$ %--- To run separ, functions ft1 should not be in the denominator %--- ?????? What if nonl. Problems? if not (pairp q and (car q='QUOTIENT) and smemberl(ft1,caddr q)) then % This exceptional case should not occure anymore <<prflag:=print_$print_:=nil$ l:=separ(q,ft1,list v,nil)$ % det. all lin. ind. factors with v for each a in reverse idx_sort for each b in l collect cons(delength car b,b) collect cdr a$ if tr_gensep then <<terpri()$write "To get rid of ",f, " we differentiate w.r.t. : ",v>>$ print_:=prflag$ %--- l is a list of dotted pairs a each representing a term in q % where car(a) is the product of v-dep. factors and cdr(a) the % product of v-independent factors %--- A list l1 of car(a) is generated for which cdr(a) depends % on f. l1 is the list of divisions to be done before differen. l1:=nil$ while l do <<a:=car l$ b:=nil$ if not freeof(cdr a,f) and (not zerop car a) then l1:=cons(car a,l1)$ l:=cdr l >>$ if tr_gensep then <<terpri()$ write v," - depending coefficients of terms containing ",f," : "$ for each ss in l1 do eqprint ss>>$ %--- Now the divisions and differentiations are done while l1 do <<b:=reval aeval car l1$ %--- b is the v-dep. coefficient l1:=cdr l1$ %--- ????? If for non-linear Problems b includes ftem functions % then the extra case 0 = b has to be considered if not zerop b then <<a:=reval aeval list('QUOTIENT,q,b)$ %--- for later backward integrations: extension of the history l:=cons(b . q ,l)$ %--- new: q is the equ. before division & diff. % formerly: l:=cons(b ,l)$ % l will be returned later by felim %--- l1 has to be updated as the coefficients % change through division and differentiation l1:=for each c in l1 collect reval list('DF,list('QUOTIENT,c,b),v)$ %--- the differentiation q:=reval list('DF,a,v)$ if tr_gensep then <<write "The new equation: "$ eqprint q$ write "The remaining factors:"$ for each ss in l1 do eqprint ss >> >> >>$ %if l then part_histy:=cons(v,l)$ %--- output if tr_gensep then <<terpri()$write "new expression (should not depend on ",f,") : "$ eqprint q$>>$ if tr_gensep and l then <<write"To invert the last steps one has to integr. wrt. ",v$ terpri()$ write "each time before multiplying with "$ for each aa in l do eqprint car aa >>$ l1:=list(q,cons(v,l)) >>$ return l1 end$ symbolic procedure backint(l4,old_histy,histy,ftem,vl)$ % l4 is a list of pairs (sep_coefficient . % sep_remainding_factor_to_be_integrated) % old_histy, histy are lists of elements % {denom,v,{(divisor . expr_before),..}} % where a sequence of divisions through factors from the % list of divisors and differentiations wrt. v and % afterwards multiplication with denom resulted in q % Integrations should only be done inverting histy, but % in choosing functions of integration, both should be used % % returns {- integrated equivalent of l4 where the cdr of each element % is the integrated expression, % - a list of check_sum conditions, % - a list of new functions} begin scalar succ,ft,q,l,v,v1,vf,s1,s2,p,f1,f2,fctr,check_sum, allfnew,new_cond,denomi$ % start of the backintegration succ:=t$ while histy and succ do <<l:=car histy$ % l={diff variable, sequence of (factor . exp before)} histy:=cdr histy$ denomi:=car l$ v:=cadr l$ l:=cddr l$ % At first putting the denominator back in % l4 is a list of pairs (sep_coefficient . sep_remainding_factor) if denomi neq 1 then l4:=for each h in l4 collect (car h . {'QUOTIENT,cdr h,denomi}); if tr_gensep then <<terpri()$ write "backward integration w.r.t. ",v>>$ % Now the sequence of integrations wrt v % l is the list of (factor . expression before division & diff) while l do << % while l and q do fctr:=caar l; check_sum:=cdar l; l:=cdr l; if tr_gensep then <<terpri()$ write "The integrals of the following partitioned subexpressions"$ terpri()$ write "added up should be equal the original expression: "$ terpri()$ eqprint check_sum >>$ %write"l4="$ %prettyprint l4; % l4 is a list of pairs (sep_coefficient . sep_remainding_factor) l4:=for each h in l4 collect if null car h then h % one integration % was not succ.ful else << ft:=smemberl(ftem,cdr h)$ fnew_:=nil$ if tr_gensep then <<terpri()$ write "Backintegration of: "$eqprint cdr h >>$ q:=integratepde(cdr h,ft,v,nil,nil)$ % genflag:=nil, potflag=nil if null q then << succ:=nil$ if print_ then << terpri()$ write "#### Back integration of "$ eqprint cdr h$ write " wrt ",v," during generalized ", "separation was not successful ####."$ terpri()$ write "The coeff. dropped in direct separation was "$ mathprint car h >> >> else << q:=reval list('TIMES,fctr,car q)$ % fctr is the next integrating factor % Neccessary: Substituting the new functions of integration by % derivatives of them such that back-integration can be made % hat fnew_ nur ein element, d.h. wird nur eine Integration gemacht % oder mehrere? for each f1 in fnew_ do <<f2:=f1$ vf:=setdiff(vl,fctargs f1)$ for each s1 in reverse append(histy,old_histy) do <<v1:=cadr s1$ % The following also if not my_freeof(f1,v1) % The reason is that divisors may contain variables which % are later integration variables s2:=reverse cddr s1$ while s2 do <<% only divisions through factors that can be swallowed by f1 if not smemberl(vf,caar s2) then f2:=list('QUOTIENT,f2,caar s2)$ if not my_freeof(f1,v1) then f2:=reval list('DF,f2,v1)$ % actually only dividing through those factors of (caar s2) % which depend on v1 and which do not contain variables % which f2 does not depent on. s2:=cdr s2 >>$ if not smemberl(vf,car s1) then f2:=list('TIMES,f2,car s1)$ >>$ % the remaining integrations in the current element of histy if histy then << s2:=reverse l$ while s2 do <<if not smemberl(vf,caar s2) then f2:=list('QUOTIENT,f2,caar s2)$ if not my_freeof(f1,v1) then f2:=reval list('DF,f2,v1)$ s2:=cdr s2 >>; >>; if f1 neq f2 then <<if tr_gensep then <<terpri()$ write f1," is replaced by "$ eqprint f2>>$ q:=subst(f2,f1,q)$ >> >>$ allfnew:=append(fnew_,allfnew)$ ftem:=union(fnew_,ftem); if succ then check_sum:={'DIFFERENCE,check_sum,{'TIMES,q,car h}}; % car h is the coefficient dropped through direct separation >>$ (car h . q) % the value to be collected to give the new l4 >>; check_sum:=reval check_sum$ if succ then new_cond:=cons(check_sum,new_cond)$ if succ and tr_gensep then <<terpri()$ write "Consistency condition: "$eqprint check_sum >>$ >> >>$ for each f in allfnew do ftem_:=fctinsert(f,ftem_)$ if tr_gensep then if succ then <<terpri()$write "yields : "$ eqprint p$>> else <<terpri()$ write "was not successful.">>$ fnew_:=nil$ return {l4,new_cond,allfnew} end$ endmodule$ %********************************************************************* module gensep_non_lin$ %********************************************************************* % Routines for generalized separation of de's % Author: Thomas Wolf since 1997 symbolic procedure my_smemberl(p,vl)$ begin scalar l,v; for each v in vl do if not my_freeof(p,v) then l:=cons(v,l); return reverse l end$ %----------- symbolic procedure stripcond(conds)$ begin scalar newconds,condi; newconds:=nil; while conds do << condi:=cdar conds; conds:=cdr conds; if length condi=1 then condi:=car condi else condi:=cons('PLUS,condi); newconds:=cons(condi,newconds) >>; return newconds end$ %----------- symbolic procedure checkli(exlist,condi)$ begin scalar ok,isincondi,isinexli,excopy,n; ok:=t; while condi and ok do << % all i in the condition car condi isincondi:=car condi; %smemberl(ilist,car condi); n:=length isincondi; % are all isincondi contained in one of the elements of exlist? excopy:=exlist; while excopy and ok do << isinexli:=smemberl(isincondi,car excopy); if isinexli then if length(isinexli) = n then ok:=nil; excopy:=cdr excopy >>; condi:=cdr condi >>; return ok end$ %----------- symbolic procedure longst(exlist)$ % returns the element of exlist which (is a list and) % has the most elements begin scalar lo; while exlist do << if not lo then lo:=car exlist else if length(lo)<length(car exlist) then lo:=car exlist; exlist:=cdr exlist >>; return lo end$ %----------- symbolic procedure starequ(n,alindep,blindep)$ % alindep is a list of lists of factors ai which are all non-zero and % are all linear independent from each other within such a list % blindep like alindep % generates all cases each with all conditions with _i representing % ai or bi, equations and new functions are not generated begin comment The equation to separate has the form 0 = sum_i ai*bi where the bi do not depend on some variable x. The procedure starequ generates cases: cases ... ( all cases ) each case ... ( list of all a-conditions, list of all b-conditions) each condition ... ( the ai,bi contained in the condition with _i representing ai and bi ) ; scalar i,j,cases,oldcases,case,ai,bi,ci,oldaconds,oldbconds, newaconds,newbcond,newbconds,newacond, ilist,cona,conb,unin,el,pri; % ,oldpri % Determine the longest union of two list, one, ai, being element of % alindep and one, bi, being from blindep %pri:=t; i:=0; if alindep then for each cona in alindep do if blindep then for each conb in blindep do if (j:=length union(cona,conb)) > i then <<ai:=cona;bi:=conb;i:=j>> else else % no blindep given if (j:=length cona) > i then <<ai:=cona;i:=j>> else else % no alindep given if blindep then for each conb in blindep do if (j:=length conb) > i then <<bi:=conb;i:=j>>; % ai, bi will now be determined % preparation of the sequence ilist of extensions ilist:=for i:=1:n collect i; if pri then <<write"222";terpri()>>$ if i neq 0 then << if ai then i:=length ai else i:=0; if bi then j:=length bi else j:=0; unin:=union(ai,bi); % extensions through bch should be done when elements from % bi are treated. This is coded in ilist through negative numbers ilist:=setdiff(ilist,unin); if i>j then << for each el in setdiff(unin,ai) do ilist:=cons(-el,ilist); for each el in ai do ilist:=cons( el,ilist) >> else << for each el in setdiff(unin,bi) do ilist:=cons( el,ilist); for each el in bi do ilist:=cons(-el,ilist) >>; ilist:=reverse ilist >>; % ilist is prepared now if pri then <<write"333 ilist=",ilist;terpri()>>$ % `cases' is a list of cases, each is a dotted pair with % the car being the a-conditions and cdr being the b-conditions % The first two cases i:=car ilist; if i<0 then i:=-i; ci:=mkid('_,i); cases:=list(list(list(list(ci)), nil ), list( nil , list(list(ci))) ); % oldpri:=print_; % print_:=nil; ilist:=cdr ilist; if pri then <<write"555";terpri()>>$ while ilist do << i:=car ilist;ilist:=cdr ilist; if pri then <<write"iii=",i;terpri()>>$ if i>0 then ci:=mkid('_, i) else ci:=mkid('_,-i); if pri then << write"666 car ilist=",i; terpri() >>$ % if i>0 then the list of cases is extended with ai else with bi oldcases:=cases; cases:=nil; while oldcases do << % for each old case do: case:=car oldcases; if pri then <<write"old case:",case;terpri()>>$ oldcases:=cdr oldcases; if i>0 then << oldaconds:=car case; if pri then <<write"888 oldaconds=",oldaconds;terpri()>>$ % at first add condition i=0 to the case cases:=cons(cons(cons(list(ci),oldaconds), cdr case) , cases); if pri then <<write"999 new case=",car cases;terpri()>>$ % then: - add to each a-condition i % - add one new b-condition with the first element `i' % and furtherelements which are the first elements of % the a-lists which have been extended newaconds:=nil; newbcond:=list(ci); while oldaconds do << j:=caar oldaconds; newaconds:=cons(cons(j,cons(ci,cdar oldaconds)), newaconds); newbcond:=cons(j,newbcond); oldaconds:=cdr oldaconds >>; if pri then <<write"newaconds=",newaconds, " rev newbcond=",reverse newbcond; terpri()>>$ cases:=cons(list(newaconds, cons(reverse newbcond,cadr case)) , cases); if pri then <<write"000 new case=",car cases;terpri()>>$ >> else << oldbconds:=cadr case; if pri then <<write"888 oldbconds=",oldbconds;terpri()>>$ % at first add condition bi=0 to the case cases:=cons(list(car case, cons(list(ci),oldbconds)), cases); if pri then <<write"999 new case=",car cases;terpri()>>$ % then: - add to each b-condition i % - add one new a-condition with the first element `i' % and further elements which are the first elements of % the b-lists which have been extended newbconds:=nil; newacond:=list(ci); while oldbconds do << j:=caar oldbconds; newbconds:=cons(cons(j,cons(ci,cdar oldbconds)), newbconds); newacond:=cons(j,newacond); oldbconds:=cdr oldbconds >>; cases:=cons(list(cons(reverse newacond,car case), newbconds), cases); if pri then <<write"000 new case=",car cases;terpri()>>$ >> >>; >>; % Throwing out cases which are forbidden by alindep and blindep alindep:=for each ci in alindep collect for each i in ci collect mkid('_,i); blindep:=for each ci in blindep collect for each i in ci collect mkid('_,i); oldcases:=nil; % ilist:=for i:=1:n collect mkid('_,i); while cases do << if checkli(alindep,caar cases%,ilist ) then if checkli(blindep,cadar cases%,ilist ) then oldcases:=cons(car cases,oldcases); cases:=cdr cases >>; % print_:=oldpri; return oldcases end$ % of starequ %----------- symbolic procedure pickfac(ex,indx)$ % returns the n'th element of ex where indx has the form _n nth(ex,compress cdr explode indx)$ %----------- symbolic procedure find_cond(bcons,ai)$ % find the element of bcons with ai as (automatically first) element % (there must be an b-condition with ai as first element % if that has not already been dropped % because ai is not the first element of the a-condition) begin while (pairp bcons) and (caar bcons neq ai) do bcons:=cdr bcons; return if pairp bcons then car bcons else nil end$ %----------- symbolic procedure starsep(exx,ex,ftem,vl,x)$ % exx is the original expression to be separated % ex is a *-expression returned from SEPAR % vl are all variables which really occur in ex or % on which ex depends on % x is a variable on which the bi do not depend on % % RETURNS a list of cases, each case having the form % {{new constants/functions of integration}, % eq1,eq2,eq3,...} % where eqi are a set of necc and suff conditions begin scalar cases,newcases,acons,bcons,acond,newca,alindep,blindep,aco,bco, ai,bi,ci,a1,avars,bvars,i,ili,cilist,ali,n,addex,bcopy,cntr, pri; % pri:=t; ili:=for i:=1:length ex collect mkid('_,i); n:=length vl; % at first determining lists of non-vanishing and linear independent % a-factors alindep and b-factors blindep % does so far only the obvious test which is useful essentially for % linear problems cntr:=0; for each ci in ex do << cntr:=add1 cntr; if pri then << write"a",cntr," = "$mathprint car ci$ write"b",cntr," = "$mathprint cdr ci$ >>$ if null smemberl(ftem,car ci) then alindep:=cons(cntr,alindep)$ if null smemberl(ftem,cdr ci) then blindep:=cons(cntr,blindep)$ >>; if alindep then alindep:=list alindep$ if blindep then blindep:=list blindep$ % generation of all logical cases with the factors ai,bi % substituted by _i cases:=starequ(length ex,alindep,blindep); if pri then <<terpri()$write"Returned from STAREQU: ",cases>>; % newcases will be the new list of all cases newcases:=nil; while cases do << % car cases is one case with alltogether n conditions which % The conditions for the a-factors are called below acons % and for the b-factors bcons. acons:= caar cases; bcons:=cadar cases; cases:= cdr cases; if pri then <<write"acons=",acons," bcons=",bcons;terpri()>>$ % The case will now be formulated with the terms of the expression ex newca:=nil; cilist:=nil; addex:=nil; bcopy:=nil; % extract at first all b-conditions with only one term as they do not % need constants of integration and are used for any grade of ex while bcons do << if length car bcons=1 then newca:=cons(cdr pickfac(ex,caar bcons),newca) else bcopy:=cons(car bcons,bcopy); bcons:=cdr bcons >>; bcons:=bcopy; % The a-factors may depend on all variables vl whereas the % b-factors do at least not depend on x. while acons do << % formulating all a-conditions aco:=car acons; % aco encodes one condition for a-factors acons:=cdr acons; if pri then <<write"aco=",aco;terpri()>>$ a1:=car aco; % e.g. a1 = _7 acond:=list car pickfac(ex,a1); if pri then <<write"acond=",acond;terpri()>>$ % acond becomes a list of all a-factors encoded by aco % adding all a-conditions to the new case newca if (length aco)=1 then newca:=cons(car acond,newca) else << ali:=for each i in aco collect car pickfac(ex,i); avars:=my_smemberl(ali,vl); if length avars neq n then << % The progress in this case is that all new equations will % have less variables. % Now determination on which variables the constants of back- % integration would depend on. This is the intersection of all % variables avars in the a-condition and the variables bvars in % the b-condition in which the constant of integration occurs. The % a-condition is known, it will be made from acond. The relevant % b-condition is determined through the current index of aco % from which the coefficient is to be determined (which is not the % first index of aco) aco:=cdr aco; % a1 already in acond while aco do << ai:=car aco;aco:=cdr aco; % find the list of variables bvers of the b-condition bco % which includes the b-factor corresponding to ai=car aco % disadvantage of this way: if bco has m elements then all % variables of bco are determined m-1 times. % determining bco as the b-condition which contains ai bco:=find_cond(bcopy,ai); % bcopy is used instead of bcons % the condition with ai may already have been dropped from % bcons because of ai depending on all variables, i.e. the % new functions in the b-conditions would depend on all % variables and be no help. if pri then <<write"bco=",bco;terpri()>>$ % bcoli:=smemberl(ili,bco); % in case bco has already had subst. % determining bvars bvars:=nil; for each bi in bco do bvars:=union(my_smemberl(cdr pickfac(ex,bi),vl),bvars); if pri then <<write"bvars=",bvars$terpri()>>$ % introducing new constants of integration ci:=newfct(fname_,intersection(avars,bvars),nfct_)$ cilist:=cons(ci,cilist); nfct_:=nfct_+1; acond:=cons(list('MINUS,list('TIMES,ci,car pickfac(ex,ai))),acond); if pri then <<write"acond=",acond;terpri()>>$ if bco:=find_cond(bcons,ai) then << bcons:=subst(subst(list('TIMES,ci,a1),a1,bco),bco,bcons); if pri then <<write"bcons=",bcons;terpri()>>$ >> >>; acond:=cons('PLUS,acond) >> else << % The condition aco is a *-condition % progress in this case is that new a-conditions % have less functions addex:=t; ali:=reverse ali; aco:=reverse aco; while length ali > 1 do << if pri then <<write"ali1=",ali;terpri()>>$ % Generate the a-condition if pri then <<write"###";terpri()>>$ ali:= if not zerop car ali then for each i in cdr ali collect reval list('DF,list('QUOTIENT,i,car ali),x) else cdr ali; if pri then <<write"ali2=",ali;terpri()>>$ % Drop that element from bcons which includes % car ali (as first element) if bco:=find_cond(bcons,car aco) then bcons:=setdiff(bcons,list bco); aco:=cdr aco >>; acond:=car ali; if (pairp acond) and (car acond = 'QUOTIENT) then acond:=cadr acond; >>; newca:=cons(acond,newca) >>; if pri then <<write"newca1=",newca;terpri()>>; >>; % all a-conditions have been dealt with if pri then <<write"newca2=",newca;terpri()>>; % completing all b-conditions for each bi in ili do bcons:=subst(cdr pickfac(ex,bi),bi,bcons); % adding all b-conditions to the new case newca while bcons do << if length car bcons = 1 then newca:=cons(caar bcons,newca) else newca:=cons(cons('PLUS,car bcons), newca); bcons:=cdr bcons >>; % if ex is an *-expression with grade>1 then possibly b-conditions % had to be dropped, so ex must be added if addex then newca:=cons(exx,newca); if pri then <<write"newca3=",newca;terpri()>>; % adding the list of constants of integeration newca:=cons(cilist,newca); if pri then <<write"cilist=",cilist;terpri()>>; newcases:=cons(newca,newcases) >>; return newcases end$ % of starsep %----------- symbolic procedure separizable(p,ftem,vl)$ begin scalar x,ft,f,ex,v,a,b,vlcp,allvarcaara,print_bak$ vlcp:=vl; repeat << x:=car vl; vl:=cdr vl; % Determining all x-dependent functions ft ft:=nil; for each f in ftem do if member(x,fctargs f) and not my_freeof(p,f) then ft:=cons(f,ft)$ f:=car reverse fctsort ft$ % sorting w.r.t. number of args v:=car setdiff(vlcp,fctargs f)$ % getting rid of f by v-differen. % preparation of the separ-call, ft are now v-indep. functions ft:=nil$ for each f in ftem do if my_freeof(f,v) then ft:=cons(f,ft)$ % ex:=separ(p,ft,list v,nil)$ % det. all lin. ind. factors print_bak:=print_; print_:=nil; ex:=separ2(p,ft,list v)$ % det. all lin. ind. factors print_:=print_bak; a:=ex; while a and << b:=vlcp; while b and not my_freeof(caar a,car b) do b:=cdr b; b >> do a:=cdr a; if a then allvarcaara:=cons(caar a,allvarcaara); >> until (null a) or (null vl); % if a then null vl then whatever x was, there is allways one % element (car a) of ex such that car of this element (caar a) % does depend on all variables --> no separability possible, % new functions would depend on all variables % if a then test whether separation would be possible by getting % rid of functions through differentiation % (this would not be the case if e.g. sin(all variables) would occur) % --> use of smemberl vl:=vlcp; while allvarcaara and not not_included(vlcp,smemberl(vlcp,car allvarcaara)) do <<allvarcaara:=cdr allvarcaara; vl:=cdr vl>>$ return if a and null allvarcaara then nil % no chance else if a then {nil,car allvarcaara,car vl} % deleting functions first else << % separation now possible if tr_gensep then <<terpri()$write "To separate directly wrt. ",x$ write " the expression : "$deprint list p$ write "will be differentiated wrt. ",v," to get rid of ",ft," ">>$ {ex,v} >> end$ %----------- symbolic procedure newgensep(p,starpro,ftem,vl)$ % ftem, vl should be correct: % ftem:=smemberl(ftem_,p)$ % vl:=varslist(p,ftem,vl)$ % starpro:=stardep(ftem,vl)$ % returns what starsep returns begin scalar pl,v,ex,a,b$ % ,gense,el1,el2,conds,newcali,l,clist$ % if pairp p and (car p = 'QUOTIENT) then <<casecheck(caddr p,vl)$ % p:=cadr p>>$ % ftem:=smemberl(ftem,p)$ % vl:=varslist(p,ftem,vl)$ % if not (starpro:=stardep(ftem,vl)) then % then no *-equation % pl:=list list(nil,p) % one case, no new functions % else % e.g. starpro=((x y z).2) % if zerop cdr starpro then pl:=nil% ############################## % %list cons(nil,separate(p,ftem,vl)) % direct sep % else % if delength(p) leq gensep_ then % generalized separation % << if print_ then <<terpri()$write "generalized separation ">>$ if tr_gensep then <<terpri()$write "de to be separated : "$ deprint list p$ terpri()$write "variable list for separation : ",car starpro$ terpri()$write "for each of these variables ",cdr starpro, " function(s) depend(s) on it.">>$ for each v in car starpro do vl:=delete(v,vl); vl:=append(car starpro,vl); a:=separizable(p,ftem,vl)$ if null a then return nil else if null car a then return << % functions to be deleted before separation are those in cadr a % ft:=smemberl(ftem,cadr a); if print_ then <<terpri()$ write"In order to be separable with this procedure at first"$terpri()$ write"one or more functions have to be eliminated through"$terpri()$ write"differentiation and algebraic elimination, for example,"$terpri()$ write"the functions: ",smemberl(ftem,cadr a)$terpri()$ >>; nil >> else <<ex:=car a;v:=cadr a>>$ for each a in reverse idx_sort for each b in ex collect cons(delength car b,b) collect cdr a$ if tr_gensep then <<terpri()$write "Return from SEPAR: "$terpri()$prettyprint ex>>$ % with v and v-dep. functions as first factors in the terms in ex pl:=starsep(p,ex,ftem,vl,v); if tr_gensep then <<terpri()$write "Return from STARSEP: "$terpri()$prettyprint pl>>$ % print_:=oldpri$ %%############################################################ % % l is a list of cases each (list of new fncts, cond1, cond2, ...) % % for each condition (neq p) in all cases calling gensep again % % if needed % pl:=nil; % the final list of cases of only non-*-equ. % while l do % checking all cases % <<clist:=caar l; % list of new functions % conds:=cdar l; % list of conditions % l:=cdr l; % newcali:=list list clist; % % newcali will be the list of new cases which starts as % % only one case and from this only the list of new functions % % but no conditions % while conds do % checking all conditions of one case % << %% if car conds = ex then %% % ex aready investigated, not again %% % append ex to the conditions of all new cases %% newcali:=for each el1 in newcali collect %% cons(car el1,cons(ex,cdr el1)) %% else << % gense:=newgensep(car conds,append(ftem,clist),vl); % newcali:=for each el1 in gense join % for each el2 in newcali collect % cons(append(car el1,car el2), % append(cdr el1,cdr el2)); %% >>; % conds:=cdr conds % >>; % pl:=append(newcali,pl) % >> % >>; return pl end$ % of newgensep %----------- symbolic procedure gen_separation2(arglist)$ % Indirect separation of a pde, 2nd version for non-linear PDEs begin scalar p,h,fl,l,l1,pdes,forg,n,result,d,contrad,newpdes$%,newfdep,bak,sol pdes:=car arglist$ forg:=cadr arglist$ if expert_mode then << l1:=selectpdes(pdes,1); if get(car l1,'starde) then flag(l1,'to_gensep) >> else l1:=pdes$ if (p:=get_gen_separ_pde(l1,nil,nil)) then if l:=newgensep(get(p,'val), get(p,'starde), get(p,'fcts), get(p,'vars)) then if cdr l then << if print_ then << terpri()$ write"The indirect separation leads to ",length l," cases."$ %terpri()$ >>$ contrad:=t$ n:=0; remflag1(p,'to_gensep)$ % bak:=backup_pdes(pdes,forg)$ % must come before drop_pde(... backup_to_file(pdes,forg,nil)$ % newfdep:=nil$ while l do << d:=car l; l:=cdr l; if not memberl(cdr d,ineq_) then << % non of the equations is an inequality if n neq 0 then << h:=restore_and_merge(l1,pdes,forg)$ pdes:= car h; forg:=cadr h; % was not assigned above as it has not changed probably % h:=restore_pdes(bak); % pdes:=car h; forg:=cadr h >>; n:=n+1$ level_:=cons(n,level_)$ if print_ then << print_level(t)$ terpri()$ write "CRACK is now called with the assumption : "$ deprint(cdr d) >>$ % formulation of new equations for each h in car d do ftem_:=fctinsert(h,ftem_); fl:=append(get(p,'fcts),car d); newpdes:=pdes$ for each h in cdr d do newpdes:=eqinsert(mkeq(h,fl,vl_,allflags_,t,list(0),nil,newpdes),newpdes); % further necessary step to call crackmain(): recycle_fcts:=nil$ % such that functions generated in the sub-call % will not clash with existing functions l1:=crackmain(newpdes,forg)$ % for each sol in l1 do % if sol then << % for each f in caddr sol do % if h:=assoc(f,depl!*) then newfdep:=cons(h,newfdep); % >>; if not contradiction_ then contrad:=nil$ if l1 and not contradiction_ then result:=union(l1,result); contradiction_:=nil$ >> >>; delete_backup()$ % % newfdep are additional dependencies of the new functions in l1 % depl!*:=append(depl!*,newfdep); contradiction_:=contrad$ if contradiction_ then result:=nil$ if print_ then << terpri()$ write"This completes the investigation of all cases of an ", "indirect separation."$ terpri()$ >>$ result:=list result % to tell crackmain that computation is completed >> else << % only one case l:=car l; for each h in car l do ftem_:=fctinsert(h,ftem_); fl:=append(get(p,'fcts),car l); pdes:=drop_pde(p,pdes,nil)$ for each h in cdr l do pdes:=eqinsert(mkeq(h,fl,vl_,allflags_,t,list(0),nil,pdes),pdes); result:=list(pdes,forg) >>$ return result$ end$ endmodule$ end$