Artifact 13c3017f9bb044cab86b300c396a846ecb1a3041a53a03e1418d0622e4ec64cb:
- Executable file
r38/packages/crack/crunder.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: 14866) [annotate] [blame] [check-ins using] [more...]
%******************************************************************** module underdetde$ %******************************************************************** % Routines for the solution of underdetermiined ODEs and PDEs % Author: Thomas Wolf % August 1998, February 1999 symbolic procedure undetlinode(arglist)$ % parametric solution of underdetermined ODEs begin scalar l,l1,p,pdes,forg,s$ pdes:=car arglist$ forg:=cadr arglist$ if expert_mode then <<l1:=selectpdes(pdes,1);flag(l1,'to_under)>> else l1:=cadddr arglist$ while l1 do if null (p:=get_ulode(l1)) then l1:=nil else << l:=underode(p,pdes)$ p:=car p$ if null l then <<remflag1(p,'to_under);l1:=delete(p,l1)>> else << if print_ then << write"Parametric solution of the underdetermined ODE ",p$ terpri(); write"giving the new ODEs "$ s:=l; while s do << write car s; s:=cdr s; if s then write "," >>$ terpri() >>$ pdes:=drop_pde(p,pdes,nil)$ for each s in l do pdes:=eqinsert(s,pdes)$ l:=list(pdes,forg)$ l1:=nil; >> >>$ return l$ end$ symbolic procedure undetlinpde(arglist)$ % parametric solution of underdetermined PDEs begin scalar l,l1,p,pdes,forg$ pdes:=car arglist$ forg:=cadr arglist$ if expert_mode then <<l1:=selectpdes(pdes,1);flag(l1,'to_under)>> else l1:=cadddr arglist$ while l1 do if null (p:=get_ulpde(l1)) then l1:=nil else << l:=underpde(p,pdes)$ % l has to be the list of new pdes p:=car p$ if null l then <<remflag1(p,'to_under);l1:=delete(p,l1)>> else << pdes:=drop_pde(p,pdes,nil)$ for each s in l do pdes:=eqinsert(s,pdes)$ l:=list(pdes,forg)$ l1:=nil; >> >>$ return l$ end$ symbolic procedure get_ulode(pdes)$ begin scalar h,best_ulode; for each p in pdes do if flagp(p,'to_under) then if null (h:=ulodep(p)) then remflag1(p,'to_under) else if ((null best_ulode) or (car h < car best_ulode)) then best_ulode:=h; return if best_ulode then cdr best_ulode else nil end$ symbolic procedure get_ulpde(pdes)$ begin scalar h,best_ulpde; for each p in pdes do if flagp(p,'to_under) then if null (h:=ulpdep(p)) then remflag1(p,'to_under) else if ((null best_ulpde) or (car h < car best_ulpde)) then best_ulpde:=h; return if best_ulpde then cdr best_ulpde else nil end$ symbolic procedure ulodep(p)$ begin scalar tr_ulode,drvs,ulode,allvarf,minord,dv,f,x,h,found,minordf,totalorder$ % minord and minordf are currently not needed later on % tr_ulode:=t; % Is it an underdetermined linear ODE for the allvarfcts? drvs:=get(p,'derivs)$ ulode:=t$ allvarf:=get(p,'allvarfcts); if tr_ulode then << write"allvarf=",allvarf$ terpri()$ >>$ minord:=1000; if not (allvarf and cdr allvarf) then ulode:=nil else % at least two allvar-fcts while ulode and drvs do << dv:=caar drvs; f:=car dv$ if tr_ulode then << write"car drvs=",car drvs," dv=",dv," f=",f, " member(f,allvarf)=",member(f,allvarf)$ terpri()$ >>$ if member(f,allvarf) then if (cdar drvs neq 1) or % not linear % x is already determined and it is not x: (cdr dv and ((x and (x neq cadr dv) ) or % there are other differentiation variables: ((cddr dv) and ((not fixp caddr dv) or cdddr dv )) ) ) then << % no ODE: ulode:=nil; if tr_ulode then << write"new ulode=",ulode$ terpri()$ >>$ >> else % can be an ODE if null cdr dv then % f has no derivatives if not member(f,found) then ulode:=nil % no ODE --> substitition else % f has already been found with a % consequently higher x-derivative else % this is an x-derivative of f if null x then << % x had not yet been determined if tr_ulode then << write"null x"$ terpri()$ >>$ found:=cons(f,found)$ x:=cadr dv; minordf:=car dv; if null cddr dv then minord:=1 else minord:=caddr dv; totalorder:=minord >> else % x has already been determined if not member(f,found) then << % only leading derivatives matter found:=cons(f,found)$ % and the first deriv. of f is leading if null cddr dv then h:=1 else h:=caddr dv; totalorder:=totalorder+h; if h<minord then << minord:=h; minordf:=car dv >> >> else else % not member(f,allvarf) if null x or % then there are only derivatives % of non-allvarfcts left member(x,fctargs f) then ulode:=nil; % no x-dependent non-allvarfcts if tr_ulode then << write"found=",found," minord=",minord," minordf=",minordf$ terpri()$ >>$ drvs:=cdr drvs; >>$ if tr_ulode then << write"ulode=",ulode$ terpri()$ >>$ return if ulode then {totalorder,p,x,minord,minordf} else nil end$ % of ulodep symbolic procedure ulpdep_(p)$ begin scalar tr_ulpde,drvs,drv,ulpde,allvarf,allvarfcop, vld,vl,v,pde,fn,f,no_of_drvs,no_of_tms,ordr,maxordr,parti$ %tr_ulpde:=t; % Is it an underdetermined linear PDE for the allvarfcts? drvs:=get(p,'derivs)$ ulpde:=t$ allvarf:=get(p,'allvarfcts); if tr_ulpde then << write"allvarf=",allvarf$ terpri()$ >>$ if not (allvarf and cdr allvarf) then ulpde:=nil else << % at least two allvar-fcts allvarfcop:=allvarf$ no_of_tms:=0; % total number of terms of all diff operators vl:=get(p,'vars)$ while ulpde and allvarfcop do << % extracting the differential operator for car allvarfcop pde:=get(p,'val); fn:=car allvarfcop; allvarfcop:=cdr allvarfcop; for each f in allvarf do if f neq fn then pde:=subst(0,f,pde); pde:=reval pde; % Is pde linear in fn? if not lin_check(pde,{fn}) then << if tr_ulpde then <<write"not linear in ",f$terpri()>>$ ulpde:=nil >> else << % add up the number of terms no_of_tms:=no_of_tms + no_of_terms(pde)$ % What are all variables in pde? This is needed to test later % whether they are disjunct from all variables from another % diff. operator vld:=nil; for each v in vl do if not freeof(pde,v) then vld:=cons(v,vld); % What is the number of derivatives of fn? % What order is the highest derivative of fn? no_of_drvs:=0; for each drv in drvs do if fn=caar drv then << ordr:=absdeg(cdar drv); if (no_of_drvs=0) or (ordr>maxordr) then maxordr:=ordr; no_of_drvs:=add1 no_of_drvs; >>; % collect the differential operators with properties in parti parti:=cons({pde,fn,vld,no_of_drvs,maxordr},parti); >> >>; if no_of_tms neq get(p,'terms) then << if tr_ulpde then << write"not a lin. homog. PDE"$terpri() >>$ ulpde:=nil; % not a lin. homog. PDE >>$ if tr_ulpde then << write"parti="$ prettyprint parti$ >>$ >>$ return if null ulpde then nil else parti end$ symbolic procedure ulpdep(p)$ begin scalar tr_ulpde,drvs,drv,ulpde,parti,pde, difop1,difop2,commu,disjun,totalcost$ %tr_ulpde:=t; % Is it an underdetermined linear PDE for the allvarfcts? drvs:=get(p,'derivs)$ ulpde:=ulpdep_(p)$ if ulpde then << parti:=ulpde$ ulpde:=t$ % Find a differential operator pde in parti pde:=nil; for each difop1 in parti do << commu:=t; % which does commute with all other diff. operators for each difop2 in parti do if (cadr difop1 neq cadr difop2) and (not zerop reval {'DIFFERENCE,subst(car difop1,cadr difop2,car difop2), subst(cadr difop1,cadr difop2, subst(car difop2,cadr difop1,car difop1))}) then << commu:=nil; if tr_ulpde then << write"no commutation of:"$terpri()$ prettyprint difop1$ write"and "$terpri()$ prettyprint difop2 >> >>$ % and is variable-disjunct with at least one other diff. operator if commu then << disjun:=nil; for each difop2 in parti do if (cadr difop1 neq cadr difop2) and freeoflist(caddr difop1,caddr difop2) then disjun:=t; if disjun then if null pde then pde:=difop1 else if ( car cddddr difop1) < (car cddddr pde) or % minimal maxord (((car cddddr difop1) = (car cddddr pde)) and % minimal number of terms ((cadddr difop1) < (cadddr pde)) ) then pde:=difop1 >> >>; if null pde then ulpde:=nil >>; if tr_ulpde then << write"ulpde=",ulpde$ terpri()$ >>$ % as a first try we take as cost for the PDE p the sum of all orders % of all derivatives of all functions in p totalcost:=0; for each drv in drvs do totalcost:=totalcost+absdeg(cdar drv); return if ulpde then {totalcost,p,pde,parti} else nil end$ % of ulpdep symbolic procedure min_ord_f(ode,allvarf,vl)$ begin scalar minord,minordf,newallvarf,f,h,tr_ulode$ % tr_ulode:=t; % does any function not occur anymore? % Which function does occur with lowest derivative: minord, minordf? minord:=1000; minordf:=nil; newallvarf:=nil; for each f in allvarf do << h:=ld_deriv_search(ode,f,vl)$ if tr_ulode then << write"ld_deriv_search(",f,")=",h$ terpri()$ >>$ if not zerop cdr h then << % otherwise f does not occur anymore in ode newallvarf:=cons(f,newallvarf)$ h:=car h$ h:=if null h then 0 else if null cdr h then 1 else cadr h$ % asuming that car h = x if h<minord then <<minord:=h;minordf:=f>> >> >>$ return {minord,minordf,newallvarf} end$ symbolic procedure underode(pchar,pdes)$ % pchar has the form {p,x,minord,minordf} begin scalar tr_ulode,p,x,allvarf,orgallvarf,ode,noallvarf,vl, minord,minordf,adj,f,h,newf,sol,sublist,rtnlist$ % tr_ulode:=t; p :=car pchar; x :=cadr pchar; minord :=caddr pchar; minordf:=cadddr pchar; allvarf:=get(p,'allvarfcts); orgallvarf:=allvarf; ode:=get(p,'val); noallvarf:=length(allvarf); vl:=get(p,'vars); while (minord>0) and (length(allvarf)=noallvarf) do << if tr_ulode then << write "x=",x," minord=",minord," minordf=",minordf, " allvarf=", allvarf$ terpri()$ >>$ repeat << adj:=intpde(ode,allvarf,vl,x,t); if tr_ulode then << write"adj=",adj$ terpri()$ >>$ h:=nil; for each f in allvarf do if not freeof(cadr adj,f) then h:=cons(f,h); if null h then % exact ode --> should do better then what is done now ode:=reval {'TIMES,x,ode}; >> until h; minordf:=cadr min_ord_f(ode,h,vl)$ % a new function (potential) is needed: newf:=newfct(fname_,vl,nfct_)$ nfct_:=add1 nfct_; if tr_ulode then << algebraic write"eqn=",{'LIST,{'PLUS,{'DF,newf,x},lisp cadr adj}}$ algebraic write"var=",{'LIST,minordf } >>$ sol:=cadr solveeval list({'LIST,{'PLUS,{'DF,newf,x},cadr adj}}, {'LIST,minordf } ); allvarf:=delete(minordf,allvarf)$ allvarf:=cons(newf,allvarf)$ % assuming that there is exacly one solution to the lin. alg. equation if tr_ulode then << terpri()$ write"sol=",sol$ terpri()$ >>$ sublist:=cons(sol,sublist)$ ode:=reval num reval {'PLUS,newf,{'MINUS,subst(caddr sol,cadr sol,car adj)}}$ if tr_ulode then algebraic(write"ode=",ode)$ h:=min_ord_f(ode,allvarf,vl)$ minord:=car h; minordf:=cadr h; allvarf:=caddr h; if minord=0 then sublist:=cons(cadr solveeval list({'LIST,ode},{'LIST,minordf}),sublist)$ if tr_ulode then << write"allvarf=",allvarf," minord=",minord," minordf=",minordf$ terpri()$ >>$ >>$ if (minord neq 0) and (not zerop ode) then rtnlist:=list ode; ode:=nil; if tr_ulode then <<write"rtnlist=", rtnlist;terpri()>>$ if tr_ulode then algebraic(write"sublist=",cons('LIST,sublist)); while sublist do << if member(cadar sublist,orgallvarf) then rtnlist:=cons(reval num reval {'PLUS,cadar sublist, {'MINUS,caddar sublist}},rtnlist)$ sublist:=cdr reval cons('LIST, subst(caddar sublist,cadar sublist,cdr sublist))$ if tr_ulode then algebraic(write"sublist=",cons('LIST,sublist)) >>$ allvarf:=smemberl(allvarf,rtnlist)$ if tr_ulode then << write"allvarf=",allvarf$ terpri()$ >>$ for each h in allvarf do ftem_:=fctinsert(h,ftem_)$ if tr_ulode then algebraic(write"rtnlist=",cons('LIST,rtnlist)); h:=for each h in rtnlist collect mkeq(h,union(get(p,'fcts),allvarf),vl,allflags_,t,list(0),nil,pdes)$ if print_ then terpri()$ return h end$ symbolic procedure underpde(pchar,pdes)$ % pchar has the form {p,difop,parti} where p is the name of the equation, % difop is the main differential operator in p and parti is a partition % of p, i.e. the list of all differential operators. Each differential % operator has the form % {pde,fn,vld,no_of_drvs,maxordr} % where pde are all terms in p with fn, vld is a list of all variables % in pde, no_of_dervs is the number of different derivatives of fn, % maxord is the maximal order of derivatives of fn (order of pde) begin scalar tr_ulpde,ldo,parti,fn,lcond,difop,h,fl,eqlist,vl$ % has to return list of new pde just like in underode % tr_ulpde:=t; ldo:=cadr pchar; parti:=caddr pchar$ vl:=get(car pchar,'vars)$ fn:=cadr ldo$ lcond:={fn}$ if tr_ulpde then << write"ldo="$prettyprint parti$ write"parti="$prettyprint parti >>$ for each difop in parti do if cadr difop neq fn then << h:=newfct(fname_,vl,nfct_)$ nfct_:=add1 nfct_$ if print_ then terpri()$ fl:=cons(h,fl)$ eqlist:=cons(cons({cadr difop,h}, reval {'DIFFERENCE,cadr difop,subst(h,fn,car ldo)}), eqlist)$ lcond:=cons(subst(h,cadr difop,car difop),lcond) >>$ eqlist:=cons(cons(append(get(car pchar,'fcts),fl), cons('PLUS,lcond)),eqlist)$ if tr_ulpde then << write"eqlist="$prettyprint eqlist$ >>$ h:=for each h in eqlist collect mkeq(cdr h,car h,get(car pchar,'vars),allflags_,t,list(0),nil,pdes)$ if print_ then terpri()$ return h end$ endmodule$ end$