Неточность чисел с плавающей точкой приводит к ошибке вычислений

Мой код решает квадратное уравнение (в игровой логике), чтобы решить задачу — найти смещение тика спутника вдоль орбиты подвижного объекта в пространстве.
И я столкнулся с ошибками в дискриминант (дальше D) расчет. Я напомню: D = b^2 - 4ac,
Как это орбита большого объекта, мой a,b & c числа порядка, как:

1E+8
1E+12
1E+16

Соответственно, b^2 номер заказа о 1E+24, & 4ac около 1E+24 тоже.
НО корень этого уравнения намного меньше чисел, потому что они просто координаты на сцене. Так корни о 1E+3 ... 1E+4,

Эта проблема (обновлено — конкретизировано): из-за плавающих значений поплавков (& двойники) b^2 & 4ac иметь неточность, которая достаточно мала (относительно этих очень больших чисел [измеренная абсолютная неточность имеет порядок около 1E+18]), Но, как D == разница среди них, поэтому, когда D является (со стороны больших значений) на значение порядка, как указано неточность (1E+18), его значение начинают колебаться в диапазоне около +1E+18 .. -1E+18 (то есть колеблющийся диапазон шире, чем [-100% .. + 100%] фактического значения!

Очевидно, что это колебание вызывает неправильные (даже неправильно направленные) смещения тиков. И мой спутник начинает качаться (и это ужасно)).

Примечание: когда я сказал «когда D приближается к нулю «на самом деле D все еще достаточно далеко от нуля, поэтому я не могу просто присвоить его нулю в этом диапазоне значений.

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

Мой вопрос: Как я могу попытаться решить мою проблему? Может быть, есть некоторые общие решения для моего случая? Большое спасибо за любые советы!

PS: все формулы хороши (я все просчитал в excel & получил правильные результаты, когда поплавки в моем коде не удалось).

PPS: я пробовал удваивать ставку с плавающей точкой (не все расчеты, но мой a, b & c двойники сейчас) & проблема не исчезла.

Обновлено: Я ошибся — запутался порядок заказов a, b & c, Так «b^2 номер заказа о 1E+16, & 4ac около 1E+28«было неправильно. Теперь это исправлено 1E+24 и то и другое. (Я написал это, чтобы уже написанные комментарии были понятны)

Обновление # 2: Раздел «Проблема» конкретизирован.

Обновление # 3: Реальный регистр значений (для справки):
Примечание: в качестве «точных значений» здесь я отмечаю значения, рассчитанные вручную в Excel.

a == 1.43963872E+8
b == 3.24884062357827E+12
c == 1.83291898112689E+16

//floats:
b^2 == 1.05549641E+25
4ac == 1.05549641E+25
D == 0.0
root:
y = -1.12835273E+4

//doubles:
b^2 == 1.0554965397412443E+25
4ac == 1.0554964543412880E+25
D == 8.5399956328598733E+17
roots:
y1 == -1.1280317962726038E+4
y2 == -1.1286737079932651E+4

//accurate values:
b^2 == 1.05549653974124E+25
4ac == 1.05549645434129E+25
D == 8.53999563285987E+17
roots:
y1 == -1.128031796E+4
y2 == -1.128673708E+4

Похоже, хорошо с двойными, но это не так, потому что здесь я дал только часть вычислений — здесь я начинаю с того же а, б & c значения, но их реальные значения в моем коде также рассчитываются. И содержат неточности, которые приводят к проблемам даже с двойниками.

1

Решение

Использование стандартной квадратичной формулы может дать «катастрофическое аннулирование», где вычитание 2 чисел одинаковой величины дает потерю точности.

Хитрость заключается в использовании альтернативной формулировки в таких случаях, см. Здесь:
https://math.stackexchange.com/a/311397

ОБНОВЛЕНИЕ: я неправильно понял ваш вопрос. Я думаю, что проблема, скорее всего, заключается в чувствительности вашего результата к вводимым числам. Давайте выберем, скажем

a = 4e8
b = -1e12
c = 6.2e14

для которых решения ~ 1138 и 1361. Теперь, если вы вычислите относительные производные. Я могу сделать это в Джулии путем автоматического дифференцирования с помощью ForwardDiff.jl пакет:

julia> import ForwardDiff.Dual

julia> function p(a,b,c)
D = sqrt(b^2-4*a*c)
(-b+D)/(2a), (-b-D)/(2a)
end

julia> p(a,Dual(b,b),c)
(Dual(1361.803398874989,15225.424859373757),Dual(1138.196601125011,-12725.424859373757))

julia> p(Dual(a,a),b,c)
(Dual(1361.803398874989,-8293.614129124373),Dual(1138.196601125011,5793.614129124373))

julia> p(a,b,Dual(c,c))
(Dual(1361.803398874989,-6931.8107302493845),Dual(1138.196601125011,6931.8107302493845))

Результатами здесь являются два решения и их масштабированные производные (то есть (df / dx) * x). Обратите внимание, что все они имеют порядок O (10000), поэтому, если вход имеет ошибку на 0,000001%, выход будет иметь ошибку на 0,1%.

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

5

Другие решения

Смотрите мой ответ на этот вопрос: Квадратное уравнение в Аде

Хитрость заключается в том, чтобы всегда использовать

x1 = (-b - sign(b) * sqrt(b^2 - 4ac)) / 2a

в качестве первого корня, и использовать

x1 * x2 = c / a

найти второй. Таким образом, вы застрахуетесь от случая, когда 4ac << b ^ 2 и -b + sqrt (дельта) демонстрирует катастрофическое аннулирование.

Если ваша предполагаемая проблема заключается в том, что b ^ 2 и 4ac имеют одинаковую величину, то дельта на самом деле невелика по сравнению с b, и у вас нет проблемы округления, и вам, возможно, следует перемасштабировать вашу проблему (оба решения очень близки к -b / 2a).

2

C ++ имеет стандартную функцию библиотеки математики fma() который предлагает простой способ сделать вычисление корней квадратного уравнения максимально точным в рамках данного типа с плавающей точкой с помощью надежного вычисления дискриминанта d = √ (b2 — 4ac):

/*
Compute a*b-c*d with error < 1.5 ulp

Claude-Pierre Jeannerod, Nicolas Louvet, and Jean-Michel Muller,
"Further Analysis of Kahan's Algorithm for the Accurate Computation of 2x2 Determinants".
Mathematics of Computation, Vol. 82, No. 284, Oct. 2013, pp. 2245-2264
*/
T diff_of_products (T a, T b, T c, T d)
{
T w = d * c;
T e = fma (-d, c, w);
T f = fma (a, b, -w);
return f + e;
}

/* George E. Forsythe, "How Do You Solve a Quadratic Equation"Stanford University Technical Report No. CS40 (June 16, 1966)
*/
T a, b, c;
T d = diff_of_products (b, b, 2*a, 2*c);
T x1 = 2*c / (-b - sqrt (d));
T x2 = 2*c / (-b + sqrt (d));

Слитая операция умножения-сложения (FMA) реализовано fma() отображается на одну аппаратную инструкцию на большинстве современных процессорных архитектур. Поскольку FMA вычисляет полный, необоснованный продукт двойной ширины перед добавлением, он используется для точного вычисления ошибки продукта.

Как намекал Саймон Бирн в его ответ, конкретная рассматриваемая проблема плохо обусловлена, и точные вычисления не могут это исправить, может только переформулировка основной математики.

1
По вопросам рекламы ammmcru@yandex.ru
Adblock
detector