Edit detail for Jet DifferentialEquation revision 1 of 1

DifferentialEquation (DE)

Basic domain for differential equations.

)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

)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

)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

)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

)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

   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