# WORKSHEET #28
#
# INFINITE SERIES
#
#
--------------------------------------------------------------------------------
> restart;
--------------------------------------------------------------------------------
# If a(n) are the terms of a series then the infinite sum is just the limit of the partial sums,
# Sum(a(n),n=1..infinity)=limit(S(N),N=infinity), where S(N)=Sum(a(n),n=1..N).
--------------------------------------------------------------------------------
> A:=Sum(a(n),n=1..infinity);
infinity
-----
\
A := ) a(n)
/
-----
n = 1
--------------------------------------------------------------------------------
> value(A);
infinity
-----
\
) a(n)
/
-----
n = 1
--------------------------------------------------------------------------------
> S:=N->Sum(a(n),n=1..N);
N
-----
\
S := N -> ) a(n)
/
-----
n = 1
--------------------------------------------------------------------------------
# EXAMPLE:
#
> a:=n->1/n^2;
1
a := n -> ----
2
n
--------------------------------------------------------------------------------
> A;
infinity
-----
\ 1
) ----
/ 2
----- n
n = 1
--------------------------------------------------------------------------------
> value(A);
2
1/6 Pi
--------------------------------------------------------------------------------
# Maple knows the exact answers to some series.
--------------------------------------------------------------------------------
> evalf(A);
1.644934067
> evalf(value(A));
1.644934068
--------------------------------------------------------------------------------
# We get slightly different answers because Maple is using different procedures to arrive at
# the answers. In the first case Maple uses some numerical procedure to add up the series,
# whereas in the second case Maple evaluates 1/6*Pi^2. To see which is more accurate we
# set Digits to 15 and redo the calculations.
--------------------------------------------------------------------------------
> Digits:=15:evalf(A);evalf(value(A));
1.64493406684823
1.64493406684823
--------------------------------------------------------------------------------
# Thus numerically evaluating the series is better than knowing its exact value.
--------------------------------------------------------------------------------
> Digits:=10:
--------------------------------------------------------------------------------
> a:='a':
# EXAMPLE:
--------------------------------------------------------------------------------
> a:=n->1/n^(2*k);
1
a := n -> ------
(2 k)
n
--------------------------------------------------------------------------------
> for k from 1 to 6 do value(A) od;k:='k':
2
1/6 Pi
4
1/90 Pi
6
1/945 Pi
8
1/9450 Pi
10
1/93555 Pi
691 12
--------- Pi
638512875
--------------------------------------------------------------------------------
# Maple knows these series.
--------------------------------------------------------------------------------
> a:='a':
--------------------------------------------------------------------------------
# EXAMPLE:
> a:=n->1/n^3;
1
a := n -> ----
3
n
--------------------------------------------------------------------------------
> A;
infinity
-----
\ 1
) ----
/ 3
----- n
n = 1
--------------------------------------------------------------------------------
> value(A);
Zeta(3)
--------------------------------------------------------------------------------
# This is the famous Riemann zeta function. The definition of Zeta(s) for s>1 is:
--------------------------------------------------------------------------------
> ?Zeta
--------------------------------------------------------------------------------
> Zeta(s)=Sum(1/n^s,n=1..infinity);
infinity
-----
\ 1
Zeta(s) = ) ----
/ s
----- n
n = 1
--------------------------------------------------------------------------------
# Putting s=1 in this series gives the harmonic series, which we know diverges. This
# definition is valid for complex numbers with real part>1.
--------------------------------------------------------------------------------
> Zeta(1);
Error, (in Zeta) singularity encountered
--------------------------------------------------------------------------------
# Zeta(s) has a singularity at s=1. This is not s surprise, but the next calculation is.
--------------------------------------------------------------------------------
> evalf(Zeta(1/2));
-1.460354509
--------------------------------------------------------------------------------
# The series in the definition of Zeta(s) diverges for s=-1/2. We can in fact extend the
# definition of Zeta(s) to other values of s, but with a different formula. The following
# equation is valid for s>0, s different from 1. The RHS converges for s>0 by the AST, and
# therefore this formula can be taken as the definition of Zeta(s) for s>0, s different from 1.
--------------------------------------------------------------------------------
> (1-2/2^(s-1))*Zeta(s)=Sum((-1)^(n-1)/n^s,n=1..infinity);
infinity
----- (n - 1)
/ 2 \ \ (-1)
|1 - --------| Zeta(s) = ) -----------
| (s - 1)| / s
\ 2 / ----- n
n = 1
--------------------------------------------------------------------------------
# There are other formulas that allow us to define Zeta(s) for all complex numbers s
# different from 1.
--------------------------------------------------------------------------------
> for k from 0 to 20 do Zeta(-k) od;k:='k':
-1/2
-1/12
0
1/120
0
-1/252
0
1/240
0
-1/132
0
691
-----
32760
0
-1/12
0
3617
----
8160
0
43867
- -----
14364
0
174611
------
6600
0
--------------------------------------------------------------------------------
> Zeta(-100);Zeta(-102);Zeta(-104);Zeta(-106);Zeta(-108);Zeta(-110);
0
0
0
0
0
0
--------------------------------------------------------------------------------
# The values of Zeta(s) at s=-2, -4, -6,... all seem to be zero, whereas the values of Zeta(s)
# at s=-1, -3, -5,... seem to be rational numbers which are negative for s=-1, -5, -9, and
# positive for s=-3, -7, -11,... The real importance of Zeta(s) concerns its values on the
# line Real(s)=1/2. Here we are thinking of s as a complex number with real part 1/2.
--------------------------------------------------------------------------------
> for k from 1 to 20 do print(Zeta(1/2+k*I)=evalf(Zeta(1/2+k*I))) od;k:='k':
Zeta(1/2 + I) = .1439364271 - .7220997435 I
Zeta(1/2 + 2 I) = .4405456503 - .3116463384 I
Zeta(1/2 + 3 I) = .5327366710 - .07889651343 I
Zeta(1/2 + 4 I) = .6067837645 + .09111213997 I
Zeta(1/2 + 5 I) = .7018123712 + .2310380084 I
Zeta(1/2 + 6 I) = .8372238081 + .3402183969 I
Zeta(1/2 + 7 I) = 1.021434231 + .3961894978 I
Zeta(1/2 + 8 I) = 1.241615106 + .3600475884 I
Zeta(1/2 + 9 I) = 1.447642452 + .1918030128 I
Zeta(1/2 + 10 I) = 1.544895220 - .1153364653 I
Zeta(1/2 + 11 I) = 1.419650152 - .4876753540 I
Zeta(1/2 + 12 I) = 1.015936651 - .7451124722 I
Zeta(1/2 + 13 I) = .4430047825 - .6554830983 I
Zeta(1/2 + 14 I) = .02224114261 - .1032581233 I
Zeta(1/2 + 15 I) = .1471099070 + .7047522416 I
Zeta(1/2 + 16 I) = .9385454085 + 1.216587816 I
Zeta(1/2 + 17 I) = 1.946654834 + .8954051920 I
Zeta(1/2 + 18 I) = 2.329154873 - .1888660058 I
Zeta(1/2 + 19 I) = 1.622495793 - 1.160117493 I
Zeta(1/2 + 20 I) = .4299138604 - 1.064291443 I
--------------------------------------------------------------------------------
# Something interesting is happening near s=1/2+14*I.
--------------------------------------------------------------------------------
> for k from 0 to 10 do print(Zeta(1/2+(14+k/10)*I)=evalf(Zeta(1/2+(14+k/10)*I))) od;k:='k':
Zeta(1/2 + 14 I) = .02224114261 - .1032581233 I
141
Zeta(1/2 + --- I) = .004698400183 - .02705828237 I
10
Zeta(1/2 + 71/5 I) = - .006816218159 + .05159699098 I
143
Zeta(1/2 + --- I) = - .01198782431 + .1322313685 I
10
Zeta(1/2 + 72/5 I) = - .01053795573 + .2143301116 I
Zeta(1/2 + 29/2 I) = - .002229099379 + .2973425897 I
Zeta(1/2 + 73/5 I) = .01313091049 + .3806853143 I
147
Zeta(1/2 + --- I) = .03568467788 + .4637454932 I
10
Zeta(1/2 + 74/5 I) = .06552127293 + .5458850989 I
149
Zeta(1/2 + --- I) = .1026725628 + .6264454436 I
10
Zeta(1/2 + 15 I) = .1471099070 + .7047522416 I
--------------------------------------------------------------------------------
> evalf(Zeta(1/2+14.134725*I));
-7 -6
.1767429836*10 - .1110202889*10 I
--------------------------------------------------------------------------------
# Zeta(s) has a zero near s=14.134725.
--------------------------------------------------------------------------------
> (1-1/(2^(s-1)))*Zeta(s)=sum((-1)^(n-1)/n^s,n=1..infinity);
infinity
----- (n - 1)
/ 1 \ \ (-1)
|1 - --------| Zeta(s) = ) -----------
| (s - 1)| / s
\ 2 / ----- n
n = 1
--------------------------------------------------------------------------------
> n^(-s)=exp(-s*ln(n));\
(- s)
n = exp(- s ln(n))
--------------------------------------------------------------------------------
> s:=sigma+I*t;
s := sigma + I t
--------------------------------------------------------------------------------
> (1-1/(2^(s-1)))*Zeta(s)=sum((-1)^(n-1)/n^s,n=1..infinity);
infinity
----- (n - 1)
/ 1 \ \ (-1)
|1 - ------------------| Zeta(sigma + I t) = ) --------------
| (sigma + I t - 1)| / (sigma + I t)
\ 2 / ----- n
n = 1
--------------------------------------------------------------------------------
> n^(-sigma-I*t)=n^(-sigma)*(cos(t*ln(n))-I*sin(t*ln(n)));
(- sigma - I t) (- sigma)
n = n (cos(t ln(n)) - I sin(t ln(n)))
--------------------------------------------------------------------------------
> (1-1/(2^(s-1)))*value(Zeta(s))=sum((-1)^(n-1)*n^(-sigma)*(cos(t*ln(n))-I*sin(t*ln(n))),n=1..infinity);
/ 1 \
|1 - ------------------| Zeta(sigma + I t) =
| (sigma + I t - 1)|
\ 2 /
infinity
-----
\ (n - 1) (- sigma)
) (-1) n (cos(t ln(n)) - I sin(t ln(n)))
/
-----
n = 1
--------------------------------------------------------------------------------
> sigma:=1/2;
sigma := 1/2
--------------------------------------------------------------------------------
> (1-1/(2^(s-1)))*value(Zeta(s))=sum((-1)^(n-1)*n^(-sigma)*(cos(t*ln(n))-I*sin(t*ln(n))),n=1..infinity);
/ 1 \
|1 - --------------| Zeta(1/2 + I t) =
| (- 1/2 + I t)|
\ 2 /
infinity
----- (n - 1)
\ (-1) (cos(t ln(n)) - I sin(t ln(n)))
) -------------------------------------------
/ 1/2
----- n
n = 1
--------------------------------------------------------------------------------
> eta:=t->sum((-1)^(n-1)*n^(-1/2)*(cos(t*ln(n))-I*sin(t*ln(n))),n=1..infinity);
infinity
----- (n - 1)
\ (-1) (cos(t ln(n)) - I sin(t ln(n)))
eta := t -> ) -------------------------------------------
/ 1/2
----- n
n = 1
--------------------------------------------------------------------------------
# The zeros of Zeta(s) on the line Real(s)=1/2 are the solutions of eta(t)=0.
--------------------------------------------------------------------------------
> s:='s':
--------------------------------------------------------------------------------
> a(n);
1
----
3
n
--------------------------------------------------------------------------------
> A;
infinity
-----
\ 1
) ----
/ 3
----- n
n = 1
--------------------------------------------------------------------------------
> S(100);
100
-----
\ 1
) ----
/ 3
----- n
n = 1
--------------------------------------------------------------------------------
> value(");
1/2 Psi(2, 101) + Zeta(3)
--------------------------------------------------------------------------------
> ?psi
--------------------------------------------------------------------------------
> ?GAMMA
--------------------------------------------------------------------------------
> GAMMA(s) = int( exp(-t)*t^(s-1), t=0..infinity );
infinity
/
| (s - 1)
GAMMA(s) = | exp(- t) t dt
|
/
0
--------------------------------------------------------------------------------
> Psi(x) = Diff( ln(GAMMA(x)), x );
d
Psi(x) = ---- ln(GAMMA(x))
dx
--------------------------------------------------------------------------------
> Psi(x)= Diff(GAMMA(x), x ) / GAMMA(x);
d
---- GAMMA(x)
dx
Psi(x) = -------------
GAMMA(x)
--------------------------------------------------------------------------------
> Psi(n,x) = Diff(Psi(x), x$n );
Psi(n, x) = Diff(Psi(x), x$n)
--------------------------------------------------------------------------------
> evalf(Psi(2,101))/2;
-.00004950249992
--------------------------------------------------------------------------------
> Zeta(3)=Sum(1/n^3,n=1..100)-(1/2)*Psi(2,101);
/ 100 \
|----- |
| \ 1 |
Zeta(3) = | ) ----| - 1/2 Psi(2, 101)
| / 3 |
|----- n |
\n = 1 /
--------------------------------------------------------------------------------
> evalf(");
1.202056903 = 1.202056903
--------------------------------------------------------------------------------
> evalf(S(100));
1.202007401
--------------------------------------------------------------------------------
> evalf(Zeta(3));
1.202056903
--------------------------------------------------------------------------------
> "-"";
.000049502
--------------------------------------------------------------------------------
> 1.0/(2*(101)^2);evalf(A-S(100));1.0/(2*(100)^2);
.00004901480247
.000049502
.00005000000000
--------------------------------------------------------------------------------
> a:='a':
--------------------------------------------------------------------------------
> a:=n->1/n;
a := n -> 1/n
--------------------------------------------------------------------------------
> A;
infinity
-----
\
) 1/n
/
-----
n = 1
--------------------------------------------------------------------------------
> value(A);
infinity
> limit(S(n)-ln(n),n=infinity);
/ n \
|----- |
| \ |
limit | ) 1/n| - ln(n)
n -> infinity | / |
|----- |
\n = 1 /
> value(");
gamma
> ?gamma
> evalf(gamma);
.5772156649
> a:='a':
> a:=n->1/n^3;
1
a := n -> ----
3
n
> A;
infinity
-----
\ 1
) ----
/ 3
----- n
n = 1
--------------------------------------------------------------------------------
> A-S(N)<1/(2*N^2);
/infinity \ / N \
| ----- | |----- |
| \ 1 | | \ 1 | 1
| ) ----| - | ) ----| < ----
| / 3 | | / 3 | 2
| ----- n | |----- n | 2 N
\ n = 1 / \n = 1 /
--------------------------------------------------------------------------------
> 1/(2*(N+1)^2)<A-S(N);
/infinity \ / N \
| ----- | |----- |
1 | \ 1 | | \ 1 |
---------- < | ) ----| - | ) ----|
2 | / 3 | | / 3 |
2 (N + 1) | ----- n | |----- n |
\ n = 1 / \n = 1 /
--------------------------------------------------------------------------------
> evalf(1+1/(2000)^3);
1.000000000
--------------------------------------------------------------------------------
> for k from 1250 to 1270 do evalf(1+1/(k^3)) od;k:='k':
1.000000001
1.000000001
1.000000001
1.000000001
1.000000001
1.000000001
1.000000001
1.000000001
1.000000001
1.000000001
1.000000000
1.000000000
1.000000000
1.000000000
1.000000000
1.000000000
1.000000000
1.000000000
1.000000000
1.000000000
1.000000000
--------------------------------------------------------------------------------
> evalf(1+1/((1259)^3));evalf(1+1/((1260)^3));
1.000000001
1.000000000
--------------------------------------------------------------------------------
> evalf((2*10^9)^(1/3));
1259.921050
--------------------------------------------------------------------------------
> (1259.921050)^3;
10
.2000000001*10
--------------------------------------------------------------------------------
> evalf(1/(1259)^(3));
-9
.5010981619*10
--------------------------------------------------------------------------------
> evalf(1/(1260)^(3));
-9
.4999060177*10
--------------------------------------------------------------------------------
> solve(1/n^3<(.5)*10^(-9),n);
{n < 0}, {1259.921050 < n}
> s[0]:=0:for n from 1 to 2000 do s[n]:=evalf(s[n-1]+1/n^3) od:n:='n':\
for n from 1 to 2000 do if modp(n,200)=0 then print(n,s[n]) fi od;n:='n':
200, 1.202044467
400, 1.202053791
600, 1.202055522
800, 1.202056121
1000, 1.202056394
1200, 1.202056594
1400, 1.202056653
1600, 1.202056653
1800, 1.202056653
2000, 1.202056653
--------------------------------------------------------------------------------
> S(10);
10
-----
\ 1
) ----
/ 3
----- n
n = 1
--------------------------------------------------------------------------------
> value(S(10));
19164113947
-----------
16003008000
--------------------------------------------------------------------------------
> value(S(20));
336658814638864376538323
------------------------
280346265322438720204800
--------------------------------------------------------------------------------
> value(S(50));
3251958144170050385476280040137248195072824750447564876132422113
----------------------------------------------------------------
2705769231561873462410773718765386619862166782990923887104000000
--------------------------------------------------------------------------------
> value(S(100));
1/2 Psi(2, 101) + Zeta(3)
--------------------------------------------------------------------------------
> convert([seq(1/n^3,n=1..100)],`+`);
8147348333074350358307418186167251193151812233617221640689414939133128970409751\
9580221863303145356050828007873151451209887
--------------------------------------------------------------------------------
6778118278309249584865634509184402157173419063091459022933216137995025717082809\
8031102950264769178652556660142954086400000
--------------------------------------------------------------------------------
> numerator(");
numerator(
8147348333074350358307418186167251193151812233617221640689414939133128970409751\
9580221863303145356050828007873151451209887
--------------------------------------------------------------------------------
6778118278309249584865634509184402157173419063091459022933216137995025717082809\
8031102950264769178652556660142954086400000
)
--------------------------------------------------------------------------------
> denominator("");
denominator(
8147348333074350358307418186167251193151812233617221640689414939133128970409751\
9580221863303145356050828007873151451209887
--------------------------------------------------------------------------------
6778118278309249584865634509184402157173419063091459022933216137995025717082809\
8031102950264769178652556660142954086400000
)
--------------------------------------------------------------------------------
> n:='n':
--------------------------------------------------------------------------------
> ?convert
--------------------------------------------------------------------------------
> ?Psi;
> a:='a':
--------------------------------------------------------------------------------
> a:=n->(-1)^(n-1)*(ln(n))^3/n;
(n - 1) 3
(-1) ln(n)
a := n -> ------------------
n
--------------------------------------------------------------------------------
> A;
infinity
----- (n - 1) 3
\ (-1) ln(n)
) ------------------
/ n
-----
n = 1
--------------------------------------------------------------------------------
> for j from 1 to 10 do s[j]:=evalf(s[j-1]+a(j)) od: j:='j':
--------------------------------------------------------------------------------
> for j from 1 to 10 do print(j,s[j]) od;j:='j':
1, 0
2, -.1665123260
3, .2754773276
4, -.3905719762
5, .4432103360
6, -.5155010264
7, .5371158036
8, -.5868423974
9, .5917966776
10, -.6290104774
--------------------------------------------------------------------------------
> for j from 1 to 1000 do s[j]:=evalf(s[j-1]+a(j)) od: j:='j':
--------------------------------------------------------------------------------
> s[999];s[1000];
.1554416297
-.1741763023
--------------------------------------------------------------------------------
> value(A);
infinity
----- (n - 1) 3
\ (-1) ln(n)
) ------------------
/ n
-----
n = 1
--------------------------------------------------------------------------------
> evalf(");
Error, (in evalf/sum1/levinu) division by zero
--------------------------------------------------------------------------------
> xi:=s->(1-1/2^(s-1))*Zeta(s);
/ 1 \
xi := s -> |1 - --------| Zeta(s)
| (s - 1)|
\ 2 /
--------------------------------------------------------------------------------
> 'xi(s)'=Sum((-1)^(n-1)/n^s,n=1..infinity);
infinity
----- (n - 1)
\ (-1)
xi(s) = ) -----------
/ s
----- n
n = 1
--------------------------------------------------------------------------------
> Diff('xi(s)',s$3)=Sum((-1)^(n-1)*diff(n^(-s),s$3),n=1..infinity);
infinity
-----
\ (n - 1) (- s) 3
Diff(xi(s), s, s, s) = ) - (-1) n ln(n)
/
-----
n = 1
--------------------------------------------------------------------------------
> A:=-Subs(s=1,Diff('xi(s)',s$3));
A := - Subs(s = 1, Diff(xi(s), s, s, s))
--------------------------------------------------------------------------------
> diff(xi(s),s$3);
3 2
ln(2) Zeta(s) ln(2) Zeta(1, s) ln(2) Zeta(2, s)
-------------- - 3 ----------------- + 3 ----------------
(s - 1) (s - 1) (s - 1)
2 2 2
/ 1 \
+ |1 - --------| Zeta(3, s)
| (s - 1)|
\ 2 /
--------------------------------------------------------------------------------
> ''A''=limit(",s=1,right);
4 3 2
'A' = - 1/4 ln(2) + ln(2) gamma + 3 ln(2) gamma(2) + 3 ln(2) gamma(1)
--------------------------------------------------------------------------------
> evalf(");
A = .0094139502
--------------------------------------------------------------------------------
>
