Artifact 5fac6ff3f8bc6c13b07e610270360630ad7b48a475545eb722a8c2441ab9903e:
- Executable file
r37/packages/crack/crint.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: 65112) [annotate] [blame] [check-ins using] [more...]
%********************************************************************* 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$