Roots Finding: Bracketing Methods

Problem 5.14.

A reversible chemical reaction 2A + B ⇄ C can be characterized by the equilibrium relationship 

 Eq.(5.14(a))

where the nomenclature c_i represents the concentration of constituent i. Suppose that we define a variable x as representing the number of moles of C that are produced. Conservation of mass can be used to reformulate the equilibrium relationship as 

 Eq.(5.14(b))

where subscript 0 designates the initial concentration of each constituent. If K = 0.016, c_a,0 = 42, c_b,0 = 28, and c_c,0 = 4, determine the value of x.

(a) Obtain the solution graphically.

(b) On the basis of (a), solve for the root with initial guesses of x_l = 0 and x_u = 20 to ε_s = 0.5%. Choose either bisection or false position to obtain your solution. Justify your choice. 

Solution:

(a) Firstly, substitute the parameters into the Eq.(5.14(b)) to obtain a single function in terms of x:

MATLAB can be used to plot the function:

>> x = linspace(0,20);

>> f = 0.016-(4+x)./(((42-2*x).^2).*(28-x));

>> plot(x,f),grid,xlabel('x'),ylabel('f(x)')

A visual inspection of the plot gives a rough estimate of the root, which is about 15.9 moles. The validity of the graphical estimate can be checked by substituting it into f(x) to yield:

>> f = 0.016-(4+15.9)./((42-2*15.9).^2.*(28-15.9))

f =

   1.9235e-04

which is close to zero. 

(b) In this case, I will choose bisection to solve the problem. This is because false-position is not suitable to find root for a function with significant curvature, as it would lead to poor and slow convergence. We can use bisect M-file below to solve for x:

function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)

%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

% varargin = variable number of arguments, can input any number of

% arguments in the 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{:});

Plug in the given initial guesses and stopping criterion:

>> f = @(x)0.016-(4+x)./((42-2*x).^2.*(28-x));

>> [root,fx,ea,iter]=bisect(f,0,20,0.5)

root =

   15.8594

fx =

   5.2493e-04

ea =

    0.4926

iter =

     8

Thus, after eight iterations, the x is computed as 15.8594 moles with an approximate relative error of 0.4926%. 

**Additional: Or else we can also use the MATLAB fzero function:

>> f = @(x)0.016-(4+x)./((42-2*x).^2.*(28-x));

>> x = fzero(f,[0,20])

x =

   15.9230

which gives the exact value of 15.9230 moles.