Generalized False Position (GFP)

The GFP Function Block type implements the Generalized False Position algorithm for the solution of systems of simultaneous linear or non-linear equations, as described below.

Quasi-Newton Methods

Quasi-Newton methods are iterative numerical procedures for finding a zero of the vector function f(x), that is, a solution to the set of simultaneous equations

f(x) = 0,
where
f(x) = (f1(x), ..., fn(x))
and
x= (x1, ..., xn).

These methods are based on the approximation

f(x) - f(k) ≈ J(k)(x - x(k)),

where x(k) is the value of x at the most recent iteration k, f(k)=f(x(k)) and J(k) is the Jacobian matrix of partial derivatives of f with respect to x at x(k):

J(k)ij = ∂fi/∂xj|x=x(k).

Quasi-Newton methods apply various techniques to estimate the value of J(k) or (J(k))-1, and then obtain the next iterated value x(k+1) by setting f(k+1) to 0 and solving for x(k+1):

x(k+1) = x(k) - (J(k))-1f(k).

As early as 1962, Rubin[ 1] noted that the Jacobian could be replaced with a global linear fit to the previous n+1 iterations. Rubin called this "Generalized False Position (GFP)" since it was a multivariable generalization of the false position method. However, Rubin's formulation requires a full matrix inversion at every step. We will use instead a formulation that allows us just to use a single Product Form of the Inverse ( PFI) pivot operation at each step (an analysis of the convergence properties of a similar method was given in 1995 by Lopes and Martínez).

Linear Estimation of the Jacobian Matrix

Let A be an n×n matrix of coefficients for the linear approximation f(x) ≈ Ax + b, that is,

fi(x) ≈ ai1x1 + ... + ainxn + bi
for i = 1, ..., n.

If the coefficients aij are determined in such a manner that the linear fit is exact at points x(k) and x(k-1), then we can write

Δf(k) = AΔx(k),
where
Δf(k) = f(k) - f(k-1)
and
Δx(k) = x(k) - x(k-1).

Let ΔX be an n×n matrix whose columns {1,...,n} serve as a circular buffer for n successive differences Δx(k), where after each iteration Δx(k)
replaces column ((k-1)mod n)+1 , and let ΔF be a similar circular buffer for values of Δf(k). Then

ΔF(k) = A(k)ΔX(k),
whence
(A(k))-1 = ΔX(k)(ΔF(k))-1
and
x(k+1) = x(k) - (A(k))-1f(k)
= x(k) - ΔX(k)(ΔF(k))-1f(k).

Since ΔF(k) only varies by a single column from ΔF(k-1), it is not necessary to do a complete matrix inversion at each iteration; rather, the matrix B(k) = (ΔF(k))-1 can be updated at each iteration via a Product Form of the Inverse ( PFI) pivot operation and the iterative step is then defined as

x(k+1) = x(k) - ΔX(k)B(k)f(k).

The GFP Algorithm

  1. Input the vector size n, initial vector x(0), perturbation vector p, convergence tolerance ε, and maximum number of iterations kMAX .
  2. Initialize B to an identity matrix of order n, and ΔX to a square zero matrix of order n.
  3. Set k = 0 and compute f(0) = f(x(0)) .
  4. If |f(k)| < ε, STOP.  x(k) is the desired solution.
  5. If k = kMAX, STOP.  The algorithm has failed to converge.
  6. Set k = k + 1.
  7. Compute Δx(k):
    • If k ≤ n, Δx(k) = pkîk, where îk is the unit vector whose j-th element is the given by the Kronecker delta function δjk.
    • Otherwise, set Δx(k) = -ΔX * Bf(k-1) .
  8. Set x(k) = x(k) + Δx(k)
  9. Let r = ((k-1) mod n) + 1.
  10. Replace column r of ΔX with Δx(k).
  11. Compute f(k) = f(x(k)) and Δf(k) = f(k) - f(k-1).
  12. Perform a PFI pivot of Δf(k) into column r of B.
  13. Repeat from Step 4.

[1] Donald I. Rubin, "Generalized Material Balance," Chemical Engineering Progress Symposium Series 37, 58 (1962), pp.54-61.