продолжаем суммировать ряды, или заменяем грубую силу на ___

пусть мы хотим посчитать сумму ряда, альтернативной формулы для которой сходу не видно

скажем, сумму обратных кубов (в отличие от суммы обратных квадратов или четвертых степеней она [гипотетически] не выражается через пи)

просуммируем первые N членов… насколько мы далеки от правильного ответа? сумма оставшихся членов примерно равна интегралу 1/x³ от N до ∞, т.е. ~1/2N²

то есть даже для N=100500 мы получаем всего 10 правильных цифр, а до 100 правильных цифр дойти практически невозможно (требуется ~10^{50} слагаемых)…


…но подождите! «ваш прогноз выполняется с вероятностью 40%? так говорите все наоборот, и будет выполняться с вероятностью 60%» — раз мы знаем, насколько наш ответ далек от правильного, так давайте соответствующим образом скорректируем ответ!

и действительно, для N=100500 сумма первых N членов +1/2N² дает на 5 больше правильных цифр


можно ли оценить сумму (в данном случае, хвост ряда) точнее, чем просто интегралом? оценка «интеграл ≈ сумма» возникает при приближении площади под графиком площадью столбиков — а ясно, что лучше заменить прямоугольники трапециями

если вдуматься, это учит, что точность повысится, если вычесть половину N-го члена

на самом деле, +1/2N²-1/2N³ — это начало серии поправок, а возникающие коэффициенты — это (примерно) числа Бернулли… писал про этот сюжет (формулу Эйлера-Маклорена и числа Б.) в статье «Алгебра, геометрия и анализ сумм степеней последовательных чисел», https://www.mathnet.ru/rus/mp882

***

посмотрим, как это работает на практике

возьмем, скажем, 57 поправочных членов… а сумму не обязательно даже считать до N=100500, ограничимся N=500 — и уже этого хватает для 100 правильных цифр после запятой!

from mpmath import *
mp.dps = 102

from textwrap import wrap

def bernoulli_upto(m):
    ans = [mpf(1)]
    for k in range(1,m+1):
        b = mpf(0)
        for i in range(k):
            b -= ans[i]*binomial(k+1,i) / (k+1)
        ans.append(chop(b,tol=1e-10))
    return ans

def improve(pre,N,m):
    bernoulli = bernoulli_upto(m-1)
    for k in range(1,len(bernoulli)+1):
        pre += (bernoulli[k-1]*k) / (2*N**(k+1))
    return pre

def partial_sum(N):
    ans = mpf(0)
    for k in range(1,N+1):
        ans += mpf(1) / k**3
    return ans

def mprint(num,prec):
    intp, frcp = nstr(num,prec).split(".")
    print(f"{intp},",*wrap(frcp,width=5))

#N = 100500
#ans = partial_sum(N)
#mprint(ans,21) # всего 10 правильных цифр

N = 500
ans = partial_sum(N)
mprint(improve(ans,N,57),101) # все 100 цифр правильные!

***

так что, может вообще зафиксировать N, а только увеличивать количество поправочных членов (аргумент m функции improve)? степени N экспоненциально убывают, ну какие-то коэффициенты стоят при них…

проблема в том, что если начало последовательности Бернулли выглядит безобидно (1, -1/2, 1/6, -1/30…), то дальше эти числа начинают быстро — почти как факториалы! — расти (и для m=57 появляются коэффициенты ~10^{30})

ситуация достаточно тонкая: какое бы (конечное) число поправочных членов мы ни написали, при достаточно больших N они увеличивают точность вычисления — но если, наоборот, зафиксировать N, то ряд из поправочных членов расходится (в частности, с какого-то момента добавление новых поправочных членов ухудшает ситуацию!)

***

в прошлый раз речь шла о версии того же для знакопеременных рядов — и там коэффициенты поправок получались не просто рациональными, а целыми — и их можно прямо увидеть в подчеркнутых цифрах

про это поведение ряда Лейбница — взял из книжки «Mathematics by Experiment» (Borwein, Bailey) — и книжку вообще, и конкретно это место в ней мне показал Сережа Маркелов

а сама математика восходит к Эйлеру — которому очень хотелось посчитать 20 знаков суммы обратных квадратов, не складывая 10^{20} чисел

(после этого Эйлер нашел и точное значение этой суммы… а также ζ(2n), ζ(-k), и даже функциональное уравнение для ζ(x) — видимо поэтому ζ(x) теперь называют дзета-функцией… Римана)