Know your tolerance

Given:

Contents

with the information given below, solve for the exit concentration. This should be simple. We use the units package to make sure we don't miss any units.

clear all; close all; clc
u = cmu.units;

Cao = 2*u.mol/u.L;
V = 10*u.L;
nu = 0.5*u.L/u.s;
k = 0.23 * u.L/u.mol/u.s;

function we need to solve

f = @(Ca) V - nu*(Cao - Ca)./(k*Ca.^2);

plot the function to estimate the solution

c = linspace(0,2)*u.mol/u.L;
plot(c/(u.mol/u.L),f(c))
ylim([-0.1 0.1]) % limit the y-range to see the zero easier
xlabel('C_A')
ylabel('f(C_A)')
% we want the value of C_A that makes f(C_A)=0. this happens around C_A=0.5

cguess = 0.5*u.mol/u.L;
c = fsolve(f,cguess)
% it is a little weird that the solution is the same as our initial guess.
% Usually the answer is at least a little different!

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.




ans =

  500.0000


Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.



500/m^3*mol

check your solution

z = f(c)
z.as(u.L)
% this is not really that close to zero, especially when you look at the
% units in L!

-0.00304348*m^3

ans =

-3.043*L

closer look at fsolve

Let's output some more information.

options = optimset('Display','iter');
[X,FVAL,EXITFLAG,OUTPUT] = fsolve(f,cguess,options)
% It should be a warning that zero iterations were taken! Our guess was
% probably not that good! The EXITFLAG=1 means the function ended ok as far
% as fsolve is concerned.

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          2    9.26276e-006                     1.85e-007               1

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.




ans =

  500.0000


                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          2    9.26276e-006                     1.85e-007               1

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.



500/m^3*mol
-0.00304348*m^3

EXITFLAG =

     1


OUTPUT = 

       iterations: 0
        funcCount: 2
        algorithm: 'trust-region dogleg'
    firstorderopt: 1.8526e-007
          message: [1x756 char]

The problem is bad scaling because the underlying units are stored in m^3, which results in very small numbers that are close to zero. So, fsolve sees a number it thinks is close to zero, and stops. One way to fix this is to decrease the tolerance on what constitutes zero.

options = optimset('TolFun',1e-9);
cguess = 0.5*u.mol/u.L;
[c2,FVAL,EXITFLAG,OUTPUT] = fsolve(f,cguess,options)
z2 = f(c2)
z2.as(u.L) % much closer to zero.

sprintf('The exit concentration is %1.2f mol/L', c2/(u.mol/u.L))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the selected value of the function tolerance, and
the problem appears regular as measured by the gradient.




ans =

  559.5545


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the selected value of the function tolerance, and
the problem appears regular as measured by the gradient.



559.554/m^3*mol
-1.24874e-006*m^3

EXITFLAG =

     1


OUTPUT = 

       iterations: 6
        funcCount: 14
        algorithm: 'trust-region dogleg'
    firstorderopt: 5.3309e-011
          message: [1x698 char]

-1.24874e-006*m^3

ans =

-0.001*L


ans =

The exit concentration is 0.56 mol/L

Rescaling the units

We can make the default length unit be dm, then the default volume unit will be 1L (dm^3). This will effectively rescale the function so that the default tolerance will be effective. Note you have to specify the unit for each unit in the order (length, time, mass, temperature, charge, mol).

u = cmu.unit.units('dm','s','kg','K','coul','mol');
Cao = 2*u.mol/u.L;
V = 10*u.L;
nu = 0.5*u.L/u.s;
k = 0.23 * u.L/u.mol/u.s;

% function we need to solve
f = @(Ca) V - nu*(Cao - Ca)./(k*Ca.^2);
cguess = 0.5*u.mol/u.L;
[c3,FVAL,EXITFLAG,OUTPUT] = fsolve(f,cguess)
sprintf('The exit concentration is %1.2f', c3)
% we were able to solve this more reasonably scaled problem without
% changing any default tolerances.

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.




ans =

    0.5596


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.



0.559584/dm^3*mol
-3.9293e-012*dm^3

EXITFLAG =

     1


OUTPUT = 

       iterations: 4
        funcCount: 10
        algorithm: 'trust-region dogleg'
    firstorderopt: 1.6772e-010
          message: [1x695 char]


ans =

The exit concentration is 0.56 (dm^3*mol)^-1

'done'
% categories: nonlinear algebra
% tags: reaction engineering, units

ans =

done

Published with MATLAB® 7.10