# 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
# quadrature routine "ccquad", avoiding the singularity-handling
# 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
>
