Horner's Rule, and plotting rounding errors

Question

I have been trying to understand the plot seen in Fig.2.3. I hope someone can explain it to me because I thought attempting to replicate it as instructed in Q12 would help me understand it. However I cannot seem to replicate this no matter what.

The text along with this image is as follows (please excuse images)

With the actual question being: "Apply Horner's rule twice (once for numerator and once for the denominator) to produce two plots, one for x = 1.606+2^{-52}i and one for x = 2.400+2^{-52}i. Your results should look like Fig.2.3".

r(x) = (4x^4 - 59x^3 + 324x^2 - 751x + 622) / (x^4 - 14x^3 + 72x^2 - 151x + 112)

so far my code is:

def horner(coeffs, x):
    result = 0
    for coeff in reversed(coeffs):
        result = result * x + coeff
    return result
numerator_coeffs = [4, -59, 324, -751, 622]
denominator_coeffs = [1, -14, 72, -151, 112]
def rational_function(x):
    num = horner(numerator_coeffs, x)
    denom = horner(denominator_coeffs, x)
    return num / denom

From here i have no idea how to plot the function. I assume it is obviously some scatter plot but I dont know what to use as the x parameter? Is it even right to use 1.606 in my plot like so..? xs = [1.606+(2**-52)*i for i in range(801)] I figure that the x axis is in fact integer multiples of i as also mentioned as a comment. However in my efforts of trying to graph this i get just a straight line as a result.

xs = [1.606+(2**-52)*i for i in range(801)]
r = []
for x in xs:
    r.append(rational_function(x))
plt.scatter(xs,r)
plt.show()

As seen here in my plot attempt of Fig.2.3

This seems close but not sure what I am missing.
Appreciate any help possible. Thanks!

Answer 1

First, your implementation of Horner's method can be improved

In [33]: def horner(coeffs, x):
    ...:     result = 0
    ...:     for coeff in coeffs[::-1]: # with reversal
    ...:         result = result * x + coeff
    ...:     return result
    ...: numerator_coeffs = [4, -59, 324, -751, 622]
    ...: x = 1.606
    ...: print(horner(numerator_coeffs, x))
    ...: print(4*x**4 - 59*x**3 + 324*x**2 - 751*x + 622)
1771.9155387629123
33.78336943078409

In [34]: def horner(coeffs, x):
    ...:     result = 0
    ...:     for coeff in coeffs: # no reversal
    ...:         result = result * x + coeff
    ...:     return result
    ...: numerator_coeffs = [4, -59, 324, -751, 622]
    ...: x = 1.606
    ...: print(horner(numerator_coeffs, x))
    ...: print(4*x**4 - 59*x**3 + 324*x**2 - 751*x + 622)
33.78336943078398
33.78336943078409

In your plot, the range of y values is (5.6-4.4)10^-11 = 1.2·10^-11 while in the original plot it is (5.507-4.245)10^-13 = 1.262·10^-13 so you cannot see the dispersion that it's present because, in your plot, the points are too close together.

The remedy is, you manually set the limits of the y axis to obviate the poor choice Matplotlib has made.

import numpy as np
import matplotlib.pyplot as plt

def horner(coeffs, x):
    result = 0
    for coeff in coeffs:
        result = result * x + coeff
    return result

num_c = [4, -59, 324, -751, 622]
den_c = [1, -14,  72, -151, 112]
dx = pow(2, -52)
x = np.array([1.606 + i*dx for i in range(801)])
ratio = horner(num_c, x)/horner(den_c, x)
dr = max(ratio)-min(ratio)
##################################################
plt.ylim((min(ratio)-0.05*dr, max(ratio)+0.05*dr))
##################################################
plt.scatter(range(801), ratio, s=1)
plt.show()

---

Why x=1.606 and x=2.400?

Probably because we have a stationary point, probably a very flat stationary point (higher derivatives close to zero) in these two points, so that the r(x) value in a neighborhood of the stationary points is influenced more from numerical errors than from some linear trend...