# 1.21 Lesson 21

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