Bracketing Methods in MATLAB | MATLAB Tutorial

MATLAB Helper Quiz

Roots problems occur when some function f can be written in terms of one or more dependent variables x, where the solutions to f(x) = 0 yields the solution to the problem. These problems often occur when a design problem presents an implicit equation for a required parameter.

Analytical solution to predict fall velocity as a function of time:

Try as you might, you cannot manipulate this equation to explicitly solve for m – that is, you cannot isolate the mass on the left side of the equation.

The alternative way of looking at the problem involves subtracting  v(t)  from both sides to give a new function:

Now we can see that the answer to the problem is the value of m that makes the function equal to zero. Hence, we call this a “roots” problem.

Bracketing

If we have a function of one variable, if we can determine that a solution to our problem lies on a particular interval [a, b], it seems reasonable that if we can select a point on that interval, say c, then if we can determine whether or not the solution lies in [a, c] or [c, b], then we have restricted the domain on which we know the solution lies. In that case, we can try again: shrink the new interval and get a tighter bound on the solution.

Bisection Method

The bisection method is a variation of the incremental search method in which the interval is always divided in half. If a function changes sign over an interval, the function value at the midpoint is evaluated. The location of the root is then determined as lying within the subinterval where the sign change occurs. The subinterval then becomes the interval for the next iteration. The process is repeated until the root is known to the required precision.

Use the graphical approach to determine the mass of the bungee jumper with a drag coefficient of 0.25 kg/m to have a velocity of 36 m/s after 4 s of free fall. The acceleration of gravity is 9.81

A graphical depiction of the method is provided in Fig 1. From fig 1, we can see that the function changes sign between values of 50 and 200. Therefore, the initial estimate of the root  lies at the midpoint of the interval.

Note that the exact value of the root is 142.7376. This means that the value of 125 calculated here has a true percent relative error of

Next we compute the product of the function value at the lower bound and at the midpoint:

which is greater than zero, and hence no sign change occurs between the lower bound and the midpoint. Consequently, the root must be located in the upper interval between 125 and 200. Therefore, we create a new interval by redefining the lower bound as 125.

Figure 1: A graphical depiction of the bisection method

At this point, the new interval extends from   to   . A revised root estimate can then be calculated as

which represents a true percent error of  . The process can be repeated to pbtain refined estimates. For example,





Therefore, the root is now in the lower interval between 125 and 162.5. The upper bound is redefined as 162.5, and the root estimate for the third iteration is calculated as

which represents a percent relative error of   . The method can be repeated until the result is accurate enough to satisfy your needs.

We ended Solution with the statement that the method could be continued to obtain a refined estimate of the root. We must now develop an objective criterion for deciding when to terminate the method.

Figure 2: Errors for the bisection method. True and approximate errors are plotted versus the number of iterations.

An initial suggestion might be to end the calculation when the error falls below some pre-specified level. For instance, in Solution, the true relative error dropped from 12.43 to 0.709% during the course of the computation. We might decide that we should terminate when the error drops below, say, 0.5%. This strategy is flawed because the error estimates in the example were based on knowledge of the true root of the function. This would not be the case in an actual situation because there would be no point in using the method if we already knew the root.

Therefore, we require an error estimate that is not contingent on foreknowledge of the root. One way to do this is by estimating an approximate percent relative error as in

where  is the root for the present iteration and  is the root from the previous iteration. When  becomes less than a prespecified stopping criterion , the computation is terminated.

The nature of the true error is due to the fact that, for bisection, the true root can lie anywhere within the bracketing interval. The true and approximate errors are far apart when the interval happens to be centred on the true root. They are close when the true root falls at either end of the interval.

Although the approximate error does not provide an exact estimate of the true error, Fig. 2 suggests that  captures the general downward trend of  . In addition, the plot exhibits the extremely attractive characteristic that  is always greater than . Thus, when  falls below  , the computation could be terminated with confidence that the root is known to be at least as accurate as the prespecified acceptable level.

function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
% Author - Gunjan Gupta, MATLAB Helper 
% Bisection Method -> Numerical Methods -> MATLABHelper.com
%% bisect: root location zeroes
% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
% uses bisection method to find the root of func

%% input:
% func = name of function
% xl, xu = lower and upper guesses
% es = desired relative error (default = 0.0001%)
% maxit = maximum allowable iterations (default = 50)
% p1,p2,... = additional parameters used by func

%% output:
% root = real root
% fx = function value at root
% ea = approximate relative error (%)
% iter = number of iterations
%%
if nargin < 3, error('At least 3 input arguments required'), end test = func(xl,varargin{:})*func(xu,varargin{:}); if test>0, error('no sign change'), end
if nargin<4 || isempty(es), es=0.0001; end
if nargin<5 || isempty(maxit), maxit=50; end
iter = 0; xr = xl; ea = 100;
while (1)
    xrold = xr;
    xr = (xl + xu)/2;
    iter = iter + 1;
    if xr ~= 0, ea = abs((xr - xrold)/xr) * 100; end
    test = func(xl,varargin{:})*func(xr,varargin{:});
    if test < 0 xu = xr; elseif test > 0
        xl = xr;
    else
        ea = 0;
    end
    if ea <= es || iter >= maxit, break, end
end
root = xr; fx = func(xr, varargin{:});

It is passed the function (func) along with lower (xl) and upper (xu) guesses. In addition, an optional stopping criterion (es) and maximum iterations (maxit) can be entered. The function first checks whether there are sufficient arguments and if the initial guesses bracket a sign change. If not, an error message is displayed and the function is terminated. It also assigns default values if maxit and es are not supplied. Then a while…break loop is employed to implement the bisection algorithm until the approximate error falls below es or the iterations exceed maxit. We need to find the root of

The bisect function can be used to determine roots as

fm=@(m) sqrt(9.81*m/0.25)*tanh(sqrt(9.81*0.25/m)*4)-36;
[mass,fx,ea,iter]= bisect(fm,40,200);
fprintf('Mass (Root)= %.2f \n',mass);
fprintf('Function Value at Root (fx) = %e \n',fx);
fprintf('Number of Interation (iter) = %d \n',iter);
fprintf('Approximate Relative Error (ea) = %e \n',ea);

Mass (Root)= 142.74
Function Value at Root (fx) = 4.608913e-07
Number of Interation (iter) = 21
Approximate Relative Error (ea) = 5.345047e-05

Thus, a result of m = 142.74 kg is obtained after 21 iterations with an approximate relative error of  0.00005345%, and a function value close to zero.

Leave a Comment