open-axiom-source 1.4.1+svn~2626-2ubuntu2 / usr / lib / open-axiom / src / algebra / kovacic.spad

--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.
 
-- Compile order for the differential equation solver:
-- oderf.spad  odealg.spad  nlode.spad  nlinsol.spad  riccati.spad
-- kovacic.spad  odeef.spad

)abbrev package KOVACIC Kovacic
++ Author: Manuel Bronstein
++ Date Created: 14 January 1992
++ Date Last Updated: 3 February 1994
++ Description:
++ \spadtype{Kovacic} provides a modified Kovacic's algorithm for
++ solving explicitely irreducible 2nd order linear ordinary
++ differential equations.
++ Keywords: differential equation, ODE
Kovacic(F, UP): Exports == Impl where
  F  : Join(CharacteristicZero, AlgebraicallyClosedField,
            RetractableTo Integer, RetractableTo Fraction Integer)
  UP : UnivariatePolynomialCategory F
 
  RF  ==> Fraction UP
  SUP ==> SparseUnivariatePolynomial RF
  LF  ==> List Record(factor:UP, exponent:Integer)
  LODO==> LinearOrdinaryDifferentialOperator1 RF
 
  Exports ==> with
    kovacic: (RF, RF, RF) -> Union(SUP, "failed")
      ++ kovacic(a_0,a_1,a_2) returns either "failed" or P(u) such that
      ++ \spad{$e^{\int(-a_1/2a_2)} e^{\int u}$} is a solution of
      ++      \spad{a_2 y'' + a_1 y' + a0 y = 0}
      ++ whenever \spad{u} is a solution of \spad{P u = 0}.
      ++ The equation must be already irreducible over the rational functions.
    kovacic: (RF, RF, RF, UP -> Factored UP) -> Union(SUP, "failed")
      ++ kovacic(a_0,a_1,a_2,ezfactor) returns either "failed" or P(u) such
      ++ that \spad{$e^{\int(-a_1/2a_2)} e^{\int u}$} is a solution of
      ++      \spad{$a_2 y'' + a_1 y' + a0 y = 0$}
      ++ whenever \spad{u} is a solution of \spad{P u = 0}.
      ++ The equation must be already irreducible over the rational functions.
      ++ Argument \spad{ezfactor} is a factorisation in \spad{UP},
      ++ not necessarily into irreducibles.
 
  Impl ==> add
    import RationalRicDE(F, UP)
 
    case2       : (RF, LF, UP -> Factored UP) -> Union(SUP, "failed")
    cannotCase2?: LF -> Boolean
 
    kovacic(a0, a1, a2) == kovacic(a0, a1, a2, squareFree)
 
    -- it is assumed here that a2 y'' + a1 y' + a0 y is already irreducible
    -- over the rational functions, i.e. that the associated Riccati equation
    -- does NOT have rational solutions (so we don't check case 1 of Kovacic's
    -- algorithm)
    -- currently only check case 2, not 3
    kovacic(a0, a1, a2, ezfactor) ==
      -- transform first the equation to the form y'' = r y
      -- which makes the Galois group unimodular
      -- this does not change irreducibility over the rational functions
    -- the following is split into 5 lines in order to save a couple of
    -- hours of compile time.
      r:RF := a1**2 
      r := r + 2 * a2 * differentiate a1
      r := r - 2 * a1 * differentiate a2
      r := r - 4 * a0 * a2
      r := r  / (4 * a2**2)
      lf := factors squareFree denom r
      case2(r, lf, ezfactor)
 
    -- this is case 2 of Kovacic's algorithm, i.e. look for a solution
    -- of the associated Riccati equation in a quadratic extension
    -- lf is the squarefree factorisation of denom(r) and is used to
    -- check the necessary condition
    case2(r, lf, ezfactor) ==
      cannotCase2? lf => "failed"
      -- build the symmetric square of the operator L = y'' - r y
      -- which is L2 = y''' - 4 r y' - 2 r' y
      l2:LODO := monomial(1, 3) - monomial(4*r, 1) - 2 * differentiate(r)::LODO
      -- no solution in this case if L2 has no rational solution
      empty?(sol := ricDsolve(l2, ezfactor)) => "failed"
      -- otherwise the defining polynomial for an algebraic solution
      -- of the Ricatti equation associated with L is
      -- u^2 - b u + (1/2 b' + 1/2 b^2 - r) = 0
      -- where b is a rational solution of the Ricatti of L2
      b := first sol
      monomial(1, 2)$SUP - monomial(b, 1)$SUP
                         + ((differentiate(b) + b**2 - 2 * r) / (2::RF))::SUP
 
    -- checks the necessary condition for case 2
    -- returns true if case 2 cannot have solutions
    -- the necessary condition is that there is either a factor with
    -- exponent 2 or odd exponent > 2
    cannotCase2? lf ==
      for rec in lf repeat
        rec.exponent = 2 or (odd?(rec.exponent) and rec.exponent > 2) =>
          return false
      true