Lesson 21: Singularities and Improper Integrals restart;

with(Student[Calculus1]): ApproximateInt(f(x), x=0..1, method = newtoncotes[8], partition=1); TR3:= unapply((2^6*ApproximateInt(f(x),x=0..1,method=newtoncotes[4],partition=n/4)-ApproximateInt(f(x),x=0..1,method=newtoncotes[4],partition=n/8))/(2^6-1),(f,n)); One thing that's not so nice about LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEsbmV3dG9uY290ZXNGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictSSNtbkdGJDYkUSI4RicvRjZRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPA== is that it has some negative coefficients. So you could have a function that's positive everywhere, but the Newton-Cotes approximation for its integral is negative. The LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEjVFJGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRicvRjNRJ25vcm1hbEYn rules don't have that problem: their coefficients are all positive. This implies that if LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbW9HRiQ2LVEifkYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0Y0LyUqc3ltbWV0cmljR0Y0LyUobGFyZ2VvcEdGNC8lLm1vdmFibGVsaW1pdHNHRjQvJSdhY2NlbnRHRjQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZDRi8=LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYtLUkjbWlHRiQ2JVEiY0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSYmbGVxO0YnRjlGO0Y+RkBGQkZERkZGSC9GS1EsMC4yNzc3Nzc4ZW1GJy9GTkZTRjUtRiw2JVEiZkYnRi9GMi1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRInhGJ0YvRjJGOUY5RjVGT0Y1LUYsNiVRImRGJ0YvRjJGOQ== on the interval LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSYmbGVxO0YnRjlGO0Y+RkBGQkZERkZGSC9GS1EsMC4yNzc3Nzc4ZW1GJy9GTkZTRjUtRiw2JVEieEYnRi9GMkY1Rk9GNS1GLDYlUSJiRidGL0YyRjk=, the TR rules will all give values between LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEiY0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRImJGJ0YvRjItRjY2LVEoJm1pbnVzO0YnRjlGO0Y+RkBGQkZERkZGSC9GS1EsMC4yMjIyMjIyZW1GJy9GTkZlbi1GLDYlUSJhRidGL0YyRjlGOS1GLDYjUSFGJ0Y5 and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiZEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1JKG1mZW5jZWRHRiQ2JC1GIzYoLUYsNiVRImJGJ0YvRjJGNS1GNjYtUSomdW1pbnVzMDtGJ0Y5RjtGPkZARkJGREZGRkgvRktRLDAuMjIyMjIyMmVtRicvRk5GZW5GNS1GLDYlUSJhRidGL0YyRjlGOUY5. Here's an example. f:= x -> exp(-100*(x-1/2)^2); evalf(ApproximateInt(f(x),x=0..1,method=newtoncotes[8],partition=1)); evalf(TR3(f,8)); evalf(int(f(x),x=0..1)); ApproximateInt(f(x),x=0..1,method=newtoncotes[8],partition=1,output=plot);

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 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEkZicnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnL0YzUSdub3JtYWxGJw== is continuous, while that for Simpson's Rule assumes a continuous 4th derivative. The higher-order rules require even more continuous derivatives. Although these rules will "work" for functions that do not have so many continuous derivatives, in the sense that the limit as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJuRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEnJnJhcnI7RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtRiw2JVEoJmluZmluO0YnRjRGN0Y+RitGPg== is correct, they may approach that limit quite slowly. We must therefore modify our approach in dealing with such functions. A typical function given by a closed-form expression will be nice and smooth almost everywhere in the interval we're integrating over, but might have a finite number of "bad" points where some derivative doesn't exist. We sometimes call these singularities. There are even more serious problems in dealing with improper integrals. Here the numerical methods we've used can't even get started (e.g. for an integral over an infinite interval, or an integrand that is undefined at some point). Of course, some improper integrals diverge, and these are outside the reach of any numerical integration method (but we would want to detect this divergence before trying numerical methods). Let's see what Maple does with an improper integral having a singularity. J := Int(1/sqrt(x^6+x), x=0..infinity); This integral is improper for two reasons: 1) the interval extends to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJmluZmluO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGMg== 2) the integrand goes to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJmluZmluO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGMg== as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RKCYjODU5NDtGJy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGUy1JI21uR0YkNiRRIjBGJ0ZARkBGK0ZARitGQA==. Maple can actually get an exact expression for this integral. value(J); evalf(%); But here we're interested in numerical integration. I'll use infolevel to take a peek at how Maple evaluates the integral numerically. infolevel[evalf]:= 1; evalf(J); It used a NAG routine for doing improper integrals. That wasn't very enlightening. I'll set Digits = 16 so it can't use NAG. Digits := 16: evalf(J); infolevel[evalf]:= 0; One of the main tools Maple uses (and we will too) for dealing with singularities is transformation (change of variables). The simplest way to transform an infinite interval to a finite one is the change of variable LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbWZyYWNHRiQ2KC1GIzYkLUkjbW5HRiQ2JFEiMUYnRj5GPi1GIzYkLUYsNiVRInhGJ0Y0RjdGPi8lLmxpbmV0aGlja25lc3NHRmZuLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmBvLyUpYmV2ZWxsZWRHRkJGPkYrRj4=, which takes [ LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiLEYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GLDYlUSgmaW5maW47RidGNEY3Rj5GK0Y+) to (0, 1/LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=] for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= > 0. In our example we want to use LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkjbW5HRiQ2JFEiMEYnRj5GPkYrRj4=, but that's only a minor annoyance. One solution is to use LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbWZyYWNHRiQ2KC1GIzYkLUkjbW5HRiQ2JFEiMUYnRj5GPi1GIzYmRistRiM2Ji1GLDYlUSJ4RidGNEY3LUY7Ni1RIitGJ0Y+RkBGQ0ZFRkdGSUZLRk0vRlBRLDAuMjIyMjIyMmVtRicvRlNGYm9GWUY+RitGPi8lLmxpbmV0aGlja25lc3NHRmZuLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmhvLyUpYmV2ZWxsZWRHRkJGPkYrRj4=, which takes [LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1JI21uR0YkNiRRIjBGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYtUSIsRidGNS8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdRJXRydWVGJy8lKXN0cmV0Y2h5R0Y+LyUqc3ltbWV0cmljR0Y+LyUobGFyZ2VvcEdGPi8lLm1vdmFibGVsaW1pdHNHRj4vJSdhY2NlbnRHRj4vJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GLDYlUSgmaW5maW47RicvJSdpdGFsaWNHRkEvRjZRJ2l0YWxpY0YnRjVGK0Y1) to (0,1]. The other, which is the one Maple uses in this case, is to break up the interval into two parts, [0,1] and [LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtSSNtbkdGJDYkUSIxRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiLEYnRjcvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHUSV0cnVlRicvJSlzdHJldGNoeUdGQC8lKnN5bW1ldHJpY0dGQC8lKGxhcmdlb3BHRkAvJS5tb3ZhYmxlbGltaXRzR0ZALyUnYWNjZW50R0ZALyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRiw2JVEoJiM4NzM0O0YnLyUnaXRhbGljR0ZDL0Y4USdpdGFsaWNGJ0Y3RitGN0YrRjc=), and use the transformation to take [LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtSSNtbkdGJDYkUSIxRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiLEYnRjcvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHUSV0cnVlRicvJSlzdHJldGNoeUdGQC8lKnN5bW1ldHJpY0dGQC8lKGxhcmdlb3BHRkAvJS5tb3ZhYmxlbGltaXRzR0ZALyUnYWNjZW50R0ZALyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRiw2JVEoJiM4NzM0O0YnLyUnaXRhbGljR0ZDL0Y4USdpdGFsaWNGJ0Y3RitGN0YrRjc=) to (0,1], leaving the other interval alone. Since we'll also have a problem at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1JI21uR0YkNiRRIjBGJ0ZARkBGK0ZARitGQA== to deal with, it's a good idea to split it up, and deal with one problem at a time. J1:= Int(1/sqrt(x^6+x), x=0..1); J2:= Int(1/sqrt(x^6+x), x=1..infinity); with(Student[Calculus1]): J2:= rhs(Rule[change,x=1/t](J2)); This one should be no problem for any of our approximate integration methods. Here's our integrand (Integrand is a command in the Student[Calculus1] package): Integrand(J2); It returns a list so it can work on expressions containing several integrals. We want the first element. integrand:= %[1]; v50:= evalf(ApproximateInt(integrand, t=1..0, method=simpson,partition=50)); ApproximateInt doesn't like integrating from 1 to 0, so we have to switch the endpoints. That introduces a minus sign. v50:= - evalf(ApproximateInt(integrand, t=0..1, method=simpson, partition=50)); v25:= - evalf(ApproximateInt(integrand, t=0..1, method=simpson, partition=25)); Here's the error estimate for v50 according to Richardson (since Simpson's Rule is of order 4) (v50 - v25)/(2^4-1); And the Richardson approximation for the integral is JR2:=(2^4*v50 - v25)/(2^4-1); Just to check, the actual error in v50 is (up to roundoff error) evalf(J2 - v50); and the actual error in JR2 is evalf(J2-JR2); Now for the other part J1. The problem here is at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JlEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIj1GJ0YyL0Y2USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQC8lKXN0cmV0Y2h5R0ZALyUqc3ltbWV0cmljR0ZALyUobGFyZ2VvcEdGQC8lLm1vdmFibGVsaW1pdHNHRkAvJSdhY2NlbnRHRkAvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZPLUkjbW5HRiQ2JVEiMEYnRjJGPEY8, where the denominator is 0. J1; What Maple seems to have done is a transformation LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNictRiw2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYkLUklbXN1cEdGJDYlLUYsNiVRInRGJ0Y2RjktSSNtbkdGJDYkUSIyRidGQC8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGQEYrRkBGK0ZARitGQA==. J3:= rhs(Rule[change, x=t^2](J1)); Notice that this has no problem at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEpJmVxdWFscztGJy9GOFEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkIvJSlzdHJldGNoeUdGQi8lKnN5bW1ldHJpY0dGQi8lKGxhcmdlb3BHRkIvJS5tb3ZhYmxlbGltaXRzR0ZCLyUnYWNjZW50R0ZCLyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUS1JI21uR0YkNiRRIjBGJ0Y+Rj5GK0Y+. How did Maple know that this transformation would work? It did a series expansion. This is essentially a Taylor series, but a little more general when necessary. series(exp(x)*cos(x), x=Pi); series(..., x=a) gives us a certain number of terms in powers of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEoJm1pbnVzO0YnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0ZRLUYsNiVRImFGJ0Y0RjdGPkYrRj4= and then LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYtUTAmQXBwbHlGdW5jdGlvbjtGJ0Y5LyUmZmVuY2VHRjgvJSpzZXBhcmF0b3JHRjgvJSlzdHJldGNoeUdGOC8lKnN5bW1ldHJpY0dGOC8lKGxhcmdlb3BHRjgvJS5tb3ZhYmxlbGltaXRzR0Y4LyUnYWNjZW50R0Y4LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGUC1JKG1mZW5jZWRHRiQ2JC1GIzYkLUklbXN1cEdGJDYlLUZUNiQtRiM2Ji1GLDYlUSJ4RicvRjdRJXRydWVGJy9GOlEnaXRhbGljRictRj02LVEoJm1pbnVzO0YnRjlGQEZCRkRGRkZIRkpGTC9GT1EsMC4yMjIyMjIyZW1GJy9GUkZkby1GLDYlUSJhRidGXG9GXm9GOUY5LUYsNiVRIm5GJ0Zcb0Zeby8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGOUY5RjlGK0Y5RitGOQ== to say that's what the remainder is. We can give it a third input to control to what order the series will be computed. series(exp(x)*cos(x), x=Pi, 8); Actually the result of series(..., x=a, n) will not necessarily have LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYtUTAmQXBwbHlGdW5jdGlvbjtGJ0Y5LyUmZmVuY2VHRjgvJSpzZXBhcmF0b3JHRjgvJSlzdHJldGNoeUdGOC8lKnN5bW1ldHJpY0dGOC8lKGxhcmdlb3BHRjgvJS5tb3ZhYmxlbGltaXRzR0Y4LyUnYWNjZW50R0Y4LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGUC1JKG1mZW5jZWRHRiQ2JC1GIzYkLUklbXN1cEdGJDYlLUZUNiQtRiM2Ji1GLDYlUSJ4RicvRjdRJXRydWVGJy9GOlEnaXRhbGljRictRj02LVEoJm1pbnVzO0YnRjlGQEZCRkRGRkZIRkpGTC9GT1EsMC4yMjIyMjIyZW1GJy9GUkZkby1GLDYlUSJhRidGXG9GXm9GOUY5LUYsNiVRIm5GJ0Zcb0Zeby8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGOUY5RjlGK0Y5RitGOQ==. Basically the n controls how the series is computed, rather than what the result will be. series((sin(x)-tan(x))/x^3, x=0, 5); series((sin(x)-tan(x))/x^3, x=0, 8); The series for our integrand around LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1JI21uR0YkNiRRIjBGJ0ZARkBGK0ZARitGQA== involves non-integer powers of x (that's why I said it was a little more general than Taylor series). series(1/sqrt(x^6+x), x=0, 10); Rule[change, x = t^2](Int(g(x),x=0..1)); If LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiw2JVEieEYnRjRGN0Y+Rj5GK0Y+ involves half-integer powers of x, LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtSSVtc3VwR0YkNiUtRiw2JVEidEYnRjRGNy1JI21uR0YkNiRRIjJGJ0Y+LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0Y+Rj5GK0Y+ will involve integer powers of t. And after multiplying by t, LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbWZyYWNHRiQ2KC1JI21uR0YkNiRRIjFGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSZtc3FydEdGJDYjLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnL0YzUSdpdGFsaWNGJy8lLmxpbmV0aGlja25lc3NHRjEvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGRS8lKWJldmVsbGVkR1EmZmFsc2VGJ0Yy becomes 1. integrand:=Integrand(J3)[1]; v50:= evalf(ApproximateInt(integrand, t=0..1, method=simpson, partition=50)); v25:= evalf(ApproximateInt(integrand, t=0..1, method=simpson, partition=25)); Here's Richardson's error estimate for v50: (v50 - v25)/(2^4-1); Maybe we'd like a bit more accuracy here. v100:= evalf(ApproximateInt(integrand, t=0..1, method=simpson, partition=100)); (v100 - v50)/(2^4-1); JR1:= (2^4*v100 - v50)/(2^4-1); So our value for the integral is JR1 + JR2; And this compares to the "exact" value: evalf(J); %-%%; Another method for removing singularities (on a finite interval) is subtraction. Let's go back to J1. J1; integrand:= Integrand(J1)[1]; S:= series(integrand, x=0, 10); The idea is that after subtracting off the "bad" terms (in particular the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJm1mcmFjR0YkNigtRiM2JC1JI21uR0YkNiRRIjFGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGOi1GIzYkLUkmbXNxcnRHRiQ2Iy1GLDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvRjtRJ2l0YWxpY0YnRjovJS5saW5ldGhpY2tuZXNzR0Y5LyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRk4vJSliZXZlbGxlZEdRJmZhbHNlRidGOkYrRjo=), the rest should be OK for something like Simpson's rule. The bad terms can be integrated exactly. badpart:= convert(S, polynom); goodpart:= integrand - badpart; badint:= int(badpart,x=0..1); ApproximateInt(goodpart,x=0..1,method=simpson, partition=50); Oops. The problem here is that two terms in goodpart are undefined at x=0. goodpart:= piecewise(x=0,0, goodpart);

Another method for removing singularities (on a finite interval) is subtraction. integrand:= Integrand(J1)[1]; S:= series(integrand, x=0, 10); The idea is that after subtracting off the "bad" terms (in particular the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJm1mcmFjR0YkNigtRiM2JC1JI21uR0YkNiRRIjFGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGOi1GIzYkLUkmbXNxcnRHRiQ2Iy1GLDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvRjtRJ2l0YWxpY0YnRjovJS5saW5ldGhpY2tuZXNzR0Y5LyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRk4vJSliZXZlbGxlZEdRJmZhbHNlRidGOkYrRjo=), the rest should be OK for something like Simpson's rule. The bad terms can be integrated exactly. badpart:= convert(S, polynom); goodpart:= integrand - badpart; badint:= int(badpart,x=0..1); ApproximateInt(goodpart,x=0..1,method=simpson, partition=50); Oops. The problem here is that two terms in goodpart are undefined at x=0. goodpart:= piecewise(x=0,0, goodpart); goodint50:= evalf(ApproximateInt(goodpart,x=0..1,method=simpson, partition=50)); v50:=evalf(badint+goodint50); evalf(J1-v50);

Here's another improper integral. In this case it's not at all obvious that it even converges. It does so because of rapid oscillation. J := Int(x*cos(x^3),x=0..infinity); plot(x*cos(x^3),x=0..5); F:=int(x*cos(x^3),x); plot(F,x=0..5); In some Maple versions there's an odd "glitch" in the graph near x = 4.3. Sometimes plots can be improved by telling Maple to use more points, with the numpoints option. plot(F,x=0..5, numpoints=1000); Maple does claim to find a value for this integral, using the Gamma function: Jtrue := value(J); evalf(%); But we're interested in numerical integration. We could try a change of variables to a finite interval: rhs(Rule[change,t=1/x](J)); But that looks awful as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEnJnJhcnI7RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSSNtbkdGJDYkUSIwRidGPkY+RitGPg==. A better transformation might be one that makes the argument of the cos simpler. rhs(Rule[change,t=x^3](J)); This actually is improper at both ends, even though the original J was fine at 0. To avoid the problems at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEpJmVxdWFscztGJy9GOFEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkIvJSlzdHJldGNoeUdGQi8lKnN5bW1ldHJpY0dGQi8lKGxhcmdlb3BHRkIvJS5tb3ZhYmxlbGltaXRzR0ZCLyUnYWNjZW50R0ZCLyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUS1JI21uR0YkNiRRIjBGJ0Y+Rj5GK0Y+, split J into two pieces. J1 := Int(x*cos(x^3), x=0..1); J2 := Int(x*cos(x^3), x=1..infinity); J1 is no problem. V1 := evalf(ApproximateInt(x*cos(x^3),x=0..1,method=simpson,partition=50)); Now use the substitution LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Jy1GLDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUYjNiQtSSVtc3VwR0YkNiUtRiw2JVEieEYnRjRGNy1JI21uR0YkNiRRIjNGJ0Y+LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0Y+RitGPkYrRj4= on J2. J3 := rhs(Rule[change,t=x^3](J2)); The next trick to use is integration by parts. J4:= Rule[parts,1/t^(1/3)/3,sin(t)](J3); Maple should know that limit. What's going on? When you want to see what's below the surface in a Maple expression, without the fancy typesetting, the lprint command is useful. lprint(%); Ah! it's using an inert Limit rather than a limit that would actually find the limit. J4 := eval(rhs(J4),Limit=limit); Notice that this integration by parts gave us a factor LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEidEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW9HRiQ2LVEqJnVtaW51czA7RicvRjZRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRlEtSSZtZnJhY0dGJDYoLUYjNiQtSSNtbkdGJDYkUSI0RidGPkY+LUYjNiQtRlo2JFEiM0YnRj5GPi8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGYW8vJSliZXZlbGxlZEdGQkY+LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0Y+ instead of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEidEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW9HRiQ2LVEqJnVtaW51czA7RicvRjZRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRlEtSSZtZnJhY0dGJDYoLUYjNiQtSSNtbkdGJDYkUSIxRidGPkY+LUYjNiQtRlo2JFEiM0YnRj5GPi8lLmxpbmV0aGlja25lc3NHRmZuLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmBvLyUpYmV2ZWxsZWRHRkJGPi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGPg==, so the integrand goes to 0 much faster as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEnJnJhcnI7RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtRiw2JVEoJmluZmluO0YnRjRGN0Y+RitGPg==. At this point we can see that the integral converges by a comparison test: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkobWZlbmNlZEdGJDYmLUYjNiQtSSZtZnJhY0dGJDYoLUYjNiUtSSNtaUdGJDYlUSRzaW5GJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUYsNiQtRiM2JC1GNjYlUSJ0RicvRjpRJXRydWVGJy9GPVEnaXRhbGljRidGPEY8RjwtRiM2Jy1JI21uR0YkNiRRIjlGJ0Y8LUkjbW9HRiQ2LVEifkYnRjwvJSZmZW5jZUdGOy8lKnNlcGFyYXRvckdGOy8lKXN0cmV0Y2h5R0Y7LyUqc3ltbWV0cmljR0Y7LyUobGFyZ2VvcEdGOy8lLm1vdmFibGVsaW1pdHNHRjsvJSdhY2NlbnRHRjsvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0Zeby1JJW1zdXBHRiQ2JUZDLUYjNiUtRjE2KC1GTTYkUSI0RidGPC1GIzYlLUZNNiRRIjNGJ0Y8RkZGSC8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGZXAvJSliZXZlbGxlZEdGO0ZGRkgvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRkZGSEZgcEZjcEZmcEZocEY8RjwvJSVvcGVuR1EifGdyRicvJSZjbG9zZUdGX3FGUC1GUTYtUSYmbGVxO0YnRjxGVEZWRlhGWkZmbkZobkZqbi9GXW9RLDAuMjc3Nzc3OGVtRicvRmBvRmZxRlAtRjE2KC1GTTYkRmJwRjxGSkZgcEZjcEZmcEZocEY8, and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYtLUkobXN1YnN1cEdGJDYnLUknbXN0eWxlR0YkNiUtSSNtb0dGJDYvUSYmaW50O0YnLyUrZm9yZWdyb3VuZEdRLlsxNDQsMTQ0LDE0NF1GJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvSSttc2VtYW50aWNzR0YkUSZpbmVydEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZALyUpc3RyZXRjaHlHUSV0cnVlRicvJSpzeW1tZXRyaWNHRkAvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGQC8lJ2FjY2VudEdGQC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlBGNUY4LUkjbW5HRiQ2JFEiMUYnRjgtSSNtaUdGJDYlUSgmaW5maW47RicvJSdpdGFsaWNHRkBGOC8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRicvJS9zdWJzY3JpcHRzaGlmdEdGaW4tSSZtZnJhY0dGJDYoLUYjNidGUy9GNlEqWzAsMCwyNTVdRicvJSlyZWFkb25seUdGRS8lMGZvbnRfc3R5bGVfbmFtZUdRKjJEfk91dHB1dEYnRjgtRiM2KS1GVDYkUSI5RidGOC1GMjYtUTEmSW52aXNpYmxlVGltZXM7RidGOEY+RkEvRkRGQEZGL0ZJRkBGSkZMRk5GUS1JJW1zdXBHRiQ2JS1GWDYlUSJ0RicvRmZuRkUvRjlRJ2l0YWxpY0YnLUZdbzYoLUZUNiRRIjRGJ0Y4LUZUNiRRIjNGJ0Y4LyUubGluZXRoaWNrbmVzc0dGVi8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0ZncS8lKWJldmVsbGVkR0ZFRmduRmFvRmNvRmVvRjhGY3FGZXFGaHEvRltyRkAtSSdtc3BhY2VHRiQ2Ji8lJ2hlaWdodEdRJjAuMGV4RicvJSZ3aWR0aEdRJjAuM2VtRicvJSZkZXB0aEdGYnIvJSpsaW5lYnJlYWtHUSVhdXRvRictRi82JS1GMjYvUTAmRGlmZmVyZW50aWFsRDtGJ0Y1RjhGO0Y+RkFGYHBGRkZhcEZKRkxGTkZRRjVGOEZlcC1GMjYtUSJ+RidGOEY+RkFGYHBGRkZhcEZKRkxGTkZRLUYyNi1RIjxGJ0Y4Rj5GQUZgcEZGRmFwRkpGTC9GT1EsMC4yNzc3Nzc4ZW1GJy9GUkZnc0Zhb0Zjb0Zlb0Y4LUkjbWlHNiMvSSttb2R1bGVuYW1lRzYiSSxUeXBlc2V0dGluZ0dJKF9zeXNsaWJHRic2JVEoJmluZmluO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRic=. We'll try this a few more times to make it converge faster. Actually, the Calculus1 package's integration by parts is a bit too picky to be convenient (it wants the exact specification of the u and v, besides having this Limit vs limit problem). There is an IntegrationTools package that has a Parts command that is easier to use, where you need only specify the expression u to be differentiated (it chooses the antiderivative of dv to use for v; that may or may not be the best choice, but you can override the choice if you want). It also has some other useful commands for manipulating integrals. We could have done the steps so far with this instead of Student[Calculus1]. with(IntegrationTools); S:=Split(J,1); J1:= op(1,S); J2:= op(2,S); J3:= Change(J2,x=t^(1/3)); J4:= Parts(J3,t^(-1/3)); And now for another integration by parts. J5:= Parts(J4,t^(-4/3)); Can we change variables to a proper integral now? Change(J5,t=1/u); The integrand here is continuous as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEidUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RKCYjODU5NDtGJy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGUy1JI21uR0YkNiRRIjBGJ0ZARkBGK0ZARitGQA== but its first derivative isn't: diff(cos(1/u)*u^(1/3),u); So we need to keep going a few more times. for k from 0 to 7 do J5 := Parts(J5,t^(-k-7/3)) end do; J6 := Change(J5,t=1/u);

Limit numpoints option for plot IntegrationTools package Split (in IntegrationTools) Change (in IntegrationTools) Parts (in IntegrationTools)

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=