[ Go to CENG301 Homepage | CENG301 Notes Table of Contents]
[Go to previous section: 1.4| Go
to next section: 1.6]
[ Go to CENG301 Homepage
| CENG301 Notes Table of Contents]
Page 2
The root finding functions ssec1 and
ssec2 were developed for a course in computer
applications to illustrate the ease of producing a general purpose
function to find a scalar root of one nonlinear equation. These
functions use the secant method discussed in Reklaitis on page 285 to
find a root in a specified interval. The functions will be found to
be very useful in this course since many of our equations require
trial and error solution for a single unknown.
Here is what help shows for ssec2:
>> help ssec2
function x=ssec2(xs,fx,c,erc,nmax) finds linear approximation of f(x)=0
xs gives an interval to seek a solution to f(x) = 0 in.
fx is a character vector that tells how to find f(x).
c is a vector of parameters in f(x)
erc gives the difference between successive values of x allowed
when convergence is deemed to be achieved.
nmax is the maximum number of trials, before quitting.
If erc and nmax are not given, then the program uses:
erc=.0001 and nmax=10
If only two arguments are given, c is set to null.
The program uses the function ssec1 for each iteration.
Example 1: >> cs=[1 -2 -4.25 7.5];
>> disp(ssec2([2 3],'polyval(c,x)',cs,.0001,10))
Example 2: >> disp(ssec2([3 4],'sin(x)'))
The program uses the function ssec1 for each
iteration.
Here is the first example suggested in the help listing for
ssec2:
>> cs=[1 -2 -4.25 7.5];
>> disp(ssec2([2 3],'polyval(c,x)',cs,.0001,10))
2.5000 <-- One root of the cubic.
>> polyval(cs,2.5)
ans =
0 <-- It checks!
>> disp(ssec2([3 4],'sin(x)'))
3.1416 <-- Very close to pi
Here is another example that shows how to find multiple roots with
ssec2
Example 3: Find all roots of a quadratic:
f(x) = -0.451*x2 + 0.2836*x + 0.5641
By using the ssec2 function:
>> ssec2([0 1],'polyval(c,x)',[-.451 .2836 .5641])
ans =
1.4761
or by:
>> ssec2([0 1],'polyval([-.451 .2836 .5641],x)')
ans =
1.4761
To find the second root, we give it a different interval:
>> ssec2([-1 0],'polyval(c,x)',[-.451 .2836 .5641]> ans = -0.8473
or:
>> ssec2([-1 0],'polyval([-.451 .2836 .5641],x)') ans = -0.8473
Since we are dealing with a polynomial, it is easier to find all the
roots at once with roots:
>> help roots
ROOTS Find polynomial roots.
ROOTS(C) computes the roots of the polynomial whose coefficients
are the elements of the vector C. If C has N+1 components,
the polynomial is C(1)*X^N + ... + C(N)*X + C(N+1).
See also POLY.
>> roots([-.451 .2836 .5641])
ans =
1.4761
-0.8473
NOTE: Sometimes SSEC2 will not return a satisfactory root,
but that does not mean that there
are no roots. Often the interval chosen does not
bracket a root you want to see or too few trials are done. In either
case, choose a new interval or increase the number of trials. If a
root is still not found, plot the function to see what is wrong or
use another root finding function such as ROOTS, ROOT3, FZERO, or
POLY.
The solution of Example Problem 5.5 in Reklaitis with these functions
is trivial. We want a root of the cubic equation:
f(x) = x - 2(1-x)3 = 0
If we look in the interval: 0 < x < 1, with
ssec2 we find:
>> ssec2(0:1,'x-2*(1-x)^3')
ans =
0.4102 <----this is the root
We can easily adapt this procedure to solve many of the trial and
error problems in process calculations. For example, we can get back
to finishing Example Problem 2.8
from the text. It is convenient to define a function that will act as
our f(x) for this problem. Er28 does this:
function er=er28(x) ti=sept(1,x,[.83333,0,.16667],5); sep(ti,5,[4 3]) er=.25-ns(3,1)./ns(3,2);
We need to go through the process of setting ns as we did up
to the point where we discovered that we needed a trial and error
solution. All we really need is for the feed stream (5) to be set and
the shape of ns to be correct. We will also need to have the
right cnms to use showm. Recall that the
ratio of benzene (component 1) to the other HC (component 2) in the
extract stream (stream 3) must be 0.25; hence the form of
er. If we execute er28 with our
first assumed value for the mass fraction of benzene in stream 3:
0.1, we get:
>> er28(.1) ans = -0.9470
While an initial guess of 0 gives:
>> er28(0)
ans =
0.2500
Obviously the answer must be between 0 and 0.1 so we can use
ssec2 to find it:
>> ssec2([0 .1],'er28(x)')
ans =
0.0231
We also have set our flow matrix to give the desired result:
>> showm(5 ,[4 3],10,2) Compound Inlet | Outlet Stream 5 | 4 3 Total benzene 700.00 | 625.00 75.00 700.00 other hc 300.00 | 0.00 300.00 300.00 SO2 3000.00 | 125.00 2875.00 3000.00 Total 4000.00 | 750.00 3250.00 4000.00
The problem asks for the percent recovery of benzene - the amount of
benzene in the raffinate (stream 4) to the amount of benzene in the
refinery stream (stream 1) is:
>> ns(4,1)/ns(1,1)
ans =
0.8929
Therefore, 89.3 percent of the benzene is recovered.
As a final example in the chapter, we will consider extracting a
solute out of an "oil" mixture with "water". One formulation of this
type of problem is in terms of two feed streams containing:
Stream Water Oil Solute
1 W 0.0 S1
2 0.0 O S2
|
|
If we mix these streams and separate the result, we can produce
two products: a light stream (of mainly oil and solute) and a heavy
stream (of mainly water and solute.) In the ideal case the "oil"
stream contains no water and vice versa. We may also be able to find
data that determines the distribution of the solute between them. One
approximation to the distribution is given by:
with D a known distribution coefficient. We can actually solve for
all flow rates in such a separation by solving a quadratic equation
involving the amount of solute in one of the phases. We will
demonstrate another way to do this by using our modular functions and
ssec2. First set globals and initialize
cnms and ns with
start301,choice (3) Other :
You have chosen to not use any compound data If you wish to enter data for a compound, exit (control-C) then execute "create", then run a "new session", choice (1) Enter the number of compounds: 3 Compound names must be less than 25 characters long! Give the name of compound # 1: water Compound names must be less than 25 characters long! Give the name of compound # 2: oil Compound names must be less than 25 characters long! Give the name of compound # 3: solute Enter the number of streams: 5
Set the feed streams:
>> ns(1:2,:)=[100 0 2;0 50 10]; >> mix(1:2,3) >> showm(1:2,3,10,2) Compound Inlet Outlet Stream 1 2 Total 3 water 100.00 0.00 100.00 100.00 oil 0.00 50.00 50.00 50.00 solute 2.00 10.00 12.00 12.00 Total 102.00 60.00 162.00 162.00
Then write a program to try various fractions of the solute sent to
the oil phase:
function out=phase(dis,x)
% Phase: calculates solute distribution in a two-phase separation.
% function out=phase(dis,x)
% Argument List:
% Argument Gives
% dis distribution coefficient
% x fraction of solute
% going into the top phase
% phase tries various fractions of solutes sent to a phase
% and returns the difference between the calculated and desired
% distribution coefficients
% See also PHASEG, PICPHASEG
% OKB, TYLC
global ns
out=[];
for k=1:length(x)
sep([0 1,x(k)],3,[4 5]);% [No water in 4, All oil in 4, x(k) of Solute in 4]
temp=sum(ns(4:5,:)');% total flow in each stream out of separator: row
xs=ns(4:5,3)./temp';% fraction solute in streams 4 and 5
subt=xs(1)/xs(2);% ratio of the fractions of solute in 4 to 5
if subt==NaN % check to see if the fractions were both zero
subt=0;
end
out=[out,dis-subt];% we want subt to equal dis
end
Note that the first argument of phase is the
distribution coefficient and the second argument is the fraction of
the solute going into the top phase. We will assume that in our case
D = 3, and for this value try a couple of values for the fraction to
send to the top:
>> phase(3,.9) <-- D = 3 and the frac. to top = .9
ans =
-11.9803
>> phase(3,.1) <-- D = 3 and the frac. = .1
ans =
2.7595
OK, there must be a solution for D = 3 somewhere between 0.1 and
0.9.
>> ssec2([.1 .9],'phase(3,x)')
ans =
0.6226
There it is. We can display our results and check them.
>> showm(3,4:5,11,4)
Compound Inlet Outlet
Stream 3 4 5 Total
water 100.0000 0.0000 100.0000 100.0000
oil 50.0000 50.0000 0.0000 50.0000
solute 12.0000 7.4708 4.5292 12.0000
Total 162.0000 57.4708 104.5292 162.0000
>> x45=[7.4708 4.5292]./[57.4708 104.5292]
x45 =
0.1300 0.0433
>> .13/.0433
ans =
3.0023 <-- Pretty close to the specification.
A more generalized version of phase is
phaseg. Help on it shows:
>> help phaseg
Phaseg: Calculates solute distribution in 2-phase separation using modules.
function er=phaseg(iw,ioil,isol,dis,in,out,x)
Argument Gives
iw index of the water in cnms.
ioil index of oil in cnms.
isol index of the solute in cnms.
dis distribution coefficient
in index of the inlet stream
out indices of the exit streams
x fraction of solute going into
the top phase
er returns the difference between dis and the
ratio actually set for the solute in the exit streams.
For additional help and picture see PICPHASEG
See also PHASE
Example: >> ssec2([0.1 0.9],'phaseg(1,2,3,3.0,3,[4 5],x)');
revised help by TYLC
Here is the image of phaseg:
Here is the same simulation that we showed for phase with D=3.0 and a guess of 0.1 for the fraction in the top stream.
>> phaseg(1,2,3,3.0,3,[4 5],0.1) <-- D = 3 and the frac. = .1
ans =
2.7595
It may be used in ssec2 just like phase
>> ssec2([.1 .9],'phaseg(1,2,3,3.0,3,4:5,x)')
ans =
0.6226