пусть мы хотим посчитать сумму ряда, альтернативной формулы для которой сходу не видно
скажем, сумму обратных кубов (в отличие от суммы обратных квадратов или четвертых степеней она [гипотетически] не выражается через пи)
просуммируем первые 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) теперь называют дзета-функцией… Римана)