The sine integral. (This is featured in problem 92 on page 400.)f:=x->int(sin(t)/t,t=0..x);plot(f(x),x=0..8*Pi);f(2);evalf(Si(2),20);Let's say we want to evaluate f(2) to many places. We could just call upon Maple to do the evaluation, and treat it as a black box. evalf(Si(2),100);Maple knows about the sine integral, and uses the name Si for the function. We could try our own hand at it. Now Riemann sums tend to be slow going. They work. They're reliable. But they're slow. How would that go? A Riemann sum for an integral LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkobXN1YnN1cEdGJDYnLUkjbW9HRiQ2MFErJkludGVncmFsO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R1EldHJ1ZUYnLyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGPC8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSVmb3JtR1EncHJlZml4RicvJSdsc3BhY2VHUSQwZW1GJy8lJ3JzcGFjZUdGSi8lKG1pbnNpemVHUSFGJy8lKG1heHNpemVHRk8tSSNtaUdGJDYnUSJhRicvJSdpdGFsaWNHRjwvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y8L0YzUSdpdGFsaWNGJy1GUzYnUSJiRidGVkZYRmVuRmduLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy8lL3N1YnNjcmlwdHNoaWZ0R0Zeby1GUzYnUSJmRidGVi9GWVErWzAsMTYwLDgwXUYnRmVuRmduLUknbXNwYWNlR0YkNiYvJSdoZWlnaHRHUSYwLjBleEYnLyUmd2lkdGhHUSYwLjVlbUYnLyUmZGVwdGhHRltwLyUqbGluZWJyZWFrR1ElYXV0b0YnLUYvNjBRMCZEaWZmZXJlbnRpYWxEO0YnRjJGNUY4L0Y7RjdGPS9GQEY3RkFGQy9GRkZPL0ZJRk8vRkxGT0ZNRlAtRlM2J1EieEYnRlZGWEZlbkZnbi1GLzYwUSJ+RidGMi9GNlEmZmFsc2VGJy9GOUZjcS9GO0ZjcS9GPkZjcS9GQEZjcS9GQkZjcS9GREZjcUZpcEZIRksvRk5RIjFGJy9GUVEpaW5maW5pdHlGJw==sets a number n, splits the interval into n pieces of equal length, picks a point from each interval at which to evaluate the integrand, and sums the value of the integrand f at that point times the common length of the intervals. That sum is an approximation to the integral, and the error, (difference between sum and actual integral) is no worse than the greatest change in the integrand over any one of the subintervals, times |b-a|, the length of the region of integration. Here, we take 200 intervals, each of length 1/100. We'd have to come up with an estimate for the derivative of sin(x)/x in order to get a rigorous take on the difference between our Riemann sum and the actual value, but here, we neglect that, partly because we trust the graph and it is telling us that our integrand is reasonably smooth, and partly because we will be getting a much better answer shortly by altogether different methods. n:=100;riemannsum:=(1/n)*sum(evalf(sin(j/n)/(j/n)),j=1..200);We could average the Riemann sums got by starting at the left, with that got by starting at the right. Happily, to do that, we need only add half the difference between the last term and what would have been the first, 1/200, in the previous calculation. This would give us a nice bump in accuracy for almost no extra computational work. riemannsum+(1/200)*(1-sin(2)/2);evalf(%);We could use the series expansion for sin(x), the one we got in class on Wednesday by running the fact that cos(x)<=1 for all x through the logic mill based on the fact that if f\342\211\244g on the interval [0,x], then JSMlP0c=LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYuLUkobXN1YnN1cEdGJDYnLUkjbW9HRiQ2MFErJkludGVncmFsO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R1EldHJ1ZUYnLyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGPC8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSVmb3JtR1EncHJlZml4RicvJSdsc3BhY2VHUSQwZW1GJy8lJ3JzcGFjZUdGSi8lKG1pbnNpemVHUSFGJy8lKG1heHNpemVHRk8tRiM2Iy1JI21uR0YkNiRRIjBGJ0YyLUYjNiMtSSNtaUdGJDYlUSJ4RicvJSdpdGFsaWNHRjwvRjNRJ2l0YWxpY0YnLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy8lL3N1YnNjcmlwdHNoaWZ0R0Zeby1GZW42J1EiZkYnRmhuLyUrZm9yZWdyb3VuZEdRK1swLDE2MCw4MF1GJy8lLHBsYWNlaG9sZGVyR0Y8RmpuLUknbXNwYWNlR0YkNiYvJSdoZWlnaHRHUSYwLjBleEYnLyUmd2lkdGhHUSYwLjVlbUYnLyUmZGVwdGhHRl5wLyUqbGluZWJyZWFrR1ElYXV0b0YnLUYvNjBRMCZEaWZmZXJlbnRpYWxEO0YnRjJGNUY4L0Y7RjdGPS9GQEY3RkFGQy9GRkZPL0ZJRk8vRkxGT0ZNRlAtRmVuNiVRInRGJ0ZobkZqbi1GLzYwUSYmbGVxO0YnRjIvRjZRJmZhbHNlRicvRjlGZnEvRjtGZnEvRj5GZnEvRkBGZnEvRkJGZnEvRkRGZnEvRkZRJmluZml4RicvRklRL3RoaWNrbWF0aHNwYWNlRicvRkxGYHIvRk5RIjFGJy9GUVEpaW5maW5pdHlGJ0YrLUZlbjYlUSJnRidGaG5Gam4tRi82MFEifkYnRjJGZXFGZ3FGaHFGaXFGanFGW3JGXHJGXHFGSEZLRmJyRmRyRmdwRl9xLUYvNjBRIi5GJ0YyRmVxRmdxRmhxRmlxRmpxRltyRlxyRl1yRkhGS0ZickZkcg==sin(x)=x-x^3/3!+x^5/5!-x^7/7!...=sum((-1)^k*x^(2*k+1)/(2*k+1)!,k=0..infinity);Of course, we cannot actually do infinitely many terms. We shall have to quit after some number of terms. Now as we showed in class, the actual value rests between the value got by quitting after n terms, and that got by quitting instead after n+1 terms. When the next term, the first one we don't use, is small enough, we know we have an answer that's good enough. The same goes for the sine integral, since the actually integrand lies always between the integrand we'd get by using n terms, and the one we'd get by using n+1. So, here goes, with twenty terms. We use the fact that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYxLUkobXN1YnN1cEdGJDYnLUkjbW9HRiQ2MFErJkludGVncmFsO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R1EldHJ1ZUYnLyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGPC8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSVmb3JtR1EncHJlZml4RicvJSdsc3BhY2VHUSQwZW1GJy8lJ3JzcGFjZUdGSi8lKG1pbnNpemVHUSFGJy8lKG1heHNpemVHRk8tRiM2Iy1JI21uR0YkNiRRIjBGJ0YyLUYjNiMtSSNtaUdGJDYlUSJ4RicvJSdpdGFsaWNHRjwvRjNRJ2l0YWxpY0YnLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy8lL3N1YnNjcmlwdHNoaWZ0R0Zeby1JKG1mZW5jZWRHRiQ2JC1GIzYlLUZVNiRRIjFGJ0YyLUYvNjBRIi9GJ0YyL0Y2USZmYWxzZUYnL0Y5Rl1wRjovRj5GXXAvRkBGXXAvRkJGXXAvRkRGXXAvRkZRJmluZml4RicvRklRLnRoaW5tYXRoc3BhY2VGJy9GTEZmcC9GTkZoby9GUVEpaW5maW5pdHlGJy1GZW42JVEidEYnRmhuRmpuRjItRi82MFEnJnNkb3Q7RidGMkZccEZecC9GO0ZdcEZfcEZgcEZhcEZicEZjcEZIRktGaHBGaXBGW3EtRi82MFEiXkYnRjJGXHBGXnBGYXFGX3BGYHBGYXBGYnBGY3AvRklRMnZlcnl0aGlubWF0aHNwYWNlRicvRkxGZnFGaHBGaXAtRmJvNiQtRiM2Ji1GVTYkUSIyRidGMi1GZW42JVEia0YnRmhuRmpuLUYvNjBRIitGJ0YyRlxwRl5wRmFxRl9wRmBwRmFwRmJwRmNwL0ZJUTBtZWRpdW1tYXRoc3BhY2VGJy9GTEZmckZocEZpcEZmb0YyRmlvRmhxLUYvNjBRIiFGJ0YyRlxwRl5wRmFxRl9wRmBwRmFwRmJwRmNwRmVxRmdxRmhwRmlwLUYvNjBRIn5GJ0YyRlxwRl5wRmFxRl9wRmBwRmFwRmJwL0ZGRk9GSEZLRmhwRmlwLUYvNjBRMCZEaWZmZXJlbnRpYWxEO0YnRjJGNUY4L0Y7RjdGPS9GQEY3RkFGQ0Zecy9GSUZPL0ZMRk9GTUZQRltxLUYvNjBRIj1GJ0YyRlxwRl5wRmFxRl9wRmBwRmFwRmJwRmNwL0ZJUS90aGlja21hdGhzcGFjZUYnL0ZMRmpzRmhwRmlwLUkmbWZyYWNHRiQ2KC1GIzYjLUklbXN1cEdGJDYlRlotRiM2I0ZocUZcby1GIzYjLUZibzYkLUYjNiZGaHFGaHJGXnFGaHFGMi8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGYXUvJSliZXZlbGxlZEdGXXAtRi82MFEiLkYnRjJGXHBGXnBGYXFGX3BGYHBGYXBGYnBGY3BGSEZLRmhwRmlwHere, x=2 because we're evaluating Si(2). a:=2;for k from 1 to 20 do a:=a+(-1)^k*2^(2*k+1)/((2*k+1)!*(2*k+1)) od:a;evalf(a,40);The next term, the one we didn't use, is k=21st, and would have given -2^(43)/(43!*43)evalf(-2^(43)/(43!*43));That's an error of about 3 times 10^(-42). This shows that taking the first 20 terms was good enough to get us 40 places. With today's machinery, it is no problem at all to take 20 terms, or even 200. a:=2;for k from 1 to 200 do a:=a+(-1)^k*2^(2*k+1)/((2*k+1)!*(2*k+1)) od:evalf(a,400);evalf(Si(2),400);We didn't really need to know the value of this integral to 4 places, much less 400. On the other hand, this tends to drive home the point that these series methods are really powerful, and permit us to say, in many cases, that our approximate answers are so good that they may as well be exact.