with(Groebner,normalf); with(Groebner,gsolve); with(Groebner,inter_reduce); with(linalg,jacobian);with(linalg,transpose);with(linalg,inverse); with(Groebner,gbasis);with(Groebner,termorder);with(Ore_algebra); ######################################################## # Maple package Indiff # # written by Elizabeth Mansfield # # University of Kent at Canterbury # # Copyright 1999 # # Manual and Information available at # #http://www.ukc.ac.uk/ims/maths/people/E.L.Mansfield.html# ########################################################$$ Buchberger2:=proc(pairnow,pairsdone,GG,allvars,termorder) local g1,g2,h,hset1,hset2,hset,hdt,hp,hc,A,Z,j; #readlib('LCD',``.`indiff`.`/`.`code`.`/`.`LCD.m`); #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`); #readlib('Idiffdivide',``.`indiff`.`/`.`code`.`/`.`Idiffdivide.m`); #we assume that pairnow is a set of 2 different integers #we assume that pairsdone is a set of sets of pairs of integers g1:=op(op(1,pairnow),GG); g2:=op(op(2,pairnow),GG); hset1:=select(has,pairsdone,op(1,pairnow)); hset2:=select(has,pairsdone,op(2,pairnow)); hset1:=map(proc(m) op(subs(op(1,pairnow)=NULL,m)) end,hset1); hset2:=map(proc(m) op(subs(op(2,pairnow)=NULL,m)) end,hset2); hset:=hset1 intersect hset2; Z:=false; for j in hset do Idiffparse2(op(j,GG),allvars,termorder,'hdt','hp','hc'); if LCD(g1,g2,allvars,termorder)<>false then Z:=Idiffdivide(hdt,LCD(g1,g2,allvars,termorder),'A'); if Z=true then break fi; fi; od; Z; end: DNAlgGB:=proc(F,allvars,mytermorder,newF) #F a list or set of idp's #allvars list of list of dpendent variables #newF name of output local A,C,N1,DTlist,m,T,n,Hmonset,Hmonset2,mat,fred,sue,sentinel,Ktmp,Kdenoms,\ j,NZ,A2,allukns; global vars,DegDiffTable; #with(Groebner,gbasis);with(Groebner,termorder);with(Ore_algebra); #with(Groebner,inter_reduce); #readlib('Idtsort',``.`indiff`.`/`.`code`.`/`.`Idtsort.m`); #readlib('degofdiff',``.`indiff`.`/`.`code`.`/`.`degofdiff.m`); #readlib('IHmon',``.`indiff`.`/`.`code`.`/`.`IHmon.m`); #readlib('sortanddiffle',``.`indiff`.`/`.`code`.`/`.`sortanddiffle.m`); #readlib('makeDTlist',``.`indiff`.`/`.`code`.`/`.`makeDTlist.m`); #readlib('makeorder',``.`indiff`.`/`.`code`.`/`.`makeorder.m`); #readlib('noKfactors',``.`indiff`.`/`.`code`.`/`.`noKfactors.m`); A:=convert(numer(F),set); if nops(A)<2 then newF:=A; RETURN(A) fi; C:=A;#initialize C #we need to collect factors of the denominators of the Kmatrix Ktmp:=convert(Kmatrix,listlist): Ktmp:={seq(op(op(j,Ktmp)),j=1..nops(Ktmp))}: Kdenoms:=select(has,map(denom,Ktmp),`In`); Kdenoms:=map(m-> if type(factor(m),`*`) then op(factor(m)) elif type(factor(m),`^`) then op(1,factor(m)) else m fi,Kdenoms); Kdenoms:=map(m-> if type(factor(m),`^`) then op(1,factor(m)) else m fi,Kdenoms); Kdenoms:=select(has,Kdenoms,`In`); #we addd in any extra we want to be nonzero if nargs>5 then NZ:=args[6]; else NZ:={} fi; sentinel:=1;allukns:=[seq(op(op(j,allvars)),j=2..nops(allvars))]; while C<>{} do #degofdiff(C,'N1','m'); #if N1 > N+1 then ERROR(`degree greater than`, N+1, `detected`) fi; sortanddiffle(C,allvars,mytermorder,'C',Kdenoms union NZ); if sentinel=1 then Hmonset:=map(IHmon,C,allvars,ttdeg,'m') fi; sentinel:=0; makeDTlist( {op(A),op(C)}, vars,allukns,'DTlist'); # organise an order for the algebraic GBasis calculation makeorder(DTlist,allvars,'T','mat'); fred:=poly_algebra(op(T)); sue:=termorder(fred,'matrix'(convert(mat,listlist),T)); A:=gbasis( {op(A),op(C)} ,sue); #we throw away factors which are in Kdenoms as these are non-zero #and interfere with termination of the algorithm #we also throw away generators which are raised to whole powers #as this is of NO interest, also any factors in NonZero A2:=noKfactors(A,Kdenoms union NZ); while A<>A2 do A:=inter_reduce(A2,sue);A2:=map(expand,noKfactors(A,Kdenoms union NZ)); od; C:={}; #pull out the new ones for n in A do degofdiff(n,allvars,mytermorder,'N1'); if op(2,N1)>DegDiffTable[op(1,N1)] then next fi; IHmon(n,allvars,ttdeg,'m'); if not member(m,Hmonset) then C:=C union {n}; Hmonset:=Hmonset union {m} fi; od; od; #need to be able to export the termordering for the simplification of Kmatrix if nargs>4 then assign(args[5],sue) fi; newF:=A; end: noKfactors:=proc(A,Kdenoms) local Atmp,m,j,k,M; Atmp:=A; for j to nops(Atmp) do m:=factor(op(j,Atmp)); if type(m,`^`) then Atmp:=subsop(j=op(1,m),Atmp); next; fi; if not type(m,`*`) then next fi; M:=map(proc(k) local k2; if type(k,`^`) then k2:=op(1,k) else k2:=k fi; if member(k2,Kdenoms) then 1 else k2 fi end,[op(m)]); M:=convert(M,`*`); Atmp:=subsop(j=M,Atmp); od; Atmp; end: #computes the correction term in D_j #computes the correction term in D_j Error:=proc(j,u) local m,ell,Ind; global GroupP,vars,ukns,Kmatrix,XiPhis,Neqs; option system; #readlib('PROL',``.`indiff`.`/`.`code`.`/`.`PROL.m`): #readlib('Invariantize',``.`indiff`.`/`.`code`.`/`.`Invariantize.m`): if nargs>2 then Ind:=args[3] fi; if member(u,ukns) then expand(-convert([seq(Kmatrix[j,m]*PROL(u,op(m,GroupP),Ind), m=1..nops(GroupP))],`+`)) elif member(u,vars,'ell') then expand(-convert([seq(Kmatrix[j,m]*Invariantize(XiPhis[m,ell]), m=1..nops(GroupP))],`+`)) else 0 fi; end: #save Error, ``.`indiff`.`/`.`code`.`/`.`Error.m`; HNI:=proc(allvars,termorder) global HNIlist,HNIdifflist; #readlib('HNIinfo',``.`indiff`.`/`.`code`.`/`.`HNIinfo.m`); HNIinfo(allvars,termorder,'HNIlist','DegDiffTable'); end: HNIinfo:=proc(allvars,termorder,HNIlist,HNIdifflist) local j,k,HDT,Hp,Hc,UNeqs,allukns,m,difflisttmp; global Neqs; #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`); allukns:=[seq(op(op(j,allvars)),j=2..nops(allvars))]; UNeqs:=select(has,Neqs,allukns); for k to nops(UNeqs) do Idiffparse2(op(k,UNeqs),allvars,termorder,'HDT[k]','Hp','Hc') od; difflisttmp:=table([seq(allukns[j]=-2,j=1..nops(allukns))]); for j to nops(UNeqs) do if nops(HDT[j])=2 then m:=nops(op(2,HDT[j])) elif nops(HDT[j])=3 then m:=convert(op(3,HDT[j]),`+`) else -2 fi; if difflisttmp[op(1,HDT[j])]{} then Idiffparse2(g,allvars,termorder,'HDT','Hp','HC'); temp:=HDT^Hp; while (has(HC,'In') and type(HC,`+`)=true) do Idiffparse2(HC,allvars,termorder,'A','B','HC'); temp:=temp*A^B; od; #we're down now to polys in the variables if type(HC,`+`)=true then M:=convert(HC,set); Top:=op(1,M); for n to nops(M) do if Idifforder(op(n,M),1,Top,1,allvars,termorder)=1 then Top:=op(n,M) fi; od; temp:=temp*Top; else temp:=temp*HC fi; #this case f has no deriv terms elif type(g,`+`) then M:=convert(g,set); Top:=op(1,M); for n to nops(M) do if Idifforder(op(n,M),1,Top,1,allvars,termorder)=1 then Top:=op(n,M) fi; od; temp:=Top; else temp:=g fi; Hmon:=temp; end: IKolRitt:=proc(F,allvars,myorder,newF) local Ftmp,A,fdeg,matorder,sentinel,probset,DTlist,j,Au,B,Bu0,Done,Done0,\ uNeqs,Xset1,Xset2,lXset,lfactor,tt,s,t,NZ,T,B2,A2,allukns,mat,\ Km1,matorder2,C,m,r,rC; global vars,Neqs,Kmatrix,DegDiffTable; #with(Groebner,normalf);with(Groebner,gbasis);with(Groebner,termorder); #with(Groebner,inter_reduce); #with(Ore_algebra); #readlib('degofdiff',``.`indiff`.`/`.`code`.`/`.`degofdiff.m`); #readlib('IKolRitt2',``.`indiff`.`/`.`code`.`/`.`IKolRitt2.m`); #readlib('Ireduce',``.`indiff`.`/`.`code`.`/`.`Ireduce.m`); #readlib('Ieqnsort',``.`indiff`.`/`.`code`.`/`.`Ieqnsort.m`); #readlib('IHmon',``.`indiff`.`/`.`code`.`/`.`IHmon.m`); #readlib('Iorthreduceall',``.`indiff`.`/`.`code`.`/`.`Iorthreduceall.m`); #readlib('sortanddiffle',``.`indiff`.`/`.`code`.`/`.`sortanddiffle.m`); #readlib('makeorder',``.`indiff`.`/`.`code`.`/`.`makeorder.m`); #readlib('makeDTlist',``.`indiff`.`/`.`code`.`/`.`makeDTlist.m`); #readlib('DNAlgGB',``.`indiff`.`/`.`code`.`/`.`DNAlgGB.m`); #we only need here to sort out the Xset,factor options, the rest are handled by #IKolRitt2 lXset:=false;lfactor:=false;NZ:={}; if nargs>4 then for tt in [args[5..nargs]] do if not type(tt,equation) then ERROR(`invalid option argument`,tt) fi; s:=op(1,tt);t:=op(2,tt); if s='info' then if member('Xset',t) then lXset:=true fi; fi; if s='NonZero' then NZ:=t fi; if s='strategy' then if member('factor',t) then lfactor:=true; #with(Groebner,gsolve); fi; fi; od; fi; Ftmp:=convert(F,set); #removes multiple copies probset:=Ftmp; #initialize the probset Done:={}; allukns:=[seq(op(op(j,allvars)),j=2..nops(allvars))]; Xset2:={}; A:=map(proc(m) global DegDiffTable; degofdiff(m,allvars,myorder,'fdeg'); if op(2,fdeg)>DegDiffTable[op(1,fdeg)]+1 then NULL else m fi; end,Ftmp); B:=Ftmp minus A; #Hmonset:=map(proc(m) global DegDiffTable; degofdiff(m,allvars,myorder,'fdeg'); #if op(2,fdeg)<>DegDiffTable[op(1,fdeg)]+1 then NULL else #IHmon(m,allvars,myorder,'k0') fi; end, A); while probset<>{} do DNAlgGB(A,allvars,myorder,'A','matorder',NZ); if not has(A,allukns) then ERROR(`system inconsistent`,A) fi; gc(); #sets the remember table for Error to NULL #next we simplify the Kmatrix so that derivatives get simpler we hope #this means that the normalisation equations are no longer inbuilt #so we keep track of whether Kmatrix has altered or not for later Km1:=op(Kmatrix); Kmatrix:=map(proc(m) local h,hd,hn; h:=normal(m); r:=traperror(normalf(numer(h),A,matorder)); if r=`polynomials not members of the algebra` or nops(A)<2 then hn:=simplify(numer(h),A) else hn:=r fi; r:=traperror(normalf(denom(h),A,matorder)); if r=`polynomials not members of the algebra` or nops(A)<2 then hd:=simplify(denom(h),A) else hd:=r fi; if hd=0 then ERROR(`denominator in Kmatrix`,denom(h), `zero on input system`) fi; hn/hd; end, op(Kmatrix)); #it can happen theat where there is more than one unknown, that A #contains idp's whose degree is higher than N+1 So we collect them #to include them in B later Au:=map(proc(m) global DegDiffTable; degofdiff(m,allvars,myorder,'fdeg'); if op(2,fdeg)<=DegDiffTable[op(1,fdeg)]+1 then NULL else m fi; end, {op(A)} ); #it can happen that Hcoefficients in A are known to be zero from things in B #when we have more than one unknown, so this next few steps gets rid of this #worrying possibility Iorthreduceall({op(B),op(A)},allvars,myorder,'B', 'algebra'=[{op(A)} minus Au,matorder]); B2:=map(proc(m) global DegDiffTable; degofdiff(m,allvars,myorder,'fdeg'); if op(2,fdeg)>DegDiffTable[op(1,fdeg)]+1 then NULL else m fi; end, {op(B)} minus {op(A)}); if B2<>{} then A:={op(A)} minus Au union B2; makeDTlist(A,vars,allukns,'DTlist'); makeorder(DTlist,allvars,'T','mat'); matorder:=termorder(poly_algebra(op(T)),'matrix'(convert(mat,listlist),T)); A:={op(inter_reduce(A,matorder))}; fi; Bu0:=map(proc(m) global DegDiffTable; degofdiff(m,allvars,myorder,'fdeg'); if op(2,fdeg)=DegDiffTable[op(1,fdeg)]+1 then NULL else m fi; end, {op(Done)} ); #once we have simplified the Kmatrix, the normalisation eqns are no #longer inbuilt, so we need to put them back in #using the sortanddiffle below uNeqs:=select(has,Neqs,allukns); sortanddiffle(uNeqs,allvars,myorder,'C'); rC:={seq(Iorthreduce(op(m,C),B,allvars,myorder,'r','algebra'=[A,matorder]),m=1..nops(C))}; probset:=probset union rC minus {0}; A:={op(A),op(probset), op(Done0)}; Done:=Done union {op(B)}; Xset2:=Xset2 union Xset1; probset:=probset union {op(B)} minus Done; B:= {op(B)} minus {op(probset)}; od; #outputting the Xset if lXset=true then Xset2:=map(proc(m) if type(m,numeric) then NULL else m fi; end,Xset2); lprint(`the Xset is generated by`,Xset2) fi; #the clean up operation B2:=map(proc(m) global DegDiffTable; degofdiff(m,allvars,myorder,'fdeg'); if op(2,fdeg)>=DegDiffTable[op(1,fdeg)]+1 then m else NULL fi; end, {op(A),op(B)}); #A2:={op(A)} union ({op(B)} minus B2); A2:= ({op(A),op(B)} minus B2); makeDTlist(A2, vars,allukns,'DTlist'); makeorder(DTlist,allvars,'T','matorder2'); matorder2:=termorder(poly_algebra(op(T)), 'matrix'(convert(matorder2,listlist),T)): A2:=inter_reduce(A2,matorder2); Iorthreduceall(B2,allvars,myorder,'B2','algebra'=[A2,matorder2]); Ieqnsort({op(A2),op(B2)},allvars,myorder,'newF'); end: IKolRitt2:=proc(G1,G2,allvars,myorder,A,matorder,F,probset,XXX) local n,t,s,tt,hdt,hp,hc,GG,pairsdone,i,j,h,top,hdtk,A2,\ hpk,hck,k,m,lshow,lfactor,otherF,allvars2,morestrat,X,XX,probset2,Kdenoms,\ Ktmp,newP,pairstodo,top2,pairnow,test,NZ,fdeg,mlist,melt;#,allukns; global vars,Kmatrix,DegDiffTable; #G2 are the new ones and G1 are the Done ones #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`); #readlib('IdSpoly',``.`indiff`.`/`.`code`.`/`.`IdSpoly.m`); #readlib('Ireduce',``.`indiff`.`/`.`code`.`/`.`Ireduce.m`); #readlib('Iorthreduceall',``.`indiff`.`/`.`code`.`/`.`Iorthreduceall.m`); #readlib('Ieqnsort',``.`indiff`.`/`.`code`.`/`.`Ieqnsort.m`); #readlib('degofdiff',``.`indiff`.`/`.`code`.`/`.`degofdiff.m`); #readlib('noKfactors',``.`indiff`.`/`.`code`.`/`.`noKfactors.m`); #readlib('Buchberger2',``.`indiff`.`/`.`code`.`/`.`Buchberger2.m`); #readlib('Ipairsort',``.`indiff`.`/`.`code`.`/`.`Ipairsort.m`); #readlib('makeDTlist',``.`indiff`.`/`.`code`.`/`.`makeDTlist.m`); #with(Groebner,normalf); #some simple syntax checking follows if type(G1,set)=false and type(G1,list)=false then ERROR(`first argument must be a list or set`) fi; if type(G2,set)=false and type(G2,list)=false then ERROR(`first argument must be a list or set`) fi; #readlib('termorders',``.`indiff`.`/`.`code`.`/`.`termorders.m`); if not member(myorder,termorders()) then ERROR(`myorder must be one of`,termorders()) fi; if nops(G1)+nops(G2)<2 then probset:={};XXX:={};F:={op(G1),op(G2)}; RETURN({op(G1),op(G2)}) fi; #if the "to do" pile is empty we exit if nops(G2)=0 then probset:={};XXX:={};F:={op(G1)}; RETURN({op(G1)}) fi; #we get the denominators of the Kmatrix; these are NON-ZERO Ktmp:=convert(Kmatrix,listlist): Ktmp:={seq(op(op(j,Ktmp)),j=1..nops(Ktmp))}: Kdenoms:=select(has,map(denom,Ktmp),`In`); Kdenoms:=map(m-> if type(factor(m),`*`) then op(factor(m)) elif type(factor(m),`^`) then op(1,factor(m)) else m fi,Kdenoms); Kdenoms:=select(has,Kdenoms,`In`); #the next few lines of code are to determine if extra options #and additional reduction are required NZ:={}; otherF:={};allvars2:=allvars; lshow:=false;lfactor:=false;#llinear:=false; if nargs>9 then for tt in [args[10..nargs]] do if not type(tt,equation) then ERROR(`invalid option argument`,tt) fi; s:=op(1,tt);t:=op(2,tt); if s='strategy' then if member('factor',t) then lfactor:=true fi; morestrat:={op(t)} minus {'factor'}; if nops(morestrat)<>0 then lprint(`these strategies are not implemented:`,morestrat) fi; fi; if s='info' then if member('shownew',t) then lshow:=true fi; fi; if s='NonZero' then NZ:=t fi; if s='extraReduction' then otherF:=t fi; if s='extraVarlist' then allvars2:=t fi; od; fi; #we use allukns to test for inconsistencies #allukns:=[seq(op(allvars[k]0,k=2..nops(allvars))]; probset2:={};A2:={}; GG:=[op(G1),op(G2)]; XX:={}; top:=nops(G1);top2:=nops(G2)+top; if top>1 then pairsdone:={seq(seq({k,j},j=k+1..top),k=1..top-1)}; pairstodo:=[seq(seq({k,j},j=top+1..top2),k=1..top), seq(seq({k,j},j=k+1..top2),k=top+1..top2-1)]; else pairsdone:={};pairstodo:=[seq(seq({k,j},j=k+1..top2),k=1..top2-1)]; fi; Ipairsort(pairstodo,GG,allvars,myorder,'newP'); while newP<>[] do pairnow:=op(1,newP); newP:=subsop(1=NULL,newP); test:=Buchberger2(pairnow,pairsdone,GG,allvars,myorder); pairsdone:=pairsdone union {pairnow}; if test=true then next fi; i:=op(1,pairnow);j:=op(2,pairnow); IdSpoly(op(i,GG),op(j,GG),allvars,myorder,'h'); if has(%,undefined) then next fi; Ireduce(h,A2 union {op(GG)},allvars,myorder,'h','Xset'='X', 'algebra'=[A,matorder]); XX:=XX union X; h:=op(numer(noKfactors([h],Kdenoms union NZ))); #the next bit is the additional reduction if otherF<>{} if otherF<>{} then for m to nops(otherF) do Ireduce(h,[op(m,otherF)],allvars2,myorder,'h') od; fi; if h<>0 then #if not has(h,allukns) then ERROR(h,`system inconsistent`) fi; #now follow the commands if the show and factor options are true if lshow=true and nops(h)<20 then lprint(h);lprint(``) elif lshow=true and nops(h)>=20 then Idiffparse2(h,allvars,myorder,'hdt','hp','hc'); lprint(`equation has`,nops(h),`summands, its HDT is`,hdt); lprint(`its Hpower is`,hp); lprint(`its Hcoeff is`,factor(hc));lprint(``) fi; if lfactor=true then m:=factor(h); if type(m,`*`) then mlist:=map(proc(v) if type(v,numeric) or member(v,NZ) or member(expand(-v),NZ) then NULL else v fi; end, [op(m)]); if mlist=[] then ERROR(m,`inconsistency found`) fi; if nops(mlist)=1 then h:=op(1,mlist) else lprint(`here are the factors`); for k to nops(mlist) do lprint(`factor number`,k); if type(op(k,mlist),`^`) then melt:=op(1,op(k,mlist)) else melt:=op(k,mlist) fi; if nops(melt)<21 then lprint(melt);lprint(``); else Idiffparse2(melt,allvars,myorder, 'hdtk','hpk','hck'); lprint(`the`,k,`th factor has`, nops(melt),`summands`); lprint(`its HDT is`,hdtk); lprint(`its Hpower is`,hpk); if nops(hck)<21 then lprint(`its Hcoeff is`,factor(hck)) fi; lprint(``); fi; od; k:=readstat(`Input factor number required >`); h:=op(k,mlist); fi; Ireduce(h,GG,allvars,myorder,'h','Xset'='X', 'algebra'=[A,matorder]);XX:=XX union X; fi; fi; #we return to the main part of the algorithm # if degree is too low we put into probset2 and leave for the algebraic part degofdiff(h,allvars,myorder,'fdeg'); if op(2,fdeg)<=DegDiffTable[op(1,fdeg)]+1 then probset2:=probset2 union {h} fi; if op(2,fdeg)op(0,HDT[f2]) or op(1,HDT[f1])<>op(1,HDT[f2])) then Spoly:=0; else Idiffcofs(HDT[f1],HDT[f2],'dcof1','dcof2'); if nops(HDT[f1])=3 then flag:=HDT[f1] else flag:='cum' fi; Idiffindex(f1temp,dcof1,'df1',flag); if has(%,undefined) then RETURN(undefined) fi; if nops(HDT[f2])=3 then flag:=HDT[f2] else flag:='cum' fi; Idiffindex(f2temp,dcof2,'df2',flag); if has(%,undefined) then RETURN(undefined) fi; if dcof1=[] then Idiffparse2(df1,allvars,termorder,'HDT[df1]','Hpower[df1]','Hcoeff[df1]'); elif nops(HDT[f1])=2 then HDT[df1]:=In[op(1,HDT[f1]),sort([op(op(2,HDT[f1])),op(dcof1)])]; Hpower[df1]:=1; coefficient(expand(df1),HDT[df1],'Hcoeff[df1]'); elif nops(HDT[f1])=3 then HDT[df1]:=In[op(1,HDT[f1]),op(2,HDT[f1]),[seq(op(j,op(3,HDT[f1]))+ op(j,dcof1))]]; Hpower[df1]:=1; coefficient(expand(df1),HDT[df1],'Hcoeff[df1]'); fi; if dcof2=[] then Idiffparse2(df2,allvars,termorder,'HDT[df2]','Hpower[df2]','Hcoeff[df2]'); elif nops(HDT[f2])=2 then HDT[df2]:=In[op(1,HDT[f2]),sort([op(op(2,HDT[f2])),op(dcof2)])]; Hpower[df2]:=1; coefficient(expand(df2),HDT[df2],'Hcoeff[df2]'); elif nops(HDT[f2])=3 then HDT[df2]:=In[op(1,HDT[f2]),op(2,HDT[f2]),[seq(op(j,op(3,HDT[f2]))+ op(j,dcof2))]]; Hpower[df2]:=1; coefficient(expand(df2),HDT[df2],'Hcoeff[df2]'); fi; algcofs(HDT[df1],Hcoeff[df1],Hpower[df1],HDT[df2],Hcoeff[df2],\ Hpower[df2],'acof1','acof2'); Spoly:=expand(acof1*df2-acof2*df1); fi; f1Sf2:=simplify(numer(Spoly)); end: Idiff:=proc(G,n) local idx,sol,j,Args2,whichU; global Neqs,HNIlist,ukns,vars; #readlib('Error',``.`indiff`.`/`.`code`.`/`.`Error.m`): #readlib('IndiffSimp',``.`indiff`.`/`.`code`.`/`.`IndiffSimp.m`); if type(G,`+`) then RETURN(map(Idiff,G,n)) fi; if type(G,`*`) then sol:=0;for j to nops(G) do sol:=sol+subsop(j=Idiff(op(j,G),n),G) od; RETURN(sol) fi; if type(G,`^`) then RETURN(op(2,G)*op(1,G)^(op(2,G)-1)*Idiff(op(1,G),n)) fi; #only cater for numerical powers for now!! if not has(G,'In') then RETURN(0) fi: if not type(G,indexed) and not type(G,function) then RETURN(IndiffSimp(diff(G,In[vars[n]]))) fi; #we need a rule for F(In[u,[]]) etc if type(G,function) then sol:=0;for j to nops(G) do sol:=sol+diff(G,op(j,G))*Idiff(op(j,G),n); od; RETURN(convert(sol,D)) fi; #we need a rule for In[f,[u1,u2],index] #note that this index is additive multi-index if nops(G)=3 then if type(args[2],integer) then idx:=op(3,G); sol:=0; for j to nops(op(2,G)) do sol:=sol+In[op(1,G),op(2,G), subsop(j=1+op(j,op(3,G)),op(3,G))]*Idiff(In[op(j,op(2,G)),[]],n) od; RETURN(sol); elif type(args[2],list) then Args2:=args[2]; whichU:=Args2[1][Args2[2]]; member(whichU,op(2,G),j); sol:=In[op(1,G),op(2,G),subsop(j=1+op(j,op(3,G)),op(3,G))]; RETURN(sol); fi; fi; if not type(args[2],integer) then ERROR(args,`Idiff undefined`) fi; #we need a rule for In[x] if nops(G)=1 then diff(G,In[vars[n]])+Error(n,op(1,G)); RETURN(simplify(%,Neqs)) fi; #we're down to the usual In terms idx:=sort([op(op(nops(G),G)),n]); #simplify(subsop(nops(G)=idx,G),select(has,Neqs,ukns),HNIlist) #+Error(sort(op(nops(G),G)),n,op(1,G)); subsop(nops(G)=idx,G)+Error(n,op(1,G),sort(op(nops(G),G))); IndiffSimp(%); end: Idiffcofs:=proc(HDT1,HDT2,dcof1,dcof2) local j,pp1,pp2,Nvars,ind1,ind2,dcof1temp,dcof2temp; option remember; #readlib('countj',``.`indiff`.`/`.`code`.`/`.`countj.m`): #procedure assumes HDT1, HDT2 are in In[] form if op(1,HDT1)<>op(1,HDT2) then ERROR(HDT1,HDT2,`diffcofs undefined`) fi; ind1:=op(nops(HDT1),HDT1): ind2:=op(nops(HDT2),HDT2): if nops(HDT1)=2 then dcof1temp:=[]: dcof2temp:=[]: Nvars:=convert(ind1,set) union convert(ind2,set); for j in Nvars do pp1:=max(countj(ind1,j),countj(ind2,j))-countj(ind1,j); pp2:=max(countj(ind1,j),countj(ind2,j))-countj(ind2,j); dcof1temp:=[op(dcof1temp),j$pp1]; dcof2temp:=[op(dcof2temp),j$pp2]; od; dcof1:=dcof1temp; dcof2:=dcof2temp; elif nops(HDT1)=3 then if not (nops(ind1)=nops(ind2) and op(2,HDT1)=op(2,HDT2)) then ERROR(HDT1,HDT2,`diffcofs undefined`) fi; #the result is being used for Idiffindex and thus needs to be #the nonadditive type dcof1:=[seq(max(op(j,ind1),op(j,ind2))-op(j,ind1),j=1..nops(ind1))]; dcof2:=[seq(max(op(j,ind1),op(j,ind2))-op(j,ind2),j=1..nops(ind1))]; fi; end: Idiffdivide:=proc(HDT1,HDT2,alpha) local Z,m,pp,ind1,ind2,j,alphatemp,vars1,vars2; option remember; # readlib('countj',``.`indiff`.`/`.`code`.`/`.`countj.m`): #calculates alpha such that D^{alpha}HDT1=HDT2 if op(1,HDT1)<>op(1,HDT2) then RETURN(false) fi; ind1:=op(nops(HDT1),HDT1): ind2:=op(nops(HDT2),HDT2): if nops(HDT1)=2 then vars1:=convert(ind1,set); vars2:=convert(ind2,set); if not(vars1 minus vars2 ={}) then RETURN(false) fi; alphatemp:=[]; for j to nops(vars2) do pp:=countj(ind2,op(j,vars2)) -countj(ind1,op(j,vars2)); if pp<0 then RETURN(false) fi; alphatemp:=[op(alphatemp),op(j,vars2)$pp]; od; alpha:=alphatemp; elif nops(HDT1)=3 then alphatemp:=[]; for j to nops(ind1) do #the lists of dependencies might be in a different order Z:=member(op(j,op(2,HDT1)),op(2,HDT2),'m'); if Z=false then ERROR(HDT1,HDT2,`inconsistent syntax`) fi; pp:=op(m,ind2)-op(j,ind1); if pp<0 then RETURN(false) fi; alphatemp:=[op(alphatemp),pp]; od; alpha:=alphatemp; fi; true end: Idiffindex:=proc(g,alpha,newg) local j,h,n; #this procedure differentiates g wrt the multi-index alpha, where #the alpha is assumed to be a multi-index #we need to do things in reverse order as not commutative #and that is how we interpret curlyD[idx] to operate #readlib('`Idiff`',``.`indiff`.`/`.`code`.`/`.`Idiff.m`): ###### if alpha=[] then newg:=g; RETURN(g) fi; h:=expand(g); if nargs=3 or ( nargs>3 and args[4]='cum') then for n from nops(alpha) by -1 to 1 do h:=Idiff(h,alpha[n]) od: else for n to nops(alpha) do if alpha[n]=0 then next fi; for j to alpha[n] do h:=Idiff(h,[op(2,args[4]),n]) od; od; fi; newg:=h; end: Idifforder:=proc (dt1, hp1, dt2, hp2, allvars, termorder) local deg1,deg2,varind,q,A,j,ind1,ind2,p1,p2,n,B,varlist,r,a1,a2,\ termvars1,termvars2,A1,A2,M,N,tt; global vars; option remember; #readlib('countj',``.`indiff`.`/`.`code`.`/`.`countj.m`): #yields 1 if dt1^hp1>dt2^hp2, 0 if equal, -1 if reverse inequality #we allow dti to be derviative terms OR monomials in the variables #note a variable may be in the form h(a) i.e. have type function if dt1 = dt2 then r:=traperror(max(hp1,hp2)); if type(r,function) then RETURN(0) elif hp1 < hp2 then RETURN(-1) elif hp2 < hp1 then RETURN(1) else RETURN(0) fi; fi; #next eliminate all the silly stuff if type(dt1,{indexed,function}) and op(0,dt1)='In' and not (type(dt2,{indexed,function}) and op(0,dt2)='In') then RETURN(1) fi; if type(dt2,{indexed,function}) and op(0,dt2)='In' and not (type(dt1,{indexed,function}) and op(0,dt1)='In') then RETURN(-1) fi; #then do case where both are derivative terms if type(dt1,indexed) and op(0,dt1)='In' and nops(dt1)>1 and type(dt2,indexed) and op(0,dt2)='In' and nops(dt2)>1 then A:=[seq(op(op(j,allvars)),j=2..nops(allvars))]; if not member(termorder,{'mtdeg','ttdeg','mrlex'}) then member(op(1,dt1),A,'p1'); member(op(1,dt2),A,'p2'); if p1 < p2 then RETURN(-1) elif p2 < p1 then RETURN(1) fi fi; if member(termorder,{'mtdeg','mrlex'}) then for j from 2 to nops(allvars) do if member(op(1,dt1),op(j,allvars)) then a1:=j fi; if member(op(1,dt2),op(j,allvars)) then a2:=j fi; od: if a1a2 then RETURN(1) fi; fi; ind1 := op(nops(dt1),dt1); ind2 := op(nops(dt2),dt2); if nops(dt1)=2 then deg1:=nops(ind1) else deg1:=convert(ind1,`+`) fi; if nops(dt2)=2 then deg2:=nops(ind2) else deg2:=convert(ind2,`+`) fi; if member(termorder,{'mtdeg','ttdeg','tdeg','rlex','mrlex'}) and deg1 deg2 then RETURN(1) fi; if member(termorder,{'mtdeg','ttdeg','mrlex'}) then member(op(1,dt1),A,'p1'); member(op(1,dt2),A,'p2'); if p1 < p2 then RETURN(-1) elif p2 < p1 then RETURN(1) fi; fi; #we should now have that op(1,hdt1)=op(1,hdt2) if nops(dt1)=2 and nops(dt2)=2 then varlist:=allvars[1]; B:=nops(varlist); if termorder ='lex' then for j from B by -1 to 1 do if countj(ind1,varlist[j])1 then RETURN(-1) fi; if nops(dt1)>1 and nops(dt2)=1 then RETURN(1) fi; #we are down to monomials in the indepts varind:=[seq(`'In'`[vars[op(j,op(1,allvars))]],j=1..nops(vars))]; #no difference here between tdeg and ttdeg if member(termorder,{'mtdeg','tdeg','ttdeg','rlex'}) then r:=traperror(max(degree(dt1,varind),degree(dt2,varind))); if type(r,function) then RETURN(0) elif degree(dt1,varind)>degree(dt2,varind) then RETURN(1) elif degree(dt1,varind)degree(dt2,{op(n,varind)}) then RETURN(1) elif degree(dt1,{op(n,varind)})degree(dt2,{op(n,varind)}) then RETURN(1) elif degree(dt1,{op(n,varind)})degree(dt2,{op(n,varind)}) then RETURN(-1) fi; od; fi; ############ fi; end: Idifforder2:=proc (dt1,hp1,hc1,dt2,hp2,hc2,allvars,termorder) local test,X,Y,A,B,L1,L2,f,g,Hh,Hg; option remember; #yields 1 if dt1^hp1>dt2^hp2, -1 if reverse inequality #if the hdt's and hp's are equal, examine the hc's #do recursively on the hc's #readlib('Idifforder',``.`indiff`.`/`.`code`.`/`.`Idifforder.m`): #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`): #readlib('IHmon',``.`indiff`.`/`.`code`.`/`.`IHmon.m`): test:=Idifforder(dt1,hp1, dt2, hp2,allvars,termorder); L1:=hc1; L2:=hc2; while test=0 do if L1=L2 then RETURN(0) fi; if has(L1,'In') and not has(L2,'In') then RETURN(1) fi; if not has(L1,'In') and has(L2,'In') then RETURN(-1) fi; if not has(L1,'In') and not has(L2,'In') then if type(L1,numeric) and type(L2,numeric) then if L1icontent(Hg) then RETURN(1) else RETURN(0) fi; fi; test:=Idifforder(f,1,g,1,allvars,termorder); break; fi; Idiffparse2(L1,allvars,termorder,'X','Y','L1'); Idiffparse2(L2,allvars,termorder,'A','B','L2'); test:=Idifforder(X,Y, A,B,allvars,termorder) od; test; end: Idiffparse:=proc(f,allvars,termorder,HDT,Hpower,Hcoeff) local j,DT,Sep, SepDT,L,L1,L2,L3,L4,L5,power,C,HDTtemp,g,n,LS,LS1,nhdt,newg, incon,v3,v4,vend1,vend2,vend3; option remember; #readlib('Idifforder',``.`indiff`.`/`.`code`.`/`.`Idifforder.m`); #readlib('termorders',``.`indiff`.`/`.`code`.`/`.`termorders.m`); if not member(termorder,termorders()) then ERROR(`termorder must be one of`,termorders()) fi; g:=f; g:=expand(g); if g=0 then Hcoeff:=1;DT:={};Hpower:=1;HDT:=0;Sep:=1;SepDT:={};RETURN(0) fi; L:={select(has,1+g,'In')} minus {0}: if (L={}) then incon:=op(g);#save incon, `inconsis.m`; ERROR(f,`equation contains no invariantised derivatives. Check input otherwise inconsistency generated. Condition saved in the file, inconsis.m, under the name incon`) fi; if type(op(L),`+`) then L:=convert(op(L),set) fi; L2:={}; for n in L do if type(n,`*`) then L2:=L2 union {op(n)} else L2:=L2 union {n} fi od; L3:=select(has,L2,'In'); for n in L3 do if type(n,`^`) then L4[n][1]:=op(1,n); L4[n][2]:=op(2,n) else L4[n][1]:=n; L4[n][2]:=1; fi od; L5:={}; for n in L3 do L5:=L5 union{L4[n][1]} od; L5:=convert(L5,list); #at this point we have a list of all derivative terms, without powers, in f HDTtemp:=op(1,L5); for n from 2 to nops(L5) do if Idifforder(HDTtemp,1,op(n,L5),1,allvars,termorder)=-1 then HDTtemp:=op(n,L5) fi; od; power:=1; for n in L3 do if (L4[n][1]=HDTtemp and L4[n][2]>power) then power:=op(2,n) fi od; #next we get the coefficient L1:={}; for n in L do if has(n,HDTtemp^power) then L1:=L1 union {n} fi od; C:=convert(L1,`+`): C:=expand(C/HDTtemp^power); if nargs>6 then assign(args[7],L3) fi; if nargs>7 then newg:=frontend(collect,[f,HDTtemp,distributed]); LS1:=convert({seq(j*HDTtemp^(j-1)*frontend(coeff,[newg,HDTtemp,j]),j=1..power)},`+`); assign(args[8],LS1); fi; if nargs >8 then if type(LS1,`+`) then LS:=convert(LS1,set) else LS:={LS1} fi; LS1:={}; for n in LS do if type(n,`*`) then LS1:=LS1 union {op(n)} else LS1:=LS1 union {n} fi; od; assign(args[9],select(has,LS1,'In')); fi; Hpower:=power; Hcoeff:=C; HDT:=HDTtemp; v3:=factor(C); vend1:=``.`HDT: `.HDTtemp.`, `.`Hpower: `.power; vend2:=``.`Hcoeff: `.v3; if nargs>7 then v4:=eval(args[8]); vend3:=``.`Sep: `.v4; linalg[matrix](3, 1, [vend1, vend2, vend3]); else linalg[matrix](2, 1, [vend1, vend2]); fi; end: Idiffparse2:=proc(f,allvars,termorder,HDT,Hpower,Hcoeff) local r,j,DT,Sep,SepDT,L,L1,L2,L3,L4,L5,power,cC,HDTtemp,g,n,LS,LS1,newg, incon; global vars,ukns; option remember; #readlib('Idifforder',``.`indiff`.`/`.`code`.`/`.`Idifforder.m`); #readlib('termorders',``.`indiff`.`/`.`code`.`/`.`termorders.m`); if not member(termorder,termorders()) then ERROR(`termorder must be one of`,termorders()) fi; g:=f; g:=expand(g); if g=0 then Hcoeff:=1;DT:={};Hpower:=1;HDT:=0;Sep:=1;SepDT:={};RETURN(0) fi; #the next bit is to solve a bit of a problem, as #select(has,In[f,[u],[1]]^2,'In') gives a strange error message r:=traperror({select(has,g,'In')}); if r=`non algebraic terms in power should be of the same type` then L:={select(has,1+g,'In')} minus {0}; else L:=r fi; if (L={}) then incon:=op(g);#save incon, `inconsis.m`; ERROR(f,`equation contains no invariantised derivatives. Check input otherwise inconsistency generated. Condition saved in the file, inconsis.m, under the name incon`) fi; if type(op(L),`+`) then L:=convert(op(L),set) fi; L2:={}; for n in L do if type(n,`*`) then L2:=L2 union {op(n)} else L2:=L2 union {n} fi od; L3:=select(has,L2,'In'); for n in L3 do if type(n,`^`) then L4[n][1]:=op(1,n); L4[n][2]:=op(2,n) else L4[n][1]:=n; L4[n][2]:=1; fi od; L5:={}; for n in L3 do L5:=L5 union{L4[n][1]} od; L5:=convert(L5,list); #at this point we have a list of all derivative terms, without powers, in f HDTtemp:=op(1,L5); for n from 2 to nops(L5) do if Idifforder(HDTtemp,1,op(n,L5),1,allvars,termorder)=-1 then HDTtemp:=op(n,L5) fi; od; power:=1; for n in L3 do if (L4[n][1]=HDTtemp and L4[n][2]>power) then power:=op(2,n) fi od; #next we get the coefficient L1:={}; for n in L do if has(n,HDTtemp^power) then L1:=L1 union {n} fi od; cC:=convert(L1,`+`): cC:=expand(cC/HDTtemp^power); if nargs>6 then assign(args[7],L3) fi; if nargs>7 then newg:=frontend(collect,[f,HDTtemp,distributed]); LS1:=convert({seq(j*HDTtemp^(j-1)*frontend(coeff,[newg,HDTtemp,j]),j=1..power)},`+`); assign(args[8],LS1); fi; if nargs >8 then if type(LS1,`+`) then LS:=convert(LS1,set) else LS:={LS1} fi; LS1:={}; for n in LS do if type(n,`*`) then LS1:=LS1 union {op(n)} else LS1:=LS1 union {n} fi; od; assign(args[9],select(has,LS1,'In')); fi; Hpower:=power; Hcoeff:=cC; HDT:=HDTtemp; end: Idtsort:=proc(F,allvars,termorder,newF) local Ftmp,X,i,v,j,sentinel,HDT,Hp,newFtemp2;#newFtemp if not type(F,list) and not type(F,set) then ERROR(F,`first arg must be a list or set`) fi; #readlib('Idifforder',``.`indiff`.`/`.`code`.`/`.`Idifforder.m`); #F is a list or set of DT's that we sort into increasing order sentinel:=1; Ftmp:={op(F)};#this removes multiple copies X:=[sentinel,op(Ftmp)];X:=convert(X,array); HDT[X[1]]:=sentinel; Hp[X[1]]:=1; for i from 2 to nops(Ftmp)+1 do if type(X[i],`^`) then HDT[X[i]]:=op(1,X[i]);Hp[X[i]]:=op(2,X[i]) else HDT[X[i]]:=X[i];Hp[X[i]]:=1 fi od; for i from 3 to nops(Ftmp)+1 do v:=X[i];j:=i; while Idifforder(HDT[X[j-1]],Hp[X[j-1]],HDT[v],Hp[v],allvars,termorder)=1 do X[j]:=X[j-1];j:=j-1 od; X[j]:=v od; #for j from 2 to nops(Ftmp)+1 do newFtemp[j-1]:=X[j] od; #newFtemp2:=[]; #for j to nops(F) do #newFtemp2:=[op(newFtemp2),newFtemp[j]] od; newFtemp2:=[seq(X[j],j=2..nops(Ftmp)+1)]: newF:=newFtemp2 end: Ieqnsort:=proc(F,allvars,termorder,newF) local X,i,v,j,nopF,sentinel,newFtemp,HDT,Hp,Hcoeff,n,newFtemp2; if not type(F,list) and not type(F,set) then ERROR(`first arg\ must be a list or set`) fi; #readlib('Idifforder2',``.`indiff`.`/`.`code`.`/`.`Idifforder2.m`); #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`); #F is a list or set of equations that we sort into increasing order #the sorting algorithm is insertion sort if nops(F)=1 then newF:=F;RETURN(F) fi; #the following is to get rid of multiple copies of elements X:=convert(F,set);nopF:=nops(X); sentinel:=1; for n in X do Idiffparse2(n,allvars,termorder,'HDT[n]','Hp[n]','Hcoeff[n]'); od; X:=[sentinel,op(X)];X:=convert(X,array); HDT[X[1]]:=sentinel; Hp[X[1]]:=1; for i from 3 to nopF+1 do v:=X[i];j:=i; while Idifforder2(HDT[X[j-1]],Hp[X[j-1]],Hcoeff[X[j-1]],HDT[v],Hp[v],Hcoeff[v],allvars,termorder)=1 do X[j]:=X[j-1];j:=j-1 od; X[j]:=v od; for j from 2 to nopF+1 do newFtemp[j-1]:=X[j] od; newFtemp2:=[]; for j to nopF do newFtemp2:=[op(newFtemp2),newFtemp[j]] od; newF:=newFtemp2 end: #file `indiff/Indiff #defines Indiff as a readlibbed package in Maple #DELETED Invariantize:=proc(G) local fn, Gd,Vars,didx,p,r; global Neqs,vars,ukns,HNIlist; if type(G,list) then RETURN(map(Invariantize,G)) fi; if type(G,set) then RETURN(map(Invariantize,G)) fi; if type(G,`+`) then RETURN(map(Invariantize,G)) fi; if type(G,`*`) then RETURN(map(Invariantize,G)) fi; if type(G,`^`) then RETURN(map(Invariantize,G)) fi; if member(G,vars) then RETURN(simplify(In[G],Neqs)) fi; #readlib('IndiffSimp',``.`indiff`.`/`.`code`.`/`.`IndiffSimp.m`); if member(G,ukns) then RETURN(IndiffSimp(In[G,[]])) fi; Gd:=convert(G,D);if not type(Gd,function) then RETURN(G) fi; if not has(op(0,Gd),D) then if not member(op(0,G),ukns) then RETURN(map(Invariantize,G)) else RETURN( IndiffSimp( In[op(0,G),[]] ) ) fi; fi; fn:=op(op(0,Gd)); Vars:=[op(Gd)]: if nops(Vars)>1 then didx:=[op(op(0,op(0,Gd)))] elif op(0,op(0,Gd))=D then didx:=[1] else r:=member(op(Vars),vars,'p'); if r=false then ERROR(op(Vars),`undeclared variable`) fi; didx:=[p$op(2,op(0,op(0,Gd)))] fi: IndiffSimp(In[fn,didx]); end: Iorthreduce:=proc(f,F,allvars,termorder,newf) local HDT,Hp,Hc,n,l,i,G,g,DTl,pl,alpha,alphabar,DT,b,A,s,t,tt,lLOWSET,lalgred,nn,r,\ AA,matermorder,m,k,j,flag,h,g1,coeffa,coeffb,acofa,acofb,myZ,exponent,allukns,nDTlist; global vars,DegDiffTable; #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`): #readlib('coefficient',``.`indiff`.`/`.`code`.`/`.`coefficient.m`): #readlib('algcofs',``.`indiff`.`/`.`code`.`/`.`algcofs.m`): #readlib('Ieqnsort',``.`indiff`.`/`.`code`.`/`.`Ieqnsort.m`): #readlib('Idtsort',``.`indiff`.`/`.`code`.`/`.`Idtsort.m`): #readlib('Idifforder',``.`indiff`.`/`.`code`.`/`.`Idifforder.m`): #readlib('Idiffindex',``.`indiff`.`/`.`code`.`/`.`Idiffindex.m`): #readlib('Idiffdivide',``.`indiff`.`/`.`code`.`/`.`Idiffdivide.m`): #readlib('makeDTlist',``.`indiff`.`/`.`code`.`/`.`makeDTlist.m`): #with(Groebner,normalf); #sort out options lalgred:=false;#lLOW:=false; for tt in [args[6..nargs]] do if not type(tt,equation) then ERROR(`invalid option argument`,tt) fi; s:=op(1,tt);t:=op(2,tt); if s='algebra' then lalgred:=true; AA:=op(1,t);matermorder:=op(2,t) fi; #if s='info' then if member('LOWSET',t) then lLOW:=true fi fi; od; lLOWSET:={}; g:=expand(f); if (g=0) then newf:=0;RETURN(0) fi; allukns:=[seq(op(op(j,allvars)),j=2..nops(allvars))]; for n in F do Idiffparse2(n,allvars,termorder,'HDT[n]','Hp[n]','Hc[n]'); if nops(HDT[n])=2 then m:=nops(op(2,HDT[n])) elif nops(HDT[n])=3 then m:=convert(op(3,HDT[n]),`+`) #else -1 fi; if m0 then newf:=g/icontent(g); RETURN(g/icontent(g)) else newf:=0; RETURN(0) fi fi; if Idifforder(HDT[G[i]],Hp[G[i]],HDT[g],Hp[g],allvars,termorder)=1 then next fi; Idtsort(DT[g],allvars,termorder,'DT[g]'); for l from nops(DT[g]) by -1 to 1 do if type(op(l,DT[g]),`^`) then DTl:=op(1,op(l,DT[g])); pl:=op(2,op(l,DT[g])) else DTl:=op(l,DT[g]); pl:=1 fi; if Idifforder(HDT[G[i]],Hp[G[i]],DTl,pl,allvars,termorder)=1\ then break fi; #if (op(0,HDT[G[i]])<>op(0,DTl)) then next fi; if (op(1,HDT[G[i]])<>op(1,DTl)) then next fi; myZ:=Idiffdivide(HDT[G[i]],DTl,'A'); if not type(alpha[HDT[G[i]],DTl],list) then alpha[HDT[G[i]],DTl]:=A fi; if (myZ=false) then next else if nops(HDT[G[i]])=3 then flag:=HDT[G[i]] else flag:='cum' fi; Idiffindex(G[i],alpha[HDT[G[i]],DTl],'h',flag); if has(%,undefined) then next fi; alphabar:=convert(alpha[HDT[G[i]],DTl],`+`); if (alphabar=0 and pl-Hp[G[i]]<0) then next fi; if alphabar=0 then coefficient(h,DTl^Hp[G[i]],'coeffa') else coefficient(h,DTl,'coeffa') fi; coefficient(g,DTl^pl,'coeffb'); algcofs(1,coeffa,1,1,coeffb,1,'acofa','acofb'); if (alphabar=0) then exponent:=pl-Hp[G[i]] else exponent:=pl-1 fi; if has(acofa,'In') then next fi; g1:=numer(map(simplify,expand(acofa*g-acofb*(DTl^exponent)*h))); if (g1=0) then newf:=0;RETURN(0) fi; fi; l:=0; if g1<>g then b:=nops(G)+1;g:=g1 fi; break; od; od; for j from nops(lLOWSET) by -1 to 1 do n:=op(j,lLOWSET); nn:=nops(nDTlist[n]); g:=numer(simplify(g,{op(j,lLOWSET)},[seq(nDTlist[n][nn-k],k=0..nn-1)])); od; if lalgred=true then r:=traperror(normalf(g,AA,matermorder)); if r=`polynomials not members of the algebra` or nops(AA)<2 then g:=simplify(g,AA) else g:=r fi; fi; if g=0 then newf:=0 else newf:=g/icontent(g) fi; end: Iorthreduceall:=proc(F,allvars,termorder,newF) local newn,Ftemp,z,j,n; #reduces F to an auto-reduced list, assumes that deg(F)>N=1 if F={} or F=[] or nops(F)=1 then newF:=F; RETURN(newF) fi; #readlib('Iorthreduce',``.`indiff`.`/`.`code`.`/`.`Iorthreduce.m`); #readlib('Ieqnsort',``.`indiff`.`/`.`code`.`/`.`Ieqnsort.m`): if not type(F,set) and not type(F,list) then \ ERROR(`first argument must be a list or set`) fi; Ftemp:=F; z:=0; while z=0 do Ieqnsort(Ftemp,allvars,termorder,'Ftemp'); if nops(Ftemp)=1 then break fi; for j from 2 to nops(Ftemp) do n:=op(j,Ftemp); Iorthreduce(n,{op(1..j-1,Ftemp)},allvars,termorder,'newn',args[5..nargs]); if newn=n then z:=1;next fi; Ftemp:=convert(Ftemp,set) minus{n}; if newn<>0 then Ftemp:=Ftemp union {newn} fi; z:=0; j:=1; break; od od; newF:=Ftemp; end: Ipairsort:=proc(P,GG,allvars,termorder,newP) local X,i,v,j,sentinel,HDT,termorder2,falsetmp,top; #readlib('Idifforder',``.`indiff`.`/`.`code`.`/`.`Idifforder.m`); #readlib('LCD',``.`indiff`.`/`.`code`.`/`.`LCD.m`); #F is a list or set of pairs of integers that we sort into increasing order #according to the LCD of the corresponding elements of GG sentinel:=1; X:=map(proc(m) if LCD(op(op(1,m),GG),op(op(2,m),GG),allvars,termorder)=false then NULL else m fi; end,P); if member(termorder,{'lex','alex'}) then termorder2:='tdeg' else termorder2:=termorder fi; for j in X do HDT[j]:=LCD(op(op(1,j),GG),op(op(2,j),GG),allvars,termorder) od; HDT[1]:=1; top:=nops(X); X:=convert([sentinel,op(X)],array); for i from 3 to top+1 do v:=X[i];j:=i; while Idifforder(HDT[X[j-1]],1,HDT[v],1,allvars,termorder2)=1 do X[j]:=X[j-1];j:=j-1 od; X[j]:=v od; newP:=subsop(1=NULL,convert(X,list)); end: Ireduce:=proc(f,F,allvars,termorder,newf) local HDT,Hp,Hc,n,l,i,G,g,DTl,pl,alpha,alphabar,DT,b,XX,A,AA,matermorder,lXset,\ flag,h,g1,coeffa,coeffb,acofa,acofb,myZ,exponent,lalgred,Xsetname,\ lLOWSET,nDTlist,allukns,k,j, m, nn, r, s, t, tt; global DegDiffTable,vars; #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`): #readlib('coefficient',``.`indiff`.`/`.`code`.`/`.`coefficient.m`): #readlib('algcofs',``.`indiff`.`/`.`code`.`/`.`algcofs.m`): #readlib('Ieqnsort',``.`indiff`.`/`.`code`.`/`.`Ieqnsort.m`): #readlib('Idtsort',``.`indiff`.`/`.`code`.`/`.`Idtsort.m`): #readlib('Idifforder',``.`indiff`.`/`.`code`.`/`.`Idifforder.m`): #readlib('Idiffindex',``.`indiff`.`/`.`code`.`/`.`Idiffindex.m`): #readlib('Idiffdivide',``.`indiff`.`/`.`code`.`/`.`Idiffdivide.m`): #readlib('degofdiff',``.`indiff`.`/`.`code`.`/`.`degofdiff.m`): #with(Groebner,normalf); lalgred:=false;lXset:=false; for tt in [args[6..nargs]] do if not type(tt,equation) then ERROR(`invalid option argument`,tt) fi; s:=op(1,tt);t:=op(2,tt); if s='algebra' then lalgred:=true; AA:=op(1,t);matermorder:=op(2,t) fi; if s='Xset' then lXset:=true; Xsetname:=t fi; od; XX:={}; g:=simplify(expand(f)); if (g=0) then if lXset=true then assign(Xsetname,{}) fi;newf:=0;RETURN(0) fi; lLOWSET:={}; allukns:=[seq(op(op(j,allvars)),j=2..nops(allvars))]; for n in F do Idiffparse2(n,allvars,termorder,'HDT[n]','Hp[n]','Hc[n]') ; if nops(HDT[n])=2 then m:=nops(op(2,HDT[n])) elif nops(HDT[n])=3 then m:=convert(op(3,HDT[n]),`+`) #else -1 fi; if m0 then newf:=g/icontent(g); RETURN(g/icontent(g)) else newf:=0; RETURN(0) fi fi; if Idifforder(HDT[G[i]],Hp[G[i]],HDT[g],Hp[g],allvars,termorder)=1 then next fi; Idtsort(DT[g],allvars,termorder,'DT[g]'); for l from nops(DT[g]) by -1 to 1 do if type(op(l,DT[g]),`^`) then DTl:=op(1,op(l,DT[g])); pl:=op(2,op(l,DT[g])) else DTl:=op(l,DT[g]); pl:=1 fi; if Idifforder(HDT[G[i]],Hp[G[i]],DTl,pl,allvars,termorder)=1\ then break fi; if (op(1,HDT[G[i]])<>op(1,DTl)) then next fi; myZ:=Idiffdivide(HDT[G[i]],DTl,'A'); if not type(alpha[HDT[G[i]],DTl],list) then alpha[HDT[G[i]],DTl]:=A fi; if (myZ=false) then next else if nops(HDT[G[i]])=3 then flag:=HDT[G[i]] else flag:='cum' fi; Idiffindex(G[i],alpha[HDT[G[i]],DTl],'h',flag); if has(%,undefined) then next fi; alphabar:=convert(alpha[HDT[G[i]],DTl],`+`); h:=simplify(h); if (alphabar=0 and pl-Hp[G[i]]<0) then next fi; if alphabar=0 then coefficient(h,DTl^Hp[G[i]],'coeffa') else coefficient(h,DTl,'coeffa') fi; coefficient(g,DTl^pl,'coeffb'); algcofs(1,coeffa,1,1,coeffb,1,'acofa','acofb'); if (alphabar=0) then exponent:=pl-Hp[G[i]] else exponent:=pl-1 fi; g1:=numer(map(simplify,expand(acofa*g-acofb*(DTl^exponent)*h))); XX:=XX union {acofa}; if (g1=0) then if lXset=true then assign(Xsetname,XX) fi; newf:=0;RETURN(0) fi; fi; l:=0; if g1<>g then b:=nops(G)+1;g:=g1 fi; break; od; od; for j from nops(lLOWSET) by -1 to 1 do n:=op(j,lLOWSET); nn:=nops(nDTlist[n]); g:=numer(simplify(g,{op(j,lLOWSET)},[seq(nDTlist[n][nn-k],k=0..nn-1)]));od; if lXset=true then assign(Xsetname,XX) fi; if lalgred=true then r:=traperror(normalf(g,AA,matermorder)); if r=`polynomials not members of the algebra` or nops(AA)<2 then g:=simplify(g,AA) else g:=r fi; fi; if g=0 then newf:=0 else newf:=g/icontent(g) fi; end: Ireduceall:=proc(F,allvars,termorder,newF) #local newn,Ftemp,z,j,n,XX,XXn,Param; local newn,Ftemp,z,j,n,XX,XXn; #Ireduces F to an auto-Ireduced list if F={} or F=[] or nops(F)=1 then newF:=F; RETURN(newF) fi; #readlib('Ireduce',``.`indiff`.`/`.`code`.`/`.`Ireduce.m`); #readlib('Ieqnsort',``.`indiff`.`/`.`code`.`/`.`Ieqnsort.m`): if not type(F,set) and not type(F,list) then \ ERROR(`first argument must be a list or set`) fi; Ftemp:=F; #if nargs>5 then Param:=args[6] else Param:={} fi; XX:={}; z:=0; while z=0 do Ieqnsort(Ftemp,allvars,termorder,'Ftemp'); if nops(Ftemp)=1 then break fi; for j from 2 to nops(Ftemp) do n:=op(j,Ftemp); Ireduce(n,{op(1..j-1,Ftemp)},allvars,termorder,'newn','XXn'); if newn=n then z:=1;next fi; XX:=XX union XXn; Ftemp:=convert(Ftemp,set) minus{n}; if newn<>0 then Ftemp:=Ftemp union {newn} fi; z:=0; j:=1; break; od od; if nargs>4 then assign(args[5],XX) fi; newF:=Ftemp; end: Kmat:=proc() local A,Jt,DTlist,J,T,PhiMatrix,j,n,nterm,idx,p; global Kmatrix,HNIlist,GroupP,vars,ukns,Neqs,XiPhis; #readlib('makeDTlist',``.`indiff`.`/`.`code`.`/`.`makeDTlist.m`); #readlib('PROL',``.`indiff`.`/`.`code`.`/`.`PROL.m`); #readlib('IndiffSimp',``.`indiff`.`/`.`code`.`/`.`IndiffSimp.m`); #with(linalg,jacobian);with(linalg,transpose);with(linalg,inverse); if nops(GroupP)<>nops(Neqs) then ERROR(`the number of normalisation equations must equal the number of group parameters`) fi; makeDTlist(Neqs,vars,ukns,'DTlist'); #we now have a list of derivative terms and indepts J:=jacobian(convert(Neqs,list),DTlist); Jt:=transpose(J); # next we make the PhiMatrix matrix #this assumes indept vars appear as In[x] #and the dept vars as In[u,[]] PhiMatrix:=matrix(nops(GroupP),nops(DTlist)); for j to nops(GroupP) do for n to nops(DTlist) do nterm:=op(n,DTlist); if member(op(1,nterm),vars,'p') then PhiMatrix[j,n]:=Invariantize(XiPhis[j,p]) else PhiMatrix[j,n]:=PROL(op(1,nterm),op(j,GroupP),op(nops(nterm),nterm)) fi; od; od; #next the T matrix T:=matrix(nops(vars),nops(DTlist)); for j to nops(vars) do for n to nops(DTlist) do nterm:=op(n,DTlist); if member(op(1,nterm),vars,'p') then if j=p then T[j,n]:=1 else T[j,n]:=0 fi; else idx:=sort([op(op(nops(nterm),nterm)),j]); T[j,n]:=IndiffSimp(subsop(nops(nterm)=idx,nterm)) fi; od; od; A:=inverse(evalm(PhiMatrix&*Jt)); Kmatrix:=evalm(T&*Jt&*A); #transpose(A); end: LCD:=proc(f1,f2,allvars,termorder) local hdt1,hdt2,hc1,hc2,hp1,hp2,outindex,lcdtmp,m,idx1,idx2,j; global vars; option remember; #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`): #readlib('countj',``.`indiff`.`/`.`code`.`/`.`countj.m`): Idiffparse2(f1,allvars,termorder,'hdt1','hp1','hc1'); Idiffparse2(f2,allvars,termorder,'hdt2','hp2','hc2'); if op(1,hdt1)<>op(1,hdt2) then RETURN(false) fi; if f1=1 then RETURN(false) fi; if hdt1=hdt2 then RETURN(hdt1) fi; lcdtmp:=false; if nops(hdt1)=3 then if op(2,hdt1)<>op(2,hdt2) then ERROR(`bad syntax`,hdt1,hdt2) fi; if nops(op(3,hdt1))<>nops(op(3,hdt2)) then ERROR(`bad syntax`,hdt1,hdt2) fi; idx1:=op(3,hdt1);idx2:=op(3,hdt2); outindex:=[seq(max(op(j,idx1),op(j,idx2)),j=1..nops(idx1))]; lcdtmp:='In'[op(1,hdt1),op(2,hdt1),outindex]; elif nops(hdt1)=2 then outindex:=NULL; for j to nops(vars) do m:=max(countj(op(2,hdt1),j),countj(op(2,hdt2),j)); if m=0 then next else outindex:=outindex,j$m fi; od; lcdtmp:='In'[op(1,hdt1),[outindex]]; fi; end: PROL:=proc(u,lambda,Ind) local xi,eta,diffseq,varslist,myf,fred; global vars,ukns,XiPhis,Neqs: option remember; #readlib('prolong',``.`indiff`.`/`.`code`.`/`.`prolong.m`): #readlib('proldata',``.`indiff`.`/`.`code`.`/`.`proldata.m`): proldata(lambda,xi,eta); fred:=op(eta[u]); varslist:=op(vars); if Ind=[] then myf:=u(varslist) else diffseq:=map(m->op(m,vars),Ind); myf:=diff(u(varslist),op(diffseq)) fi; Invariantize(eval(prolong(myf,vars,ukns,op(xi),op(eta)))) end: algcofs:=proc(HDT1,Hcoeff1,Hpower1,HDT2,Hcoeff2,Hpower2,\ algcof1,algcof2) local MM,NN,KK,Hlist1,Hlist2,Hset,Hset1,Hset2,n,temp1,temp2,j,j1,j2,HT1,HT2,m; Hlist1:=table(); Hlist2:=table(); #readlib(select); temp1:=factor(Hcoeff1); if type(temp1,`*`) then HT1:={op(temp1)} else HT1:={temp1} fi; temp2:=factor(Hcoeff2); if type(temp2,`*`) then HT2:={op(temp2)} else HT2:={temp2} fi; HT1:=HT1 union {HDT1^Hpower1} minus {1}; HT2:=HT2 union {HDT2^Hpower2} minus {1}; for n in HT1 do if type(n,`^`) then Hlist1[n]:=op(1,n) else Hlist1[n]:=n fi od; for n in HT2 do if type(n,`^`) then Hlist2[n]:=op(1,n) else Hlist2[n]:=n fi od; Hset1:=convert(Hlist1,set); Hset2:=convert(Hlist2,set); Hset:=Hset1 intersect Hset2; temp1:=1; temp2:=1; for n in HT1 do if not member(Hlist1[n],Hset) then temp1:=temp1*n else j1:=0; j2:=0; if type(n,`^`) then NN:=op(1,n);j1:=op(2,n) else NN:=n; j1:=1 fi; for m in HT2 do if Hlist2[m]=Hlist1[n] then if type(m,`^`) then MM:=op(1,m);j2:=op(2,m) else MM:=m; j2:=1; break fi; fi od; j:=max(j1,j2); KK:=simplify((j-j2)/j1); temp1:=temp1*NN^(j-j2); fi od; for n in HT2 do if not member(Hlist2[n],Hset) then temp2:=temp2*n else j1:=0; j2:=0; if type(n,`^`) then NN:=op(1,n);j2:=op(2,n) else NN:=n;j2:=1 fi; for m in HT1 do if Hlist1[m]=Hlist2[n] then if type(m,`^`) then j1:=op(2,m);MM:=op(1,m); else MM:=m;j1:=1;break fi; fi od; j:=max(j1,j2); KK:=simplify((j-j1)/j2); temp2:=temp2*NN^(j-j1); fi od; algcof1:=temp1; algcof2:=temp2; end: coefficient:=proc(f,term,coeff) local g,L2,n,z; g:=expand(f); if type(g,`+`) then g:=convert(g,set) else g:={g} fi; L2:=[]; for n in g do if has(n,term) then L2:=[op(L2),n] fi od; z:=convert(L2,`+`); z:=expand(z/term); coeff:=z end: countj:=proc(index,k) local j,n; j:=0; for n to nops(index) do if index[n]=k then j:=1+j fi; od; j; end: degofdiff:=proc(f,allvars,termorder,DL) #f is an idp, #N1= [highest unknown, degree of differentiation] local degf,HDT,Hp,Hc; #readlib('Idiffparse2',``.`indiff`.`/`.`code`.`/`.`Idiffparse2.m`); Idiffparse2(f,allvars,termorder,'HDT','Hp','Hc'); if nops(HDT)=2 then degf:=nops(op(2,HDT)) elif nops(HDT)=3 then degf:=convert(op(3,HDT),`+`) else degf:=-1 fi; DL:=[op(1,HDT),degf]; end: makeDTlist:=proc(F,vars,ukns,LL) local g,n,Ltmp,j,L1,L2,L3; Ltmp:={}; for j to nops(F) do g:=expand(op(j,F)); if g=0 then RETURN([]) fi; L1:={select(has,g,{op(vars),op(ukns)})}; if type(op(L1),`+`) then L1:=convert(op(L1),set) fi; L2:={}; for n in L1 do if type(n,`*`) then L2:=L2 union {op(n)} else L2:=L2 union {n} fi; od; L3:=select(has,L2,{op(vars),op(ukns)}); L3:=map(w->if type(w,`^`) then op(1,w) else w fi,L3); Ltmp:=Ltmp union L3; od; LL:=convert(Ltmp,list); end: makeorder:=proc(DTlist,allvars,T,C) local hdeg,ldeg,nuknlists,LL1,newlist,j,k,DTLtmp,countt,countt1,B,x,y,Mlist,K,p,Mlist1,S,V; #DTlist is a list of invariantized derivative terms #allvars is a list of lists of ukns # T is the output name if nops(DTlist)<2 then T:=`lexdeg`(DTlist); RETURN(`lexdeg`(DTlist)) fi; #readlib('diffcountn',``.`indiff`.`/`.`code`.`/`.`diffcountn.m`); LL1:=op(map(m->if nops(m)=2 then nops(op(2,m)) elif nops(m)=3 then convert(op(3,m),`+`) else NULL fi,DTlist)); if LL1=NULL then ERROR(`bad DTlist`,DTlist) fi; hdeg:=max(LL1);ldeg:=min(LL1); nuknlists:=nops(allvars)-1; if nuknlists<1 then ERROR(`bad allvars list`,allvars) fi; newlist:=NULL;countt:=NULL;DTLtmp:=DTlist; for j to nuknlists do for k from ldeg to hdeg do LL1:=select((f,L)->if member(op(1,f),L) then true else false fi,DTLtmp,allvars[j+1]); if select(diffcountn,LL1,k)<>[] then newlist:=op(select(diffcountn,LL1,k)),newlist; countt:=nops(select(diffcountn,LL1,k)),countt fi; DTLtmp:=remove((f,L)->if member(op(1,f),L) and diffcountn(f,k) then true else false fi,DTLtmp,allvars[j+1]); od; od; T:=[newlist]: K:=nops([countt]); #unfortunately, groebner was written by someone who #uses square matrices only!! if K=1 then C:=linalg[matrix](countt[1],countt[1], (x,y)->if x=1 or x-1=y then 1 else 0 fi); RETURN(op(C)) fi; Mlist:=linalg[blockmatrix](1,2,[linalg[matrix](countt[1],countt[1], (x,y)->if x=1 or x-1=y then 1 else 0 fi) , linalg[matrix](countt[1],sum('countt[p]','p'=2..K),(x,y)->0)]); for j from 2 to K-1 do S:=linalg[matrix](countt[j],sum('countt[p]','p'=1..j-1),(x,y)->0); V:=linalg[matrix](countt[j],sum('countt[p]','p'=j+1..K),(x,y)->0); B:=linalg[matrix](countt[j],countt[j],(x,y)->if x=1 or x-1=y then 1 else 0 fi); Mlist1:=linalg[blockmatrix](1,3,[S, B,V]); Mlist:=linalg[blockmatrix](2,1,[Mlist,Mlist1]); od; C:=linalg[blockmatrix](2,1,[Mlist, linalg[blockmatrix](1,2, [linalg[matrix](countt[K],sum('countt[p]','p'=1..K-1),(x,y)->0), linalg[matrix](countt[K],countt[K],(x,y)->if x=1 or x-1=y then 1 else 0 fi) ])]); op(C); end: diffcountn:=proc(f,n) if nops(f)=2 and nops(op(2,f))=n then true elif nops(f)=3 and convert(op(3,f),`+`)=n then true else false fi; end: prolong := proc(eqn,ivar,dvar,xi,eta) #local lfn,ii,lfind,eqnind,prol,expr,result,xi,eta: local lfn,ii,lfind,lfind2,eqnind,prol,result: #readlib('`indiff/desolv/gendl`',``.`indiff`.`/`.`code`.`/`.`gendl.m`); #readlib('`indiff/desolv/isderiv`',``.`indiff`.`/`.`code`.`/`.`isderiv.m`); #readlib('`indiff/desolv/transform1`',``.`indiff`.`/`.`code`.`/`.`transform1.m`); #readlib('`indiff/desolv/transform2`',``.`indiff`.`/`.`code`.`/`.`transform2.m`); #readlib('`indiff/desolv/orderd`',``.`indiff`.`/`.`code`.`/`.`orderd.m`); #readlib('`indiff/desolv/sortder`',``.`indiff`.`/`.`code`.`/`.`sortder.m`); #readlib('`indiff/desolv/map2ind`',``.`indiff`.`/`.`code`.`/`.`map2ind.m`); #readlib('`indiff/desolv/infin`',``.`indiff`.`/`.`code`.`/`.`infin.m`); #xi:=op(Fxi);eta:=op(Feta); lfn := select(`indiff/desolv/isderiv`,indets(eqn,function),dvar): lfn := `indiff/desolv/sortder`(convert(lfn,list),-1): lfind := map(`indiff/desolv/map2ind`,lfn): lfind2:=map(convert,map(convert,lfind,string),name): eqnind := `indiff/desolv/transform1`(eqn,lfn,lfind2,ivar,dvar): prol := seq([ZETA_(ivar[ii],ivar,dvar,xi,eta),ivar[ii]],ii=1..nops(ivar)): prol := prol,seq([ZETA_(dvar[ii],ivar,dvar,xi,eta),dvar[ii]],ii=1..nops(dvar)): #prol := [prol,seq([ZETA_(lfind[ii],ivar,dvar,xi,eta), #convert(lfind[ii],string)],ii=1..nops(lfind))]: prol := [prol,seq([ZETA_(lfind[ii],ivar,dvar,xi,eta),lfind2[ii]],ii=1..nops(lfind))]: result := map(proc(oper,ex) oper[1]*diff(ex,oper[2]) end,prol,eqnind): result := convert(result,`+`): result := expand(eval(subs(ZETA_ = `indiff/desolv/infin`,result))): result := `indiff/desolv/transform2`(result,lfn,lfind2,ivar,dvar): end: `indiff/desolv/isderiv` := proc(func,dvar) local der: if not type(func,function) then false elif has(op(0,func),D) then member(op(op(0,func)),dvar) elif has(op(0,func),diff) then der := func: while has(der,diff) do der := op(1,der) od: member(op(0,der),dvar) else false fi end: `indiff/desolv/sortder` := proc(lder) local num_der,lord,ii,num_: num_der := nops(lder): if num_der < 2 then RETURN(lder) fi: num_ := 10^length(num_der): lord := [seq(`indiff/desolv/orderd`(lder[ii])*num_+ii,ii=1..num_der)]: lord := sort(lord): lord := map(proc(a,b) irem(a,b) end,lord,num_): if nargs > 1 and args[2] = -1 then [seq(lder[lord[num_der-ii]],ii=0..num_der-1)] else [seq(lder[lord[ii]],ii=1..num_der)] fi end: `indiff/desolv/map2ind` := proc(der) local ldvar,func,lindex,lvar,ii: if op(0,der) = `diff` then ldvar := NULL: func := der: while op(0,func) = `diff` do ldvar := ldvar,op(2,func): func := op(1,func) od: func := op(0,func) else lindex := [op(op(0,op(0,der)))]: lvar := [op(der)]: ldvar := seq(lvar[lindex[ii]],ii=1..nops(lindex)): func := op(op(0,der)): fi: func[ldvar] end: `indiff/desolv/transform1` := proc(eqn,lfn,lstr,ivar,dvar) local result,lfdvar,lsubfn,lfd: #local result,ii,lder,der,tder,lfdvar,fun,lsubfn,lfd: # Replace functions/derivatives in 'lfn' by strings in 'lstr' respectively. # lsubfn := zip((x,y)-> x = convert(y,string),lfn,lstr): lsubfn := zip((x,y)-> x = y,lfn,lstr): result := subs(lsubfn,eqn): lfd := indets(result,function): lfdvar := select((f,ld)->member(op(0,f),ld),lfd,dvar): lfd := lfd minus lfdvar: lsubfn := map((f)-> f = op(0,f),lfdvar): result := subs(lsubfn,result): # convert 'D' into diff operator. lfd := subs(lsubfn,lfd): lfd := select((f)->has(op(0,f),`D`) and evalb(nops(f)=1) and type(op(f),name),lfd): lsubfn := map(proc(der) local ordder,func: if type(op(0,op(0,der)),function) then ordder := op(2,op(0,op(0,der))) else ordder := 1 fi: func := op(op(0,der))(op(der)): der = diff(func,op(func)$ordder) end,lfd): subs(lsubfn,result) end: `indiff/desolv/orderd` := proc(func) # deriv = a derivative. # Calculate the order of derivative of "deriv". #_______________________________________________________________________ local ord_der,der; if has(func,diff) then # func is a derivative in "diff" form. ord_der := 0: der := func; while has(der,diff) do ord_der := ord_der + 1: der := op(1,der) od: ord_der elif has(op(0,func),D) then # func is a derivative in 'D' form. if type(op(0,op(0,func)),indexed) then nops(op(0,op(0,func))) elif type(op(0,op(0,func)),function) then op(2,op(0,op(0,func))) else 1 fi else 0 fi # func is not a derivative. end: `indiff/desolv/transform2` := proc(eqn,lfn,lstr) local result,lder,der,mem,func1,func2,lsubfn,fvars,fvar,flag: # Replace strings in 'lstr' by functions/derivatives in 'lfn' respectively. # lsubfn := zip((x,y)->convert(x,string)=y,lstr,lfn): lsubfn := zip((x,y)->x=y,lstr,lfn): result := subs(lsubfn,eqn): lder := indets(result,function): lder := select((f)->has(op(0,f),`diff`),lder): lder := `indiff/desolv/sortder`(lder): # convert diff into 'D' operator. lsubfn := zip((x,y)->x=y,lstr,lfn): for der in lder do func1 := der: while has(op(0,func1),diff) do func1 := op(1,func1) od: fvars := [op(func1)]: func2 := der: flag := false: for fvar in fvars do if member(fvar,lfn,'mem') then func2 := subs(fvar=lstr[mem],func2): flag := true fi od: if flag then func2 := convert(func2,D): func2 := subs(lsubfn,func2): result := subs(der=func2,result) fi: od: result end: `indiff/desolv/infin` := proc(func,ivar,dvar,xi,eta) local ii,result,zeta_km1,fname,lvar,ldvar: lvar := op(`indiff/desolv/gendl`(ivar,dvar,0)): if not type(func,indexed) then if member(func,ivar) then xi[func](lvar) else eta[func](lvar) fi elif nops(func) = 1 then fname := op(0,func): ldvar := [op(func)]: result := [seq(diff(fname(op(ivar)),ivar[ii])* diff(xi[ivar[ii]](lvar),op(func)),ii=1..nops(ivar))]: diff(eta[fname](lvar),op(func)) - expand(convert(result,`+`)) else fname := op(0,func): ldvar := [op(func)]: zeta_km1 := `indiff/desolv/infin`(subsop(1=NULL,func),ivar,dvar,xi,eta): result := [seq(diff(fname(op(ivar)),op(subsop(1=ivar[ii],ldvar)))* diff(xi[ivar[ii]](lvar),op(1,ldvar)),ii=1..nops(ivar))]: diff(zeta_km1,op(1,ldvar)) - convert(result,`+`) fi end: `indiff/desolv/gendl` := proc(ivar,dvar,ordder) local xvar,deplist,lder1,lder2,tlder,ii: deplist := NULL: lder1 := map(proc(fname,lvar) fname(op(lvar)) end,dvar,ivar): for ii to ordder do lder2 := NULL: for xvar in ivar do tlder := map(diff,lder1,xvar): tlder := select(proc(a,l) if member(a,l) then false else true fi end, tlder,[lder2]): lder2 := lder2,op(tlder) od: deplist := deplist,op(lder1): lder1 := [lder2] od: [op(ivar),deplist,op(lder1)] end: proldata:=proc(lambda,xi,eta) local nvars,j, varslist,p; global vars,ukns,XiPhis,GroupP; member(lambda,GroupP,'p'); varslist:=(op(vars),op(ukns)); nvars:=nops(vars); for j to nvars do xi[op(j,vars)]:=subs(DUMMY1=p,DUMMY2=j, proc(varslist) XiPhis[DUMMY1,DUMMY2] end): od; for j to nops(ukns) do eta[op(j,ukns)]:=subs(DUMMY1=p,DUMMY2=nvars+j,DUMMY3=varslist, proc(DUMMY3) XiPhis[DUMMY1,DUMMY2] end): od; op(xi);op(eta); end: #Written by E.Mansfield sortanddiffle:=proc(F,allvars,mytermorder,C) # F a list or set of idp's # N1 is supposed to be maximum of the deg(Neqs)+1 # C an output name local n,nset,fdeg,m,j,ftmp,probset,endset,N1,lKuZ,KuZ; global vars,DegDiffTable,Pout; #readlib('degofdiff',``.`indiff`.`/`.`code`.`/`.`degofdiff.m`); #readlib('Idiff',``.`indiff`.`/`.`code`.`/`.`Idiff.m`); #readlib('noKfactors',``.`indiff`.`/`.`code`.`/`.`noKfactors.m`); if nargs>4 then lKuZ:=true; KuZ:=args[5] fi; endset:={};Pout:={}; probset:=convert(numer(F),set); N1:=max(op(map(rhs,op(op(DegDiffTable)))))+1; while probset<>{} do for n from 0 to N1 do nset[n]:={} od; for n to nops(probset) do degofdiff(probset[n],allvars,mytermorder,'fdeg'); nset[op(2,fdeg)]:=nset[op(2,fdeg)] union {probset[n]} od; probset:={}; for n from 0 to N1-1 do if nset[n]={} then next fi; for m in nset[n] do for j to nops(vars) do ftmp:=numer(Idiff(m,j)); if lKuZ=true then ftmp:=noKfactors([ftmp],KuZ)[1] fi; ftmp:=simplify(ftmp,F); if ftmp=0 then next fi; degofdiff(ftmp,allvars,mytermorder,'fdeg'); if op(2,fdeg)>DegDiffTable[op(1,fdeg)]+1 then Pout:=Pout union {ftmp}; next fi; if op(2,fdeg)<=n then probset:=probset union {ftmp}; next; fi; nset[op(2,fdeg)]:=nset[op(2,fdeg)] union {ftmp}; od; od; od; endset:=endset union {seq(op(nset[j]),j=0..N1)}; od; C:=endset union Pout; end: termorders:=proc() {'mrlex','mtdeg','ttdeg','tdeg','lex','alex','rlex'}; end: #save termorders,``.`indiff`.`/`.`code`.`/`.`termorders.m`: IndiffSimp:=proc(f) local Xset,Uset,Simpf; global Neqs,vars,ukns,HNIlist; Uset:=select(has,Neqs,ukns); Xset:=convert(Neqs,set) minus convert(Uset,set); Simpf:=f; if Uset<>[] then Simpf:=simplify(Simpf,Uset,HNIlist) fi; if Xset<>{} then Simpf:=simplify(Simpf,Xset) fi; Simpf; end: