Ускорение сходимости

Недавно пришлось разобраться в вопросе об ускорении сходимости последовательностей. До этого я был в полной уверенности, что знаю достаточно об этом, а чего не знаю — и не надо никогда. Но оказалось, что бывает надо. Расскажу про один простой метод, описанный в статье Дэвида Бродхерста (Broadhurst). Статья, конечно, не про метод.

Рассмотрим в качестве подопытной свинки $\zeta$-функцию Римана:

$$\zeta(n)=\sum_{k=1}^{\infty}\frac 1{k^n}$$

Просуммировав $N$ членов суммы, получим примерно $P\approx (n-1)\log_{10}N$ десятичных знаков. Можно, конечно, "пришить хвост", т.е., аппроксимировать остаток суммы интегралом, но это не может изменить того факта, что сходимость — логарифмическая $P\propto \log N$. Меняется только коэффициент пропорциональности. Справедливости ради заметим, что для многих практических целей достаточно, чтобы этот коэффициент был достаточно большим. Но что, если нам все-же хочется уметь вычислять такие суммы с высокой точностью?

Замечательный метод ускорения сходимости для таких последовательностей описан в вышеупомянутой статье. Я сначала не особенно сильно на него надеялся. Интуитивно мне казалось, что если мы вычислили $N$ первых членов суммы, то лучшее, что мы с ними можем сделать — просто сложить. Фокус в том, что мы не просто знаем эти члены, но еще и знаем кое-что об остатке $r_k$. А именно, функция $k^{n-1}r_k$ раскладывается в ряд вблизи точки $k=\infty$. Это позволяет построить некоторую взвешенную сумму первых $N$ членов так, что $P\propto N$!

В качестве иллюстрации прилагаю маленький скрипт, написанный по мотивам вышеупомянутой статьи, который вычисляет $\zeta(n)$. Внимательный взгляд найдет пару отличий от описанного в статье. Поскольку написан он на javascript, то точность в какой-то момент теряется (точнее, не растет), но это проблемы вычислений с фиксированной точностью. Вбейте в форму ниже $N$ и $n$ и нажмите кнопку.

---

N= членов суммы для ζ() дают
при тупом суммировании, и
после ускорения.

---

Если кому-то стало интересно, предлагаю выяснить как выглядит вектор весов и определить коэффициент пропорциональности между P и N (экспериментально). Результат может удивить. Еще один вопрос — какой должна быть рабочая точность? Для удобства ниже — упрощенный (без учета потери точности) код скрипта
 function Zeta(N,n) {
  var terms= [0];
  //accumulating from above
  for(var i=N;i>0;i--) {
  terms.push(terms[terms.length-1]+1/Math.pow(i,n));
  };
  //acceleration starts
  for(var m=0;m<N;m++) {
  for(var k=0;k<(N-m);k++) {
  terms[k]=((N-k)*terms[k]-(N-k-m-n+1)*terms[k+1])/(m+n-1)
  };
  };
  //direct sum
  document.getElementById('direct').value=terms[N];
  //accelerated
  document.getElementById('accelerated').value=terms[N]-terms[0];
 }