[Go to previous section: 1.4| Go
to next section: 1.6]
[ Go to CENG301 Homepage | CENG301 Notes Table of Contents]
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
>> 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
>> 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
>> 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