Initial population.

1 files changed, 255 insertions, 0 deletions

diff --git a/src/algebra/exprode.spad.pamphlet b/src/algebra/exprode.spad.pamphlet
new file mode 100644
index 00000000..3aa24347
--- /dev/null
+++ b/src/algebra/exprode.spad.pamphlet
@@ -0,0 +1,255 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra exprode.spad}
+\author{Manuel Bronstein}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package EXPRODE ExpressionSpaceODESolver}
+<<package EXPRODE ExpressionSpaceODESolver>>=
+)abbrev package EXPRODE ExpressionSpaceODESolver
+++ Taylor series solutions of ODE's
+++ Author: Manuel Bronstein
+++ Date Created: 5 Mar 1990
+++ Date Last Updated: 30 September 1993
+++ Description: Taylor series solutions of explicit ODE's;
+++ Keywords: differential equation, ODE, Taylor series
+ExpressionSpaceODESolver(R, F): Exports == Implementation where
+  R: Join(OrderedSet, IntegralDomain, ConvertibleTo InputForm)
+  F: FunctionSpace R
+
+  K   ==> Kernel F
+  P   ==> SparseMultivariatePolynomial(R, K)
+  OP  ==> BasicOperator
+  SY  ==> Symbol
+  UTS ==> UnivariateTaylorSeries(F, x, center)
+  MKF ==> MakeUnaryCompiledFunction(F, UTS, UTS)
+  MKL ==> MakeUnaryCompiledFunction(F, List UTS, UTS)
+  A1  ==> AnyFunctions1(UTS)
+  AL1 ==> AnyFunctions1(List UTS)
+  EQ  ==> Equation F
+  ODE ==> UnivariateTaylorSeriesODESolver(F, UTS)
+
+  Exports ==> with
+    seriesSolve: (EQ, OP, EQ, EQ) -> Any
+      ++ seriesSolve(eq,y,x=a, y a = b) returns a Taylor series solution
+      ++ of eq around x = a with initial condition \spad{y(a) = b}.
+      ++ Note: eq must be of the form
+      ++ \spad{f(x, y x) y'(x) + g(x, y x) = h(x, y x)}.
+    seriesSolve: (EQ, OP, EQ, List F) -> Any
+      ++ seriesSolve(eq,y,x=a,[b0,...,b(n-1)]) returns a Taylor series
+      ++ solution of eq around \spad{x = a} with initial conditions
+      ++ \spad{y(a) = b0}, \spad{y'(a) = b1},
+      ++ \spad{y''(a) = b2}, ...,\spad{y(n-1)(a) = b(n-1)}
+      ++ eq must be of the form
+      ++ \spad{f(x, y x, y'(x),..., y(n-1)(x)) y(n)(x) +
+      ++ g(x,y x,y'(x),...,y(n-1)(x)) = h(x,y x, y'(x),..., y(n-1)(x))}.
+    seriesSolve: (List EQ, List OP, EQ, List EQ) -> Any
+      ++ seriesSolve([eq1,...,eqn],[y1,...,yn],x = a,[y1 a = b1,...,yn a = bn])
+      ++ returns a taylor series solution of \spad{[eq1,...,eqn]} around
+      ++ \spad{x = a} with initial conditions \spad{yi(a) = bi}.
+      ++ Note: eqi must be of the form
+      ++ \spad{fi(x, y1 x, y2 x,..., yn x) y1'(x) +
+      ++ gi(x, y1 x, y2 x,..., yn x) = h(x, y1 x, y2 x,..., yn x)}.
+    seriesSolve: (List EQ, List OP, EQ, List F) -> Any
+      ++ seriesSolve([eq1,...,eqn], [y1,...,yn], x=a, [b1,...,bn])
+      ++ is equivalent to
+      ++ \spad{seriesSolve([eq1,...,eqn], [y1,...,yn], x = a,
+      ++ [y1 a = b1,..., yn a = bn])}.
+    seriesSolve: (List F, List OP, EQ, List F) -> Any
+      ++ seriesSolve([eq1,...,eqn], [y1,...,yn], x=a, [b1,...,bn])
+      ++ is equivalent to
+      ++ \spad{seriesSolve([eq1=0,...,eqn=0], [y1,...,yn], x=a, [b1,...,bn])}.
+    seriesSolve: (List F, List OP, EQ, List EQ) -> Any
+      ++ seriesSolve([eq1,...,eqn], [y1,...,yn], x = a,[y1 a = b1,..., yn a = bn])
+      ++ is equivalent to
+      ++ \spad{seriesSolve([eq1=0,...,eqn=0], [y1,...,yn], x = a,
+      ++ [y1 a = b1,..., yn a = bn])}.
+    seriesSolve: (EQ, OP, EQ, F) -> Any
+      ++ seriesSolve(eq,y, x=a, b) is equivalent to
+      ++ \spad{seriesSolve(eq, y, x=a, y a = b)}.
+    seriesSolve: (F, OP, EQ, F) -> Any
+      ++ seriesSolve(eq, y, x = a, b) is equivalent to
+      ++ \spad{seriesSolve(eq = 0, y, x = a, y a = b)}.
+    seriesSolve: (F, OP, EQ, EQ) -> Any
+      ++ seriesSolve(eq, y, x = a, y a = b) is equivalent to
+      ++ \spad{seriesSolve(eq=0, y, x=a, y a = b)}.
+    seriesSolve: (F, OP, EQ, List F) -> Any
+      ++ seriesSolve(eq, y, x = a, [b0,...,bn]) is equivalent to
+      ++ \spad{seriesSolve(eq = 0, y, x = a, [b0,...,b(n-1)])}.
+
+  Implementation ==> add
+    checkCompat: (OP, EQ, EQ) -> F
+    checkOrder1: (F, OP, K, SY, F) -> F
+    checkOrderN: (F, OP, K, SY, F, NonNegativeInteger) -> F
+    checkSystem: (F, List K, List F) -> F
+    div2exquo  : F -> F
+    smp2exquo  : P -> F
+    k2exquo    : K -> F
+    diffRhs    : (F, F) -> F
+    diffRhsK   : (K, F) -> F
+    findCompat : (F, List EQ) -> F
+    findEq     : (K, SY, List F) -> F
+    localInteger: F -> F
+
+    opelt := operator("elt"::Symbol)$OP
+    --opex  := operator("exquo"::Symbol)$OP
+    opex  := operator("fixedPointExquo"::Symbol)$OP
+    opint := operator("integer"::Symbol)$OP
+
+    Rint? := R has IntegerNumberSystem
+
+    localInteger n == (Rint? => n; opint n)
+    diffRhs(f, g) == diffRhsK(retract(f)@K, g)
+
+    k2exquo k ==
+      is?(op := operator k, "%diff"::Symbol) =>
+        error "Improper differential equation"
+      kernel(op, [div2exquo f for f in argument k]$List(F))
+
+    smp2exquo p ==
+      map(k2exquo,#1::F,p)$PolynomialCategoryLifting(IndexedExponents K,
+                                                             K, R, P, F)
+
+    div2exquo f ==
+--      one?(d := denom f) => f
+      ((d := denom f) = 1) => f
+      opex(smp2exquo numer f, smp2exquo d)
+
+-- if g is of the form a * k + b, then return -b/a
+    diffRhsK(k, g) ==
+      h := univariate(g, k)
+      (degree(numer h) <= 1) and ground? denom h =>
+        - coefficient(numer h, 0) / coefficient(numer h, 1)
+      error "Improper differential equation"
+
+    checkCompat(y, eqx, eqy) ==
+      lhs(eqy) =$F y(rhs eqx) => rhs eqy
+      error "Improper initial value"
+
+    findCompat(yx, l) ==
+      for eq in l repeat
+        yx =$F lhs eq => return rhs eq
+      error "Improper initial value"
+
+    findEq(k, x, sys) ==
+      k := retract(differentiate(k::F, x))@K
+      for eq in sys repeat
+        member?(k, kernels eq) => return eq
+      error "Improper differential equation"
+
+    checkOrder1(diffeq, y, yx, x, sy) ==
+      div2exquo subst(diffRhs(differentiate(yx::F,x),diffeq),[yx],[sy])
+
+    checkOrderN(diffeq, y, yx, x, sy, n) ==
+      zero? n => error "No initial value(s) given"
+      m     := (minIndex(l := [retract(f := yx::F)@K]$List(K)))::F
+      lv    := [opelt(sy, localInteger m)]$List(F)
+      for i in 2..n repeat
+        l  := concat(retract(f := differentiate(f, x))@K, l)
+        lv := concat(opelt(sy, localInteger(m := m + 1)), lv)
+      div2exquo subst(diffRhs(differentiate(f, x), diffeq), l, lv)
+
+    checkSystem(diffeq, yx, lv) ==
+      for k in kernels diffeq repeat
+        is?(k, "%diff"::SY) =>
+          return div2exquo subst(diffRhsK(k, diffeq), yx, lv)
+      0
+
+    seriesSolve(l:List EQ, y:List OP, eqx:EQ, eqy:List EQ) ==
+      seriesSolve([lhs deq - rhs deq for deq in l]$List(F), y, eqx, eqy)
+
+    seriesSolve(l:List EQ, y:List OP, eqx:EQ, y0:List F) ==
+      seriesSolve([lhs deq - rhs deq for deq in l]$List(F), y, eqx, y0)
+
+    seriesSolve(l:List F, ly:List OP, eqx:EQ, eqy:List EQ) ==
+      seriesSolve(l, ly, eqx,
+                  [findCompat(y rhs eqx, eqy) for y in ly]$List(F))
+
+    seriesSolve(diffeq:EQ, y:OP, eqx:EQ, eqy:EQ) ==
+      seriesSolve(lhs diffeq - rhs diffeq, y, eqx, eqy)
+
+    seriesSolve(diffeq:EQ, y:OP, eqx:EQ, y0:F) ==
+      seriesSolve(lhs diffeq - rhs diffeq, y, eqx, y0)
+
+    seriesSolve(diffeq:EQ, y:OP, eqx:EQ, y0:List F) ==
+      seriesSolve(lhs diffeq - rhs diffeq, y, eqx, y0)
+
+    seriesSolve(diffeq:F, y:OP, eqx:EQ, eqy:EQ) ==
+      seriesSolve(diffeq, y, eqx, checkCompat(y, eqx, eqy))
+
+    seriesSolve(diffeq:F, y:OP, eqx:EQ, y0:F) ==
+      x      := symbolIfCan(retract(lhs eqx)@K)::SY
+      sy     := name y
+      yx     := retract(y lhs eqx)@K
+      f      := checkOrder1(diffeq, y, yx, x, sy::F)
+      center := rhs eqx
+      coerce(ode1(compiledFunction(f, sy)$MKF, y0)$ODE)$A1
+
+    seriesSolve(diffeq:F, y:OP, eqx:EQ, y0:List F) ==
+      x      := symbolIfCan(retract(lhs eqx)@K)::SY
+      sy     := new()$SY
+      yx     := retract(y lhs eqx)@K
+      f      := checkOrderN(diffeq, y, yx, x, sy::F, #y0)
+      center := rhs eqx
+      coerce(ode(compiledFunction(f, sy)$MKL, y0)$ODE)$A1
+
+    seriesSolve(sys:List F, ly:List OP, eqx:EQ, l0:List F) ==
+      x      := symbolIfCan(kx := retract(lhs eqx)@K)::SY
+      fsy    := (sy := new()$SY)::F
+      m      := (minIndex(l0) - 1)::F
+      yx     := concat(kx, [retract(y lhs eqx)@K for y in ly]$List(K))
+      lelt   := [opelt(fsy, localInteger(m := m+1)) for k in yx]$List(F)
+      sys    := [findEq(k, x, sys) for k in rest yx]
+      l      := [checkSystem(eq, yx, lelt) for eq in sys]$List(F)
+      center := rhs eqx
+      coerce(mpsode(l0,[compiledFunction(f,sy)$MKL for f in l])$ODE)$AL1
+
+@
+\section{License}
+<<license>>=
+--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+--All rights reserved.
+--
+--Redistribution and use in source and binary forms, with or without
+--modification, are permitted provided that the following conditions are
+--met:
+--
+--    - Redistributions of source code must retain the above copyright
+--      notice, this list of conditions and the following disclaimer.
+--
+--    - Redistributions in binary form must reproduce the above copyright
+--      notice, this list of conditions and the following disclaimer in
+--      the documentation and/or other materials provided with the
+--      distribution.
+--
+--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
+--      names of its contributors may be used to endorse or promote products
+--      derived from this software without specific prior written permission.
+--
+--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+@
+<<*>>=
+<<license>>
+
+<<package EXPRODE ExpressionSpaceODESolver>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}