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