# 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
>
