1)abbrev package DEFINTEF ElementaryFunctionDefiniteIntegration 2++ Definite integration of elementary functions. 3++ Author: Manuel Bronstein 4++ Date Created: 14 April 1992 5++ Description: 6++ \spadtype{ElementaryFunctionDefiniteIntegration} 7++ provides functions to compute definite 8++ integrals of elementary functions. 9ElementaryFunctionDefiniteIntegration(R, F) : Exports == Implementation where 10 R : Join(PolynomialFactorizationExplicit, Comparable, CharacteristicZero, 11 RetractableTo Integer, LinearlyExplicitOver Integer) 12 F : Join(TranscendentalFunctionCategory, PrimitiveFunctionCategory, 13 AlgebraicallyClosedFunctionSpace R) 14 15 B ==> Boolean 16 SE ==> Symbol 17 Z ==> Integer 18 P ==> SparseMultivariatePolynomial(R, K) 19 K ==> Kernel F 20 UP ==> SparseUnivariatePolynomial F 21 OFE ==> OrderedCompletion F 22 U ==> Union(f1:OFE, f2:List OFE, fail:"failed", pole:"potentialPole") 23 24 Exports ==> with 25 integrate : (F, SegmentBinding OFE) -> U 26 ++ integrate(f, x = a..b) returns the integral of 27 ++ \spad{f(x)dx} from a to b. 28 ++ Error: if f has a pole for x between a and b. 29 integrate : (F, SegmentBinding OFE, String) -> U 30 ++ integrate(f, x = a..b, "noPole") returns the 31 ++ integral of \spad{f(x)dx} from a to b. 32 ++ If it is not possible to check whether f has a pole for x 33 ++ between a and b (because of parameters), then this function 34 ++ will assume that f has no such pole. 35 ++ Error: if f has a pole for x between a and b or 36 ++ if the last argument is not "noPole". 37 innerint : (F, SE, OFE, OFE, B) -> U 38 ++ innerint(f, x, a, b, ignore?) should be local but conditional 39 40 Implementation ==> add 41 import from ElementaryFunctionSign(R, F) 42 import from DefiniteIntegrationTools(R, F) 43 import from FunctionSpaceIntegration(R, F) 44 45 polyIfCan : (P, K) -> Union(UP, "failed") 46 int : (F, SE, OFE, OFE, B) -> U 47 nopole : (F, SE, K, OFE, OFE) -> U 48 checkFor0 : (P, K, OFE, OFE) -> Union(B, "failed") 49 checkSMP : (P, SE, K, OFE, OFE) -> Union(B, "failed") 50 checkForPole : (F, SE, K, OFE, OFE) -> Union(B, "failed") 51 posit : (F, SE, K, OFE, OFE) -> Union(B, "failed") 52 negat : (F, SE, K, OFE, OFE) -> Union(B, "failed") 53 moreThan : (OFE, Fraction Z) -> Union(B, "failed") 54 55 if R has Join(ConvertibleTo Pattern Integer, PatternMatchable Integer) 56 and F has SpecialFunctionCategory then 57 import from PatternMatchIntegration(R, F) 58 59 innerint(f, x, a, b, ignor?) == 60 ((u := int(f, x, a, b, ignor?)) case f1) or (u case f2) 61 or ((v := pmintegrate(f, x, a, b)) case "failed") => u 62 [v::F::OFE] 63 64 else 65 innerint(f, x, a, b, ignor?) == int(f, x, a, b, ignor?) 66 67 integrate(f : F, s : SegmentBinding OFE) == 68 innerint(f, variable s, low(segment(s)), high(segment(s)), false) 69 70 integrate(f : F, s : SegmentBinding OFE, str : String) == 71 innerint(f, variable s, low(segment(s)), high(segment(s)), ignore? str) 72 73 int(f, x, a, b, ignor?) == 74 a = b => [0$F::OFE] 75 k := kernel(x)@Kernel(F) 76 (z := checkForPole(f, x, k, a, b)) case "failed" => 77 ignor? => nopole(f, x, k, a, b) 78 ["potentialPole"] 79 z::B => error "integrate: pole in path of integration" 80 nopole(f, x, k, a, b) 81 82 checkForPole(f, x, k, a, b) == 83 ((u := checkFor0(d := denom f, k, a, b)) case "failed") or (u::B) => u 84 ((u := checkSMP(d, x, k, a, b)) case "failed") or (u::B) => u 85 checkSMP(numer f, x, k, a, b) 86 87-- true if p has a zero between a and b exclusive 88 checkFor0(p, x, a, b) == 89 (u := polyIfCan(p, x)) case UP => checkForZero(u::UP, a, b, false) 90 (v := isTimes p) case List(P) => 91 for t in v::List(P) repeat 92 ((w := checkFor0(t, x, a, b)) case "failed") or (w::B) => return w 93 false 94 (z := isExpt p) case "failed" => "failed" 95 k := z.var 96-- functions with no real zeros 97 is?(k, 'exp) or is?(k, 'acot) or is?(k, 'cosh) => false 98-- special case for log 99 is?(k, 'log) => 100 (w := moreThan(b, 1)) case "failed" or not(w::B) => w 101 moreThan(-a, -1) 102 "failed" 103 104-- returns true if a > b, false if a < b, "failed" if can't decide 105 moreThan(a, b) == 106 (r := retractIfCan(a)@Union(F, "failed")) case "failed" => -- infinite 107 whatInfinity(a) > 0 108 (u := retractIfCan(r::F)@Union(Fraction Z, "failed")) case "failed" => 109 "failed" 110 u::Fraction(Z) > b 111 112-- true if p has a pole between a and b 113 checkSMP(p, x, k, a, b) == 114 (u := polyIfCan(p, k)) case UP => false 115 (v := isTimes p) case List(P) => 116 for t in v::List(P) repeat 117 ((w := checkSMP(t, x, k, a, b)) case "failed") or (w::B) => return w 118 false 119 (v := isPlus p) case List(P) => 120 n := 0 -- number of summand having a pole 121 for t in v::List(P) repeat 122 (w := checkSMP(t, x, k, a, b)) case "failed" => return w 123 if w::B then n := n + 1 124 zero? n => false -- no summand has a pole 125 (n = 1) => true -- only one summand has a pole 126 "failed" -- at least 2 summands have a pole 127 (z := isExpt(p)) case "failed" => "failed" 128 kk := z.var 129 -- nullary operators have no poles 130 nullary? operator kk => false 131 ak := argument kk 132 f := first ak 133 -- functions which are defined over all the reals: 134 is?(kk, 'exp) or is?(kk, 'sin) or is?(kk, 'cos) 135 or is?(kk, 'sinh) or is?(kk, 'cosh) or is?(kk, 'tanh) 136 or is?(kk, 'sech) or is?(kk, 'atan) or is?(kk, 'acot) 137 or is?(kk, 'asinh) or is?(kk, 'erf) or is?(kk, 'erfi) 138 or is?(kk, 'fresnelC) or is?(kk, 'fresnelS) 139 or is?(kk, 'Si) or is?(kk, 'Shi) => checkForPole(f, x, k, a, b) 140 -- functions which are defined on (-1, +1): 141 is?(kk, 'asin) or is?(kk, 'acos) or is?(kk, 'atanh) => 142 ((w := checkForPole(f, x, k, a, b)) case "failed") or (w::B) => w 143 ((w := posit(f - 1, x, k, a, b)) case "failed") or (w::B) => w 144 negat(f + 1, x, k, a, b) 145 -- functions which are defined on (+1, +infty): 146 is?(kk, 'acosh) => 147 ((w := checkForPole(f, x, k, a, b)) case "failed") or (w::B) => w 148 negat(f - 1, x, k, a, b) 149 -- functions which are defined on (0, +infty): 150 is?(kk, 'log) or is?(kk, 'nthRoot) or is?(kk, 'Ei) or is?(kk, 'Ci) 151 or is?(kk, 'Chi) or is?(kk, 'dilog) => 152 ((w := checkForPole(f, x, k, a, b)) case "failed") or (w::B) => w 153 negat(f, x, k, a, b) 154 is?(kk, 'Gamma) and #ak = 2 and D(f, x) = 0 => 155 f2 := ak(2) 156 ((w := checkForPole(f2, x, k, a, b)) case "failed") or (w::B) => w 157 negat(f2, x, k, a, b) 158 -- defined on (-infty, 1) 159 is?(kk, 'polylog) and D(f, x) = 0 => 160 f2 := ak(2) 161 ((w := checkForPole(f2, x, k, a, b)) case "failed") or (w::B) => w 162 posit(f - 1, x, k, a, b) 163 "failed" 164 165-- returns true if it is certain that f takes at least one strictly positive 166-- value for x in (a, b), false if it is certain that f takes no strictly 167-- positive value in (a,b), "failed" otherwise 168-- f must be known to have no poles in (a, b) 169 posit(f, x, k, a, b) == 170 z := 171 (r := retractIfCan(a)@Union(F, "failed")) case "failed" => sign(f, x, a) 172 sign(f, x, r::F, "right") 173 (b1 := z case Z) and z::Z > 0 => true 174 z := 175 (r := retractIfCan(b)@Union(F, "failed")) case "failed" => sign(f, x, b) 176 sign(f, x, r::F, "left") 177 (b2 := z case Z) and z::Z > 0 => true 178 b1 and b2 => 179 ((w := checkFor0(numer f, k, a, b)) case "failed") or (w::B) => "failed" 180 false 181 "failed" 182 183-- returns true if it is certain that f takes at least one strictly negative 184-- value for x in (a, b), false if it is certain that f takes no strictly 185-- negative value in (a,b), "failed" otherwise 186-- f must be known to have no poles in (a, b) 187 negat(f, x, k, a, b) == 188 z := 189 (r := retractIfCan(a)@Union(F, "failed")) case "failed" => sign(f, x, a) 190 sign(f, x, r::F, "right") 191 (b1 := z case Z) and z::Z < 0 => true 192 z := 193 (r := retractIfCan(b)@Union(F, "failed")) case "failed" => sign(f, x, b) 194 sign(f, x, r::F, "left") 195 (b2 := z case Z) and z::Z < 0 => true 196 b1 and b2 => 197 ((w := checkFor0(numer f, k, a, b)) case "failed") or (w::B) => "failed" 198 false 199 "failed" 200 201-- returns a UP if p is only a poly w.r.t. the kernel x 202 polyIfCan(p, x) == 203 q := univariate(p, x) 204 ans : UP := 0 205 while q ~= 0 repeat 206 member?(x, tower(c := leadingCoefficient(q)::F)) => return "failed" 207 ans := ans + monomial(c, degree q) 208 q := reductum q 209 ans 210 211-- integrate f for x between a and b assuming that f has no pole in between 212 nopole(f, x, k, a, b) == 213 (u := integrate(f, x)) case F => 214 (v := computeInt(k, u::F, a, b, false)) case "failed" => ["failed"] 215 [v::OFE] 216 ans := empty()$List(OFE) 217 for g in u::List(F) repeat 218 (v := computeInt(k, g, a, b, false)) case "failed" => return ["failed"] 219 ans := concat!(ans, [v::OFE]) 220 [ans] 221 222--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. 223--All rights reserved. 224-- 225--Redistribution and use in source and binary forms, with or without 226--modification, are permitted provided that the following conditions are 227--met: 228-- 229-- - Redistributions of source code must retain the above copyright 230-- notice, this list of conditions and the following disclaimer. 231-- 232-- - Redistributions in binary form must reproduce the above copyright 233-- notice, this list of conditions and the following disclaimer in 234-- the documentation and/or other materials provided with the 235-- distribution. 236-- 237-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the 238-- names of its contributors may be used to endorse or promote products 239-- derived from this software without specific prior written permission. 240-- 241--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 242--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 243--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 244--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 245--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 246--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 247--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 248--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 249--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 250--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 251--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 252