Artifact cdc3505a63c699070e160134d91df20a6be41d632b6eaf70b70267dffb921795:
- Executable file
r37/packages/eds/edsexptl.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: 7490) [annotate] [blame] [check-ins using] [more...]
module edsexptl; % Experimental (algebraic mode) operators % Author: David Hartley % These procedures need the other packages loaded during compilation load_package 'xideal; %%%% Characteristic variety, symbol relations and symbol matrix Comment. At present, algebraic routines. endcomment; fluid '(!*varopt !*arbvars xvars!* !*allbranch); symbolic operator indexnames; symbolic procedure indexnames u; begin u := makelist uniqids foreach k in getrlist u collect !*a2k k; apply1('indexrange,{{'equal,gensym(),u}}); return u; end; algebraic procedure symbol_relations(S,name); % S:eds, name:id -> symbol_relations:list of 1-form begin scalar tbl,ix,sys,pis,!*varopt,!*arbvars; pform name(i,j) = 1; tbl := tableau S; ix := indexnames independence S; for i:=1:first length tbl do indexrange !!symbol!!index=i; pis := for i:=1:first length tbl collect foreach j in ix collect name(i,-j); sys := for i:=1:first length tbl join for j:=1:length ix collect (tbl(i,j) - part(pis,i,j)); pis := foreach l in pis join l; sys := first solve(sys,append(cobasis s,pis)); sys := foreach x in sys join if lhs x member pis then {lhs x - rhs x} else {}; return sys; end; algebraic procedure symbol_matrix(S,name); % S:eds, name:id -> symbol_matrix:matrix of 0-form begin scalar sys,wlist,n; pform name(i) = 0,{!!symbol!!pi(i,j),!!symbol!!w(i)}=1; n := first length tableau S; wlist := for i:=1:n collect !!symbol!!w(i); sys := symbol_relations(S,!!symbol!!pi); rl := for i:=1:n join foreach j in indexnames independence S collect make_rule(!!symbol!!pi(i,-j),!!symbol!!w(i)*name(-j)); let rl; % sys := (sys where rl); sys := sys; % write showrules !!symbol!!pi; clearrules rl; matrix !!symbol!!mat(length sys,length wlist); for i:=1:length sys do for j:=1:length wlist do !!symbol!!mat(i,j) := coeffn(part(sys,i),part(wlist,j),1); return !!symbol!!mat; end; algebraic procedure characteristic_variety(S,name); % S:eds, name:id -> characteristic_variety:{list of 0-form,list of % variable} begin scalar ix,m,sys; scalar xvars!*; % make all 0-forms coefficients ix := indexnames independence S; m := symbol_matrix(S,name); if first length m > second length m then m := tp m; for i:=1:second length m do indexrange symbol!!index!!=i; wlist := for i:=1:second length m collect !!symbol!!w(i); www := 1; for i:=1:first length m do www := (for j:=1:length wlist sum m(i,j)*!!symbol!!w(j))^www; return {excoeffs www,foreach i in ix collect name(-i)}; end; algebraic procedure make_rule(lh,rh); lh => rh; %%% Invariants, or first integrals. fluid '(!*edsdebug print_ fname_ time_ xvars!* !*allbranch !*arbvars); mkform!*('eds!:t,0); algebraic procedure edsorderp(x,y); % Just a hook for sort if ordp(x,y) then 1 else 0; put('invariants,'psopfn,'invariants); symbolic procedure invariants u; if length u = 2 then (algebraic invariants0(x,y)) where x=car u, y=cadr u else if length u = 1 then (algebraic invariants0(x,y)) where x=car u, y=makelist nil else rederr "Wrong number of arguments to invariants"; algebraic procedure invariants0(S,C); begin scalar ans,inv,cfrm,Z,xvars!*; load_package odesolve,crack; cfrm := coframing(); if part(S,0) = eds then << set_coframing S; if C = {} then C := coordinates S; S := systemeds S >> % Use systemeds rather than system for % compiler. else S := xauto S; if C = {} then C := reverse sort(coordinates S,edsorderp); Z := for a:=1:length S collect lisp mkform!*(mkid('eds!:u,a),0); ans := foliation(S,C,Z); inv := solve(ans,Z); if length Z = 1 then inv := foreach x in inv collect {x}; if lisp !*edsdebug then write "Constants"; if lisp !*edsdebug then write inv; if length inv neq 1 then rederr "Not a unique solution"; set_coframing cfrm; return foreach x in first inv collect rhs x; end; algebraic procedure foliation(S,C,Z); begin scalar r,n,x,S0,Z0,g,Q,f,f0; scalar print_,fname_,time_,!*allbranch,!*arbvars,xvars!*; % Constants r := length S; n := length C; fname_ := 'eds!:c; % Deal with errors and end case if r > n then rerror(eds,000,"Not enough coordinates in foliation"); if r neq length Z then rerror(eds,000,"Wrong number of invariant labels in foliation"); if r = n then << g := for a:=1:r collect part(C,a) = part(Z,a); lisp edsdebug("Intermediate result",g,'prefix); return g >>; % Choose truncation S0 := {}; Z0 := {}; while length S0 < r do << x := first C; C := rest C; Z0 := x . Z0; S0 := xauto(S xmod {d x}) >>; C := append(C,rest Z0); lisp edsdebug("Truncating coordinate : ",x,'prefix); % Compute foliation for truncation g := foliation(S0,C,Z); % Calculate ODE foreach y in Z do << lisp(y := !*a2k y); fdomain y=y(eds!:t) >>; S := pullback(S,g); S := pullback(S,{x = eds!:t}); Q := foreach f in S collect @eds!:t _| f; Q := solve(Q,foreach y in Z collect @(y,eds!:t)); if r neq 1 then Q := first Q; Q := foreach f in Q collect (lhs f - rhs f); Q := sub(partdf=df,Q); lisp edsdebug("CRACK ODE",Q,'prefix); % Solve ODE Q := crack(Q,{},Z,{}); lisp edsdebug("CRACK solution",Q,'prefix); % Analyse result for the general solution f := {}; while Q neq {} do << f := first Q; Q := rest Q; Z0 := third f; if first f = {} and length Z0 = r then Q := {} else if length Z0 > r then if length(f0 := solve(first f,Z)) = 0 then f := {} else << if r = 1 then f0 := {{first f0}}; Z0 := foreach v in Z0 join if v member Z then {} else {v}; f := {{},append(second f,first f0),Z0}; Q := {} >> else f := {} >>; foreach y in Z do << lisp(y := !*a2k y); remfdomain y >>; if f = {} then rerror(eds,000,"Intermediate ODE general solution not found"); % Compose general solution with truncated foliation g := sub(second f,g); f := (eds!:t = x) . for a := 1:r collect part(Z0,a) = part(Z,a); g := sub(f,g); lisp edsdebug("Intermediate result",g,'prefix); return g; end; %%% Homotopy operator algebraic procedure poincare df; % with df a closed p-form POINCARE returns a (p-1)-form f % satisfying df=d f. begin scalar f; pform !!lambda!! = 0; f := sub(for each c in coordinates df collect c = c * !!lambda!!,df); % f := sub(for each c in allcoords df collect c = c * !!lambda!!,df); f := @(!!lambda!!) _| f; f := int(f,!!lambda!!); f := sub(!!lambda!! = 1,f) - sub(!!lambda!! = 0,f); % if d f - df neq 0 then write "error in POINCARE"; return reval f; end; %%% Integrability conditions put('integrability,'rtypefn,'quotelist); put('integrability,'listfn,'evalintegrability); symbolic procedure evalintegrability(s,v); % s:eds|rlist, v:bool -> evalintegrability:rlist if edsp(s := reval car s) then !*sys2a1(nonpfaffpart eds_sys edscall closure s,v) else algebraic append(s xmod one_forms s,d s xmod one_forms s); endmodule; end;