# Lesson 15. Symbolic and Numerical Integration # ----------------------------------- # # Last time we saw that there can be more to a definite integral than F(b) - F(a). > restart:\ J:= int(sqrt(tan(x)), x); 1/2 1/2 1/2 1/2 1/2 2 tan(x) 1/2 tan(x) + 2 tan(x) + 1 J := 1/2 2 arctan(--------------) - 1/2 2 ln(---------------------------) 1 - tan(x) 2 1/2 (1 + tan(x) ) > int(sqrt(tan(x)), x=0 .. Pi/3); 1/2 1/4 1/2 2 3 1/2 1/2 1/2 1/4 1/2 - 1/2 2 arctan(----------) - 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 ) 1/2 - 1 + 3 1/2 + 1/2 2 Pi > evalf("); .787779048 > subs(x=Pi/3, J) - subs(x = 0, J); 1/2 1/2 1/2 2 tan(1/3 Pi) 1/2 2 arctan(-------------------) 1 - tan(1/3 Pi) 1/2 1/2 1/2 tan(1/3 Pi) + 2 tan(1/3 Pi) + 1 - 1/2 2 ln(-------------------------------------) 2 1/2 (1 + tan(1/3 Pi) ) 1/2 1/2 1/2 2 tan(0) - 1/2 2 arctan(--------------) 1 - tan(0) 1/2 1/2 1/2 tan(0) + 2 tan(0) + 1 + 1/2 2 ln(---------------------------) 2 1/2 (1 + tan(0) ) > "; 1/2 1/4 1/2 2 3 1/2 1/2 1/2 1/4 1/2 1/2 2 arctan(---------) - 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 ) 1/2 1 - 3 > evalf("); -1.433662421 # The discontinuities of the antiderivative are found with the "discont" command. This must first be read # in to Maple (although that's not necessary if you evaluate a definite integral first). > readlib(discont): > discont(J,x); {1/4 Pi + Pi _Z~, Pi _Z7~ + 1/2 Pi} # A word about "arctan". The discontinuity at Pi/4 in this case could be removed with a better choice of # formula for the antiderivative, based on arccot instead of arctan. > J; 1/2 1/2 1/2 1/2 1/2 2 tan(x) 1/2 tan(x) + 2 tan(x) + 1 1/2 2 arctan(--------------) - 1/2 2 ln(---------------------------) 1 - tan(x) 2 1/2 (1 + tan(x) ) > Jp:= sqrt(2)/2*arccot((1-tan(x))/(sqrt(2*tan(x)))) + op(2,J); 1/2 1/2 (1 - tan(x)) 2 Jp := 1/2 2 arccot(1/2 -----------------) 1/2 tan(x) 1/2 1/2 1/2 tan(x) + 2 tan(x) + 1 - 1/2 2 ln(---------------------------) 2 1/2 (1 + tan(x) ) # This has no discontinuity when tan(x) = 1. It does have one when tan(x)=0, but that's at an endpoint of # our interval. > subs(x=Pi/3, Jp) - limit(Jp, x=0,right); 1/2 1/2 (1 - tan(1/3 Pi)) 2 1/2 2 arccot(1/2 ----------------------) 1/2 tan(1/3 Pi) 1/2 1/2 1/2 tan(1/3 Pi) + 2 tan(1/3 Pi) + 1 - 1/2 2 ln(-------------------------------------) 2 1/2 (1 + tan(1/3 Pi) ) > "; 1/2 1/2 1/2 (1 - 3 ) 2 1/2 2 arccot(1/2 ---------------) 1/4 3 1/2 1/2 1/2 1/4 1/2 - 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 ) > evalf("); .7877790480 # There is an option to "int" that omits the search for discontinuities in the interval of integration. This # may sometimes speed things up, but use it at your own risk! > int(sqrt(tan(x)), x=0 .. Pi/3, continuous); 1/2 1/4 1/2 2 3 1/2 1/2 1/2 1/4 1/2 - 1/2 2 arctan(----------) - 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 ) 1/2 - 1 + 3 # What if we had an example where Maple couldn't identify the discontinuities? > g:= normal(diff(1/(2*x-tan(x)),x));\ 2 - 1 + tan(x) g := ----------------- 2 (- 2 x + tan(x)) > J:= int(g,x); 1 J := - -------------- - 2 x + tan(x) > discont(J,x); {Pi _Z11~ + 1/2 Pi, 0} # It didn't see the discontinuity at the point between 0 and Pi/2 where # tan(x) = 2*x. > plot({2*x, tan(x)}, x=0..Pi/2, 0..Pi); > int(g, x=Pi/4 .. 3*Pi/8) = subs(x=3*Pi/8, J) - subs(x=Pi/4, J); 4 2 1 1 ----------------- - ------ = - ---------------------- + ---------------------- 1/2 Pi - 2 - 3/4 Pi + tan(3/8 Pi) - 1/2 Pi + tan(1/4 Pi) 3 Pi - 4 2 - 4 > evalf("); -18.98765027 = -18.98765042 # Why is this answer absurd? What should the answer be? # Apart from such mishaps, Maple can do improper integrals. It checks for convergence or divergence. # Here are some you probably met in Math 101. > int(1/x^2, x= 1 .. infinity); 1 > int(1/x, x= -1 .. 1); 1 / | | 1/x dx | / -1 > int(exp(-a*x^2), x=-infinity .. infinity); infinity / | 2 | exp(- a x ) dx | / - infinity # Why didn't that one work? Because nobody told Maple that a > 0. If a < 0 the integral doesn't converge. # Don't expect Maple to read your mind! And remember that variables are allowed to be complex unless # otherwise specified, so replacing a with a^2 wouldn't have worked either. When you want to make an # assumption on a variable, you can use "assume". > assume(a > 0); > int(exp(-a*x^2), x=-infinity .. infinity); 1/2 Pi ----- 1/2 a~ # Maple reminds you that there is an assumption on the variable "a" by printing it as "a~". # Here are some that you may not have met. Maple evaluates # them correctly. > int(sin(x)/x, x=-infinity .. infinity); Pi > int((ln(x))^2/(1+x^2), x= 0 .. infinity); 3 1/8 Pi # In this case the answer didn't come from a formula for the antiderivative, since in fact Maple can't find # an antiderivative. # In Math 300 and 301 you learn how to evaluate integrals such as this using residues. > int((ln(x))^2/(1+x^2), x); / 2 | ln(x) | ------ dx | 2 / 1 + x # On the other hand... > int(sin(x)/(1+x^2), x=-infinity .. infinity); I sinh(1) Pi # Oops again. This is another case where Maple's formula for the antiderivative has a discontinuity, this # time involving complex numbers. It's actually something you'll meet in Math 301, called a branch cut. # # Conclusion: don't put too much trust in "int" for definite integrals. It's good to be able to check its # answers. This can often be done using numerical integration, which we'll look at next. # We can use "evalf" on an integral, even if Maple can't do the integration symbolically. > int(exp(sin(x)), x=0..1); 1 / | | exp(sin(x)) dx | / 0 > evalf("); 1.631869608 # This is using numerical techniques (which we will discuss). In fact we can avoid the attempt at symbolic # integration: > evalf(Int(exp(sin(x)), x=0 .. 1)); 1.631869608 # evalf(Int(...)) can be used for improper integrals as well. > evalf(Int(1/x, x=-1 .. 1)); Error, (in evalf/int) unable to handle singularity > evalf(Int(1/x^2, x= 1 .. infinity)); 1. > evalf(Int(ln(x)^2/(1+x^2), x=0 .. infinity)); 3.875784585 # Here is the answer we obtained for that one with "int". > evalf(Pi^3/8); 3.875784586 # One thing evalf(Int(...)) can't handle is a symbolic variable other than the variable of integration. > evalf(Int(exp(-a*x^2), x = -infinity .. infinity)); Error, (in evalf/int) function does not evaluate to numeric > evalf(Int(exp(-x^2), x = -infinity .. infinity)) = evalf(sqrt(Pi)); 1.772453851 = 1.772453851 # Warning!! If you try evalf(Int(...)) on an improper integral that diverges, but Maple doesn't detect that it # diverges, the result could be a VERY long attempt by Maple to approximate the integral numerically. # Since the desired accuracy can never be attained, Maple might keep trying forever. In some cases, # especially on your home PC, the usual methods of interrupting Maple's computation might not work, and # you could be forced to reboot your computer, losing any unsaved work. For example, DON'T try the # following (unless you're prepared for a long wait, have saved your work, and the lab isn't busy). > # evalf(Int(g, x = Pi/4 .. 3*Pi/8)); # In some cases, even when there is an "exact" formula for the antiderivative, that formula could be so # complicated that numerical integration may be faster and more accurate. > timer:= time():\ int(exp(x)/(x^2 - x + 1)^2, x = 0 .. 1);\ evalf(");\ time()-timer; 1/2 1/3 exp(1) + 1/3 exp(1/2) %1 Ei(1, - 1/2 - 1/2 I 3 ) 1/2 1/2 + 2/9 I 3 exp(1/2) %2 Ei(1, - 1/2 + 1/2 I 3 ) 1/2 + 1/3 exp(1/2) %2 Ei(1, - 1/2 + 1/2 I 3 ) 1/2 1/2 - 2/9 I 3 exp(1/2) %1 Ei(1, - 1/2 - 1/2 I 3 ) + 1/3 1/2 - 1/3 exp(1/2) %1 Ei(1, 1/2 - 1/2 I 3 ) 1/2 1/2 - 2/9 I 3 exp(1/2) %2 Ei(1, 1/2 + 1/2 I 3 ) 1/2 - 1/3 exp(1/2) %2 Ei(1, 1/2 + 1/2 I 3 ) 1/2 1/2 + 2/9 I 3 exp(1/2) %1 Ei(1, 1/2 - 1/2 I 3 ) 1/2 %1 := exp(- 1/2 I 3 ) 1/2 %2 := exp(1/2 I 3 ) -9 2.515476400 - .2*10 I 2.550 > timer:= time():\ evalf(Int(exp(x)/(x^2 - x + 1)^2, x = 0 .. 1));\ time()-timer; 2.515476400 .716 # In general, evalf(Int(...)) is quite reliable and reasonably fast, except for that one difficulty with # undetected singularities, and some other circumstances we'll see where it has difficulties.