{"id":1280,"date":"2015-10-26T12:00:53","date_gmt":"2015-10-26T17:00:53","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=1280"},"modified":"2016-04-13T14:45:08","modified_gmt":"2016-04-13T19:45:08","slug":"zeroin-part-2-brents-version","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2015\/10\/26\/zeroin-part-2-brents-version\/","title":{"rendered":"Zeroin, Part 2: Brent&#8217;s Version"},"content":{"rendered":"<div class=\"content\">\r\n\r\n<!--introduction-->Richard Brent's improvements to Dekker's <i>zeroin<\/i> algorithm, published in 1971, made it faster, safer in floating point arithmetic, and guaranteed not to fail.\r\n\r\n<!--\/introduction-->\r\n<h3>Contents<\/h3>\r\n<div>\r\n<ul>\r\n\t<li><a href=\"#cc67c076-9561-49b0-b334-0dc29d2e4c57\">Richard Brent<\/a><\/li>\r\n\t<li><a href=\"#ac81b842-7fbd-4190-8ab3-ee104044bd1c\">Weakness of Zeroin<\/a><\/li>\r\n\t<li><a href=\"#fe8667d3-af9a-43de-80e3-fb75fe1386b4\">Two improvements<\/a><\/li>\r\n\t<li><a href=\"#4f35e8ec-f7d4-4b50-a53b-7e17f222fbb4\">Muller's method<\/a><\/li>\r\n\t<li><a href=\"#8368e5f6-42c7-4b46-ad76-093f494a02fe\">Inverse quadratic interpolation<\/a><\/li>\r\n\t<li><a href=\"#f9b0f31b-1d22-478c-a8a6-c1684ec12f1a\">Brent's algorithm<\/a><\/li>\r\n\t<li><a href=\"#d354b12f-235d-474e-be05-d4dba9fa6fcc\">Fzero<\/a><\/li>\r\n\t<li><a href=\"#f4631cb3-c086-46df-abee-85d128db8d75\">References<\/a><\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>Richard Brent<a name=\"cc67c076-9561-49b0-b334-0dc29d2e4c57\"><\/a><\/h4>\r\nRichard Brent was a graduate student in computer science at Stanford in 1968-71. He wrote a Ph. D. thesis under George Forsythe's direction titled <i>Algorithms for finding zeros and extrema of functions without calculating derivatives<\/i>. I happened to be a visiting professor at Stanford and was on Brent's committee. Forsythe was editor of an important series of books in computer science published by Prentice-Hall. He was so impressed with Brent's thesis that he published it almost unchanged in the series in 1973 with the title <i>Algorithms for Minimization without Derivatives<\/i>.\r\n\r\nThree important algorithms were described in Brent's book. The algorithm for finding a zero of a function of one real variable, which is a successor to Dekker's <i>zeroin<\/i>, is still in use in MATLAB today. This is the algorithm that I want to write about here. The other two algorithms, one for minimizing a function of one real variable and one for minimizing a function of several real variables are no longer in wide spread use.\r\n<h4>Weakness of Zeroin<a name=\"ac81b842-7fbd-4190-8ab3-ee104044bd1c\"><\/a><\/h4>\r\nThe version of <i>zeroin<\/i> that was presented by Dekker and that I described in <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/10\/12\/zeroin-part-1-dekkers-algorithm\/\">Part 1 of this series<\/a> is elegant and concise, but not foolproof. It has trouble with multiple roots. A simple example is provided by\r\n\r\n$$ f(x) = x^3 $$\r\n\r\nIf we start with a symmetric interval spanning the origin, the initial secant step hits the zero exactly, but the stopping criterion fails to recognize the good luck. The code takes a minimal step, which in this case is <tt>eps*realmin<\/tt>, the smallest possible denormal normal floating point number. It will continue essentially forever, taking impossibly small steps.\r\n<pre class=\"language-matlab\">f =\r\n    @(x)x.^3;\r\n<\/pre>\r\n<pre class=\"language-matlab\">zeroin(f,-1,1)\r\n<\/pre>\r\n<pre>     1  initial  -1.000000000000000e+00  -1.000000000000000e+00\r\n     2  initial   1.000000000000000e+00   1.000000000000000e+00\r\n     3  secant    0.000000000000000e+00   0.000000000000000e+00\r\n     4  minimal -4.940656458412465e-324  -0.000000000000000e+00\r\n     5  minimal -9.881312916824931e-324  -0.000000000000000e+00\r\n     6  minimal -1.482196937523740e-323  -0.000000000000000e+00\r\n     7  minimal -1.976262583364986e-323  -0.000000000000000e+00\r\n     8  minimal -2.470328229206233e-323  -0.000000000000000e+00\r\n         .\r\n         .\r\n         .\r\n 12341  minimal -6.095781938389300e-320  -0.000000000000000e+00\r\n 12342  minimal -6.096276004035141e-320  -0.000000000000000e+00\r\n 12343  minimal -6.096770069680982e-320  -0.000000000000000e+00\r\n 12344  minimal -6.097264135326824e-320  -0.000000000000000e+00\r\n         .\r\n         .\r\n         .<\/pre>\r\nA starting interval that is not symmetric about the origin does not work any better. We don't hit the zero exactly, so we take some small secant steps for a while. But eventually we get bogged down repeatedly taking the minimal step, and just making roundoff in roundoff.\r\n<pre class=\"language-matlab\">zeroin(f,-.5,1)\r\n<\/pre>\r\n<pre>     1  initial  -5.000000000000000e-01  -1.250000000000000e-01\r\n     2  initial   1.000000000000000e+00   1.000000000000000e+00\r\n     3  secant   -3.333333333333334e-01  -3.703703703703705e-02\r\n     4  secant   -2.631578947368422e-01  -1.822423093745445e-02\r\n     5  secant   -1.951779563719863e-01  -7.435193904825083e-03\r\n     6  secant   -1.483300289942098e-01  -3.263527261310824e-03\r\n     7  secant   -1.116805310361257e-01  -1.392940003647091e-03\r\n     8  secant   -8.438934127192455e-02  -6.009838348927866e-04\r\n         .\r\n         .\r\n         .\r\n   124  secant   -5.752411316818776e-16  -1.903486478002247e-46\r\n   125  minimal  -5.752411316818775e-16  -1.903486478002246e-46\r\n   126  secant   -4.143774132574865e-16  -7.115218233323205e-47\r\n   127  minimal  -4.143774132574865e-16  -7.115218233323202e-47\r\n   128  minimal  -4.143774132574864e-16  -7.115218233323200e-47\r\n   129  minimal  -4.143774132574864e-16  -7.115218233323197e-47\r\n   130  minimal  -4.143774132574863e-16  -7.115218233323194e-47\r\n   131  minimal  -4.143774132574863e-16  -7.115218233323192e-47\r\n   132  minimal  -4.143774132574862e-16  -7.115218233323189e-47\r\n         .\r\n         .\r\n         .<\/pre>\r\n<h4>Two improvements<a name=\"fe8667d3-af9a-43de-80e3-fb75fe1386b4\"><\/a><\/h4>\r\nTo avoid these difficulties, Brent made two important improvements to Dekker's algorithm. One is a check to see if an exact zero has accidentally been found. Another is to introduce the variable <tt>e<\/tt> in the code below that keeps track of the step size ratios in two successive steps. This variable is involved in the choice about whether to use bisection. As a result Brent's algorithm never takes more than twice as many steps as bisection alone.\r\n\r\nLet's see how these improvements affect our example in <tt>fzero<\/tt>.\r\n<pre class=\"language-matlab\">opt = optimset(<span class=\"string\">'display'<\/span>,<span class=\"string\">'iter'<\/span>)\r\n<\/pre>\r\nStarting with an interval that is centered about the origin, the exact zero results in an early termination.\r\n<pre class=\"language-matlab\">fzero(f,[-1,1],opt)\r\n<\/pre>\r\n<pre>  Func-count    x          f(x)             Procedure\r\n     2               1             1        initial\r\n     3               0             0        bisection<\/pre>\r\n<pre class=\"language-matlab\">Zero <span class=\"string\">found<\/span> <span class=\"string\">in<\/span> <span class=\"string\">the<\/span> <span class=\"string\">interval<\/span> <span class=\"string\">[-1, 1]<\/span>\r\nans =\r\n    0\r\n<\/pre>\r\nStarting with an interval that is not centered about the origin, the iteration terminates in a reasonable number of steps.\r\n<pre class=\"language-matlab\">fzero(f,[-.5,1],opt)\r\n<\/pre>\r\n<pre>  Func-count    x          f(x)             Procedure\r\n     2            -0.5        -0.125        initial\r\n     3       -0.333333     -0.037037        interpolation\r\n     4       -0.265664    -0.0187499        interpolation\r\n     5       -0.197929     -0.007754        interpolation\r\n     6       -0.197929     -0.007754        bisection\r\n     7       -0.133649   -0.00238724        interpolation\r\n         .\r\n         .\r\n         .\r\n   151    -4.51944e-16   -9.2311e-47        interpolation\r\n   152    -7.85458e-18  -4.84584e-52        interpolation\r\n   153    -7.85458e-18  -4.84584e-52        bisection\r\n   154    -7.85458e-18  -4.84584e-52        interpolation<\/pre>\r\n<pre class=\"language-matlab\">Zero <span class=\"string\">found<\/span> <span class=\"string\">in<\/span> <span class=\"string\">the<\/span> <span class=\"string\">interval<\/span> <span class=\"string\">[-0.5, 1]<\/span>\r\nans =\r\n    -7.854580142952130e-18\r\n<\/pre>\r\n<h4>Muller's method<a name=\"4f35e8ec-f7d4-4b50-a53b-7e17f222fbb4\"><\/a><\/h4>\r\nDuring the course of the iterations we often have three distinct values of $x$, namely $a$, $b$, and $c$, and corresponding function values, $f(a)$, $f(b)$, and $f(c)$. We could interpolate these three values with a quadratic in $x$, find the roots of this quadratic, and take the root that is closest to our best approximation so far to be the next approximation to the zero of $f(x)$. This is Muller's method. The trouble is the roots of the quadratic might be complex. This would be an advantage if we were looking for complex solutions, and Muller's method is often used in this situation. But we don't want them here and Brent chose to do something else.\r\n<h4>Inverse quadratic interpolation<a name=\"8368e5f6-42c7-4b46-ad76-093f494a02fe\"><\/a><\/h4>\r\nInstead of a quadratic in $x$, we can interpolate the three points with a quadratic function in $y$. That's a ``sideways'' parabola, $P(y)$, determined by the interpolation conditions\r\n\r\n$$ a = P(f(a)), \\ b = P(f(b)), \\ c = P(f(c)) $$\r\n\r\nThis parabola always intersects the $x$-axis, which is $y = 0$. So $x = P(0)$ is the next iterate.\r\n\r\nThis method is known as <i>Inverse Quadratic Interpolation<\/i>, abbreviated IQI. Here is a naive implementation that illustrates the idea. It uses <tt>polyinterp<\/tt>, taken from <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/37976-numerical-computing-with-matlab\"><i>Numerical Computing with MATLAB<\/i><\/a>\r\n<pre class=\"codeinput\">   type <span class=\"string\">iqi<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">  function x = iqi(f,a,b,c)\r\n     while abs(c-b) &gt; eps(abs(b))\r\n        x = polyinterp([f(a),f(b),f(c)],[a,b,c],0);\r\n        disp(x)\r\n        a = b;\r\n        b = c;\r\n        c = x;\r\n     end\r\n  end\r\n<\/pre>\r\nOne difficulty with this simple code is that polynomial interpolation requires the abscissae, which in this case are $f(a)$, $f(b)$, and $f(c)$, to be distinct. We have no guarantee that they are. For example, if we try to compute $\\sqrt{2}$ using $f(x) = x^2 - 2$ and start with $a = -2, b = 0, c = 2$, we are starting with $f(a) = f(c)$ and the first step is undefined. If we start near this singular situation, say with $a = -2.001, b = 0, c = 1.999$, the next iterate is near $x = 500$.\r\n\r\nSo IQI converges faster than secant once it is near the zero, but its global behavior is unpredictable. It needs to be combined with bisection to have a reliable algorithm.\r\n<h4>Brent's algorithm<a name=\"f9b0f31b-1d22-478c-a8a6-c1684ec12f1a\"><\/a><\/h4>\r\nLike Dekker's original version, Brent's version of <i>zeroin<\/i> starts with an interval $[a,b]$ on which the function $f(x)$ changes sign. The objective is to reduce the interval to a tiny subinterval on which the function still changes sign. The variables $a$, $b$, and $c$ play the same role:\r\n<div>\r\n<ul>\r\n\t<li>$b$ is the best zero so far, in the sense that $f(b)$ is the smallest value of $f(x)$ so far.<\/li>\r\n\t<li>$a$ is the previous value of $b$, so $a$ and $b$ produce the secant.<\/li>\r\n\t<li>$c$ and $b$ bracket the sign change, so $b$ and $c$ provide the midpoint.<\/li>\r\n<\/ul>\r\n<\/div>\r\nAt each iteration, the choice is made from four possible steps:\r\n<div>\r\n<ul>\r\n\t<li>If $a \\ne c$ , try inverse quadratic interpolation.<\/li>\r\n\t<li>If $a = c$, try secant interpolation.<\/li>\r\n\t<li>If the interpolation step is near the endpoint, or outside the interval, use bisection.<\/li>\r\n\t<li>If the step is smaller than the tolerance, use the tolerance.<\/li>\r\n<\/ul>\r\n<\/div>\r\nHere is a somewhat cluttered plot that shows most of these possibilities. The three values $a$, $b$, and $c$ are plotted in black. The function\r\n\r\n$$ f(x) = x^3 - 3x - 2 $$\r\n\r\nis plotted in light grey to emphasize that the algorithm so far has only evaluated it three times to give the three values shown by the circles.\r\n<pre class=\"codeinput\">   zeroin_plot(1,2.4)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/zeroin_blog2_01.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nThe green line is the secant determined by $a$ and $b$. A secant step would go to the green $\\times$ labeled $t$, but that's outside the interval so <i>zeroin<\/i> would not go there. Actually Brent's <i>zeroin<\/i> would not even be considering the secant in this situation because $a$ and $c$ are distinct.\r\n\r\nThe blue curve is the quadratic in $y$ that goes through the three circles. It crosses the $x$ -axis at the blue $\\times$ labeled $q$, so that would be the IQI step. However, in this case $q$ is more than three-quarters of the way from $b$ to $c$, so Brent would choose bisection instead.\r\n\r\nThe red $\\times$ labeled $m$ is the bisection point halfway between $b$ and $c$. With a value of 2.025 it happens to be really close to the zero at 2. Lucky shot.\r\n<h4>Fzero<a name=\"d354b12f-235d-474e-be05-d4dba9fa6fcc\"><\/a><\/h4>\r\nBrent included Algol and Fortran implementations of this algorithm in his book. Forsythe, Mike Malcolm, and I made the Fortran program the basis for the zero finding chapter of <i>Computer Methods for Mathematical Computations<\/i>. The code is <a href=\"http:\/\/www.netlib.org\/fmm\/zeroin.f\">still available from netlib<\/a>.\r\n\r\nWe named the MATLAB version <tt>fzero<\/tt>. Here is the textbook version, <tt>fzerotx<\/tt>, from <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/37976-numerical-computing-with-matlab\"><i>Numerical Computing with MATLAB<\/i><\/a>. This is not all of the MATLAB <tt>fzero<\/tt>. In my next blog post I will deal with the need to know an interval on which the function changes sign.\r\n<pre class=\"codeinput\">   type <span class=\"string\">fzerotx<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">function b = fzerotx(F,ab,varargin)\r\n%FZEROTX  Textbook version of FZERO.\r\n%   x = fzerotx(F,[a,b]) tries to find a zero of F(x) between a and b.\r\n%   F(a) and F(b) must have opposite signs.  fzerotx returns one \r\n%   end point of a small subinterval of [a,b] where F changes sign.\r\n%   Arguments beyond the first two, fzerotx(F,[a,b],p1,p2,...),\r\n%   are passed on, F(x,p1,p2,..).\r\n%\r\n%   Examples:\r\n%      fzerotx(@sin,[1,4])\r\n%      F = @(x) sin(x); fzerotx(F,[1,4])\r\n\r\n%   Copyright 2014 Cleve Moler\r\n%   Copyright 2014 The MathWorks, Inc.\r\n\r\n% Initialize.\r\na = ab(1);\r\nb = ab(2);\r\nfa = F(a,varargin{:});\r\nfb = F(b,varargin{:});\r\nif sign(fa) == sign(fb)\r\n   error('Function must change sign on the interval')\r\nend\r\nc = a;\r\nfc = fa;\r\nd = b - c;\r\ne = d;\r\n\r\n% Main loop, exit from middle of the loop\r\nwhile fb ~= 0\r\n\r\n   % The three current points, a, b, and c, satisfy:\r\n   %    f(x) changes sign between a and b.\r\n   %    abs(f(b)) &lt;= abs(f(a)).\r\n   %    c = previous b, so c might = a.\r\n   % The next point is chosen from\r\n   %    Bisection point, (a+b)\/2.\r\n   %    Secant point determined by b and c.\r\n   %    Inverse quadratic interpolation point determined\r\n   %    by a, b, and c if they are distinct.\r\n\r\n   if sign(fa) == sign(fb)\r\n      a = c;  fa = fc;\r\n      d = b - c;  e = d;\r\n   end\r\n   if abs(fa) &lt; abs(fb)\r\n      c = b;    b = a;    a = c;\r\n      fc = fb;  fb = fa;  fa = fc;\r\n   end\r\n   \r\n   % Convergence test and possible exit\r\n   m = 0.5*(a - b);\r\n   tol = 2.0*eps*max(abs(b),1.0);\r\n   if (abs(m) &lt;= tol) | (fb == 0.0)\r\n      break\r\n   end\r\n   \r\n   % Choose bisection or interpolation\r\n   if (abs(e) &lt; tol) | (abs(fc) &lt;= abs(fb))\r\n      % Bisection\r\n      d = m;\r\n      e = m;\r\n   else\r\n      % Interpolation\r\n      s = fb\/fc;\r\n      if (a == c)\r\n         % Linear interpolation (secant)\r\n         p = 2.0*m*s;\r\n         q = 1.0 - s;\r\n      else\r\n         % Inverse quadratic interpolation\r\n         q = fc\/fa;\r\n         r = fb\/fa;\r\n         p = s*(2.0*m*q*(q - r) - (b - c)*(r - 1.0));\r\n         q = (q - 1.0)*(r - 1.0)*(s - 1.0);\r\n      end;\r\n      if p &gt; 0, q = -q; else p = -p; end;\r\n      % Is interpolated point acceptable\r\n      if (2.0*p &lt; 3.0*m*q - abs(tol*q)) &amp; (p &lt; abs(0.5*e*q))\r\n         e = d;\r\n         d = p\/q;\r\n      else\r\n         d = m;\r\n         e = m;\r\n      end;\r\n   end\r\n   \r\n   % Next point\r\n   c = b;\r\n   fc = fb;\r\n   if abs(d) &gt; tol\r\n      b = b + d;\r\n   else\r\n      b = b - sign(b-a)*tol;\r\n   end\r\n   fb = F(b,varargin{:});\r\nend\r\n<\/pre>\r\n<h4>References<a name=\"f4631cb3-c086-46df-abee-85d128db8d75\"><\/a><\/h4>\r\nRichard P. Brent, \"An algorithm with guaranteed convergence for finding a zero of a function,\" <i>Computer Journal<\/i> 14, 422-25, 1971.\r\n\r\nRichard P. Brent, <i>Algorithms for Minimization without Derivatives<\/i>, Prentice-Hall, Englewood Cliffs, New Jersey, 1973, 195 pp. Reprinted by Dover, New York, 2002 and 2013. <a href=\"http:\/\/maths-people.anu.edu.au\/~brent\/pub\/pub011.html\">&lt;http:\/\/maths-people.anu.edu.au\/~brent\/pub\/pub011.html<\/a>&gt;\r\n\r\nGeorge Forsythe, Michael Malcolm, and Cleve Moler, <i>Computer Methods for Mathematical Computations<\/i>, Prentice-Hall, 259pp, 1977.\r\n\r\n<script>\/\/ <![CDATA[\r\nfunction grabCode_a8127df70bd8474491ad3dec2945b6ac() {\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='a8127df70bd8474491ad3dec2945b6ac ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' a8127df70bd8474491ad3dec2945b6ac';\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('\r\n\r\n\r\n\r\n\r\n\r\n<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>\r\n\r\n\r\n\r\n\r\n\r\n\r\n\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }\r\n\/\/ ]]><\/script>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\"><a><span style=\"font-size: x-small; font-style: italic;\">Get\r\nthe MATLAB code<noscript>(requires JavaScript)<\/noscript><\/span><\/a><\/p>\r\nPublished with MATLAB\u00ae R2015a\r\n\r\n<\/div>\r\n<!--\r\na8127df70bd8474491ad3dec2945b6ac ##### SOURCE BEGIN #####\r\n%% ZEROIN, Part 2: Brent's Version\r\n% Richard Brent's improvements to Dekker's _zeroin_ algorithm, published in\r\n% 1971, made it faster, safer in floating point arithmetic, and guaranteed\r\n% not to fail.\r\n\r\n%% Richard Brent\r\n% Richard Brent was a graduate student in computer science at Stanford in\r\n% 1968-71.  He wrote a Ph. D. thesis under George Forsythe's direction\r\n% titled _Algorithms for finding zeros and extrema of functions without\r\n% calculating derivatives_.  I happened to be a visiting professor at\r\n% Stanford and was on Brent's committee.  Forsythe was editor of an\r\n% important series of books in computer science published by Prentice-Hall.\r\n% He was so impressed with Brent's thesis that he published it almost\r\n% unchanged in the series in 1973 with the title _Algorithms for Minimization\r\n% without Derivatives_.\r\n\r\n%%\r\n% Three important algorithms were described in Brent's book.\r\n% The algorithm for finding a zero of a function of one real variable,\r\n% which is a successor to Dekker's _zeroin_, is still in use in MATLAB today.\r\n% This is the algorithm that I want to write about here.\r\n% The other two algorithms, one for minimizing a function of one real\r\n% variable and one for minimizing a function of several real variables\r\n% are no longer in wide spread use.\r\n\r\n%% Weakness of Zeroin\r\n% The version of _zeroin_ that was presented by Dekker and that I described in\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2015\/10\/12\/zeroin-part-1-dekkers-algorithm\/ % Part 1 of this series> is elegant and concise, but not foolproof.\r\n% It has trouble with multiple roots.  A simple example is provided by\r\n%\r\n% $$ f(x) = x^3 $$\r\n%\r\n% If we start with a symmetric interval spanning the origin,\r\n% the initial secant step hits the zero exactly, but the stopping\r\n% criterion fails to recognize the good luck.  The code takes a minimal\r\n% step, which in this case is |eps*realmin|, the smallest possible\r\n% denormal normal floating point number.  It will continue essentially\r\n% forever, taking impossibly small steps.\r\n%\r\n%   f =\r\n%       @(x)x.^3;\r\n%\r\n%   zeroin(f,-1,1)\r\n%\r\n%       1  initial  -1.000000000000000e+00  -1.000000000000000e+00\r\n%       2  initial   1.000000000000000e+00   1.000000000000000e+00\r\n%       3  secant    0.000000000000000e+00   0.000000000000000e+00\r\n%       4  minimal -4.940656458412465e-324  -0.000000000000000e+00\r\n%       5  minimal -9.881312916824931e-324  -0.000000000000000e+00\r\n%       6  minimal -1.482196937523740e-323  -0.000000000000000e+00\r\n%       7  minimal -1.976262583364986e-323  -0.000000000000000e+00\r\n%       8  minimal -2.470328229206233e-323  -0.000000000000000e+00\r\n%           .\r\n%           .\r\n%           .\r\n%   12341  minimal -6.095781938389300e-320  -0.000000000000000e+00\r\n%   12342  minimal -6.096276004035141e-320  -0.000000000000000e+00\r\n%   12343  minimal -6.096770069680982e-320  -0.000000000000000e+00\r\n%   12344  minimal -6.097264135326824e-320  -0.000000000000000e+00\r\n%           .\r\n%           .\r\n%           .\r\n\r\n%%\r\n% A starting interval that is not symmetric about the origin does not\r\n% work any better.  We don't hit the zero exactly, so we take some\r\n% small secant steps for a while.\r\n% But eventually we get bogged down repeatedly taking the minimal step,\r\n% and just making roundoff in roundoff.\r\n%\r\n%   zeroin(f,-.5,1)\r\n%\r\n%       1  initial  -5.000000000000000e-01  -1.250000000000000e-01\r\n%       2  initial   1.000000000000000e+00   1.000000000000000e+00\r\n%       3  secant   -3.333333333333334e-01  -3.703703703703705e-02\r\n%       4  secant   -2.631578947368422e-01  -1.822423093745445e-02\r\n%       5  secant   -1.951779563719863e-01  -7.435193904825083e-03\r\n%       6  secant   -1.483300289942098e-01  -3.263527261310824e-03\r\n%       7  secant   -1.116805310361257e-01  -1.392940003647091e-03\r\n%       8  secant   -8.438934127192455e-02  -6.009838348927866e-04\r\n%           .\r\n%           .\r\n%           .\r\n%     124  secant   -5.752411316818776e-16  -1.903486478002247e-46\r\n%     125  minimal  -5.752411316818775e-16  -1.903486478002246e-46\r\n%     126  secant   -4.143774132574865e-16  -7.115218233323205e-47\r\n%     127  minimal  -4.143774132574865e-16  -7.115218233323202e-47\r\n%     128  minimal  -4.143774132574864e-16  -7.115218233323200e-47\r\n%     129  minimal  -4.143774132574864e-16  -7.115218233323197e-47\r\n%     130  minimal  -4.143774132574863e-16  -7.115218233323194e-47\r\n%     131  minimal  -4.143774132574863e-16  -7.115218233323192e-47\r\n%     132  minimal  -4.143774132574862e-16  -7.115218233323189e-47\r\n%           .\r\n%           .\r\n%           .\r\n\r\n%% Two improvements\r\n% To avoid these difficulties, Brent made two important improvements to\r\n% Dekker's algorithm.  One is a check to see if an exact zero has accidentally\r\n% been found.  Another is to introduce the variable |e| in the code below\r\n% that keeps track of the step size ratios in two successive steps.  This\r\n% variable is involved in the choice about whether to use bisection.  As a\r\n% result Brent's algorithm never takes more than twice as many steps as\r\n% bisection alone.\r\n\r\n%%\r\n% Let's see how these improvements affect our example in |fzero|.\r\n%\r\n%   opt = optimset('display','iter')\r\n%\r\n% Starting with an interval that is centered about the origin, the exact zero\r\n% results in an early termination.\r\n%\r\n%   fzero(f,[-1,1],opt)\r\n%\r\n%    Func-count    x          f(x)             Procedure\r\n%       2               1             1        initial\r\n%       3               0             0        bisection\r\n%\r\n%   Zero found in the interval [-1, 1]\r\n%   ans =\r\n%       0\r\n\r\n%%\r\n% Starting with an interval that is not centered about the origin,\r\n% the iteration terminates in a reasonable number of steps.\r\n%\r\n%   fzero(f,[-.5,1],opt)\r\n%\r\n%    Func-count    x          f(x)             Procedure\r\n%       2            -0.5        -0.125        initial\r\n%       3       -0.333333     -0.037037        interpolation\r\n%       4       -0.265664    -0.0187499        interpolation\r\n%       5       -0.197929     -0.007754        interpolation\r\n%       6       -0.197929     -0.007754        bisection\r\n%       7       -0.133649   -0.00238724        interpolation\r\n%           .\r\n%           .\r\n%           .\r\n%     151    -4.51944e-16   -9.2311e-47        interpolation\r\n%     152    -7.85458e-18  -4.84584e-52        interpolation\r\n%     153    -7.85458e-18  -4.84584e-52        bisection\r\n%     154    -7.85458e-18  -4.84584e-52        interpolation\r\n%\r\n%   Zero found in the interval [-0.5, 1]\r\n%   ans =\r\n%       -7.854580142952130e-18\r\n\r\n%% Muller's method\r\n% During the course of the iterations we often have three distinct values\r\n% of $x$, namely $a$, $b$, and $c$, and corresponding function values,\r\n% $f(a)$, $f(b)$, and $f(c)$.  We could interpolate these three values with\r\n% a quadratic in $x$, find the roots of this quadratic, and take the root that\r\n% is closest to our best approximation so far to be the next approximation\r\n% to the zero of $f(x)$.  This is Muller's method.  The trouble is the\r\n% roots of the quadratic might be complex.  This would be an advantage if\r\n% we were looking for complex solutions, and Muller's method is often used\r\n% in this situation.  But we don't want them here and Brent chose to do\r\n% something else.\r\n\r\n%% Inverse quadratic interpolation\r\n% Instead of a quadratic in $x$, we can interpolate the three points\r\n% with a quadratic function in $y$.  That's a ``sideways'' parabola,\r\n% $P(y)$, determined by the interpolation conditions\r\n%\r\n% $$ a = P(f(a)), \\ b = P(f(b)), \\ c = P(f(c)) $$\r\n%\r\n% This parabola always intersects the $x$-axis, which is $y = 0$.\r\n% So $x = P(0)$ is the next iterate.\r\n%\r\n% This method is known as _Inverse Quadratic Interpolation_,\r\n% abbreviated IQI.  Here is a naive implementation that illustrates\r\n% the idea.  It uses |polyinterp|, taken from\r\n% <https:\/\/www.mathworks.com\/moler.htmlncmfilelist.html % _Numerical Computing with MATLAB_>\r\n\r\ntype iqi\r\n\r\n%%\r\n% One difficulty with this simple code is that polynomial\r\n% interpolation requires the abscissae, which in this case are $f(a)$,\r\n% $f(b)$, and $f(c)$, to be distinct.  We have no guarantee that they\r\n% are.  For example, if we try to compute $\\sqrt{2}$ using $f(x) = x^2 - 2$\r\n% and start with $a = -2, b = 0, c = 2$, we are starting with $f(a) = f(c)$\r\n% and the first step is undefined.  If we start near this singular situation,\r\n% say with $a = -2.001, b = 0, c = 1.999$, the next iterate is near $x = 500$.\r\n\r\n%%\r\n% So IQI converges faster than secant once it is near the zero, but its\r\n% global behavior is unpredictable.  It needs to be combined with bisection\r\n% to have a reliable algorithm.\r\n\r\n%% Brent's algorithm\r\n% Like Dekker's original version,\r\n% Brent's version of _zeroin_ starts with an interval $[a,b]$ on which\r\n% the function $f(x)$ changes sign.  The objective\r\n% is to reduce the interval to a tiny subinterval on which the function\r\n% still changes sign.  The variables $a$, $b$, and $c$ play the same role:\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 four possible steps:\r\n%\r\n% * If $a \\ne c$ , try inverse quadratic interpolation.\r\n% * If $a = c$, try secant interpolation.\r\n% * If the interpolation step is near the endpoint, or outside the interval,\r\n% use bisection.\r\n% * If the step is smaller than the tolerance, use the tolerance.\r\n\r\n%%\r\n% Here is a somewhat cluttered plot that shows most of these possibilities.\r\n% The three values $a$, $b$, and $c$ are plotted in black.\r\n% The function\r\n%\r\n% $$ f(x) = x^3 - 3x - 2 $$\r\n%\r\n% is plotted in light grey to emphasize that the algorithm so far has only\r\n% evaluated it three times to give the three values shown by the circles.\r\n\r\nzeroin_plot(1,2.4)\r\n\r\n%%\r\n% The green line is the secant determined by $a$ and $b$.  A secant step\r\n% would go to the green $\\times$ labeled $t$, but that's outside the interval\r\n% so _zeroin_ would not go there.  Actually Brent's _zeroin_ would not even\r\n% be considering the secant in this situation because $a$ and $c$ are\r\n% distinct.\r\n\r\n%%\r\n% The blue curve is the quadratic in $y$ that goes through the three circles.\r\n% It crosses the $x$ -axis at the blue $\\times$ labeled $q$, so that would be\r\n% the IQI step.  However, in this case $q$ is more than three-quarters of the\r\n% way from $b$ to $c$, so Brent would choose bisection instead.\r\n\r\n%%\r\n% The red $\\times$ labeled $m$ is the bisection point halfway between $b$ and\r\n% $c$.  With a value of 2.025 it happens to be really close to the zero at 2.\r\n% Lucky shot.\r\n\r\n%%\r\n\r\n%% Fzero\r\n% Brent included Algol and Fortran implementations of this algorithm\r\n% in his book.  Forsythe, Mike Malcolm, and I made the Fortran program\r\n% the basis for the zero finding chapter of _Computer Methods for\r\n% Mathematical Computations_.  The code is\r\n% <http:\/\/www.netlib.org\/fmm\/zeroin.f still available from netlib>.\r\n\r\n%%\r\n% We named the MATLAB version |fzero|.   Here is the textbook version,\r\n% |fzerotx|, from <https:\/\/www.mathworks.com\/moler.htmlncmfilelist.html % _Numerical Computing with MATLAB_>.\r\n% This is not all of the MATLAB |fzero|.\r\n% In my next blog post I will deal with the need to know an interval on\r\n% which the function changes sign.\r\n\r\ntype fzerotx\r\n\r\n%% References\r\n% Richard P. Brent, \"An algorithm with guaranteed convergence for finding\r\n% a zero of a function,\" _Computer Journal_ 14, 422-25, 1971.\r\n\r\n%%\r\n% Richard P. Brent, _Algorithms for Minimization without Derivatives_,\r\n% Prentice-Hall, Englewood Cliffs, New Jersey, 1973, 195 pp.\r\n% Reprinted by Dover, New York, 2002 and 2013.\r\n% <http:\/\/maths-people.anu.edu.au\/~brent\/pub\/pub011.html % http:\/\/maths-people.anu.edu.au\/~brent\/pub\/pub011.html>\r\n\r\n%%\r\n% George Forsythe, Michael Malcolm, and Cleve Moler,\r\n% _Computer Methods for Mathematical Computations_,\r\n% Prentice-Hall, 259pp, 1977.\r\n\r\n##### SOURCE END ##### a8127df70bd8474491ad3dec2945b6ac\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/zeroin_blog2_01.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction-->Richard Brent's improvements to Dekker's <i>zeroin<\/i> algorithm, published in 1971, made it faster, safer in floating point arithmetic, and guaranteed not to fail.\r\n\r\n<!--\/introduction-->... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/10\/26\/zeroin-part-2-brents-version\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11,4,16,8],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1280"}],"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=1280"}],"version-history":[{"count":8,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1280\/revisions"}],"predecessor-version":[{"id":1439,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1280\/revisions\/1439"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=1280"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=1280"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=1280"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}