# WORKSHEET #29
#
# ASYMPTOTIC SERIES
#
#
> restart;
--------------------------------------------------------------------------------
# In this worksheet we consider the problem of computing values of the Sine integral,
# Si(x)=Int(sin(t)/t,t=0..x).
--------------------------------------------------------------------------------
> Si(x)=Int(sin(t)/t,t=0..x);
x
/
| sin(t)
Si(x) = | ------ dt
| t
/
0
--------------------------------------------------------------------------------
# One possible idea is to substitute the Taylor series of sin(t), about t=0, into the integral
# and then integrate term-by-term.
--------------------------------------------------------------------------------
> sin(t)=Sum(((1/(2*k+1)!)*(-1)^k)*(t^(2*k+1)),k=0..infinity);
infinity
----- k (2 k + 1)
\ (-1) t
sin(t) = ) ----------------
/ (2 k + 1)!
-----
k = 0
--------------------------------------------------------------------------------
> sin(t)/t=Sum(((1/(2*k+1)!)*(-1)^k)*(t^(2*k)),k=0..infinity);
infinity
----- k (2 k)
sin(t) \ (-1) t
------ = ) ------------
t / (2 k + 1)!
-----
k = 0
--------------------------------------------------------------------------------
> Int(sin(t)/t,t=0..x)=Sum(((-1)^k)/(2*k+1)!*int(t^(2*k),t=0..x),k=0..infinity);
/ (2 k + 1) (2 k + 1)\
x infinity k | t x |
/ ----- (-1) |limit ---------- + ----------|
| sin(t) \ \t -> 0+ 2 k + 1 2 k + 1 /
| ------ dt = ) ---------------------------------------
| t / (2 k + 1)!
/ -----
0 k = 0
--------------------------------------------------------------------------------
# The limit term is 0, and so
--------------------------------------------------------------------------------
> Si(x)=Sum((-1)^k/((2*k+1)*((2*k+1)!))*x^(2*k+1),k=0..infinity);
infinity
----- k (2 k + 1)
\ (-1) x
Si(x) = ) --------------------
/ (2 k + 1) (2 k + 1)!
-----
k = 0
--------------------------------------------------------------------------------
# This series converges by the alternating series test since the terms do alternate in sign and
# eventually their absolute values decrease monotonically to zero. For small values of x the
# convergence is quite rapid, but it's another matter for large values of x. We will set
# A(k,x) equal to the k^(th) term of the series and then consider the problem of computing
# Si(10) and Si(100) from this alternating series.
--------------------------------------------------------------------------------
> A:=(k,x)->(-1)^k/((2*k+1)*((2*k+1)!))*x^(2*k+1);
k (2 k + 1)
(-1) x
A := (k,x) -> --------------------
(2 k + 1) (2 k + 1)!
--------------------------------------------------------------------------------
# The exact values of the first 21 terms of the series for Si(10) are:
--------------------------------------------------------------------------------
> for k from 0 to 20 do print(k,A(k,10)) od;k:='k':
0, 10
1, -500/9
2, 500/3
125000
3, - ------
441
1562500
4, -------
5103
15625000
5, - --------
68607
390625000
6, ---------
3162159
781250000
7, - ---------
15324309
24414062500
8, -----------
1476241767
1220703125000
9, - -------------
282135852999
6103515625000
10, -------------
6548521640661
305175781250000
11, - ----------------
1814564163190779
6103515625000
12, ---------------
236682282155319
7629394531250000
13, - -------------------
2243037987985958163
190734863281250000
14, ---------------------
489065356861975396503
1907348632812500000
15, - -----------------------
48619842201140519590281
5960464477539062500
16, -------------------------
1707968005065871801090839
11920928955078125000
17, - --------------------------
43113252976359733645717239
1490116119384765625000
18, -----------------------------
75885484274532611178411728817
74505805969238281250000
19, - --------------------------------
59270665136478862984997852731959
186264514923095703125000
20, ----------------------------------
2554717643446691504558497190831361
--------------------------------------------------------------------------------
# This is not too helpful, so we change the above loop by replacing 10 by 10.0. This will
# replace the exact value of A(k,10) by a 10 place approximation.
--------------------------------------------------------------------------------
> for k from 0 to 20 do print(k,A(k,10.0)) od;k:='k':
0, 10.0
1, -55.55555556
2, 166.6666667
3, -283.4467120
4, 306.1924358
5, -227.7464399
6, 123.5311064
7, -50.98109155
8, 16.53798385
9, -4.326650130
10, .9320448125
11, -.1681813118
12, .02578780114
13, -.003401366616
14, .0003899987202
15, -.00003922984005
-5
16, .3489798673*10
-6
17, -.2765026560*10
-7
18, .1963637886*10
-8
19, -.1257043527*10
-10
20, .7291002017*10
--------------------------------------------------------------------------------
# Notice that the absolute values of the terms of the series are initially increasing, but after
# term 4 they begin to decrease, and in fact tend to zero. By the alternating series test the
# sum of the exact values of the terms A(k,10), k=0..20, should approximate Si(10) with an
# error no more than the absolute value of the next term in the sequence, that is
# abs(A(21,10)).
--------------------------------------------------------------------------------
> A(21,10.0);
-11
-.3849327599*10
--------------------------------------------------------------------------------
# Maple knows the exact answer to Si(10). Let's use the exact answer to see if the
# predicted error is reasonable. A speedy way to add up the terms of a series is to use the
# convert command.
--------------------------------------------------------------------------------
> Si(10.0);
1.658347594
--------------------------------------------------------------------------------
> ?convert
--------------------------------------------------------------------------------
# The exact value of Si(10) to 10 places, possibly rounded off in the last place, is
# 1.658347594. The exact value of the sum of the first 21 terms in the series is:
--------------------------------------------------------------------------------
> convert([seq(A(k,10),k=0..20)],`+`);k:='k':
9733366047998806860644040501494565352596170
-------------------------------------------
5869315987738947410279534073197070730772067
--------------------------------------------------------------------------------
> evalf("");
1.658347594
--------------------------------------------------------------------------------
# This is in perfect agreement with the alternating series test. The actual sum of the 10
# place decimal approximations to the first 21 terms is given by:
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,10.0),k=0..20)],`+`);k:='k':
s := 1.658347506
--------------------------------------------------------------------------------
> s-Si(10.0);
-7
-.88*10
--------------------------------------------------------------------------------
# This is the actual error made when we add up the 10 place decimal approximations of the
# first 21 terms, an error that is not too close to the predicted error -.3849327599*10^(-11).
# To see why there is such a large discrepancy consider the terms of the series and the
# partial sums out to k=20. The first term is A(0,10)=10.
--------------------------------------------------------------------------------
> S[0]:=A(0,10);
S[0] := 10
--------------------------------------------------------------------------------
> for k from 1 to 10 do S[k]:=S[k-1]+A(k,10.0) od:k:='k':
--------------------------------------------------------------------------------
> for k from 1 to 10 do print(k,A(k,10.0),S[k]) od;k:='k':\
1, -55.55555556, -45.55555556
2, 166.6666667, 121.1111111
3, -283.4467120, -162.3356009
4, 306.1924358, 143.8568349
5, -227.7464399, -83.8896050
6, 123.5311064, 39.6415014
7, -50.98109155, -11.33959015
8, 16.53798385, 5.19839370
9, -4.326650130, .871743570
10, .9320448125, 1.803788383
--------------------------------------------------------------------------------
# Notice that 5 of the terms, namely those for k=2,3,4,5, and 6, all have absolute values
# more than 100, and only 7 places after the decimal point. The largest absolute value is
# about 306. We can not expect to add these terms up and have an answer which is
# accurate to 9 places after the decimal point. The actual error might be more like
# 5*10^(-7), and in fact it is. One solution is to compute to more decimal places to begin
# with. Calculations to 12 places might be sufficient to give us 10 significant places in our
# answer because 12 places allows for 3 before the decimal point and 9 after in the 5 large
# terms in the sum. The exact value of Si(10) has 1 place before the decimal point and we
# are trying to get 9 places after.
--------------------------------------------------------------------------------
> Digits:=12:
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,10.0),k=0..20)],`+`);k:='k':
s := 1.65834759372
--------------------------------------------------------------------------------
> Digits:=10:
--------------------------------------------------------------------------------
> evalf(""");
1.658347594
--------------------------------------------------------------------------------
> s-Si(10.0);
0
--------------------------------------------------------------------------------
# Rounding off to 10 significant places thus produces the correct value. If we want to
# compute Si(10) to within the truncation error 0.3849327599*10^(-11) then we should be
# carrying 12 places after the decimal point and therefore do our calculations to 15
# significant figures, 3 before the decimal point and 12 after. First we try 14 digits and then
# 15.
--------------------------------------------------------------------------------
> Digits:=14:
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,10.0),k=0..20)],`+`);k:='k':
s := 1.6583475942293
--------------------------------------------------------------------------------
> s-Si(10.0);
-10
.104*10
--------------------------------------------------------------------------------
# Not quite enough, just as we suspected.
--------------------------------------------------------------------------------
> Digits:=15:
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,10.0),k=0..20)],`+`);k:='k':
s := 1.65834759422126
--------------------------------------------------------------------------------
> s-Si(10.0);
-11
.239*10
--------------------------------------------------------------------------------
# This agrees perfectly with our predicted truncation error.
--------------------------------------------------------------------------------
> Digits:=10:
--------------------------------------------------------------------------------
# If we are to use the alternating series for Si(x) as a computational tool we should first
# determine how large k must be in order that the terms abs(A(k,x)) are decreasing. this
# we do by comparing A(k+1,x) with A(k,x).
--------------------------------------------------------------------------------
> abs(A(k+1,x)/A(k,x));
(k + 1) (2 k + 3)
(-1) x (2 k + 1) (2 k + 1)!
abs(-------------------------------------------)
k (2 k + 1)
(2 k + 3) (2 k + 3)! (-1) x
--------------------------------------------------------------------------------
> simplify(");
2
x (2 k + 1)
1/2 abs(------------------)
2
(2 k + 3) (k + 1)
--------------------------------------------------------------------------------
# Thus the terms are decreasing if 2*k+3>x. Now we try x=100. The maximum absolute
# value should occur for k=49.
--------------------------------------------------------------------------------
> for k from 46 to 52 do print(k,A(k,100.0)) od;k:='k':
40
46, .9295421620*10
41
47, -.1019006555*10
41
48, .1071731218*10
41
49, -.1082333624*10
41
50, .1050397302*10
40
51, -.9803933355*10
40
52, .8806952167*10
--------------------------------------------------------------------------------
# From k=49 onward the terms are decreasing. The largest absolute value is in fact
# abs(A(49,100))=0.1082333624*10^(41), a very large term. It is not surprising that adding
# up an alternating series with such large terms is not accurate, unless one does the
# calculations to a great many places. Moreover it is easy to see that the maximum value of
# abs(Si(x)) occurs at x=Pi, and that this value is Si(Pi)=1.851937052.
--------------------------------------------------------------------------------
> plot(sin(t)/t,t=0..20*Pi);
--------------------------------------------------------------------------------
#
** Maple V Graphics **
--------------------------------------------------------------------------------
> evalf(Si(Pi));
1.851937052
--------------------------------------------------------------------------------
# The individual terms in the sum for Si(100) have a maximum absolute value of about
# 10^(40), so if we want to approximate Si(100) to 10 places with this sum we should do
# our calculations to at least 50 places, and perhaps more. Next we estimate how many
# terms are necessary to make the truncation error less than 0.5*10^(-10). We do it by
# experiment.
--------------------------------------------------------------------------------
> A(100,100.0);
23
.3138479130*10
--------------------------------------------------------------------------------
> A(200,100.0);
-72
.9711740670*10
--------------------------------------------------------------------------------
> A(150,100.0);
-17
.3606317547*10
--------------------------------------------------------------------------------
> A(140,100.0);
-8
.7550840789*10
--------------------------------------------------------------------------------
# Thus the vakue of k we are looking for is somewhere between k=141 and k=149.
--------------------------------------------------------------------------------
> for k from 141 to 149 do print(k,A(k,100.0)) od;k:='k':
-9
141, -.9394629422*10
-9
142, .1152545367*10
-10
143, -.1394354050*10
-11
144, .1663668471*10
-12
145, -.1957855579*10
-13
146, .2272770322*10
-14
147, -.2602746120*10
-15
148, .2940689792*10
-16
149, -.3278287386*10
--------------------------------------------------------------------------------
# Thus k=143 should suffice.
--------------------------------------------------------------------------------
> Digits:=50:
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,100.0),k=0..142)],`+`);k:='k':
s := 1.5622254671018448150728932832037677330244317383092
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;
elapsedtime := 2.567
--------------------------------------------------------------------------------
# The convert command is indeed quite fast.
--------------------------------------------------------------------------------
> evalf(Si(100));
1.5622254668890562933523451388045026772278249805411
--------------------------------------------------------------------------------
> s-evalf(Si(100));
-9
.2127885217205481443992650557966067577681*10
--------------------------------------------------------------------------------
# This is not quite accurate to 10 places. The problem must be that there are too many
# terms in the series with absolute value near 10^(40).
--------------------------------------------------------------------------------
> Digits:=51:
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,100.0),k=0..142)],`+`);k:='k':
s := 1.56222546696625795842933532621285898862481574590685
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;
elapsedtime := 2.467
--------------------------------------------------------------------------------
> s-evalf(Si(100));
-10
.7720166507699018740835631139699076536577*10
--------------------------------------------------------------------------------
# Still not there.
--------------------------------------------------------------------------------
> Digits:=52:
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,100.0),k=0..142)],`+`);k:='k':
s := 1.562225466904483402648495792398277959940529764588085
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;
elapsedtime := 3.184
--------------------------------------------------------------------------------
> s-evalf(Si(100));
-10
.15427109296150653593775282712704784047002*10
#
--------------------------------------------------------------------------------
> Digits:=53:
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,100.0),k=0..142)],`+`);k:='k':
s := 1.5622254669021758132370188522007506652086266569193657
--------------------------------------------------------------------------------
> s-evalf(Si(100));
-10
.131195198846737133962479879808016763782822*10
--------------------------------------------------------------------------------
# This is now in quite good agreement with the truncation error A(143,100).
--------------------------------------------------------------------------------
> A(143,100.0);
-10
-.13943540504088827334043593535271574839854399819848964*10
--------------------------------------------------------------------------------
> Digits:=10:
--------------------------------------------------------------------------------
# If we do not carry 53 significant digits then the actual error can be very far removed from
# the predicted error, that is from the truncation error. Consider the next calculation. Here
# we are using only 10 places for the terms in the series.
--------------------------------------------------------------------------------
> s:=convert([seq(A(k,100.0),k=0..142)],`+`);k:='k':
31
s := -.3538632006*10
--------------------------------------------------------------------------------
# What a ridiculous answer!
--------------------------------------------------------------------------------
# It is quite clear that alternating series have limited use if the the series has large terms
# and the limit is relatively small. Fortunately there are other ways of estimating the Sine
# integral.
--------------------------------------------------------------------------------
> Si(x)=Int(sin(t)/t,t=0..infinity)-Int(sin(t)/t,t=x..infinity);
infinity infinity
/ /
| sin(t) | sin(t)
Si(x) = | ------ dt - | ------ dt
| t | t
/ /
0 x
--------------------------------------------------------------------------------
# The exact value of Si(infinity)=Int(sin(t)/t,t=0..infinity) is Pi/2, so the problem of
# computing Si(x) becomes the problem of computing the complementary Sine integral
# Int(sin(t)/t,t=x..infinity).
--------------------------------------------------------------------------------
> Si(infinity);
1/2 Pi
--------------------------------------------------------------------------------
> Sic:=x->Int(sin(t)/t,t=x..infinity);
infinity
/
| sin(t)
Sic := x -> | ------ dt
| t
/
x
--------------------------------------------------------------------------------
# We can substitute the Taylor series for sin(t)/t into the integral, but we can no longer
# integrate term-by-term. We need another idea.
--------------------------------------------------------------------------------
> with(student):
--------------------------------------------------------------------------------
> B[1]:=Sic(x);
infinity
/
| sin(t)
B[1] := | ------ dt
| t
/
x
--------------------------------------------------------------------------------
# In the next calculation we integrate repeatedly by parts.
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> for k from 2 to 150 do B[k]:=intparts(B[k-1],1/t^(k-1)) od:k:='k':
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;
elapsedtime := 16.150
--------------------------------------------------------------------------------
> B[10];
cos(x) sin(x) cos(x) sin(x) cos(x) sin(x) cos(x)
------ + ------ - 2 ------ - 6 ------ + 24 ------ + 120 ------ - 720 ------
x 2 3 4 5 6 7
x x x x x x
infinity
/
sin(x) cos(x) | cos(t)
- 5040 ------ + 40320 ------ - | 362880 ------ dt
8 9 | 10
x x / t
x
--------------------------------------------------------------------------------
# Notice that the first 9 terms in B[10] can easily be computed for any value of x, but that
# the tenth term, the integral that is, can only be estimated. Thus the idea is to use the
# non-integral terms to approximate Sic(x) and then estimate the integral to determine the
# error. This means we have to be able to select the non-integral terms from the
# expressions B[k], and also the integral term. This is done by the operand command.
--------------------------------------------------------------------------------
> ?op
--------------------------------------------------------------------------------
# Looking at B[10] we see that the number of operands is 10 and that the tenth operand is
# the integral term.
--------------------------------------------------------------------------------
> nops(B[10]);
10
--------------------------------------------------------------------------------
> op(10,B[10]);
infinity
/
| cos(t)
- | 362880 ------ dt
| 10
/ t
x
--------------------------------------------------------------------------------
> B[10]-op(10,B[10]);
cos(x) sin(x) cos(x) sin(x) cos(x) sin(x) cos(x)
------ + ------ - 2 ------ - 6 ------ + 24 ------ + 120 ------ - 720 ------
x 2 3 4 5 6 7
x x x x x x
sin(x) cos(x)
- 5040 ------ + 40320 ------
8 9
x x
--------------------------------------------------------------------------------
# This last expression is our approximation to Sic(x) and the expression before this one is
# the error. Unfortunately Maple does not always give the terms of integration by parts in
# the same order so we have to sort them first. For example consider B[25].
--------------------------------------------------------------------------------
> B[25];
sin(x) cos(x) cos(x) sin(x) sin(x)
120 ------ - 720 ------ + 40320 ------ - 5040 ------ + 362880 ------
6 7 9 8 10
x x x x x
cos(x) sin(x) cos(x)
- 3628800 ------ - 39916800 ------ + 479001600 ------
11 12 13
x x x
sin(x) cos(x) sin(x)
+ 51090942171709440000 ------ - 87178291200 ------ + 6227020800 ------
22 15 14
x x x
sin(x) sin(x) cos(x)
- 1307674368000 ------ + 355687428096000 ------ + 20922789888000 ------
16 18 17
x x x
sin(x) cos(x) cos(x) sin(x)
+ ------ - 2 ------ - 6402373705728000 ------ - 121645100408832000 ------
2 3 19 20
x x x x
sin(x) cos(x) cos(x)
- 6 ------ + 2432902008176640000 ------ - 1124000727777607680000 ------
4 21 23
x x x
infinity
/
cos(x) | sin(t) cos(x)
+ ------ + | 620448401733239439360000 ------ dt + 24 ------
x | 25 5
/ t x
x
sin(x)
- 25852016738884976640000 ------
24
x
--------------------------------------------------------------------------------
> ?sort
--------------------------------------------------------------------------------
> sort(B[25]);
infinity
/
| sin(t) cos(x) sin(x) cos(x)
| 620448401733239439360000 ------ dt + ------ + ------ - 2 ------
| 25 x 2 3
/ t x x
x
sin(x) cos(x) sin(x) cos(x) sin(x)
- 6 ------ + 24 ------ + 120 ------ - 720 ------ - 5040 ------
4 5 6 7 8
x x x x x
cos(x) sin(x) cos(x) sin(x)
+ 40320 ------ + 362880 ------ - 3628800 ------ - 39916800 ------
9 10 11 12
x x x x
cos(x) sin(x) cos(x)
+ 479001600 ------ + 6227020800 ------ - 87178291200 ------
13 14 15
x x x
sin(x) cos(x) sin(x)
- 1307674368000 ------ + 20922789888000 ------ + 355687428096000 ------
16 17 18
x x x
cos(x) sin(x)
- 6402373705728000 ------ - 121645100408832000 ------
19 20
x x
cos(x) sin(x)
+ 2432902008176640000 ------ + 51090942171709440000 ------
21 22
x x
cos(x) sin(x)
- 1124000727777607680000 ------ - 25852016738884976640000 ------
23 24
x x
> for k from 1 to 150 do Err[k]:=op(1,sort(B[k])); C[k]:=B[k]-Err[k] od:k:='k':
--------------------------------------------------------------------------------
# Sorting the terms of B[k] puts the integral first which we can then select by the operand
# command. Err[k] is the integral term, that is the error when we approximate Sic(x), and
# C[k] is the approximation.
--------------------------------------------------------------------------------
> Err[25];
infinity
/
| sin(t)
| 620448401733239439360000 ------ dt
| 25
/ t
x
--------------------------------------------------------------------------------
> 24!;
620448401733239439360000
--------------------------------------------------------------------------------
> sort(C[25]);
cos(x) sin(x) cos(x) sin(x) cos(x) sin(x) cos(x)
------ + ------ - 2 ------ - 6 ------ + 24 ------ + 120 ------ - 720 ------
x 2 3 4 5 6 7
x x x x x x
sin(x) cos(x) sin(x) cos(x)
- 5040 ------ + 40320 ------ + 362880 ------ - 3628800 ------
8 9 10 11
x x x x
sin(x) cos(x) sin(x)
- 39916800 ------ + 479001600 ------ + 6227020800 ------
12 13 14
x x x
cos(x) sin(x) cos(x)
- 87178291200 ------ - 1307674368000 ------ + 20922789888000 ------
15 16 17
x x x
sin(x) cos(x)
+ 355687428096000 ------ - 6402373705728000 ------
18 19
x x
sin(x) cos(x)
- 121645100408832000 ------ + 2432902008176640000 ------
20 21
x x
sin(x) cos(x)
+ 51090942171709440000 ------ - 1124000727777607680000 ------
22 23
x x
sin(x)
- 25852016738884976640000 ------
24
x
--------------------------------------------------------------------------------
# The error term will have the form Err[k]=(k-1)!*Int(trig(t)/t^k,t=x..infinity), where
# trig(t) is one of the 4 functions sin(t), cos(t), -sin(t) or -cos(t). This allows us to find an
# easy bound on the error.
--------------------------------------------------------------------------------
> Err[k]=(k-1)!*Int(trig(t)/t^k,t=x..infinity);
infinity
/
| trig(t)
Err[k] = (k - 1)! | ------- dt
| k
/ t
x
--------------------------------------------------------------------------------
> abs(Err[k])<(k-1)!*int(t^(-k),t=x..infinity);
/ (- k + 1) (- k + 1)\
| t x |
abs(Err[k]) < (k - 1)! |limit ---------- + ----------|
\t -> infinity- - k + 1 k - 1 /
--------------------------------------------------------------------------------
# The limit term is zero for k>1 and therefore we have.
--------------------------------------------------------------------------------
> abs(Err[k])<(k-2)!/x^(k-1);
(k - 2)!
abs(Err[k]) < --------
(k - 1)
x
--------------------------------------------------------------------------------
> Bound:=(k,x)->(k-2)!/x^(k-1);
(k - 2)!
Bound := (k,x) -> --------
(k - 1)
x
--------------------------------------------------------------------------------
> Bound(k,x)/Bound(k-1,x);
(k - 2)
(k - 2)! x
-----------------
(k - 1)
x (k - 3)!
--------------------------------------------------------------------------------
> simplify(");
k - 2
-----
x
--------------------------------------------------------------------------------
# Keeping x fixed and letting k vary we see that the bound on the error,Bound(k,x), is the
# smallest when k=floor(x+2). For larger or smaller values of k the bound is bigger.
# Moreover, Limit(Bound(k,x),k=infinity)=infinity, and therefore there is definitely a limit
# on the amount of accuracy obtainable this way. In fact Limit(Err[k],k=infinity)=infinity
# and the approximating series, that is, Limit(C[k],k=infinity), actually diverges. We
# conclude with a few examples.
--------------------------------------------------------------------------------
> for k from 10 to 14 do print(k,Bound(k,10.0)) od;k:='k':
10, .00004032000000
11, .00003628800000
12, .00003628800000
13, .00003991680000
14, .00004790016000
--------------------------------------------------------------------------------
# Thus the bound on the error estimate is smallest for k=11 or 12. There is clearly no point
# in using more than 11 terms. Thus we predict that C[10] should be the best
# approximation.
--------------------------------------------------------------------------------
> evalf(Sic(10));
-.08755126742
--------------------------------------------------------------------------------
> for k from 2 to 20 do print(k,evalf(subs(x=10.0,C[k])),abs(evalf(subs(x=10.0,C[k]))-evalf(Sic(10)))) od;k:='k':
2, -.08390715291, .00364411451
3, -.08934736402, .00179609660
4, -.08766922096, .00011795354
5, -.08734280829, .00020845913
-5
6, -.08754418546, .708196*10
7, -.08760946800, .00005820058
-5
8, -.08754905485, .221257*10
9, -.08752163618, .00002963124
-5
10, -.08755546754, .420012*10
11, -.08757520898, .00002394156
-5
12, -.08754476076, .650666*10
13, -.08752304518, .00002822224
14, -.08756323687, .00001196945
15, -.08759711318, .00004584576
16, -.08752396435, .00002730307
17, -.08745282410, .00009844332
18, -.08762838127, .00007711385
19, -.08782188274, .00027061532
20, -.08728467777, .00026658965
--------------------------------------------------------------------------------
# From these results we see that the best approximation is actually C[8] and not C[10]. The
# reason for this is that in obtaining the bound we replaced abs(trig(t)) by 1, thereby
# making the error approximation bigger.
--------------------------------------------------------------------------------
> evalf(Sic(100));
.008570859906
--------------------------------------------------------------------------------
> for k from 92 to 112 do print(k,Bound(k,100.0)) od;k:='k':
-43
92, .1485715964*10
-43
93, .1352001528*10
-43
94, .1243841405*10
-43
95, .1156772507*10
-43
96, .1087366157*10
-43
97, .1032997849*10
-44
98, .9916779349*10
-44
99, .9619275968*10
-44
100, .9426890449*10
-44
101, .9332621544*10
-44
102, .9332621544*10
-44
103, .9425947760*10
-44
104, .9614466715*10
-44
105, .9902900716*10
-43
106, .1029901675*10
-43
107, .1081396758*10
-43
108, .1146280564*10
-43
109, .1226520203*10
-43
110, .1324641819*10
-43
111, .1443859583*10
-43
112, .1588245542*10
--------------------------------------------------------------------------------
> Digits:=45:
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> for k from 90 to 112 do print(k,evalf(subs(x=100.0,C[k])),abs(evalf(subs(x=100.0,C[k]))-evalf(Sic(100)))) od;k:='k':
-43
90, .00857085990584032587897655283524876487075973083, .1168*10
-44
91, .00857085990584032587897655283524876487075972247, .332*10
-44
92, .00857085990584032587897655283524876487075970960, .955*10
-44
93, .00857085990584032587897655283524876487075971648, .267*10
-44
94, .00857085990584032587897655283524876487075972722, .807*10
-44
95, .00857085990584032587897655283524876487075972129, .214*10
-44
96, .00857085990584032587897655283524876487075971200, .715*10
-44
97, .00857085990584032587897655283524876487075971723, .192*10
-44
98, .00857085990584032587897655283524876487075972579, .664*10
-44
99, .00857085990584032587897655283524876487075972091, .176*10
-44
100, .00857085990584032587897655283524876487075971279, .636*10
-44
101, .00857085990584032587897655283524876487075971749, .166*10
-44
102, .00857085990584032587897655283524876487075972555, .640*10
-44
103, .00857085990584032587897655283524876487075972077, .162*10
-44
104, .00857085990584032587897655283524876487075971248, .667*10
-44
105, .00857085990584032587897655283524876487075971749, .166*10
-44
106, .00857085990584032587897655283524876487075972638, .723*10
-44
107, .00857085990584032587897655283524876487075972090, .175*10
-44
108, .00857085990584032587897655283524876487075971116, .799*10
-44
109, .00857085990584032587897655283524876487075971735, .180*10
-44
110, .00857085990584032587897655283524876487075972866, .951*10
-44
111, .00857085990584032587897655283524876487075972135, .220*10
-43
112, .00857085990584032587897655283524876487075970766, .1149*10
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;
elapsedtime := 59.933
--------------------------------------------------------------------------------
# This is a much better way of computing Si(100) then our previous approach using Taylor
# series. Si(100) to high accuracy is given by:
--------------------------------------------------------------------------------
> evalf(subs(x=100.0,C[102]));
.00857085990584032587897655283524876487075972555
--------------------------------------------------------------------------------
