Page 1

Page 2

1.5.4 The Root Finding "Module": ssec2

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.

Example 1: Find one root of a cubic polynomial

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!```

Example 2: Find a root of sin(x)=0

```>> 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).

>> 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.

1.5.5 Extraction of a Solute with: phase

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
% 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.```

1.5.6 Extraction of a Solute with: phaseg

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
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```

Go back to Page 1

[Go to previous section: 1.4| Go to next section: 1.6]