# 1.15 Lesson 15

```# 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).
> 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
# 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.

```