|
DifferentialEquation (DE)
Basic domain for differential equations.
axiom )lib JBC JBC-
JetBundleCategory is now explicitly exposed in frame initial
JetBundleCategory will be automatically loaded when needed from
/var/zope2/var/LatexWiki/JBC.NRLIB/code
JetBundleCategory& is now explicitly exposed in frame initial
JetBundleCategory& will be automatically loaded when needed from
/var/zope2/var/LatexWiki/JBC-.NRLIB/code
axiom )lib JBFC JBFC-
JetBundleFunctionCategory is now explicitly exposed in frame initial
JetBundleFunctionCategory will be automatically loaded when needed
from /var/zope2/var/LatexWiki/JBFC.NRLIB/code
JetBundleFunctionCategory& is now explicitly exposed in frame
initial
JetBundleFunctionCategory& will be automatically loaded when needed
from /var/zope2/var/LatexWiki/JBFC-.NRLIB/code
axiom )lib SEM
SparseEchelonMatrix is now explicitly exposed in frame initial
SparseEchelonMatrix will be automatically loaded when needed from
/var/zope2/var/LatexWiki/SEM.NRLIB/code
axiom )lib VF DIFF
VectorField is now explicitly exposed in frame initial
VectorField will be automatically loaded when needed from
/var/zope2/var/LatexWiki/VF.NRLIB/code
Differential is now explicitly exposed in frame initial
Differential will be automatically loaded when needed from
/var/zope2/var/LatexWiki/DIFF.NRLIB/code
spad )abb domain DE DifferentialEquation
B ==> Boolean
Sy ==> Symbol
I ==> Integer
PI ==> PositiveInteger
NNI ==> NonNegativeInteger
EQ ==> Equation
L ==> List
V ==> Vector
M ==> Matrix
VD ==> V D
MD ==> M D
OUT ==> OutputForm
JBC ==> JetBundleCategory
JBFC ==> JetBundleFunctionCategory JB
SEM ==> SparseEchelonMatrix(JB,D)
DIFF ==> Differential(JB,D)
MVREC ==> Record(Rank:NNI, NumMultVar:NNI, Betas:L NNI)
SREC ==> Record(SDe:$,IC:L D)
-- ------------------------- --
-- DifferentialEquation (DE) --
-- ------------------------- --
++ Description:
++ \axiomType{DifferentialEquation} provides the basic data structures and
++ procedures for differential equations as needed in the geometric theory.
++ Differential equation means here always a submanifold in the jet bundle.
++ The concrete equations which define this submanifold are called system.
++ In an object of the type \axiomType{DifferentialEquation} much more than
++ only the system is stored. \axiom{D} denotes the class of functions allowed
++ as equations. It is assumed that the \axiom{simplify} procedure of \axiom{D}
++ returns only independent equations and a system with symbol in row echelon
++ form.
DifferentialEquation(JB:JBC,D:JBFC) : Cat == Def where
Cat ==> with
order : % -> NNI
++ \axiom{order(de)} yields the order of the differential equation
++ \axiom{de}.
coerce : % -> OUT
++ \axiom{coerce(de)} transforms the differential equation \axiom{de}
++ to \axiomType{OutputForm}.
printSys : L D -> OUT
++ \axiom{printSys(sys)} writes a list of functions as a vector of
++ equations (with right hand side 0) and coerces the result to
++ \axiomType{OutputForm}.
display : % -> Void
++ \axiom{display(de)} prints all information stored about the
++ differential equation \axiom{de}. This comprises the system
++ ordered by the order of the equations, the Jacobi matrices
++ separately for each order and the index of the independent
++ variable with respect to which the equation was lastly
++ differentiated (1 for not prolonged equations).
copy : % -> %
++ \axiom{copy(De)} returns a copy of the equation \axiom{De}.
retract : % -> L D
++ \axiom{retract(de)} returns the system defining the differential
++ equation \axiom{de}.
jacobiMatrix : % -> L SEM
++ \axiom{jacobiMatrix(De)} returns a list of Jacobi matrices sorted
++ by the order of the equations.
generate : L D -> %
++ \axiom{generate(sys)} generates a differential equation from a system.
join : (%,%) -> %
++ \axiom{join(de1,de2)} combines \axiom{de1} and \axiom{de2}
++ to a single differential equation.
insert : (L D,%) -> %
++ \axiom{insert(sys,de)} adds the system \axiom{sys=0}
++ to the differential equation \axiom{de}.
dimension : (%,NNI) -> NNI
++ \axiom{dimension(de,q)} computes the dimension of the differential
++ equation \axiom{de} as a submanifold of the \axiom{q}-th order jet
++ bundle. The result is correct only, if \axiom{de} is simplified.
setSimpMode : NNI -> NNI
++ \axiom{setSimpMode(i)} sets the flag controlling the used
++ simplifications and returns the old value. Current values are:
++ \axiom{i=0} -> No simplification modulo lower order equations.
++ \axiom{i=1} -> Simplification modulo lower order equations.
++ Default is 0.
simplify : % -> SREC
++ \axiom{simplify(de)} simplifies the equations of each order separately
++ using the procedure \axiom{simplify} from \axiom{D}. Found
++ integrability conditions are also returned separately.
extractSymbol : (%,B) -> SEM
++ \axiom{extractSymbol(de,solved?)} computes the symbol of the differential
++ equation \axiom{de}. If \axiom{solved?} is true, the row echelon form of
++ the symbol is computed at once.
analyseSymbol : SEM -> MVREC
++ \axiom{analyseSymbol(symb)} computes the multiplicative variables of
++ the symbol \axiom{symb}.
prolongSymbol : SEM -> SEM
++ \axiom{prolongSymbol(symb)} prolongs directly the symbol \axiom{symb}.
prolongMV : MVREC -> MVREC
++ \axiom{prolongMV(mv)} calculates the number of multiplicative variables
++ for the prolongation of an involutive symbol.
project : (%,NNI) -> %
++ \axiom{project(de,q)} projects the differential equation \axiom{de}
++ of order higher than \axiom{q} into the \axiom{q}-th order jet bundle.
prolong : % -> SREC
++ \axiom{prolong(de)} prolongs the differential equation \axiom{de}.
++ Additionally the arising integrability conditions are returned.
prolong : (%,NNI) -> SREC
++ \axiom{prolong(de,q)} is like \axiom{prolong(de)}. However, only
++ equations of lower order than \axiom{q} are prolonged.
tableau : (SEM,DIFF) -> SEM
++ \axiom{tableau(symb,chi)} computes the tableau parametrized by a given
++ one-form.
tableau : (SEM,L DIFF) -> SEM
++ \axiom{tableau(symb,lchi)} computes the extended tableau parametrized
++ by a given list of one-forms.
Def ==> add
nn:PI := numIndVar()$JB
mm:PI := numDepVar()$JB
-- global variables for the number of independent and dependent
-- variables in JB
simpMode:NNI := 0 -- global flag for simplification mode
setSimpMode(i:NNI):NNI ==
j := simpMode
simpMode := i
j
adapt(der:L NNI, pro?:L B, dep:Union("failed",L L NNI)) : _
Record(Der:L NNI, Pro?:L B) ==
-- Adapts Deriv and Prolonged? after simplification.
-- Local function.
dep case "failed" => [[1 for i in der], [false for i in der]]
resDer:L NNI := empty
resPro?:L B := empty
for d in dep::L L NNI repeat
if one?(#d) then
resDer := cons(qelt(der, first d), resDer)
resPro? := cons(qelt(pro?, first d), resPro?)
else
j:NNI := reduce(min, [qelt(der,i) for i in d])
b:B := reduce("and", [qelt(pro?,i) for i in d])
resDer := cons(j,resDer)
resPro? := cons(b,resPro?)
[reverse! resDer, reverse! resPro?]
-- -------------- --
-- Representation --
-- -------------- --
SysRec := Record(Eqs:L D, JM:SEM, Deriv:L NNI, Prolonged?:L B,_
Simp?:B, Dim?:B, Dim:NNI)
Rep := Record(Sys:L SysRec, Order:L NNI)
-- The equations in Sys are stored order by order. Order contains for
-- each sublist the order. The first list contains the equations of
-- highest order. The Jacobi matrix for each subsystem is in JM.
-- Deriv tells for each equation the index of the independent variable
-- with resprect to which it was last differentiated. For new equations
-- that is always one. Simp? contains flags whether the subsystems have
-- already been simplified, Dim? whether the dimension at this order has
-- already been computed. Dim contains the dimension. Prolonged? tells
-- whether or not an equation has already been prolonged.
copy(De:%):% ==
newSys:L SysRec := [[copy sys.Eqs, copy sys.JM, copy sys.Deriv,_
copy sys.Prolonged?, sys.Simp?, sys.Dim?, sys.Dim] _
for sys in De.Sys]
newOrd := copy De.Order
[newSys,newOrd]
order(De:%):NNI ==
empty? De.Order => 0
first De.Order
retract(De:%):L D ==
LSys := [sys.Eqs for sys in De.Sys]
reduce(append,LSys,empty)
jacobiMatrix(De:%):L SEM == [sys.JM for sys in De.Sys]
printSys(sys:L D):OUT ==
-- Prints system as "vector" of equations
empty? sys => empty
leq:L EQ D := [eq=0 for eq in sys]
tmp:L OUT := empty
for eq in leq repeat
tmp := cons(eq::OUT, cons(" ",tmp))
vconcat reverse tmp
coerce(De:%):OUT == printSys retract De
display(De:%):Void ==
for sys in De.Sys for ord in De.Order repeat
print(hconcat("Order: ",ord::OUT))$OUT
print(" System:")$OUT
print(hconcat(" ",printSys(sys.Eqs)))$OUT
if sys.Simp? then
print(" (system simplified)")$OUT
if sys.Dim? then
print(hconcat(" Dimension: ",sys.Dim::OUT))$OUT
print(" Jacobi matrix:")$OUT
print(hconcat(" ",sys.JM::OUT))$OUT
print(hconcat(" ",allIndices(sys.JM)::OUT))$OUT
print(" Last derivations:")$OUT
print(hconcat(" ",sys.Deriv::OUT))$OUT
void
-- --------------------------- --
-- Basic Operations on Systems --
-- --------------------------- --
generate2(sys:L D,jm:SEM,der:L NNI):% ==
-- Equations are sorted by order.
-- No check for empty system.
-- Local function.
lord:L NNI := [order first row(jm,i).Indices for i in 1..nrows(jm)]
resOrd := reverse sort removeDuplicates lord
nord := #resOrd
inds := allIndices jm
ljm:L SEM := empty
for q in resOrd repeat
while order(first inds)>q repeat
inds := rest inds
ljm := cons(new(inds,1),ljm)
vsys:V L D := new(nord, empty)
vder:V L NNI := new(nord, empty)
vjm:V SEM := construct reverse! ljm
for eq in reverse sys for i in reverse der for q in reverse lord _
for j in nrows(jm)..1 by -1 repeat
pos := position(q,resOrd)+1-minIndex(resOrd)
if empty? qelt(vsys,pos) then
qsetelt!(vsys,pos,[eq])
setRow!(qelt(vjm,pos),1,row(jm,j))
qsetelt!(vder,pos,[i])
else
qsetelt!(vsys,pos,cons(eq,qelt(vsys,pos)))
consRow!(qelt(vjm,pos),row(jm,j))
qsetelt!(vder,pos,cons(i,qelt(vder,pos)))
--for j in minIndex(vjm)..maxIndex(vjm) repeat
-- elimZeroCols! qelt(vjm,i)
resSys:L SysRec := empty
for ord in resOrd for i in minIndex(vsys).. repeat
rec:SysRec := [qelt(vsys,i), qelt(vjm,i), qelt(vder,i), _
new(#qelt(vder,i),false), false, false, 0]
resSys := cons(rec,resSys)
[reverse! resSys, resOrd]
generate(sys:L D):% ==
empty? sys => [[],[]]
nsys := [numerator eq for eq in sys]
der:L NNI := [1 for eq in nsys]
jm:SEM := jacobiMatrix nsys
generate2(nsys,jm,der)
join(De1:%,De2:%):% ==
cDe1 := copy De1; cDe2 := copy De2
sys1 := cDe1.Sys; sys2 := cDe2.Sys
ord1 := cDe1.Order; ord2 := cDe2.Order
resSys:L SysRec := empty
resOrd:L NNI := empty
while not(empty? ord1 and empty? ord2) repeat
if empty? ord1 then
resSys := concat!(reverse! sys2, resSys)
resOrd := concat!(reverse! ord2, resOrd)
ord2 := empty
else if empty? ord2 then
resSys := concat!(reverse! sys1, resSys)
resOrd := concat!(reverse! ord1, resOrd)
ord1 := empty
else
o1 := first ord1
o2 := first ord2
if o1>o2 then
resSys := cons(first sys1, resSys)
resOrd := cons(o1, resOrd)
sys1 := rest sys1
ord1 := rest ord1
else if o2>o1 then
resSys := cons(first sys2, resSys)
resOrd := cons(o2, resOrd)
sys2 := rest sys2
ord2 := rest ord2
else
rec1 := first sys1; rec2 := first sys2
rec:SysRec := [concat!(rec1.Eqs,rec2.Eqs), join(rec1.JM,rec2.JM),_
concat!(rec1.Deriv,rec2.Deriv), _
concat!(rec1.Prolonged?,rec2.Prolonged?), _
false, false, 0]
resSys := cons(rec,resSys)
resOrd := cons(o1,resOrd)
sys1 := rest sys1; sys2 := rest sys2
ord1 := rest ord1; ord2 := rest ord2
[reverse! resSys, reverse! resOrd]
insert(sys:L D,De:%):% ==
newDe:$ := generate sys
join(De,newDe)
dimension(De:%,q:NNI):NNI ==
-- The dimension is computed order by order
-- Caution: Dim and Dim? are changed in De!
empty? De.Order => dimJ(q)$JB
simp? := true$B
tsys := copy De.Sys
tord := copy De.Order
resSys := empty()$L(SysRec)
while first(tord)>q repeat
resSys := cons(first tsys, resSys)
tsys := rest tsys
tord := rest tord
qq := q::I
res := 0$NNI
for sys in tsys for ord in tord repeat
for j in ord+1..qq repeat
res := res + dimS(j)$JB
qq := ord-1
simp? := simp? and sys.Simp?
if sys.Dim? then
res := res + sys.Dim
else
d := orderDim(sys.Eqs,sys.JM,ord)$D
res := res + d
sys.Dim? := true
sys.Dim := d
resSys := cons(sys,resSys)
if not simp? then
print("***** Warning: system not simplified in dimension")$OUT
if qq>=0 then
res := res + dimJ(qq::NNI)$JB
De.Sys := reverse! resSys
res
-- -------------- --
-- Simplification --
-- -------------- --
simplify(De:%):SREC ==
-- Simplification is performed order by order, starting with the highest
-- order. If equations of lower order are generated, they are moved into
-- the corresponding subsystem.
resSys := empty()$L(SysRec)
resOrd:= empty()$L(NNI)
ICs := empty()$L(D)
cDe := copy De
tsys := cDe.Sys
tord := cDe.Order
AllEqs := empty()$L(D)
if simpMode>0 then
AllEqs := retract cDe
while not empty? tord repeat
q := first tord
sys := first tsys
if sys.Simp? then -- already simplified
resSys := cons(sys,resSys)
resOrd := cons(q,resOrd)
else -- not yet simplified
if simpMode>0 then
while not(empty?(AllEqs) or (order(first AllEqs)<q)) repeat
AllEqs := rest AllEqs
-- simplify equations of highest order
if simpMode>0 then
tmp := simpMod(sys.Eqs,sys.JM,AllEqs)$D
tmp := simplify(tmp.Sys,tmp.JM)$D
else
tmp := simplify(sys.Eqs,sys.JM)$D
newEqs := tmp.Sys
newJM := tmp.JM
ad := adapt(sys.Deriv, sys.Prolonged?, tmp.Depend)
newDer := ad.Der
newPro? := ad.Pro?
-- check for equations of lower order
j:NNI := 0
for eq in newEqs for pro? in newPro? for i in 1.. repeat
o := order first row(newJM,i-j).Indices
o>q => error "order raised in simplify"
if o<q then
ICs := cons(eq,ICs)
j := j+1
pos1 := i-j+1
pos2 := i-j+minIndex(newEqs)
newEqs := delete(newEqs,pos2)
newDer := delete(newDer,pos2)
newPro? := delete(newPro?,pos2)
djm:SEM := extract(newJM,pos1,pos1)
sortedPurge!(djm,order(#1)>o)
deleteRow!(newJM,pos1)
pos := position(o,tord)
if pos>=minIndex(tord) then
rec:SysRec := qelt(tsys,pos)
concat!(rec.Eqs,eq)
appendRow!(rec.JM,row(djm,1))
concat!(rec.Deriv,1)
concat!(rec.Prolonged?,pro?)
rec.Simp? := false
rec.Dim? := false
rec.Dim := 0
qsetelt!(tsys,pos,rec)
else
rec:SysRec := [[eq],djm,[1],[pro?],false,false,0]
hord:L NNI := empty
pos := minIndex(tord)-1
while not empty? tord and first(tord)>o repeat
hord := cons(first tord, hord)
tord := rest tord
pos := pos+1
if empty? tord then
tord := reverse! cons(o,hord)
concat!(tsys,rec)
else
tord := concat!(reverse! hord, cons(o,tord))
tsys := insert!(rec,tsys,pos)
rec:SysRec := [newEqs,newJM,newDer,newPro?,true,false,0]
resSys := cons(rec,resSys)
resOrd := cons(q,resOrd)
tsys := rest tsys
tord := rest tord
-- check for inconsistencies or conditions on independent variables
if zero? q then -- algebraic equations
jm0 := first(resSys).JM
for eq in first(resSys).Eqs for i in 1.. repeat
lj := row(jm0,i).Indices
empty? lj => error "inconsistent system"
u?:B := false
for j in 1..mm until u? repeat
u? := member?(U(j::PI)$JB,lj)
not u? => error "independent variables not independent"
[[reverse! resSys, reverse! resOrd], reverse! ICs]
-- --------------------------- --
-- Prolongation and Projection --
-- --------------------------- --
project(De:%,q:NNI):% ==
cDe := copy De
q>=order De => cDe
resSys := cDe.Sys
resOrd := cDe.Order
check:B := true
while not empty? resOrd and first(resOrd)>q repeat
check := check and first(resSys).Simp?
resSys := rest resSys
resOrd := rest resOrd
if not check then
print("***** Warning: projection of not simplified system")$OUT
[resSys,resOrd]
prolong(De:%):SREC ==
-- Prolonged equations are at once simplified. This yields the
-- integrability conditions.
-- prolong equations of highest order
pEqs:L D := empty
pDer:L NNI := empty
pJV:L L JB := empty
pIC:L D := empty
rec := first De.Sys
q := first De.Order
for eq in rec.Eqs for j in rec.Deriv for k in 1.. repeat
jmeq := extract(rec.JM,k,k)
for i in nn..j by -1 repeat
FDiff := formalDiff2(eq,i::PI,jmeq)
pEqs := cons(FDiff.DPhi,pEqs)
pDer := cons(i::NNI,pDer)
pJV := cons(FDiff.JVars,pJV)
pEqs := reverse! pEqs
pJV := reverse! pJV
pDer := reverse! pDer
pJM := jacobiMatrix(pEqs,pJV)
pRec:SysRec:= [pEqs, pJM, pDer, [false for i in pDer], false, false, 0]
pSys:L SysRec := [pRec]
pOrd:L NNI := [q+1]
-- prolong remaining equations, if necessary
-- this yields additional integrability conditions
lastRec := copy rec
lastRec.Prolonged? := [true for j in rec.Deriv]
lastOrd := q
for rec in rest De.Sys for ord in rest De.Order repeat
pEqs := empty
pDer := empty
pJV := empty
for eq in rec.Eqs for j in rec.Deriv for pro? in rec.Prolonged? _
for k in 1.. repeat
if not pro? then
jmeq := extract(rec.JM,k,k)
for i in nn..j by -1 repeat
FDiff := formalDiff2(eq,i::PI,jmeq)
pEqs := cons(FDiff.DPhi,pEqs)
pDer := cons(i::NNI,pDer)
pJV := cons(FDiff.JVars,pJV)
if empty? pEqs then
pSys := cons(lastRec,pSys)
pOrd := cons(lastOrd,pOrd)
else
pIC := append(pIC,pEqs)
pJM := jacobiMatrix(pEqs,pJV)
if ord+1<lastOrd then
pRec := [pEqs, pJM, pDer, [false for i in pDer], false, false, 0]
pSys := cons(pRec,cons(lastRec,pSys))
pOrd := cons(ord+1,cons(lastOrd,pOrd))
else
pRec := [append(lastRec.Eqs,pEqs), join(lastRec.JM,pJM), _
append(lastRec.Deriv,pDer), _
append(lastRec.Prolonged?, [false for i in pDer]), _
false, false, 0]
pSys := cons(pRec,pSys)
pOrd := cons(lastOrd,pOrd)
lastRec := copy rec
lastRec.Prolonged? := [true for j in rec.Deriv]
lastOrd := ord
pSys := cons(lastRec,pSys)
pOrd := cons(lastOrd,pOrd)
res:$ := [reverse! pSys, reverse! pOrd]
tmp := simplify res
[tmp.SDe, concat!(pIC,tmp.IC)]
prolong(De:%,q:NNI):SREC ==
-- Prolongs only equations of order lower than q.
-- Prolonged equations are at once simplified. This yields the
-- integrability conditions.
cDe := copy De
tsys := cDe.Sys
tord := cDe.Order
pSys:L SysRec := empty
pOrd:L NNI := empty
pIC:L D := empty
while first(tord)>q repeat
pSys := cons(first tsys, pSys)
pOrd := cons(first tord, pOrd)
tsys := rest tsys
tord := rest tord
if not first(tord)=q then
-- prolong equations of highest order (<q)
pEqs:L D := empty
pDer:L NNI := empty
pJV:L L JB := empty
rec := first tsys
ord := first tord
for eq in rec.Eqs for j in rec.Deriv for k in 1.. repeat
jmeq := extract(rec.JM,k,k)
for i in nn..j by -1 repeat
FDiff := formalDiff2(eq,i::PI,jmeq)
pEqs := cons(FDiff.DPhi,pEqs)
pDer := cons(i::NNI,pDer)
pJV := cons(FDiff.JVars,pJV)
pEqs := reverse! pEqs
pJV := reverse! pJV
pDer := reverse! pDer
pJM := jacobiMatrix(pEqs,pJV)
pIC := pEqs
pRec:SysRec:= [pEqs, pJM, pDer, [false for i in pDer], false, false, 0]
pSys := cons(pRec,pSys)
pOrd := cons(ord+1,pOrd)
lastRec := first tsys
lastOrd := first tord
-- prolong remaining equations, if necessary
for rec in rest tsys for ord in rest tord repeat
pEqs:L D := empty
pDer:L NNI := empty
pJV:L L JB := empty
for eq in rec.Eqs for j in rec.Deriv for pro? in rec.Prolonged? _
for k in 1.. repeat
if not pro? then
jmeq := extract(rec.JM,k,k)
for i in nn..j by -1 repeat
FDiff := formalDiff2(eq,i::PI,jmeq)
pEqs := cons(FDiff.DPhi,pEqs)
pDer := cons(i::NNI,pDer)
pJV := cons(FDiff.JVars,pJV)
if empty? pEqs then
pSys := cons(lastRec,pSys)
pOrd := cons(lastOrd,pOrd)
else
pEqs := reverse! pEqs
pJV := reverse! pJV
pDer := reverse! pDer
pJM := jacobiMatrix(pEqs,pJV)
pIC := append(pIC,pEqs)
if (ord+1)<lastOrd then
pRec:SysRec := [pEqs, pJM, pDer, [false for i in pDer], _
false, false, 0]
pSys := cons(pRec,cons(lastRec,pSys))
pOrd := cons(ord+1,cons(lastOrd,pOrd))
else
pRec:SysRec := [concat!(lastRec.Eqs,pEqs), join(lastRec.JM,pJM), _
concat!(lastRec.Deriv,pDer), _
concat!(lastRec.Prolonged?, [false for i in pDer]), _
false, false, 0]
pSys := cons(pRec,pSys)
pOrd := cons(lastOrd,pOrd)
rec.Prolonged? := [true for j in rec.Deriv]
lastRec := rec
lastOrd := ord
pSys := cons(lastRec,pSys)
pOrd := cons(lastOrd,pOrd)
res:$ := [reverse! pSys, reverse! pOrd]
tmp := simplify res
[tmp.SDe, concat!(pIC,tmp.IC)]
-- --------------------- --
-- Operations on Symbols --
-- --------------------- --
extractSymbol(De:%,solved?:B):SEM ==
res := extractSymbol(first(De.Sys).JM)$D
if solved? then
res := rowEchelon(res).Ech
res
analyseSymbol(Symb:SEM):MVREC ==
tmp := rowEchelon Symb
ech := tmp.Ech
pivs := pivots ech
MSum := 0$NNI
BetaI := 0$NNI
LastClass:NNI := nn
LBeta:L NNI := empty
for jv in pivs.Indices repeat
CurClass := class jv
if CurClass=LastClass then
BetaI := BetaI + 1
else
LBeta := cons(BetaI,LBeta)
MSum := MSum + BetaI*LastClass
for k in 2..LastClass-CurClass repeat
LBeta := cons(0,LBeta)
BetaI := 1
LastClass := CurClass
LBeta := cons(BetaI,LBeta)
MSum := MSum + BetaI*LastClass
for k in 2..LastClass repeat
LBeta := cons(0,LBeta)
[tmp.Rank, MSum, LBeta]
prolongSymbol(Symb:SEM):SEM ==
-- Direct prolongation of a symbol without differentiation.
oldInds := allIndices Symb
newInds:L JB := empty
for jv in reverse oldInds repeat
for i in 1..nn repeat
newInds := cons(differentiate(jv,i::PI)::JB, newInds)
newInds := sort!(#2<#1, removeDuplicates! newInds)
res:SEM := new(newInds, nn*nrows Symb)
for j in 1..nrows(Symb) repeat
r := row(Symb,j)
for i in nn..1 by -1 repeat
ninds := [differentiate(jv,i::PI)::JB for jv in r.Indices]
setRow!(res, nn*j-i+1, ninds, r.Entries)
res
prolongMV(mv:MVREC):MVREC ==
oldBeta := reverse mv.Betas
newBeta:L NNI := empty
sum:NNI := 0
rank:NNI := 0
msum:NNI := 0
for beta in oldBeta for k in nn..1 by -1 repeat
sum := sum + beta
rank := rank + sum
msum := msum + k*sum
newBeta := cons(sum,newBeta)
[rank, msum, reverse! newBeta]
-- ---------------------- --
-- Operations on Tableaux --
-- ---------------------- --
power(lc:L D, mu:L NNI, mask:L PI):D ==
-- local function to compute the monomial generated by the
-- one-form with coefficients lc and the multi-index mu.
-- mask tells which coefficients are not zero.
res:D := 1
k:PI := 1
while not empty? mask repeat
while k<first(mask) repeat
mu := rest mu
k := k + 1
res := res * first(lc)**first(mu)
lc := rest lc
mask := rest mask
mu := rest mu
k := k + 1
res
extPower(llc:MD, mu:L NNI, nu:L NNI):D ==
-- the rows of the matrix llc contain the coefficients of the one-forms
-- nu determines the repeated indices; mu the exponents
snu := allRepeated(nu)$JB
rmu := m2r(mu)$JB
q := #first(snu)
res:D := 0
for s in snu repeat
prod:D := 1
for si in s for mi in rmu repeat
prod := prod * qelt(llc,nn-si+1,mi)
res := res + prod
res
tableau(Symb:SEM,chi:DIFF):SEM ==
diffs := differentials chi
last(diffs)>X(nn)$JB => error "illegal differential in tableau"
coeffs := coefficients chi
cinds := [index d for d in diffs]
res:SEM := new([U(i::PI) for i in mm..1 by -1], nrows(Symb))
for k in 1..nrows(Symb) repeat
r := row(Symb,k)
sum:VD := new(mm,0)
for jv in r.Indices for ent in r.Entries repeat
a := index jv
mu := multiIndex jv
qsetelt!(sum, a, qelt(sum,a) + ent*power(coeffs,mu,cinds))
li:L JB := empty
le:L D := empty
for i in 1..mm for s in entries sum repeat
if not zero? s then
li := cons(U(i::PI),li)
le := cons(s,le)
setRow!(res,k,li,le)
res
tableau(Symb:SEM,lchi:L DIFF):SEM ==
q := order first allIndices Symb
inds := variables(q,(nn-#lchi+1)::PI)$JB
mco:MD := new(#lchi,nn,0)
for i in 1.. for chi in lchi repeat
for j in 1..nn repeat
qsetelt!(mco,i,j,coefficient(chi,X(j::PI)$JB))
res:SEM := new(inds,nrows(Symb))
for vv in reverse inds repeat
a := index vv
nu := multiIndex vv
for k in 1..nrows(Symb) repeat
r := row(Symb,k)
s:D := 0
for jv in r.Indices for ent in r.Entries repeat
if index(jv)=a then
mu := multiIndex jv
s := s + ent*extPower(mco,mu,nu)
if not zero? s then
rres := row(res,k)
rres.Indices := cons(vv,rres.Indices)
rres.Entries := cons(s,rres.Entries)
setRow!(res,k,rres)
res
spad Compiling FriCAS source code from file
/var/zope2/var/LatexWiki/7564794135111312932-25px002.spad using
old system compiler.
DE abbreviates domain DifferentialEquation
processing macro definition B ==> Boolean
processing macro definition Sy ==> Symbol
processing macro definition I ==> Integer
processing macro definition PI ==> PositiveInteger
processing macro definition NNI ==> NonNegativeInteger
processing macro definition EQ ==> Equation
processing macro definition L ==> List
processing macro definition V ==> Vector
processing macro definition M ==> Matrix
processing macro definition VD ==> V D
processing macro definition MD ==> M D
processing macro definition OUT ==> OutputForm
processing macro definition JBC ==> JetBundleCategory
processing macro definition JBFC ==> JetBundleFunctionCategory JB
processing macro definition SEM ==> SparseEchelonMatrix(JB,D)
processing macro definition DIFF ==> Differential(JB,D)
processing macro definition MVREC ==> Record(Rank: NNI,NumMultVar: NNI,Betas: L NNI)
processing macro definition SREC ==> Record(SDe: $,IC: L D)
processing macro definition Cat ==> -- the constructor category
processing macro definition Def ==> -- the constructor capsule
------------------------------------------------------------------------
initializing NRLIB DE for DifferentialEquation
compiling into NRLIB DE
compiling exported setSimpMode : NonNegativeInteger -> NonNegativeInteger
Time: 0.12 SEC.
compiling local adapt : (List NonNegativeInteger,List Boolean,Union(failed,List List NonNegativeInteger)) -> Record(Der: List NonNegativeInteger,Pro?: List Boolean)
Time: 0.20 SEC.
compiling exported copy : $ -> $
Time: 0.14 SEC.
compiling exported order : $ -> NonNegativeInteger
Time: 0.01 SEC.
compiling exported retract : $ -> List D
Time: 0.04 SEC.
compiling exported jacobiMatrix : $ -> List SparseEchelonMatrix(JB,D)
Time: 0.02 SEC.
compiling exported printSys : List D -> OutputForm
Time: 0.15 SEC.
compiling exported coerce : $ -> OutputForm
Time: 0.01 SEC.
compiling exported display : $ -> Void
Time: 0.12 SEC.
compiling local generate2 : (List D,SparseEchelonMatrix(JB,D),List NonNegativeInteger) -> $
Time: 0.87 SEC.
compiling exported generate : List D -> $
Time: 0.06 SEC.
compiling exported join : ($,$) -> $
Time: 0.16 SEC.
compiling exported insert : (List D,$) -> $
Time: 0 SEC.
compiling exported dimension : ($,NonNegativeInteger) -> NonNegativeInteger
Time: 0.15 SEC.
compiling exported simplify : $ -> Record(SDe: $,IC: List D)
Time: 0.65 SEC.
compiling exported project : ($,NonNegativeInteger) -> $
Time: 0.04 SEC.
compiling exported prolong : $ -> Record(SDe: $,IC: List D)
Time: 0.39 SEC.
compiling exported prolong : ($,NonNegativeInteger) -> Record(SDe: $,IC: List D)
Time: 0.49 SEC.
compiling exported extractSymbol : ($,Boolean) -> SparseEchelonMatrix(JB,D)
Time: 0.02 SEC.
compiling exported analyseSymbol : SparseEchelonMatrix(JB,D) -> Record(Rank: NonNegativeInteger,NumMultVar: NonNegativeInteger,Betas: List NonNegativeInteger)
Time: 0.05 SEC.
compiling exported prolongSymbol : SparseEchelonMatrix(JB,D) -> SparseEchelonMatrix(JB,D)
Time: 0.07 SEC.
compiling exported prolongMV : Record(Rank: NonNegativeInteger,NumMultVar: NonNegativeInteger,Betas: List NonNegativeInteger) -> Record(Rank: NonNegativeInteger,NumMultVar: NonNegativeInteger,Betas: List NonNegativeInteger)
Time: 0.02 SEC.
compiling local power : (List D,List NonNegativeInteger,List PositiveInteger) -> D
Time: 0.06 SEC.
compiling local extPower : (Matrix D,List NonNegativeInteger,List NonNegativeInteger) -> D
Time: 0.08 SEC.
compiling exported tableau : (SparseEchelonMatrix(JB,D),Differential(JB,D)) -> SparseEchelonMatrix(JB,D)
Time: 0.34 SEC.
compiling exported tableau : (SparseEchelonMatrix(JB,D),List Differential(JB,D)) -> SparseEchelonMatrix(JB,D)
Time: 0.25 SEC.
(time taken in buildFunctor: 0)
;;; *** |DifferentialEquation| REDEFINED
;;; *** |DifferentialEquation| REDEFINED
Time: 0 SEC.
Warnings:
[1] generate: not known that (Ring) is of mode (CATEGORY domain (SIGNATURE order ((NonNegativeInteger) $)) (SIGNATURE coerce ((OutputForm) $)) (SIGNATURE printSys ((OutputForm) (List D))) (SIGNATURE display ((Void) $)) (SIGNATURE copy ($ $)) (SIGNATURE retract ((List D) $)) (SIGNATURE jacobiMatrix ((List (SparseEchelonMatrix JB D)) $)) (SIGNATURE generate ($ (List D))) (SIGNATURE join ($ $ $)) (SIGNATURE insert ($ (List D) $)) (SIGNATURE dimension ((NonNegativeInteger) $ (NonNegativeInteger))) (SIGNATURE setSimpMode ((NonNegativeInteger) (NonNegativeInteger))) (SIGNATURE simplify ((Record (: SDe $) (: IC (List D))) $)) (SIGNATURE extractSymbol ((SparseEchelonMatrix JB D) $ (Boolean))) (SIGNATURE analyseSymbol ((Record (: Rank (NonNegativeInteger)) (: NumMultVar (NonNegativeInteger)) (: Betas (List (NonNegativeInteger)))) (SparseEchelonMatrix JB D))) (SIGNATURE prolongSymbol ((SparseEchelonMatrix JB D) (SparseEchelonMatrix JB D))) (SIGNATURE prolongMV ((Record (: Rank (NonNegativeInteger)) (: NumMultVar (NonNegativeInteger)) (: Betas (List (NonNegativeInteger)))) (Record (: Rank (NonNegativeInteger)) (: NumMultVar (NonNegativeInteger)) (: Betas (List (NonNegativeInteger)))))) (SIGNATURE project ($ $ (NonNegativeInteger))) (SIGNATURE prolong ((Record (: SDe $) (: IC (List D))) $)) (SIGNATURE prolong ((Record (: SDe $) (: IC (List D))) $ (NonNegativeInteger))) (SIGNATURE tableau ((SparseEchelonMatrix JB D) (SparseEchelonMatrix JB D) (Differential JB D))) (SIGNATURE tableau ((SparseEchelonMatrix JB D) (SparseEchelonMatrix JB D) (List (Differential JB D)))))
[2] simplify: Depend has no value
[3] prolong: pEqs has no value
[4] prolong: pJV has no value
[5] prolong: pDer has no value
[6] prolong: DPhi has no value
[7] prolong: JVars has no value
[8] prolong: pSys has no value
[9] prolong: pOrd has no value
[10] prolong: pIC has no value
[11] analyseSymbol: LBeta has no value
[12] analyseSymbol: LastClass has no value
[13] tableau: li has no value
[14] tableau: le has no value
[15] tableau: s has no value
Cumulative Statistics for Constructor DifferentialEquation
Time: 4.51 seconds
finalizing NRLIB DE
Processing DifferentialEquation for Browser database:
--------(order (NNI %))---------
--------(coerce (OUT %))---------
--------(printSys (OUT (L D)))---------
--------(display ((Void) %))---------
--------(copy (% %))---------
--------(retract ((L D) %))---------
--------(jacobiMatrix ((L SEM) %))---------
--------(generate (% (L D)))---------
--------(join (% % %))---------
--------(insert (% (L D) %))---------
--------(dimension (NNI % NNI))---------
--------(setSimpMode (NNI NNI))---------
--------(simplify (SREC %))---------
--------(extractSymbol (SEM % B))---------
--------(analyseSymbol (MVREC SEM))---------
--------(prolongSymbol (SEM SEM))---------
--------(prolongMV (MVREC MVREC))---------
--------(project (% % NNI))---------
--------(prolong (SREC %))---------
--------(prolong (SREC % NNI))---------
--------(tableau (SEM SEM DIFF))---------
--------(tableau (SEM SEM (L DIFF)))---------
--------constructor---------
------------------------------------------------------------------------
DifferentialEquation is now explicitly exposed in frame initial
DifferentialEquation will be automatically loaded when needed from
/var/zope2/var/LatexWiki/DE.NRLIB/code
|