钍。J. Dekker零点1969年的算法是我最喜欢的算法之一。一种优雅的技术，结合了一分的分配和SECANT方法来查找真实变量函数的零，它已成为Fzero.今天在Matlab。这是三部分系列中的第一个。

内容

Dirk Dekker.

我刚刚来自于Fortieth Woudschoten数值分析会议由荷兰/弗拉芒举办Werkgemeenschap科学计算团体。其中一个特殊的客人是th。J.“Dirk”Dekker，阿姆斯特丹大学的数学和计算机科学的退休教授。

1968年，Dekker与同事沃尔特霍夫曼一起发布了来自阿姆斯特丹的数学中心的两零零份报告，描述了一个综合的AlGOL图书馆，用于矩阵计算。本报告前三年的威尔金森和雷恩斯的出版物自动计算手册，线性代数这为eispack形成了基础，最终产生了matlab。

ZHILIN IN ALGOL.

DEKKER用于计算实际对称三角形矩阵的几个特征值的程序涉及Sturm序列并呼叫实用程序，零点。这是一个扫描零点从Dekker的1968年AlGOL报告中获取。这是我即将描述的MATLAB代码的起点。

Dekker提出零点在会议上，由B. dejon和Peter Henrici编辑了诉讼程序。Jim Wilkinson在1967年的斯坦福报表中描述了一种类似的算法。斯坦福博士。D.学生理查德布伦特提出了重要的改进，我将在我的下一个博客文章中描述。Forsythe，Malcolm和我制作了Fortran Zero Finder的基础对于数学计算的计算机方法。它随后演变为Fzero.在matlab。

测试功能

这是我将用来演示的功能零点。

$$ f（x）= \ frac {1} {x-3} -6 $$

f = @（x）1./( x-3 )-6;Ezplot（f，[3,4]）持有在情节（3 + 1 / 6,0，'ro'）设置（GCA，'xtick'，[3：1/6：4]）保持离开

此功能是有意棘手。有一个极点$ X = 3 $，我们正在努力寻找零附近是在$ X = 3 \压裂{1} {6} $。其中函数是正的实线的唯一部分是左侧的六分之一的上述$ X $这两个点之间轴摆动。

对分

DEKKER的算法需要从给定函数$ f（x）$更改符号的间隔$ [a，b] $开始。在浮点算术中，真实变量的连续功能的概念变得有点难以捉摸，因此我们将我们的目标设置为找到一个更小的子因素，其中功能更改符号。我通过消除指定收敛间隔长度的公差来简化呼叫序列。以下所有功能迭代直到终端点$ B $的间隔数为$ [A，B] $的长度为圆形循环错误。

可靠，故障安全部分零点是二分算法。这个想法是重复削减一半的时间间隔[a，b] $，同时继续跨越标志。不使用实际函数值，只使用标志。这是自身分娩的代码。

类型一分钟

函数b =平分（F，A，B）S = '％5.0F％19.15f％19.15f \ n';FPRINTF（S，1，A，F（A））FPRINTF（S，2，B，F（b））k = 2;而ABS（B-A）> EPS（ABS（B））X =（A + B）/ 2;如果标志（F（X））==标志（F（B））B = X;否则A = X;端k = k + 1;fprintf中（S，K，X，F（X））结束端

让我们看看如何一分钟执行我们的测试功能。间隔[3,4]提供令人满意的起始间隔，因为IEEE浮点算术产生正确签名的无限远，+ INF，在杆子。

二分（F，3,4）

1 3.000000000000000 Inf文件2 4.000000000000000 -5.000000000000000 3 3.500000000000000 -4.000000000000000 4 3.250000000000000 -2.000000000000000 5 3.125000000000000 2.000000000000000 6 3.187500000000000 -0.666666666666667 7 3.156250000000000 0.400000000000000 8 3.171875000000000 -0.181818181818182 9 3.164062500000000 0.095238095238095 10 3.167968750000000 -0.046511627906977 11 3.166015625000000 0.023529411764706 12 3.166992187500000 -0.011695906432749 13 3.166503906250000 0.005865102639296 14 3.166748046875000 -0.00292825768667615 3.166625976562500 0.001465201465201 16 3.166687011718750 -0.000732332478946 17 3.166656494140625 0.000366233290606 18 3.166671752929688 -0.000183099880985 19 3.166664123535156 0.000091554131380 20 3.166667938232422 -0.000045776017944 21 3.166666030883789 0.000022888270905 22 3.166666984558106 -0.000011444069969 23 3.166666507720947 0.000005722051355 24 3.166666746139526 -0.000002861021585 25 3.166666626930237 0.000001430511816 26 3.166666686534882 -0.000000715255.6.5.227.3.。16.6.6.6.6.656732559 0.000000357627890 28 3.166666671633720 -0.000000178813929 29 3.166666664183140 0.000000089406969 30 3.166666667908430 -0.000000044703484 31 3.166666666045785 0.000000022351742 32 3.166666666977108 -0.000000011175871 33 3.166666666511446 0.000000005587935 34 3.166666666744277 -0.000000002793968 35 3.166666666627862 0.000000001396984 36 3.166666666686069 -0.000000000698492 37 3.166666666656965 0.000000000349246 38 3.166666666671517 -0.000000000174623 39 3.166666666664241 0.000000000087311 40 3.166666666667879 -0.000000000043656 41 3.166666666666060 0.000000000021828 42 3.166666666666970 -0.000000000010914 43 3.166666666666515 0.000000000005457 44 3.166666666666743 -0.000000000002728 45 3.166666666666629 0.000000000001364 46 3.166666666666686 -0.000000000000682 47 3.166666666666657 0.000000000000341 48 3.166666666666671 -0.000000000000171 49 3.166666666666664 0.000000000000085 50 3.166666666666668 -0.000000000000043 51 3.166666666666666 0.000000000000021 52 3.166666666666667 -0.000000000000011 53 3.166666666666667 0.000000000000005 ans = 3.166666666666667

您可以看到间隔的长度在早期阶段的一半被切割成一半X值，然后是值f（x）在后续阶段切成两半。这需要53次迭代。事实上，我的停止标准涉及EPS（b）中，如果开始一种和B.大小相当，一分钟需要53次迭代任何功能是因为双精度浮点数有52位。

除了在我的下一篇文章中描述的终点附近的某些病理行为案例，保证了小分割融合。如果您有一个具有符号更改的开始间隔，则此算法几乎肯定会找到一个带有符号更改的微小子宫。但它太慢了。53迭代通常太多了。我们需要一些少步的东西。

割割方法

如果我们有两个迭代$ a $和$ b $，则具有相应的函数值，他们是否表现出签署变化，我们可以将下一个迭代拿到直线通过$ [a，f（a）] $和$ [b，f（b）] $与$ x $ -axis相交。这是SECANT方法。这是代码，通过计算致密地完成，以避免不必要的下溢或溢出。麻烦是没有什么可以阻止下一个点在间隔之外。

类型割线

功能b = secant（f，a，b）s ='％5.0f％19.15f％23.15e \ n';FPRINTF（S，1，A，F（A））FPRINTF（S，2，B，F（b））k = 2;ABS（B-A）> EPS（ABS（B））C = a;a = b;B = B +（B  -  C）/（F（c）/ f（b） -  1）;k = k + 1;fprintf（s，k，b，f（b））结束

透过无穷大的点没有意义，所以我们必须从稍微较小的间隔开始，但即使这很快就会陷入困境。

割线（F，3.1,3.5）

1 3.100000000000000 3.999999999999991e + 00 2 3.500000000000000 -4.000000000000000e + 00 3 3.300000000000000 -2.666666666666665e + 00 4 2.900000000000000 -1.600000000000004e + 01 5 3.379999999999999 -3.368421052631575e + 00 6 3.507999999999999 -4.031496062992121e + 00 7 2.729760000000002 -9.700414446418032e + 00 8 4.061451519999991-5.057893854634069e + 00 9 5.512291472588765 -5.601957013781701e + 00 10 -9.426310620985419 -6.080474408736508e + 00 11 180.397066104648590 -5.994362928192905e + 00 12 13394.317035490772000 -5.999925324746076e + 00 13 -14239910.406137096000000 -6.000000070225146e + 00 14 1144132943372.700400000000000 -5.999999999999126e + 00 1597754126073011126000.000000000000000 -6.000000000000000e + 00 16 -671105854578598490000000000000000.000000000000000 -6.000000000000000e + 00 17 -Inf -6.000000000000000e + 00 ANS = -Inf

这说明割线方法是不可靠的。在步骤4中通过先前的两个点的割线符合$ X - 轴$为$ X = 2.9 $这已经是初始区间之外。通过该步骤10中的迭代已经跃升到功能的其他分支。在步骤12中的计算值超过我的输出格式。

有趣的是，如果我们扭转了$ a $和$ b $的角色，割线没有逃脱间隔。

割线（F，3.5,3.1）

1 3.500000000000000 -4.000000000000000e + 00 2 3.100000000000000 3.999999999999991e + 00 3 3.300000000000000 -2.666666666666665e + 00 4 3.220000000000000 -1.454545454545450e + 00 5 3.124000000000000 2.064516129032251e + 00 6 3.180320000000000 -4.543034605146419e-01 7 3.170161920000000 -1.232444955957233e-01 8 3.166380335513600 1.032566086067455E-02 9 3.166672671466170 -2.161649939527166e-04 10 3.166666676982834 -3.713819847206423e-07 11 3.166666666666295 1.338662514172029e-11 12 3.166666666666667 5.329070518200751e-15 13 3.166666666666667 5.329070518200751e-15 ANS = 3.166666666666667

Zeroin算法

Dekker的算法保持三个变量，$ a $，$ b $和$ c $：最初，$ c $设置等于$ a $。

在每次迭代中，选择有可能由三个步骤做：

托in在matlab.

类型零点

功能B = Zeroin（F，A，B）％零点（F，A，B）将间隔[A，B]还原为Tiny％Interal，其中函数f（x）改变符号。％返回该间隔的一个结束点。fa = f（a）;fb = f（b）;如果符号（fa）== sign（fb）错误（'抱歉，f（x）必须更改间隔[a，b]。'）结束s ='％5.0f％8s％19.15f％23.15e \n';FPRINTF（S，1，'初始'，a，fa）fprintf（s，2，'初始'，b，fb）k = 2;％A是B和[B，C]的先前值始终包含零。c = a;FC = FA;虽然真正的IF符号（fb）==符号（fc）c = a; fc = fa; end % Swap to insure f(b) is the smallest value so far. if abs(fc) < abs(fb) a = b; fa = fb; b = c; fb = fc; c = a; fc = fa; end % Midpoint. m = (b + c)/2; if abs(m - b) <= eps(abs(b)) return % Exit from the loop and the function here. end % p/q is the the secant step. p = (b - a)*fb; if p >= 0 q = fa - fb; else q = fb - fa; p = -p; end % Save this point. a = b; fa = fb; k = k+1; % Choose next point. if p <= eps(q) % Minimal step. b = b + sign(c-b)*eps(b); fb = f(b); fprintf(s,k,'minimal',b,fb) elseif p <= (m - b)*q % Secant. b = b + p/q; fb = f(b); fprintf(s,k,'secant ',b,fb) else % Bisection. b = m; fb = f(b); fprintf(s,k,'bisect ',b,fb) end end end

这里是我们测试功能的性能。

零点（f，3​​,4）

1初始3.000000000000000 Inf文件2初始4.000000000000000 -5.000000000000000e + 00 3正割4.000000000000000 -5.000000000000000e + 00 4最小3.999999999999999 -4.999999999999999e + 00 5平分3.500000000000000 -3.999999999999998e + 00 6平分3.250000000000000 -2.000000000000000e + 00 7平分3.125000000000000 2.000000000000000e +0.0.8.割线3.。18.7.5.0.0.000000000 -6.666666666666670e-01 9 secant 3.171875000000000 -1.818181818181817e-01 10 secant 3.166015625000000 2.352941176470580e-02 11 secant 3.166687011718750 -7.323324789458852e-04 12 secant 3.166666746139526 -2.861021584976697e-06 13 secant 3.166666666656965 3.492459654808044e-10 14 secant 3.166666666666667 5.329070518200751e-15 15 minimal 3.166666666666667 -1.065814103640150e-14 ans = 3.166666666666667

最小的一步甚至有助于让东西远离杆子。几个小分阶段将间隔降至下降到

$$ [3 \ frac {1} {8}，3 \ frac {1} {4}] $$

然后，割线可以安全地接管并以六个步骤获得零。

参考

T. J. Dekker和W. Hoffmann（第2部分），Algol 60个是程序在数值代数那第1部分和2，炸药22和23，Mathematisch Centrum Amsterdam，1968年。

T. J. DEKKER，“通过连续的线性插值找到零”，在B. dejon和P. henrici（EDS）中，代数基本定理的建设性方面，纽约，1969年赛道。

获取matlab代码

发布与Matlab®R2015A