霍纳规则，并绘制舍入误差

Question

我一直在尝试理解图 2.3 中的情节。我希望有人可以向我解释它，因为我认为尝试按照问题 12 中的指示复制它会帮助我理解它。然而无论如何我似乎都无法复制这个。

图片文字如下（请原谅图片）

实际问题是：“应用霍纳规则两次（一次用于分子，一次用于分母）以生成两个图，一个用于 x = 1.606+2^{-52}i，另一个用于 x = 2.400+2^ {-52}i. 您的结果应如图 2.3".

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

到目前为止我的代码是：

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

从这里我不知道如何绘制函数。我认为这显然是一些散点图，但我不知道使用什么作为 x 参数？像这样在我的情节中使用 1.606 是否正确？ xs = [1.606+(2**-52)*i 对于范围 (801) 内的 i] 我认为 x 轴实际上是 i 的整数倍，正如评论中提到的那样。然而，在我尝试绘制图表的过程中，我得到的结果只是一条直线。

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

正如我在图 2.3 的绘图尝试中所见

这看起来很接近，但不确定我错过了什么。
感谢任何可能的帮助。谢谢！

Answer 1

首先，你对霍纳方法的实现可以改进

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

在您的图中，y值的范围是 (5.6-4.4)10^-11 = 1.2·10^-11，而在原始图中它是 (5.507-4.245)10^-13 = 1.262 ·10^-13，所以你看不到它存在的色散，但在你的图中，这些点太靠近了。补救措施是，您手动设置 y 轴的限制，以避免 Matplotlib 做出的错误选择。

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