# 2.15 REMOVING SINGULARITIES

```# WORKSHEET # 26
#
# REMOVING SINGULARITIES
#
#
--------------------------------------------------------------------------------
> restart;
--------------------------------------------------------------------------------
# In this worksheet we consider numerical evaluation of improper integals.
--------------------------------------------------------------------------------
> J1:=Int(1/sqrt(x^5+1),x=0..infinity);

infinity
/
|           1
J1 :=    |      ----------- dx
|        5     1/2
/       (x  + 1)
0
--------------------------------------------------------------------------------
# This integral is improper because of the infinite range of integration.  Numerical
# techniques such as the Trapezoidal rule and Simpson's rule do not work, so we try to
# change the integral into a proper one.  In the student package you will find the change of
# variables command.
--------------------------------------------------------------------------------
> with(student);

[D, Doubleint, Int, Limit, Lineint, Sum, Tripleint, changevar, combine,

completesquare, distance, equate, extrema, integrand, intercept, intparts,

isolate, leftbox, leftsum, makeproc, maximize, middlebox, middlesum,

midpoint, minimize, powsubs, rightbox, rightsum, showtangent, simpson,

slope, trapezoid, value]
--------------------------------------------------------------------------------
> ?changevar
--------------------------------------------------------------------------------
# The following change of variables changes the range of integration to a finite interval.
--------------------------------------------------------------------------------
> 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               /
--------------------------------------------------------------------------------
# You would think that Maple would be smart enough to cancel some t's, but it's not!  The
# simplify command has an  option called "symbolic" which allows one to replace sqrt(t^2)
# by t.  In this case t is non-negative, so this is OK.
--------------------------------------------------------------------------------
> J2:=simplify(J2,symbolic);

1
/                  1/2
|                  t
J2 :=  |  ----------------------------------- dt
|                 2       3      4 1/2
/   (1 - 5 t + 10 t  - 10 t  + 5 t )
0
--------------------------------------------------------------------------------
# This integral is no longer improper (the integrand behaves like sqrt(t) near t=0), but
# numerical techniques such as the trapezoidal rule do not work too well since the
# integrand is not differentiable at t=0.  Recall that the error in the trapezoidal rule is a
# constant times some second derivative, assuming that the integrand f(x) is continuously
# differentiable.  Let T[n] denote the trapezoidal rule.  Then
--------------------------------------------------------------------------------
>  Int(f(x),x=a..b)=T[n]-D(D(f))(u)*(b-a)^3/(12*n^2);\

b
/                        (2)              3
|                        D   (f)(u) (b - a)
|  f(x) dx = T[n] - 1/12 -------------------
|                                  2
/                                  n
a
--------------------------------------------------------------------------------
# Simpson's rule also does not work too well.  The next change of variables together with a
# simplification will make the integrand infinitely differentiable.
--------------------------------------------------------------------------------
> J3:=changevar(t=s^2,J2,s);

1
/                   2 1/2
|                  (s )    s
J3 :=  |  2 ------------------------------------ ds
|            2       4       6      8 1/2
/     (1 - 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
--------------------------------------------------------------------------------
> 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 )
--------------------------------------------------------------------------------
# This function actually has infinitely many derivatives on the interval t=0..1.  For example,
# the Maclaurin series expansion about s=0 is given by the following command.
--------------------------------------------------------------------------------
> ?series
--------------------------------------------------------------------------------
> series(f(s),s=0);

2      4      6
2 s  + 5 s  + O(s )
--------------------------------------------------------------------------------
# We can get more terms by including an option in the series command.
--------------------------------------------------------------------------------
> series(f(s),s=0,10);

2      4         6          8      10
2 s  + 5 s  + 35/4 s  + 105/8 s  + O(s  )
--------------------------------------------------------------------------------
# We can apply any of our numerical integration techniques to J3.
--------------------------------------------------------------------------------
>  h:=n->(b-a)/n;\

b - a
h := n -> -----
n
--------------------------------------------------------------------------------
>  X:=(i,n)->a+i*(b-a)/n;

i (b - a)
X := (i,n) -> a + ---------
n
--------------------------------------------------------------------------------
> Trap:=n->(h(n)/2)*sum(f(X(i,n))+f(X(i+1,n)),i=0..n-1);

/n - 1                              \
|-----                              |
| \                                 |
Trap := n -> 1/2 h(n) |  )   (f(X(i, n)) + f(X(i + 1, n)))|
| /                                 |
|-----                              |
\i = 0                              /
--------------------------------------------------------------------------------
> Simpson:=n->evalf((h(n)/3)*sum(f(X(2*i,n))+4*f(X(2*i+1,n))+f(X(2*i+2,n)),i=0..n/2-1));

Simpson := n -> evalf(

/1/2 n - 1                                                       \
|  -----                                                         |
|   \                                                            |
1/3 h(n) |    )     (f(X(2 i, n)) + 4 f(X(2 i + 1, n)) + f(X(2 i + 2, n)))|)
|   /                                                            |
|  -----                                                         |
\  i = 0                                                         /
--------------------------------------------------------------------------------
> a:=0;b:=1;

a := 0

b := 1
--------------------------------------------------------------------------------
> Simpson(20);

1.549619113
--------------------------------------------------------------------------------
> Simpson(40);

1.549696107
--------------------------------------------------------------------------------
# The integral J1 =J3 can be computed exactly in terms of the Gamma function.
--------------------------------------------------------------------------------
> 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
--------------------------------------------------------------------------------
# The true error in Simpson[40] can be computed.
--------------------------------------------------------------------------------
> "-Simpson(40);

-6
.168*10
--------------------------------------------------------------------------------
# By the theory of Richardson extrapolation (Simpson[40]-Simpson[20])/15 equals the first
# term in the error when J3 is approximated by Simpson[40].
--------------------------------------------------------------------------------
> (Simpson(40)-Simpson(20))/15;

-5
.51329*10
--------------------------------------------------------------------------------
# The first term is bigger than the actual error, which is what we should expect.
--------------------------------------------------------------------------------
> Simpson(40)+(Simpson(40)-Simpson(20))/15;

1.549701240
--------------------------------------------------------------------------------
# We see that Simpson(40)+(Simpson(40)-Simpson(20))/15 is greater than the actual
# value and Simpson(40)  is less.
--------------------------------------------------------------------------------
# Another method for removing singularities, on a finite interval, is by subtracting the
# singular part.  We will illustrate this technique by reconsidering J2.
--------------------------------------------------------------------------------
> J2;

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

1/2
t
g := -----------------------------------
2       3      4 1/2
(1 - 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
--------------------------------------------------------------------------------
# In other words in a neighbourhood of t=0 the function g(t) has this expansion.  Actually,
# this expansion is valid for t=0..1.  For Simpson's rule we want to have 4 continuous
# derivatives.  Thus the terms from t^(9/2) onward are OK, but the earlier terms are not.
--------------------------------------------------------------------------------
> 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
--------------------------------------------------------------------------------
# The idea is to split off the singular part, integrate this exactly, and then integrate the
# non-singular part (i.e. g(t)-the singular part) by some numerical technique.  This idea
# works quite well since the integral of the singular part is quite easy to detemine.
--------------------------------------------------------------------------------
> int(singularpart,t=0..1);

35/8
--------------------------------------------------------------------------------
> Int(g-singularpart,t=0..1);

Int(

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

t = 0 .. 1)
--------------------------------------------------------------------------------
# Thus the integral J2 is the sum of these last 2 integrals.  We will call it J4.
--------------------------------------------------------------------------------
> J4:=int(singularpart,t=0..1)+Int(g-singularpart,t=0..1);

J4 := 35/8 + Int(

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

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

f := t ->

1/2
t                      1/2        3/2         5/2           7/2
----------------------------------- - t    - 5/2 t    - 35/8 t    - 105/16 t
2       3      4 1/2
(1 - 5 t + 10 t  - 10 t  + 5 t )
--------------------------------------------------------------------------------
> 35/8+Simpson(40);

1.549696021
--------------------------------------------------------------------------------
> 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
--------------------------------------------------------------------------------
# Another method of dealing with singularities is integration by parts.  It is not even clear
# that the next integral converges.
--------------------------------------------------------------------------------
> J5:=Int(x*sin(x^3),x=0..infinity);

infinity
/
|             3
J5 :=    |      x sin(x ) dx
|
/
0
--------------------------------------------------------------------------------
# The next change of variables at least makes the sine term easier, but still leaves open the
# question of convergence.
--------------------------------------------------------------------------------
> J5:=changevar(x^3=t,J5,t);

infinity
/
|          sin(t)
J5 :=    |      1/3 ------ dt
|            1/3
/            t
0
--------------------------------------------------------------------------------
# This last integral is improper for 2 reasons: the range of integration is infinite and there
# exists a singularity at t=0.  This suggest that we split the integral into 2 pieces.
--------------------------------------------------------------------------------
> Int((1/3)*sin(t)/t^(1/3),t=0..infinity)=Int((1/3)*sin(t)/t^(1/3),t=0..1)+Int((1/3)*sin(t)/t^(1/3),t=1..infinity);

infinity                   1                  infinity
/                       /                     /
|          sin(t)       |      sin(t)         |          sin(t)
|      1/3 ------ dt =  |  1/3 ------ dt +    |      1/3 ------ dt
|            1/3        |        1/3          |            1/3
/            t          /        t            /            t
0                       0                     1
--------------------------------------------------------------------------------
# The integral from t=0..1 presents no real problem, we can handle it in the same way we
# did for J2 above.  Actually we will just evaluate it directly, the default method of
# numerical integration in Maple knows how to handle the singularity at t=0.
--------------------------------------------------------------------------------
> J5a:=evalf(Int((1/3)*sin(t)/t^(1/3),t=0..1));

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

infinity
/
|          sin(t)
J5b :=    |      1/3 ------ dt
|            1/3
/            t
1
--------------------------------------------------------------------------------
# Now we do repeated integration by parts.
--------------------------------------------------------------------------------
> ?intparts
--------------------------------------------------------------------------------
> Int(u(x)*diff(v(x),x),x=A..B);

B
/
|       /  d      \
|  u(x) |---- v(x)| dx
|       \ dx      /
/
A
--------------------------------------------------------------------------------
> Int(u(x)*diff(v(x),x),x=A..B)=intparts(",u(x));

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

infinity
/
|          cos(t)
J5b := 1/3 cos(1) -    |      1/9 ------ dt
|            4/3
/            t
1
--------------------------------------------------------------------------------
# The new integral has t^(4/3) in the denominator and so clearly converges.  We repeat
# integration by parts.
--------------------------------------------------------------------------------
> J5b:=intparts(J5b,t^(-4/3));

infinity
/
|             sin(t)
J5b := 1/3 cos(1) + 1/9 sin(1) +    |      - 4/27 ------ dt
|               7/3
/               t
1
--------------------------------------------------------------------------------
# Now we have t^(7/3) in the denominator, thus faster convergence.  We can even do a loop
# consisting of integration by parts.
--------------------------------------------------------------------------------
> for j from 1 to 8 do J5b:=intparts(J5b,t^(-j-4/3)) od;j:='j':

infinity
/
|         28  cos(t)
J5b := 5/27 cos(1) + 1/9 sin(1) -    |      - ---- ------ dt
|         81    10/3
/               t
1

infinity
/
19              |      280 sin(t)
J5b := 5/27 cos(1) - ---- sin(1) +    |      --- ------ dt
81              |      243   13/3
/            t
1

infinity
/
325           19              |      3640 cos(t)
J5b := --- cos(1) - ---- sin(1) -    |      ---- ------ dt
243           81              |       729   16/3
/             t
1

infinity
/
325          3469             |        58240 sin(t)
J5b := --- cos(1) + ---- sin(1) +    |      - ----- ------ dt
243           729             |         2187   19/3
/                t
1

infinity
/
55315          3469             |        1106560 cos(t)
J5b := - ----- cos(1) + ---- sin(1) -    |      - ------- ------ dt
2187           729             |          6561    22/3
/                  t
1

infinity
/
55315          1075339             |      24344320 sin(t)
J5b := - ----- cos(1) - ------- sin(1) +    |      -------- ------ dt
2187            6561              |        19683    25/3
/                 t
1

infinity
/
23846485          1075339             |      608608000 cos(t)
J5b := -------- cos(1) - ------- sin(1) -    |      --------- ------ dt
19683             6561              |        59049     28/3
/                  t
1

infinity
/
23846485          598929949             |        17041024000 sin(t)
J5b := -------- cos(1) + --------- sin(1) +    |      - ----------- ------ dt
19683             59049               |           177147     31/3
/                      t
1
--------------------------------------------------------------------------------
# Imagine what it would be like to have to do these calculations by hand, and get the
# correct answer!  Each of the integrals here still has the same problem, namely the range
# of integration is infinite.  So we make the substitution s=1/t.
--------------------------------------------------------------------------------
> J5b:=changevar(s=1/t,J5b,s);

1
/
23846485          598929949           |    17041024000  25/3
J5b := -------- cos(1) + --------- sin(1) +  |  - ----------- s     sin(1/s) ds
19683             59049             |       177147
/
0
--------------------------------------------------------------------------------
> f:=integrand(J5b,s);

17041024000  25/3
f := - ----------- s     sin(1/s)
177147
--------------------------------------------------------------------------------
# This new integrand is now 4 times continuously differentiable at s=0, but not 5 times.
# Thus we can apply Simpson's rule.
--------------------------------------------------------------------------------
> diff(",s\$4);

2849259212800000  13/3            455881474048000  10/3
- ---------------- s     sin(1/s) + --------------- s     cos(1/s)
14348907                          4782969

10360942592000  7/3            1090625536000  4/3
+ -------------- s    sin(1/s) - ------------- s    cos(1/s)
531441                         531441

17041024000  1/3
- ----------- s    sin(1/s)
177147
--------------------------------------------------------------------------------
> diff(",s);

37040369766400000  10/3            7408073953280000  7/3
- ----------------- s     sin(1/s) + ---------------- s    cos(1/s)
43046721                          14348907

673461268480000  4/3            35445329920000  1/3
+ --------------- s    sin(1/s) - -------------- s    cos(1/s)
4782969                         1594323

1107666560000 sin(1/s)   17041024000 cos(1/s)
- ------------- -------- + ----------- --------
531441       2/3        177147      5/3
s                      s
--------------------------------------------------------------------------------
# I experimented with the "integration by parts loop" so that I finally got an integrand that
# was 4 times continuously differentiable.
--------------------------------------------------------------------------------
> J5c:=op(1,J5b)+op(2,J5b);

23846485          598929949
J5c := -------- cos(1) + --------- sin(1)
19683             59049
--------------------------------------------------------------------------------
> evalf(J5c);

9189.573179
--------------------------------------------------------------------------------
> f:=unapply(f,s);

25/3
f := s -> - 17041024000/177147 s     sin(1/s)
--------------------------------------------------------------------------------
> J5d:=Simpson(40);
Error, (in sum) division by zero

--------------------------------------------------------------------------------
# The integrand has a removable singularity at s=0 that Maple can not handle, even though
# the integrand is 4 times continuously differentiable there.  One way around this is to
# integrate  over a slightly smaller interval, say s=(10)^(-10)..1.
--------------------------------------------------------------------------------
> a:=(10)^(-10);

a := 1/10000000000
--------------------------------------------------------------------------------
> J5d:=Simpson(40);

J5d := -9189.416599
--------------------------------------------------------------------------------
# The original integral is now just the sum of J5a, J5c and J5d.
--------------------------------------------------------------------------------
> Int(x*sin(x^3),x=0..infinity)=evalf(J5a+J5c+J5d);

infinity
/
|             3
|      x sin(x ) dx = .341910
|
/
0
--------------------------------------------------------------------------------

```