forked from fricas/fricas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
kovacic.spad
124 lines (114 loc) · 5.67 KB
/
kovacic.spad
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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
)abbrev package KOVACIC Kovacic
++ Author: Manuel Bronstein
++ Date Created: 14 January 1992
++ Description:
++ \spadtype{Kovacic} provides a modified Kovacic's algorithm for
++ solving explicitly 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 : NonNegativeInteger)
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 from 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
--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