Section 3: Steady State Systems and Algebraic Equations

This section discusses algebraic equations and how they may be solved. It starts of by showing how one single equation may be solved iteratively, and then shows how this may be extended to sets of equations.


Section 3.1: General Introduction to Algebraic Equations

Section 3.1: General Introduction to Algebraic Equations

We will consider a particular case where sets of algebraic equations may be solved very simply. However, we will first discuss briefly the case of single algebraic equations.

Single Algebraic Equations

Consider the single equation:
x -3 = 0

We solve this equation by rearranging it to the form:
x = 3

This is also an equation, because it contains an `=' sign, but it is additionally a formula because it is written so as to enable an unknown quantity, here x to be calculated explicitly. The original form of the equation did not have this property.

Not all single equations in single unknowns may be easily rearranged to provide formulas. For example:
x² - 3x - 4 = 0
can be rewritten as a formula, but not particularly easily, although the form is a standard one. It also yields more than one possible value for x which can sometimes complicate matters.

x³ - x² - 3x - 4 = 0
can be rewritten only as a complicated formula, and so is seldom solved in this way.

x - log x = 0
This cannot be rewritten as a formula and so must be solved by a trial and error method.

In general, we can rewrite an algebraic equation as a formula only when one of the following is true:

Programming languages, including spreadsheet programs, are able to evaluate formulas and so can solve single algebraic equations which can be suitably rearranged. In this introduction to modelling we will not consider problems where equations cannot be rearranged in this way.

Sets of Equations Suitable for `Direct' Solution

Consider the pair of equations:
x - 3 = 0
y - 4 x + 1 = 0

This is a set of two equations in two unknowns ( x and y) and can therefore, in principle, be solved for both its unknowns.

For any set of equations to be solvable, there must be exactly as many unknowns as equations.

It is also fairly obvious how the above equations can be solved. The first can be rearranged to give a formula for x:
x = 3

The second can be rearranged to give a formula for y, so that when x is known, y can also be calculated:
y = 4 x - 1

The multiequation/multivariable equivalent of the condition for solving a single equation by rearrangement to a formula, which we will describe as direct solution , is in fact two conditions. The first of these should be obvious from the above example.

For direct solution, it must be possible to rearrange the equations so that a formula can be written for every unknown.

The second condition is less obvious but may be illustrated by a case where it is not satisfied. consider the two equations:
y - 4 x + 1 = 0
- = 0

Since neither equation contains either x or y on its own, it is not possible to solve immediately for either of the unknowns. To be able to do this it would be necessary for one equation to contain only one of the unknowns, it would not matter which, and the other could contain both. In general, a set of n equations to be solved this way must be triangular , i.e. one must contain no more than one unknown, another no more than two, another no more than 3 and so forth.

The occurrence of these must also be in the correct pattern so that it is possible for the equations to be both rearranged and reordered so that they can be written, for example, for a set of n equations in unknowns x1 to xn:

x1 = constant expression
x2 = f(x1, constants)
x3 = f(x1, x2, constants)
x4 = f(x1, x2, x3, constants)
.
.
.
xn = f(x1 ... xn, constants)

There are a number of computer algorithms for testing and reordering sets of equations for solution in this way. Most modern spreadsheet systems will employ them automatically. There are also algorithms which will carry out the rearrangement, but these are not widely available. When using a simple program or spreadsheet the user will thus need to carry out the rearrangement, which normally will involve checking that a suitable order for direct solution does in fact exist.

General Methods for Sets of Algebraic Equations

Many, indeed probably most, serious modelling problems give rise to sets of equations which cannot be solved by the above methods and so iterative, trial-and-error solution is required, usually for several, sometimes for all of the unknowns. Discussion of these methods is beyond the scope of this introductory course, but modern computers in fact have little difficulty in handling very large sets of algebraic equations.

As two special cases, we have provided in the modelling laboratory a simple iterative solver for single algebraic equations (L1.1) and a standard solver for sets of linear equations (L1.2).


Section 3.2: Simple Ordered Sets of Equations

Section 3.2: Simple Ordered Sets of Equations

As has been shown, there are many cases when a set of algebraic equations in several unknowns can can be reordered and/or rearranged to enable all variables to be solved for in sequence rather than simultaneously.

We can create an incidence or occurrence table (or matrix) for the set of equations, see example (1.1.2.2) with a row for each equation and a column for each variable.

By trial and error, it is possible to try rearranging the rows and columns to see if a `lower triangular' structure can be arrived at, with all the non-zero elements (i.e. x's in the table) arranged into the bottom left hand half of the table.

Here are some examples to try (section 3.2.1).

If you are using a modern spreadsheet to solve a problem involving algebraic equations of this sort, then the spreadsheet will probably sort the equations out in this way using one of the standard algorithms which have been developed for this task. You will still have to perform some analysis of your equations to ensure that the individual equations are written, after rearrangement if necessary:

The last of these is quite a subtle requirement. If more than one rearrangement appears to be possible, use the incidence matrix to help determine an order.

If a correct order has not been determined, or does not exist, then the spreadsheet system will probably complain and give an error message about `circular references'.

If a lower triangular arrangement of the equations is not possible, then a iterative trail-and-error procedure must be used. This will be discussed after we have consider iterative solution of single equations in the next section.


Section 3.2.1: Ordering Equations

Section 3.2.1: Ordering Equations - Initial Examples

The sets of equations given below can all be solved directly after suitable reordering and rearrangement. For each set of equations:


Section 3.2.2: Ordering Equations Solutions

Section 3.2.2: Ordering Equations Solutions

1.
Unknowns: x1, x2, t, p
There are four unknowns and four equations, therefore the set may be solved.

Incidence Table
Equationx1x2tp
(1)11
(2)11
(3)1
(4)1111

Rearranging and reordering gives the following incidence matrix:
Equationtx1x2p
(3)1
(2)11
(1)11
(4)1111

This matrix is indeed lower triangular and so the equations may be solved in the following order:

(3) t = 4/a
(2) x1 = -3/(t + t2)
(1) x2 = 1 - x1
(4) p = (-t x2)/x1


Section 3.3: Iterative Solution of Single Algebraic Equations

Section 3.3: Iterative Solution of Single Algebraic Equations

Introduction

To solve a single equation analytically it must be possible to rewrite it as a formula with its unknown on the left hand side. As already noted, this is in general possible only if either: A very few other special equations, such quadratics and cubics, have analytical solutions. These do not, however, generally appear in physical systems and so they are not treated here.

In practice there are also very few interesting problem in chemical engineering which can conveniently be turned into just one equation. However, we will later see how the methods used for a single equation can often be extended to more realistic problems described by quite large numbers of equations by analysing the structure of the set of equations. This uses essentially the same technique we have just used to solve certain sets of equations `directly' by rearrangement.

Rather than looking for trivial chemical engineering problems with single equations we shall deal with entirely abstract examples. Consider our single equation to involve a single unknown x, and may be written as:

f( x ) = 0

Here f implies any functional form which can be expressed using any of the algebraic operators with which we are already familiar.

We will also assume that our equations have a single solution at one value of x. In principle we know that this is not the case, e.g. a quadratic has two solutions, a cubic has three, etc. In practice, however, it is unusual for an engineering problem to have multiple meaningful solutions. There may be additional solutions which are mathematically valid, but they usually have values which are physically unacceptable, e.g. at mole fractions greater than one, or negative Kelvin temperatures. (Note: exceptions do occur, but they are outside the scope of this course.) Since the methods used to find the solution normally involve searching within predefined bounds, these physically unrealistic solutions may often be excluded.

Methods

The notes below describe three methods used to solve single algebraic equations:

Single Equation Solvers

Most spreadsheet systems now provide a solve or `goal seek' facility which can be programmed to adjust a cell, containing the value of the unknown (e.g. x) until another cell containing a formula involving the unknown (e.g. f( x )) is reduced to zero. The methods used in different spreadsheets vary in both power and sophistication. Try some of the examples below.

A WWW based solver using bisection and written in JavaScript is available here. this has a single step facility and it is worth using this on a few examples to watch the progress towards the solution.

Example Equations

  1. x - 3 = 0
    A linear equation is actually a special type of algebraic equation.

  2. x2 - x -12 = 0
    There are two solutions, one positive and one negative. Run the solver with different initial bounds (e.g. 0 to 12 and -12 to 0) and find both solutions separately. What happens if a set of bounds encompassing both solutions is given?

  3. x + log (x) = 0
    The solution lies between 0 and 1.0 but log(0) of course cannot be calculated.

  4. x² - e-x = 0 ; solution is between -1 and +1

  5. 1 - exp(-x) - exp(-2x) = 0 ; solution is between -1 and +1

  6. 1 - (p1 + p2) = 0
    where p1 = exp(11(1-300/x))
    and p2 = exp(11(1-250/x))
    Solution is between 100 and 400.

Algebraic Equations

Algebraic Equations

3.3.1 Solving Algebraic Equations

3.3.2 Bisection Method

3.3.3.Secant Method

3.3.4 Newton's Method

 

The following are all algebraic equations.

x - 3 = 5 (1)

sin 3x = y cos x (2)


a x2 + b x + c = 0 (3)

Analytical Solution

Only equations which are linear in their unknown have a guaranteed analytical solution:

a x = b (4)

Or:

a x - b =0 (5)


\begin{displaymath}x = \frac{b}{a} \end{displaymath}

However, consider:

\begin{displaymath}\log x - 4 = 3 \end{displaymath} (6)


\begin{displaymath}\sin(3\pi y) + 4 - \cos \frac{3 \pi}{2} = 0 \end{displaymath} (7)

Solving a Single Algebraic Equation

Where is the solution ?






Here:

It clearly lies on the x-axis at the point where the function f(x) value is zero.

Given an equation

f(x) = 0

find the value of x such that the function f(x) is zero to within some specified tolerance.

Engineering Problems have Bounded Solutions

Intrinsic:

$ 0 \le mol \; fraction \le 1.0 $

Operational: cooling water temperature, vessel pressure, etc

3.3.2 Simplest Method: Bisection

First bisection gives a new, better x

Next discard redundant bound, here xmin

Replace it by xm which becomes `new' xmin

Repeat

Bisection Algorithm

1.
Specify bounds xmin and xmax
2.
Evaluate fl=f(xmin) and fu=f(xmax)
3.
Check that sign $f_l \neq $ sign fu. Abandon with error message if not so.
4.
Start iteration.
5.
Let xm = (xmin+xmax)/2
6.
Evaluate fm = f(xm)
7.
If sign fm = sign fl then make xmin = xm and fl=fm
8.
Else make xmax = xm and fu = fm
9.
If $ f_m \approx 0 $ then stop
10.
Repeat from step 4

Accuracy:

\begin{displaymath}r = \frac{1}{2^n} \end{displaymath}

3.3.3 Secant Method

Method 2 - Secant

Same algorithm as bisection, but a different way of choosing the new value.

Secant Update

If f(x) is nearly linear between xmin and xmax, then a straight line between the point (xmin , fl) and (xmax , fu) will cross the x-axis at xs , close to the true solution of f(x)=0.

From similar triangles A and B:

\begin{displaymath}\frac{x_{max}-x_s}{-f_u} = \frac{x_s-x_{min}}{f_l} \end{displaymath}

\begin{displaymath}x_s = \frac{f_u x_{min} - f_l x_{max}}{f_u - f_l}
\end{displaymath}

New step 4:
1.
Let xm= (fu xmin - fl xmax)/(fu - fl)

3.3.4 Newton's Method

Method 3 - Newton-Raphson

Bisection uses only sign,

Secant uses value,

Newton also uses slope.

Newton Update


\begin{displaymath}\frac{f(x_0)}{x_0 - x_1} = S = f'(x_0)
\end{displaymath}

Hence:

\begin{displaymath}x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}
\end{displaymath}

Possible problems:

Numerical Differentiation

By perturbing x by a small amount $\delta x$


\begin{displaymath}f'(x) = \frac{df}{dx} \approx \frac{f(x+\delta x) - f(x)}{\delta x}
\end{displaymath}

Then:

   \begin{displaymath}x_1 = x_0 - \frac{f(x_0) \delta x}{f(x+\delta x) - f(x)}
\end{displaymath}

<

NB, if we set:

x0 = xmin


\begin{displaymath}x_0+\delta x = x_{max}\end{displaymath}

\begin{displaymath}\delta x = x_{max} - x_{min}\end{displaymath}

...this gets back to the secant formula.

A Method which Usually Doesn't Work

Methods involving rearrangement of:
f(x)=0

to the form:
   \begin{displaymath}x = \phi(x)
\end{displaymath}

E.g. consider:

x3 +b x + c = 0

By alternative rearrangements:
\begin{displaymath}x = - ( b x + c) ^{\frac{1}{3}}
\end{displaymath}

\begin{displaymath}x=-\frac{(x^3 + c)}{b}
\end{displaymath}

Try with:
x3 - 0.01 x + 3 = 0

x3 -100 x - 1 =0

First equation can be solved by first rearrangement but not second. Second equation solves by second method but not first.

Sufficient condition for convergence is:

\begin{displaymath}\vert\frac{d\phi(x)}{dx}\vert \leq 1 \end{displaymath}

at the solution

In general this class of repeated substitution method cannot be guaranteed to work. An exception is a situation where the physics of the system ensures that the above condition is always met. One example of this is in solving process recycle problems.

A Chemical Process with Recycle


Consider the recycle of a reactant.

If c is the fractional conversion in the reactor and s the fractional recovery of the recycled reactant in the separator, then the classic procedure of `guessing a recycle stream' is equivalent to solving for x , the amount of reactant entering the reactor:

x = f + r

Here r=s(1-c)x) so:

x = f = s(1-c) x

Thus$\phi$ is the RHS expression, i.e.:

\begin{displaymath}\phi \equiv s(1-c) x \end{displaymath}

And:

\begin{displaymath}\frac{d\phi(x)}{dx} = c(1-c) \end{displaymath}

Since c and s are both fractions, this must be less than unity.

Repeated substitution methods should be avoided unless it is certain that this condition is fulfilled. Even when it is, other methods are usually better.


Section 3.4: Iterative Solution of Several Algebraic Equations using Single Equation Methods

Section 3.4: Iterative Solution of Several Algebraic Equations using Single Equation Methods

This is intended to be an informal introduction to material which is formally introduced in section 3.5. If you understand what this section is about, then all that the further treatment will really add is some terminology and notation.

The iterative methods which have been described for single nonlinear equations cannot easily be extended to handle problems where it is necessary to iterate on several variables at once, although such methods are available.

However, as we have seen, we can have sets of equations which can be solved directly by rearrangement and reordering without any iteration. It should thus be clear that we would expect to find sets of several equations where iteration is required, but only in a single unknown.

Here are two cases where this is obvious, and a solution method can easily be devised.

Example 1

x + log (x) = 0

y = x²

Clearly an iterative solution in the single variable serves to determine x. y is then calculated from the second equation which is already written appropriately.

We can refer to the second equation as the `tail' of the equation set.

Example 2

y² = 4
x + log (x) + y= 0

The first equation involves only y and so can be solved. once y is known the second equation contains only x, although it must be solved iteratively.

We could call the first equation the `head' of the equation set.

Iterative Solution of `Quasi-single' Equations

The above cases are trivial and obvious. The following is slightly less obvious, although it has already been used as an example.

1 - (p1 + p2) = 0
where p1 = exp(11(1-300/x))
and p2 = exp(11(1-250/x))

This is really a set of 3 equations in 3 unknowns, no one of which, as written, may be solved on its own for any single variable.

However, it is obvious that by simple substitution we can arrive at one equation in x, namely:

1 - ( exp(11(1-300/x)) + exp(11(1-250/x)) ) = 0

Once x is determined iteratively then the second and third equations, which form a `tail', may be used to determine p1 and p2 directly.

When a numerical solution is to be determined there are two ways of setting the problem out. One is essentially as above, but the other avoids any algebraic substitution with the attendant risk of errors. In each case we show the Fortran instructions necessary to evaluate both the l.h.s. of equation to be solved iteratively, i.e.
1 - (p1 + p2)

and the `tail variables' p1 and p2.

Method 1: bad

f = 1.0 - (exp(11.0*(1.0-300.0/x) + exp(11.0*(1.0-250.0/x))
! This evaluates the function
! when it is zero x will be correct, and so will...
p1 = exp(11.0*(1.0-300.0/x))
p2 = exp(11.0*(1.0-250.0/x))

Method 2: good

! At the current x evaluate:
p1 = exp(11.0*(1.0-300.0/x))
p2 = exp(11.0*(1.0-250.0/x))
! and hence ..
f = 1.0 - (p1 + p2)
The overall effect in each case is exactly the same. But this second method avoids more algebra.

We call a group of equations requiring iterative solution in this way a `partition' of the set of equations.


Section 3.5: Iterative Solution of Sets of Equations

Section 3.5: Iterative Solution of Sets of Equations

We have now looked at: Many problems in practice do not fall into the above two categories; indeed there are very few problems which are described by a single equation. There are two other factors which may prevent us from using the two simple techniques. In fact we have already shown how to deal with the first situation. The individual equation or equations involve only a single unknown, all others having been calculated, and so may be solved using an appropriate iterative method as single equations. This thus presents no new problem.

The second situation is the more usual one. There are several possible variants.

  1. There is set of several nonlinear equations which have to be solved simultaneously.
  2. The problem involves a set of nonlinear equations but it can be reduced to iteration in one variable, i.e. not all the equations need to be solved simultaneously.
  3. There is a set of equations which do have to be solved simultaneously, but they are all linear in the unknowns for which they must be solved.
In the first, most general situation, it is most convenient to use on of a range of packaged procedures. We will not discuss these here.

In practice a number of important chemical engineering problems can be reduced to the second case. The incidence matrix analysis introduced earlier enables us to distinguish between the first and second cases, and is dealt with in section 3.5.1, systems reducible to iteration in a single unknown.

In the third case, because the equations are are linear, it is possible to obtain what is in effect an analytical solution. This is covered in section 3.6, systematic methods for sets of linear equations.


Section 3.5.1: Systems Reducible to Iteration in a Single Unknown

Section 3.5.1: Systems Reducible to Iteration in a Single Unknown

It is frequently possible, by analysing the structure of a set of equations, typically by creating an incidence or occurrence table to see which equations appear in which unknowns, to discover that what looks like a large set of simultaneous equations can in fact be solved by a combination of rearrangement, and by iterating in a single variable. Both of these techniques have already been covered individually; we will now see how they can be combined to solve more complicated problems.

Illustrative Example

Consider the 4 equations in unknowns w, x, y, z:

log(w) + 3 = 0

w + 5 x2 - y - 4 = 0

x + log( y) = 0

w + x + y - z0.5 + 10 = 0

Without performing any analysis of the equation, we might suppose that it would be necessary to guess the values of all 4 unknowns and iterate on all of them simultaneously, using techniques which we have not yet discussed.

However, experience with problems involving direct solution where we could reorder and rearrange sets of equations to solve without any iteration whatsoever, tells us that we should look more closely at the equations.

Here is an incidence matrix for the equations as written above in the rows and the variables in order w, x, y, z in the columns.

It is obvious that the first equation contains only the unknown w and so can immediately be solved after rearranging to the formula:

w = exp(-3) = 0.0498

w is now known and so is no longer a variable. We can if we wish remove the w column from the incidence table and look at the remaining 3 variable problem. In particular, w is no longer an unknown in the second equation, which now effectively contains only x and y. However, all the other equations also contain at least both x and y and so it is not possible to go any further with `direct' solution.

In fact, it can be seen that the second and third equations now contain both x and y and only these unknowns. They must be solved simultaneously, and for this purpose could be rewritten with the value of w substituted in as:

0.0498 + 5 x2 - y - 4 = 0

x + log( y) = 0

Although the above equations will have to be solved simultaneously somehow, once we have done this and determined values for x and y it can be seen that in the last equation there now remains only the single unknown, z. This can then be determined by rearrangement:

z = (10 + w + x + y

Summary so Far

We can divide the equations in a set of algebraic equations into three categories: We can use the incidence matrix to help identify these subsets.

Solving the Partition

Consider the two `partition' equations: 0.0498 + 5 x2 - y - 4 = 0

x + log( y) = 0

Suppose we knew x and rearranged the first equation to give y:

y = 0.0498 + 5 x2 - 4

If we then substituted the values of x and y into the LHS of the second equation then it should equal zero.

However, if we only guess the value of x, but nonetheless calculate y and substitute it into the second equation, the latter will evaluate to zero if we have chosen the correct value of x. In effect the two equations if rearranged appropriately and treated in sequence:

y = 0.0498 + 5 x2 - 4

f = x + log( y)

effectively result in the evaluation of a single function of x:

f(x) = x + log[0.0498 + 5 x2 - 4]

This is just a single equation in x and can be solved using any suitable single equation solving method. However the function evaluation is carried out in two steps rather than in one. Note that we do not actually make the substitution shown above. Doing the evaluation in two steps saves us the need to carry out any algebra other than the rearrangement of the first equation.

This method of solving sets of equations is called `tearing' and the single variable x for which we solve is called the `tear variable'. In a more complicated system of equations it may be necessary to tear, i.e. guess and iterate on, more than one tear variable, and so it will not be possible to use the simple solution schemes which have been described for single unknowns. However, many practically useful problems may be reduced to iteration in a single variable.

Inspection of the incidence matrix for the partition will show how many tear variables are needed for a given problem: this is equal to the number of variables (columns) which have nonzero table entries above the diagonal of the matrix. For this example it can be seen that there is only one.


Section 3.5.2: Additional Notes - Reducing the number of equations

Section 3.5.2: Additional Notes - Reducing the number of equations



Consider:

\begin{eqnarray*}f_1(x_1) & = & 0 \\
f_2(x_1,x_2) & = & 0 \\
f_3(x_2,x_3,x_5...
... & 0 \\
f_5(x_1,x_4,x_5) & = & 0 \\
f_6(x_5,x_6) & = & 0 \\
\end{eqnarray*}


We note that:

We can call the subset (f1,f2) the head of the equation set, andf6 the tail. (In general, of course, the equations need not have been numbered and ordered so that these appeared at the top and bottom of the list of equations!)


Partitions

The group (f3,f4,f6) is called a partition of the set of equations and represents a subset that must be solved simultaneously. If we delete x1 andx2, which will be known after solvingf1 andf2, the partition is:

\begin{eqnarray*}f_3(x_3,x_5) & = & 0 \\
f_4(x_3,x_4) & = & 0 \\
f_5(x_4,x_5) & = & 0
\end{eqnarray*}


There are three equations in three unknowns, so our `work factor' is 32 = 9, only a quarter of that for the complete set.

Incidence Matrix


\begin{displaymath}\left[
\begin{array}{c c c c c c}
1 & & & & & \\
1 & 1 & & &...
... \\
1 & & & 1 & 1 & \\
& & & & 1 & 1 \\
\end{array}\right]
\end{displaymath}

The partition is identifiable as the the group of equations with above diagonal elements.


Tearing

The partition alone is 3 equations, f3, f4, f5 in $\{x_3, x_4,x_5\}$

These can be solved as follows:

Rearrange f3 to give x5 , in terms of x3 , i.e:

\begin{displaymath}x_5 = \phi_3(x_3) \end{displaymath}

Similarly f5 to give x4 in terms of x5 :
\begin{displaymath}x_4 = \phi_5(x_5) \end{displaymath}

And finally either: Evaluate f4 ; if $f_4 \approx 0$ terminate; otherwise adjust x3

This procedure is referred to as tearing the set of equations, reducing them here to solution for a single unknown x3 , called the tear variable.

The incidence matrix:

\begin{displaymath}\left[
\begin{array}{c c c}
1 & & 1 \\
1 & 1 & \\
& 1 & 1 \\
\end{array}\right]
\end{displaymath}

has one variable with a supra-diagonal element; this tells us that we can reduce the equations to solving for a single variable.

Solution by Tearing

The procedure can be represented:

\begin{displaymath}\begin{array}{ccccc}
\phi_3 & \rightarrow & x_5 & \rightarrow...
...ow \\
x_3 & \leftarrow & \phi_4 & \leftarrow & x_4
\end{array}\end{displaymath}

Clearly any one of x3 , x4 or x5 could have been chosen as the unknown.

Summarizing ...

Solve essentially by rearranging the equations from the form:

\begin{displaymath}{\bf f}({\bf x}) = {\bf0} \end{displaymath}

into:

\begin{displaymath}\tilde{{\bf x}} = {\bf\phi}({\bf x}) \end{displaymath}

$\tilde{{\bf x}}$ , the output set of the equations contains all the xi from ${\bf x}$ but not necessarily in the same order. For the example above ... Head:

\begin{eqnarray*}x_1 & = & \phi_1 \\
x_2 & = & \phi_2(x_1) \\
\end{eqnarray*}


Partition:

\begin{eqnarray*}x_5 & = & \phi_3(x_2,x_3) \\
x_4 & = & \phi_5(x_1,x_5) \\
f_4(x_3,x_4) & = & 0 \\
\end{eqnarray*}


Tail:

\begin{eqnarray*}x_6 & = & \phi_6(x_5) \\
\end{eqnarray*}


Equation set as a Graph

One and only one edge may leave an equation node. One and only one edge may enter a variable node. Partition cycle shown bold.


Bubble Point Calculation Revisited

Given a liquid mixture of n components with mol fraction composition $x_i, \ i=1 \ldots n$ , at pressure P , determine the temperature T , and the composition yi of the vapour in equilibrium with the liquid.
P*i - P*i(T) = 0 (1)
$\displaystyle k_i - \frac{P^*_i}{P}$ = 0   (2)
yi - ki xi = 0   (3)
$\displaystyle \Sigma y_i - 1$ = 0   (4)

The range of i and the summation is from 1 to n , the number of components.

With the equations in the above sequence and the output set, say for a binary, chosen to be:

\begin{displaymath}\{P^*_1, k_1, y_1, P^*_2, k_2, y_2, T \}\end{displaymath}

The incidence matrix is:

\begin{displaymath}\left[
\begin{array}{c c c c c c c}
1 & & & & & & 1 \\
1 & 1...
...1 \\
& & & & 1 & 1 \\
& & 1 & & & 1 \\
\end{array}\right]
\end{displaymath}

The final, summation, equation does not involve T explicitly. However, for any value of T the value of the function can be determined, and will be zero when T satisfies the equations.


Section 3.5.4: Solutions to Example Questions

Section 3.5.4: Solutions to Example Questions

1. The incidence table for the equation set is:

Equationx1x2 x3x4x5
(1)111
(2)11
(3)111
(4)11
(5)1

Rearranging and reordering the incidence table leaves us with:
Equationx5x1x3x4x2
(5)1
(4)11
(3)11
(1)111
(2)11
For this set of equations, equation (5) is the head, equations (4), (3) and (1) the partition, and equation (2) the tail. Either of x1, x3 or x4 would be suitable for tear variables here.


Section 3.6: Systematic Methods for Sets of Linear Equations

Section 3.6: Systematic Methods for Sets of Linear Equations

A set of equations which is linear in all its unknowns has the special property, like a single linear equation, that it has what is effectively an analytical solution.

It is not usual to attempt to write out such a solution, for a system of more than two equations it becomes very cumbersome, but the systematic procedure which enables an analytical solution to be found can more conveniently be used to compute the numerical values of the unknowns.

Definitions

These properties formally exist only for a set of equations which is linear in all its unknowns. A set of equations is linear if and only if all its unknowns appear only as first powers.

In general, if any unknown appears to a power other than unity, or as the argument of a function such as log or exp, then the whole set of equations becomes nonlinear and in general cannot be solved in this way. However, we will identify some special cases where a subset of the equations, in fact a partition (section 3.5.1) of the set, is entirely linear then it will be possible to solve a set of equations with a `head' and `tail' which are nonlinear.

We will describe a set of equations with an independently solvable linear subset as partially linear. These will be dealt with in connection with a process case study (section 3.7.2).

Example

If x, y, z are unknowns, then the following comprise a set of 3 linear equations:

x - y + 2 z = 2

2 x + 2 y = 4

x + 3 y + z = 0

Each equation as written is of the form:

(coeff. of 1st variable) (1st variable) +/- (coeff. of 2nd variable) (2ndvariable) +/- ...... = constant

This is the only form possible if the equation is linear in all its variables. (Although the constant could have been written on the LHS with zero on the RHS as has been the usual practice with previous equations, it is conventional to write linear equations in this way.)

The regular arrangement of the equations and variables lends itself to a tabular representation as shown below.
x y z   Constant
Equation 1 1 -1 2   2
Equation 2 2 2 0 = 4
Equation 3 1 3 1   0

We can think of the set of equations as being represented by:

If we let a single symbol stand for each of these entities, e.g.

Then the set of equations can be written as symbols in matrix-vector notation:

A u = b

Written in full this becomes:

Or more usually with square brackets throughout (see notes on vectors (section 1.1.1) ) as:

The individual matrix coefficients may be referred to if required, ( notes on matrices (section 1.2.1.1) as:

ar,c

where r is the row number and c is the column number.

Solution of Linear Equations

The solution of sets of linear equations is a standard problem and many computer programs are available to perform it. It is extremely unlikely that you should ever have to write such a program, and so we will not spend time on describing the methods in detail.

However, a summary of the method is available in section 3.6.1 along with an algorithm in the form of a Fortran program.

To solve a set of linear equations with specified numerical coefficients and constant values we have provided access to such a program in section L1.2

One of the things that spreadsheet systems cannot conveniently do is solve sets of linear equations. However, we have developed ways of doing this. If you would like to download a linear solver for up to 12 equations which can be loaded into an Excel spreadsheet, then save this link in your local memory and load it. The file is in Microsoft Sylk ASCII format which is acceptable to a number of spreadsheets.

There are a number of more sophisticated options available for making spreadsheet based linear solvers which you can find out about in section L2.4 .

Failure of Linear Solution Procedures

Two things may go wrong in the course of attempts to solve a set of linear equations.

The set of equations may not have a solution. For example, consider the following, all of which can be written as sets of equations, and indeed have coefficients and constants entered into one of the solvers, but none of which have valid solutions. You can check this by trying them with a solver.

x = 3
2x = 4
The two equations are are inconsistent.

x1 = 3
2x1 = 6
The equations are consistent, but if intended for solution in two unknowns, i.e. x1 and x2, the second equation provides no information additional to the first, and so there is no means of solving for x2.

x1 + x2 = ..
x1 + x1 = ..
Depending on the values on the RHS the two equations are either inconsistent or one is redundant. In neither case can a solution be found.

x1 + x2 = ..
2x1 + 2x1 = ..
This is the same situation, although slightly less obvious.

In all the above situations the equations are said to be linearly dependent. This occurs whenever one LHS is a multiple of that of another equation. More generally this occurs if any LHS can be constructed from any linear combination of another, i.e. by adding together multiples of right hand sides. The solution procedure will fail for such a set of equations regardless of the values on the RHS.

The problem is not always obvious:
x1 + x2 + x3 = 1
3 x1 -2 x2 + x3 = 2
5 x1 + 3 x3 = 2
The 3rd LHS is equal to twice the first plus the second.

The above are all examples of equations sets which are said to be singular . This depends only on the structure or numerical values of the coefficient matrix, which is also said to be singular.

Another problem can occur during the solution procedure (section 3.6.1), because this involves adding and subtracting coefficients, if the coefficients are are very different sizes. Because computers work to a limited numerical precision, if a very small number is added to a very large one. This can also happen if the values of the solutions are very different, e.g. one is a small mole fraction (say 0.0001) and another is an enthalpy (say 100,000).

In these circumstances it is a good idea to scale the problem variables to get them within the same 2 or 3 orders of magnitude.

Numerical accuracy also becomes significant if the equations are nearly singular.


Section 3.6.1: Elimination Methods for Linear Equations

Section 3.6.1 Elimination Methods for Linear Equations

Solve:

4 x1 + 2 x2 =  1

x1 - 4 x2 = 5

  1. Eliminate x1 from both equations by making the coefficients of x1in each the same and subtracting:
  2. Now substitute x2 back into the first equation as reduced in the first step of 1. above, i.e.:

    x1 + 0.5 (-1.055556) = 0.25

    This is in a convenient form for obtaining x1:

    x1 = 0.25 + 0.5 (1.055556) = 0.77778

In symbols:

\begin{eqnarray*}a_{1,1} x_1 + a_{1,2} x_2 & = & b_1 \\
a_{2,1} x_1 + a_{2,2} x_2 & = & b_2 \\
\end{eqnarray*}

 

Eliminate x1

( a2,2 - a2,1 a1,2 / a1,1 ) x2 = b2 - a2,1 b1 / a1,1


\begin{displaymath}x_2 = \frac{b_2 - a_{2,1} b_1 / a_{1,1}}
{ a_{2,2} - a_{2,1} a_{1,2} / a_{1,1} } \end{displaymath}


x1 = b1 / a1,1 - (a1,2 / a1,1) x2

Substituting:


\begin{displaymath}x_1 = b_1 / a_{1,1} - a_{1,2} / a_{1,1}
\frac{b_2 - a_{2,1} b_1 / a_{1,1}}
{ a_{2,2} - a_{2,1} a_{1,2} / a_{1,1} } \end{displaymath}

Set of n equations

\begin{eqnarray*}a_{1,1} x_1 + a_{1,2} x_2 + %
a_{2,1} x_1 + a_{2,2} x_2 + %
\vdots \ \ \ \ & = & \vdots \\
a_{n,1} x_1 + a_{n,2} x_2 + %
\end{eqnarray*}

\begin{displaymath}\left[
\begin{array}{c c c c}
a_{1,1} & a_{1,2} & & a_{1,n} \...
...egin{array}{c}
b_1 \\ b_2 \\ \vdots \\ b_n
\end{array}\right]
\end{displaymath} (1)

This may be written in matrix-vector form:

\begin{displaymath}{\bf A x = b}
\end{displaymath}
     !Start the elimination loop.....  
       do i=1,n-1  
     !... go down the matrix by rows i...  
     !.. reducing all the other rows k.. 
       do k=i+1,n 
       pivot=(a(k,i)/a(i,i))  
     !This is the number to multiply each  
     !row of equn i by before subtracting it..  
     !.. along each row k ..  
       do j=1,n
       a(k,j)=a(k,j)-pivot*a(i,j)  
       end do  
     ! .. and also the constant.  
       b(k)=b(k)-pivot*b(i)  
     !  
       end do  
     !.. go on to next row with current elimination 
     !end do  
     !.. and then eliminate next equation  
     !  
     !.. and finally solve the one remaining  
     !equation for x(n) 
       x(n)=b(n)/a(n,n) 
     !--- Elimination Complete --------- 
     !Back substitute for other x() ...
       do l=n-1,1,-1 
       x(l)=b(l) 
       do j=l+1,n  
       x(l)=x(l)-a(l,j)*x(j) end do  
       x(l)=x(l)/a(l,l)  
       end do  
     !----- Finished ------------ 
    

Using the Library Solver

       program lintest  
       implicit none 
     !test linear solver  
       integer :: n,i,j,ifg  
       integer, parameter :: max=50  
       real :: a(max,max),x(max),c(max) 
       print *, 'Solver for small sets of linear equations...'
       print *, 'Enter number of equations, less than ',max  
       read *, n  
       print *, 'Now enter coefficients ...'  
       print *, 'i.e. ',n+1, ' numbers per row...'  
       do i=1,n  
       print *, 'Elements for row ', i, ':'  
       read *, (a(i,j), j=1,n), c(i) 
       end do  
       ....... 
       print *, 'Matrix and constants are:'  
       do i=1,n  
     !check print  
       write(6,'(7f10.4)') (a(i,j),j=1,n), c(i)  
       end do 
       call solve_linear (a,max,x,n,c,ifg) 
       write(6,*) 'flag =',ifg  
       write(6,*) 'solution..'  
       write(6,*) (x(i),i=1,n) 
       end program 


Section 3.6.2: Example Questions

Section 3.6.2: Example Questions

To solve these examples you can either use the web-based linear solver (Section L1.2) or you can make a spreadsheet solver (Section L2.4) for yourself.

If you are using Fortran:

You should copy the linear solver demonstration program and use it for the first seven examples. For the final example, modify the program so that it sets up the equation coefficients by means of assignment instructions rather than asking for them from the keyboard.

Problems

In the following examples, all symbols represent unknowns unless otherwise stated. Some problems do not lead to solutions. The web solver will in this case return all unknowns as zero and and print an error message. The spreadsheet solver will produce error messages in the solution cells. if this happens, look carefully at the equations you have suppled and see why a solution could not be found.


Section 3.6.3: Solutions to Questions

Section 3.6.3: Solutions to Questions

This section contains the answers to the linear equations question in section 3.6.2.

1. x1 = 2, x2 = 1

2. x1 = 0.4, x2 = 2

3. x1 = -0.5, x2 = 5

4. x1 = 2, x2 = 1

5. x1 = -0.15, x2 = 310.6

6. x1 = 5, x2 = -3, x3 = 2

7. x1 = -3.27, x2 = 14.25, x3 = 17.52, x4 = -18.50, x5 = -1.95


Section 3.7: Case Study

Section 3.7: Case Study

There are two case studies here. The first involves vapour liquid calculations. The theory is discussed and several example problems are given with solutions. This section finishes off with some flash calculations which will be completed as a handin.

The second case study looks at a complete process with recycle. This will also be completed as a handin.


Section 3.7.1: Phase Equilibrium Calculations

Section 3.7.1: Phase Equilibrium Calculations

Before looking at any of the following sections, please read over the revision material in section 1.2.2, Binary Vapour Liquid Equilibrium.


Section 3.7.1.1: Note on Relative Volatility

Section 3.7.1.1: Note on Relative Volatility

A convenient measure for use particularly in approximate vapour-liquid equilibrium calculations is the relative volatility of one species with respect to an arbitrarily chosen reference or key component.

Relative volatility $\alpha_{i,r}$ of species i to a reference component r is defined as:


\begin{displaymath}\alpha_{i,r} = \frac{K_i(T,P)}{ K_r(T,P)}
\end{displaymath}

For an ideal system:

\begin{displaymath}K_i = \frac{P^*_i(T)}{ P} \end{displaymath}

Hence:

\begin{displaymath}\alpha_{i,r} = \frac{P^*_i(T)}{ P^*_r(T)}
\end{displaymath}

Although relative volatility is thus a function of temperature T, much of its usefulness results from it being a rather weak function of temperature, particularly for nearly ideal mixtures. This is because the slope of vapour pressure versus curves tend to be similar for similar components.

Thus over a reasonably limited range of temperatures relative volatility may be estimated at a suitable `average' or reference temperature and then regarded as a constant.

For an n component mixture there are of course n relative volatilities. However if the reference component is taken as on e of those in the mixture, which would be the usual practice, then the volatility of this component w.r.t. itself will be one.


Section 3.7.1.2: Example Questions

Section 3.7.1.2: Example Questions

  1. It is required to determine the pressure at which a liquid mixture of propane and butane will be at its bubble point at a temperature of 293K. Prepare an incidence table for this system of equations and show that the problem can be reduced to iteration on a single variable. (Most of this problem has already been covered in the course material.)

    Calculate the bubble point pressure for a an equimolar mixture of the species and the composition of the vapour phase formed.

  2. Repeat the above analysis and calculation for the dewpoint pressure calculation at 293K.

  3. Repeat the analysis of question 1 for the case of calculating the bubble point temperature given the pressure. Show that the problem can be reduced to iteration on a single variable. (Most of this has also been described in the course material.) Calculate the bubble point temperature of an equimolar butane-pentane mixture at 2 bara.

  4. Water and pentane form two immiscible liquid phases. In these circumstances it can be assumed that each liquid phase contains only a single pure component and the equilibrium condition is given by:

    (P*1(T) + P*2(T)) - P = 0

    Set up the equations, incidence table, etc. for this problem and show how it can be solved by iteration in a single variable. Show also that the composition of the vapour can be calculated. Note that this is technically a minimum boiling heterogeneous azeotrope. Determine the azeotrope temperature and composition at 2 atm.

  5. The vapour pressures of butane and pentane at 300K are 2.73 atm and 0.713 atm respectively. Show that if the pressure is specified the problem can be described by 4 linear equations in the 4 mole fractions (two liquid and two vapour). find the two phase compositions at 2 atm and 300K by setting up and solving the 4 equations. (Here is a link to the linear equation solver (Section L1.2)).
N.B.: Create your spreadsheet models using the spreadsheet model generator (Section L2.5).

Data

Boiling points:
Section 3.7.1.3: Solutions to Questions

Section 3.7.1.3: Solutions to Questions

1.
If we let propane be component 1, and let butane be component 2 then the equations for the system are as follows:

(1) P1* = exp((11(1 - (TB1/T)))
(2) P2* = exp((11(1 - (TB2/T)))
(3) Py1 = x1P1*
(4) Py2 = x2P2*
(5) y1 + y2 = 1

where:
Pi* = vapour pressure of component i
TBi = atmospheric boiling point of component i
T = temperature
yi = vapour mole fraction of component i
xi = liquid mole fraction of component i
P = pressure

The unknowns here are: P1*, P2* y1, y2 and P. There are five equations and five unknowns, therefore in principle the system may be solved.

TB1, TB2 and T have been given. We are told that there is an equilmolar mixture of the species. As we are calculating the bubble point, this means that there is only liquid present. Thus x1 = x2 = 0.5

The incidence table for the above set of equations is:
EquationP1* P2*y1 Py2
(1)1
(2)1
(3)111
(4)11
(5)11

Looking at this incidence table, it should be obvious that it will not be possible to rearrange this into a lower triangular form. There is a head, equations (1) and (2), and a partition, equations (3), (4) and (5). The spreadsheet solver may only solve equations arranged into a lower triangular format, thus some more work must be done.

There are two methods for solving this. The first involves rewriting equations (3) and (4) using partial pressures instead of vapour mole fractions and pressure. This would leave the set of equations as used in the revision section 1.1.2.2. This would mean extra manipulation to answer the question. This approach will be left as an exercise for the reader.

The other approach is to define a tear variable, and to give a value to one of the unknowns. The "solver" function in the spreadsheet may then be used to obtain the result. This approach will be used here.

The model for the spreadsheet may be viewed here, and is used to generate this spreadsheet.

model question1

parameter
 Tb1 = 231    !1 atm boiling point of propane
 Tb2 = 272.6  !1 atm boiling point of butane
 x1 = 0.5     !liquid mole fraction of propane
 x2 = 0.5     !liquid mole fraction of butane
 T = 293      !temperature of system
 P = 1        !pressure of system
end parameter

variable
 P1star  !vapour pressure of propane 
 P2star  !vapour pressure of butane
 y1      !vapour mole fraction of propane
 y2      !vapour mole fraction of butane
 tear    !tear varible
end variable

equation
 P1star = exp(11*(1-(Tb1/T)))
 P2star = exp(11*(1-(Tb2/T)))
 y1 = x1*P1star/P
 y2 = x2*P2star/P
 tear = y1 + y2 - 1
end equation

steadystate

end model question1

In order to have solved the set of equations, the tear variable must be equal to zero. This is achieved by adjusting the pressure. It may either be altered manually, or the solver function of the spreadsheet (if present) may be used. On Excel this may be found under the Tools tab.

Using this gives the following results:
Pressure = 6.2 bar
y1 = 0.827
y2 = 0.173


Section 3.7.1.4: Flash Calculations with Fixed Relative Volatility

Section 3.7.1.4: Flash Calculations with Fixed Relative Volatility

If the phase equilibrium properties of a mixture can be represented adequately by fixed relative volatilities, see notes here, then the equations describing a vapour-liquid `flash' system become particularly convenient to handle.

The situation is as illustrated in the diagram.



Section 3.7.2: Recycle Process Case Study

Section 3.7.2: Process Case Study

You should first look at some revision material on material balances in section 1.2.1.

We describe here a simple recycle chemical process which will be used in a case study to illustrate a number of aspects of steady-state modelling and algebraic equation solving.

Process Description

The process is to make product C from reactants A and B by the reaction

A + B -> C

It is fed with a mixture of A and excess B containing none of the product material C.

The conversion in the reactor is R based on A.

Following the reactor all unreacted A and B are recovered in the overheads of a distillation column and recycled. The product C is not fully separated, a fraction S of that fed to the separator remains in the overheads.

To avoid the buildup of unreacted B, a fraction P of the overhead stream is removed as a purge.

Simple Process Model

The parameters of the model are: 18 variables x1 to x18 can be associated with the 3 molar flows of each of the species A, B and C in the 6 unknown process streams as shown in the flowsheet above. Here x1 is the flowrate of A entering the reactor, x2 that of B, and x3 that of C. Similarly x4 is the flowrate of A out of the reactor, and so forth.

The operations are all described by linear equations, as discussed in section 1.2.1. hence the whole process amy be described by a set of 18 linear equations as shown below.


Section 3.7.2.1: Process Case Study Assignments

Section 3.7.2.1: Process Case Study Assignments

Preliminary Exercise

Set up the material balance equations in matrix-vector form and write them out.

Using the `typical' values given for the parameters, solve using one of the solvers with which you were supplied.

Linear Process Modelling Exercise

The above method of describing a process is not very convenient, and the output obtained from the material balance calculation is not related to the way in which an engineer would think about the plant, i.e. in terms of streams and component flows.

You are required to produce a program (or spreadsheet if you are not working with Fortran) in which you need specify only the three process parameters of conversion, recovery and purge fraction. the program should be designed so that these can easily be set by the user.

This program should also produce a process-oriented output rather than a list of values. It is suggested that this should be a table of streams and component flowrates, e.g. with a column for each component and a row for each of the six streams (or vice-versa), all suitably labelled.

Partially Linear Process Model

1. It is required to present the results in terms of total component flows and mol fraction compositions. The latter are of course calculated from the ratio of one component flow to the total flow. This involves a nonlinear equation; however, this nonlinear subset constitutes a `tail' of the equation set and can also be calculated directly, i.e. without iteration.

Adapt your model to calculate this information.

2. To maintain reactor performance it is necessary to ensure that no more than 1% of C is present in the reactor inlet stream. This is done by adjusting the purge to maintain that value.

The set of equations now becomes nonlinear because P is now an unknown equations like the following occur:
x13 - P x10 = 0

Because there is an additional unknown another equation is needed. This is provided by the mol fraction specification at the reactor inlet:
x1 / (x1 + x2 + x3) - 0.01 = 0

(As written this is also nonlinear, although unlike the preceding equation, it may be rearranged to linear form if required.)

However, if we regard P as a tear variable we may guess a value for it and see if the mol fraction specification equations is satisfied. This is in effect iteration on a single variable and so may be handled as was, e.g. the bubble point calculation. The whole linear equation solution procedure is `nested' within this iteration.

Adapt your model as necessary to find the required purge rate at 0.8 fractional conversion and 0.1 fractional recovery of C in the column overheads. If you are using bisection you will appreciate that P must lie between 0.0 and 1.0. However, a value of zero will cause a `divide by zero' in the linear solver, so use a small positive value instead.


Section 3.7.2.2: PC Linear Solver

Section 3.7.2.2: PC Linear Equation Solver

The following linear solver can be run on a PC from an MS-DOS window. The command to run it is just:
linprog

You need to copy the executable , the memory manager and a library file. These can most conveniently all be in the same directory.