Dekker 的数值分析方法收敛不正确

一些可能有帮助的更改：

1. 更改为整个代码。fnonfunc
2. 通过更改为 来更正收敛检查中的条件。abs(bk - b_minus_1)abs(bk_plus_1 - bk)

你有很多问题。

首先，您只是在使用二分法时进行计算。你需要在每次迭代时进行计算，因为你需要每次都进行比较。mkmkskmk

接下来，要前进和结尾的两行代码顺序不正确，因此在保存值之前将其丢弃。bkbk_plus_1bk

接下来，您正在执行计算以检查是否需要移动，但您从未使用过该计算的结果。这似乎有效：ak

import numpy as np

def dekker(func, x_bounds, tolerance, max_iterations):
    ak, bk = x_bounds
    b_minus_1 = ak # previous bk but start with a
    
    print(' k xk        f(xk)')

    for k in range(max_iterations):
        fk = func(bk)
        f_minus_1 = func(b_minus_1)

        # check for division by zero when attempting secant method

        mk = (ak + bk) / 2 # bisection 
        if fk == f_minus_1:
            bk_plus_1 = mk
        else:
            sk = bk - (bk - b_minus_1) / (fk - f_minus_1) * fk #secant
            
            if mk <= sk <= bk or bk <= sk <= mk:
                bk_plus_1 = sk
            else:
                bk_plus_1 = mk

        fk_plus_1 = func(bk_plus_1)

        b_minus_1 = bk

        if fk * fk_plus_1 < 0:
            ak = bk

        bk = bk_plus_1

        if abs(func(ak)) < abs(func(bk)):
            best_solution = ak
        else:
            best_solution = bk

        print('{0:2.0f}  {1:2.8f}  {2:2.2e}'.format(k, bk_plus_1, abs(fk_plus_1)), best_solution)

        if(abs(fk - fk_plus_1) < tolerance):
            print('Converged', bk)
            return bk

    print('Not converged')
    return None


def function(x):
    return -np.cos(x) + x*x*x + 2*x*x + 1

def plotme():
    import matplotlib.pyplot as plt
    x = np.arange( -3, 1, .01 )
    y = function(x)
    plt.plot(x,y)
    plt.show()

#plotme()
print(dekker(function, [-3,1], 1e-6, 100))

输出：

 k xk        f(xk)
 0  -0.32179374  2.25e-01 -0.3217937387377827
 1  -0.41378380  3.56e-01 -0.41378379538544446
 2  -1.70689190  1.99e+00 -1.7068918976927223
 3  -2.35344595  2.52e-01 -2.3534459488463613
 4  -2.28064072  1.92e-01 -2.280640719265942
 5  -2.31209175  6.87e-03 -2.3120917484530077
 6  -2.31325946  2.00e-04 -2.313259456426466
 7  -2.31322651  1.97e-07 -2.313226508522713
 8  -2.31322654  5.65e-12 -2.313226541057924
Converged -2.313226541057924
-2.313226541057924

这是一个工作版本。

def dekker(fnon, x_bounds, tolerance, max_iterations):
    # Initial bounds
    a0, b0 = x_bounds
    ak = a0
    bk = b0
    bk_minus_1 = a0 # previous bk but start with a
    
    k = 0

    print(' k xk        f(xk)')

    while k < max_iterations:
        fk = fnon(bk)
        fk_minus_1 = fnon(bk_minus_1)

        # check for division by zero when attemping secant method
        if fk - fk_minus_1 == 0:
            # implement bisection method when division by zero 
            bk_plus_1 = (a0 + bk) / 2
        else:
            sk = bk - (bk - bk_minus_1) / (fk - fk_minus_1) * fk #secant
            mk = (a0 + bk) / 2 # bisection 

            if sk >= mk and sk <= bk:
                bk_plus_1 = sk
            else:
                bk_plus_1 = mk

        fk_plus_1 = fnon(bk_plus_1)

        if fk * fk_plus_1 < 0:
            ak_plus_1 = bk
            bk_plus_1 = bk_plus_1
        else:
            ak_plus_1 = ak
            bk_plus_1 = bk_plus_1

        if abs(fnon(ak_plus_1)) < abs(fnon(bk_plus_1)):
            best_solution = ak_plus_1
        else:
            best_solution = bk_plus_1

        print('{0:2.0f}  {1:2.8f}  {2:2.2e}'.format(k, bk_plus_1, abs(fk_plus_1)))

        if(abs(bk - bk_plus_1) < tolerance):
            print('Converged')
            break

        bk_minus_1 = bk
        bk = bk_plus_1
        ak = ak_plus_1
        k += 1

    if k == max_iterations:
        print('Not converged')

和一些验证

def f2(x):
    return x*x - 4
dekker(f2, [1,4], 1e-6, 100)

指纹

 k xk        f(xk)
 1  2.50000000  2.25e+00
 2  2.15384615  6.39e-01
 3  2.01652893  6.64e-02
 4  2.00060976  2.44e-03
 5  2.00000251  1.00e-05
 6  2.00000000  1.53e-09
 7  2.00000000  1.78e-15
Converged

有许多问题以前没有提到，但一个重要的问题是在循环结束时重新初始化。如果你设置了 ，然后你实际上也设置了等于，并且丢失了你想设置为等于的值。您尚未从计算的 .bk = bk_plus_1bk_minus_1 = bkbk_minus_1bk_plus_1bkbk_minus_1akak_plus_1

您可以考虑进行一些进一步的改进。例如，您应该让该方法返回您找到的值。这样，函数的结果就可以用于进一步的计算。xk

来源

Dekker 的算法可能是什么有几个来源。我使用

• Cleve Moler （2015），*“Zeroin，第 1 部分，Dekker 算法

它使用来自

• T. J. Dekker 和 W. Hoffmann （1968），“数值代数中的 Algol 60 程序”

并引用了会议报告

• Dekker （1969），“通过连续线性插值找到零点”

混合了割线和二分步骤的寻根算法也作为算法 A 进行讨论。

• Bus， Dekker （1975）， “Two Efficient Algorithms with Guaranteed 求函数零点的收敛性”

该文章的主旨是讨论使用函数的非线性近似根作为迭代点对基本思想的修改。

维基百科的布伦特文章仅引用了会议报告Dekker（1969）和其他从布伦特算法角度进行的评论。因此，仅对于 Dekker 来说，这是一个可疑的来源，因为这两个来源都可能将算法简化为其基本思想。

算法

在平分法和其他包围法的扩展中，戴克斯法的状态有 3 个点 a、b、c。它们的功能是从 a，b 计算下一个割线根，并且 [b，c] 是一个括号区间，也就是说，b，c 处的值具有不同的符号，确保区间中存在根。此外，b是“最佳”点，在三个点中具有最小的函数值。该算法强制要求新点位于 [b，c] 内，更接近 b。

在方法开始时，给出区间 [a，b] 并将其扩展为“三角形” [a，b，c] = [a，b，a]。对点进行排序，使 b 的值较小。选择 [b，c] 内的下一个点 v 后，割线序列移动到 * a，b = b，v*。然后选择新的包围间隔。新 a、b、c 处值的符号序列有一个交替。如果它再次位于 b 和 c 之间，则无需执行任何操作，序列 a，b，c 保持“线性”。否则交替在 a 和 b 之间，因此将 c=a 设置为新的对位点。

中点选择有 3 个变体：

• 在“割线”模式下，割线根 v=s 位于括号间隔内，而且位于中点 m 和 b 之间。在大多数情况下，s 处的值与 b 处的值具有相同的符号，并且大小明显更小，因此另一个区间端 c 保持不变，割线根序列可以不间断地继续，a，b = b，s。 b 和 v=s 处的值之间的符号变化不太可能，但仍然很有可能。

• 在“中点”模式下，检测到割线根序列偏离包围间隔，要么落在外面，要么从b到m的另一侧意外地走了一大步。在这种情况下，更有可能的是中点的值不小于以前的值，或者符号在 b 和 v=m 的值之间变化

• 执行“最小”步骤来解决数值不明确的情况。例如，如果割线根序列收敛到低于容差的步长，则使用一个或多个最小步长将包围间隔也缩小到此最小大小。在没有符号变化的多根处收敛时，这在数值上是模糊的，一个或多个最小步骤会将此根移出括号间隔，并重新开始搜索具有符号变化的根。

应用于测试问题

此方程的 Dekker 方法的迹线应如下所示

a=x[-1], b=x[-1], c=x[-1] -> initial   : x[ 0] =  -3.000000000000000,  f(x[ 0]) =  -7.010007503399553
a=x[ 0], b=x[-1], c=x[-1] -> initial   : x[ 1] =   1.000000000000000,  f(x[ 1]) =   3.459697694131860
a=x[ 0], b=x[ 1], c=x[ 0] -> secant    : x[ 2] =  -0.321793738737783,  f(x[ 2]) =   0.225110648393170
a=x[ 1], b=x[ 2], c=x[ 0] -> secant    : x[ 3] =  -0.413783795385445,  f(x[ 3]) =   0.355981221390157
a=x[ 2], b=x[ 3], c=x[ 0] -> midpoint  : x[ 4] =  -1.706891897692722,  f(x[ 4]) =   1.989640412053239
a=x[ 3], b=x[ 4], c=x[ 0] -> midpoint  : x[ 5] =  -2.353445948846361,  f(x[ 5]) =  -0.252473245321568
a=x[ 4], b=x[ 5], c=x[ 4] -> secant    : x[ 6] =  -2.280640719265942,  f(x[ 6]) =   0.192012983736426
a=x[ 5], b=x[ 6], c=x[ 5] -> secant    : x[ 7] =  -2.312091748453008,  f(x[ 7]) =   0.006873812772492
a=x[ 6], b=x[ 7], c=x[ 5] -> secant    : x[ 8] =  -2.313259456426466,  f(x[ 8]) =  -0.000199582032968
a=x[ 7], b=x[ 8], c=x[ 7] -> secant    : x[ 9] =  -2.313226508522713,  f(x[ 9]) =   0.000000197276950
a=x[ 8], b=x[ 9], c=x[ 8] -> secant    : x[10] =  -2.313226541057924,  f(x[10]) =   0.000000000005652
a=x[ 9], b=x[10], c=x[ 8] -> minimal   : x[11] =  -2.313226541058924,  f(x[11]) =  -0.000000000000412
exit by interval length
a=x[10], b=x[11], c=x[10] -> best value: x[11] =  -2.313226541058924,  f(x[11]) =  -0.000000000000412

的使用取自 Dekker/Bus 论文。这里是与当前最佳点值相反的对位点，如果最后一步是割线步，则通常是前一个最佳点。a,b,ccba

请注意函数值中的常规符号开关，它表示包围间隔的大幅缩短。[b,c]

---

该代码具有打印出算法轨迹的所有功能，是

from math import sin, cos, exp

def dekker(f_in, a, b, tol=1e-7):
    mode = "initial"
    A=B=C=None
    
    class Point(object):
        k=0
        def __init__(self, x):
            self.x = x
            self.k = Point.k
            self.f = f_in(x); 
            a_str = "       " if not A else f"a=x[{A.k:2d}]"
            b_str = "       " if not B else f"b=x[{B.k:2d}]"
            c_str = "       " if not C else f"c=x[{C.k:2d}]"
            print(f"{a_str}, {b_str}, {c_str} -> {mode:10s}:", end='')
            print(f" x[{self.k:2d}] = {self.x:19.15f},  f(x[{self.k:2d}]) = {self.f:19.15f}"); 
            Point.k+=1; 
 
    A = Point(a)
    B = Point(b)
    if A.f*B.f > 0: 
        error("need sign change in initial interval")
    # [a,b] is bracketing interval, establish triangle formation, [b,c] now also bracketing
    C = A
    while 1>0:
        # make b the best. In superlinear mode nothing should happen here
        if abs(C.f) < abs(B.f): 
            A, B, C = B, C, B    # triangle formation
        # atomic step size
        db = tol
        if abs(C.x-B.x)<2*db: print("exit by interval length"); break
        # minimal change in case of stalling
        h = B.x + (db if C.x>B.x else -db);
        # elements of the secant step
        p = B.f*(A.x-B.x); q = (B.f-A.f);
        if p<0: p,q = -p,-q
        # default for midpoint step
        v = m = 0.5*(B.x+C.x); mode = "midpoint"
        # secant step if step update is large enough
        if p > abs(db*q):
            s = B.x+p/q
            # step goes not beyond midpoint, else remain at midpoint
            if p < q*(v-B.x): v = s; mode = "secant"
        else: v=h; mode = "minimal"
        # shift the secant sequence
        A, B = B, Point(v);
        # Now a < b=v <= m < c or in the other direction. Between a and c there is a sign change
        # If the sign also changes between b and c, do nothing, linear formation remains.
        # Else copy a to c, establishing a triangle formation
        if B.f*C.f > 0:
            C = A
        # Now [b,c] is ensured to be bracketing and (a,b) is ready for the next secant root computation
    mode = "best value"
    Point.k = B.k; Point(B.x);    # this is just for some nice printing
    return B.x, C.x
        
dekker(lambda x: -cos(x) + x*x*x + 2*x*x + 1, -3, 1, tol=1e-12)