aboutsummaryrefslogtreecommitdiff
path: root/src/algebra/exprode.spad.pamphlet
diff options
context:
space:
mode:
authordos-reis <gdr@axiomatics.org>2007-08-14 05:14:52 +0000
committerdos-reis <gdr@axiomatics.org>2007-08-14 05:14:52 +0000
commitab8cc85adde879fb963c94d15675783f2cf4b183 (patch)
treec202482327f474583b750b2c45dedfc4e4312b1d /src/algebra/exprode.spad.pamphlet
downloadopen-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz
Initial population.
Diffstat (limited to 'src/algebra/exprode.spad.pamphlet')
-rw-r--r--src/algebra/exprode.spad.pamphlet255
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}