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

Question

I'm following Landau's A Survey of Computational Physics and I am doing this problem (part b):

I've written code to generate the 4 solutions below. I'm confused as to what the author wants me to do here.

1. When testing for subtractive cancellation errors, is it as simple as subtracting result 1 from result 2? The point is to see how a machine screws up subtraction due to rounding, so there is no way I can get the exact answer from a machine, right?

2. For my solution to this, I found that using c=10^-1 results in 1.38777878078145e-17 for method 1 and 1.11022302462516e-16 for method 2. When I crank it up to c=10^-16 (without getting inf as a result), I get -1.10223024625156e-17 for method 1 and -0.0992800745259007 for method 2. This tells me that for method 2 there was a very large change in value as I modified c. But I don't think I'm grasping the essence of what the author wants me to understand here.

If you want to see the code, here it is. Otherwise, you can see some sample output below:

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))

Sample output:

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

First let us simplify your code. You have two methods to compute the solutions of the quadratic equation 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)
}

(We are only considering the case a != 0, c != 0 and b^2 - 4ac >= 0 here.)

Both methods would return the same results if the computations were done exactly, but due to the floating point rounding errors they are different.

Applied to the polynomial

x^2 + x + 1.0E-10

you get the results:

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)

If we compare those numbers with the results from 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

then we see that qsolve1() computes the solution near -1.0 more precisely, and qsolve2() computes the solution near 0.0 more precisely.

Why does that happen? Let us inspect the intermediate results:

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

A binary floating point number has a limited precision, about 16 decimal digits for IEEE 754 double precision numbers.

When computing x1 in qsolve1() or x2 in qsolve2(), those two numbers are subtracted, and the difference is approximately

0.00000000020000002

which has only 8 significant digits, and not 16 anymore. That is the "subtractive cancellation" and it happens when two numbers of about the same value are subtracted from each other.

It is therefore numerically better to add the numbers. Whether x1 or x2 are more precise depends on the sign of b. In any case, the other solution can be computes from the relation x1 * x2 = c without loss of precision.

That gives the following method:

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)
}

and for our test case we get

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