2.18 ASYMPTOTIC SERIES

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

```