Matlab Ridder's Method Code

Question

I need to write a proper implementation of the Ridder's method in Matlab. I must define the function as

function [x sol, f at x sol, N iterations] = Ridders(f, x1, x2, eps f, eps x)

The explanation I was given is:

1. bracket the roots (x1, x2)

2. evaluate the midpoint (x1 + x2)/2

3. find new approximation for the root

   x4 = x3 + sign(f1 - f2) [f3/((f3)^2 - f1f2)^1/2)](x3 - x1)

4. check if x4 satisfies the convergence condition. if yes, stop. if not...

5. rebracket the root using x4 and whichever of x1, x2, or x3 is closer to the root

6. loop back to 1

I have no idea how to implement this in matlab. Help?

Answer 1

In order to find a zero of some function f(x) on a given interval [x0,x2], most implementations of the Ridders algorithm make use of the signal function in the iterative formula

x = x1 + signum(f0) * f1 * d / sqrt(f1^2 - f0 * f2)

where x1 = (x0 + x2)/2, d = |x2 - x0|/2 and the next x0 or x2 are substituted by x according to the signal of f(x). This formula was simplified by Ridders himself in his original paper to

x = x1 + (f1/f0) * d / sqrt((f1/f0)^2 - (f2/f0)).

Put in that way, the algorithm can be written in a more compact code.

function [x,e,it]=ridders(f,a,b,tol)
% x     solution
% e     approx. error in f(x)
% it    iterations done
% f     function
% a,b   interval bounds
% tol   tolerance
itmax = 200;
it = 0;
if f(a) < 0
    x0 = a;
    x2 = b;
else
    x0 = b;
    x2 = a;
end
if f(a) * f(b) > 0
    error(strcat('The function does not change signal on the extremes of the interval [',num2str(a),',',num2str(b),']. Redefine the bounds.'));
end
conve = false;
while (~conve)
    it = it + 1;
    x1 = (x0 + x2) / 2;
    d = abs(x2 - x0) / 2;
    f0 = f(x0);
    f1 = f(x1);
    f2 = f(x2);
    x = x1 + d * (f1/f0) / sqrt((f1/f0)^2 - f2/f0);
    e = abs(f(x));
    if f(x) > 0
        x2 = x;
    else
        x0 = x;
    end
    conve = (e<tol) | (it>itmax);
end

end

Let's test some functions:

• f = @(x) (x^3 - x - 5) on [-1, 3] results

  [x,e,it]=ridders(f,-1,3,1e-7)
  x = 1.9042
  e = 3.9639e-08
  it = 3
  

• g = @(x) (tan(x)^tan(x) - 1000) on [1.3, 1.4] results

  [x,e,it]=ridders(g,1.3,1.4,1e-7)
  x = 1.3547
  e = 4.6612e-12
  it = 5
  

Answer 2

I wrote a Matlab implementation of Ridder's method on the Matlab file exchange: submission 54458. I've copied the code below for reference:

function xZero = rootSolve(func,xLow,xUpp)
% XZERO = ROOTSOLVE(FUNC, XLOW, XUPP)
%
% FUNCTION: This function uses Ridder's Method to return a root, xZero,
%     of func on the interval [xLow,xUpp]
%
% INPUTS:
%   func = a function for a SISO function: y = f(x)
%   xLow = the lower search bound
%   xUpp = the upper search bound
%
% OUTPUTS:
%   xZero = the root of the function on the domain [xLow, xUpp]
%
% NOTES:
%   1) The function must be smooth
%   2) sign(f(xLow)) ~= sign(f(xUpp))
%   3) This function will return a root if one exists, and the function is
%   not crazy. If there are multiple roots, it will return the first one
%   that it finds.

maxIter = 50;
fLow = feval(func,xLow);
fUpp = feval(func,xUpp);
xZero = [];

tol = 10*eps;

if (fLow > 0.0 && fUpp < 0.0) || (fLow < 0.0 && fUpp > 0.0)
    for i=1:maxIter
        xMid = 0.5*(xLow+xUpp);
        fMid = feval(func,xMid);
        s = sqrt(fMid*fMid - fLow*fUpp);
        if s==0.0, break; end
        xTmp = (xMid-xLow)*fMid/s;
        if fLow >= fUpp
            xNew = xMid + xTmp;
        else
            xNew = xMid - xTmp;
        end
        xZero = xNew;
        fNew = feval(func,xZero);
        if abs(fNew)<tol, break; end

        %Update
        if sign(fMid) ~= sign(fNew)
            xLow = xMid;
            fLow = fMid;
            xUpp = xZero;
            fUpp = fNew;
        elseif sign(fLow) ~= sign(fNew)
            xUpp = xZero;
            fUpp = fNew;
        elseif sign(fUpp) ~= sign(fNew)
            xLow = xZero;
            fLow = fNew;
        else
            error('Something bad happened in riddersMethod!');
        end

    end
else
    if fLow == 0.0
        xZero = xLow;
    elseif fUpp == 0.0
        xZero = xUpp;
    else
        error('Root must be bracketed in Ridder''s Method!');
    end
end

Answer 3

Some principal concepts that might help:

• function handles (you need to provide f in this format)
• valid variable names (many variable names in your definition are not valid)

Answer 4

In Matlab, you would define your function as:

function [list of outputs] = myfunc(list of input variables)

    %function definition to compute outputs using the input variables.

In your case if you would like x4 (i.e root) to be your output, you would do:

function root  = riddler(func, x1, x2, xaccuracy, N)


    xl = x1;
    xh = x2;
    fl=func(x1)
    fh=func(x2)

   for i = 1:N
       xm = 0.5*(xl+xh);
       fm = func(xm);
       s = sqrt(fm*fm - fl*fh)
       if s == 0
         return;
       end
       xnew = xm + (xm - xl)*sign(fl - fh)*fm/s %update formula
       .
       .
       . % extra code to check convergence and assign final answer for root 
       . % code to update xl, xh, fl, fh, etc. (i.e rebracket)
       .

    end % end for