当我更改变量 `c` 的值时，如何测试机器精度的错误？

Question

我正在关注 Landau 的计算物理概览并且我正在做这个问题（b 部分）：

我已经编写了生成以下 4 个解决方案的代码。我对作者要我在这里做什么感到困惑。

在测试减法抵消错误时，是否就像结果 2 减去结果 1 一样简单？关键是要看看机器是如何由于舍入而搞砸减法的，所以我无法从机器那里得到准确的答案，对吧？
对于我的解决方案，我发现使用 c=10^-1 会导致方法 1 的 1.38777878078145e-17 和方法 2 的 1.11022302462516e-16 。当我启动它时到 c=10^-16（结果没有得到 inf），我得到方法 1 的 -1.10223024625156e-17 和方法 2 的 -0.0992800745259007。这告诉我，对于方法 2，有一个非常当我修改 c 时，值发生了很大的变化。但是我觉得我没有抓住作者想让我在这里理解的本质。

如果你想看代码，就在这里。否则，您可以在下面看到一些示例输出：

import Foundation

var a: Double = 1.0
var b: Double = 1.0
var c: Double = pow(10,-16)
var x: Double = 0.0

//The following variables are just used to run the calculations within 
//the functions, i.e. they can be reused within each function definition
var x1: Double = 0.0
var x2: Double = 0.0

func x_func(a_var: Double, b_var: Double, c_var: Double) -> Double {

    x1 = 0.0
    x2 = 0.0
    x = 0.0

    x1 = pow(b,2.0)-4.0*a*c
    x2 = -b + x1.squareRoot()
    x = x2/(2.0*a)  

    return x

}
func x_func2(a_var: Double, b_var: Double, c_var: Double) -> Double {

    x1 = 0.0
    x2 = 0.0
    x = 0.0

    x1 = pow(b,2.0)-4.0*a*c
    x2 = -b - x1.squareRoot()
    x = x2/(2.0*a)  

    return x

}
func xp_func(a_var: Double, b_var: Double, c_var: Double) -> Double {

    x1 = 0.0
    x2 = 0.0
    x = 0.0

    x1 = pow(b,2.0)-4.0*a*c
    x2 = x1.squareRoot() + b
    x = -2.0*c/x2

    return x

}
func xp_func2(a_var: Double, b_var: Double, c_var: Double) -> Double {

    x1 = 0.0
    x2 = 0.0
    x = 0.0

    x1 = pow(b,2.0)-4.0*a*c
    x2 = b - x1.squareRoot()
    x = -2.0*c/x2

    return x

}

print("Method 1: Positive")
print(x_func(a_var: a, b_var: b, c_var: c))
print(" \n")



print("Method 1: Negative")
print(x_func2(a_var: a, b_var: b, c_var: c))
print(" \n")



print("Method 2: Positive")
print(xp_func(a_var: a, b_var: b, c_var: c))
print(" \n")



print("Method 2: Negative")
print(xp_func2(a_var: a, b_var: b, c_var: c))

print(" \n")

print("Subtractive cancellation error for Method 1:")
print(" \n")
print("Method 1 (positive) minus Method 2 (positive)",x_func(a_var: a, b_var: b, c_var: c) - xp_func(a_var: a, b_var: b, c_var: c))

print(" \n")
print("Subtractive cancellation error for Method 2:")
print(" \n")
print("Method 1 (negative) minus Method 2 (negative)",x_func2(a_var: a, b_var: b, c_var: c) - xp_func2(a_var: a, b_var: b, c_var: c))

示例输出：

Method 1: Positive
-1.11022302462516e-16

Method 1: Negative
-1.0

Method 2: Positive
-1e-16

Method 2: Negative
-0.900719925474099

Subtractive cancellation error for Method 1:

Method 1 (positive) minus Method 2 (positive) -1.10223024625156e-17

Subtractive cancellation error for Method 2:

Method 1 (negative) minus Method 2 (negative) -0.0992800745259007

Answer 1

首先让我们简化您的代码。你有两种计算方法二次方程的解 a x^2 + b x + c == 0:

func qsolve1(a: Double, b: Double, c: Double) -> (Double, Double) {
    let x1 = (-b - (b*b - 4*a*c).squareRoot())/(2*a)
    let x2 = (-b + (b*b - 4*a*c).squareRoot())/(2*a)
    return (x1, x2)
}

func qsolve2(a: Double, b: Double, c: Double) -> (Double, Double) {
    let x1 = -2*c/(b - (b*b - 4*a*c).squareRoot())
    let x2 = -2*c/(b + (b*b - 4*a*c).squareRoot())
    return (x1, x2)
}

（这里只考虑a != 0、c != 0和b^2 - 4ac >= 0的情况。）

如果计算完成，两种方法都会return相同的结果完全正确，但由于浮点舍入误差，它们是不同的。

应用于多项式

x^2 + x + 1.0E-10

你得到结果：

print(qsolve1(a: 1.0, b: 1.0, c: 1.0E-10))
// (-0.99999999989999999, -1.000000082740371e-10)

print(qsolve2(a: 1.0, b: 1.0, c: 1.0E-10))
// (-0.99999991725963588, -1.0000000001000001e-10)

如果我们将这些数字与 PARI/GP

的结果进行比较

? solve(x=-1.5, -0.5, x^2+x+1.0e-10)
%2 = -0.99999999989999999998999999999800000000
? solve(x=-0.5, +0.5, x^2+x+1.0e-10)
%3 = -1.0000000001000000000200000000050000000 E-10

然后我们看到 qsolve1() 计算了 -1.0 附近的解更准确地说，qsolve2() 计算出 0.0 附近的解更准确地说。

为什么会这样？让我们检查一下中间结果：

print("b=", b.debugDescription, " sqrt(b^2-4ac)=", (b*b - 4*a*c).squareRoot().debugDescription)
// b= 1.0  sqrt(b^2-4ac)= 0.99999999979999998

二进制浮点数的精度有限，大约 IEEE 754 双精度数的 16 位十进制数字。

计算qsolve1()中的x1或qsolve2()中的x2时，那些两个数相减，差值约为

0.00000000020000002

只有 8 位有效数字，不再是 16 位。那就是 "subtractive cancellation" 当两个数字大约相同的值相互减去。

因此，加上数字在数值上更好。无论 x1 或 x2 更精确取决于 b 的符号。在任何情况下，另一个解决方案可以从关系计算 x1 * x2 = c 不损失精度。

给出以下方法：

func qsolve(a: Double, b: Double, c: Double) -> (Double, Double) {
    let x1: Double
    if b >= 0 {
        x1 = (-b - (b*b - 4*a*c).squareRoot())/(2*a)
    } else {
        x1 = (-b + (b*b - 4*a*c).squareRoot())/(2*a)
    }
    let x2 = c/x1
    return (x1, x2)
}

对于我们的测试用例，我们得到

print(qsolve(a: 1.0, b: 1.0, c: 1.0E-10))
// (-0.99999999989999999, -1.0000000001000001e-10)

当我更改变量 `c` 的值时，如何测试机器精度的错误？

How can I test for the error in machine precision as I change the value of the variable `c`?

precision

rounding-error

subtraction

swift