Популярно о науке

Previous Entry Share Next Entry
Update
ahiin
К задаче.

Уменьшил требуемую абсолютную погрешность до Лимит времени тот же, 10 секунд.

Во избежание.

  • 1

Самый интересный для меня момент

Что суммирования с конца оказалось достаточно при вычислении рядов и боле-менее разумной точности.

Это далеко не тривиальный момент, особенно если учесть, что в FP и абсолютная и относительная погрешности расположены неравномерно (грубо говоря в 3ке подряд идущих чисел могут быть два разных интервала).

Re: Самый интересный для меня момент

Я, кстати, подумал тут, что по сути ваш подход и мой (2002 года) вполне близки в итоге, просто мой вариант выражен более аналитически.
В пятницу выложу и теорию и код.

Re: Самый интересный для меня момент

Мой какой?
С квадратом точности и рядом эйлера или с кубом точности и интегралом?


ПС
С рядом эйлера -- это вроде то же самое, что вы закрыли вчера. Только у меня "погрешность" в явном виде прибавляется в коде, а в закрытом комменте она внесена в член ряда.



Edited at 2018-03-28 01:13 pm (UTC)

Re: Самый интересный для меня момент

Тот вариант, что в виде исходника.

Общие соображения №2 (да простит меня автор)

Вообще я не занимаюсь научными рассчётами (очень редко -- исключительно редко надо было что-то посчитать по работе, ну и в студенчестве). Поэтому если я напишу явную чушь -- прошу меня понять и простить поправить.

Лично для меня эти наблюдения были интересны:
1. Для подсчёта суммы ряда в FP достаточно складывать его элементы с конца. Как ни странно это обеспечивает достаточную точность (при всех других фокусах FP).

2. Методы повышения точности "пропорцией" -- не работают.
Условно у нас есть ряды: 1/(x+1)^2 < 1/(x*(x+a)) < 1/x^2.
Мы можем, в зависимости от значения "a" сказать -- к какой из границ ближе среднее значение:
- а = 1 -- посредине
- а = 0 -- на верхней границе
- 0 < a < 1 -- пропорция
Этот метод увеличивает точность не более, чем вдвое. Это мало.

3. Если нас интересуют какие-то обычные научно-инженерные рассчёты, с не растущей непредсказуемо ошибкой (т.е. относительная точность < 10^-17, прикидываться должно "за пару минут на обычном процессоре),То для прикидок:
- если ошибка обратно пропорциональна квадрату числа посчитанных членов -- слишком долго
- если ошибка обратно пропорциональна кубу числа посчитанных членов, то это в самый раз.

ПС
Спасибо за внимание, извините, если пишу нубство или ввожу в заблуждение ;)

Re: Общие соображения №2 (да простит меня автор)

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

Re: Общие соображения №2 (да простит меня автор)

Кстати,
код ваш опубликовать разрешаете?

Re: Общие соображения №2 (да простит меня автор)

Да конечно.
Мне вообще странны такие вопросы.


ПС
А какая ассимптота в итоговом решении, если не секрет (по-большому счёту интересует это какой-то вариант вложенных сходящихся рядов и скорость схождения там показательная? Или всё-таки 1 раз оценён остаток и что-то степенное)?

Потому, как в кубическом (то, что я назвал "в лоб", а вы "не в лоб" -- с явной квадратичной сходимостью) -- время для 10^-14 странно похоже на ваше 3миллисекунды.

Edited at 2018-03-29 10:50 am (UTC)

Re: Общие соображения №2 (да простит меня автор)

Люди бывают разные, всегда лучше спросить.

Я остановился на восьмой степени, но можно было и продолжить. Так как возник интерес, я сделал более точную оценку времени, выкинув затраты на вывод на экран. Реальное полное время счета около 80 микросекунд.


Re: Общие соображения №2 (да простит меня автор)

>> я остановился на 8й степени

Хм.. я правильно понимаю, что это всё-таки вложенные (всё более точные) ассимптоты, но процесс перехода от n_ой ассимптоты к n+1_ой -- требует привлечение аналитического аппарата?

Re: Общие соображения №2 (да простит меня автор)

Технически да, но структура перехода абсолютна идентична и требует лишь скопипастить константу из вольфрама. Собственно, на восьмой степени копипаст меня задолбал.

В принципе, формула длиной в 25 означенных констант даст точность 1е-15 используя всего 4 члена итогового ряда.

Re: Общие соображения №2 (да простит меня автор)

1/(x+const_k)^2k -- такая ассимптота?


ПС
В смысле сумма рядя 1/(x+const)^2k -- аналитическая?

Edited at 2018-03-29 12:26 pm (UTC)

Re: Общие соображения №2 (да простит меня автор)

Не. Все в итоге сводится к суммрованию рада 1/((k^m)*(k+x))
где m можно сделать, технически, сколь угодно большим.
Подробности завтра)

Edited at 2018-03-29 12:37 pm (UTC)

Re: Общие соображения №2 (да простит меня автор)

Сложение с конца (с маленьких слагаемых) хорошо борется с потерей значимости при сложении чисел слишком разных порядков. Есть еще более хардкорные методы, например https://en.wikipedia.org/wiki/Kahan_summation_algorithm

Re: Общие соображения №2 (да простит меня автор)

Спасибо за ссылку.


Я скорее о том, что:
а) понятно откуда может взяться ошибка при суммировании ряда.
б) на практике в данной задаче оказалось достаточным суммирование с конца.


ПС
Я сначала наколхозил что-то подобное (попроще -- уже просуммированное не разбивал пользуясь тем, что ряд растёт достаточно "гладко"), но потом обнаружил, что при простом суммировании с конца ошибка только в 16м знаке появляется.

Ну, раз ты значения выложил, то я схалявил на Питоне.



from __future__ import print_function

import math
import timeit

# https://ahiin.livejournal.com/271654.html

N = 1000000

truth = [
  1.6449340668482264,
  1.5346072449045607,
  1.4408788415467229,
  1.360082586782444,
  1.2895778007910417,
  1.2274112777602189,
  1.1721051961250153,
  1.1225193425357527,
  1.0777588727442429,
  1.0371109178506583,
  0.99999999999999989
]

def f(x):
     a = math.pi ** 2 / 6.0
     b = 0
     for k in xrange(N, 0, -1):
         b += 1.0 / (k + x) / k / k

     return a - x * b
  
def calc(i):
    p = truth[i]
    x = i / 10.0
    q = f(x)
    print('{:.1f} {:.16f} {:g}'.format(x, q, abs(q - p)))

def main():
    t1 = timeit.default_timer()
    print(' x ', '{:18}'.format('f(x)'), 'err')
    for n in xrange(11):
        calc(n)
    t2 = timeit.default_timer()

    print(t2 - t1, 'sec')

if __name__ == '__main__':
    main()




Результат


 x  f(x)               err
0.0 1.6449340668482264 0
0.1 1.5346072449046106 4.996e-14
0.2 1.4408788415468228 9.99201e-14
0.3 1.3600825867825941 1.50102e-13
0.4 1.2895778007912417 2.00062e-13
0.5 1.2274112777604689 2.50022e-13
0.6 1.1721051961253153 2.99982e-13
0.7 1.1225193425361026 3.49942e-13
0.8 1.0777588727446430 4.00124e-13
0.9 1.0371109178511084 4.50084e-13
1.0 1.0000000000005000 5.00155e-13
1.60336350154 sec




Это ровно то, что надо) Практически дословно мое решение 2002 года, с поправкой на язык.

Прикрыл до завтра.

  • 1
?

Log in

No account? Create an account