# 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
