用修正函数确定n的Wilkinson多项式

Question

An example by Wilkinson[1963] shows that minute alterations in the coefficients of a polynomial may have massive effects on the roots. Let

f(x) = (x − 1)(x − 2) · · · (x − 20)

which has become known as the Wilkinson polynomial. Try to determine what happens to the zero r = 20 when the function is altered to f(x) − 10^(−8) * x^19 . Hint: use secant method in double precision

所以我想做的是首先f的零显然是整数1, 2, . . . , 20。但是我认为我需要使用割线方法来定位(20, 21)之间的零？我该怎么走？你知道吗

这是我的割线代码：

def SecantMethod(f,x0,x1,nmax,error):

    for i in range(nmax):

        if f(x1)-f(x0) == 0 or abs(f(x1)-f(x0)) <= error:
            return x1

        x2 = x1 - (f(x1)*(x1-x0))/(f(x1)-f(x0))
        x0 = x1
        x1 = x2
    return x1

Answer 1

试着解释你的公式，我得到以下解释

def f(x):
    res = 1;
    for k in range(20): res *= x-k-1
    return res

def f1(x): return f(x)-1e-8*x**19

x1, x2 = 20, 20.1
while abs(x2-x1)>1e-13: 
    print "x=%.16f,  f(x)=%.12g"%(x1, f1(x1)); 
    x1,x2 = x2, x2 - (f1(x2)*(x2-x1))/(f1(x2)-f1(x1))
print "x=%.16f,  f(x)=%.12g"%(x1, f1(x1)); 
print "x=%.16f,  f(x)=%.12g"%(x2, f1(x2)); 

迭代序列的结果输出

x=20.0000000000000000,  f(x)=-5.24288e+16
x=20.1000000000000014,  f(x)=-4.0426485012e+16
x=20.4368223967820448,  f(x)=1.4159426847e+17
x=20.1748076541532662,  f(x)=-2.31966692518e+16
x=20.2116899573180113,  f(x)=-1.12078103606e+16
x=20.2461694572670474,  f(x)=2.5395983862e+15
x=20.2397999600010081,  f(x)=-2.01198072116e+14
x=20.2402675359738211,  f(x)=-3.21814911103e+12
x=20.2402751363869271,  f(x)=4179299808
x=20.2402751265293332,  f(x)=-86592
x=20.2402751265295358,  f(x)=-776
x=20.2402751265295393,  f(x)=752

这是最好的实现作为情节

plt.plot(I,f1(20.24027512652954+1e-14*I));plt.grid(); plt.show()

演示在浮点精度中，-776和752之间没有中间值。你知道吗