Точность вычислений с плавающей точкой в окрестности деления на ноль

Question

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

C(s,t) = P - Q/(s-1) + (Q + Q/(s-1))/pow(s,t/T)           (1)

где P, Q, T - константы.

У этого выражения есть особенность: при s==1 получаем 0 в знаменателе. Однако, при abs(s-1) <= EPSILON оно имеет предел и вырождается в выражение:

C(s->1, t) = P + Q*(1 - t/T)                              (2)

Собственно мой вопрос: как максимально точно посчитать выражение (1) в окрестностях, близких к s==1 ? Использовать его непосредственно в виде (1), или лучше как-то преобразовать его? Например в форме (3):

C(s,t) = P + Q*((s - pow(s,t/T)/(s - 1))/pow(s,t/T)       (3)

В последнем случае, чем ближе знаменатель к нулю, тем ближе и числитель к нулю, а значит, интуитивно полагаю, погрешность вычислений должна быть меньше, чем если считать по формуле (1). Но я в этом не уверен, по-этому решил спросить совета. Просьба к математикам и продвинутым знатокам FPU не судить меня строго и не спешить минусовать.

Answer 1

Это не ответ, а расширенный комментарий.

Вы уверены, что вам нужна "точная" формула? Я тут прикинул, разложив функцию возведения в степень в ряд Тейлора.

Ваша функция представляется таким рядом:

Здесь в скобках стоит обобщённый биномиальный коэффициент:

Это разложение сходится при |s-1| < 1

Я сделал прототип на Python

P=2.0
Q=1.0
T=1.0

def C(s,t):
    return P - Q/(s-1) + (Q + Q/(s-1))/pow(s,t/T)

# Обобщённый биномиальный коэффициент
def binom(a,n):
    res = 1
    for k in range(1,n+1):
        res *= (a-k+1)/k
    return res

# Вычисление разложением в ряд. По умолчанию ряд вычисляется до 4 степени
def C2(s,t, n=4):
    d = s - 1
    f = 0
    for i in range(1,n+1):
        f += binom(-t/T, i)*(d**(i-1))
    return P + Q + Q*(1+d)*f

# Сравнение исходной формулы и разложения в ряд до 20 степени (оверкил)
# Перебираем epsilon от 0.1 до 1e-14
for p in range(1, 15):
    s = 1+pow(10,-p)
    c1 = C(s,0.5)
    c2 = C2(s,0.5, n=20)
    print(f"e = 1e-{p}: s = {s}, c1 = {c1}, c2 = {c2}, abs(c2-c1) = {abs(c2-c1)}")

Вот что получилось

e = 1e-1: c1 = 2.4880884817015154, c2 = 2.4880884817015154, abs(c2-c1) = 0.0
e = 1e-2: c1 = 2.4987562112089137, c2 = 2.4987562112089026, abs(c2-c1) = 1.1102230246251565e-14
e = 1e-3: c1 = 2.4998750624608874, c2 = 2.4998750624609647, abs(c2-c1) = 7.72715225139109e-14
e = 1e-4: c1 = 2.4999875006251386, c2 = 2.499987500624961, abs(c2-c1) = 1.7763568394002505e-13
e = 1e-5: c1 = 2.4999987500195857, c2 = 2.4999987500062497, abs(c2-c1) = 1.333599897179738e-11
e = 1e-6: c1 = 2.499999874853529, c2 = 2.4999998750000625, abs(c2-c1) = 1.4653345203896606e-10
e = 1e-7: c1 = 2.499999986961484, c2 = 2.4999999875000007, abs(c2-c1) = 5.385167867189011e-10
e = 1e-8: c1 = 2.5, c2 = 2.49999999875, abs(c2-c1) = 1.2500001034254637e-09
e = 1e-9: c1 = 2.5, c2 = 2.499999999875, abs(c2-c1) = 1.2500001034254637e-10
e = 1e-10: c1 = 2.5, c2 = 2.4999999999874998, abs(c2-c1) = 1.2500223078859563e-11
e = 1e-11: c1 = 2.5, c2 = 2.49999999999875, abs(c2-c1) = 1.2501111257279263e-12
e = 1e-12: c1 = 2.5, c2 = 2.4999999999998748, abs(c2-c1) = 1.2523315717771766e-13
e = 1e-13: c1 = 2.5, c2 = 2.4999999999999876, abs(c2-c1) = 1.2434497875801753e-14
e = 1e-14: c1 = 2.515625, c2 = 2.4999999999999987, abs(c2-c1) = 0.015625000000001332

Видно, что почти для всех эпсилонов расхождение между исходной формулой и разложением в ряд не хуже чем 1e-9

Заметная разница только для 1e-14, но это и так предел точности представления вещественных чисел в double.

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

UPD: сравнение с приближённой формулой

Сравнил приближённую формулу C(s->1, t) = P + Q*(1 - t/T) с результатом разложения в ряд. В выводе ниже c2 результат разложения в ряд, c3 - значение приближенной формулы.

e = 1e-1: c3 = 2.5, c2 = 2.4880884817015154, abs(c2-c3) = 0.01191151829848458
e = 1e-2: c3 = 2.5, c2 = 2.4987562112089026, abs(c2-c3) = 0.0012437887910974332
e = 1e-3: c3 = 2.5, c2 = 2.4998750624609647, abs(c2-c3) = 0.00012493753903530802
e = 1e-4: c3 = 2.5, c2 = 2.499987500624961, abs(c2-c3) = 1.2499375038999005e-05
e = 1e-5: c3 = 2.5, c2 = 2.4999987500062497, abs(c2-c3) = 1.2499937502852276e-06
e = 1e-6: c3 = 2.5, c2 = 2.4999998750000625, abs(c2-c3) = 1.2499993751191596e-07
e = 1e-7: c3 = 2.5, c2 = 2.4999999875000007, abs(c2-c3) = 1.2499999257897798e-08
e = 1e-8: c3 = 2.5, c2 = 2.49999999875, abs(c2-c3) = 1.2500001034254637e-09
e = 1e-9: c3 = 2.5, c2 = 2.499999999875, abs(c2-c3) = 1.2500001034254637e-10
e = 1e-10: c3 = 2.5, c2 = 2.4999999999874998, abs(c2-c3) = 1.2500223078859563e-11
e = 1e-11: c3 = 2.5, c2 = 2.49999999999875, abs(c2-c3) = 1.2501111257279263e-12
e = 1e-12: c3 = 2.5, c2 = 2.4999999999998748, abs(c2-c3) = 1.2523315717771766e-13
e = 1e-13: c3 = 2.5, c2 = 2.4999999999999876, abs(c2-c3) = 1.2434497875801753e-14
e = 1e-14: c3 = 2.5, c2 = 2.4999999999999987, abs(c2-c3) = 1.3322676295501878e-15

Получается, что ошибки округления начинают сказываться только для epsilon порядка 1e-14, а при выбранной вами точности 0.00001 ошибка приближенной формулы на пять порядков хуже, чем у точной. Поэтому я могу лишь повторить "посыл" - не беспокойтесь о потери точности из-за малости отклонения от 1.

БЛОГ НА HUSL

Точность вычислений с плавающей точкой в окрестности деления на ноль

Ответы (1 шт):