ON THE APPLICATION OF A NEWTON

RAPHSON’S ITERATIVE METHOD OF THE

FIXED POINT THEORY TO THE SOLUTION OF
A CHEMICAL EQUILIBRIUM PROBLEM
Okereke, C. Emmanuel, Ph.D. †

ABSTRACT

In this work, I discussed the solution of a chemical equilibrium

problem aiming to obtain its fixed point. To do this, the preliminary and basic
ideas introducing the fixed point theory was x-rayed and the Newton Raphson’s
iterative method for solving the system of non-linear equations discussed, then
the problem of the chemical equilibrium involving principal reactions in the
production of synthesis gas by partial oxidation of methane with oxygen were
stated.
Using a computer program the O reactant ratio that produce an
adiabatic equilibrium temperature were obtained by developing a system of
seven simultaneous nonlinear equations that has the form which we now solve
using the Newton Raphson’s method described in section 2.2 and hence the
desired fixed point of the chemical equilibrium problem.

THE NETWON’S METHOD, A PRELLIMINARY TO THE NEWTON

RAPHSON’S ITERATIVE METHOD
Let T be an operator mapping a set X into itself, a point x ∈ X is called
a fixed point of T if
x = T (x) … (1.1)
By (1.1), we achieve a natural construction of the method of successive
approximations.
x n+1 = T (xn), n ≥ 0 ∈ x … (1.2)
and if the sequence (xn), n ≥ 0 converges to some point x = x* ∈ X for some
initial guess x0 ε x, where T is a continuous operator in a Banach space X, we
have
x* = lim x
n → ∞ n +1 n→ ∞ n (n → ∞ n )
= lim T ( x ) = T lim x
that is x* is a fixed point of operations T. Hence we now state without proof the
following important results that make easy the understanding of the Newton
Raphson’s method used in this work
Theorem (1.1): If T is a continuous operator in a Banach X, {xn}, (n ≥ 0)
generated by (1.2) converges to some point x* ∈X for some initial guess x0 ∈ X
and we say that x* is a fixed point of the operator T [1].
To investigate the uniqueness property, we introduce the concept of
contraction mapping as follows. Let (x,|| ||) be a metric space and T a mapping of
X into itself. The operator T is said to be a contraction if there exists a real
number k, 0 ≤ k < I such that

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 1

 F ( x) − F ( y ) ≤ k x − y , for all x, y ε x … (1.3)
Hence, every contraction mapping T is uniformly continuous. Indeed T is
lipschitz continuous with a lipschitz constant k which may also be called the
contraction constant for T. With the above; we now discuss the Banach fixed
point extensively as related to the target of this research.
Theorem 1.1A [9] (Banach fixed point theorem (1922). Suppose that
we are given an operator T:M ⊆ X → M, i.e., M is mapped into itself by T;
M is a closed nonempty set in a complete metric space (x,d);
T is k – contractive, i.e. d(Tx, Ty) ≤ k .
Then the following hold:
Existence and uniqueness:- T has exactly one fixed point on M;
Convergence of the iteration: the sequence {xn} of successive approximations
converges to the solution, x. for an arbitrary choice of initial point x0 in M;
Error estimate: for all n = 0, 1, 2… we have the a prior error estimate d(xn, x) ≤
kn (1-k)-1 d(x0, x), and the a posteriori error estimate d(xn+1, x) ≤ k(1-k)-1
d(xn,xn+1);
Rate of convergence; for all n =0,1,2,…We have d(xn+1,x) ≤ kd(xn, x)
Definition 1.1 [9] An operator T:M ⊆ X → X on a metric space (X,d) is called
k- contractive if f (1.3) holds for all x, y ε M with fixed k, 0 ≤ k < 1,
T is called Lipschitz continuous and if.
d(Tx,Ty) < d(x,y) for all x, y ε M with x ≠ y, … (1.4)
T is called contractive for T and we obviously have the implications:-
k- Contractive → contractive→ no n expansive → Lipschitz continuous.
Every Banach space called the (x,|| ||) also is a complete metric space as
(x,d) under d(x,y) = ||x-y||.
On a B-space, (1.3) therefore becomes
||Tx-Ty|| ≤ k|| x-y||.
On the above the following follows
{xn} is a Cauchy sequence. This follows from
d(xn,xn+1) = d(Txn-1,Txn) ≤ kd(xn-1.xn)
≤ k2d(xn-2,xn-1) ≤… ≤ knd(x0,x1).
(1.5)
Repeated application of the triangle inequality and finally summing the formula
for a geometric series yields
d(xn,xn+m) ≤ d(xn,xn+1) + d(xn+1,xn+2)+…+ d(xn+m-1, xn+m)
≤ (kn +kn+1 +…+kn+m-1)d(x0,x1)
≤ kn(1-k)-1d(x0,x1)…
(1.6)
Since X is complete, the Cauchy sequence converges, i.e., xn→x as n → ∞ [3].
Equation (1.5) follows by letting m → ∞.
(II) The error estimate (1.6) follows by letting m →∞ in
d(xn+1,xn+m+1) ≤ d(xn+1,xn+2) +…+ d(xn+m-1, xn+m+1)
≤ (k +k2 +…+km)d(xn,xn+1)

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 2

 ≤ k(1-k)-1d(xn,xn+m)…
(1.7)
The point x is a solution of (1.1) for T is continuous by (1.4). Since T(M) ⊆ M
and x0 ∈ M, we have xn ∈ M also, for all n. since M is closed and xn →x as n
→ ∞, we get x ∈ M. Equation (1.2) implies that Tx = x for n → ∞.
Equation (1.6) follows d (xn+1, x) = d(Txn,Tx) ≤ kd(xn, x).
Uniqueness of solution. Suppose x = Tx and y = Ty; then d(x,y)
= d(Tx, Ty) ≤ kd(x,y), which forces d(x,y) = 0, i.e. x =y.

Continuous Dependence On A Parameter

It is important to note that in many applications, T depends on an additional
parameter P. then, (1.1) is replaced by the equation.
xp = Tpxp, xp∈M, … (1.8)
where p ∈ P.
Proposition 1.2 (Corollary to Theorem 1.1A.) let the following conditions be
P is a metric space, called the parameter space.
For each P, the operator Tp satisfies the hypotheses of Theorem (1.A) but with k
in (1.3) independent of p.
For a fixed p0 ∈P, and for all x ∈ M, limp→p0 TpX=Tp0X.
Then, for each p∈P, (1.8) has exactly one solution
xp∈M, and limp→po xp = xp0 [3].
Proof: later xp be the solution of (1.8) given by theorem 1.1 A, then
d(xp, xp0) = d(Tpxp, Tp0xp0)
≤ d(Tpxp,Tpxpo) + d(Tpxpo, Tpoxpo)
≤ kd(xp,xpo) + d(Tpxp0, Tp0xp0),
And therefore,
d(xp, xp0) < (1-k)-1d(Tpxpo,Tpoxpo) → 0 as p →po, by (iii). [9]

1.5. Accelerated Convergence and Newton’s method [5]

We begin with the insight which underlines the acceleration of iterative
methods. Let x be a solution of the real equation X = F (x), and suppose the
sequence of iterations (x, where
xn+1 = f (xn) …
(1.9)
And xn €[a, b] for all n, converges to x as n → ∞.
Now for the key: Suppose further that f is m – times differentiable on [a,b], with
f′(x) = f(2)(x) =… = f(m-1) = 0
(1.10)

(a)
x1 Figure 1.4 (b) x1
X0 x0

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 3

Since xn+1 = f(xn) and x = f(x), we have
|xn+1-x| ≤ sup f(m) (ξ)|| xn – x/m/m! - - -
a <ξ <b
(1.11)
If the supremum in (1.11) is finite, we obtain the convergence of order m, as
opposed to the linear convergence (m=1) of (1.9)
Example 1.1[10]. The trick to Newton’s method consists of rewriting the
f ( x)
equation f(x) = 0 in the equivalent form. X = F(x), where f(x) = x -
f ( x)
f ( xn )
Then the iterative method becomes xn+1 = xn -
f ( xn )
We assume that f′ (xn) ≠ 0 for all n. the, f′(x) = f(x) f n (x) / f′(x)2
So that if x is a solution of f (x) = 0 with f′(x) ≠ 0, then f′(x) = 0. Thus we have a
method with m = 2 in (1.1.9) i.e. we have quadratic convergence.
We apply this to the equation x = T(x) in (1.1.1) i.e. f(x) = T(x).
Computing, we obtain the iterative values xn with linear convergences.
The geometric interpretation of Newton’s methods is seen in figure 1.4
(a). To find a zero, x, of f, take the initial value, x0, and determine the
corresponding functional value, f(x0). The next iterative value, x1 is the
intersection of the tangent line at (x0,f(x0) and the x- axis. Keep repeating the
process, it is typical of Newton’s method that it converges very rapidly if the
initial value x0 is already in the vicinity of the zero, but figure 1.4(b) shows a
better of this.
Q0
D0

Figure 1.5
t
t0 -c t0 t0 +c
However, we know that the above discussed fixed point method is just the
traditional fixed point method that is restricticted to the solution of only linear
systems and for the purpose of this research we advance onto the modified
Newton’s method which is the Newton Raphson’s iterative method here below
generated for use in section three

2. NEWTON’S RAPHSON’S METHOD

Sections 2 and 3 are concerned with finding the solution, or solutions, of the
system.

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 4

 f1(x1, x2 ,...,xn ) = 0,
f 2 ( x1, x2 ,...,xn ) = 0, M (2.1)

f n ( x1, x2 ,...,xn ) = 0,

Involving n real functions of the n real variables x1,x2,…,xn. Following the

previous notation, x = {x1,x2,…,xn}t, we shall write fi (x) = fi (x1,x2,…,xn) here,
and in the subsequent development, 1 ≤ i ≤n. Then let α = [α1,α2,…, αn]t be a
solution of (2.1), that is, let fi(α) = 0.
Let the n functions fi (x) be such that xi = Fi (x) (2.2)
Implies fi (x) = 0, 1 ≤ j ≤ n. Basically, the n equations (2.2) will constitute a
suitable rearrangement of the original system (2.1). In particular,
let αi = fi (α) … (2.3)
Let the starting vector x0 = [x10, x20,…,xn0]t be an approximation to α. Define
successive new estimates of the solution vector, xk = [x1k,x2k…,xnk]t, k = 1, 2,…,
by computing the individual elements from the recursion relations.
xik = Fi (x1,k-1, x2,k-1,…,xn,k-1). (2.4)
suppose there is starting R describable as |xj - αj| ≤ h, 1 ≤ j ≤ n, and for x in R
there is a positive number µ, less than one, such that
n
∂Fi ( x)
Σ
j =1 ∂x j
≤ µ ( 2 .5 )

Then, if the starting vector x0 lies in R, we show that the iterative method
expressed by (2.4)
converges to a solution of the system (2.1), that is,
lim x k = α ( 2 .6 )
k → ∞

Using the mean-value theorem, the truth of (2.1) is established by first noting
from (2.3) and (2.4), that
xik- αI = Fi (xk-1)-Fi (α).

n ∂ F i [α + ξ i , k − 1 ( x k − 1 − α )]
= Σ (x j ,k −1 −α j ) , ( 2 .7 )
j =1 ∂x j

In which 0 < ξi,k-1 < 1. that is,

n ∂Fi
xik − α ≤ h Σ ≤ µh < h, (2.8)
j =1 ∂x j

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 5

Showing that the points xk lie in R. also, by induction, from (2.5) and (2.7),
|xik -α| ≤ µ max (|xj,k-1- αj|) ≤ µkh. (2.9)
Therefore, (2.6) is true, and the procedure converges to a solution of
(2.1). Note that if the Fi (x) are linear, we have the Newton’s method, and the
sufficient conditions of (2.5) are the same as the second set of sufficient
conditions controlling the Newton’s iteration.
For the nonlinear equations, there is also a counterpart to the Newton’s
method, previously discussed for the linear case. We proceed as before, except
that some replacements are made by
Xik=Fi (xik,x2k,…,xi-1,k,xi,k-1,…xn,k-1). (2.10)
That is, the most recently computed elements of the solution vector are always
used in evaluating the Fi. The proof of convergence according to (2.10) is much
the same as for the jacobi-type iteration. We have
n ∂Fi (Σ )
x ik − α i = Σ ( x j , k −1 − α j ) ik
,
j =1 ∂x j
Where ∑Xik - α = [α1 + ξik (x1,k-1 - α),…,αn + ξik (xn,k -1- αn]t,
It will appear inductively that the above is true, because the various points
concerned remain in R. If ek-1 is the largest of the numbers |xj,k-1 - αj|, then
xik- α1| ≤ µek-1 <ek-1 < h.
it follows that
∂F (Σ ) n ∂F (Σ )
− α 2 = (x -α1) 2 2k + Σ  x − α j  2 2k
x 
2k ik ∂x j j = 2  j ,k −1  ∂ xj

where
ε2k = [α1+ ξ2k(x1k-α1), α2+ξ2k(x2,k-1-α2),… αn + ξ2k (xn,k-1- αn)]t. that is, |x2k-α2| ≤
µek-1 <ek-1 < h.
Therefore, |xik- αi| ≤ µkh, and convergence according to (2.1) is again
established.
Observe that the first of the sufficiency conditions of the same (2.10) has been
reaffirmed under slightly general circumstance.

2.2 Newton- Raphson’s Iteration for Nonlinear Equations.

The equations to be solved are again those of (2.1), and we retain the
nomenclature of the previous section. The Newton-Raphson process, to be
described, is once more iterative in character. We first define.

∂ f i (x )
f ij ( x ) = ... ( 2 . 11 )
∂x j
Next define the matrix φ (x) as
φ (x) = (fi(x)), 1 ≤ I ≤n, 1 ≤ j ≤ n. … (2.12)

Thus det (φ (x)) is the Jacobian of the system (2.1) for the vector x = [x1,
x2,…,xn)t. now define the vector
f (x) as f(x) = [fi (x), f2(x),…, fn(x)]t. (2.13)

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 6

With these definitions in mind, and with the starting vector
x0 = [x10,x20,…,xn0]t, let xk+1 = xk+δk, … (2.14)
The fundamental theorem concerning convergence is much less restrictive than
those of the previous sections. We have the result that if the components of φ (x)
are continuous in a neighborhood of a point α such that
f (α) = 0, if det (φ(α)) ≠ 0, and if x0 is “near” α, then lim k →∞ x k = α .
… (2.15)
An outline for a method of proof follows. By (2.13) and (2.14), since
fi (α) = 0, δk = φ-1 (xk) [f(α)]. …
(2.15)
By the mean- value theorem,
n
f i (x k ) − f i (α ) = Σ ∫ ij (α + ξ ik ( x k − α ))( x jk − α ),
j =1

where 0 < ξ ik < 1.

For the ith row of a matrix ψ use [f11 (α + ξik (xk - α)),…, fin(α + ξik(xk- α))].
Then xk+1-α = xk - α + δk = φ-1 (xk) φ(xk)-ψ] (xk- α).
Since the entries in the matrix φ (xk)- ψ are differences of the type fij (xk) - fij
(α+ξik (xk- α)), they can be kept uniformly small if the starting vector X0 lies in
an initially chosen region R describable as |xi-αi| ≤ h, 1 ≤ i ≤ n. concurrent with
this is the fact that since det (φ(xk)) can be bounded from zero. The net result is
that, for 0 < µ < 1, |xik-αi| ≤ h µk,1, ≤ i ≤ n. thus the sequence [xk] converges to α
[6].
Example (2.1) [7]. To illustrate the procedure, we use the example below
namely.
1 x x
f1 (x1 , x2 ) = sin(x1 , x2 ) − 2 − 1 = 0
2 4π 2
 1  ex
f 2 ( x1 , x2 ) = 1−  (e 2x1 − e) + 2 − 2ex1 = 0. ... (2.16)
 4π  π
it is readily seen that
∂f1 1 x cos(x1 , x2 ) ∂f1 1 x cos(x1 , x2 )
=− 2 , =− + 1 ,
∂x1 2 2 ∂x1 4π 2
∂f 2  1  ∂f 2 e
= −2e +  2 − e 2 x1 , = ,
∂x1  2π  ∂x2 π

The increments ∆x1 and ∆x2 in x1 and x2 are determined by

∂f 1 ∂f
∆x1 + 1 ∆x 2 = − f 1 ,
∂x1 ∂x 2
∂f 2 ∂f 2
∆x1 + ∆x 2 = − f 2 .
∂x1 ∂x 2
Or, writing the determinant D of the coefficient matrix (the Jacobian),

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 7

 ∂f1 ∂f 2 ∂f1 ∂f 2
D= − ,
∂x1 ∂x 2 ∂x 2 ∂x1
then
 ∂f ∂f   ∂f 2 ∂f 
 f 2 1 − f1 2   f1 − f1 1 
∂ ∂ ∂ ∂
∆x1  2 
, ∆x 2 
x 2 x x1 x1 
.
 D   D 
   
   
for case in verification, detailed results are tabulated in Table 3.1 and moreover,
calculations were carried out using slide rule and the entries 0.00000 showed
tiny negative values.

APPLICATION OF NEWTON- RAPHSON’S METHOD IN SOLVING

THE CHEMICAL EQUILIBRUIM PROBLEM
The principal reactions in the production of synthesis gas by partial oxidation of
methane with oxygen are:
CH4 + ½ O2 → CO + 2H2 → - - -
(3.1)
CH4 + H2O → CO + 3H2 → - - -
(3.2)
H2 + CO2 → C0 +H2O → - - -
(3.3)
Write a program that finds the 0 reactant ratio that will produce an adiabatic
equilibrium temperature of 22000 F at an operating pressure of 20 atmospheres,
when the reactant gases are preheated to an entering temperature of 10000F.
Assuming that the gases behave ideally, so that the component
activities are identical with component partial pressures, the equilibrium
constants at 22000F for the three equations are respectively:

PcoP 2 H 2
K1 = = 1 .3 x10 11 ↓ → ... (3.4 )
PCH P 1 / 2
02

PcoP 3 H 2
K2 =
P P
= 1 .7837 x10
5
→ ... (3.5)
CH 4 H 2O

PCO PH
20
K3 = = 2.6058 → ( 3 .6 )
PCO PH
2 2

Here PCO, Pco2, PH2O, PCH4 and PO2 are the partial pressured of CO (carbon
monoxide), CO2 (carbon dioxide), H2O (water vapor), H2 (hydrogen), CH4
(methane), and O2 (oxygen), respectively. Enthalpies of the various components
at 10000F and 22000F are listed in Table (3.1)

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 8

 Table (3.1) Component Enthalpies in BTU/b mole
Component 10000F 22000F
CH4 -13492 8427
H 2O -90546 -78213
CO2 -154958 -139009
CO -38528 -28837
H2 10100 18927
O2 10690 20831

A fourth reaction may also occur at high temperatures:

C+ CO2 2CO (3.1) at 22000F, any carbon formed would be deposited as a solid;
the equilibrium constant is given by

P 2 co
K 4 = = 1 . 7837 x 10 5
(3 . 8 )
ac Pco 2

where ac is the activity of carbon in the solid state. Do not include reaction (3.7)
in the equilibrium analysis. After establishing the equilibrium composition,
considering only the homogeneous gaseous reactions given by (3.1), (3.2), and
(3.3), determine the thermodynamic likelihood that solid carbon would appear as
a result of reaction (3.7). Assume that the activity of solid carbon is unaffected
by pressure and equals unity.
Use the Newton- Raphson method to solve the system of simultaneous
nonlinear equations developed as the result of the equilibrium analysis.

3.2 METHOD OF SOLUTION

Because of the magnitude of K, the equilibrium constant for reactions, the first
reaction can be assumed to go to completion at 22000F, that is virtually no
unrelated oxygen will remain in the product gases at equilibrium.
Let the following nomenclature be used.
x1 mole fraction of CO in the equilibrium mixture
x2 mole fraction of CO2 in the equilibrium mixture
x3 mole fraction of H2O in the equilibrium mixture
x4 mole fraction of H2 in the equilibrium mixture
x5 mole fraction of CH4 in the equilibrium mixture
x6 number of moles of O2 per mole of CH2 in the feed gas
x7 number of moles of product gases in the equilibrium mixture per
mole of CH4 in the feed gases.
Then a system of seven simultaneous equations may be generated from three
atom balances an energy balance, a mole fraction constraint and two equilibrium
relations.
Atom conservation balances: the number of atoms of each element
entering equals the number of atoms of each clement in the equilibrium mixture.
Oxygen: x6 = (1/2x1 + x2 + 1/2x3 → … (3.9)
Hydrogen: 4 = (2x3 + 2x4 + 4x2 →
… (3.10)
Carbon: 1 = (x1 + x2 + x5 → … (3.11)

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 9

Since the reaction is to be conducted adiabatically, that is, no energy is added to
or removed from the reacting gases, the enthalpy (H) of the reactants must equal
the enthalpy of the products.
[HCH4 + x6H02]1000’F = x7[x1HC02 + x3H20 + x4H2 + x5HCH2]2200’F … (3.12)
mole fraction constraint.
x1 + x2 + x4 + x5 + = 1 … (3.13)
Equilibrium relations.
2 3
p X 1X 4
K 2 = = 1 . 7837 x 10 5
→ ... ( 3 . 14 )
x3 x5
X 1X
K 3 = 3
= 2 . 6058 → ... ( 3 . 15 )
X 2X 4

The relationships (3.14) and (3.15) follow directly from (3.5) and (3.6),
respectively, where P is the total pressure and Pco = Px1, etc. in addition, there
are five side conditions.
x1 ≥ 0, I = 1, 2, 5 … (3.16)
These C ions more that all mole fractions in the equilibrium mixture are
nonnegative, that is, any solution of equation (3.9) to 3.15) that contains
negative mole fractions is physically meaningless from physical-chemical
principle there is one and only one solution of the equation that satisfies
conditions (3.16). Any irrelevant solutions may be detected easily.

The seven equations may be rewritten in the form.

f1 (x) = 0, I = 1, 2,…, 7 (5.5.17)
where x = [x1x2 x3 x4 x5 x6 x7]t, as follows:

1 1 x
f1 ( x ) = x1 + x 2 + x3 − 6 = 0 → ( 3 . 19 a )
2 2 x7
2
f2 (x) = x3 + x4 + 2x5 − = 0 → ( 3 . 19 b )
x7
1
f 3 ( x ) = x1 + x 2 + x 5 − = 0 → ( 3 . 19 c )
x7
f 4 ( x ) = − 28837 x 1 − 139009 x 2 − 78213 x 3 + 18927 x4
13492 x6
+ 8427 x5 + − 10690 = 0 → ( 3 . 19 d )
x7 x7
f 5 ( x) = x1 x 2 x 3 x 4 x 5 − 1 = 0 → ( 3 . 19 e )

f6 (x) = P 2
x 1 x 4 − 1 . 7837 x 10 5
x3x5 = 0 → ( 3 . 19 f )

f 7 ( x ) = x 1 x 3 − 2 . 6058 x2x4 = 0 → ( 3 . 19 g )

The system of simultaneous nonlinear equations has the form (2.1), and
will be solved using the Newton- Raphson method, described in section 2.2.
The partial derivatives of above may be found by partial differentiation of the
seven functions, f1 (x), with respect to each of the seven variables. For example,

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 10

 ∂f1 1 ∂f1 ∂f1 x
= , = 0, = 62
∂ x1 2 ∂x 4 ∂x7 x7
∂f1 ∂f1
= 1, = 0,
∂ x 21 ∂x5
∂f1 1 ∂f1 1
= , = − ,
∂x3 2 ∂x 6 x7
The Newton- Raphson method may be summarized as follows:
Choose a starting vector xK = x0 = [x10, x20, …, x70], where x0 is hopefully near
a solution
solve the system of linear equations (2.14),
φ (xk)δk = - f(xk),
where
φij (xk) = ∂ f 1 ( x ) i = 1 , 2 ,..., 7 ,
( 3 . 20 )
∂x j
k
j = 1 , 2 ,..., 7 ,
f ( x k ) = lf 1 ( x k ), f 2 ( x k ),....., f7 ( xk )t , ( 2 . 21 )
for the increment vector
t
 
δ k = δ δ  ( 3 . 22 )
1k 2 k ,......, δ
 7 k 
update the approximation to the root for the next iteration. xk+1 = xk + δk.
check for possible convergence to a root α. One such test might be
| δ ik| < ε2, i = 1,2,…,7. … (3.23)
if (3.23) is true for all i, then xk +1 is taken to be the root. Is test (3.23) is failed
for any 1, then the process is repeated starting with step 2. The iterative process
is continued until test (3.23) is passed for some k, or when k exceeds some
specified upper limit. In the programs that follow, the elements of the augments
matrix A=[φ(xk): - f(xk)] … (3.24) Are evaluated by a subroutine named
CALCN. The system of linear equations (3.24) is solved by calling on the
function SIMUL, described in detail in example (2.1)
The main program is a general one, in that it is not specifically written to solve
only the seven equations of interest. By properly defining the subroutine
CALCN, the main program could be used to solve any system of n simultaneous
nonlinear equations. The main program reads data values for itmax, iprint, n, Σ1,
Σ2, and x1, x2,…, xn here, itmax is the maximum number of Newton- Raphson’s
iterations, print is a variable that controls printing of intermediate output, n is the
number of nonlinear equations, Σ1, is the minimum pivot magnitude allowed in
the Gauss-Jordan reduction algorithm, Σ2, is a small positive number used in test
(3.23), and x1,x20,….,xno, that is, the elements of x0.

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 11

 FLOW DIAGRAM
Main Program

Itmax, iprint Evaluate elements

K=1,2,…, itmax,
Beg n, Σ1, Σ2 aij of augmented
in x10,…,xno matrix
A (see (5.5.24)).
(Subroutine

“Matrix ill-
End T Solve system of n
conditioned
linear equations
T d=
(5.44) for the
F increments
i= Itcon δ ik, δ 2k,…, δ nk
and determinant, d.
(function SIMUL)

| δ ik| <T Σ2 Itcon 0

“Convergen
ce” k, d, n, End
xi,k+1,....,xn,k +1
5
xi,k+1 T

F “No
itcon= 9 End
End Convergenc
e”

Subroutine CALC (Arguments: xk, A, N)

Calculate elements
End aij i = 1, 2,…., n Retur
j = 1, 2,…..n+1 of
matrix A
(see(5.5.24)).

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 12

Fortran Implementation
List of Principle various program Symbol
(Main)
A Augmented matrix of coefficients, A (see (3.22).
DETER d, determinant of the matrix φ (the jacobian).
EPS1 Σ1, minimum pivot magnitude permitted in subroutine SIMUL.
EPS2 Σ2, small positive number, used in convergence test (3.23).
subscript, i.
IPRINT Print control variable, if iprint = 1, intermediate solutions are
printed after each iteration.
ITCON used in convergence test (3.23). ITCON 1 if (3.23) is passed
for all i, i = 1, 2,…..,n: otherwise ITCON = 0.
ITER Iteration counter, k.
ITMAX maximum number of iterations permitted, itmax.
N number of nonlinear equations, n.
XINC vector of increments, ik, i = 1, 2, ….,n.
XOLD vector of approximations to the solution, xik.
SIMUL function developed in Example (2.1) solves the system of n
linear equations (2.15) for the increments,. ik i = 1, 2,…,n.

(subroutine CALCN)
DXOLD same as XOLD. Used to avoid an excessive number of
reference to subroutine arguments in CALCN.
I, J, i and j, row and column subscript, respectively.
NRC N, dimension of the matrix A in the calling program. A
is
assumed to have the same number of rows and
columns.
P pressure, P, atm.

Program Listing(Chemical Equilibrium by Newton- Raphson Method)

Main Program
c APPLIED NUMERICAL METHODS, EXAMPLES 5.5
C CHEMICAL EQUILIBRIUM- NEWTON-RAPSION METHOD
C
C THIS PROGRAM SOLVES N SIMULTANEOUS NON-LINEAR
EQUATIONS
C IN N UNKNOWNS BY THE NEWTON-RAPSON ITERATIVE
PROCEDURE
C INITIAL GUESSES FOR VALUES OF THE UNKNOWNS ARE
READ INTO
C XOLD (1)…….XOLD (N). THE PROGRAM FIRST CALLS ON
THE SUBROUTINE
C CALCN TO COMPUTE THE ELEMENTS OF A. THE
AUGMENTED MATRIX OF

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 13

C PARTIAL DERIVATIVES, THEN ON FINCTION SIMUL (SEE
PROBLEM 5.2)
C TO SOLVE THE GENERATED SET OF LINEAR EQUATIONS
FOR THE CHANGES
C IN THE SOLUTION VALUES XINC(1)……XINC(N). DETER
IS THE
C JOCABIAN COMPUTED BY SIMUL. THE SOLUTIONS ARE
UPDATED AND THE
C PROCESS CONTINUED UNTIL ITER, THE NUMBER OF
ITERATIONS,
C EXCEEDS ITMAX OR UNTIL THE CHANGE IN EACH OF
THE N VARIABLES
C IS SMALLER IN MAGNITUDE THAN EPS2 (ITCON =1
UNDER THESE
C CONDITIONS). EPS1 IS THE MINIMUM PIVOT MAGNITUDE
PERMITTED
C IN SIMUL. WHEN IPRINT =1, INTERMEDIATE SOLUTION
VALUES ARE
C PRINTED AFTER EACH ITERATION.
C
DIMENSION XOLD (21), XINC (21), A (21,21)
C
C ….READ AND PRINT DATA…
1 READ (5,100) ITMAX, IPRINT, EPS1, EPS2, (XOLD(1), 1=1, N)
WRITE(6,200) ITMAX, IPRINT, N, EPS1, EPS2, N, (XOLD)(1),
1=1,N)
C
C … CALL ON CALCN TO SET UP THE A MATRIX…
CALL CALCN (XOLD, A, 21)
C
C … CALL SIMUL TO COMPUTE JACOBIAN AND CORRECTIONS
IN XINC…
DETER = SIMUL (N, A, XINC, EPS1, 1, 21)
IF ( DETER. NE. D. ) GO TO 3
WRITE (6 GO TO
C
C …CHECK FOR CONVERGENCE AND UPDATE XOLD VALUES…
3 ITCON = 1
DO 5 1 – 1, N
IF (ABS(XINC(1)) . GT. EPS 2 ) ITCON = 0
5 XOLD (1) = XOLD(1) + XINC(1)
I F ( IPRINT. EQ. 1 ) WRITE (6, 202) ITER, DETER, N, (XOLD (1), 1 =
1,N)
I F ( ITCON. EQ. 0 ) ITER, N, (XOLD(1), 1 = 1, N)
WRITE (6,203) ITER, N, (XOLD(1), 1 =1,N)
GO TO 1
9 CONTINUE
C
WRITE 96, 204)

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 14

PROGRAM Listing (continued)
Subroutine CALCN
SUBROUTINE CALCN D XOLD, A, NRC )
C THIS SUBROUTINE SETS UP THE AUGMENTED MATRIX OF
PARTIAL
C DERIVATIVES REQUIRED FOR THE SOLUTION OF THE
NON-LINEAR
C EQUATIONS WHICH DESCRIBE THE EQUILIBRIUM
CONCENTRATIONS
C OF CHEMICAL CONSTITUENTS RESULTS FROM PARTIAL
OXIDATION
C OF METHANE WITH OXYGEN TO PRODUCE SYNTHESIS
GAS. THE PRESSURE
C IS 20 ATMOSPHERES. SEE TEXT FOR MENINGS OF
XOLD(1)…XOLD(N)
C AND A LISTING OF THE EQUATIONS. DXOLD HAS BEEN
USED AS THE
C DUMMY ARGUMENT FOR XOLD TO AVOID AN
EXCESSIVE NUMBER OF
C REFENCENCES TO ELEMENTS IN THE ARGUMENT LIST.
C
DIMENSION XOLD(20), DXOLD(NRC), A(NRC , NRC)
C
DATA P / 20. /
C
C …SHIFT ELEMENTS OF DXOLD TO XOLD AND CLEAR A
ARRAY…
DO 1 I =1, 7
XOLD(1) = DXOLD(1)
DO 1 J = 1, 8
1 A(1 , J) = 0.
C
C …COMPUTE NON-ZERO ELEMENTS OF A…
A (1,1) = 0.5
A (1,2) = 1.0
A (1,3) = 0.9
A (1,6) = 1.0 / XOLD (7)
A (1,7) = XOLD (6) / XOLD (7)**2
A (1,8) = - XOLD(1) / 2. – XOLD (2) - XOLD (3) /2, * XOLD (6) /
XOLD(7)
A (2,3) = 1.0
A (2,4) = 1.0
A (2,5) = 2.0
A (2,7) = 2.0 / XOLD (7) **2
A (2,8) = XOLD(1) – XOLD(2) – XOLD(5) + 1.0 / XOLD(7

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 15

 A (3,1) = 1.0
A (3,2) = 1.0
A (3,5) = 1.0
A (3,7) = 1.0 / XOLD(7) **2
A (3,8) = - XOLD(1) – XOLD(2) – XOLD(5) + 1.0 / XOLD(7)
A (4,1) = - 28837.
A (4,2) = - 139009.
A (4,3) = - 78213.
A (4,4) = 18927
A (4,5) = 8427
A (4,6) = - 10690. / XOLD (7)
A (4,7) = (- 13492. + 10690. * XOLD (6) ) / XOLD (7)**2
A (4,8) = 28837. * XOLD(4) – 8427. * XOLD(5) – 13492. / XOLD 97) + 10690.
1 -18927. * XOLD(4) - 8427. * XOLD(5) – 13492. / XOLD 97) + 10690.
2 * XOLD(6) / XOLD(7)
A (5,1) = 1.0
A (5,2) = 1.0
A (5,3) = 1.0
A (5,4) = 1.0
A (5,5) = 1.0
A (5,8) = 1.0 – XOLD(1) – XOLD(2) - XOLD(3) - XOLD(4) - XOLD(5)
A (6,1) = P*P* XOLD(4)**3
A (6,3) = - 1.7837E5 * XOLD (5)
A (6,4) = 3.0 *P*P* XOLD(1) * XOLD(4)**2
A (6,5) = - 1.7837E5* XOLD (3)
A (6,8) = 1.7837E5*XOLD(5)*XOLD(3) – P*P*XOLD(1)*XOLD(4)**3
A (7,1) = XOLD (3)
A (7,2) = - 2.6058*XOLD(4)
A (7,3) = XOLD(1)
A (7,4) = - 2.6058* XOLD(2)
A (7,8) = 2.6058* XOLD(4) *XOLD(2) - XOLD(1)*XOLD(3)
RETURN
C
END

Chemical Equilibrium (Newton- Raphson Method)

Program listing (continued)
Data
ITMAX = 50 IPRINT = 1 N =
7
EPS1 = 1.0E – 10 EPS2 = 1. 0E – 05
XOLD (1)… XOLD(5) = 0.500 0.000 0.000
0.500 0.000
XOLD(6)….XOLD(7) = 0.500 2.000
ITMAX = 50 IPRINT = 0 N
= 7
EPS1 = 1.0E-10 EPS2 = 1.0E – 05
XOLD(1)….XOLD(5) = 0.200 0.200 0.200
0.200 0.200
XOLD(6)….XOLD(7) = 0.500 2.000
ITMAX = 50 IPRINT = 0 N
= 7

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 16

 EPS1 = 1.0E - 10 EPS2 = 1.0E -05
XOLD(1)….XOLD(5) = 0.220 0.075
0.001 0.580 0.125
XOLD(6)….XOLD(7) = 0.436 2.350

Computer Output
Results for the 1st Data Set

ITMAX = 50
IPRINT = 1
N = 7
EPS1 = 1.0E -10
EPS2 = 1. 0E -05
XOLD(1)....XOLD ( 7)
5.000000E -01 0.0 0.0
5.000000E-01
0.0 5.000000E -01 2.000000E 00
ITER = 1
DETER = -0.97077E 07
XOLD(1)...XOLD( 7)
2.210175E -01 2.592762E -02 6.756210E -02
4.263276E -01
2.591652E -01 3.3432350E -01 1.975559E 00
ITER = 2
DETER = -0.10221E 10
XOLD(1)…XOLD( 7)
3.101482E -01 7.142063E -03 5.538273E -02
5.791981E-01
4.812878E -02 4.681466E -01 2.524948E 00
ITER = 3
DETER = -0.41151E 09
XOLD(1)…XOLD( 7)
3.202849E -01 9.554777E -03 4.671279E -02
6.129664E -01
1.048106E-02 5.533223E -01 2.880228E 00
ITER = 4
DETER = -0.22807E 09
XOLD(1)…XOLD( 7)
3.228380E -01 9.22480E -03 4.603060E -02
6.180951E-01
3.811378E -03 5.758237E -01 2.974139E 00
ITER = 5
DETER = -0.20218E 09
XOLD(1)…XOLD( 7)
3.228708E -01 9.223551E -03 4.601710E -02
6.181716E -01

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 17

 3.716873E -03 5.767141E -01 2.977859E 00
ITER = 6
DETER = -0.20134E 09
XOLD(1)…XOLD( 7)
3.228708E -01 9.223547E -03 4.601710E -02
6.181716E -01
3.716847E -03 5.767153E -01 2.977863E 00

Computer Output
SUCCESSFUL CONVERGENCE
ITER = 6
XOLD(1)…XOLD( 7)
3.228708E -01 9.223547E -03
4.601710E -02 6.181716E -01
3.716847E -03 5.767153E-01 2.97863E
00
Results for the 3rd Data Set
ITMAX = 50
IPRINT = 1
N = 7
EPS1 = 1.0E -10
EPS2 = 1. 0E -05
XOLD(1)...XOLD ( 7)
2.200000E -01 7.499999e -02
9.999999e -04 5.800000E-01
1.250000e -01 4.360000e -01 2.349999e
00
ITER = 1
DETER = -0.61808E 08
XOLD(1)...XOLD( 7)
6.9514955E -01 -8.022028E -02 1.272939E
-02 1.217132E 00
-8.447912E -01 1.314754E 00 5.969404E
00
ITER = 2
DETER = 0.12576E 09
XOLD(1)…XOLD( 7)
4.958702E -01 -1.698154E -02 5.952045E
-03 9.518250E -01
-3.65007E -01 2.379797E 00
1.043425E 01
ITER = 3
DETER = 0.77199E 07
XOLD(1)…XOLD( 7)

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 18

 4.559822E -01 -9.799302E -04 -7.583648E -04
9.107630E -01
-3.650070E -01 2.509821E 00 1.107038E 01
ITER = 4
DETER = 0.53378 07
XOLD(1)…XOLD( 7)
4.569673E -01 -4.071472E -04 -2.142648E -03
9.152630E -01
-3.696806E -01 2.608933E 00 1.149338E 01

ITER = 5
DETER = 0.49739E 07
XOLD(1)… XOLD( 7)
4.569306E -01 -4.071994E -04 -2.125205E -03
9.151721E -01
-3.695704E -01 2.610552E 00 1.150046E 01

ITER = 6
DETER = 0.49611E 07
XOLD(1)…XOLD( 7)
4.569306E -01 -4.071984E -04 -2.125199E -03
9.151720E-01
-3.695703E -01 2.610549E 00 1.150045E 01
SUCCESSFUL CONVERHENCE
ITER = 6
XOLD(1)…XOLD( 7)
4.569306E -01 -4.071984E -04 -2.125199E -03
9.151720E -01
-3.695703R -01 2.610549E 00 1.150045E 01

Tables (3.2) Equilibrium Gas mixture

x1 Mole fraction CO 0.322871
x2 Mole fraction CO2 0.009224
x3 Mole fraction H2O 0.046017
x4 Mole fraction H2 0.618172
x5 Mole fraction CH4 0.003717
x6 Mole fraction O2 / CH4 0.576715
x7 Total moles of product 2.977863

In the feed gases, and total number of moles of product per mole of HC4 in the
feed are tabulated in Table (3.2). Thus the required feed ratio is 0.5767 moles of
oxygen per moles of methane in the feed gases.
To establish if carbon is likely to be formed according to reaction (5.5.7) at
22000F for a gas of the computed composition, it is necessary to calculate the
magnitude of
v P co Px2 2

K=a P =a x
1
... (3.25)
c co2 c 2

If
v
K is larger than k4 from (3.25), then there will be a tendency for reaction
(3.24) to shift toward the left; carbon will be formed. Assuming that ac = 1,

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 19

v 20x(0.322871) 2

k = 1x0.009224 = 226.03 < k4 = 1329.5. − − − (3.26)

Therefore there will be no tendency for carbon to form.

† Okereke, C. Emmanuel, Ph.D., Department of Mathematics, Michael Okpara

University of Agriculture, Nigeria

REFERNCES

[1] Alefeld, G., Potra, F. A and Shen, z, on the existence theorems of

kontorovich, more and mirands, preprint 2001.
[2] Ames, W., Nonlinear ordinary differential equations in transport
processes. Academic press, New York 1968.
[3] Agyros, I. K. on an interative algorithon for solving non linear
equations beritrage zir numerschen math, 10, 1 (1991) 83-92.
[4] Agyros, I. K., some generalized projection methods for solving
operator equations, J. comp. Appl math 39, 1 (192) 1-6.
[5] A Ralston and H. S., wilf, Mathematical Methods for digital computers,
Wiley, New York. (1960).
[6] Brice, carnahan, luther, H. A., Wilkes O. James. Applied numerical
methods John Wiley and sons in cooperated, New York London. Sydney-
Toronto.
[7] F. A. Jenkins and H. E. White fundamentals of optics 2nd Edition.
[8] J. F. traub, iterative methods for the solution of equations. Prentice-
Hall, eagle wood cliffa. New Jessey 1964.
[9] Mc Graw – Hill, New York 1951.
[10] Ortega, J. m, and Rheinboldt, W. C. iterative solution of non linear
equations in several variable academic press, New York, 1970.
[11] Rall, L. B. computational solution of non linear operator egns. Wiley
New York 1968.

Journal of Mathematical Sciences & Mathematics Education Vol. 9 No. 2 20