Changes http://wiki.axiom-developer.org/SsolveSpad/diff
--
\begin{spad}
-- Copyright (C) 2006 Martin Rubey <[EMAIL PROTECTED]>
--
-- This program is free software; you can redistribute it and/or
-- modify it under the terms of the GNU General Public License as
-- published by the Free Software Foundation; either version 2 of
-- the License, or (at your option) any later version.
--
-- This program is distributed in the hope that it will be
-- useful, but WITHOUT ANY WARRANTY; without even the implied
-- warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
-- PURPOSE. See the GNU General Public License for more details.
)abbrev domain SUPEXPR SparseUnivariatePolynomialExpressions
SparseUnivariatePolynomialExpressions(R: Ring): Exports == Implementation where
Exports == UnivariatePolynomialCategory R with
if R has TranscendentalFunctionCategory
then TranscendentalFunctionCategory
Implementation == SparseUnivariatePolynomial R add
if R has TranscendentalFunctionCategory then
exp(p: %): % ==
ground? p => coerce exp ground p
output(hconcat("exp p for p= ", p::OutputForm))$OutputPackage
error "SUPTRAFUN: exp only defined for elements of the
coefficient ring"
sin(p: %): % ==
ground? p => coerce sin ground p
output(hconcat("sin p for p= ", p::OutputForm))$OutputPackage
error "SUPTRAFUN: sin only defined for elements of the
coefficient ring"
asin(p: %): % ==
ground? p => coerce asin ground p
output(hconcat("asin p for p= ", p::OutputForm))$OutputPackage
error "SUPTRAFUN: asin only defined for elements of the
coefficient ring"
cos(p: %): % ==
ground? p => coerce cos ground p
output(hconcat("cos p for p= ", p::OutputForm))$OutputPackage
error "SUPTRAFUN: cos only defined for elements of the
coefficient ring"
acos(p: %): % ==
ground? p => coerce acos ground p
output(hconcat("acos p for p= ", p::OutputForm))$OutputPackage
error "SUPTRAFUN: acos only defined for elements of the
coefficient ring"
)abbrev package UTSSOL TaylorSolve
TaylorSolve(F, UTSF, UTSSUPF): Exports == Implementation where
F: Field
SUP ==> SparseUnivariatePolynomialExpressions
UTSF: UnivariateTaylorSeriesCategory F
UTSSUPF: UnivariateTaylorSeriesCategory SUP F
NNI ==> NonNegativeInteger
Exports == with
seriesSolve: (UTSSUPF -> UTSSUPF, List F) -> UTSF
Implementation == add
seriesSolve(f, l) ==
c1 := map(#1::(SUP F), l)$ListFunctions2(F, SUP F)::(Stream SUP F)
coeffs: Stream SUP F := concat(c1, generate(monomial(1$F,1$NNI)))
-- coeffs: Stream SUP F := concat(c1, monomial(1$F,1$NNI))
st: List Stream SUP F := [coeffs, coeffs]
next: () -> F :=
nr := st.1
res: F
if ground?(coeff: SUP F := nr.1)$(SUP F)
then
res := ground coeff
st.1 := rest nr
else
ns := st.2
eqs: Stream SUP F := coefficients f series ns
while zero? first eqs repeat eqs := rest eqs
eq: SUP F := first eqs
if degree eq > 1 then
if monomial? eq then res := 0
else
output(hconcat("The equation is: ", eq::OutputForm))
$OutputPackage
error "seriesSolve: equation for coefficient not
linear"
else res := (-coefficient(eq, 0$NNI)$(SUP F)
/coefficient(eq, 1$NNI)$(SUP F))
nr.1 := res::SUP F
-- concat!(st.2, monomial(1$F,1$NNI))
st.1 := rest nr
res
series generate next
)abbrev package EXPRSOL ExpressionSolve
ExpressionSolve(R, F, UTSF, UTSSUPF): Exports == Implementation where
R: Join(OrderedSet, IntegralDomain, ConvertibleTo InputForm)
F: FunctionSpace R
UTSF: UnivariateTaylorSeriesCategory F
SUP ==> SparseUnivariatePolynomialExpressions
UTSSUPF: UnivariateTaylorSeriesCategory SUP F
OP ==> BasicOperator
SY ==> Symbol
NNI ==> NonNegativeInteger
MKF ==> MakeBinaryCompiledFunction(F, UTSSUPF, UTSSUPF, UTSSUPF)
Exports == with
seriesSolve: (F, OP, SY, List F) -> UTSF
Implementation == add
opelt := operator("elt"::Symbol)$OP
opdiff := operator("D"::Symbol)$OP
opcoerce := operator("coerce"::Symbol)$OP
replaceDiffs: (F, OP, Symbol) -> F
replaceDiffs (expr, op, sy) ==
lk := kernels expr
for k in lk repeat
if freeOf?(coerce k, sy) then
expr := subst(expr, [k], [opcoerce [coerce k]])
if is?(k, op) then
arg := first argument k
if arg = sy::F
then expr := subst(expr, [k], [(name op)::F])
else expr := subst(expr, [k], [opelt [(name op)::F,
replaceDiffs(arg, op,
sy)]])
-- => "iterate"
if is?(k, %diff) then
args := argument k
differentiand := replaceDiffs(subst(args.1, args.2 =
args.3), op, sy)
expr := subst(expr, [k], [opdiff differentiand])
-- => "iterate"
expr
seriesSolve(expr, op, sy, l) ==
ex := replaceDiffs(expr, op, sy)
f := compiledFunction(ex, name op, sy)$MKF
seriesSolve(f(#1, monomial(1,1)$UTSSUPF), l)$TaylorSolve(F, UTSF,
UTSSUPF)
\end{spad}
--
forwarded from http://wiki.axiom-developer.org/[EMAIL PROTECTED]
