# Lesson 21. Series
# ------------------
#
# Series of various sorts are very useful in many types of
# computation, and Maple has a number of facilities for dealing
# with them.
#
# First of all, there is the "sum" function, which can be used
# for both finite sums and infinite series. The inert version
# is "Sum".
> Sum(1/k^2, k=1..6) = sum(1/k^2, k=1..6);
6
-----
\ 1 5369
) ---- = ----
/ 2 3600
----- k
k = 1
> Sum(1/k^2, k=1..infinity) =
> sum(1/k^2, k=1..infinity);
infinity
-----
\ 1 2
) ---- = 1/6 Pi
/ 2
----- k
k = 1
# Maple can find the sums of some infinite series in "closed
# form", i.e. as exact expressions. If it can't do this, "sum"
# just returns the series.
> sum(k^2/(k^4+k+1), k=1..infinity);
infinity
----- 2
\ k
) ----------
/ 4
----- k + k + 1
k = 1
# Even if there is no exact expression for the sum, "evalf" may
# be able to approximate it by numerical techniques.
> evalf(");
.9316843051
# If you just want to use the numerical techniques and not have
# Maple look for an exact answer, you can use "evalf" with "Sum"
# (in the same way as we used "evalf" with "Int" for integrals).
> evalf(Sum(k^2/(k^4+k+1), k=1..infinity));
.9316843051
# The main theoretical question about a series is whether or not
# it converges. Maple can sometimes figure this out, but is not
# very good at it.
> sum(a^n, n=1..infinity);
a
- -----
a - 1
# Not a good answer: it's valid if |a| < 1, but not otherwise.
> sum(2^n, n=1..infinity);
infinity
> sum((-1)^n, n=1..infinity);
infinity
-----
\ n
) (-1)
/
-----
n = 1
# Well, at least it didn't apply the formula in those cases.
# "evalf" is even less clever about convergence.
> evalf(Sum(2^n, n=1..infinity));
-2.000000000
# And even for convergent series, "evalf" doesn't always succeed.
> evalf(Sum(ln(k+1)/k^3, k=1..infinity));
FAIL(.9594129036)
# So how do we tell whether a series converges or diverges?
# There are a number of "convergence tests" that can be used:
# the integral test, comparison test, ratio test, root test and
# alternating series test are the ones that tend to be found in
# calculus texts. Maple can be used as a tool in applying any
# of them. For example, the ratio test says the following:
#
# If |x[n+1]/x[n]| has a limit L < 1 as n -> infinity, then
# sum(x[n], n=1..infinity)
# converges. If it has a limit L > 1 (possibly +infinity), then
# the series diverges.
#
# Maple can find limits with the "limit" command.
> S:= Sum(3^(k^2)/k!, k=1..infinity);
infinity 2
----- (k )
\ 3
S := ) -----
/ k!
-----
k = 1
> x:= k -> 3^(k^2)/k!;
2
(k )
3
x := k -> -----
k!
> ratio:= x(k+1)/x(k);
2
((k + 1) )
3 k!
ratio := --------------
2
(k )
(k + 1)! 3
> limit(ratio, k=infinity);
infinity
# So this one diverges. It's not hard to see why.
> simplify(ratio);
(2 k + 1)
3
----------
k + 1
# Let's try another one.
> x:= k -> k^k/(k!)^2;
k
k
x := k -> -----
2
(k!)
> ratio:= x(k+1)/x(k);
(k + 1) 2
(k + 1) (k!)
ratio := --------------------
2 k
((k + 1)!) k
> limit(ratio, k=infinity);
0
> ratio:= simplify(ratio);
/ k \(- k)
|-----|
\k + 1/
ratio := ------------
k + 1
# This one converges. What is the sum?
> S:= Sum(x(k), k=1..infinity);
infinity
----- k
\ k
S := ) -----
/ 2
----- (k!)
k = 1
> evalf(S);
3.548128226
# The sum of an infinite series is, by definition, the limit of
# the partial sums.
# Let's calculate some partial sums for this series.
> seq( evalf(Sum(x(k), k=1 .. nn)), nn = 1 .. 10);
1., 2., 2.750000000, 3.194444444, 3.411458333, 3.501458333, 3.533879244,
3.544199224, 3.547141317, 3.547900723
# This makes Maple's answer look plausible. The partial sums
# appear to be approaching
# a limit that is approximately "evalf"'s answer. What n should
# we use to get the
# sum S correct to 10 significant digits?
#
# We want to estimate R(n) = S - sum(x(k), k=1..n) = sum(x(k),
# k=n+1 .. infinity), the "tail" of the series.
# The ratio test gives us a clue for this: the ratio x(k+1)/x(k)
# was
> ratio;
/ k \(- k)
|-----|
\k + 1/
------------
k + 1
# Note that (k/(k+1))^(-k) = (1+1/k)^k goes to a limit of E as k
# -> infinity. It is in fact less than E for every k:
# ln((1+1/k)^k) = k ln(1+1/k) <= k * (1/k) = 1.
# So the ratio is less than E/(k+1). In particular, for k >= 10
# it is less than E/11.
# This means that for k >= N >= 10, we have the bound
> 'x(k)' <= (E/11)^(k-N) *'x(N)';
(k - N)
x(k) <= (1/11 E) x(N)
# And so for the tail of the series, we have an estimate:
> R(N) <= sum((E/11)^(k-N)*x(N), k=N+1 .. infinity);
(- N) N (N + 1)
(1/11 E) N (1/11 E)
R(N) <= - 11 --------------------------------
2
(N!) (E - 11)
> EstR:= expand(rhs("));
N
N E
EstR := - --------------
2
(N!) (E - 11)
# We want the right side to be less than, say, .5*10^(-9). For
# N=10 we have
> R10:= evalf(subs(N=10,EstR));
R10 := .0002492573473
# The ratio between one EstR and the next is bounded by about
# E/10, since it's the same
# ratio as x(N+1)/x(N). So we want
> .5*10^(-9) = R10 * (E/10)^(N-10);
-9 (N - 10)
.5000000000*10 = .0002492573473 (1/10 E)
> solve(",N);
20.07180907
> evalf(subs(N=20,EstR));
-11
.5814665202*10
# It's actually considerably better than our estimate. That
# might be expected.
> ApproxS:= evalf(sum(x(k), k=1..20), 15);
ApproxS := 3.54812822600477
# And this agrees with "evalf"'s answer for the infinite sum.
#
# That was a series that converged very quickly. For series
# that converge slowly,
# approximating the sum can be more difficult.
#
# Let's look at the one that "evalf" couldn't do.
> x:= k -> ln(k+1)/k^3;
ln(k + 1)
x := k -> ---------
3
k
# The ratio test is inconclusive on this one: x(k+1)/x(k) -> 1
# as k -> infinity.
# What does work to prove convergence would be a comparison
# test: compare it, e.g.
# to sum(1/k^2, k=1..infinity). We have 0 < x(k) < 1/k^2 since
# ln(k+1) < k for all k.
# But what does that say about the tail of the series?
> Rest := Sum(1/k^2, k=N+1 .. infinity);
infinity
-----
\ 1
Rest := ) ----
/ 2
----- k
k = N + 1
> value(Rest);
Psi(1, N + 1)
# Well, Maple actually has a function for this one. What N
# would we need for Rest <= .5*10^(-10)?
> solve("=.5*10^(-10));
# I was a bit optimistic to ask Maple to solve this exactly.
# What about a numerical answer?
> fsolve(value(Rest)=.5*10^(-10));
11
.2000000000*10
# Not a very encouraging result: I really don't want to add up
# 20 billion terms of the
# series! Fortunately there are other ways to do this that will
# require far fewer terms.
