| blob | 697c83303d5c7b86f2b1ebcc175529792a517419 |
1 /* Original version of this file copyright 1999 by Michael Wester,
2 * and retrieved from http://www.math.unm.edu/~wester/demos/DefiniteIntegrals/problems.macsyma
3 * circa 2006-10-23.
4 *
5 * Released under the terms of the GNU General Public License, version 2,
6 * per message dated 2007-06-03 from Michael Wester to Robert Dodier
7 * (contained in the file wester-gpl-permission-message.txt).
8 *
9 * See: "A Critique of the Mathematical Abilities of CA Systems"
10 * by Michael Wester, pp 25--60 in
11 * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
12 * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
13 */
14 /* ----------[ M a c s y m a ]---------- */
15 /* ---------- Initialization ---------- */
16 showtime: all$
17 prederror: false$
18 /* ---------- Definite Integrals ---------- */
19 /* The following two functions have a pole at a. The first integral has a
20 principal value of zero; the second is divergent */
21 integrate(1/(x - a), x, a - 1, a + 1);
22 errcatch(integrate(1/(x - a)^2, x, a - 1, a + 1));
23 /* Different branches of the square root need to be chosen in the intervals
24 [0, 1] and [1, 2]. The correct results are 4/3, [4 - sqrt(8)]/3,
25 [8 - sqrt(8)]/3, respectively */
26 integrate(sqrt(x + 1/x - 2), x, 0, 1);
27 integrate(sqrt(x + 1/x - 2), x, 1, 2);
28 integrate(sqrt(x + 1/x - 2), x, 0, 2);
29 /* => sqrt(2) [a modification of a problem due to W. Kahan] */
30 integrate(sqrt(2 - 2*cos(2*x))/2, x, -3*%pi/4, -%pi/4);
31 /* Contour integrals => pi/a e^(-a) for a > 0. See Norman Levinson and
32 Raymond M. Redheffer, _Complex Variables_, Holden-Day, Inc., 1970, p. 198.
33 */
34 assume(a > 0)$
35 'integrate(cos(x)/(x^2 + a^2), x, -inf, inf);
36 ev(%, integrate);
37 /* Integrand with a branch point => pi/sin(pi a) for 0 < a < 1
38 [Levinson and Redheffer, p. 212] */
39 assume(a < 1)$
40 declare(a, noninteger)$
41 'integrate(t^(a - 1)/(1 + t), t, 0, inf);
42 ev(%, integrate);
43 remove(a, noninteger)$
44 forget(a > 0, a < 1)$
45 /* Integrand with a residue at infinity => -2 pi [sin(pi/5) + sin(2 pi/5)]
46 (principal value) [Levinson and Redheffer, p. 234] */
47 errcatch(integrate(5*x^3/(1 + x + x^2 + x^3 + x^4), x, -inf, inf));
48 integrate(5*x^3/(1 + x + x^2 + x^3 + x^4), x, -inf, inf), intanalysis = false;
49 /* integrate(1/[1 + x + x^2 + ... + x^(2 n)], x = -infinity..infinity)
50 = 2 pi/(2 n + 1) [1 + cos(pi/[2 n + 1])] csc(2 pi/[2 n + 1])
51 [Levinson and Redheffer, p. 255] => 2 pi/5 [1 + cos(pi/5)] csc(2 pi/5) */
52 q: (2*atan(43*sqrt(3)/(9*sqrt(257))) + %pi)/6;
53 assume(equal(cos(q)*sin(2*q) - sin(q)*cos(2*q) - sin(q), 0))$
54 integrate(1/(1 + x + x^2 + x^4), x, -inf, inf);
55 forget(equal(cos(q)*sin(2*q) - sin(q)*cos(2*q) - sin(q), 0))$
56 remvalue(q)$
57 /* Integrand with a residue at infinity and a branch cut => pi [sqrt(2) - 1]
58 [Levinson and Redheffer, p. 234] */
59 integrate(sqrt(1 - x^2)/(1 + x^2), x, -1, 1);
60 factor(%);
61 /* This is a common integral in many physics calculations
62 => q/p sqrt(pi/p) e^(q^2/p) (Re p > 0) [Gradshteyn and Ryzhik 3.462(6)]
63 */
64 assume(p > 0)$
65 integrate(x*exp(-p*x^2 + 2*q*x), x, -inf, inf);
66 forget(p > 0)$
67 /* => 2 Euler's_constant [Gradshteyn and Ryzhik 8.367(5-6)] */
68 integrate(1/log(t) + 1/(1 - t) - log(log(1/t)), t, 0, 1);
69 /* This integral comes from atomic collision theory => 0 [John Prentice] */
70 integrate(sin(t)/t*exp(2*%i*t), t, -inf, inf);
71 /* => 1/12 [Gradshteyn and Ryzhik 6.443(3)] */
72 integrate(log(gamma(x))*cos(6*%pi*x), x, 0, 1);
73 /* => 36/35 [Gradshteyn and Ryzhik 7.222(2)] */
74 integrate((1 + x)^3*legendre_p(1, x)*legendre_p(2, x), x, -1, 1);
75 /* => 1/sqrt(a^2 + b^2) (a > 0 and b real)
76 [Gradshteyn and Ryzhik 6.611(1)] */
77 integrate(exp(-a*x)*bessel_j[0](b*x), x, 0, inf);
78 /* Integrand contains a special function => 4/(3 pi) [Tom Hagstrom] */
79 integrate((bessel_j[1](x)/x)^2, x, 0, inf);
80 /* => (cos 7 - 1)/7 [Gradshteyn and Ryzhik 6.782(3)] */
81 integrate(cos_int(x)*bessel_j[0](2*sqrt(7*x)), x, 0, inf);
82 /* This integral comes from doing a two loop Feynman diagram for a QCD problem
83 => - [17/3 + pi^2]/36 + log 2/9 [35/3 - pi^2/2 - 4 log 2 + log(2)^2]
84 + zeta(3)/4 = 0.210883... [Rolf Mertig] */
85 integrate(x^2*li[3](1/(x + 1)), x, 0, 1);
86 romberg(x^2*li[3](1/(x + 1)), x, 0, 1);
87 sfloat(- (17/3 + %pi^2)/36 + log(2)/9*(35/3 - %pi^2/2 - 4*log(2) + log(2)^2)
88 + zeta(3)/4);
89 /* Integrate a piecewise defined step function s(t) multiplied by cos t, where
90 s(t) = 0 (t < 1); 1 (1 <= t <= 2); 0 (t > 2)
91 => 0 (u < 1); sin u - sin 1 (1 <= u <= 2); sin 2 - sin 1 (u > 2)
92 */
93 s(t):= if 1 <= t and t <= 2 then 1 else 0$
94 integrate(s(t)*cos(t), t, 0, u);
95 s(t):= unit_step(t - 1) - unit_step(t - 2)$
96 integrate(s(t)*cos(t), t, 0, u);
97 ratsimp(%);
98 remfunction(s)$
99 /* Integrating first with respect to y and then x is much easier than
100 integrating first with respect to x and then y
101 => (|b| - |a|) pi [W. Kahan] */
102 assume(a > 0, b > 0)$
103 integrate(integrate(x/(x^2 + y^2), y, -inf, inf), x, a, b);
104 integrate(integrate(x/(x^2 + y^2), x, a, b), y, -inf, inf);
105 (forget(a > 0, b > 0), assume(a < 0, b > 0))$
106 integrate(integrate(x/(x^2 + y^2), y, -inf, inf), x, a, b);
107 integrate(integrate(x/(x^2 + y^2), x, a, b), y, -inf, inf);
108 (forget(a < 0, b > 0), assume(a < 0, b < 0))$
109 integrate(integrate(x/(x^2 + y^2), y, -inf, inf), x, a, b);
110 integrate(integrate(x/(x^2 + y^2), x, a, b), y, -inf, inf);
111 forget(a < 0, b < 0)$
112 /* => [log(sqrt(2) + 1) + sqrt(2)]/3 [Caviness et all, section 2.10.1] */
113 assume(not(equal(y, 0)))$
114 integrate(integrate(sqrt(x^2 + y^2), x, 0, 1), y, 0, 1);
115 ratsimp(logarc(%));
116 factor(log(sqrtdenest(sqrtdenest(sqrtdenest(radcan(exp(logcontract(%))))))));
117 tldefint(integrate(sqrt(x^2 + y^2), x, 0, 1), y, 0, 1);
118 forget(not(equal(y, 0)))$
119 /* => (pi a)/2 [Gradshteyn and Ryzhik 4.621(1)] */
120 assume((sin(a)*sin(y) - 1)*(sin(a)*sin(y) + 1) < 0)$
121 integrate(integrate(sin(a)*sin(y)/sqrt(1 - sin(a)^2*sin(x)^2*sin(y)^2),
122 x, 0, %pi/2),
123 y, 0, %pi/2);
124 forget((sin(a)*sin(y) - 1)*(sin(a)*sin(y) + 1) < 0)$
125 /* => 46/15 [Paul Zimmermann] */
126 assume(not(equal(x, 0)), x^2 < 2)$
127 integrate(integrate(abs(y - x^2), y, 0, 2), x, -1, 1);
128 forget(not(equal(x, 0)), x^2 < 2)$
129 /* Multiple integrals: volume of a tetrahedron => a b c / 6 */
130 'integrate('integrate('integrate(1, z, 0, c*(1 - x/a - y/b)),
131 y, 0, b*(1 - x/a)),
132 x, 0, a);
133 ev(%, integrate);