Artifact 4368087e7b7333cd9a4a231be9b72c6234bb0464401392b86f3a0b9f165f6515:
- Executable file
r38/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: 91100) [annotate] [blame] [check-ins using] [more...]
%********************************************************************* module integration$ %********************************************************************* % Routines for integration of pde's % Authors: Andreas Brand 1993 1995 % Thomas Wolf since 1993 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 = 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 = 1 then f:=caar h else if length lds = 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 = 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: v, 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 if freeabs_ and (null freeof(qih,'ABS)) 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 pri 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 print_ and null pih then << terpri()$ write"Inhomogeneous part: "$ mathprint qih$ write"can not be integrated explicitly wrt. ",v$ >>$ 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)$ begin scalar ft,ip,h,itgl1,itgl2$ if potint then return intpde_(p,ftem,vl,x,potint)$ % ft are functions of x for each h in ftem do if not freeof(assoc(h,depl!*),x) then ft:=cons(h,ft); ip:=int_partition(p,ft,x)$ if null cdr ip then return intpde_(p,ftem,vl,x,potint)$ while ip do << h:=intpde_(car ip,ftem,vl,x,potint)$ if null h then << ip:=nil; itgl1:=nil; itgl2:=nil >> else << % itgl1:=cons(car h,itgl1); % itgl2:=cons(cadr h,itgl2); itgl1:=nconc(list(car h),itgl1); itgl2:=nconc(list(cadr h),itgl2); ip:=cdr ip >>; >>$ if itgl1 then <<itgl1:=reval cons('PLUS,itgl1); itgl2:=reval cons('PLUS,itgl2) >>$ return if null itgl1 then nil else {itgl1,itgl2} end$ symbolic procedure drop_x_dif(der,x)$ begin scalar dv,newdv$ % der is a derivative like {'DF,f,u,2,x,3,y,z,2} or {'DF,f,u,2,x,y,z,2} % drops the x-derivative(s) dv:=cddr der; while dv do << if car dv=x then if cdr dv and fixp cadr dv then dv:=cdr dv else else newdv:=cons(car dv,newdv); dv:=cdr dv >>; return if newdv then cons('DF,cons(cadr der,reverse newdv)) else cadr der end$ symbolic procedure strip_x(ex,ft,x)$ begin scalar tm,cex$ % ex is a term % ft is a list of all functions of x which possibly occur in ex return if freeoflist(ex,ft) then 1 else if not pairp ex then ex else if car ex='MINUS then strip_x(cadr ex,ft,x) else if car ex='DF then drop_x_dif(ex,x) else if car ex='EXPT then if not pairp cadr ex then ex else if caadr ex='DF then {'EXPT,drop_x_dif(cadr ex,x),caddr ex} else 1 % strange else if car ex='TIMES then << ex:=cdr ex; while ex do << cex:=car ex; ex:=cdr ex; if not freeoflist(cex,ft) then if not pairp cex then tm:=cons(cex,tm) else if car cex='DF then tm:=cons(drop_x_dif(cex,x),tm) else if car cex='EXPT then if not pairp cadr cex then tm:=cons(cex,tm) else if caadr cex='DF then tm:=cons({'EXPT,drop_x_dif(cadr cex,x),caddr cex},tm) % else strange - no polynomial in ft >>; if null tm then 1 % strange else if length tm > 1 then reval cons('TIMES,tm) % product of factors else car tm % single factor >> else 1 % strange end$ symbolic procedure sort_partition(pname,p,ft,x)$ % The equation is either given by its name pname or its value p % if keep_parti=t then the partitioning will be stored begin scalar stcp,pcop,parti; % parti will be the list of partial sums if pname then if get(pname,'partitioned) then return get(pname,'partitioned) else p:=get(pname,'val)$ if (not pairp p) or (car p neq 'PLUS) then p:=list p else p:=cdr p; while p do << % sort each term into a partial sum stcp:=strip_x(car p,ft,x); % first strip off x_dependent details pcop:=parti; % search for the label in parti while pcop and caar pcop neq stcp do pcop:=cdr pcop; if null pcop then parti:=cons({stcp,1,{car p}},parti) % open a new partial sum else rplaca(pcop,{stcp,add1 cadar pcop, cons(car p,caddar pcop)}); % add the term to an existing partial sum p:=cdr p >>; if pname and keep_parti then put(pname,'partitioned,parti)$ return parti end$ % of sort_partition symbolic procedure int_partition(p,ft,x)$ begin scalar parti,ft,pcop; % the special case of a quotient if (pairp p) and (car p='QUOTIENT) then return if not freeoflist(caddr p,ft) then list p else << pcop:=int_partition(cadr p,ft,x)$ for each h in pcop collect {'QUOTIENT,h,caddr p} >>$ parti:=sort_partition(nil,p,ft,x)$ parti:=idx_sort for each h in parti collect cdr h; return for each h in parti collect if car h = 1 then caadr h else cons('PLUS,cadr h) end$ 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 <<if (pairp p) and (car p = 'QUOTIENT) then << l1:=lderiv(cadr p,f,vl)$ % numerator if cdr l1 neq 'INFINITY then << l2:=lderiv(caddr p,f,vl)$ % denomiator if cdr l2 = 'INFINITY then l1:=l2 else << if car l1 and (car l1=car l2) then l1:=(car l1 . 2) % nonlinearity else << l:=lderiv({'PLUS,car l1,car l2},f,vl)$ % comparison of both if car l2 % changed from car l1 as car l1 gives endless loops with % call intpde((quotient (minus 1) (df u y)),(u),(y x t),y,nil) and (car l = car l2) then l1:=(car l2 . -1) >> >> >>; l:=l1; >> else l1:=l:=lderiv(p,f,vl)$ while not (flag or << iflag:=intlintest(l,x); if (iflag='NOXDRV) or (iflag='NODRV) then << l2:=start_let_rules()$ p:=reval aeval p$ stop_let_rules(l2)$ l:=lderiv(p,f,vl)$ iflag:=intlintest(l,x) >>; iflag >> ) do <<if (pairp p) and (car p = 'QUOTIENT) then k:=reval {'QUOTIENT,coeffn(cadr p,car l,cdr l),caddr p} else 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 diffrelp(l1,(l:=lderiv(p,f,vl)),vl) 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 if freeabs_ and (null freeof(s,'ABS)) then nil else << k:=start_let_rules()$ !#if (equal version!* "REDUCE 3.6") l2:=reval aeval list('DF,l1,x)$ if 0 neq reval reval aeval list('DIFFERENCE,l,l2) then << !#else l2 := reval {'DF,l1,x} where !*precise=nil; if 0 neq (reval {'DIFFERENCE,l,l2} where !*precise=nil) then << !#endif 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,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 (null flagp(q,'to_fullint)) then nil else if get(q,'starde) then nil else if (null get(q,'allvarfcts)) % or (cdr get(q,'allvarfcts)) % if more than one allvar-function 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)$ if fullint and (null dl) then remflag1(q,'to_fullint)$ return dl end else t$ symbolic procedure integrate(q,genintflag,fullint,pdes)$ % 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 % parameter pdes only for drop_pde_from_idties(), drop when pdes_ global % and for mkeqlist() for adding inequalities begin scalar l,fli,fnew_old$ if fli:=intflagtest(q,fullint) then <<if fullint then <<fnew_old:=fnew_;fnew_:=nil>>$ if (l:=integratede(get(q,'val),get(q,'fcts),genintflag)) then if fullint and not null car ldiffp(car l,car fli) then <<remflag1(q,'to_fullint); for each f in fnew_ do drop_fct(f)$ fnew_:=fnew_old; l:=nil; if print_ then << terpri()$write"Not enough integrations to solve for a function. " >> >> else <<fnew_:=union(fnew_old,fnew_)$ for each f in fnew_ do ftem_:=fctinsert(f,ftem_)$ fnew_:=nil$ flag1(q,'to_eval)$ update(q,car l,ftem_,get(q,'vars),t,list(0),nil)$ drop_pde_from_idties(q,pdes,nil)$ drop_pde_from_properties(q,pdes)$ l:=cons(q,mkeqlist(cdr l,ftem_,get(q,'vars), allflags_,t,get(q,'orderings),pdes))$ put(q,'dec_with,nil); % 23.3.99 added --> cycling? put(q,'dec_with_rl,nil); % " added --> cycling? if print_ then << terpri()$ if cdr l then if get(q,'nvars)=get(cadr l,'nvars) then write "Potential integration of ",q," yields ",l else write "Partially potential integration of ",q," yields ",l else write "Integration of ",q$ terpri()>>$ remflag1(q,'to_fullint)$ remflag1(q,'to_int) >> else << remflag1(q,'to_fullint)$ remflag1(q,'to_int) >> >>$ % if print_ and null l and fullint then terpri()$ % prints unnecc. nl return l$ end$ symbolic procedure quick_integrate_one_pde(pdes)$ begin scalar q,p,r,nv,nvc,miordr,minofu,minodv,ordr,nofu,nodv$ % ,nvmax$ % nvmax:=0; % for each q in ftem_ do if (r:=fctlength q)>nvmax then nvmax:=r; nv:=no_fnc_of_v()$ % the number of functions for each variable % find the lowest order derivative wrt. only one variable miordr:=10000; % the order of the currently best equation minofu:=10000; % the number of functions depending on the % variable wrt. which shall be integrated minodv:=10000; % the number of differentiation variables of % the so far best equation while pdes and (get(car pdes,'length) = 1) do << % only 1 term q:=get(car pdes,'derivs)$ if q and % (get(car pdes,'nvars) = nvmax) cdaar q % any differentiations at all then << q:=caar q$ nodv:=0$ % no of differentiation variables ordr:=0$ % total order of the derivative r:=cdr q; while r do << if fixp car r then ordr:=ordr-1+car r else <<ordr:=add1 ordr;nodv:=add1 nodv>>; r:=cdr r >>$ if nodv>1 then nofu:=10000 % nodv = no of functions depending else << % on the integration variable nvc:=nv; while cadr q neq caar nvc do nvc:=cdr nvc; nofu:=cdar nvc; >>$ % no of fncs of v if nodv=1 then if (ordr<miordr) or ((ordr=miordr) and (nodv<minodv) or ((nodv=minodv) and (nofu<minofu))) then <<p:=car pdes; minofu:=nofu; miordr:=ordr; minodv:=nodv >>; >>$ pdes:=cdr pdes >>$ if p then p:=integrate(p,nil,t,pdes)$ return p end$ symbolic procedure integrate_one_pde(pdes,genintflag,fullint)$ % trying to integrate one pde begin scalar l,l1,m,p,pdescp$ % ,nvmax,h,f$ % nvmax:=0; % for each f in ftem_ do if (h:=fctlength f)>nvmax then nvmax:=h; % at first selecting all eligible de's m:=-1$ pdescp:=pdes$ while pdescp do << if flagp(car pdescp,'to_int) and not(get(car pdescp,'starde)) then << l:=cons(car pdescp,l); if get(car pdescp,'nvars)>m then m:=get(car pdescp,'nvars)$ >>; pdescp:=cdr pdescp >>; l:=reverse l; if mem_eff then % try the shortest equation first while l do if p:=integrate(car l,genintflag,fullint,pdes) then l:=nil else l:=cdr l else % find an equation to be integrated with as many as possible variables % if (m=nvmax) or (null fullint) then while m>=0 do << l1:=l$ while l1 do if (get(car l1,'nvars)=m) and (p:=integrate(car l1,genintflag,fullint,pdes)) then << m:=-1$ l1:=nil >> else l1:=cdr l1$ % if fullint then m:=-1 else 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(if q1=1 then car l1 else cons('TIMES, append(if pairp q1 and car q1='TIMES then cdr q1 else list q1, if pairp car l1 and caar l1='TIMES then cdar l1 else list 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; err_catch_solve(eqs,hh) >>; 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:=err_catch_solve(eqs,hh) >>; 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 procedure add_factors(gk,p)$ % gk is a list of factors and p anything but a quotient if null p then gk else if ( not pairp p ) or (( pairp p ) and (car p neq 'TIMES) ) then union(list if (pairp p) and (car p = 'MINUS) then cadr p else p,gk) else <<for each h in cdr p do gk:=union(list if (pairp h) and (car h = 'MINUS) then cadr h else h,gk); gk >>$ %------------------ 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,precise_old, carftr; % exfactors is the list of factors extracted at the beginning % pri:=t; !#if (equal version!* "REDUCE 3.6") !#else precise_old:=!*precise$ !*precise:=nil$ !#endif 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); carftr:=car ftr; if not evalnumberp carftr then if (pairp carftr) and (car carftr='QUOTIENT) then gk:=add_factors(add_factors(gk,cadr carftr), caddr carftr ) else gk:=add_factors(gk,carftr); >> 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 err_catch_fac(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: "$ >>; !#if (equal version!* "REDUCE 3.6") !#else !*precise:=precise_old$ !#endif 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)$ begin scalar newde,newnewde,l,h,newcond$ h:= % the integrated equation if not xnew then << % Integr. einer alg. Gl. fuer eine Abl. newde:=cadr solveeval list(de,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(de,fold,xnew,ordr,ftem) % --> ode fuer eine Abl. von f else << newde:=odeconvert(de,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 >> >>; return if not h then nil else cons(h,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 % length get(p,'derivs) ? 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:= reval 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()$ if !*rational then << off rational; newde:=reval newde; on rational>> else 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 freeabs_ and null freeof(newde,'ABS) 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 pairp newde and (car newde = 'EQUAL) then if (pairp cadr newde) and (caadr newde = 'QUOTIENT) and (zerop caddr newde) then newde:={'EQUAL,cadadr newde,0} else if (pairp caddr newde) and (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$ %******************************************************************** module divintegration$ %******************************************************************** % Routines to write a given expression as divergence % Author: Thomas Wolf % 1998 symbolic operator intcurrent1$ % used in conlaw2,4 symbolic procedure intcurrent1(divp,ulist,xlist,dulist, nx,eqord,densord)$ % computes a list in reverse order from which the conserved % current is computed through integration begin scalar ii,il,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11,contrace,u, nega,pii,mo,pl,nu$ %contrace:=t; xlist:=cdr xlist; ulist:=cdr ulist; nu:=length ulist; mo:=if eqord>densord then eqord-1 else densord-1$ pl:=for ii:=1:nx collect << % the components W^i il:=nil; pii:=nil; repeat << h11:=cons(ii,il); h1:=derili(h11); h11:=combi(sortli(h11)); if contrace then <<write"==== ii=",ii," il=",il," h1=",h1," h11=",h11;terpri()>>; h2:=il; h3:=nil; if null h2 then h4:=list(nil . nil) else << h4:=list(car h2 . nil); while h2 do << h3:=cons(car h2,h3);h2:=cdr h2; h4:=cons((if h2 then car h2 else nil ) . derili(h3),h4)$ >>; >>; if contrace then <<write"h3=",h3," h4=",h4;terpri()>>; for e1:=1:nu do << % for each function u u:=nth(ulist,e1); h5:=mkid(u,h1); if contrace then <<write"h5=",h5;terpri()>>; % h6:=nil; if contrace then <<write"h6-1=", list('DF,divp,h5); terpri()>>; % h6:=cons(reval list('QUOTIENT,list('DF,divp,h5),h11), h6); % if nde=1 then h6:=car h6 % else h6:=cons('PLUS,h6); h6:=reval list('QUOTIENT,list('DF,divp,h5),h11); if contrace then <<write"h6=",h6;terpri()>>; nega:=nil; h7:=h4; % h7 is in such an order that the first term is always positive while h7 and not zerop h6 do << h8:=car h7; h7:=cdr h7; h9:=car h8; h8:=cdr h8; if contrace then <<write if nega then "--" else "++", " h8=",h8," h9=",h9;terpri()>>; if contrace then <<write"h9-2=",h9;terpri()>>; if h9 then h6:=totdif(h6,nth(xlist,h9),h9,dulist); if contrace then <<write"h6-3=",h6;terpri()>>; h10:=if h8 then mkid(u,h8) else u; if contrace then <<write"h10=",h10;terpri()>>; h10:=list('TIMES,h6,h10); if nega then h10:=list('MINUS,h10); if contrace then algebraic write">>>> h10=",h10; pii:=cons(h10,pii)$ nega:=not nega; >> >>; % for each function u il:=newil(il,mo,nx); >> until null il; pii:=reval if null pii then 0 else if length pii=1 then car pii else cons('PLUS,pii); if contrace then algebraic write"pii-end=",pii; pii >>; % for all ii return cons('LIST,pl) end$ % of intcurrent1 %------------- symbolic operator intcurrent2$ % used in conlaw2,4, crdec symbolic procedure intcurrent2(divp,ulist,xlist)$ % computes the conserved current P_i from the divergence through % partial integration % potential improvement: one could substitute low order derivatives % by high order derivatives using remaining conditions and try again begin scalar h2,h3,h4,h5,h6,h7,h8,e2,e3; % repeated partial integration to compute P_i ulist :=cdr reval ulist; xlist :=cdr reval xlist; h4:=list xlist; % dequ is here a list containing only the conserved density % and flux to drop commen factors repeat << e3:=divp; h3:=car h4; % h3 is the list of variables is a spec. order h4:=cdr h4; h5:=for e2:=1:length h3 collect 0; % h5 is old list of the conserved quantities h8:=0; % h8 counts integrations wrt. all variables repeat << % integrate several times wrt. all variables h8:=h8+1; h2:=h3; % h2 is a copy of h3 h6:=nil; % h6 is new list of the conserved quantities h7:=nil; % h7 is true if any integration was possible while h2 neq nil do << % integrating wrt. each variable e2:=intpde(e3,ulist,h3,car h2,t); if null e2 then e2:=list(nil,e3) else e3:=cadr e2; if (car e2) and (not zerop car e2) then h7:=t; h6:=cons(list('PLUS,car e2,car h5),h6); h5:=cdr h5; h2:=cdr h2 >>; h5:=reverse h6; >> until (h7=nil) % no integration wrt. no variable was possible or (e3=0) % complete integration or (h8=10); % max. 10 integrations wrt. all variables >> until (e3=0) or (h4=nil); return {'LIST,reval cons('LIST, h5),e3} % end of the computation of the conserved current P % result is a {{P_i},remainder} % was successful if 0 = remainder (=e3) end$ % of intcurrent2 %------------- symbolic operator intcurrent3$ % crident symbolic procedure intcurrent3(divp,ulist,xlist)$ % computes the conserved current P_i from the divergence through % partial integration with restriction of maximal 2 terms begin scalar xl,h1,h2,h3,h4,h5,resu1,resu2,resu3,succ; % repeated partial integration to compute P_i ulist :=cdr reval ulist; xlist :=cdr reval xlist; xl:=xlist; resu1:=nil; succ:=nil; % try all possible different pairs of variables while (cdr xl) and not succ do << h1:=intpde(divp,ulist,xlist,car xl,t); if h1 and not zerop car h1 then << resu2:=cons(car h1,resu1); h2:=cdr xl; repeat << h3:=intpde(cadr h1,ulist,xlist,car h2,nil); if h3 and zerop cadr h3 then << h4:=cons(car h3,resu2); for each h5 in cdr h2 do h4:=cons(0,h4); succ:=t; resu3:= {'LIST,cons('LIST,reverse h4),0} >>; h2:=cdr h2; resu2:=cons(0,resu2) >> until succ or null h2; >>; resu1:=cons(0,resu1); xl:=cdr xl >>$ return if succ then resu3 else {'LIST,cons('LIST,cons(0,resu1)),divp} end$ % of intcurrent3 endmodule$ %******************************************************************** module quasilinpde_integration$ %******************************************************************** % Routines to solve a quasilinear first order PDE % Author: Thomas Wolf % summer 1995 %---------------------------- algebraic procedure select_indep_var(clist,xlist)$ begin scalar s,scal,notfound,cq,x,xx,xanz,sanz,xcop,ok; % Find the simplest non-zero element of clist notfound:=t; xcop:=xlist; while notfound and (xcop neq {}) do << x :=first xcop ; xcop :=rest xcop ; cq:=first clist; clist:=rest clist; if cq neq 0 then << xanz:=0; for each xx in xlist do if %df(cq,xx) not freeof(cq,xx) then xanz:=xanz+1; ok:=nil; if not s then ok:=t else % to guarantee s neq nil if xanz<sanz then ok:=t else % as few dependencies as poss. if xanz=sanz then if length(cq)=1 then ok:=t else if length(scal)>1 then if part(scal,0)=PLUS then % if possible not a sum if (part(cq,0) neq PLUS) or (length(cq)<length(scal)) then ok:=t; if ok then <<s:=x; scal:=cq; sanz:=xanz>>; if scal=1 then notfound:=nil >> >>; return {s,scal} end$ % of select_indep_var %---------------------------- algebraic procedure charsyscrack(xlist,cqlist,s,scal,u,ode)$ % formulation and solution of the characteristic ODE-system begin scalar lcopy,x,cS,flist,soln,printold,timeold,facintold, adjustold,safeintold,freeintold,odesolveold, e1,e2,e3,e4,n,nmax,dff,proclistold; % formulating the ode/pde - system . lcopy := cqlist ; cS := {} ; for each x in xlist do << cq := first lcopy ; lcopy := rest lcopy ; if x neq s then << depend x,s; cS := .(scal*df(x,s)-cq,cS) >> >>; if s neq u then <<flist:={u}; depend u,s; cS := .(scal*df(u,s)+ode,cS) >> else flist:={}; for each x in xlist do if s neq x then flist:=.(x,flist); lisp << timeold := time_; time_ :=nil; facintold:=facint_; facint_:=1000; adjustold:=adjust_fnc; adjust_fnc:=t; safeintold:=safeint_; safeint_:=nil; freeintold:=freeint_; freeint_:=nil; odesolveold:=odesolve_; odesolve_:=50; proclistold:=proc_list_; proc_list_:=delete('alg_solve_single,proc_list_) >>$ % solving the odes using crack. if lisp(print_ neq nil) then lisp << write "The equivalent characteristic system:";terpri(); deprint cdr algebraic cS$ %terpri()$ write "for the functions: "; fctprint( cdr reval algebraic flist);write"."; >>; soln:=crack(cS,flist,flist,{}); lcopy:={}; for each x in soln do << e1:=first x; if (e1 neq {}) and (length e1 + length second x = length flist) then << % all remaining conditions are algebraic (non-differential) in all the % functions of flist? e2:={}; e3:={}; for each e4 in third x do if freeof(flist,e4) then e3:=e4 . e3 else e2:=e4 . e2; if (length cS) = (length e3) then << % sufficiently many integrations done for each e4 in e2 do for each e5 in e1 do if lisp(not freeof(lderiv(reval algebraic e5, reval algebraic e4, list(reval algebraic s)),'DF)) then <<e2:={};e1:={}>>; if e2 neq {} then << % It may be possible that derivatives of unsolved functions % occur in the expressions of the evaluated functions: second x nmax:=0; for each e4 in e2 do << % for each unsolved function for each e5 in second x do << % for each solved expression lisp << n:=lderiv(reval algebraic rhs e5,reval algebraic e4, list(reval algebraic s)); n:=if (car n) = nil then 0 else if (length car n) = 3 then 1 else cadddr car n >>; n:=lisp n; if n>nmax then nmax:=n; >>; if nmax>0 then << % adding further conditions e5:=e1; while freeof(first e5,e4) do e5:=rest e5; e5:=first e5; dff:=e4; for n:=1:nmax do << e5 :=df(e5,s); e1:= e5 . e1; dff:=df(dff,s); e3:=dff . e3 >> >> >>; lcopy:=cons({append(second x,e1),e3},lcopy); >> >> >> else if (first x = {}) and (length cS = length third x) then lcopy:=cons({second x,third x},lcopy) >>; lisp << time_:=timeold; facint_:=facintold; adjust_fnc:=adjustold; safeint_:=safeintold; freeint_:=freeintold; odesolve_:=odesolveold; proc_list_:=proclistold >>; return if lcopy={} then <<for each x in flist do if not my_freeof(x,s) then nodepend x,s; {}>> else s . lcopy % { {{x=..,y=..,u=..,..,0=..},{df(z,s),df(z,s,2),..,c1,c2,c3,..}},..} % df(z,s,..) only if df(z,s,..) turns up in x, y, u. .. . end$ % of charsyscrack %---------------------------- procedure charsyspsfi(xlist,cqlist,u,ode,denf); begin scalar h,s; h:=cqlist; cqlist:={}; while h neq {} do <<cqlist:=cons(first h*denf,cqlist);h:=rest h>>; cqlist:=cons(-ode*denf,cqlist); xlist:=cons(u,reverse xlist); s:=lisp gensym(); for each h in xlist do depend h,s; h:=psfi(cqlist,xlist,s); for each h in xlist do if not my_freeof(h,s) then nodepend h,s; return h end$ % of charsyspsfi %---------------------------- algebraic procedure storedepend(xlist)$ % stores the dependencies of the elements of xlist in the returned % list and clears the dependencies begin scalar q,e1,e2$ return for each e1 in xlist collect << q:=fargs e1; for each e2 in q do nodepend e1,e2; cons(e1,q)>> end$ % of storedepend %---------------------------- algebraic procedure restoredepend(prev_depend)$ % assigns the dependencies stored in the argument begin scalar q,s,x; for each s in prev_depend do <<q:=first s; s:=rest s; for each x in s do depend q,x>> end$ % of restoredepend %---------------------------- symbolic procedure simplifiziere(q,fl)$ begin scalar n; return if not pairp q then q else if member(car q,ONE_ARGUMENT_FUNCTIONS_) or (member(car q,{'EXPT,'QUOTIENT}) and not smemberl(fl,caddr q) ) then simplifiziere(cadr q,fl) else if car q = 'QUOTIENT and not smemberl(fl,cadr q) then simplifiziere(caddr q,fl) else << n:=ncontent(q); if n=1 then q else simplifiziere(reval list('QUOTIENT, q, reval n),fl) >> end$ %---------------------------- algebraic procedure quasilinpde(f,u,xlist); % f ... PDE, u ... unknown function, xlist ... variable list begin scalar i, q, qlist, cq, cqlist, ode, soln, tran, xcop, s, s1, x, xx, h1, h2, scal, qlin, prev_depend, tr_qlp,xlist_cp1,xlist_cp2; symbolic put('ff,'simpfn,'simpiden)$ tr_qlp:=t; if lisp print_ then write"The quasilinear PDE: 0 = ",f,"."; % changing the given pde into a quasi-linear ode . i := 0 ; ode := f; qlist := {}; cqlist:={}; for each x in xlist do <<i := i+1 ; q := mkid(p_,i) ; qlist := .(q,qlist) ; f := sub(df(u,x) = q , f) ; cq:=df(f,q); cqlist:=.(df(f,q),cqlist) ; ode := ode - df(u,x)*df(f,q) ; %if not my_freeof(u,x) then nodepend u, x ; >> ; lisp(depl!*:=delete(assoc(u,depl!*),depl!*))$ lisp(depl!*:=delete(assoc(mkid(u,'_),depl!*),depl!*))$ qlist := reverse qlist ; cqlist := reverse cqlist ; % checking for linearity. qlin:=t; for each cq in cqlist do for each q in qlist do if df(cq,q) neq 0 then qlin:=nil; if not qlin then return {}; % soln:=charsyspsfi(xlist,cqlist,u,ode,den f); % Determination of the independent variable for the ODEs pcopy:=cons(-ode,cqlist);xcop:=cons(u,xlist); scal:=select_indep_var(pcopy,xcop)$ s1:=first scal; prev_depend:=storedepend(xlist)$ soln:=charsyscrack(xlist,cqlist,s1,second scal,u,ode); if soln={} then << % try all other coefficients as factors repeat << repeat <<s :=first xcop ;xcop :=rest xcop; scal:=first pcopy;pcopy:=rest pcopy>> until (pcopy={}) or ((scal neq 0) and (s neq s1)); if (s neq s1) and (scal neq 0) then << if lisp print_ then lisp <<terpri()$ write"New attempt with a different independent variable:"; terpri()>>$ soln:=charsyscrack(xlist,cqlist,s,scal,u,ode) >> >> until (soln neq {}) or (xcop={}) >>; % solving for the constants(eg..c1,c2 etc.) and put them in a % linear form. % The first element of soln is the ODE-variable tran:={}; if soln neq {} then << s1:=first soln; for each s in rest soln do <<cq:=first solve(first s,second s); x:={}; for each xx in cq do if lisp atom algebraic lhs xx then x:=.(lisp <<q:=reval algebraic rhs xx; simplifiziere(q,cdr xlist)>>,x); lisp << x:=cdr x; xlist_cp1:=cdr xlist; xlist_cp2:=xlist_cp1; repeat << for each h1 in xlist_cp1 do if member(h1,x) then xlist_cp2:=delete(h1,xlist_cp2); if xlist_cp1 neq xlist_cp2 then << xlist_cp1:=xlist_cp2; x:=for each h1 in x collect simplifiziere(h1,xlist_cp1) >> >> until xlist_cp1=xlist_cp2; x:=cons('LIST,x); >>$ xx:=tran; while (xx neq {}) and <<h1:=x;h2:=first xx; while (h1 neq {}) and (h2 neq {}) and (first h1 = first h2) do <<h1:=rest h1; h2:=rest h2>>; if (h1={}) and (h2={}) then nil else t >> do xx:=rest xx; if xx={} then tran:=.(x,tran); >>; for each s in xlist do if (s neq s1) and (not my_freeof(s,s1)) then nodepend s,s1 >>; for each x in xlist do depend u,x; if lisp print_ then if tran neq {} then << write"The general solution of the PDE is given through"; for each x in tran do write"0 = ", lisp( cons('ff,cdr reval algebraic x)); if length tran>1 then write"with arbitrary function(s) ff(..)." else write"with arbitrary function ff(..)." >>; % restoring the dependencies of the variables of the PDE restoredepend(prev_depend)$ return tran; end$ % of quasilinpde endmodule$ end$