## File: Example-of-accelerating-a-series.html

package info (click to toggle)
gsl-ref-html 2.3-1
• area: non-free
• in suites: bullseye, buster, sid
• size: 6,876 kB
• ctags: 4,574
• sloc: makefile: 35
 file content (162 lines) | stat: -rw-r--r-- 6,774 bytes parent folder | download
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162  GNU Scientific Library – Reference Manual: Example of accelerating a series

31.3 Examples

The following code calculates an estimate of \zeta(2) = \pi^2 / 6 using the series,

\zeta(2) = 1 + 1/2^2 + 1/3^2 + 1/4^2 + ...

After N terms the error in the sum is O(1/N), making direct summation of the series converge slowly.

#include <stdio.h> #include <gsl/gsl_math.h> #include <gsl/gsl_sum.h>  #define N 20  int main (void) {   double t[N];   double sum_accel, err;   double sum = 0;   int n;      gsl_sum_levin_u_workspace * w      = gsl_sum_levin_u_alloc (N);    const double zeta_2 = M_PI * M_PI / 6.0;      /* terms for zeta(2) = \sum_{n=1}^{\infty} 1/n^2 */    for (n = 0; n < N; n++)     {       double np1 = n + 1.0;       t[n] = 1.0 / (np1 * np1);       sum += t[n];     }      gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err);    printf ("term-by-term sum = % .16f using %d terms\n",            sum, N);    printf ("term-by-term sum = % .16f using %zu terms\n",            w->sum_plain, w->terms_used);    printf ("exact value      = % .16f\n", zeta_2);   printf ("accelerated sum  = % .16f using %zu terms\n",            sum_accel, w->terms_used);    printf ("estimated error  = % .16f\n", err);   printf ("actual error     = % .16f\n",            sum_accel - zeta_2);    gsl_sum_levin_u_free (w);   return 0; }

The output below shows that the Levin u-transform is able to obtain an estimate of the sum to 1 part in 10^10 using the first eleven terms of the series. The error estimate returned by the function is also accurate, giving the correct number of significant digits.

\$ ./a.out
term-by-term sum =  1.5961632439130233 using 20 terms term-by-term sum =  1.5759958390005426 using 13 terms exact value      =  1.6449340668482264 accelerated sum  =  1.6449340669228176 using 13 terms estimated error  =  0.0000000000888360 actual error     =  0.0000000000745912

Note that a direct summation of this series would require 10^10 terms to achieve the same precision as the accelerated sum does in 13 terms.