Как я могу проверить ошибку в точности машины, поскольку я изменяю значение переменной `c`?

Question

Я следую опрос Ландау вычислительной физики и я делаю эту проблему (часть б):

Я написал код для создания 4 решения ниже. Я смущен тем, что автор хочет, чтобы я сделал здесь.

1. При тестировании на ошибки сокращения вычитания это так же просто, как вычесть результат 1 из результата 2? Дело в том, чтобы увидеть, как машина затягивает вычитание из-за округления, поэтому я не могу получить точный ответ от машины, не так ли?

2. Для моего решения это, я обнаружил, что использование c=10^-1 результатов в 1.38777878078145e-17 для метода 1 и 1.11022302462516e-16 для метода 2. Когда я провернуть его до c=10^-16 (без получения inf в результате), я получаю -1.10223024625156e-17 для метода 1 и -0.0992800745259007 для метода 2. Это говорит о том, что для метода 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: (. Мы только рассматриваем случай a != 0, c != 0 и b^2 - 4ac >= 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) 
} 

Оба метода будет возвращать те же результаты, если вычисления были выполнены точно, но из-за ошибок округления с плавающей запятой они различны.

Применительно к полиномиальному

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 

Двоичное число с плавающей точкой имеет ограниченную точность, около 16 десятичных цифр для IEEE 754 двойной точности.

При вычислении x1 в qsolve1() или x2 в qsolve2(), эти два числа вычитали, и разница составляет около

0.00000000020000002 

, который имеет только 8 значащих цифр, а не 16 больше. То есть «субтрактивное аннулирование», и это происходит, когда два числа примерно одинаковое значение вычитаются друг из друга.

Поэтому численно лучше добавить цифры. 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)