\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra kovacic.spad} \author{Manuel Bronstein} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package KOVACIC Kovacic} <>= )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 @ \section{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. @ <<*>>= <> -- Compile order for the differential equation solver: -- oderf.spad odealg.spad nlode.spad nlinsol.spad riccati.spad -- kovacic.spad odeef.spad <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}