Thu Feb 25 21:52:43 MST 1999 aquarius% maple |\^/| Maple V Release 5.1 (WMI Campus Wide License) ._|\| |/|_. Copyright (c) 1981-1998 by Waterloo Maple Inc. All rights \ MAPLE / reserved. Maple and Maple V are registered trademarks of <____ ____> Waterloo Maple Inc. | Type ? for help. # ----------[ M a p l e ]---------- #interface(echo = 3); # ---------- Initialization ---------- > readlib(forget): > readlib(showtime): > _Envadditionally:= true: > on; # ---------- Definite Integrals ---------- # The following two functions have a pole at a. The first integral has a # principal value of zero; the second is divergent O1 := int(1/(x - a), x = (a - 1)..(a + 1)); bytes used=1000084, alloc=851812, time=0.62 a + 1 / | 1 | ----- dx | x - a / a - 1 time = 0.78, bytes = 1098650 O2 := int(1/(x - a), x = (a - 1)..(a + 1), CauchyPrincipalValue); 0 time = 0.08, bytes = 72206 O3 := int(1/(x - a)^2, x = (a - 1)..(a + 1)); infinity time = 0.05, bytes = 37854 # Different branches of the square root need to be chosen in the intervals # [0, 1] and [1, 2]. The correct results are 4/3, [4 - sqrt(8)]/3, # [8 - sqrt(8)]/3, respectively O4 := int(sqrt(x + 1/x - 2), x = 0..1); bytes used=2000304, alloc=1507052, time=1.70 4/3 time = 0.71, bytes = 607358 O5 := int(sqrt(x + 1/x - 2), x = 1..2); 1/2 - 2/3 2 + 4/3 time = 0.13, bytes = 129926 O6 := int(sqrt(x + 1/x - 2), x = 0..2); 1/2 - 2/3 2 + 8/3 time = 0.06, bytes = 57942 # => sqrt(2) [a modification of a problem due to W. Kahan] O7 := int(sqrt(2 - 2*cos(2*x))/2, x = -3*Pi/4..-Pi/4); bytes used=3000608, alloc=1834672, time=2.99 1/2 1/2 1/2 2 4 time = 1.87, bytes = 1743726 O8 := simplify(%); 1/2 2 time = 0.02, bytes = 21094 # Contour integrals => pi/a e^(-a) for a > 0. See Norman Levinson and # Raymond M. Redheffer, _Complex Variables_, Holden-Day, Inc., 1970, p. 198. O9 := assume(a > 0): time = 0.02, bytes = 15122 bytes used=4000856, alloc=2096768, time=4.08 O10 := Int(cos(x)/(x^2 + a^2), x = -infinity..infinity); infinity / | cos(x) | -------- dx | 2 2 / x + a~ -infinity time = 0.26, bytes = 4154 O11 := value(%); Pi (-sinh(a~) + cosh(a~)) ------------------------- a~ time = 0.92, bytes = 976990 O12 := convert(%, exp); bytes used=5001392, alloc=2424388, time=5.33 Pi ---------- exp(a~) a~ time = 0.33, bytes = 11054 # Integrand with a branch point => pi/sin(pi a) for 0 < a < 1 # [Levinson and Redheffer, p. 212] O13 := assume(a < 1): time = 0.03, bytes = 17658 O14 := Int(t^(a - 1)/(1 + t), t = 0..infinity); infinity / (a~ - 1) | t | --------- dt | 1 + t / 0 time = 0.01, bytes = 4046 O15 := value(%); bytes used=6004372, alloc=2686484, time=6.87 infinity / (a~ - 1) | t | --------- dt | 1 + t / 0 time = 1.89, bytes = 1353582 O16 := a:= 'a': time = 0.00, bytes = 3590 # Integrand with a residue at infinity => -2 pi [sin(pi/5) + sin(2 pi/5)] # (principal value) [Levinson and Redheffer, p. 234] O17 := int(5*x^3/(1 + x + x^2 + x^3 + x^4), x = -infinity..infinity, O17 := CauchyPrincipalValue); bytes used=7005344, alloc=3014104, time=8.27 bytes used=8005708, alloc=3210676, time=9.80 bytes used=9005952, alloc=3341724, time=11.41 infinity / 3 | x | 5 -------------------- dx | 2 3 4 / 1 + x + x + x + x -infinity time = 4.17, bytes = 2643838 # integrate(1/[1 + x + x^2 + ... + x^(2 n)], x = -infinity..infinity) # = 2 pi/(2 n + 1) [1 + cos(pi/[2 n + 1])] csc(2 pi/[2 n + 1]) # [Levinson and Redheffer, p. 255] => 2 pi/5 [1 + cos(pi/5)] csc(2 pi/5) O18 := int(1/(1 + x + x^2 + x^4), x = -infinity..infinity); bytes used=10006528, alloc=3341724, time=12.86 0 time = 2.02, bytes = 1621310 # Integrand with a residue at infinity and a branch cut => pi [sqrt(2) - 1] # [Levinson and Redheffer, p. 234] O19 := int(sqrt(1 - x^2)/(1 + x^2), x = -1..1); bytes used=11006696, alloc=3341724, time=14.30 bytes used=12006904, alloc=3472772, time=16.16 bytes used=13007076, alloc=3538296, time=18.01 bytes used=14007412, alloc=3538296, time=19.84 1/2 2 Pi - Pi time = 7.55, bytes = 4189470 O20 := factor(%); 1/2 Pi (2 - 1) time = 0.02, bytes = 11746 # This is a common integral in many physics calculations # => q/p sqrt(pi/p) e^(q^2/p) (Re p > 0) [Gradshteyn and Ryzhik 3.462(6)] O21 := assume(Re(p) > 0): time = 0.01, bytes = 19226 O22 := int(x*exp(-p*x^2 + 2*q*x), x = -infinity..infinity); bytes used=15007616, alloc=3538296, time=21.67 Definite integration: Can't determine if the integral is convergent. Need to know the sign of --> -2*q Will now try indefinite integration and then take limits. bytes used=16007788, alloc=3538296, time=23.43 infinity / | 2 | x exp(-p~ x + 2 q x) dx | / -infinity time = 3.40, bytes = 1912150 O23 := p:= 'p': time = 0.00, bytes = 3418 O24 := assume(p > 0): time = 0.01, bytes = 6230 O25 := int(x*exp(-p*x^2 + 2*q*x), x = -infinity..infinity); 2 q 1/2 exp(----) q Pi p~ ----------------- 3/2 p~ time = 0.07, bytes = 56578 O26 := p:= 'p': time = 0.01, bytes = 3242 # => 2 Euler's_constant [Gradshteyn and Ryzhik 8.367(5-6)] O27 := int(1/log(t) + 1/(1 - t) - log(log(1/t)), t = 0..1); bytes used=17007984, alloc=3603820, time=25.17 bytes used=18011652, alloc=3603820, time=27.07 bytes used=19011844, alloc=3734868, time=28.98 bytes used=20012240, alloc=3734868, time=30.72 bytes used=21013348, alloc=3800392, time=32.57 1 / | 1 1 | ----- + ----- - ln(ln(1/t)) dt | ln(t) 1 - t / 0 time = 8.92, bytes = 4791842 # This integral comes from atomic collision theory => 0 [John Prentice] O28 := int(sin(t)/t*exp(2*I*t), t = -infinity..infinity); bytes used=22013576, alloc=3865916, time=34.36 infinity / | sin(t) exp(2 I t) | ----------------- dt | t / -infinity time = 1.58, bytes = 977294 # => 1/12 [Gradshteyn and Ryzhik 6.443(3)] O29 := int(log(gamma(x))*cos(6*Pi*x), x = 0..1); 1 / | | ln(gamma(x)) cos(6 Pi x) dx | / 0 time = 0.20, bytes = 171610 # => 36/35 [Gradshteyn and Ryzhik 7.222(2)] O30 := int((1 + x)^3*orthopoly[P](1, x)*orthopoly[P](2, x), x = -1..1); 36 -- 35 time = 0.03, bytes = 25486 # => 1/sqrt(a^2 + b^2) (a > 0 and b real) # [Gradshteyn and Ryzhik 6.611(1)] O31 := assume(a > 0): time = 0.02, bytes = 12702 O32 := assume(b, real): time = 0.01, bytes = 5418 O33 := int(exp(-a*x)*BesselJ(0, b*x), x = 0..infinity); bytes used=23013828, alloc=3865916, time=36.01 Definite integration: Can't determine if the integral is convergent. Need to know the sign of --> b Will now try indefinite integration and then take limits. bytes used=24014004, alloc=3931440, time=38.05 bytes used=25014544, alloc=3931440, time=39.88 infinity / | | exp(-a~ x) BesselJ(0, b~ x) dx | / 0 time = 4.57, bytes = 2143022 O34 := b:= 'b': time = 0.00, bytes = 3678 O35 := q:= int(exp(-a*x)*BesselJ(0, b*x), x = 0..infinity); Definite integration: Can't determine if the integral is convergent. Need to know the sign of --> b Will now try indefinite integration and then take limits. 2 a~ hypergeom([1/2], [], - ---) 2 b --------------------------- b time = 0.85, bytes = 636818 O36 := assume(b, real): time = 0.01, bytes = 4994 O37 := simplify(q); signum(b~) -------------- 2 2 1/2 (b~ + a~ ) time = 0.20, bytes = 164014 O38 := q:= 'q': time = 0.00, bytes = 3318 O39 := a:= 'a': time = 0.00, bytes = 3406 O40 := b:= 'b': time = 0.01, bytes = 3270 # Integrand contains a special function => 4/(3 pi) [Tom Hagstrom] O41 := int((BesselJ(1, x)/x)^2, x = 0..infinity); 1 4/3 ---- Pi bytes used=26017720, alloc=3931440, time=41.67 time = 0.14, bytes = 125594 # => (cos 7 - 1)/7 [Gradshteyn and Ryzhik 6.782(3)] O42 := int(Ci(x)*BesselJ(0, 2*sqrt(7*x)), x = 0..infinity); infinity / | 1/2 1/2 | Ci(x) BesselJ(0, 2 7 x ) dx | / 0 time = 0.82, bytes = 650774 # This integral comes from doing a two loop Feynman diagram for a QCD problem # => - [17/3 + pi^2]/36 + log 2/9 [35/3 - pi^2/2 - 4 log 2 + log(2)^2] # + zeta(3)/4 = 0.210883... [Rolf Mertig] */ O43 := int(x^2*polylog(3, 1/(x + 1)), x = 0..1); bytes used=27017960, alloc=3931440, time=43.44 1 / | 2 1 | x polylog(3, -----) dx | x + 1 / 0 time = 1.45, bytes = 794602 O44 := evalf(Int(x^2*polylog(3, 1/(x + 1)), x = 0..1)); bytes used=28018272, alloc=3931440, time=45.04 bytes used=29018920, alloc=4128012, time=46.52 bytes used=30019428, alloc=4324584, time=48.22 bytes used=31019608, alloc=4324584, time=49.94 bytes used=32019796, alloc=4324584, time=51.53 bytes used=33020212, alloc=4324584, time=53.08 bytes used=34020388, alloc=4324584, time=54.64 bytes used=35020696, alloc=4324584, time=56.22 bytes used=36020908, alloc=4324584, time=57.79 bytes used=37021308, alloc=4324584, time=59.37 bytes used=38021660, alloc=4324584, time=60.95 .2108828596 time = 17.30, bytes = 10841322 O45 := evalf(- (17/3 + Pi^2)/36 + log(2)/9*(35/3 - Pi^2/2 - 4*log(2) + log(2)^2) O45 := + Zeta(3)/4); .2108828597 time = 0.01, bytes = 12394 # Integrate a piecewise defined step function s(t) multiplied by cos t, where # s(t) = 0 (t < 1); 1 (1 <= t <= 2); 0 (t > 2) # => 0 (u < 1); sin u - sin 1 (1 <= u <= 2); sin 2 - sin 1 (u > 2) O46 := s:= proc(t) if 1 <= t and t <= 2 then 1 else 0 fi end: time = 0.03, bytes = 14078 O47 := int(s(t)*cos(t), t = 0..u); Error, cannot evaluate boolean time = 0.01, bytes = 3634 O48 := s:= t -> piecewise(1 <= t and t <= 2, 1, 0): time = 0.00, bytes = 3626 O49 := int(s(t)*cos(t), t = 0..u); -sin(u) Heaviside(1 - u) + sin(u) Heaviside(-u + 2) + sin(2) + sin(1) Heaviside(1 - u) - sin(1) - sin(2) Heaviside(-u + 2) time = 0.59, bytes = 501010 O50 := convert(%, piecewise); { 0 u <= 1 { { sin(u) - sin(1) u <= 2 { { -sin(1) + sin(2) 2 < u time = 0.12, bytes = 100426 O51 := s:= 's': time = 0.00, bytes = 3662 # Integrating first with respect to y and then x is much easier than # integrating first with respect to x and then y # => (|b| - |a|) pi [W. Kahan] O52 := assume(a, real, b, real): time = 0.02, bytes = 16534 O53 := int(int(x/(x^2 + y^2), y = -infinity..infinity), x = a..b); bytes used=39022756, alloc=4324584, time=62.67 b~ Pi signum(b~) - a~ Pi signum(a~) time = 1.50, bytes = 652742 O54 := factor(simplify(%)); Pi (| b~ | - | a~ |) time = 0.06, bytes = 43102 O55 := int(int(x/(x^2 + y^2), x = a..b), y = -infinity..infinity); bytes used=40023092, alloc=4324584, time=64.63 bytes used=41023400, alloc=4324584, time=66.57 bytes used=42023604, alloc=4455632, time=68.43 b~ Pi signum(b~) - a~ Pi signum(a~) + 1/2 %2 ln(-%2) - 1/2 %2 ln(%2) - 1/2 %1 ln(-%1) + 1/2 %1 ln(%1) 2 2 %1 := RootOf(a~ + _Z ) 2 2 %2 := RootOf(b~ + _Z ) time = 5.99, bytes = 3240298 O56 := simplify(%); bytes used=43024312, alloc=4455632, time=70.31 Pi | b~ | - Pi | a~ | + 1/2 %2 ln(-%2) - 1/2 %2 ln(%2) - 1/2 %1 ln(-%1) + 1/2 %1 ln(%1) 2 2 %1 := RootOf(a~ + _Z ) 2 2 %2 := RootOf(b~ + _Z ) time = 0.74, bytes = 124154 O57 := a:= 'a': time = 0.01, bytes = 3422 O58 := b:= 'b': time = 0.01, bytes = 3306 # => [log(sqrt(2) + 1) + sqrt(2)]/3 [Caviness et all, section 2.10.1] O59 := int(int(sqrt(x^2 + y^2), x = 0..1), y = 0..1); 1 / | 1 | csgn(y) y hypergeom([1/2, -1/2], [3/2], - ----) dy | 2 / y 0 time = 0.67, bytes = 497790 # => (pi a)/2 [Gradshteyn and Ryzhik 4.621(1)] O60 := int(int(sin(a)*sin(y)/sqrt(1 - sin(a)^2*sin(x)^2*sin(y)^2), O60 := x = 0..Pi/2), O60 := y = 0..Pi/2); bytes used=44024484, alloc=4455632, time=72.16 bytes used=45024684, alloc=4521156, time=73.95 bytes used=46024948, alloc=4586680, time=75.91 bytes used=47025124, alloc=4586680, time=77.85 2 1/2 sin(a) Pi hypergeom([1/2, 1/2], [3/2], sin(a) ) time = 7.18, bytes = 3713470 O61 := convert(%, StandardFunctions); 1/2 Pi arcsin(sin(a)) time = 0.09, bytes = 102482 # => 46/15 [Paul Zimmermann] O62 := int(int(abs(y - x^2), y = 0..2), x = -1..1); bytes used=48025320, alloc=4586680, time=79.76 46 -- 15 time = 1.61, bytes = 738574 # Multiple integrals: volume of a tetrahedron => a b c / 6 O63 := Int(Int(Int(1, z = 0..c*(1 - x/a - y/b)), O63 := y = 0..b*(1 - x/a)), O63 := x = 0..a); a b (1 - x/a) c (1 - x/a - y/b) / / / | | | | | | 1 dz dy dx | | | / / / 0 0 0 time = 0.03, bytes = 13722 O64 := value(map(value, %)); 1/6 c b a time = 0.20, bytes = 157730 # ---------- Quit ---------- O65 := quit bytes used=48370524, alloc=4586680, time=80.21 real 82.67 user 80.27 sys 2.30