{"id":1269,"date":"2015-10-12T12:00:35","date_gmt":"2015-10-12T17:00:35","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=1269"},"modified":"2016-07-05T13:51:54","modified_gmt":"2016-07-05T18:51:54","slug":"zeroin-part-1-dekkers-algorithm","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2015\/10\/12\/zeroin-part-1-dekkers-algorithm\/","title":{"rendered":"Zeroin, Part 1: Dekker&#8217;s Algorithm"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>Th. J. Dekker's <i>zeroin<\/i> algorithm from 1969 is one of my favorite algorithms. An elegant technique combining bisection and the secant method for finding a zero of a function of a real variable, it has become <tt>fzero<\/tt> in MATLAB today.  This is the first of a three part series.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#a6ecbf14-02ae-491e-b956-18af08e52e9d\">Dirk Dekker<\/a><\/li><li><a href=\"#466bc11f-42d0-4edc-9cf9-9fe9049d5d3e\">Zeroin in Algol<\/a><\/li><li><a href=\"#5cbc4ee9-2180-4542-8aa6-2c9cbda8e02e\">The test function<\/a><\/li><li><a href=\"#e996c443-bc81-4bd3-82a2-5e7a7016e989\">Bisection<\/a><\/li><li><a href=\"#7e654aa5-e93a-4189-9712-06a4b9bc42c1\">Secant method<\/a><\/li><li><a href=\"#cf04caae-7244-461a-80a6-446ecb6d3e21\">Zeroin algorithm<\/a><\/li><li><a href=\"#6977f5dc-1980-4a0d-bb63-f67c61d81b93\">Zeroin in MATLAB<\/a><\/li><li><a href=\"#1fb0d597-b7a3-4e61-9584-3e263b3e2070\">References<\/a><\/li><\/ul><\/div><h4>Dirk Dekker<a name=\"a6ecbf14-02ae-491e-b956-18af08e52e9d\"><\/a><\/h4><p>I have just come from the <a href=\"https:\/\/wsc.project.cwi.nl\/woudschoten-conferences\/2015-woudschoten-conference\">Fortieth Woudschoten Numerical Analysis Conference<\/a>, organized by the Dutch\/Flemish <i>Werkgemeenschap Scientific Computing<\/i> group. One of the special guests was Th. J. \"Dirk\" Dekker, a retired professor of mathematics and computer science at the University of Amsterdam.<\/p><p>In 1968 Dekker, together with colleague Walter Hoffmann, published a two-part report from the Mathematical Centre in Amsterdam, that described a comprehensive library of Algol 60 procedures for matrix computation. This report preceded by three years the publication of the Wilkinson and Reinsch <i>Handbook for Automatic Computation, Linear Algebra<\/i> that formed the basis for EISPACK and eventually spawned MATLAB.<\/p><h4>Zeroin in Algol<a name=\"466bc11f-42d0-4edc-9cf9-9fe9049d5d3e\"><\/a><\/h4><p>Dekker's program for computing a few of the eigenvalues of a real symmetric tridiagonal matrix involves Sturm sequences and calls a utility procedure, <i>zeroin<\/i>.  Here is a scan of <i>zeroin<\/i> taken from Dekker's 1968 Algol report.  This was my starting point for the MATLAB code that I am about to describe.<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/zeroin_scan.jpg\" alt=\"\"> <\/p><p>Dekker presented <i>zeroin<\/i> at a conference whose proceedings were edited by B. Dejon and Peter Henrici.  Jim Wilkinson described a similar algorithm in a 1967 Stanford report.  Stanford Ph. D. student Richard Brent made important improvements that I will describe in my next blog post. Forsythe, Malcolm and I made Brent's work the basis for the Fortran zero finder in <i>Computer Methods for Mathematical Computations<\/i>. It subsequently evolved into <tt>fzero<\/tt> in MATLAB.<\/p><h4>The test function<a name=\"5cbc4ee9-2180-4542-8aa6-2c9cbda8e02e\"><\/a><\/h4><p>Here is the function that I will use to demonstrate <i>zeroin<\/i>.<\/p><p>$$ f(x) = \\frac{1}{x-3}-6 $$<\/p><pre class=\"codeinput\">   f = @(x) 1.\/(x-3)-6;\r\n   ezplot(f,[3,4])\r\n   hold <span class=\"string\">on<\/span>\r\n   plot(3+1\/6,0,<span class=\"string\">'ro'<\/span>)\r\n   set(gca,<span class=\"string\">'xtick'<\/span>,[3:1\/6:4])\r\n   hold <span class=\"string\">off<\/span>\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/zeroin_blog1_01.png\" alt=\"\"> <p>This function is intentionally tricky.  There is a pole at $x = 3$ and the zero we are trying to find is nearby at $x = 3 \\frac{1}{6}$.  The only portion of the real line where the function is positive is on the left hand one-sixth of the above $x$-axis between these two points.<\/p><h4>Bisection<a name=\"e996c443-bc81-4bd3-82a2-5e7a7016e989\"><\/a><\/h4><p>Dekker's algorithm needs to start with an interval $[a,b]$ on which the given function $f(x)$ changes sign.  The notion of a continuous function of a real variable becomes a bit elusive in floating point arithmetic, so we set our goal to be finding a much smaller subinterval on which the function changes sign. I have simplified the calling sequence by eliminating the tolerance specifying the length of the convergence interval.  All of the following functions iterate until the length of the interval $[a,b]$ is of size roundoff error in the endpoint $b$.<\/p><p>The reliable, fail-safe portion of <i>zeroin<\/i> is the bisection algorithm. The idea is to repeatedly cut the interval $[a,b]$ in half, while continuing to span a sign change.  The actual function values are not used, only the signs.  Here is the code for bisection by itself.<\/p><pre class=\"codeinput\">   type <span class=\"string\">bisect<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction b = bisect(f,a,b)\r\n   s = '%5.0f %19.15f %19.15f\\n';\r\n   fprintf(s,1,a,f(a))\r\n   fprintf(s,2,b,f(b))\r\n   k = 2;\r\n   while abs(b-a) &gt; eps(abs(b))\r\n      x = (a + b)\/2;\r\n      if sign(f(x)) == sign(f(b))\r\n         b = x;\r\n      else \r\n         a = x;\r\n      end\r\n      k = k+1;\r\n      fprintf(s,k,x,f(x))\r\n   end\r\nend\r\n<\/pre><p>Let's see how <tt>bisect<\/tt> performs on our test function. The interval <tt>[3,4]<\/tt> provides a satisfactory starting interval because IEEE floating point arithmetic generates a properly signed infinity, <tt>+Inf<\/tt>, at the pole.<\/p><pre class=\"codeinput\">   bisect(f,3,4)\r\n<\/pre><pre class=\"codeoutput\">    1   3.000000000000000                 Inf\r\n    2   4.000000000000000  -5.000000000000000\r\n    3   3.500000000000000  -4.000000000000000\r\n    4   3.250000000000000  -2.000000000000000\r\n    5   3.125000000000000   2.000000000000000\r\n    6   3.187500000000000  -0.666666666666667\r\n    7   3.156250000000000   0.400000000000000\r\n    8   3.171875000000000  -0.181818181818182\r\n    9   3.164062500000000   0.095238095238095\r\n   10   3.167968750000000  -0.046511627906977\r\n   11   3.166015625000000   0.023529411764706\r\n   12   3.166992187500000  -0.011695906432749\r\n   13   3.166503906250000   0.005865102639296\r\n   14   3.166748046875000  -0.002928257686676\r\n   15   3.166625976562500   0.001465201465201\r\n   16   3.166687011718750  -0.000732332478946\r\n   17   3.166656494140625   0.000366233290606\r\n   18   3.166671752929688  -0.000183099880985\r\n   19   3.166664123535156   0.000091554131380\r\n   20   3.166667938232422  -0.000045776017944\r\n   21   3.166666030883789   0.000022888270905\r\n   22   3.166666984558106  -0.000011444069969\r\n   23   3.166666507720947   0.000005722051355\r\n   24   3.166666746139526  -0.000002861021585\r\n   25   3.166666626930237   0.000001430511816\r\n   26   3.166666686534882  -0.000000715255652\r\n   27   3.166666656732559   0.000000357627890\r\n   28   3.166666671633720  -0.000000178813929\r\n   29   3.166666664183140   0.000000089406969\r\n   30   3.166666667908430  -0.000000044703484\r\n   31   3.166666666045785   0.000000022351742\r\n   32   3.166666666977108  -0.000000011175871\r\n   33   3.166666666511446   0.000000005587935\r\n   34   3.166666666744277  -0.000000002793968\r\n   35   3.166666666627862   0.000000001396984\r\n   36   3.166666666686069  -0.000000000698492\r\n   37   3.166666666656965   0.000000000349246\r\n   38   3.166666666671517  -0.000000000174623\r\n   39   3.166666666664241   0.000000000087311\r\n   40   3.166666666667879  -0.000000000043656\r\n   41   3.166666666666060   0.000000000021828\r\n   42   3.166666666666970  -0.000000000010914\r\n   43   3.166666666666515   0.000000000005457\r\n   44   3.166666666666743  -0.000000000002728\r\n   45   3.166666666666629   0.000000000001364\r\n   46   3.166666666666686  -0.000000000000682\r\n   47   3.166666666666657   0.000000000000341\r\n   48   3.166666666666671  -0.000000000000171\r\n   49   3.166666666666664   0.000000000000085\r\n   50   3.166666666666668  -0.000000000000043\r\n   51   3.166666666666666   0.000000000000021\r\n   52   3.166666666666667  -0.000000000000011\r\n   53   3.166666666666667   0.000000000000005\r\nans =\r\n   3.166666666666667\r\n<\/pre><p>You can see the length of the interval being cut in half in the early stages of the <tt>x<\/tt> values and then the values of <tt>f(x)<\/tt> being cut in half in the later stages.  This takes 53 iterations.  In fact, with my stopping criterion involving <tt>eps(b)<\/tt>, if the starting <tt>a<\/tt> and <tt>b<\/tt> are of comparable size, <tt>bisect<\/tt> takes 53 iterations for <i>any<\/i> function because the double precision floating point fraction has 52 bits.<\/p><p>Except for some cases of pathological behavior near the end points that I will describe in my next post, bisection is guaranteed to converge. If you have a starting interval with a sign change, this algorithm will almost certainly find a tiny subinterval with a sign change.  But it is too slow.  53 iterations is usually too many.  We need something that takes many fewer steps.<\/p><h4>Secant method<a name=\"7e654aa5-e93a-4189-9712-06a4b9bc42c1\"><\/a><\/h4><p>If we have two iterates $a$ and $b$, with corresponding function values, whether or not they exhibit a sign change, we can take the next iterate to be the point where the straight line through $[a,f(a)]$ and $[b,f(b)]$ intersects the $x$ -axis.  This is the secant method.  Here is the code, with the computation of the secant done carefully to avoid unnecessary underflow or overflow.  The trouble is there is nothing to prevent the next point from being outside the interval.<\/p><pre class=\"codeinput\">   type <span class=\"string\">secant<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction b = secant(f,a,b)\r\n   s = '%5.0f %19.15f %23.15e\\n';\r\n   fprintf(s,1,a,f(a))\r\n   fprintf(s,2,b,f(b))\r\n   k = 2;\r\n   while abs(b-a) &gt; eps(abs(b))\r\n      c = a;\r\n      a = b;\r\n      b = b + (b - c)\/(f(c)\/f(b) - 1);\r\n      k = k+1;\r\n      fprintf(s,k,b,f(b))\r\n   end\r\n<\/pre><p>A secant through the point at infinity does not make sense, so we have to start with a slightly smaller interval, but even this one soon gets into trouble.<\/p><pre class=\"codeinput\">   secant(f,3.1,3.5)\r\n<\/pre><pre class=\"codeoutput\">    1   3.100000000000000   3.999999999999991e+00\r\n    2   3.500000000000000  -4.000000000000000e+00\r\n    3   3.300000000000000  -2.666666666666665e+00\r\n    4   2.900000000000000  -1.600000000000004e+01\r\n    5   3.379999999999999  -3.368421052631575e+00\r\n    6   3.507999999999999  -4.031496062992121e+00\r\n    7   2.729760000000002  -9.700414446418032e+00\r\n    8   4.061451519999991  -5.057893854634069e+00\r\n    9   5.512291472588765  -5.601957013781701e+00\r\n   10  -9.426310620985419  -6.080474408736508e+00\r\n   11 180.397066104648590  -5.994362928192905e+00\r\n   12 13394.317035490772000  -5.999925324746076e+00\r\n   13 -14239910.406137096000000  -6.000000070225146e+00\r\n   14 1144132943372.700400000000000  -5.999999999999126e+00\r\n   15 97754126073011126000.000000000000000  -6.000000000000000e+00\r\n   16 -671105854578598490000000000000000.000000000000000  -6.000000000000000e+00\r\n   17                -Inf  -6.000000000000000e+00\r\nans =\r\n  -Inf\r\n<\/pre><p>This shows the secant method can be unreliable.  At step 4 the secant through the previous two points meets the $x$ -axis at $x = 2.9$ which is already outside the initial interval.  By step 10 the iterate has jumped to the other branch of the function.  And at step 12 the computed values exceed my output format.<\/p><p>Interestingly, if we reverse the roles of $a$ and $b$ in this case, <tt>secant<\/tt> does not escape the interval.<\/p><pre class=\"codeinput\">   secant(f,3.5,3.1)\r\n<\/pre><pre class=\"codeoutput\">    1   3.500000000000000  -4.000000000000000e+00\r\n    2   3.100000000000000   3.999999999999991e+00\r\n    3   3.300000000000000  -2.666666666666665e+00\r\n    4   3.220000000000000  -1.454545454545450e+00\r\n    5   3.124000000000000   2.064516129032251e+00\r\n    6   3.180320000000000  -4.543034605146419e-01\r\n    7   3.170161920000000  -1.232444955957233e-01\r\n    8   3.166380335513600   1.032566086067455e-02\r\n    9   3.166672671466170  -2.161649939527166e-04\r\n   10   3.166666676982834  -3.713819847206423e-07\r\n   11   3.166666666666295   1.338662514172029e-11\r\n   12   3.166666666666667   5.329070518200751e-15\r\n   13   3.166666666666667   5.329070518200751e-15\r\nans =\r\n   3.166666666666667\r\n<\/pre><h4>Zeroin algorithm<a name=\"cf04caae-7244-461a-80a6-446ecb6d3e21\"><\/a><\/h4><p>Dekker's algorithm keeps three variables, $a$, $b$, and $c$: Initially, $c$ is set equal to $a$.<\/p><div><ul><li>$b$ is the best zero so far, in the sense that $f(b)$ is the smallest value of $f(x)$ so far.<\/li><li>$a$ is the previous value of $b$, so $a$ and $b$ produce the secant.<\/li><li>$c$ and $b$ bracket the sign change, so $b$ and $c$ provide the midpoint.<\/li><\/ul><\/div><p>At each iteration, the choice is made from three possible steps:<\/p><div><ul><li>a minimal step if the secant step is within roundoff error of an interval endpoint.<\/li><li>a secant step if that step is within the interval and not within roundoff error of an endpoint.<\/li><li>a bisection step otherwise.<\/li><\/ul><\/div><h4>Zeroin in MATLAB<a name=\"6977f5dc-1980-4a0d-bb63-f67c61d81b93\"><\/a><\/h4><pre class=\"codeinput\">   type <span class=\"string\">zeroin<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction b = zeroin(f,a,b)\r\n   % Zeroin(f,a,b) reduces the interval [a,b] to a tiny\r\n   % interal on which the function f(x) changes sign.\r\n   % Returns one end point of that interval.\r\n\r\n   fa = f(a);\r\n   fb = f(b);\r\n   if sign(fa) == sign(fb)\r\n      error('Sorry, f(x) must change sign on the interval [a,b].')\r\n   end\r\n   s = '%5.0f %8s %19.15f %23.15e\\n';\r\n   fprintf(s,1,'initial',a,fa)\r\n   fprintf(s,2,'initial',b,fb)\r\n   k = 2;\r\n\r\n   % a is the previous value of b and [b, c] always contains the zero.\r\n   c = a; fc = fa;\r\n   while true\r\n      if sign(fb) == sign(fc)\r\n         c = a; fc = fa;\r\n      end\r\n\r\n      % Swap to insure f(b) is the smallest value so far.\r\n      if abs(fc) &lt; abs(fb)\r\n         a = b; fa = fb;\r\n         b = c; fb = fc;\r\n         c = a; fc = fa;\r\n      end\r\n\r\n      % Midpoint.\r\n      m = (b + c)\/2;\r\n      if abs(m - b) &lt;= eps(abs(b))\r\n         return   % Exit from the loop and the function here.\r\n      end\r\n\r\n      % p\/q is the the secant step.\r\n      p = (b - a)*fb;\r\n      if p &gt;= 0\r\n         q = fa - fb;\r\n      else\r\n         q = fb - fa;\r\n         p = -p;\r\n      end\r\n      % Save this point.\r\n      a = b; fa = fb;\r\n      k = k+1;\r\n\r\n      % Choose next point.\r\n      if p &lt;= eps(q)\r\n         % Minimal step.\r\n         b = b + sign(c-b)*eps(b);\r\n         fb = f(b);\r\n         fprintf(s,k,'minimal',b,fb)\r\n\r\n      elseif p &lt;= (m - b)*q\r\n         % Secant.\r\n         b = b + p\/q;\r\n         fb = f(b);\r\n         fprintf(s,k,'secant ',b,fb)\r\n\r\n      else\r\n         % Bisection.\r\n         b = m;\r\n         fb = f(b);\r\n         fprintf(s,k,'bisect ',b,fb)\r\n      end\r\n   end\r\nend\r\n<\/pre><p>And here is the performance on our test function.<\/p><pre class=\"codeinput\">   zeroin(f,3,4)\r\n<\/pre><pre class=\"codeoutput\">    1  initial   3.000000000000000                     Inf\r\n    2  initial   4.000000000000000  -5.000000000000000e+00\r\n    3  secant    4.000000000000000  -5.000000000000000e+00\r\n    4  minimal   3.999999999999999  -4.999999999999999e+00\r\n    5  bisect    3.500000000000000  -3.999999999999998e+00\r\n    6  bisect    3.250000000000000  -2.000000000000000e+00\r\n    7  bisect    3.125000000000000   2.000000000000000e+00\r\n    8  secant    3.187500000000000  -6.666666666666670e-01\r\n    9  secant    3.171875000000000  -1.818181818181817e-01\r\n   10  secant    3.166015625000000   2.352941176470580e-02\r\n   11  secant    3.166687011718750  -7.323324789458852e-04\r\n   12  secant    3.166666746139526  -2.861021584976697e-06\r\n   13  secant    3.166666666656965   3.492459654808044e-10\r\n   14  secant    3.166666666666667   5.329070518200751e-15\r\n   15  minimal   3.166666666666667  -1.065814103640150e-14\r\nans =\r\n   3.166666666666667\r\n<\/pre><p>The minimal step even helps get things away from the pole. A few bisection steps get the interval down to<\/p><p>$$ [3 \\frac{1}{8}, 3 \\frac{1}{4}] $$<\/p><p>Then secant can safely take over and obtain a zero in half a dozen steps.<\/p><h4>References<a name=\"1fb0d597-b7a3-4e61-9584-3e263b3e2070\"><\/a><\/h4><p>T. J. Dekker and W. Hoffmann (part 2), <i>Algol 60 Procedures in Numerical Algebra<\/i>, <i>Parts 1 and 2<\/i>, Tracts 22 and 23, Mathematisch Centrum Amsterdam, 1968.<\/p><p>T. J. Dekker, \"Finding a zero by means of successive linear interpolation\", in B. Dejon and P. Henrici (eds.), <i>Constructive aspects of the fundamental theorem of algebra<\/i>, Interscience, New York, 1969.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_57b24e5247914d2882c155812963d252() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='57b24e5247914d2882c155812963d252 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 57b24e5247914d2882c155812963d252';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2015 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_57b24e5247914d2882c155812963d252()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2015a<br><\/p><\/div><!--\r\n57b24e5247914d2882c155812963d252 ##### SOURCE BEGIN #####\r\n%% ZEROIN, Part 1: Dekker's Algorithm\r\n% Th. J. Dekker's _zeroin_ algorithm from 1969 is one of my favorite algorithms.\r\n% An elegant technique combining bisection and the secant method for finding\r\n% a zero of a function of a real variable, it has become |fzero| in MATLAB\r\n% today.  This is the first of a three part series.\r\n\r\n%% Dirk Dekker\r\n% I have just come from the \r\n% <https:\/\/wsc.project.cwi.nl\/woudschoten\/2015\/conferentieE.php\r\n% Fortieth Woudschoten Numerical Analysis Conference>, organized by the\r\n% Dutch\/Flemish _Werkgemeenschap Scientific Computing_ group.\r\n% One of the special guests was Th. J. \"Dirk\" Dekker, a retired professor\r\n% of mathematics and computer science at the University of Amsterdam.\r\n\r\n%%\r\n% In 1968 Dekker, together with colleague Walter Hoffmann, published a two-part\r\n% report from the Mathematical Centre in Amsterdam, that described a\r\n% comprehensive library of Algol 60 procedures for matrix computation.  \r\n% This report preceded by three years the publication of the Wilkinson and\r\n% Reinsch _Handbook for Automatic Computation, Linear Algebra_ that\r\n% formed the basis for EISPACK and eventually spawned MATLAB.\r\n\r\n%% Zeroin in Algol\r\n% Dekker's program for computing a few of the eigenvalues of a real\r\n% symmetric tridiagonal matrix involves Sturm sequences and calls a utility\r\n% procedure, _zeroin_.  Here is a scan of _zeroin_ taken from Dekker's\r\n% 1968 Algol report.  This was my starting point for the MATLAB code\r\n% that I am about to describe.\r\n% \r\n% <<zeroin_scan.jpg>>\r\n%\r\n\r\n%%\r\n% Dekker presented _zeroin_ at a conference whose proceedings were edited\r\n% by B. Dejon and Peter Henrici.  Jim Wilkinson described a similar algorithm\r\n% in a 1967 Stanford report.  Stanford Ph. D. student Richard Brent made\r\n% important improvements that I will describe in my next blog post.\r\n% Forsythe, Malcolm and I made Brent's work the basis for the Fortran\r\n% zero finder in _Computer Methods for Mathematical Computations_.\r\n% It subsequently evolved into |fzero| in MATLAB.\r\n\r\n%% The test function\r\n% Here is the function that I will use to demonstrate _zeroin_.\r\n%\r\n% $$ f(x) = \\frac{1}{x-3}-6 $$\r\n%\r\n\r\n   f = @(x) 1.\/(x-3)-6;\r\n   ezplot(f,[3,4])\r\n   hold on\r\n   plot(3+1\/6,0,'ro')\r\n   set(gca,'xtick',[3:1\/6:4])\r\n   hold off\r\n\r\n%%\r\n% This function is intentionally tricky.  There is a pole at $x = 3$ and\r\n% the zero we are trying to find is nearby at $x = 3 \\frac{1}{6}$.  The\r\n% only portion of the real line where the function is positive is on the\r\n% left hand one-sixth of the above $x$-axis between these two points.\r\n\r\n%% Bisection\r\n% Dekker's algorithm needs to start with an interval $[a,b]$ on which\r\n% the given function $f(x)$ changes sign.  The notion of a continuous\r\n% function of a real variable becomes a bit elusive in floating point\r\n% arithmetic, so we set our goal to be finding a much smaller subinterval\r\n% on which the function changes sign.\r\n% I have simplified the calling sequence by eliminating the tolerance\r\n% specifying the length of the convergence interval.  All of the following\r\n% functions iterate until the length of the interval $[a,b]$ is of size\r\n% roundoff error in the endpoint $b$.\r\n\r\n%%\r\n% The reliable, fail-safe portion of _zeroin_ is the bisection algorithm.\r\n% The idea is to repeatedly cut the interval $[a,b]$ in half,\r\n% while continuing to span a sign change.  The actual function values are\r\n% not used, only the signs.  Here is the code for bisection by itself.\r\n\r\n   type bisect\r\n\r\n%%\r\n% Let's see how |bisect| performs on our test function.\r\n% The interval |[3,4]| provides a satisfactory starting interval\r\n% because IEEE floating point arithmetic generates a properly signed infinity,\r\n% |+Inf|, at the pole.\r\n\r\n   bisect(f,3,4)\r\n\r\n%%\r\n% You can see the length of the interval being cut in half in the\r\n% early stages of the |x| values and then the values of |f(x)| being cut\r\n% in half in the later stages.  This takes 53 iterations.  In fact, with\r\n% my stopping criterion involving |eps(b)|, if the starting |a| and |b| are\r\n% of comparable size, |bisect| takes 53 iterations for _any_ function\r\n% because the double precision floating point fraction has 52 bits.\r\n\r\n%%\r\n% Except for some cases of pathological behavior near the end points that\r\n% I will describe in my next post, bisection is guaranteed to converge.\r\n% If you have a starting interval with a sign change, this algorithm will\r\n% almost certainly find a tiny subinterval with a sign change.  But it is\r\n% too slow.  53 iterations is usually too many.  We need something that takes\r\n% many fewer steps.\r\n\r\n%% Secant method\r\n% If we have two iterates $a$ and $b$, with corresponding function values,\r\n% whether or not they exhibit a sign change, we can take the next iterate\r\n% to be the point where the straight line through $[a,f(a)]$ and $[b,f(b)]$\r\n% intersects the $x$ -axis.  This is the secant method.  Here is the code,\r\n% with the computation of the secant done carefully to avoid unnecessary\r\n% underflow or overflow.  The trouble is there is nothing to prevent the next\r\n% point from being outside the interval.\r\n\r\n   type secant\r\n\r\n%%\r\n% A secant through the point at infinity does not make sense, so we have to\r\n% start with a slightly smaller interval, but even this one soon gets into\r\n% trouble.\r\n\r\n   secant(f,3.1,3.5)\r\n\r\n%%\r\n% This shows the secant method can be unreliable.  At step 4 the secant\r\n% through the previous two points meets the $x$ -axis at $x = 2.9$ which\r\n% is already outside the initial interval.  By step 10 the iterate has\r\n% jumped to the other branch of the function.  And at step 12 the computed\r\n% values exceed my output format.\r\n\r\n%%\r\n% Interestingly, if we reverse the roles of $a$ and $b$ in this case,\r\n% |secant| does not escape the interval.\r\n\r\n   secant(f,3.5,3.1)\r\n\r\n%% Zeroin algorithm\r\n% Dekker's algorithm keeps three variables, $a$, $b$, and $c$:\r\n% Initially, $c$ is set equal to $a$.\r\n%\r\n% * $b$ is the best zero so far, in the sense that $f(b)$ is the smallest\r\n% value of $f(x)$ so far.\r\n% * $a$ is the previous value of $b$, so $a$ and $b$ produce the secant. \r\n% * $c$ and $b$ bracket the sign change, so $b$ and $c$ provide the midpoint.\r\n\r\n%%\r\n% At each iteration, the choice is made from three possible steps:\r\n%\r\n% * a minimal step if the secant step is within roundoff error of an interval\r\n% endpoint.\r\n% * a secant step if that step is within the interval and not within\r\n% roundoff error of an endpoint.\r\n% * a bisection step otherwise.\r\n\r\n%% Zeroin in MATLAB\r\n\r\n   type zeroin\r\n  \r\n%%\r\n% And here is the performance on our test function.\r\n\r\n   zeroin(f,3,4)\r\n\r\n%%\r\n% The minimal step even helps get things away from the pole.\r\n% A few bisection steps get the interval down to\r\n%\r\n% $$ [3 \\frac{1}{8}, 3 \\frac{1}{4}] $$\r\n%\r\n% Then secant can safely take over and obtain a zero in half a dozen steps.\r\n\r\n\r\n\r\n\r\n%% References\r\n% T. J. Dekker and W. Hoffmann (part 2), _Algol 60 Procedures in Numerical\r\n% Algebra_, _Parts 1 and 2_, Tracts 22 and 23, Mathematisch Centrum\r\n% Amsterdam, 1968.\r\n%\r\n% T. J. Dekker, \"Finding a zero by means of successive linear interpolation\",\r\n% in B. Dejon and P. Henrici (eds.), _Constructive aspects of the fundamental\r\n% theorem of algebra_, Interscience, New York, 1969.\r\n\r\n##### SOURCE END ##### 57b24e5247914d2882c155812963d252\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/feature_image\/zeroin_blog1_01.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div><!--introduction--><p>Th. J. Dekker's <i>zeroin<\/i> algorithm from 1969 is one of my favorite algorithms. An elegant technique combining bisection and the secant method for finding a zero of a function of a real variable, it has become <tt>fzero<\/tt> in MATLAB today.  This is the first of a three part series.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/10\/12\/zeroin-part-1-dekkers-algorithm\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":1272,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11,13,4,16,8],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1269"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=1269"}],"version-history":[{"count":10,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1269\/revisions"}],"predecessor-version":[{"id":1835,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1269\/revisions\/1835"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media\/1272"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=1269"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=1269"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=1269"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}