axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / ODESYS.spad

)abbrev package ODESYS SystemODESolver
++ Author: Manuel Bronstein
++ Date Created: 11 June 1991
++ Date Last Updated: 13 April 1994
++ Description: 
++ SystemODESolver provides tools for triangulating
++ and solving some systems of linear ordinary differential equations.

SystemODESolver(F, LO) : SIG == CODE where
  F : Field
  LO : LinearOrdinaryDifferentialOperatorCategory F

  N   ==> NonNegativeInteger
  Z   ==> Integer
  MF  ==> Matrix F
  M   ==> Matrix LO
  V   ==> Vector F
  UF  ==> Union(F, "failed")
  UV  ==> Union(V, "failed")
  REC ==> Record(mat: M, vec: V)
  FSL ==> Record(particular: UF, basis: List F)
  VSL ==> Record(particular: UV, basis: List V)
  SOL ==> Record(particular: F, basis: List F)
  USL ==> Union(SOL, "failed")
  ER  ==> Record(C: MF, g: V, eq: LO, rh: F)

  SIG ==> with

    triangulate : (MF, V) -> Record(A:MF, eqs: List ER)
      ++ triangulate(M,v) returns
      ++ \spad{A,[[C_1,g_1,L_1,h_1],...,[C_k,g_k,L_k,h_k]]}
      ++ such that under the change of variable \spad{y = A z}, the first
      ++ order linear system \spad{D y = M y + v} is uncoupled as
      ++ \spad{D z_i = C_i z_i + g_i} and each \spad{C_i} is a companion
      ++ matrix corresponding to the scalar equation \spad{L_i z_j = h_i}.

    triangulate : (M, V) -> REC
      ++ triangulate(m, v) returns \spad{[m_0, v_0]} such that \spad{m_0}
      ++ is upper triangular and the system \spad{m_0 x = v_0} is equivalent
      ++ to \spad{m x = v}.

    solve : (MF,V,(LO,F)->USL) -> Union(Record(particular:V, basis:MF),"failed")
      ++ solve(m, v, solve) returns \spad{[[v_1,...,v_m], v_p]} such that
      ++ the solutions in \spad{F} of the system \spad{D x = m x + v} are
      ++ \spad{v_p + c_1 v_1 + ... + c_m v_m} where the \spad{c_i's} are
      ++ constants, and the \spad{v_i's} form a basis for the solutions of
      ++ \spad{D x = m x}.
      ++ Argument \spad{solve} is a function for solving a single linear
      ++ ordinary differential equation in \spad{F}.

    solveInField : (M, V, (LO, F) -> FSL) -> VSL
      ++ solveInField(m, v, solve) returns \spad{[[v_1,...,v_m], v_p]} such
      ++ that the solutions in \spad{F} of the system \spad{m x = v} are
      ++ \spad{v_p + c_1 v_1 + ... + c_m v_m} where the \spad{c_i's} are
      ++ constants, and the \spad{v_i's} form a basis for the solutions of
      ++ \spad{m x = 0}.
      ++ Argument \spad{solve} is a function for solving a single linear
      ++ ordinary differential equation in \spad{F}.

  CODE ==> add

    import PseudoLinearNormalForm F

    applyLodo   : (M, Z, V, N) -> F
    applyLodo0  : (M, Z, Matrix F, Z, N) -> F
    backsolve   : (M, V, (LO, F) -> FSL) -> VSL
    firstnonzero: (M, Z) -> Z
    FSL2USL     : FSL -> USL
    M2F         : M -> Union(MF, "failed")

    diff := D()$LO

    solve(mm, v, solve) ==
      rec  := triangulate(mm, v)
      sols:List(SOL) := empty()
      for e in rec.eqs repeat
          (u := solve(e.eq, e.rh)) case "failed" => return "failed"
          sols := concat(u::SOL, sols)
      n := nrows(rec.A)    -- dimension of original vectorspace
      k:N := 0             -- sum of sizes of visited companionblocks
      i:N := 0             -- number of companionblocks
      m:N := 0             -- number of Solutions
      part:V := new(n, 0)
      -- count first the different solutions
      for sol in sols repeat 
        m := m + count((f1:F):Boolean +-> f1 ^= 0, sol.basis)$List(F)
      SolMatrix:MF := new(n, m, 0)
      m := 0
      for sol in reverse_! sols repeat
          i := i+1
          er := rec.eqs.i
          nn := #(er.g)           -- size of active companionblock
          for s in sol.basis repeat
              solVec:V := new(n, 0)
              -- compute corresponding solution base with recursion (24)
              solVec(k+1) := s
              for l in 2..nn repeat solVec(k+l) := diff solVec(k+l-1)
              m := m+1
              setColumn!(SolMatrix, m, solVec)
          -- compute with (24) the corresponding components of the part. sol.
          part(k+1) := sol.particular
          for l in 2..nn repeat part(k+l) := diff part(k+l-1) - (er.g)(l-1)
          k := k+nn
      -- transform these values back to the original system
      [rec.A * part, rec.A * SolMatrix]

    triangulate(m:MF, v:V) ==
      k:N := 0       -- sum of companion-dimensions
      rat := normalForm(m, 1, (f1:F):F +-> - diff f1)
      l   := companionBlocks(rat.R, rat.Ainv * v)
      ler:List(ER) := empty()
      for er in l repeat
        n := nrows(er.C)         -- dimension of this companion vectorspace
        op:LO := 0               -- compute homogeneous equation
        for j in 0..n-1 repeat op := op + monomial((er.C)(n, j + 1), j)
        op := monomial(1, n) - op
        sum:V := new(n::N, 0)    -- compute inhomogen Vector (25)
        for j in 1..n-1 repeat sum(j+1) := diff(sum j) + (er.g) j
        h0:F := 0                 -- compute inhomogenity (26)
        for j in 1..n repeat h0 := h0 - (er.C)(n, j) * sum j
        h0 := h0 + diff(sum n) + (er.g) n
        ler := concat([er.C, er.g, op, h0], ler)
        k := k + n
      [rat.A, ler]

-- like solveInField, but expects a system already triangularized
    backsolve(m, v, solve) ==
      part:V
      r := maxRowIndex m
      offset := minIndex v - (mr := minRowIndex m)
      while r >= mr and every?(zero?, row(m, r))$Vector(LO) repeat r := r - 1
      r < mr => error "backsolve: system has a 0 matrix"
      (c := firstnonzero(m, r)) ^= maxColIndex m =>
        error "backsolve: undetermined system"
      rec := solve(m(r, c), v(r + offset))
      dim := (r - mr + 1)::N
      if (part? := ((u := rec.particular) case F)) then
        part := new(dim, 0)                           -- particular solution
        part(r + offset) :=  u::F
-- hom is the basis for the homogeneous solutions, each column is a solution
      hom:Matrix(F) := new(dim, #(rec.basis), 0)
      for i in minColIndex hom .. maxColIndex hom for b in rec.basis repeat
        hom(r, i) := b
      n:N := 1                 -- number of equations already solved
      while r > mr repeat
        r := r - 1
        c := c - 1
        firstnonzero(m, r) ^= c => error "backsolve: undetermined system"
        degree(eq := m(r, c)) > 0 => error "backsolve: pivot of order > 0"
        a := leadingCoefficient(eq)::F
        if part? then
           part(r + offset) := (v(r + offset) - applyLodo(m, r, part, n)) / a
        for i in minColIndex hom .. maxColIndex hom repeat
          hom(r, i) := - applyLodo0(m, r, hom, i, n)
        n := n + 1
      bas:List(V) := [column(hom,i) for i in minColIndex hom..maxColIndex hom]
      part? => [part, bas]
      ["failed", bas]

    solveInField(m, v, solve) ==
      ((n := nrows m) = ncols m) and
         ((u := M2F(diagonalMatrix [diff for i in 1..n] - m)) case MF) =>
             (uu := solve(u::MF, v, 
               (l1:LO,f2:F):USL +-> FSL2USL solve(l1, f2))) case "failed" =>
                  ["failed", empty()]
             rc := uu::Record(particular:V, basis:MF)
             [rc.particular, [column(rc.basis, i) for i in 1..ncols(rc.basis)]]
      rec := triangulate(m, v)
      backsolve(rec.mat, rec.vec, solve)

    M2F m ==
        mf:MF := new(nrows m, ncols m, 0)
        for i in minRowIndex m .. maxRowIndex m repeat
            for j in minColIndex m .. maxColIndex m repeat
               (u := retractIfCan(m(i, j))@Union(F, "failed")) case "failed" =>
                     return "failed"
               mf(i, j) := u::F
        mf

    FSL2USL rec ==
        rec.particular case "failed" => "failed"
        [rec.particular::F, rec.basis]

-- returns the index of the first nonzero entry in row r of m
    firstnonzero(m, r) ==
      for c in minColIndex m .. maxColIndex m repeat
        m(r, c) ^= 0 => return c
      error "firstnonzero: zero row"

-- computes +/[m(r, i) v(i) for i ranging over the last n columns of m]
    applyLodo(m, r, v, n) ==
      ans:F := 0
      c := maxColIndex m
      cv := maxIndex v
      for i in 1..n repeat
        ans := ans + m(r, c) (v cv)
        c := c - 1
        cv := cv - 1
      ans

-- computes +/[m(r, i) mm(i, c) for i ranging over the last n columns of m]
    applyLodo0(m, r, mm, c, n) ==
      ans := 0
      rr := maxRowIndex mm
      cc := maxColIndex m
      for i in 1..n repeat
        ans := ans + m(r, cc) mm(rr, c)
        cc := cc - 1
        rr := rr - 1
      ans

    triangulate(m:M, v:V) ==
      x := copy m
      w := copy v
      nrows := maxRowIndex x
      ncols := maxColIndex x
      minr  := i := minRowIndex x
      offset := minIndex w - minr
      for j in minColIndex x .. ncols repeat
        if i > nrows then leave x
        rown := minr - 1
        for k in i .. nrows repeat
          if (x(k, j) ^= 0) and ((rown = minr - 1) or
                              degree x(k,j) < degree x(rown,j)) then rown := k
          rown = minr - 1 => "enuf"
          x := swapRows_!(x, i, rown)
          swap_!(w, i + offset, rown + offset)
        for k in i+1 .. nrows | x(k, j) ^= 0 repeat
          l := rightLcm(x(i,j), x(k,j))
          a := rightQuotient(l, x(i, j))
          b := rightQuotient(l, x(k, j))
          -- l = a x(i,j) = b x(k,j)
          for k1 in j+1 .. ncols repeat
            x(k, k1) :=  a * x(i, k1) - b * x(k, k1)
          x(k, j) := 0
          w(k + offset) := a(w(i + offset)) - b(w(k + offset))
        i := i+1
      [x, w]