diff options
author | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
---|---|---|
committer | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
commit | ab8cc85adde879fb963c94d15675783f2cf4b183 (patch) | |
tree | c202482327f474583b750b2c45dedfc4e4312b1d /src/algebra/nlode.spad.pamphlet | |
download | open-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz |
Initial population.
Diffstat (limited to 'src/algebra/nlode.spad.pamphlet')
-rw-r--r-- | src/algebra/nlode.spad.pamphlet | 207 |
1 files changed, 207 insertions, 0 deletions
diff --git a/src/algebra/nlode.spad.pamphlet b/src/algebra/nlode.spad.pamphlet new file mode 100644 index 00000000..2f7f672d --- /dev/null +++ b/src/algebra/nlode.spad.pamphlet @@ -0,0 +1,207 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra nlode.spad} +\author{Manuel Bronstein} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{package NODE1 NonLinearFirstOrderODESolver} +<<package NODE1 NonLinearFirstOrderODESolver>>= +)abbrev package NODE1 NonLinearFirstOrderODESolver +++ Author: Manuel Bronstein +++ Date Created: 2 September 1991 +++ Date Last Updated: 14 October 1994 +++ Description: NonLinearFirstOrderODESolver provides a function +++ for finding closed form first integrals of nonlinear ordinary +++ differential equations of order 1. +++ Keywords: differential equation, ODE +NonLinearFirstOrderODESolver(R, F): Exports == Implementation where + R: Join(OrderedSet, EuclideanDomain, RetractableTo Integer, + LinearlyExplicitRingOver Integer, CharacteristicZero) + F: Join(AlgebraicallyClosedFunctionSpace R, TranscendentalFunctionCategory, + PrimitiveFunctionCategory) + + N ==> NonNegativeInteger + Q ==> Fraction Integer + UQ ==> Union(Q, "failed") + OP ==> BasicOperator + SY ==> Symbol + K ==> Kernel F + U ==> Union(F, "failed") + P ==> SparseMultivariatePolynomial(R, K) + REC ==> Record(coef:Q, logand:F) + SOL ==> Record(particular: F,basis: List F) + BER ==> Record(coef1:F, coefn:F, exponent:N) + + Exports ==> with + solve: (F, F, OP, SY) -> U + ++ solve(M(x,y), N(x,y), y, x) returns \spad{F(x,y)} such that + ++ \spad{F(x,y) = c} for a constant \spad{c} is a first integral + ++ of the equation \spad{M(x,y) dx + N(x,y) dy = 0}, or + ++ "failed" if no first-integral can be found. + + Implementation ==> add + import ODEIntegration(R, F) + import ElementaryFunctionODESolver(R, F) -- recursive dependency! + + checkBernoulli : (F, F, K) -> Union(BER, "failed") + solveBernoulli : (BER, OP, SY, F) -> Union(F, "failed") + checkRiccati : (F, F, K) -> Union(List F, "failed") + solveRiccati : (List F, OP, SY, F) -> Union(F, "failed") + partSolRiccati : (List F, OP, SY, F) -> Union(F, "failed") + integratingFactor: (F, F, SY, SY) -> U + + unk := new()$SY + kunk:K := kernel unk + + solve(m, n, y, x) == +-- first replace the operator y(x) by a new symbol z in m(x,y) and n(x,y) + lk:List(K) := [retract(yx := y(x::F))@K] + lv:List(F) := [kunk::F] + mm := eval(m, lk, lv) + nn := eval(n, lk, lv) +-- put over a common denominator (to balance m and n) + d := lcm(denom mm, denom nn)::F + mm := d * mm + nn := d * nn +-- look for an integrating factor mu + (u := integratingFactor(mm, nn, unk, x)) case F => + mu := u::F + mm := mm * mu + nn := nn * mu + eval(int(mm,x) + int(nn-int(differentiate(mm,unk),x), unk),[kunk],[yx]) +-- check for Bernoulli equation + (w := checkBernoulli(m, n, k1 := first lk)) case BER => + solveBernoulli(w::BER, y, x, yx) +-- check for Riccati equation + (v := checkRiccati(m, n, k1)) case List(F) => + solveRiccati(v::List(F), y, x, yx) + "failed" + +-- look for an integrating factor + integratingFactor(m, n, y, x) == +-- check first for exactness + zero?(d := differentiate(m, y) - differentiate(n, x)) => 1 +-- look for an integrating factor involving x only + not member?(y, variables(f := d / n)) => expint(f, x) +-- look for an integrating factor involving y only + not member?(x, variables(f := - d / m)) => expint(f, y) +-- room for more techniques later on (e.g. Prelle-Singer etc...) + "failed" + +-- check whether the equation is of the form +-- dy/dx + p(x)y + q(x)y^N = 0 with N > 1 +-- i.e. whether m/n is of the form p(x) y + q(x) y^N +-- returns [p, q, N] if the equation is in that form + checkBernoulli(m, n, ky) == + r := denom(f := m / n)::F + (not freeOf?(r, y := ky::F)) + or (d := degree(p := univariate(numer f, ky))) < 2 + or degree(pp := reductum p) ^= 1 or reductum(pp) ^= 0 + or (not freeOf?(a := (leadingCoefficient(pp)::F), y)) + or (not freeOf?(b := (leadingCoefficient(p)::F), y)) => "failed" + [a / r, b / r, d] + +-- solves the equation dy/dx + rec.coef1 y + rec.coefn y^rec.exponent = 0 +-- the change of variable v = y^{1-n} transforms the above equation to +-- dv/dx + (1 - n) p v + (1 - n) q = 0 + solveBernoulli(rec, y, x, yx) == + n1 := 1 - rec.exponent::Integer + deq := differentiate(yx, x) + n1 * rec.coef1 * yx + n1 * rec.coefn + sol := solve(deq, y, x)::SOL -- can always solve for order 1 +-- if v = vp + c v0 is the general solution of the linear equation, then +-- the general first integral for the Bernoulli equation is +-- (y^{1-n} - vp) / v0 = c for any constant c + (yx**n1 - sol.particular) / first(sol.basis) + +-- check whether the equation is of the form +-- dy/dx + q0(x) + q1(x)y + q2(x)y^2 = 0 +-- i.e. whether m/n is a quadratic polynomial in y. +-- returns the list [q0, q1, q2] if the equation is in that form + checkRiccati(m, n, ky) == + q := denom(f := m / n)::F + (not freeOf?(q, y := ky::F)) or degree(p := univariate(numer f, ky)) > 2 + or (not freeOf?(a0 := (coefficient(p, 0)::F), y)) + or (not freeOf?(a1 := (coefficient(p, 1)::F), y)) + or (not freeOf?(a2 := (coefficient(p, 2)::F), y)) => "failed" + [a0 / q, a1 / q, a2 / q] + +-- solves the equation dy/dx + l.1 + l.2 y + l.3 y^2 = 0 + solveRiccati(l, y, x, yx) == +-- get first a particular solution + (u := partSolRiccati(l, y, x, yx)) case "failed" => "failed" +-- once a particular solution yp is known, the general solution is of the +-- form y = yp + 1/v where v satisfies the linear 1st order equation +-- v' - (l.2 + 2 l.3 yp) v = l.3 + deq := differentiate(yx, x) - (l.2 + 2 * l.3 * u::F) * yx - l.3 + gsol := solve(deq, y, x)::SOL -- can always solve for order 1 +-- if v = vp + c v0 is the general solution of the above equation, then +-- the general first integral for the Riccati equation is +-- (1/(y - yp) - vp) / v0 = c for any constant c + (inv(yx - u::F) - gsol.particular) / first(gsol.basis) + +-- looks for a particular solution of dy/dx + l.1 + l.2 y + l.3 y^2 = 0 + partSolRiccati(l, y, x, yx) == +-- we first do the change of variable y = z / l.3, which transforms +-- the equation into dz/dx + l.1 l.3 + (l.2 - l.3'/l.3) z + z^2 = 0 + q0 := l.1 * (l3 := l.3) + q1 := l.2 - differentiate(l3, x) / l3 +-- the equation dz/dx + q0 + q1 z + z^2 = 0 is transformed by the change +-- of variable z = w'/w into the linear equation w'' + q1 w' + q0 w = 0 + lineq := differentiate(yx, x, 2) + q1 * differentiate(yx, x) + q0 * yx +-- should be made faster by requesting a particular nonzero solution only + (not((gsol := solve(lineq, y, x)) case SOL)) + or empty?(bas := (gsol::SOL).basis) => "failed" + differentiate(first bas, x) / (l3 * first bas) + +@ +\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>> + +-- Compile order for the differential equation solver: +-- oderf.spad odealg.spad nlode.spad nlinsol.spad riccati.spad odeef.spad + +<<package NODE1 NonLinearFirstOrderODESolver>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |