# 1.20 Lesson 20

```# Lesson 20.  Removing Singularities
# =======================
#
# The various numerical integration methods we have considered
# all depend on the integrand being a "smooth" function.  The
# error estimates for the Trapezoid and Midpoint Rules, for
# example, assume that f'' is continuous, while that for
# Simpson's Rule assumes a continuous 4th derivative.  Although
# these rules will "work" for functions that do not have so many
# continuous derivatives, in the sense that the limit as n ->
# infinity is correct, they may approach that limit quite
# slowly.  We must therefore modify our approach in dealing with
# such functions.
#
# There are even more serious problems in dealing with improper
# integrals.  Here we often can't even get started (e.g. for an
# integral over an infinite interval, or an integrand that is
# undefined at an endpoint).
# Of course, some improper integrals diverge, and these are
# outside the reach of any numerical integration method (but we
# want to be sure to detect this divergence before trying
# numerical methods).
#
# Our first tool for dealing with such integrals is a change of
# variables, by which we try to transform the integral into a
# proper integral with a smooth integrand.  In some cases we
# break up the interval of integration into several pieces to
# make this simpler.
#
> restart:
> J1:= Int(1/sqrt(x^5+1), x=0..infinity);

infinity
/
|           1
J1 :=    |      ----------- dx
|        5     1/2
/       (x  + 1)
0
# This is a "type I" improper integral, i.e. the integral of a
# continuous function on an infinite interval.  We know it
# converges since the integrand is less than x^(-5/2) and
# int(x^(-5/2), x=1..infinity) converges.
# The simplest way to transform an infinite interval to a finite
# one is the change of variable t = 1/x, which takes
# [a,infinity) to (0, 1/a] for a > 0.  Here we want to use a=0,
# but that's only a minor annoyance: we can use t = 1/(x+1),
# which takes [0,infinity) to (0,1].
#
> with(student):
> J2:=changevar(1/(x+1)=t, J1, t);

1
/
|            1
J2 :=  |  -------------------- dt
|  /       5    \1/2
/   |(1 - t)     |     2
0   |-------- + 1|    t
|    5       |
\   t        /
> J2:=simplify(J2);

1
/
|                     1
J2 :=  |  -------------------------------------- dt
|  /              2       3      4\1/2
/   |1 - 5 t + 10 t  - 10 t  + 5 t |     2
0   |------------------------------|    t
|               5              |
\              t               /
# Why won't it take the t^5 out of the square root?  Release 2
# would: it would simplify sqrt(t^2) to t, but that would run
# into trouble when t wasn't positive.  In this case t is in the
# interval 0 to 1, but Maple doesn't take that into account.
# The "simplify" command has an option "symbolic", that is
# useful when you don't need to worry about signs.
> J2:=simplify(J2,symbolic);

1
/
|                    t
J2 :=  |  ------------------------------------ dt
|          2       3       4      5 1/2
/   (t - 5 t  + 10 t  - 10 t  + 5 t )
0
# We might take out a sqrt(t) from numerator and denominator.
# Then the integral is a proper integral.  But it's not smooth
# at 0: it's approximately sqrt(t) there.  To smooth it out, try
# another change of variables.
> J3:=changevar(t=s^2, J2, s);

1
/                      3
|                      s
J3 :=  |  2 -------------------------------------- ds
|      2      4       6       8      10 1/2
/     (s  - 5 s  + 10 s  - 10 s  + 5 s  )
0
> J3:=simplify(J3,symbolic);

1
/                   2
|                   s
J3 := 2  |  ------------------------------------ ds
|          2       4       6      8 1/2
/   (1 - 5 s  + 10 s  - 10 s  + 5 s )
0
# No problem now: any of our numerical methods will handle this.
> h:= n ->(b-a)/n:
> X:= (i,n) -> a + i*(b-a)/n:
> a:=0: b:= 1:
> f:= unapply(2*integrand(J3),s);

2
s
f := s -> 2 ------------------------------------
2       4       6      8 1/2
(1 - 5 s  + 10 s  - 10 s  + 5 s )
# For convenience, I'm making a slight change in "simpson", by
# putting in an "evalf".
> simpson:= n ->
>    evalf(h(n) * sum(1/3*f(X(2*i,n))+4/3*f(X(2*i+1,n))+1/3*f(X(2*i+2,n)), i=0..n/2-1));

simpson := n -> evalf(h(n)

/1/2 n - 1                                                                 \
|  -----                                                                   |
|   \                                                                      |
|    )     (1/3 f(X(2 i, n)) + 4/3 f(X(2 i + 1, n)) + 1/3 f(X(2 i + 2, n)))|
|   /                                                                      |
|  -----                                                                   |
\  i = 0                                                                   /

)
> simpson(40);

1.549696108
> (simpson(40)-simpson(20))/15;

-5
.51330*10
> Jexact:= value(J1);

3/2  1/2       1/2 1/2          3/2  1/2  1/2       1/2 1/2
Pi    2    (5 + 5   )           Pi    5    2    (5 + 5   )
Jexact := 1/5 ------------------------ - 1/25 -----------------------------
GAMMA(4/5) GAMMA(7/10)             GAMMA(4/5) GAMMA(7/10)
> evalf(");

1.549696275
> evalf(Jexact-simpson(40));

-6
.167*10
# Another method for removing singularities (on a finite
# interval) is subtraction.  Let's go back to J2.
> J2;

1
/
|                    t
|  ------------------------------------ dt
|          2       3       4      5 1/2
/   (t - 5 t  + 10 t  - 10 t  + 5 t )
0
> g:= integrand(J2);

t
g := ------------------------------------
2       3       4      5 1/2
(t - 5 t  + 10 t  - 10 t  + 5 t )
> series(g,t=0);

1/2        3/2         5/2   105  7/2   1155  9/2   2875  11/2      13/2
t    + 5/2 t    + 35/8 t    + --- t    + ---- t    + ---- t     + O(t    )
16         128         256
# For Simpson's Rule, we want the integrand to have a continuous
# fourth derivative.  This means we must remove the t^(1/2),
# t^(3/2), t^(5/2) and t^(7/2) terms, but the t^(9/2) should be
# OK.
> singularpart:= t^(1/2)+5/2*t^(3/2)+35/8*t^(5/2)+105/16*t^(7/2);

1/2        3/2         5/2   105  7/2
singularpart := t    + 5/2 t    + 35/8 t    + --- t
16
> J4:= int(singularpart,t=0..1)+Int(g-singularpart,t=0..1);

J4 := 35/8 + Int(

t                     1/2        3/2         5/2   105  7/2
------------------------------------ - t    - 5/2 t    - 35/8 t    - --- t   ,
2       3       4      5 1/2                                  16
(t - 5 t  + 10 t  - 10 t  + 5 t )

t = 0 .. 1)
> f:= unapply(g-singularpart, t);

f := t ->

t                     1/2        3/2         5/2           7/2
------------------------------------ - t    - 5/2 t    - 35/8 t    - 105/16 t
2       3       4      5 1/2
(t - 5 t  + 10 t  - 10 t  + 5 t )
> 35/8 + simpson(40);
Error, (in sum) division by zero

# Where could we have a division by 0?
# In f(0), which is not defined.  Well, the easy way out is to
# let "a" be slightly greater than 0, not enough to make a real
# difference.

> a:= 10^(-10);

a := 1/10000000000
> 35/8 + simpson(40);

1.549696021
> evalf(Jexact)-";

-6
.254*10
# A third method for dealing with singularities is integration
# by parts.
> J5:= Int(x*sin(x^3), x=0..infinity);

infinity
/
|             3
J5 :=    |      x sin(x ) dx
|
/
0
# It's not obvious that this even converges.  We could first
# transform to a finite interval.
> changevar(t=1/(x+1),J5,t);

3
(1 - t)
1 (1 - t) sin(--------)
/                 3
|                 t
|  --------------------- dt
|             3
/             t
0
# Still not obvious that it converges.  A different change of
# variables will be more useful to us, one that makes the "sin"
# part simpler.
> changevar(x^3=t,J5,t);

infinity
/
|          sin(t)
|      1/3 ------ dt
|            1/3
/            t
0
# This actually is improper at both ends, even though the
# original J5 was fine at 0.  To avoid the problems at 0, split
# J5 into two pieces.
# The integral from 0 to 1, say, is no problem.  I won't even
# bother using "simpson" on it, just "evalf".
> J5a:= evalf(Int(x*sin(x^3),x=0..1));

J5a := .1853301486
> J5b:= Int(x*sin(x^3),x=1..infinity);

infinity
/
|             3
J5b :=    |      x sin(x ) dx
|
/
1
> J:= changevar(x^3=t,J5b,t);

infinity
/
|          sin(t)
J :=    |      1/3 ------ dt
|            1/3
/            t
1
# The "intparts" command in the "student" package does an
# integration by parts.
> Int(u(x)*diff(v(x),x), x=A..B);

B
/
|       /  d      \
|  u(x) |---- v(x)| dx
|       \ dx      /
/
A
> intparts(", u(x));

B
/
|  /  d      \
u(B) v(B) - u(A) v(A) -  |  |---- u(x)| v(x) dx
|  \ dx      /
/
A
> J:=intparts(J,t^(-1/3));

infinity
/
|          cos(t)
J := 1/3 cos(1) -    |      1/9 ------ dt
|            4/3
/            t
1
# This is progress: t^(-4/3) goes to 0 as t-> infinity faster
# than t^(-1/3).  Fast enough, in fact, that we can now see that
# the integral converges by the Comparison Test.
> J:=intparts(J,t^(-4/3));

infinity
/
|             sin(t)
J := 1/3 cos(1) + 1/9 sin(1) +    |      - 4/27 ------ dt
|               7/3
/               t
1
> J:=intparts(J,t^(-7/3));

infinity
/
|         28  cos(t)
J := 5/27 cos(1) + 1/9 sin(1) -    |      - ---- ------ dt
|         81    10/3
/               t
1
> J:=intparts(J,t^(-10/3));

infinity
/
19              |      280 sin(t)
J := 5/27 cos(1) - ---- sin(1) +    |      --- ------ dt
81              |      243   13/3
/            t
1
# Can we change variables to a proper integral now?
> J:= changevar(s=1/t,J,s);

1
/
19            |  280  7/3
J := 5/27 cos(1) - ---- sin(1) +  |  --- s    sin(1/s) ds
81            |  243
/
0
# The integrand here has one continuous derivative but not two
# (each time you differentiate sin(1/s) or cos(1/s), you get a
# 1/s^2 factor).  So it's not yet time to call in "simpson".
# We'd need s^(25/3) instead of s^(7/3).  Another six
# integrations by parts should do it.
> J:=intparts(J,s^(13/3));

1
/
325           19            |  3640  10/3
J := --- cos(1) - ---- sin(1) -  |  ---- s     cos(1/s) ds
243           81            |   729
/
0
> for jj from 2 to 6 do J:=intparts(J,s^(jj+10/3)) od:
> J;

1
/
23846485          598929949           |    17041024000  25/3
-------- cos(1) + --------- sin(1) +  |  - ----------- s     sin(1/s) ds
19683             59049             |       177147
/
0
> J5c:= op(1,J)+op(2,J);

23846485          598929949
J5c := -------- cos(1) + --------- sin(1)
19683             59049
> f:= unapply(integrand(J),s);

25/3
f := s -> - 17041024000/177147 s     sin(1/s)
> J5d:= simpson(40);

J5d := -9189.416601
> eest:= (J5d-simpson(20))/15;

eest := .0488144
# If we want, say, error bounded less than 10^(-6), what n do we
# need?
> nest:= evalf((eest/10^(-6))^(1/4)*40);

nest := 594.5617691
> J5d:= simpson(600);

J5d := -9189.367611
> myJ:= J5a+evalf(J5c)+J5d;

myJ := .390898
> infolevel[evalf]:= 3; infolevel[int]:= 3;

infolevel[evalf] := 3

infolevel[int] := 3
> evalf(J5);
evalf/int:   entering
evalf/int/improper:   integrating on interval   0   ..   infinity
evalf/int/improper:   applying transformation   x   =   1/x
evalf/int/improper:   interval is   0   ..   1
for the two integrands:

1
sin(----)
3
3        x
x sin(x ), ---------
3
x
evalf/int/control:   integrating on   0   ..   1   the integrand

3
x sin(x )
evalf/int/control:   procedure for evaluation is
proc (x) evalf(x*sin(x^3)) end
evalf/int/control:   from ccquad, error =   .6599428644151505e-11
error tolerance =   .4633253715030e-10
integrand evals =   19
result =   .1853301486012
evalf/int/control:   integrating on   0   ..   1   the integrand

1
sin(----)
3
x
---------
3
x
evalf/int/control:   singularity at left end-point
evalf/int/transform:   series contains   {sin(1/x^3)}
evalf/int/findtransform:   no transform found
evalf/int/series:   integrating on   0 .. 1.   the series:

1
sin(----)
3
x
---------
3
x
int/indef:   first-stage indefinite integration
int/indef2:   second-stage indefinite integration
int/indef2:   applying change of variables
int/indef:   first-stage indefinite integration
int/indef:   first-stage indefinite integration
int/indef2:   second-stage indefinite integration
int/trigon:   case of integrand containing trigs
int/prptrig:   case ratpoly*trig(arg)
int/rischnorm:   enter Risch-Norman integrator
int/risch:   enter Risch integration
int/risch/algebraic1:
RootOfs should be algebraic numbers and functions
int/risch:   the field extensions are

2
RootOf(_Z  + 1)
[_X, exp(---------------)]
3
_X
int/risch:   Introduce the namings:

2
RootOf(_Z  + 1)
{_th[1] = exp(---------------)}
3
_X
unknown:   integrand is

1
_th[1] - ------
_th[1]
---------------
3
_X
unknown:   integrand expressed as

_th[1]        1
------ - ----------
3      3
_X     _X  _th[1]
int/risch/exppoly:   integrating

_th[1]        1
------ - ----------
3      3
_X     _X  _th[1]
int/risch/diffeq:   solving Risch d.e.  y' + f y = g  where f,g are:

2
RootOf(_Z  + 1)   1
- 3 ---------------, ---
4          3
_X         _X
int/risch/DEratpoly:   solving Risch d.e.  y' + f y = g  where f,g are:

2
RootOf(_Z  + 1)   1
- 3 ---------------, ---
4          3
_X         _X
int/risch/exppoly:   Risch d.e. has no solution
int/indef:   first-stage indefinite integration
int/indef:   first-stage indefinite integration
int/indef2:   second-stage indefinite integration
int/indef2:   applying change of variables
int/indef:   first-stage indefinite integration
int/indef:   first-stage indefinite integration
int/indef2:   second-stage indefinite integration
int/exp:   case of integrand containing exp
int/prpexp:   case ratpoly*exp(arg)
int/indef:   first-stage indefinite integration
int/indef:   first-stage indefinite integration
int/indef2:   second-stage indefinite integration
int/indef2:   applying change of variables
int/indef:   first-stage indefinite integration
int/indef:   first-stage indefinite integration
int/indef2:   second-stage indefinite integration
int/exp:   case of integrand containing exp
int/prpexp:   case ratpoly*exp(arg)
int/risch:   exit Risch integration
int/def:   definite integration
int/contour:   contour integration
evalf/int/series:   integration failed, result was:   int(1/x^3*sin(1/x
^3),x = 0 .. 1.)
evalf/int/control:   integrating on   0   ..   1   the integrand

1
sin(----)
3
x
---------
3
x
evalf/int/control:   procedure for evaluation is
proc (x) evalf(1/x^3*sin(1/x^3)) end
evalf/int/control:   applying double-exponential method
evalf/int/control:   integrating on   0   ..   1
evalf/int/control:   integrating on   0   ..   .5000000000000000
evalf/int/control:   integrating on   0   ..   .2500000000000000
evalf/int/control:   integrating on   0   ..   .1250000000000000
evalf/int/control:   integrating on   0   ..   .6250000000000000e-1
evalf/int/control:   integrating on   0   ..   .3125000000000000e-1
evalf/int/control:   integrating on   0   ..   .1562500000000000e-1
evalf/int/control:   integrating on   0   ..   .7812500000000000e-2
evalf/int/control:   integrating on   0   ..   .3906250000000000e-2
evalf/int/control:   integrating on   0   ..   .1953125000000000e-2
evalf/int/control:   integrating on   0   ..   .9765625000000000e-3
evalf/int/control:   integrating on   0   ..   .4882812500000000e-3
evalf/int/control:   integrating on   0   ..   .2441406250000000e-3
evalf/int/control:   integrating on   0   ..   .1220703125000000e-3
evalf/int/control:   integrating on   0   ..   .6103515625000000e-4
evalf/int/control:   integrating on   0   ..   .3051757812500000e-4
evalf/int/control:   integrating on   0   ..   .1525878906250000e-4

> ?int[numeric]
# - The default numerical method used is Clenshaw-Curtis
# quadrature.  If singularities are detected in or near the
# interval of integration then some techniques of symbolic
# analysis are used to deal with the singularities.  For
# problems with non-removable singularities, an adaptive
# double-exponential quadrature method is used.
#
# - The optional fourth parameter is a flag; its value may be
# any of:
#
#   _CCquad   indicates to use only the Clenshaw-Curtis
# routines.
#
#   _Dexp     indicates to use only the adaptive
# double-exponential routine "quadexp", avoiding Clenshaw-Curtis
# quadrature and the singularity-handling routines.
#
#   _NCrule   indicates to use only the adaptive Newton-Cotes
# rule "quanc8", which is a fixed-order method and hence is not
# effective for high precisions (e.g. Digits > 15).
#
> evalf(J);
evalf/int:   entering
evalf/int/control:   integrating on   0   ..   1   the integrand

17041024000  25/3
- ----------- s     sin(1/s)
177147
evalf/int/control:   singularity at left end-point
evalf/int/transform:   series contains   {sin(1/s), s^(25/3)}
evalf/int/singleft:   applying transformation   s    =   s^3
evalf/int/singleft:   interval is   0   ..   1.   for the integrand:

17041024000  27       1
- ----------- s   sin(----)
59049              3
s
evalf/int/control:   integrating on   0   ..   1.   the integrand

17041024000  27       1
- ----------- s   sin(----)
59049              3
s
evalf/int/control:   singularity at left end-point
evalf/int/transform:   series contains   {}
evalf/int/findtransform:   no transform found
evalf/int/series:   integrating on   0 .. .2261972975167   the series:

27
O(s  )
int/indef:   first-stage indefinite integration
evalf/int/series:   result =   0
evalf/int/control:   integrating on   .2261972975167   ..   1.
the integrand

17041024000  27       1
- ----------- s   sin(----)
59049              3
s
evalf/int/control:   procedure for evaluation is
proc (s) evalf(-17041024000/59049*s^27*sin(1/s^3)) end
evalf/int/control:   applying double-exponential method
evalf/int/control:   integrating on   .2261972975167000   ..   1
evalf/int/control:   integrating on   .2261972975167000   ..   .
6130986487583500
evalf/int/control:   integral on   .2261972975167000   ..   .
6130986487583500   =   .1012353861726033e-1
evalf/int/control:   integrating on   .6130986487583500   ..   1
evalf/int/control:   integral on   .6130986487583500   ..   1   =   -
9189.377730242926
evalf/int/control:   errest =   .2651421953474889e-10   abserr =   0
relerr =   .5000000000000000e-9   int(abs(f)) =   9189.518481224533
evalf/int/control:   from quadexp, error =   .2651421953475e-10
error tolerance =   .4594759240613e-5
integrand evals =   254
result =   -9189.367606704
evalf/int:   exiting

.205572
> J5a+";

.3909021486
>

```