# Approximating the sum of an infinite series

I started thing about this as I was preparing some material on the use of Matlab for teaching engineering mathematics – material about sequences and series. It’s a long time since I did much numerical analysis, so i was wondering how you would best approximate the sum of an infinite alternating series. We know by the alternating series test that if is a sequence of positive values for which

then the series

converges. However, simply adding terms one by one is, in general very slow and inefficient. For example, take the well known series

which converges to . Here’s a table of the sum of the first terms, and multiplied by 4, with correct digits in bold:

As you see, obtaining correct digits requires terms to be summed.

There are many different methods for obtaining an approximate sum of an infinite series; most of these are known as techniques for accelerating convergence or series acceleration.

#### Euler’s method

One of the first documented methods to accelerate the convergence of a series was given by Euler (but of course); he noticed that an alternating series

can be written as

This replaces the original series with a new alternating series which should converge faster, because the terms are smaller.

Now for some experiments. For this blog post I’m going to use the software Pari/GP which is designed primarily for computational number theory, but contains a whole lot of other goodies, in particular some routines for approximating series sums.

To obtain the examples above, for example with :

? s = sum(i=0,10^4,(-1)^i/(2*i+1);
? \p 12
? 4.0*s
%1 = 3.14169264359


Here the first line computes the sum. Since GP works in symbolic, and accurate representation mode as a default, the sum is computed as a massive fraction, which we don;t want displayed. Hence the semi-colon at the end of the line. The second line sets the precision, and the last line displays our approximation to .

We are going to work with the first 11 terms, and apply Euler’s trick. First, create all the terms, noting that we are dealing with the positive values; we don’t need to include the alternating signs:

? a = vector(11,n,1/(2*n-1))
%2 = [1, 1/3, 1/5, 1/7, 1/9, 1/11, 1/13, 1/15, 1/17, 1/19, 1/21]


Note that GP indexing starts at one, rather than zero. The alternating sum of these terms from above is 3.232315809405594, which does not even have one decimal place accuracy. We can create the new sequence by taking differences:

? b = a/2;
? for(i = 2, length(a), b[i] = (a[i]-a[i-1])/2)


and we can compare the values of terms in the original and new series:

? 1.0*matconcat(a;b)~

%3 =
[  1.00000000000    0.500000000000]

[ 0.333333333333   -0.333333333333]

[ 0.200000000000  -0.0666666666667]

[ 0.142857142857  -0.0285714285714]

[ 0.111111111111  -0.0158730158730]

[0.0909090909091  -0.0101010101010]

[0.0769230769231 -0.00699300699301]

[0.0666666666667 -0.00512820512821]

[0.0588235294118 -0.00392156862745]

[0.0526315789474 -0.00309597523220]

[0.0476190476190 -0.00250626566416]


The terms in the new series are an order of magnitude smaller than in the original series. Now add them as an alternating series:

? sgns = vector(11,i,(-1)^(i-1));
>> 4.0 * sum(i = 1, 11, b[i]*sgns[i])

%4 = 3.13707771417


This value is a great deal closer to than the first result!

This process can be repeated any number of times, depending on the number of terms of the original series we have at our disposal:

? c = b/2;
? for(i = 2, length(b), c[i] = (b[i]-b[i-1])/2)
? 4.0 * sum(i = 1, 11, c[i]*sgns[i])

%5 = 3.14209024550

? d = c/2;
? for(i = 2, length(c), d[i] = (c[i]-c[i-1])/2)
? 4.0 * sum(i = 1, 11, c[i]*sgns[i])

%6 = 3.14150053593


and already we have far greater accuracy than at the start, without computing any further terms of the series.

Another way of implementing this acceleration technique is to write the first term as

and consider partial sums instead of individual terms. The -th partial sum of the original series is

Now consider the -th partial sum of the series of differences:

Breaking this up according to the first and second terms in each of the
numerators, we obtain

For the series above, we start with the terms including the sign changes:

? N = 11;
? a = vector(N,i,(-1)^(i-1)/(2*i-1));


Now it will be convenient to create a function called cumsum, which produces the cumulative sum of elements of a vector:

cumsum(x) = {
s = x;
for(i = 2, length(x), s[i] = s[i-1]+x[i])
return(s)
}


Now we can use this to experiment with Euler’s method again:

? s = cumsum(a);
? t = concat(0, vector(N-1, i, (s[i]+s[i+1])/2);
? 4.0*matconcat(s;t)~

%7 =
[4.00000000000             0]

[2.66666666667 3.33333333333]

[3.46666666667 3.06666666667]

[2.89523809524 3.18095238095]

[3.33968253968 3.11746031746]

[2.97604617605 3.15786435786]

[3.28373848374 3.12989232989]

[3.01707181707 3.15040515041]

[3.25236593472 3.13471887590]

[3.04183961893 3.14710277682]

[3.23231580941 3.13707771417]


And again we can take the pairwise averages of all the elements of , and again:

? u = concat(0,vector(N-1, i, (t[i]+t[i+1])/2));
? v = concat(0,vector(N-1, i, (u[i]+u[i+1])/2));
? \p 6
? 4.0 * matconcat(s;t;u;v)~

%8 =
[4.00000       0       0        0]

[2.66667 3.33333 1.66667 0.833333]

[3.46667 3.06667 3.20000  2.43333]

[2.89524 3.18095 3.12381  3.16190]

[3.33968 3.11746 3.14921  3.13651]

[2.97605 3.15786 3.13766  3.14343]

[3.28374 3.12989 3.14388  3.14077]

[3.01707 3.15041 3.14015  3.14201]

[3.25237 3.13472 3.14256  3.14136]

[3.04184 3.14710 3.14091  3.14174]

[3.23232 3.13708 3.14209  3.14150]


We can write a small program to do this for the initial values:

euler_trans(x) = {
local(s,N,p,t);
s = x;
N = length(s);
p = Vec(s[N]);
for(i = 1, N,
t = concat(0,vector(N-1,i,(s[i]+s[i+1])/2));
p = concat(p,t[N]);
s = t;
);
return(p);
}


And now use this program:

? \r euler_trans.gp
? \p 10
? s = cumsum(a*1.0);
? Mat(4*euler_trans(s))~

%9 =
[3.232315809]

[3.137077714]

[3.142090245]

[3.141500536]

[3.141618478]

[3.141582188]

[3.141598683]

[3.141587686]

[3.141598683]

[3.141581088]

[3.141633874]

[3.139152897]


Note that better accuracy is in fact obtained not at the end, but about two-thirds of the way down. This observation (and the proof that it does indeed have higher accuracy), is the contribution of Adriaan van Wijngaarden, who published it in 1965.

So one last check, for 120 terms, stopping our loop at 80:

? \p 20
? N = 120;
? a = vector(N,i,(-1)^(i-1)/(2*i-1));
? s = cumsum(a*1.0);
? p = euler_trans(s);
? forstep(i = 1,floor(N*2/3),10,print(p[i]))

3.1332594649198298159
3.1415926535897932384
3.1415926535897932385
3.1415926535897932385
3.1415926535897932385
3.1415926535897932385
3.1415926535897932385
3.1415926535897932385


Our accuracy is at the limit of our set precision:

? print(abs(4*p[80]-Pi))

1.5281426560689737604 E-37


Notice that all of this was produced with just the first 120 terms of the series; at no time in the computation of the loop did we compute any further terms.

For comparison, the direct sum of these terms is a very poor approximation indeed:

? print(abs(4*s[N]-Pi))

0.0083331886699634225838


One last example, to compute

Before we attempt it ourselves, we can check its value with the sumalt command of Pari/GP:

? altsqrt = sumalt(i=2,(-1)^i/sqrt(i))
%10 = 0.39510135657836962975273408576404450022


As above,:

? N = 120;
? a = vector(N,i,(-1)^(i+1)./sqrt(i+1));
? s = cumsum(a); s(N)

%11 = 0.34974072346962855572717264600532619101

? abs(s[N]-altsqrt)

%12 = 0.045360633108741074025561439758718309212


So the initial sum of 120 elements is not very good – as we’d expect. So, let’s see how Euler-van Wijngaarden manages:

? p = euler_trans(s);
? forstep(i = 1, floor(N*2/3), 10, print(p[i]))

0.34974072346962855572717264600532619101
0.39510135657836962968976358865551977318
0.39510135657836962975273408576388436451
0.39510135657836962975273408576404450023
0.39510135657836962975273408576404450022
0.39510135657836962975273408576404450022
0.39510135657836962975273408576404450022
0.39510135657836962975273408576404450022


And in fact we have reached accuracy to our set precision even earlier:

? forstep(i = 1, floor(N*2/3), 10, print([i,abs(p[i]-altsqrt)]))

[1, 0.045360633108741074025561439758718309212]
[11, 6.2970497108524727047861082221352448351 E-20]
[21, 1.6013571021345509597 E-31]
[31, 4.408103815583578155 E-39]
[41, 0.E-38]
[51, 1.4693679385278593850 E-39]
[61, 2.938735877055718770 E-39]
[71, 4.408103815583578155 E-39]


Note again that the power of this method is that without computing any extra terms of the series, we have increased the precision of the sum from one decimal place to nearly 40.