Погрешности расчета, накопленные при расчете частичной суммы гармонического ряда - PullRequest
Новая
       18

Погрешности расчета, накопленные при расчете частичной суммы гармонического ряда

0 голосов
/ 06 июня 2019

В настоящее время я использую два подхода для вычисления частичной суммы ряда гармоник:

  1. Грубое принуждение (итерация)
  2. Использование определенных формул, которые предлагают ярлык для решения

Проблема, которая меня беспокоит, заключается в следующем: двойное с точностью до 15 десятичных знаков и длинное двойное с точностью до 19 десятичных знаков.Я имею дело с частичными суммами от 5 до 50. Например, чтобы достичь суммы 5, вам потребуется около 83 членов гармонического ряда, а для достижения суммы 50 вам потребуется около 2911068851283175669760 членовгармонический ряд.В настоящее время я использую double.

При вычислении больших сумм (например,> 20) теряю ли я точность при вычислении таких сумм с использованием итераций и формул?

Если это так, будет ли этовыгодно использовать некоторые высокоточные библиотеки?

1 Ответ

1 голос
/ 07 июня 2019

Используя любой язык, который обеспечивает точную арифметику для (неограниченного) большого целого и дробного числа, мы можем попытаться легко оценить ошибку округления с помощью очень наивного кода, по крайней мере, для небольшого n (точная арифметика может тогда стать запретительной)

Я произвольно выберу здесь Squeak Smalltalk, чтобы проиллюстрировать это.

Сначала давайте определим замыкание, которое будет возвращать точную частичную сумму 

hn := [:n | (1 to: n) detectSum: [:i | i reciprocal]].

Затем еще одно, которое возвратит Float(двойная точность) по итерации, начиная с меньших членов: 

fn := [:n | (n to: 1 by: -1) detectSum: [:i | i asFloat reciprocal]].

А теперь еще одно замыкание, которое оценит ошибку в терминах ulp: 

un := [:n | | f | ((f := fn value: n) asFraction - (hn value: n)) abs / f ulp].

Теперь давайте найдемнаихудший случай для первых 1000 неполных сумм: 

i := (1 to: 1000) detectMax: [:n | un value: n ].
{i. un value: i}.

Результат - #(807 6.101840327204507), ошибка 6 ulp для H807.Это не так много, как вы можете видеть.

Если мы попытаемся сначала суммировать с большей гармоникой: 

fn := [:n | (1 to: n) detectSum: [:i | i asFloat reciprocal]].

Тогда результат будет #(471 5.714808765246509), 5,7 ulp, удивите немного меньше.

Мы можем попытаться взглянуть немного дальше с более поздними значениями fn: 

un value: 10000.
-> 1.8242087614993667

Используя мемоизацию, я нахожу, что максимальная ошибка составляет около 54 ulp для 100 000 первых членов (H ~ 7) и около26 ulp, если суммирование более аккуратное, начиная с меньших терминов.

Таким образом, до нескольких тысяч терминов нет смысла начинать суммирование с меньших терминов (которые позволяют запоминать, если вы собираетесь сканировать увеличение), и это не такнеобходимо использовать суммирование по Кахану, если только вы действительно не заботитесь о последних нескольких ульпах.

Кахан станет интересным для миллионов до миллиардов терминов, но для большей суммы формулы будут намного быстрее, а также более точными, еслиу вас есть хорошее приближение Euler Mascheroni, натуральный логарифм никогда не будет больше, чем несколько уловов.Но вы должны показать свой код, в противном случае это всего лишь предположение.

Например, эта наивная формулировка второго порядка 

fn := [:n | n ln + (0.5772156649015328606 + ((n*2.0) reciprocal - (n squared*12.0) reciprocal))].

даст ошибку около 5,5 ulp для малых значений n, но меньшечем 1,5 ulp выше 100 000 с приличным MacOs libm.

...