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