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