# Lesson 15. Symbolic and Numerical Integration
# -----------------------------------
#
# Last time we saw that there can be more to a definite integral than F(b) - F(a).
> restart:\
J:= int(sqrt(tan(x)), x);
1/2 1/2 1/2 1/2
1/2 2 tan(x) 1/2 tan(x) + 2 tan(x) + 1
J := 1/2 2 arctan(--------------) - 1/2 2 ln(---------------------------)
1 - tan(x) 2 1/2
(1 + tan(x) )
> int(sqrt(tan(x)), x=0 .. Pi/3);
1/2 1/4
1/2 2 3 1/2 1/2 1/2 1/4 1/2
- 1/2 2 arctan(----------) - 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 )
1/2
- 1 + 3
1/2
+ 1/2 2 Pi
> evalf(");
.787779048
> subs(x=Pi/3, J) - subs(x = 0, J);
1/2 1/2
1/2 2 tan(1/3 Pi)
1/2 2 arctan(-------------------)
1 - tan(1/3 Pi)
1/2 1/2
1/2 tan(1/3 Pi) + 2 tan(1/3 Pi) + 1
- 1/2 2 ln(-------------------------------------)
2 1/2
(1 + tan(1/3 Pi) )
1/2 1/2
1/2 2 tan(0)
- 1/2 2 arctan(--------------)
1 - tan(0)
1/2 1/2
1/2 tan(0) + 2 tan(0) + 1
+ 1/2 2 ln(---------------------------)
2 1/2
(1 + tan(0) )
> ";
1/2 1/4
1/2 2 3 1/2 1/2 1/2 1/4 1/2
1/2 2 arctan(---------) - 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 )
1/2
1 - 3
> evalf(");
-1.433662421
# The discontinuities of the antiderivative are found with the "discont" command. This must first be read
# in to Maple (although that's not necessary if you evaluate a definite integral first).
> readlib(discont):
> discont(J,x);
{1/4 Pi + Pi _Z~, Pi _Z7~ + 1/2 Pi}
# A word about "arctan". The discontinuity at Pi/4 in this case could be removed with a better choice of
# formula for the antiderivative, based on arccot instead of arctan.
> J;
1/2 1/2 1/2 1/2
1/2 2 tan(x) 1/2 tan(x) + 2 tan(x) + 1
1/2 2 arctan(--------------) - 1/2 2 ln(---------------------------)
1 - tan(x) 2 1/2
(1 + tan(x) )
> Jp:= sqrt(2)/2*arccot((1-tan(x))/(sqrt(2*tan(x)))) + op(2,J);
1/2
1/2 (1 - tan(x)) 2
Jp := 1/2 2 arccot(1/2 -----------------)
1/2
tan(x)
1/2 1/2
1/2 tan(x) + 2 tan(x) + 1
- 1/2 2 ln(---------------------------)
2 1/2
(1 + tan(x) )
# This has no discontinuity when tan(x) = 1. It does have one when tan(x)=0, but that's at an endpoint of
# our interval.
> subs(x=Pi/3, Jp) - limit(Jp, x=0,right);
1/2
1/2 (1 - tan(1/3 Pi)) 2
1/2 2 arccot(1/2 ----------------------)
1/2
tan(1/3 Pi)
1/2 1/2
1/2 tan(1/3 Pi) + 2 tan(1/3 Pi) + 1
- 1/2 2 ln(-------------------------------------)
2 1/2
(1 + tan(1/3 Pi) )
> ";
1/2 1/2
1/2 (1 - 3 ) 2
1/2 2 arccot(1/2 ---------------)
1/4
3
1/2 1/2 1/2 1/4 1/2
- 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 )
> evalf(");
.7877790480
# There is an option to "int" that omits the search for discontinuities in the interval of integration. This
# may sometimes speed things up, but use it at your own risk!
> int(sqrt(tan(x)), x=0 .. Pi/3, continuous);
1/2 1/4
1/2 2 3 1/2 1/2 1/2 1/4 1/2
- 1/2 2 arctan(----------) - 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 )
1/2
- 1 + 3
# What if we had an example where Maple couldn't identify the discontinuities?
> g:= normal(diff(1/(2*x-tan(x)),x));\
2
- 1 + tan(x)
g := -----------------
2
(- 2 x + tan(x))
> J:= int(g,x);
1
J := - --------------
- 2 x + tan(x)
> discont(J,x);
{Pi _Z11~ + 1/2 Pi, 0}
# It didn't see the discontinuity at the point between 0 and Pi/2 where
# tan(x) = 2*x.
> plot({2*x, tan(x)}, x=0..Pi/2, 0..Pi);
> int(g, x=Pi/4 .. 3*Pi/8) = subs(x=3*Pi/8, J) - subs(x=Pi/4, J);
4 2 1 1
----------------- - ------ = - ---------------------- + ----------------------
1/2 Pi - 2 - 3/4 Pi + tan(3/8 Pi) - 1/2 Pi + tan(1/4 Pi)
3 Pi - 4 2 - 4
> evalf(");
-18.98765027 = -18.98765042
# Why is this answer absurd? What should the answer be?
# Apart from such mishaps, Maple can do improper integrals. It checks for convergence or divergence.
# Here are some you probably met in Math 101.
> int(1/x^2, x= 1 .. infinity);
1
> int(1/x, x= -1 .. 1);
1
/
|
| 1/x dx
|
/
-1
> int(exp(-a*x^2), x=-infinity .. infinity);
infinity
/
| 2
| exp(- a x ) dx
|
/
- infinity
# Why didn't that one work? Because nobody told Maple that a > 0. If a < 0 the integral doesn't converge.
# Don't expect Maple to read your mind! And remember that variables are allowed to be complex unless
# otherwise specified, so replacing a with a^2 wouldn't have worked either. When you want to make an
# assumption on a variable, you can use "assume".
> assume(a > 0);
> int(exp(-a*x^2), x=-infinity .. infinity);
1/2
Pi
-----
1/2
a~
# Maple reminds you that there is an assumption on the variable "a" by printing it as "a~".
# Here are some that you may not have met. Maple evaluates
# them correctly.
> int(sin(x)/x, x=-infinity .. infinity);
Pi
> int((ln(x))^2/(1+x^2), x= 0 .. infinity);
3
1/8 Pi
# In this case the answer didn't come from a formula for the antiderivative, since in fact Maple can't find
# an antiderivative.
# In Math 300 and 301 you learn how to evaluate integrals such as this using residues.
> int((ln(x))^2/(1+x^2), x);
/ 2
| ln(x)
| ------ dx
| 2
/ 1 + x
# On the other hand...
> int(sin(x)/(1+x^2), x=-infinity .. infinity);
I sinh(1) Pi
# Oops again. This is another case where Maple's formula for the antiderivative has a discontinuity, this
# time involving complex numbers. It's actually something you'll meet in Math 301, called a branch cut.
#
# Conclusion: don't put too much trust in "int" for definite integrals. It's good to be able to check its
# answers. This can often be done using numerical integration, which we'll look at next.
# We can use "evalf" on an integral, even if Maple can't do the integration symbolically.
> int(exp(sin(x)), x=0..1);
1
/
|
| exp(sin(x)) dx
|
/
0
> evalf(");
1.631869608
# This is using numerical techniques (which we will discuss). In fact we can avoid the attempt at symbolic
# integration:
> evalf(Int(exp(sin(x)), x=0 .. 1));
1.631869608
# evalf(Int(...)) can be used for improper integrals as well.
> evalf(Int(1/x, x=-1 .. 1));
Error, (in evalf/int) unable to handle singularity
> evalf(Int(1/x^2, x= 1 .. infinity));
1.
> evalf(Int(ln(x)^2/(1+x^2), x=0 .. infinity));
3.875784585
# Here is the answer we obtained for that one with "int".
> evalf(Pi^3/8);
3.875784586
# One thing evalf(Int(...)) can't handle is a symbolic variable other than the variable of integration.
> evalf(Int(exp(-a*x^2), x = -infinity .. infinity));
Error, (in evalf/int) function does not evaluate to numeric
> evalf(Int(exp(-x^2), x = -infinity .. infinity)) = evalf(sqrt(Pi));
1.772453851 = 1.772453851
# Warning!! If you try evalf(Int(...)) on an improper integral that diverges, but Maple doesn't detect that it
# diverges, the result could be a VERY long attempt by Maple to approximate the integral numerically.
# Since the desired accuracy can never be attained, Maple might keep trying forever. In some cases,
# especially on your home PC, the usual methods of interrupting Maple's computation might not work, and
# you could be forced to reboot your computer, losing any unsaved work. For example, DON'T try the
# following (unless you're prepared for a long wait, have saved your work, and the lab isn't busy).
> # evalf(Int(g, x = Pi/4 .. 3*Pi/8));
# In some cases, even when there is an "exact" formula for the antiderivative, that formula could be so
# complicated that numerical integration may be faster and more accurate.
> timer:= time():\
int(exp(x)/(x^2 - x + 1)^2, x = 0 .. 1);\
evalf(");\
time()-timer;
1/2
1/3 exp(1) + 1/3 exp(1/2) %1 Ei(1, - 1/2 - 1/2 I 3 )
1/2 1/2
+ 2/9 I 3 exp(1/2) %2 Ei(1, - 1/2 + 1/2 I 3 )
1/2
+ 1/3 exp(1/2) %2 Ei(1, - 1/2 + 1/2 I 3 )
1/2 1/2
- 2/9 I 3 exp(1/2) %1 Ei(1, - 1/2 - 1/2 I 3 ) + 1/3
1/2
- 1/3 exp(1/2) %1 Ei(1, 1/2 - 1/2 I 3 )
1/2 1/2
- 2/9 I 3 exp(1/2) %2 Ei(1, 1/2 + 1/2 I 3 )
1/2
- 1/3 exp(1/2) %2 Ei(1, 1/2 + 1/2 I 3 )
1/2 1/2
+ 2/9 I 3 exp(1/2) %1 Ei(1, 1/2 - 1/2 I 3 )
1/2
%1 := exp(- 1/2 I 3 )
1/2
%2 := exp(1/2 I 3 )
-9
2.515476400 - .2*10 I
2.550
> timer:= time():\
evalf(Int(exp(x)/(x^2 - x + 1)^2, x = 0 .. 1));\
time()-timer;
2.515476400
.716
# In general, evalf(Int(...)) is quite reliable and reasonably fast, except for that one difficulty with
# undetected singularities, and some other circumstances we'll see where it has difficulties.
