Having trouble with root finding problems that involve both bracketing and open methods.

The company you now work for just bought a horizontal cylindrical tank to store a liquid that your company sells. The expression relating the volume of the liquid, V , the radius of the tank, R, the length of the tank, L, and the height of the liquid in the tank, h, is given by

V = [R^2*cos^-1(R – h/R) − (R − h)*sqr(2Rh – h^2) ]*L

Your company’s tank has a radius of 3 feet and a length of 7 feet, and you currently have 140 cubic feet of liquid that needs to go into the tank.

Your boss wants you to implement an algorithm for finding the liquid height, h, given all the other parameters. Your boss, however, is old-school and doesn’t trust Matlab’s built-in wizardry (i.e., functions. So you “decide” to do the following.

1. Recast Equation 1 as a root finding problem. In other words, identify a function f(h) such that f(h) = 0 when h is the height of the liquid in your tank.

2. This is your first real engineering job, so to be sure that you get the right answer you use the bisection and secant functions. You also use a tolerance equal to machine precision. (Maybe overkill, but perhaps you can use this tidbit to impress your boss.)

3. Using both the bisection and the secant functions serves as double-check (i.e., verification) that you really have the right answer. To show your boss that you could have used Matlab to generate the same result, you also solve for the liquid height using Matlab’s fzero function. 4. You tabulate your results, showing (1) the liquid height, h, that you calculated for each of the three methods, (2) a calculation of liquid volume using the value of h that you calculated for each of the three methods (this is another verification check), and (3) the value of the function f(h) using the value of h from each of the three methods (yet another verification check). You are now really confident that you have the right answer!

Secant Method:

function [xr,k] = secant(f,xa,xb,tol,M, plotConv)

% bisection finds a root of a given function using a bracketing method

% Input

%   f...........the function

%   fp..........

%   x0

%   tol.........the tolerance

%   M...........the total number of iterations

%   plotConv....will generate plots if greater than zero

% Make sure that at least one root is bracketed

if (abs(f(xa)) < tol)

      warning('The value f(x0) is already within the search tolerance!')

      xr = xa;

      k = 0;

end

% Find a root!

xOld = xb;

xOlder = xa;

xp = linspace(0,5,100);

for k = 1:M

      % Use NR method

      xNew = xOld - f(xOld)*(xOld - xOlder)/(f(xOld) - f(xOlder));

      if (plotConv > 0),

          xp = linspace(min([0,xOld,xOlder,xNew]),max([5,xOld,xOlder,xNew]),100);

          plot(xp,f(xp))

          hold on

          plot([xNew xOld xOlder],[0 f(xOld) f(xOlder)],'r-o');

          grid on

          pause

          stem(xNew,f(xNew))

          hold off

      end

      % If xr is a root (within tolerance) then quit looking further

      if (abs(f(xNew)) < tol),

          break;

      else

          xOlder = xOld;

          xOld = xNew;

      end

      if (plotConv > 0), pause; end

end

% xr is either a root, or the closest we could get within the number of

% iterations that we performed

xr = xNew;

end

Bisection Method:

function [x,k] = bisection(f,a,b,tol,M, plotConv)

% bisection finds a root of a given function using a bracketing method

% Input

%   f...........the function

%   a...........the initial left side of the bracket

%   b...........the initial right side of the bracket

%   tol.........the tolerance

%   M...........the total number of iterations

%   plotConv....will generate plots if greater than zero

% Set the values of f(a) and f(b) to fa and fb, respectively, so we don't

% need multiple function calls (potentially expensive)

fa = f(a);

fb = f(b);

% Do some basic logic tests

% 1) Make sure that at least one root is bracketed

if (fa*fb > 0),

      error('The a and b values do not bracket a root!')

% 2) Make sure that a is not already a root

elseif (abs(fa) < tol)

      warning('The value f(a) is already within the search tolerance!')

      x = a;

      k = 0;

% 3) Make sure that b is not arleady a root

elseif (abs(fb) < tol)

      warning('The value f(b) is already within the search tolerance!')

      x = b;

      k = 0;

end

% Find a root!

for k = 1:M

      % Use bisection to approximately locate the root using the midpoint of

      % the bracket

      xr = a + (b-a)/2;

      fr = f(xr);

      if (plotConv > 0),

          if (mod(k,3)==1), xp = linspace(a,b,100); end

          plot(xp,f(xp))

          hold on

          stem(a,fa,'Color','red','LineWidth',2)

          stem(b,fb,'Color','red','LineWidth',2)

          stem(xr,fr,'Color','green','LineWidth',2)

          grid on

          hold off

      end

      % Reset the bracket

      if (fa*fr < 0),

          b = xr;

          fb = fr;

      else

          a = xr;

          fa = fr;

      end

      % If xr is a root (within tolerance) then quit looking further

      if (abs(fr) < tol),

          break;

      end

      if (plotConv > 0), pause; end

end

% xr is either a root, or the closest we could get within the number of

% iterations that we performed

x = xr;

end

Fzero Method

f (x) = e −ax cos(bx) + c can be defined as

>> a=0.1;

>> b=1.0;

>> c=0.5;

>> anonFunc = @(x) exp(-a*x).*cos(b*x)+c;

Variables a, b, and c in anonFunc are taken from the workspace. Then we can use

>> fzero(anonFunc, 2) % With anonymous functions, no @

>> fzero(@myFunc, 2) % With functions in a file, need @

[Geoff Hayes]

Steve - you've posted your homework assignment and a lot of code, but you have indicated what the problem is that you are experiencing besides stating that you are having trouble. Please clarify what you need help with.