# 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
--------------------------------------------------------------------------------
