/usr/share/axiom-20170501/src/algebra/KOVACIC.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 | )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.
Kovacic(F, UP) : SIG == CODE 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
SIG ==> 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.
CODE ==> 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, 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, 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
|