# 1.22 Lesson 22

```# Lesson 22.  Tails of Series
# ---------------------------
> restart;
# In the last lesson, we wanted to find a numerical value for
# the sum of the following series, accurate to .5*10^(-10).
> x:= k -> ln(k+1)/k^3:
> S:= Sum(x(k), k=1..infinity);

infinity
-----
\      ln(k + 1)
S :=    )     ---------
/           3
-----       k
k = 1
> evalf(S);

FAIL(.9594129036)
> PartialSum:= N -> sum(x(k), k=1..N);

N
-----
\
PartialSum := N ->   )   x(k)
/
-----
k = 1
# The rather unhelpful result we began with was that the partial
# sum S[N] would work if N was 2*10^10.  But we can do much
# better.
# One improvement is to use the following estimate, based on the
# Integral Test:
# If f(t) is a positive, decreasing function of t, then
> Int(f(t), t=n-1 .. n) >= f(n);

n
/
|
f(n) <=    |  f(t) dt
|
/
n - 1
> f(n) >= Int(f(t), t=n .. n+1);

n + 1
/
|
|    f(t) dt <= f(n)
|
/
n
# And so if Int(f(t), t=N .. infinity) converges,
> Int(f(t),t=N .. infinity) >= Sum(f(n), n=N+1 .. infinity);

infinity          infinity
-----               /
\                 |
)     f(n) <=    |      f(t) dt
/                 |
-----             /
n = N + 1           N
> Sum(f(n), n=N+1 .. infinity) >= Int(f(t), t=N+1 .. infinity);

infinity             infinity
/                  -----
|                    \
|      f(t) dt <=     )     f(n)
|                    /
/                    -----
N + 1                n = N + 1
# Is our x(t) decreasing?  It's certainly positive.
> normal(D(x)(t));

- t + 3 ln(t + 1) t + 3 ln(t + 1)
- ---------------------------------
4
(t + 1) t
# For t >= 1, 3 ln(t+1) = ln((t+1)^3) >= ln(8) > 1, so 3 ln(t+1)
# t > t and D(x)(t) < 0.
# So the estimates apply.  The "tail" R(N) of the series:
> R:= N -> Sum(x(n), n=N+1 .. infinity);

infinity
-----
\
R := N ->     )     x(n)
/
-----
n = N + 1
# lies between the following bounds:
> Rupper:= unapply(int(x(t), t=N .. infinity), N);

2              2
ln(N + 1) N  - N - ln(N) N  - ln(N + 1)
Rupper := N -> - 1/2 ---------------------------------------
2
N
> Rlower:= unapply(Rupper(N+1),N);

2                            2
ln(N + 2) (N + 1)  - N - 1 - ln(N + 1) (N + 1)  - ln(N + 2)
Rlower := N -> - 1/2 -----------------------------------------------------------
2
(N + 1)
# For example, for N=100:
> evalf([Rlower(100), Rupper(100)]);

[.000251040, .000255591]
# Thus the true sum is between the following bounds:
> evalf([PartialSum(100)+Rlower(100), PartialSum(100)+Rupper(100)]);

[.9594041692, .9594087199]
# A good approximation to the sum would be the average of these:
> Approx:= evalf(PartialSum(100)+(Rlower(100)+Rupper(100))/2);

Approx := .9594064442
# And the absolute value of the error in this approximation
# would be no greater than
> Ebound:= evalf((Rupper(100)-Rlower(100))/2);

-5
Ebound := .2276*10
# This is a much better approximation than the partial sum
# itself, whose error is more than 2.5*10^(-4).
# How large an N would we need to get the error bound to be less
# than .5*10^(-10)?
> fsolve((Rupper(N)-Rlower(N))/2 = .5*10^(-10), N);

4376.287743
> evalf((Rupper(4377)-Rlower(4377))/2);

0
> evalf((Rupper(4377)-Rlower(4377))/2, 15);

-10
.4998*10
# So 4377 terms ought to do it.  This would be feasible to
# calculate, though it might take a little while.
> timer:= time():\
A4377:= evalf(PartialSum(4377)+(Rlower(4377)+Rupper(4377))/2, 20);\
time()-timer;

A4377 := .95940642259401994639

44.034
# Let's see if we can do better.
# The bounds we used were essentially comparing int(f(x),x=n ..
# n+1) to the values of f at the two endpoints: f(n) and f(n+1).
#  We're using an integral to approximate a sum, but we can
# relate this to what we've done before in approximating
# integrals by sums.  And using the value of f at the left or
# right endpoint is a very crude way to approximate the
# integral.  For an improvement, we might try the Trapezoid Rule:
> Int(f(t), t=n .. n+1) = (f(n) + f(n+1))/2 - (D@@2)(f)(c)/12;

n + 1
/
|                                              (2)
|    f(t) dt = 1/2 f(n) + 1/2 f(n + 1) - 1/12 D   (f)(c)
|
/
n
# where c is between n and n+1.  Adding these up for n from N to
# infinity, we get
> Int(f(t), t=N .. infinity) \
= f(N)/2 + Sum(f(n), n=N+1 .. infinity)  \
- Sum((D@@2)(f)(c[n])/12, n=N .. infinity);

infinity                      / infinity     \   /infinity                   \
/                          |  -----       |   | -----                     |
|                           |   \          |   |  \            (2)         |
|      f(t) dt = 1/2 f(N) + |    )     f(n)| - |   )     1/12 D   (f)(c[n])|
|                           |   /          |   |  /                        |
/                            |  -----       |   | -----                     |
N                            \n = N + 1     /   \ n = N                     /
# so the tail of the series is
> Sum(f(n), n=N+1 .. infinity) = Int(f(t), t=N .. infinity) - f(N)/2 + Sum((D@@2)(f)(c[n])/12, n=N .. infinity);

infinity         infinity                      /infinity                   \
-----              /                          | -----                     |
\                |                           |  \            (2)         |
)     f(n) =    |      f(t) dt - 1/2 f(N) + |   )     1/12 D   (f)(c[n])|
/                |                           |  /                        |
-----            /                            | -----                     |
n = N + 1          N                            \ n = N                     /
# Now if (D@@2)(f)(t) is a decreasing function of t, we have
# (D@@2)(f)(n) >= (D@@2)(f)(c[n]) >= (D@@2)(f)(n+1).
> Sum(f(n), n=N+1 .. infinity) \
>= Int(f(t), t=N .. infinity) - f(N)/2 \
+ Sum((D@@2)(f)(n)/12, n=N+1 .. infinity);

infinity                      / infinity                \     infinity
/                          |  -----                  |      -----
|                           |   \            (2)      |       \
|      f(t) dt - 1/2 f(N) + |    )     1/12 D   (f)(n)| <=     )     f(n)
|                           |   /                     |       /
/                            |  -----                  |      -----
N                            \n = N + 1                /    n = N + 1
> Sum(f(n), n=N+1 .. infinity) \
<= Int(f(t), t=N .. infinity) - f(N)/2 \
+ Sum((D@@2)(f)(n)/12, n=N .. infinity);

infinity          infinity                      /infinity                \
-----               /                          | -----                  |
\                 |                           |  \            (2)      |
)     f(n) <=    |      f(t) dt - 1/2 f(N) + |   )     1/12 D   (f)(n)|
/                 |                           |  /                     |
-----             /                            | -----                  |
n = N + 1           N                            \ n = N                  /
# Is this going to be helpful?  Do we have another infinite
# series to calculate?  Well, compare again to integrals.  Again
# if (D@@2)(f)(t) is decreasing,
# (D@@2)(f)(n) >= int((D@@2)(f)(t), t=n .. n+1)
#              >= (D@@2)(f)(n+1).
# And note that (if D(f)(t) -> 0 as t -> infinity),
> Int((D@@2)(f)(t), t=N .. infinity) = -D(f)(N);

infinity
/
|       (2)
|      D   (f)(t) dt = - D(f)(N)
|
/
N
> Sum(f(n), n=N+1 .. infinity) \
<= Int(f(t), t=N .. infinity) - f(N)/2 - D(f)(N-1)/12;

infinity          infinity
-----               /
\                 |
)     f(n) <=    |      f(t) dt - 1/2 f(N) - 1/12 D(f)(N - 1)
/                 |
-----             /
n = N + 1           N
> Sum(f(n), n=N+1 .. infinity) \
>= Int(f(t), t=N .. infinity) - f(N)/2 - D(f)(N+1)/12;

infinity                                           infinity
/                                                -----
|                                                  \
|      f(t) dt - 1/2 f(N) - 1/12 D(f)(N + 1) <=     )     f(n)
|                                                  /
/                                                  -----
N                                                n = N + 1
# So a good approximation to the sum would be the average of the
# two bounds, and a bound on the error in this approximation
# would be half the difference between the bounds.
> Approx:= (f,N) -> Sum(f(n),n=1..N) \
+ Int(f(t),t=N .. infinity) - f(N)/2 \
- (D(f)(N+1)+D(f)(N-1))/24;

/  N       \    infinity
|-----     |       /
| \        |      |
Approx := (f,N) -> |  )   f(n)| +    |      f(t) dt - 1/2 f(N)
| /        |      |
|-----     |     /
\n = 1     /     N

- 1/24 D(f)(N + 1) - 1/24 D(f)(N - 1)
> ErrEst:= (f,N) -> (D(f)(N+1)-D(f)(N-1))/24;

ErrEst := (f,N) -> 1/24 D(f)(N + 1) - 1/24 D(f)(N - 1)
# Does our x(t) meet the requirements?
> normal((D@@3)(x)(t));

3       2                        3                  2
- (- 47 t  - 81 t  - 36 t + 60 ln(t + 1) t  + 180 ln(t + 1) t  + 180 ln(t + 1) t

/         3  6
+ 60 ln(t + 1))  /  ((t + 1)  t )
/
# This is clearly negative at least when ln(t+1) >= 1 (true for
# t >= 2).  So (D@@2)(x)(t) is decreasing.
> ErrEst(x,N);

1                ln(N + 2)         1               ln(N)
------------------- - 1/8 --------- - ------------- + 1/8 --------
3               4               3              4
24 (N + 2) (N + 1)         (N + 1)    24 N (N - 1)        (N - 1)
> E10:=evalf(ErrEst(x,10));

E10 := .00001954662216
# What N would we need to make the error estimate less than
# .5*10^(-10)?
# We could try "fsolve" again, but here's a different approach.
# The "asympt" command gives an asymptotic formula for an
# expression as a variable goes to infinity.
> asympt(ErrEst(x,N),N);

ln(N) - 7/12       1
------------ + O(----)
5            6
N            N
# Ignoring the ln(N), if ErrEst(x,N) was proportional to 1/N^5,
# this is the N we would need:
> (E10/(.5*10^(-10)))^(1/5)*10;

131.3470559
# Trying it out:
> E131:= evalf(ErrEst(x,131));

-9
E131 := .111604983*10
# Not quite, try it again.
> (E131/(.5*10^(-10)))^(1/5)*131;

153.8204267
# This time add a little bit of a safety or "fudge factor".
> evalf(ErrEst(x,155));

-10
.49981479*10
# Success!  Here's the approximation.
> A155:= evalf( Approx(x,155),20);

A155 := .95940642259480509054
# To check how well this agrees with our previous approximation:
> evalf(A155-A4377,20);

-12
.78514415*10
>

```