Tue Mar 10 22:18:15 MST 1998 aquarius% mupad *----* MuPAD 1.4.0 -- Multi Processing Algebra Data Tool /| /| *----* | Copyright (c) 1997 - 98 by SciFace Software GmbH | *--|-* All rights reserved. |/ |/ *----* Licensed to: Michael Wester >> # ----------[ M u P A D ]---------- # >> # ---------- Initialization ---------- # >> TEXTWIDTH:= 80: >> read("../../Time.mupad"): >> # ---------- Definite Integrals ---------- # >> # The following two functions have a pole at a. The first integral has a &> principal value of zero; the second is divergent # >> int(1/(x - a), x = (a - 1)..(a + 1)); undefined Time: 34170 msec Type: undefined >> int(1/(x - a), x = (a - 1)..(a + 1), PrincipalValue); 0 Time: 920 msec Type: DOM_INT >> int(1/(x - a)^2, x = (a - 1)..(a + 1)); infinity Time: 1070 msec Type: stdlib::Infinity >> # 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 # >> int(sqrt(x + 1/x - 2), x = 0..1); 4/3 Time: 17880 msec Type: DOM_RAT >> int(sqrt(x + 1/x - 2), x = 1..2); 1/2 2 2 4/3 - ------ 3 Time: 1590 msec Type: "_plus" >> int(sqrt(x + 1/x - 2), x = 0..2); 1/2 2 2 8/3 - ------ 3 Time: 2470 msec Type: "_plus" >> # => sqrt(2) [a modification of a problem due to W. Kahan] # >> int(sqrt(2 - 2*cos(2*x))/2, x = -3*PI/4..-PI/4); 1/2 2 Time: 2840 msec Type: "_power" >> # Contour integrals => pi/a e^(-a) for a > 0. See Norman Levinson and &> Raymond M. Redheffer, _Complex Variables_, Holden-Day, Inc., 1970, p. 198. # >> assume(a > 0): >> int(cos(x)/(x^2 + a^2), x = -infinity..infinity); PI exp(-a) ---------- a Time: 5560 msec Type: "_mult" >> # Integrand with a branch point => pi/sin(pi a) for 0 < a < 1 &> [Levinson and Redheffer, p. 212] # >> assume(a < 1, _and): >> int(t^(a - 1)/(1 + t), t = 0..infinity); gamma(a) gamma(1 - a) Time: 2920 msec Type: "_mult" >> unassume(a): >> # Integrand with a residue at infinity => -2 pi [sin(pi/5) + sin(2 pi/5)] &> (principal value) [Levinson and Redheffer, p. 234] # >> int(5*x^3/(1 + x + x^2 + x^3 + x^4), x = -infinity..infinity, PrincipalValue); / 3 \ | 5 x | int| --------------------, x = -infinity..infinity | | 2 3 4 | \ x + x + x + x + 1 / Time: 15130 msec Type: "int" >> # 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) # >> int(1/(1 + x + x^2 + x^4), x = -infinity..infinity); / 1 \ int| ---------------, x = -infinity..infinity | | 2 4 | \ x + x + x + 1 / Time: 3340 msec Type: "int" >> # Integrand with a residue at infinity and a branch cut => pi [sqrt(2) - 1] &> [Levinson and Redheffer, p. 234] # >> int(sqrt(1 - x^2)/(1 + x^2), x = -1..1); / 2 1/2 \ | (1 - x ) | int| -----------, x = -1..1 | | 2 | \ x + 1 / Time: 4080 msec Type: "int" >> # 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)] &> # >> assume(Re(p) > 0): >> int(x*exp(-p*x^2 + 2*q*x), x = -infinity..infinity); / / 2 \ / 3/2 1/2 \ | 1/2 | q | | p x - p q | 3/2 2 | p q PI exp| -- | erf| --------------- | - p exp(2 q x - p x ) | \ p / \ p / limit| -------------------------------------------------------------------, | 5/2 \ 2 p \ / | | | | | | x = infinity | - limit| | | / \ / 2 \ / 3/2 1/2 \ 1/2 | q | | p x - p q | 3/2 2 p q PI exp| -- | erf| --------------- | - p exp(2 q x - p x ) \ p / \ p / -------------------------------------------------------------------, 5/2 2 p \ | | | x = -infinity | | / Time: 11020 msec Type: "_plus" >> unassume(Re(p) > 0): >> # => 2 Euler's_constant [Gradshteyn and Ryzhik 8.367(5-6)] # >> int(1/ln(t) + 1/(1 - t) - ln(ln(1/t)), t = 0..1); / 1 / / 1 \ \ 1 \ int| ----- - ln| ln| - | | + -----, t = 0..1 | \ ln(t) \ \ t / / 1 - t / Time: 23790 msec Type: "int" >> # This integral comes from atomic collision theory => 0 [John Prentice] # >> int(sin(t)/t*exp(2*I*t), t = -infinity..infinity); / sin(t) exp(2 I t) \ int| -----------------, t = -infinity..infinity | \ t / Time: 7380 msec Type: "int" >> # => 1/12 [Gradshteyn and Ryzhik 6.443(3)] # >> int(ln(gamma(x))*cos(6*PI*x), x = 0..1); int(ln(gamma(x)) cos(6 x PI), x = 0..1) Time: 4730 msec Type: "int" >> # => 36/35 [Gradshteyn and Ryzhik 7.222(2)] # >> int((1 + x)^3*orthpoly::legendre(1, x)*orthpoly::legendre(2, x), x = -1..1); Error: Illegal indeterminate [poly]; in procedure 'intlib::definite' >> # => 1/sqrt(a^2 + b^2) (a > 0 and b real) &> [Gradshteyn and Ryzhik 6.611(1)] # >> assume(a > 0): >> int(exp(-a*x)*besselJ(0, b*x), x = 0..infinity); 1 ------------ 2 2 1/2 (a + b ) Time: 7720 msec Type: "_power" >> unassume(a > 0): >> # Integrand contains a special function => 4/(3 pi) [Tom Hagstrom] # >> int((besselJ(1, x)/x)^2, x = 0..infinity); Error: Argument out of range [DIGITS]; in procedure 'funcattr(besselJ, "float")' >> # => (cos 7 - 1)/7 [Gradshteyn and Ryzhik 6.782(3)] # >> #int(Ci(x)*besselJ(0, 2*sqrt(7*x)), x = 0..infinity);# >> # 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] # >> int(x^2*polylog(3, 1/(x + 1)), x = 0..1); / 2 / 1 \ \ int| x polylog| 3, ----- |, x = 0..1 | \ \ x + 1 / / Time: 14270 msec Type: "int" >> float(hold(int(x^2*polylog(3, 1/(x + 1)), x = 0..1))); / 2 / 1 \ \ int| x polylog| 3, ----- |, x = 0..1 | \ \ x + 1 / / Time: 7030 msec Type: "int" >> float(- (17/3 + PI^2)/36 + ln(2)/9*(35/3 - PI^2/2 - 4*ln(2) + ln(2)^2) &> + zeta(3)/4); 0.2108828595 Time: 500 msec Type: DOM_FLOAT >> # 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) &> # >> s:= proc(t) begin if 1 <= t and t <= 2 then 1 else 0 end_if end_proc: >> int(s(t)*cos(t), t = 0..u); Error: Can't evaluate to boolean [_less]; in procedure 's' >> s:= func(heaviside(t - 1) - heaviside(t - 2), t): >> int(s(t)*cos(t), t = 0..u); int(cos(t) (heaviside(t - 1) - heaviside(t - 2)), t = 0..u) Time: 7610 msec Type: "int" >> s:= NIL: >> # 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] # >> assume(a > 0): >> assume(b > 0): >> int(int(x/(x^2 + y^2), y = -infinity..infinity), x = a..b); 2 2 b PI a PI ------- - ------- 2 1/2 2 1/2 (b ) (a ) Time: 6070 msec Type: "_plus" >> simplify(%); b PI - a PI Time: 760 msec Type: "_plus" >> int(int(x/(x^2 + y^2), x = a..b), y = -infinity..infinity); / / x \ \ int| int| -------, x = a..b |, y = -infinity..infinity | | | 2 2 | | \ \ x + y / / Time: 11510 msec Type: "int" >> unassume(a > 0): >> unassume(b > 0): >> assume(a < 0): >> assume(b > 0): >> int(int(x/(x^2 + y^2), y = -infinity..infinity), x = a..b); 2 2 b PI a PI ------- - ------- 2 1/2 2 1/2 (b ) (a ) Time: 6520 msec Type: "_plus" >> simplify(%); a PI + b PI Time: 780 msec Type: "_plus" >> int(int(x/(x^2 + y^2), x = a..b), y = -infinity..infinity); / / x \ \ int| int| -------, x = a..b |, y = -infinity..infinity | | | 2 2 | | \ \ x + y / / Time: 6960 msec Type: "int" >> unassume(a < 0): >> unassume(b > 0): >> assume(a < 0): >> assume(b < 0): >> int(int(x/(x^2 + y^2), y = -infinity..infinity), x = a..b); 2 2 b PI a PI ------- - ------- 2 1/2 2 1/2 (b ) (a ) Time: 5880 msec Type: "_plus" >> simplify(%); a PI - b PI Time: 780 msec Type: "_plus" >> int(int(x/(x^2 + y^2), x = a..b), y = -infinity..infinity); / / x \ \ int| int| -------, x = a..b |, y = -infinity..infinity | | | 2 2 | | \ \ x + y / / Time: 7280 msec Type: "int" >> unassume(a < 0): >> unassume(b < 0): >> # => [log(sqrt(2) + 1) + sqrt(2)]/3 [Caviness et all, section 2.10.1] # >> int(int(sqrt(x^2 + y^2), x = 0..1), y = 0..1); / 2 1/2 2 2 1/2 2 2 1/2 \ | (y + 1) y ln((y ) ) y ln((y + 1) + 1) | int| ----------- - -------------- + ----------------------, y = 0..1 | \ 2 2 2 / Time: 30510 msec Type: "int" >> # => (pi a)/2 [Gradshteyn and Ryzhik 4.621(1)] # >> int(int(sin(a)*sin(y)/sqrt(1 - sin(a)^2*sin(x)^2*sin(y)^2), &> x = 0..PI/2), &> y = 0..PI/2); / 2 2 1/2 PI \ int| sin(a) sin(y) ellipticK((sin(a) sin(y) ) ), y = 0..-- | \ 2 / Time: 20000 msec Type: "int" >> # => 46/15 [Paul Zimmermann] # >> int(int(abs(y - x^2), y = 0..2), x = -1..1); 2 int(int(abs(y - x ), y = 0..2), x = -1..1) Time: 940 msec Type: "int" >> # Multiple integrals: volume of a tetrahedron => a b c / 6 # >> int(int(int(1, z = 0..c*(1 - x/a - y/b)), &> y = 0..b*(1 - x/a)), &> x = 0..a); a b c ----- 6 Time: 550 msec Type: "_mult" >> # ---------- Quit ---------- # >> quit real 298.83 user 277.36 sys 2.11