f
(
x
)
=
1
x
x
+a
1
x
+a
2
x
+a
3
x
+a
4
x
8
+b
1
x
6
+b
2
x
4
+b
3
x
2
+b
4
+ ✒
(
x
)
✒
(
x
)
< 5 $ 10
−7
a
1
= 38.027264 b
1
= 40.021433
a
2
= 265.187033 b
2
= 322.624911
a
3
= 335.677320 b
3
= 570.236280
a
4
= 38.102495 b
4
= 157.105423
g
(
x
)
=
1
x
2
x
8
+a
1
x
6
+a
2
x
4
+a
3
x
2
+a
4
x
8
+b
1
x
6
+b
2
x
4
+b
3
x
2
+b
4
+ ✒
(
x
)
✒
(
x
)
< 3 $ 10
−7
a
1
= 42.242855 b
1
= 48.196927
a
2
= 302.757865 b
2
= 482.485984
a
3
= 352.018498 b
3
= 1114.978885
a
4
= 21.821899 b
4
= 449.690326
All of these equations are from Handbook of Mathematical Functions, Abramowitz and Stegun, Dover,
1965.
The results from both of these routines are only accurate to about seven significant digits. Below ,
x = ✜
the accuracy may be as good as 11 or 12 significant digits. Notice that the maximum argument for both
routines is 1E12 radians. This limit results from the use of the built-in sin() and cos() functions.
The evaluation of ci(x) and si(x) is also complicated by the fact that nInt() evaluates the integrals very
slowly, and with poor accuracy, for arguments less than about 1E-4. To get around this, I used Taylor
series expansions for the integrands near x=0, then symbolically integrated those. For the cosine
integral, I use
cos(t)−1
t
=−
t
2
+
t
3
24
−
t
5
720
+
t
7
40320
+ ...
and integrated this expansion to find
Ci(x) j −
x
2
4
+
x
4
96
−
x
6
4320
+
x
8
322560
Similarly, for the sine integral:
sin(t)
t
= 1 −
t
2
6
+
t
4
120
−
t
6
5040
+
t
8
362880
−
t
10
39916800
+ ...
and integrating gives
Si(x) j x −
x
3
18
+
x
5
600
−
x
7
35280
+
x
9
3265920
−
x
11
439084800
These approximations give good accuracy for x < 0.01.
Code listing for si(x):
si(x)
Func
©(x) return si(x), x real, |x|<112
©Must be in folder \spfn
©24jan01/dburkett@infinet.com
6 - 63