# 2.16 TABULATING FUNCTIONS NUMERICALLY

```# WORKSHEET # 27
#
# TABULATING FUNCTIONS NUMERICALLY
#
#
--------------------------------------------------------------------------------
>  restart;
--------------------------------------------------------------------------------
# In this worksheet we will estimate the Sine integral, Si(x)=Int(sin(t)/t,t=0..x), for a
# sequence of values of x.  This is a convenient test case for the general procedure of
# building tables for functions since Maple already knows the values of the Sine integral. In
# other words we can tabulate the values of Si(x) using some numerical procedure and then
# compare with the known values.  But first we use the help command to get some
# information on the Sine integral.  I have edited the help page and only copied what is
# relevant.
--------------------------------------------------------------------------------
> ?Si
--------------------------------------------------------------------------------
#
# FUNCTION: Si - The Sine Integral
#
# CALLING SEQUENCE:
#    Si(x)
#
# PARAMETERS:
#    x - an expression
#
# SYNOPSIS:
# - This integral is defined as follows:
#
#       Si(x)   = int(sin(t)/t, t=0..x)
#
# EXAMPLES:
#
# > Si(3.14159+7.6*I);
#                           63.60695388 - 123.1816272 I
#
#
# > Si(12345.67890);
#                                   1.570738760
#
#
--------------------------------------------------------------------------------
> Si(x)=Int(sin(t)/t,t=0..x);

x
/
|  sin(t)
Si(x) =  |  ------ dt
|     t
/
0
--------------------------------------------------------------------------------
# Now suppose we want to tabulate the values of Si(x) for x=0.0, 0.1, 0.2,...,10.0, that is the
# values of Si(x[n]), where x[n]=n*(0.1), n=0..100.  Each of the integrals
# Int(f(t),t=x[n]..x[n+1]) will be estimated by Simpson's rule S[2] and then added up to give
# an estimate of Int(f(t),t=0..x) for x=0.0, 0.1, 0.2,...,10.0.  Thus the step size h=0.05.
--------------------------------------------------------------------------------
> h:=n->(b-a)/n;

b - a
h := n -> -----
n
--------------------------------------------------------------------------------
>  X:=(i,n)->a+i*(b-a)/n;

i (b - a)
X := (i,n) -> a + ---------
n
--------------------------------------------------------------------------------
> Simpson:=n->(h(n)/3)*sum(f(X(2*i,n))+4*f(X(2*i+1,n))+f(X(2*i+2,n)),i=0..n/2-1);

Simpson :=

/1/2 n - 1                                                       \
|  -----                                                         |
|   \                                                            |
n -> 1/3 h(n) |    )     (f(X(2 i, n)) + 4 f(X(2 i + 1, n)) + f(X(2 i + 2, n)))|
|   /                                                            |
|  -----                                                         |
\  i = 0                                                         /
--------------------------------------------------------------------------------
# Let S[n] denote Simpson's rule with n subdivisions.  Assume that the integrand f(t) is 4
# times continuously differentiable.  Then the error in Simpson's rule is given by a constant
# times a fourth order derivative.
--------------------------------------------------------------------------------
> Int(f(t),t=a..b)=S[n]-(D@@4)(f)(u)*(b-a)^5/(180*n^4);

b
/                         (4)              5
|                         D   (f)(u) (b - a)
|  f(t) dt = S[n] - 1/180 -------------------
|                                   4
/                                   n
a
--------------------------------------------------------------------------------
# Here u is some unknown point between a and b.
--------------------------------------------------------------------------------
> f:=t->sin(t)/t;

sin(t)
f := t -> ------
t
--------------------------------------------------------------------------------
> Diff(f(u),u\$4)=(D@@4)(f)(u);

sin(u)                sin(u)     cos(u)      sin(u)      cos(u)      sin(u)
Diff(------, u, u, u, u) = ------ + 4 ------ - 12 ------ - 24 ------ + 24 ------
u                     u          2           3           4           5
u           u           u           u
--------------------------------------------------------------------------------
# We do not know the value of u so, if we are to use this in the error formula then we must
# estimate the maximum absolute value on the interval u=0..x.  This can be quite tricky.
# Fortunately there is another way.
--------------------------------------------------------------------------------
> int(cos(t*x),x=0..1)=Int(cos(t*x),x=0..1);

1
/
sin(t)    |
------ =  |  cos(t x) dx
t      |
/
0
#
--------------------------------------------------------------------------------
# In fact let's redefine our function f(t) by this formula.  This has the advantage that Maple
# will now easily compute the value of f(t) at t=0.  Another great advantage is that we can
# easily compute derivatives of f(t) because we can differentiate with respect to x inside the
# integral sign.
--------------------------------------------------------------------------------
> f:=t->Int(cos(t*x),x=0..1);

1
/
|
f := t ->  |  cos(t x) dx
|
/
0
--------------------------------------------------------------------------------
> f(0);

1
/
|
|  1 dx
|
/
0
--------------------------------------------------------------------------------
> f(t);

1
/
|
|  cos(t x) dx
|
/
0
--------------------------------------------------------------------------------
> Diff(f(t),t\$4)=Int(diff(cos(t*x),t\$4),x=0..1);

1                              1
/                              /
|                              |            4
Diff( |  cos(t x) dx, t, t, t, t) =  |  cos(t x) x  dx
|                              |
/                              /
0                              0
--------------------------------------------------------------------------------
# From this we easily get an estimate for the fourth derivative.
--------------------------------------------------------------------------------
> abs((D@@4)(f)(u))<int(t^4,t=0..1);

1
/
|            4
abs( |  cos(u x) x  dx) < 1/5
|
/
0
--------------------------------------------------------------------------------
# Applying this to the Sine integral for positive x we have
--------------------------------------------------------------------------------
> abs(Si(x)-S[n])<(x^(5)/900)/n^4;

5
x
abs(Si(x) - S[n]) < 1/900 ----
4
n
--------------------------------------------------------------------------------
> for n from 0 to 100 do x[n]:=n*(0.1) od:n:='n':
--------------------------------------------------------------------------------
> abs(Int(sin(t)/t,t=x[n]..x[n+1])-S(2))<(0.1)^(5)/(900*16);

x[n + 1]
/
|      sin(t)                            -9
abs(   |      ------ dt - S(2)) < .6944444444*10
|         t
/
x[n]
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> for n from 0 to 99 do a:=x[n]:  b:=x[n+1]: J[n]:=Simpson(2) od:n:='n':
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;

elapsedtime := 19.533
--------------------------------------------------------------------------------
> J[0];

.09994446182
--------------------------------------------------------------------------------
> evalf(Int(sin(t)/t,t=0..0.1));

.09994446111
--------------------------------------------------------------------------------
> Si(0.1);

.09994446111
--------------------------------------------------------------------------------
> J[1];

.09961162812
--------------------------------------------------------------------------------
> evalf(Int(sin(t)/t,t=0.1..0.2));

.09961162742
--------------------------------------------------------------------------------
> J[0]+J[1];

.1995560899
--------------------------------------------------------------------------------
> evalf(Int(sin(t)/t,t=0.0..0.2));

.1995560885
--------------------------------------------------------------------------------
> Si(0.2);

.1995560885
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> for n from 0 to 99 do sum(J[k],k=0..n):Int(sin(t)/t,t=0..x[n+1])=Si(x[n+1]):sum(J[k],k=0..n)-Si(x[n+1]):(n+1)*(0.1)^(5)/(900*16) od:n:='n':
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;

elapsedtime := 5.000
--------------------------------------------------------------------------------
> for n from 0 to 99 do if modp(n,5)=4 then print(sum(J[k],k=0..n),Int(sin(t)/t,t=0..x[n+1])=Si(x[n+1]),sum(J[k],k=0..n)-Si(x[n+1]),(n+1)*(0.1)^(5)/(900*16)) fi od;n:='n':

.5
/
|   sin(t)                         -8                -8
.4931074214,  |   ------ dt = .4931074180, .34*10  , .3472222222*10
|      t
/
0

1.0
/
|   sin(t)                         -8                -8
.9460830765,  |   ------ dt = .9460830704, .61*10  , .6944444444*10
|      t
/
0

1.5
/
|   sin(t)                        -8                -7
1.324683539,  |   ------ dt = 1.324683531, .8*10  , .1041666667*10
|      t
/
0

2.0
/
|   sin(t)                        -8                -7
1.605412986,  |   ------ dt = 1.605412977, .9*10  , .1388888888*10
|      t
/
0

2.5
/
|   sin(t)                        -8                -7
1.778520181,  |   ------ dt = 1.778520173, .8*10  , .1736111111*10
|      t
/
0

3.0
/
|   sin(t)                        -8                -7
1.848652534,  |   ------ dt = 1.848652528, .6*10  , .2083333333*10
|      t
/
0

3.5
/
|   sin(t)                        -8                -7
1.833125402,  |   ------ dt = 1.833125399, .3*10  , .2430555555*10
|      t
/
0

4.0
/
|   sin(t)                                    -7
1.758203139,  |   ------ dt = 1.758203139, 0, .2777777777*10
|      t
/
0

4.5
/
|   sin(t)                         -8                -7
1.654140412,  |   ------ dt = 1.654140414, -.2*10  , .3124999999*10
|      t
/
0

5.0
/
|   sin(t)                         -8                -7
1.549931242,  |   ------ dt = 1.549931245, -.3*10  , .3472222222*10
|      t
/
0

5.5
/
|   sin(t)                         -8                -7
1.468724070,  |   ------ dt = 1.468724073, -.3*10  , .3819444444*10
|      t
/
0

6.0
/
|   sin(t)                         -8                -7
1.424687549,  |   ------ dt = 1.424687551, -.2*10  , .4166666666*10
|      t
/
0

6.5
/
|   sin(t)                                    -7
1.421794274,  |   ------ dt = 1.421794274, 0, .4513888888*10
|      t
/
0

7.0
/
|   sin(t)                        -8                -7
1.454596616,  |   ------ dt = 1.454596614, .2*10  , .4861111110*10
|      t
/
0

7.5
/
|   sin(t)                        -8                -7
1.510681535,  |   ------ dt = 1.510681531, .4*10  , .5208333333*10
|      t
/
0

8.0
/
|   sin(t)                        -8                -7
1.574186828,  |   ------ dt = 1.574186822, .6*10  , .5555555555*10
|      t
/
0

8.5
/
|   sin(t)                        -8                -7
1.629597107,  |   ------ dt = 1.629597100, .7*10  , .5902777777*10
|      t
/
0

9.0
/
|   sin(t)                        -8                -7
1.665040084,  |   ------ dt = 1.665040076, .8*10  , .6249999999*10
|      t
/
0

9.5
/
|   sin(t)                        -8                -7
1.674463350,  |   ------ dt = 1.674463342, .8*10  , .6597222221*10
|      t
/
0

10.0
/
|    sin(t)                        -8                -7
1.658347600,  |    ------ dt = 1.658347594, .6*10  , .6944444444*10
|       t
/
0

```