浮点数不准确会导致计算错误

Question

我的代码解决二次方程(在游戏逻辑滴答中)来解决任务 - 找到沿着空间中可移动物体轨道的卫星刻度线偏移.我在判别(更远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仍然远离零时,所以我不能只在这个值范围内将它分配给zerro.

我考虑过使用定点计算(这可以帮助我解决问题).但是,建议不要在刻度逻辑中使用(因为它们的优化程度要低得多,而且可能会非常慢).

我的问题:我怎样才能解决我的问题？对我的案子可能有一些常见的解决方案吗？非常感谢任何建议!

PS:所有公式都很好(我在excel中计算了所有结果,并且在我的代码中浮动时失败了).

PPS:我尝试双打insted的花车(不是所有的计算,但我的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

看起来好像有双打,但事实并非如此,'我在这里只给出了部分计算 - 这里我从相同的a,b&c值开始,但我们的代码中的实际值也是计算出来的.并且包含不准确性,即使使用双打也会产生问题.

Answer 1

使用标准二次公式可以给出"灾难性的消除",其中减去相同幅度的2个数字会导致精度损失.

诀窍是在这种情况下使用替代配方,请参见此处:https: //math.stackexchange.com/a/311397

更新:我误解了你的问题.我认为问题更可能是结果对输入数字的敏感性.我们选择吧

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

解决方案是~1138和1361.现在,如果你计算相对导数.我可以使用ForwardDiff.jl包通过自动区分在Julia中执行此操作:

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%.

这里唯一的解决方案是重新设计您的问题,使其对输入值不那么敏感.