How to find out root of a function

Question

The equation 𝑒^−0.2𝑥sin(𝑥+2) = 0.1 has three solutions in the interval 0<𝑥<10. Find these three solutions. (Answers: 𝑥 = 1.0187, 4.5334, 7.0066).

How I can get all roots? I used the fzero function like so:

fun = @(x) (exp(-0.2.*x).*sin(x+2));
x   = 0:.01:10;
z0  = fzero(fun, 0)

Answer 1

I'm going to provide a canonical answer to this question, because it is asked quite frequently in one way or another -- and not just here on Stack Overflow.

If you know beforehand that your function is smooth and continuous, you can use the method used by Stephen Morris in this File Exchange submission, which I reworked a little bit for efficiency and clarity.

His method approximates the function on the interval of interest by its Chebyshev polynomial. This gives a close-to-optimal approximation, with minimal function evaluations. The roots of this polynomial can be found easily with a method akin to MATLAB's own roots function.

Here is the reworked function:

% FINDREALROOTS     Find approximations to all real roots of any function 
%                   on an interval [a, b].
%
% USAGE:
%  Roots = FindRealRoots(funfcn, a, b, n, vectorized, make_plot)
%
% FINDREALROOTS() approximates all the real roots of the function 'funfcn' 
% in the interval [a,b]. It does so by finding the roots of an [n]-th degree 
% Chebyshev polynomial approximation, via the eigenvalues of the associated 
% companion matrix. 
%
% When the argument [vectorized] is [true], FINDREALROOTS() will evaluate 
% the function 'funfcn' at all [n] required points in the interval 
% simultaneously. Otherwise, it will use ARRAYFUN() to calculate the [n] 
% function values one-by-one. [vectorized] defaults to [false]. 
%
% When the argument [make_plot] is true, FINDREALROOTS() plots the 
% original function and the Chebyshev approximation, and shows any roots on
% the given interval. Also [make_plot] defaults to [false]. 
%
% All [Roots] (if any) will be sorted.
% 
% See also roots, eig.
function Roots = FindRealRoots(funfcn, a, b, n, vectorized, make_plot)

% First version 26th May 2007 by Stephen Morris, 
% Nightingale-EOS Ltd., St. Asaph, Wales.
%
% Modified 2017/January/02 by Rody Oldenhuis, LuxSpace sarl

    % parse input and initialize.
    inarg = nargin; 
    if n <= 2, n = 3; end                    % Minimum [n] is 3:    
    if (inarg < 5), vectorized = false; end  % default: function isn't vectorized
    if (inarg < 6), make_plot = false; end   % default: don't make plot

    % some convenient variables
    bma   = (b-a)/2;  
    bpa   = (b+a)/2; 
    Roots = [];

    % Obtain the Chebyshev coefficients for the function
    %
    % Based on the routine given in Numerical Recipes (3rd) section 5.8;
    % calculates the Chebyshev coefficients necessary to approximate some
    % function over the interval [a,b]

    % initialize 
    c = zeros(1,n);  
    k = (1:n)';  
    y = cos(pi*((1:n)-1/2)/n); 

    % evaluate function on Chebychev nodes
    if vectorized
        f = feval(funfcn, y*bma + bpa);
    else
        f = arrayfun(@(x) feval(funfcn,x), y*bma + bpa);
    end

    % compute the coefficients
    for jj = 1:n
        c(jj) = ( f(:).' * (cos((pi*(jj-1)) * ((k-0.5)/n))) ) * (2-(jj==1))/n; 
    end       

    % Coefficients may be [NaN] if [inf]
    if any(~isfinite(c(:))) || (c(n) == 0)
        return; end

    % Define [A] as the Frobenius-Chebyshev companion matrix. This is based
    % on the form given by J.P. Boyd, Appl. Num. Math. 56 pp.1077-1091 (2006).
    one          = ones(n-3,1);
    A            = diag([one/2; 0],-1) + diag([1; one/2],+1);
    A(end, :)    = -c(1:n-1)/2/c(n);
    A(end,end-1) = A(end,end-1) + 0.5;

    % Now we have the companion matrix, we can find its eigenvalues using the
    % MATLAB built-in function. We're only interested in the real elements of
    % the matrix:
    eigvals  = eig(A);  
    realvals = eigvals(imag(eigvals)==0);

    % if there aren't any real roots, return
    if isempty(realvals)
        return; end

    % Of course these are the roots scaled to the canonical interval [-1,1]. We
    % need to map them back onto the interval [a, b]; we widen the interval just
    % a tiny bit to make sure that we don't miss any that are right on the 
    % boundaries.
    rangevals = nonzeros(realvals(abs(realvals) <= 1+1e-5));

    % also sort the roots
    Roots = sort(rangevals*bma + bpa);

    % As a sanity check we'll plot out the original function and its Chebyshev
    % approximation: if they don't match then we know to call the routine again
    % with a larger 'n'.
    if make_plot

        % simple grid
        grid = linspace(a,b, max(25,n));

        % evaluate function like before
        if vectorized
            fungrid = feval(funfcn, grid);
        else
            fungrid = arrayfun(@(x) feval(funfcn,x), grid);
        end        
        % corresponding Chebychev-grid (more complicated but essentially
        % identical)
        y  = (2.*grid-a-b) ./ (b-a); 
        d  = zeros(1,length(grid)); 
        dd = d;
        for j = length(c):-1:2
            sv=d;
            d=(2*y.*d)-dd+c(j);
            dd=sv;
        end
        chebgrid = (y.*d) - dd + c(1);

        % Now make plot
        figure(1), clf,  hold on
        plot(grid, fungrid ,'color' , 'r');
        line(grid, chebgrid,'color' , 'b'); 
        line(grid, zeros(1,length(grid)), 'linestyle','--')
        legend('function', 'interpolation')

    end % make plot

end % FindRealRoots

The following script uses this function to solve your problem. Since your function is rather trivial to take derivatives of, I'm going to show how fzero compares to Newton's method, Halley's method and the 4^th order Householder method:

% Parameters
N    = 9;  % Degree of the Chebyshev polynomial to use
xmin = 0;  % LB
xmax = 10; % UB

% Function of interest, and its first 3 derivatives
f    = @(x) +exp(-0.2*x) .*       sin(x+2)                  - 0.1;
fp   = @(x) +exp(-0.2*x) .* (     cos(x+2)  - 0.200*sin(x+2));
fpp  = @(x) -exp(-0.2*x) .* ( 0.40*cos(x+2) + 0.960*sin(x+2));
fppp = @(x) +exp(-0.2*x) .* (-0.88*cos(x+2) + 0.592*sin(x+2));

% Approximate the roots 
R = FindRealRoots(f, xmin, xmax, N, true, false);

% Refine the roots
% Compare fzero against the first 3 Householder methods
A = R;  fevals_fzero  = N;
B = R;  fevals_Newton = N;
C = R;  fevals_Halley = N;
D = R;  fevals_4th    = N;

options = optimset('TolX'   , 1e-10,...
                   'display', 'off');

for ii = 1:numel(A)

    % FZERO
    [A(ii), ~,~, out] = fzero(f, R(ii), options);       
    fevals_fzero = fevals_fzero + out.funcCount;

    % Newton's method
    xprev = inf;
    x     = R(ii);    
    while ( abs(xprev - x) > 1e-10 )
        xprev = x;
        x     = x - f(x)/fp(x);
        fevals_Newton = fevals_Newton + 2;
    end    
    B(ii) = x;

    % Halley's method
    xprev = inf;
    x     = R(ii);    
    while ( abs(xprev - x) > 1e-10 )

        fx   = f(x);
        fpx  = fp(x);
        fppx = fpp(x);

        xprev = x;
        x     = x - 2*fx*fpx/(2*fpx*fpx - fx*fppx);

        fevals_Halley = fevals_Halley + 3;
    end    
    C(ii) = x;

    % 4th order method 
    xprev = inf;
    x     = R(ii);    
    while ( abs(xprev - x) > 1e-10 )

        fx    = f(x);    fx2  = fx*fx;
        fpx   = fp(x);   fpx2 = fpx*fpx;
        fppx  = fpp(x);
        fpppx = fppp(x);

        xprev = x;
        x     = x -  3*(2*fx*fpx2 - fx2*fppx) / ...
                    (6*(fpx2*fpx - fx*fpx*fppx) + fx2*fpppx);

        fevals_4th = fevals_4th + 4;
    end    
    D(ii) = x;

end

% Check for convergence
all_equal = all(abs(A-B) < 1e-8) && ...
            all(abs(A-C) < 1e-8) && ...
            all(abs(A-D) < 1e-8)

% function evaluations per method
fevals_fzero 
fevals_Newton
fevals_Halley
fevals_4th 

Figure generated when the last argument to FindRealRoots is set to true:

Command line output:

all_equal =
 1
fevals_fzero =
    29
fevals_Newton =
    27
fevals_Halley =
    36
fevals_4th =
    41

Conclusions:

• For smooth, continuous functions, the Chebyshev polynomial is quite efficient and suffices for most intents and purposes. The example here used only 9 function evaluations to find (1) the number of roots, and (2) initial estimates for each to about 2 digits accuracy. These numbers only improve for increased polynomial order. For this particular function, the accuracy improves with less computational cost than any of the more "standard" refinement methods demonstrated here -- setting the polynomial order to about 19 (implying 19 function evaluations) gives comparable accuracy and performance, so no refinement is actually needed.
• fzero is easy to use, but certainly not the best under all circumstances. In this particular case, Newton's method is trivial to set up, and has better performance.
• Don't rely too heavily on statements on the rate of convergence of a numerical method; that is a theoretical limit that does not by itself prove the superiority of a method used in practice. It is shown in the above that the higher the order of the refinement method, the overall performance gets worse.

Answer 2

As a general numerical problem, this is indeed hard.

You need to make some assumptions to help fsolve() (implicit in M_Tornack’s answer). For example, assume that a coarse grid search of some tolerance finds all zero crossings. The tolerance defines the minimum separation of the zeroes. Grid search is slow and inexact, so therefore we refine the answers with fsolve afterwards. As follows:

fun = @(x) (exp(-0.2.*x).*sin(x+2))-0.1;
x = 0:0.1:10; % interval start : minimum separation of zeroes : end
zeros_approx = x(find(diff(fun(x)>0))) % grid search
for i = 1:length(zeros_approx)
    zeros_refined(i) = fsolve(fun,zeros_approx(i))
end

Depending on your problem, you might want to check your boundaries as well for zero solutions. Another way is to perform grid search iteratively, that is, to start with a coarse grid and then refine around approximate zeroes.

Answer 3

A good results is also given in: https://es.mathworks.com/matlabcentral/answers/143643-fzero-function-calculating-all-zeros-within-interval

The solution is composed in 3 steps:

1. step: multiply function with itself shifted by 1 element to get zero crossings and get negative values

   fun =@(x) exp(-0.2.*x).*sin(x+2) - 0.1;
   x = linspace(0,10);
   y = fun(x);
   zx = x(fun(x).*circshift(fun(x),[0 -1]) <= 0);  % Estimate zero crossings
   

2. step: Delete last entry

   zx = zx(1:end-1);  % Eliminate any due to ‘wrap-around’ effect
   

3. step:calculate all roots at the zero crossings

   for k1 = 1:length(zx)
       fz(k1) = fzero(fun, zx(k1));
   end
   

Answer 4

The function fzero stops when it finds the first zero. In the way you coded, this zero will be the closest one to '0' (second parameter).

You will find these links useful, as the questions are very similar. The solution is to loop around the function fzero.

https://www.mathworks.com/matlabcentral/answers/103169-how-do-i-find-all-the-zeros-of-a-function

https://www.mathworks.com/matlabcentral/answers/18398-getting-3-zeros-from-a-function-using-fzero