%% This document created by Scientific Word (R) Version 2.5


\documentclass[12pt,thmsa]{article}
\usepackage{amssymb}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{sw20jart}

%TCIDATA{TCIstyle=article/art4.lat,jart,sw20jart}

%TCIDATA{Created=Sun Feb 13 21:16:24 2000}
%TCIDATA{LastRevised=Fri Feb 25 11:41:16 2000}

\input{tcilatex}
\begin{document}


\parindent0cm

%\Large

\begin{center}
{\Huge \textbf{Nonlinear Equations}}
\end{center}

\vspace{0.5cm}

\begin{center}
W. E. Schiesser\\[0pt]
Iacocca D307\\[0pt]
111 Research Drive\\[0pt]
Lehigh University\\[0pt]
Bethlehem, PA 18015\\[0pt]

\vspace{0.5cm}

(610) 758-4264 (office)\\[0pt]
(610) 758-5057 (fax)\\[0pt]

\vspace{0.5cm}

wes1@lehigh.edu \\[0pt]
\vspace{0.2cm} http://www.lehigh.edu/\~{\thinspace}wes1/wes1.html
\end{center}

\newpage

\begin{center}
{\Huge \textbf{Nonlinear Equations}}
\end{center}

\vspace{0.2cm}

\begin{center}
http://www.lehigh.edu/\~{\thinspace }wes1/apci/24feb00-b.ps\\[0pt]
http://www.lehigh.edu/\~{\thinspace }wes1/apci/24feb00-b.pdf
\end{center}

\vspace{0.5cm}

\begin{center}
http://www.lehigh.edu/\~{\thinspace }wes1/apci/24feb00-s.tex\\[0pt]
http://www.lehigh.edu/\~{\thinspace }wes1/apci/24feb00-s.ps\\[0pt]
http://www.lehigh.edu/\~{\thinspace }wes1/apci/24feb00-s.pdf
\end{center}

\vspace{0.5cm}

Previous lectures are in:

\begin{center}
http://www.lehigh.edu/\~{\thinspace }wes1/apci/28jan00.tex\\[0pt]
http://www.lehigh.edu/\~{\thinspace }wes1/apci/28jan00.ppt\\[0pt]
http://www.lehigh.edu/\~{\thinspace }wes1/apci/28jan00.ps\\[0pt]
http://www.lehigh.edu/\~{\thinspace }wes1/apci/28jan00.pdf\\[0pt]
\end{center}

\vspace{0.5cm}

\begin{center}
http://www.lehigh.edu/\~{\thinspace }wes1/apci/18feb00.ps\\[0pt]
http://www.lehigh.edu/\~{\thinspace }wes1/apci/18feb00.pdf
\end{center}

A Table of Contents is in:

\begin{center}
http://www.lehigh.edu/\~{\thinspace }wes1/apci/files.txt
\end{center}

\vspace{0.5cm}

\begin{enumerate}
\item  Newton's method and quasilinearization (WES)

\item  Obtaining derivatives (LTB)

\item  Broyden's method (LTB)

\item  Continuation methods (WES)

\item  Davidenko's equation (WES)

\item  Levenberg-Marquardt method (WES)

\item  Tearing algorithms (LTB)

\item  Direct substitution (LTB)

\item  Wegstein method (LTB)

\item  References and software (LTB/WES)
\end{enumerate}

\newpage

\begin{center}
{\Huge \textbf{The General Nonlinear Problem}}
\end{center}

\vspace{0.2cm}The n x n nonlinear problem is: 
\begin{equation}
\begin{array}{c}
f_{1}(x_{1},x_{2},\cdots ,x_{n})=0 \\ 
f_{2}(x_{1},x_{2},\cdots ,x_{n})=0 \\ 
\vdots \\ 
f_{n}(x_{1},x_{2},\cdots ,x_{n})=0
\end{array}
\tag{1}
\end{equation}

which can be summarized as

\begin{equation}
\mathbf{f}(\mathbf{x})=\mathbf{0}  \tag{2}
\end{equation}

where we use bold face to denote vectors:

\[
\mathbf{f}=\left[ 
\begin{array}{c}
f_{1} \\ 
f_{2} \\ 
\vdots \\ 
f_{n}
\end{array}
\right] ,\mathbf{x}=\left[ 
\begin{array}{c}
x_{1} \\ 
x_{2} \\ 
\vdots \\ 
x_{n}
\end{array}
\right] ,\mathbf{0}=\left[ 
\begin{array}{c}
0 \\ 
0 \\ 
0 \\ 
0
\end{array}
\right] 
\]

Two points:

\vspace{0.5cm}

\begin{itemize}
\item  Nonlinear equations are generally solved numerically by the iterative
solution of linear equations

\vspace{0.5cm}

\item  John Rice states (approximately) that the solution of nonlinear
equations is the most difficult problem in scientific computation
\end{itemize}

\newpage

\begin{center}
{\Huge \textbf{Newton's Method and Quasilinearization}}
\end{center}

\vspace{0.5cm}

To obtain a form of Newton's method for systems of equations, we start with
the Taylor series in $n$ variables, e.g., $n=2$, we want to obtain the
values of $x_{1},x_{2}$ that satisfy the 2 x 2 system

\begin{equation}
\begin{array}{c}
f_{1}(x_{1},x_{2})=0 \\ 
f_{2}(x_{1},x_{2})=0
\end{array}
\tag{7}
\end{equation}

If $f_{1}$ and $f_{2}$ are expanded in two Taylor series in two dimensions
(with respect to the iteration counter $k$)

\begin{equation}
\begin{array}{c}
f_{1}(x_{1}^{k+1},x_{2}^{k+1})=f_{1}(x_{1}^{k},x_{2}^{k})+\dfrac{\partial
f_{1}(x_{1}^{k},x_{2}^{k})}{\partial x_{1}}(x_{1}^{k+1}-x_{1}^{k})+\dfrac{%
\partial f_{1}(x_{1}^{k},x_{2}^{k})}{\partial x_{2}}(x_{2}^{k+1}-x_{2}^{k})+%
\cdots \\ 
f_{2}(x_{1}^{k+1},x_{2}^{k+1})=f_{2}(x_{1}^{k},x_{2}^{k})+\dfrac{\partial
f_{2}(x_{1}^{k},x_{2}^{k})}{\partial x_{1}}(x_{1}^{k+1}-x_{1}^{k})+\dfrac{%
\partial f_{2}(x_{1}^{k},x_{2}^{k})}{\partial x_{2}}(x_{2}^{k+1}-x_{2}^{k})+%
\cdots
\end{array}
\tag{8}
\end{equation}

truncation after the first derivative terms, with the resulting equations
expressed in terms of Newton corrections

\[
\begin{array}{c}
\Delta x_{1}^{k}=(x_{1}^{k+1}-x_{1}^{k}) \\ 
\Delta x_{2}^{k}=(x_{2}^{k+1}-x_{2}^{k})
\end{array}
\]

and with just the right Newton corrections so that the LHS's of eqs. (8) are
zero, gives

\begin{equation}
\begin{array}{c}
\dfrac{\partial f_{1}(x_{1}^{k},x_{2}^{k})}{\partial x_{1}}\Delta x_{1}^{k}+%
\dfrac{\partial f_{1}(x_{1}^{k},x_{2}^{k})}{\partial x_{2}}\Delta
x_{2}^{k}=-f_{1}(x_{1}^{k},x_{2}^{k}) \\ 
\dfrac{\partial f_{2}(x_{1}^{k},x_{2}^{k})}{\partial x_{1}}\Delta x_{1}^{k}+%
\dfrac{\partial f_{2}(x_{1}^{k},x_{2}^{k})}{\partial x_{2}}\Delta
x_{2}^{k}=-f_{2}(x_{1}^{k},x_{2}^{k})
\end{array}
\tag{9}
\end{equation}

Eqs. (9) are a 2 x 2 linear system in the two Newton corrections, $\Delta
x_{1}^{k},\Delta x_{2}^{k}$ then can be solved by any standard method for
linear equations, e.g., Gauss row reduction. The unknowns can then be
updated iteratively by

\begin{equation}
\begin{array}{c}
x_{1}^{k+1}=\Delta x_{1}^{k}+x_{1}^{k} \\ 
x_{2}^{k+1}=\Delta x_{2}^{k}+x_{2}^{k}
\end{array}
k=0,1,2,\cdots  \tag{10}
\end{equation}

Eq. (9) can be generalized to the n x n case by writing it in matrix form

\begin{equation}
J\mathbf{\Delta x}=-\mathbf{f}  \tag{11}
\end{equation}

where

\begin{equation}
J=\left[ 
\begin{array}{cccc}
\dfrac{\partial f_{1}(x_{1}^{k},x_{2}^{k},\cdots ,x_{n}^{k})}{\partial x_{1}}
& \dfrac{\partial f_{1}(x_{1}^{k},x_{2}^{k},\cdots ,x_{n}^{k})}{\partial
x_{2}} & \cdots  & \dfrac{\partial f_{1}(x_{1}^{k},x_{2}^{k},\cdots
,x_{n}^{k})}{\partial x_{n}} \\ 
\dfrac{\partial f_{2}(x_{1}^{k},x_{2}^{k},\cdots ,x_{n}^{k})}{\partial x_{1}}
& \dfrac{\partial f_{2}(x_{1}^{k},x_{2}^{k},\cdots ,x_{n}^{k})}{\partial
x_{2}} & \cdots  & \dfrac{\partial f_{2}(x_{1}^{k},x_{2}^{k},\cdots
,x_{n}^{k})}{\partial x_{n}} \\ 
\vdots  & \vdots  & \ddots  & \vdots  \\ 
\dfrac{\partial f_{n}(x_{1}^{k},x_{2}^{k},\cdots ,x_{n}^{k})}{\partial x_{1}}
& \dfrac{\partial f_{n}(x_{1}^{k},x_{2}^{k},\cdots ,x_{n}^{k})}{\partial
x_{2}} & \cdots  & \dfrac{\partial f_{n}(x_{1}^{k},x_{2}^{k},\cdots
,x_{n}^{k})}{\partial x_{n}}
\end{array}
\right] ,\Delta x=\left[ 
\begin{array}{c}
\Delta x_{1}^{k} \\ 
\Delta x_{2}^{k} \\ 
\vdots  \\ 
\Delta x_{n}^{k}
\end{array}
\right]   \tag{12)(13}
\end{equation}

Eq. (11) is Newton's method for an n x n system (a remarkably compact and
useful equation).

\vspace{0.5cm}

We now consider three small nonlinear problems:

\vspace{0.5cm}

(1) The third order polynomial

\[
f(x)=x^{3}-2x-5=0 
\]

As noted on 28jan00, this polynomial may be the first application of
Newton's method, since Newton used it to illustrate his method

\vspace{0.5cm}

If we apply Newton's method, eq. (9), for the scalar ($n=1$) case

\[
\frac{df(x^{k})}{dx}(x^{k+1}-x^{k})=-f(x^{k}) 
\]

or the familiar form

\[
x^{k+1}=x^{k}-\frac{f(x^{k})}{df(x^{k})/dx},k=0,1,\cdots 
\]

where $k$ is an iteration index or counter.

\vspace{0.5cm}

We can make some observations about the application of Newton's method to
Newton's polynomial and the generalization of Newton's method to the n x n
problem (eq. (11)):

\vspace{0.5cm}

\begin{itemize}
\item  The polynomial has three roots (three reals or one real and two
complex conjugates - conjugates are required for the polynomial to have real
coefficients). This conclusion is from a famous proof by C. F. Gauss; in
general, \textit{nonlinear problems will have multiple (nonunique) roots.}

\vspace{0.5cm}

\item  If we make the choice of the initial starting point at a root of

\[
\frac{df(x)}{dx}=3x^{2}-2=0
\]

or

\[
x^{0}=\pm \sqrt{2/3}
\]

the method will fail because this is a \textit{singular point} (the \textit{%
Jacobian matrix}, or in this case, the first derivative, is singular).

\vspace{0.5cm}

\item  Also, if we choose a starting point near the singular point, the
method will most likely fail, i.e., the system is \textit{near-singular}
(this is clear if we consider Newton's method as a projection along the
tangent to $f(x)$).

\vspace{0.5cm}

\item  The choice of an initial point that is not singular or near-singular
becomes increasingly difficult as the size of the problem increases (more
equations), and we should have a test in our software for the \textit{\
condition} of the Jacobian matrix.

\vspace{0.5cm}

\item  The Jacobian matrix, \textbf{$\mathbf{J}$} (in eq. (11)), is n x n,
and therfore grows quickly with the problem size.

\vspace{0.5cm}

\item  The Jacobian matrix is not constant (a property of nonlinear
systems), and therefore must be updated as the iterations proceed; updating
at every iteration may not be required since all that is required is
convergence of the iterations.

\vspace{0.5cm}

\item  The right hand side (RHS) vector, $\mathbf{f}$ in eq. (11), must also
be updated as the iterations proceed, but this will be less costly than for
the Jacobian matrix because of the smaller size of $\mathbf{f}$.

\vspace{0.5cm}

\item  A stopping criterion is required for the iterations, e.g.,

\[
\left| \Delta x_{i}^{k}\right| <\varepsilon _{i},i=1,2,\cdots ,n\text{
(relative or absolute)}
\]
\[
\left| \Delta x_{i}^{k+1}-\Delta x_{i}^{k}\right| <\varepsilon
_{i},i=1,2,\cdots ,n\text{ (relative or absolute)}
\]

\[
\left| f_{i}(\mathbf{x}_{i}^{k})\right| <\varepsilon _{i},i=1,2,\cdots ,n
\]
\[
\left| f_{i}(\mathbf{x}_{i}^{k+1})-f_{i}(\mathbf{x}_{i}^{k})\right|
<\varepsilon _{i},i=1,2,\cdots ,n
\]

\vspace{0.5cm}

\item  Quasilinearization is the characterization of eq. (11) denoting the
linearization with respect to the iteration index or counter, $k$.

\vspace{0.5cm}

\item  Eq. (11) is almost never solved as

\[
\mathbf{\Delta x}=-J^{-1}\mathbf{f}
\]

Rather, $J$ is factored (LU, QR factorization), and eq. (11) is then solved
as a system of linear algebraic equations (for the vector of Newton
corrections \textbf{$\Delta x$}$\mathbf{)}$.

\vspace{0.5cm}

\item  For large nonlinear systems, e.g., $n>100$, the structure of the
Jacobian matrix in eq. (11) can often be exploited to good advantage, e.g.,
as a banded or sparse matrix.
\end{itemize}

\vspace{0.5cm}

(2) A 2 x 2 system provided by L. T. Biegler (1991):

\begin{equation}
\begin{array}{c}
f_{1}(x_{1},x_{2})=x_{1}^{2}+x_{2}^{2}-17 \\ 
f_{2}(x_{1},x_{2})=2x_{1}^{1/3}+x_{2}^{1/2}-4
\end{array}
\tag{14}
\end{equation}

Eqs. (14) have the solutions $\mathbf{x}=(1,4)^{T}$ and $\mathbf{x}%
=(4.07150,0.65027)^{T}$.

\vspace{0.5cm}

Note the nonlinearity of this system. This is manifest in the Jacobian
matrix,

\begin{equation}
\mathbf{J}=\left[ 
\begin{array}{cc}
\partial f_{1}/\partial x_{1}=2x_{1} & \partial f_{1}/\partial x_{2}=2x_{2}
\\ 
\partial f_{2}/\partial x_{1}=(2/3)x_{1}^{-2/3} & \partial f_{2}/\partial
x_{2}=(1/2)x_{2}^{-1/2}
\end{array}
\right]  \tag{15}
\end{equation}

since the elements of $J$ are functions of $x_{1},x_{2}$ (for a linear
system, the Jacobian matrix is a constant matrix).

\vspace{0.5cm}

$J$ is singular if

\[
\left| \mathbf{J}\right| =\left| 
\begin{array}{cc}
2x_{1} & 2x_{2} \\ 
(2/3)x_{1}^{-2/3} & (1/2)x_{2}^{-1/2}
\end{array}
\right| =0 
\]

or

\begin{equation}
\left( 2x_{1}\right) (1/2)x_{2}^{-1/2}-\left( 2x_{2}\right)
(2/3)x_{1}^{-2/3}=0  \tag{16}
\end{equation}

Eq. (16) is a relationship between $x_{1}$ and $x_{2}$ that defines a locus
of singular points for system (14). Later we will start a numerical solution
at a singular point given by eq. (16) and observe what happens with the
calculation of the solution.

\vspace{0.5cm}

Applying eq. (11) to eqs. (14), we have

\[
\left[ 
\begin{array}{cc}
2x_{1} & 2x_{2} \\ 
(2/3)x_{1}^{-2/3} & (1/2)x_{2}^{-1/2}
\end{array}
\right] ^{k}\left[ 
\begin{array}{c}
\Delta x_{1} \\ 
\Delta x_{2}
\end{array}
\right] ^{k}=-\left[ 
\begin{array}{c}
x_{1}^{2}+x_{2}^{2}-17 \\ 
2x_{1}^{1/3}+x_{2}^{1/2}-4
\end{array}
\right] ^{k} 
\]

or

\begin{equation}
\begin{array}{c}
2x_{1}^{k}\Delta x_{1}^{k}+2x_{2}^{k}\Delta
x_{2}^{k}=-(x_{1}^{k})^{2}-(x_{2}^{k})^{2}+17 \\ 
(2/3)(x_{1}^{k})^{-2/3}\Delta x_{1}^{k}+(1/2)(x_{2}^{k})^{-1/2}\Delta
x_{2}^{k}=-2(x_{1}^{k})^{1/3}-(x_{2}^{k})^{1/2}+4
\end{array}
\tag{17}
\end{equation}

A basic algorithm for solving eqs. (14) is therefore

\[
\begin{array}{cc}
1 & \text{Select }x_{1}^{k}=x_{1}^{0},x_{2}^{k}=x_{2}^{0} \\ 
2 & \text{Substitute }x_{1}^{k},x_{2}^{k}\text{ in eqs. (17) } \\ 
3 & \text{Solve eqs. (17) for }\Delta x_{1}^{k},\Delta x_{2}^{k} \\ 
4 & \left| \Delta x_{1}^{k}\right| <\varepsilon _{1}\text{ and }\left|
\Delta x_{2}^{k}\right| <\varepsilon _{2}\text{ ?; yes - output, quit; no -
continue} \\ 
5 & x_{1}^{k+1}=x_{1}^{k}+\Delta x_{1}^{k},x_{2}^{k+1}=x_{2}^{k}+\Delta
x_{2}^{k} \\ 
6 & \text{Let }x_{1}^{k}=x_{1}^{k+1},x_{2}^{k}=x_{2}^{k+1};\text{go to 2}
\end{array}
\]

Note that all that is required to implement this algorithm is a linear
equation solver (to compute the Newton corrections, $\Delta x_{1}^{k},
\Delta x_{2}^{k}$).

\vspace{0.5cm}

(3) A 4 x 4 system from Riggs, J. B. (1994), ``An Introduction to Numerical
Methods for Chemical Engineers'', 2nd ed., Texas Tech University Press, pp
82-89:

\vspace{0.5cm}

A continuous stirred tank reactor (CSTR) is used for the following reaction
scheme

\[
A\stackrel{r_{1}}{\rightarrow }2B 
\]
\[
A\stackrel{r_{2}}{\stackunder{r_{3}}{\rightleftarrows }}C 
\]
\[
B\stackrel{r_{4}}{\rightarrow }C+D 
\]

where

\[
\begin{array}{c}
r_{1}=k_{1}C_{A} \\ 
r_{2}=k_{2}C_{A}^{3/2} \\ 
r_{3}=k_{3}C_{C}^{2} \\ 
r_{4}=k_{4}C_{B}^{2} \\ 
k_{1}=1.0\text{ sec}^{-1} \\ 
k_{2}=0.2\text{ liter}^{1/2}\text{/gm mol}^{1/2}\text{-sec} \\ 
k_{3}=0.05\text{ liter/gm mol-sec} \\ 
k_{4}=0.4\text{ liter/gm mol-sec}
\end{array}
\]

\vspace{0.5cm}

Each reaction rate ($r_{1}$ to $r_{4}$) is in gm mols/liter-sec. If the
volume of the reactor is $V_{r}$ $=100$ liters, the volumetric flow rate
through the reactor is $Q=50$ liters/sec, and the feed to the reactor is $%
C_{A0}=1.0$ gm mol/liter, the steady state material balances for the reactor
are:

\begin{equation}
f_{1}(C_{A},C_{B},C_{C},C_{D})=-QC_{A}+QC_{A0}+V_{R}(-k_{1}C_{A}-k_{2}C_{A}^{3/2}+k_{3}C_{C}^{2})=0
\tag{17a}
\end{equation}

\begin{equation}
f_{2}(C_{A},C_{B},C_{C},C_{D})=-QC_{B}+V_{R}(2k_{1}C_{A}-k_{4}C_{B}^{2})=0 
\tag{17b}
\end{equation}

\begin{equation}
f_{3}(C_{A},C_{B},C_{C},C_{D})=-QC_{C}+V_{R}(k_{2}C_{A}^{3/2}-k_{3}C_{C}^{2}+k_{4}C_{B}^{2})=0
\tag{17c}
\end{equation}

\begin{equation}
f_{4}(C_{A},C_{B},C_{C},C_{D})=-QC_{D}+V_{R}(k_{4}C_{B}^{2})=0  \tag{17d}
\end{equation}

We will subsequently consider a solution to eqs. (17) to three significant
figures, i.e., the values of $C_{A},C_{B},C_{C},C_{D}$ that satisfy eqs.
(17).

\vspace{0.5cm}

One possible approach to computing a solution would be to convert eqs. (17)
to a system of ordinary differential equations (ODEs):

\begin{equation}
V\frac{dC_{A}}{dt}%
=-QC_{A}+QC_{A0}+V_{R}(-k_{1}C_{A}-k_{2}C_{A}^{3/2}+k_{3}C_{C}^{2})=0,C_{A}(0)=C_{A}^{0}
\tag{18a}
\end{equation}

\begin{equation}
V\frac{dC_{B}}{dt}%
=-QC_{B}+V_{R}(2k_{1}C_{A}-k_{4}C_{B}^{2})=0,C_{B}(0)=C_{B}^{0}  \tag{18b}
\end{equation}

\begin{equation}
V\frac{dC_{C}}{dt}%
=-QC_{C}+V_{R}(k_{2}C_{A}^{3/2}-k_{3}C_{C}^{2}+k_{4}C_{B}^{2})=0,C_{C}(0)=C_{C}^{0}
\tag{18c}
\end{equation}

\begin{equation}
V\frac{dC_{D}}{dt}=-QC_{D}+V_{R}(k_{4}C_{B}^{2})=0,C_{D}(0)=C_{D}^{0} 
\tag{18d}
\end{equation}

where: (a) $V$ is the volume of the reactor (if we are interested in the
true dynamic response of the CSTR) or (b) $V$ is a ``pseudo'' time constant
if all we are interested in is the final steady state with $\dfrac{dC_{A}}{dt%
}\approx \dfrac{dC_{B}}{dt}\approx \dfrac{dC_{C}}{dt}\approx \dfrac{dC_{D}}{%
dt}\approx 0$.

\vspace{0.5cm}

Method (b) is a perfectly satisfactory way of solving eqs. (14) called the
method of ``false transients''. However, care must be given to adding the
derivatives to the nonlinear equations with the correct sign to ensure a
stable steady state or equilibrium solution (which may not be obvious for
large systems of equations). Later we will consider a systematic procedure
for adding the derivatives.

\newpage

\begin{center}
{\Huge \textbf{Continuation Methods}}
\end{center}

\vspace{0.5cm}

Consider again the 2 x 2 system of eq. (7). If we assume two related
functions, 
\[
\begin{array}{c}
g_{1}(x_{1},x_{2}) \\ 
g_{2}(x_{1},x_{2})
\end{array}
\]

with known solutions based on arbitrary values $%
x_{1}=x_{1}^{0},x_{2}=x_{2}^{0},$i.e., 
\[
\begin{array}{c}
g_{1}(x_{1}^{0},x_{2}^{0})=0 \\ 
g_{2}(x_{1}^{0},x_{2}^{0})=0
\end{array}
\]

then we can form a vector of \textit{homotopy functions} 
\begin{equation}
\begin{array}{c}
h_{1}(x_{1},x_{2},t)=tf_{1}(x_{1},x_{2})+(1-t)g_{1}(x_{1},x_{2}) \\ 
h_{2}(x_{1},x_{2},t)=tf_{2}(x_{1},x_{2})+(1-t)g_{2}(x_{1},x_{2})
\end{array}
\tag{19}
\end{equation}

where $t$ is an \textit{embedded homotopy continuation parameter}.

\vspace{0.5cm}

Now we vary $t$ over the interval $0\leq t\leq 1$ while always setting the
homotopy functions to zero. This can be easily visualized using a table

\[
\begin{array}{cccc}
t & h_{1}(x_{1},x_{2},t)=0 & x_{1,}x_{2} & \text{Comments} \\ 
& h_{2}(x_{1},x_{2},t)=0 &  &  \\ 
&  &  &  \\ 
0 & g_{1}(x_{1},x_{2})=0 & x_{1}^{0} & \text{Assumed initial} \\ 
& g_{2}(x_{1},x_{2})=0 & x_{2}^{0} & \text{solution} \\ 
&  &  &  \\ 
0.1 & 0.1f_{1}(x_{1},x_{2})+0.9g_{1}(x_{1},x_{2})=0 & x_{1}^{0.1} & 
x_{1}^{0.1}\text{ close to }x_{1}^{0} \\ 
& 0.1f_{2}(x_{1},x_{2})+0.9g_{2}(x_{1},x_{2})=0 & x_{2}^{0.1} & x_{2}^{0.1}%
\text{ close to }x_{2}^{0} \\ 
&  &  &  \\ 
0.2 & 0.2f_{1}(x_{1},x_{2})+0.8g_{1}(x_{1},x_{2})=0 & x_{1}^{0.2} & \text{
Root finder required} \\ 
& 0.2f_{2}(x_{1},x_{2})+0.8g_{2}(x_{1},x_{2})=0 & x_{2}^{0.2} & \text{ at
each step} \\ 
&  &  &  \\ 
0.3 & 0.3f_{1}(x_{1},x_{2})+0.7g_{1}(x_{1},x_{2})=0 & x_{1}^{0.3} &  \\ 
& 0.3f_{2}(x_{1},x_{2})+0.7g_{2}(x_{1},x_{2})=0 & x_{2}^{0.3} &  \\ 
&  &  &  \\ 
\vdots & \vdots & \vdots &  \\ 
&  &  &  \\ 
0.9 & 0.9f_{1}(x_{1},x_{2})+0.9g_{1}(x_{1},x_{2})=0 & x_{1}^{0.9} &  \\ 
& 0.9f_{2}(x_{1},x_{2})+0.9g_{2}(x_{1},x_{2})=0 & x_{2}^{0.9} &  \\ 
&  &  &  \\ 
1 & f_{1}(x_{1},x_{2})=0 & x_{1}^{1} & \text{Required solution} \\ 
& f_{2}(x_{1},x_{2})=0 & x_{2}^{1} & 
\end{array}
\]

\vspace{0.5cm}

The essential features of this continuation method are:

\vspace{0.5cm}

\begin{itemize}
\item  The solution is \textit{continued} from the known solution

\vspace{0.5cm}

\item  The closer $g_{1}(x_{1},x_{2})$ and $g_{1}(x_{1},x_{2})$ are to $%
f_{1}(x_{1},x_{2})$ and $f_{2}(x_{1},x_{2})$, the better chances of this
procedure working.

\vspace{0.5cm}

\item  However, the success of the method is generally not particularly
sensitive to the selection of $g_{1}(x_{1},x_{2})$ and $g_{1}(x_{1},x_{2})$.
For example, we could try linear functions

\[
\begin{array}{c}
g_{1}(x_{1},x_{2})=(x_{1}-x_{1}^{0})+(x_{2}-x_{2}^{0}) \\ 
g_{2}(x_{1},x_{2})=(x_{1}-x_{1}^{0})+(x_{2}-x_{2}^{0})
\end{array}
\]

\vspace{0.5cm}

\item  A root finder is required at each step (each value of $t$).

\vspace{0.5cm}

\item  However, the initial guess for the root finder can be the solution of
the preceding step (since the change in the solution will be small from one
step to the next).

\vspace{0.5cm}

\item  If the method fails, the increment in $t$ can be reduced and the
calculation repeated (the calculations in the preceding table are easily
programmed in a loop, and the loop can be be repeated more often to improve
the successful continuation from the assumed initial condition to the final
solution).

\vspace{0.5cm}

\item  The extension to the n x n problem is straightforward.
\end{itemize}

\vspace{0.5cm}

The homotopy functions of eq. (19) are just one possibility with one
particular form of embedding of the continuation parameter $t$. Here are
some other possibilities:

\vspace{0.5cm}

(1)

\begin{equation}
\begin{array}{c}
h_{1}(x_{1},x_{2,}t)=f_{1}(x_{1},x_{2})-(1-t)f_{1}(x_{1}^{0},x_{2}^{0}) \\ 
h_{2}(x_{1},x_{2},t)=f_{2}(x_{1},x_{2})-(1-t)f_{2}(x_{1}^{0},x_{2}^{0})
\end{array}
\tag{20}
\end{equation}

with the properties

\[
\begin{array}{c}
h_{1}(x_{1},x_{2,}0)=0 \\ 
h_{2}(x_{1},x_{2,}0)=0
\end{array}
\]

\[
\begin{array}{c}
h_{1}(x_{1},x_{2,}1)=f_{1}(x_{1},x_{2}) \\ 
h_{2}(x_{1},x_{2},1)=f_{2}(x_{1},x_{2})
\end{array}
\]

as required.

\vspace{0.5cm}

(2)

\begin{equation}
\begin{array}{c}
h_{1}(x_{1},x_{2,}t)=f_{1}(x_{1},x_{2})-e^{-t}f_{1}(x_{1}^{0},x_{2}^{0}) \\ 
h_{2}(x_{1},x_{2},t)=f_{2}(x_{1},x_{2})-e^{-t}f_{2}(x_{1}^{0},x_{2}^{0})
\end{array}
\tag{21}
\end{equation}

with the properties

\[
\begin{array}{c}
h_{1}(x_{1},x_{2,}0)=0 \\ 
h_{2}(x_{1},x_{2,}0)=0
\end{array}
\]

\[
\begin{array}{c}
h_{1}(x_{1},x_{2,}\infty )=f_{1}(x_{1},x_{2}) \\ 
h_{2}(x_{1},x_{2},\infty )=f_{2}(x_{1},x_{2})
\end{array}
\]

as required (note that the interval in $t$ is now $0\leq t\leq \infty $).

\newpage

\begin{center}
{\Huge \textbf{Davidenko's Equation}}
\end{center}

\vspace{0.5cm}

If we consider again the embedding of eqs. (21), we can analyze it by
considering the differentials of the homotopy functions

\begin{equation}
\begin{array}{c}
dh_{1}=\dfrac{\partial h_{1}}{\partial x_{1}}dx_{1}+\dfrac{\partial h_{1}}{%
\partial x_{2}}dx_{2}+\dfrac{\partial h_{1}}{\partial t}dt=0 \\ 
dh_{2}=\dfrac{\partial h_{2}}{\partial x_{1}}dx_{1}+\dfrac{\partial h_{2}}{%
\partial x_{2}}dx_{2}+\dfrac{\partial h_{2}}{\partial t}dt=0
\end{array}
\tag{22}
\end{equation}

(these differentials are zero since $h_{1}(x_{1},x_{2},t)$ $%
=0,h_{2}(x_{1},x_{2},t)$ $=0)$. From eqs. (21) and (22), we have

\[
\begin{array}{c}
dh_{1}=\dfrac{\partial f_{1}}{\partial x_{1}}dx_{1}+\dfrac{\partial f_{1}}{%
\partial x_{2}}dx_{2}+e^{-t}f_{1}(x_{1}^{0},x_{2}^{0})dt=0 \\ 
dh_{2}=\dfrac{\partial f_{2}}{\partial x_{1}}dx_{1}+\dfrac{\partial f_{2}}{%
\partial x_{2}}dx_{2}+e^{-t}f_{2}(x_{1}^{0},x_{2}^{0})dt=0
\end{array}
\]

or

\begin{equation}
\begin{array}{c}
\dfrac{\partial f_{1}}{\partial x_{1}}\dfrac{dx_{1}}{dt}+\dfrac{\partial
f_{1}}{\partial x_{2}}\dfrac{dx_{2}}{dt}%
=-e^{-t}f_{1}(x_{1}^{0},x_{2}^{0})=-f_{1}(x_{1},x_{2}) \\ 
\dfrac{\partial f_{2}}{\partial x_{1}}\dfrac{dx_{1}}{dt}+\dfrac{\partial
f_{2}}{\partial x_{2}}\dfrac{dx_{2}}{dt}%
=-e^{-t}f_{2}(x_{1}^{0},x_{2}^{0})=-f_{2}(x_{1},x_{2})
\end{array}
\tag{23}
\end{equation}

Eqs. (23), a differential form of Newton's method, are the \textit{Davidenko
differential equation}, that can be written in matrix form for the n x n
problem as

\begin{equation}
J\frac{d\mathbf{x}}{dt}=-\mathbf{f}  \tag{24}
\end{equation}

where again $J$ and $\mathbf{f}$\textbf{\ }are the Jacobian matrix and
function vector of the nonlinear system, respectively, and $\dfrac{d\mathbf{x%
}}{dt}$ is the derivative of the solution vector with respect to the
continuation parameter, $t$ (compare eq. (11) with eq. (24) to see why
Davidenko's ODE is a differential form of Newton's method).

\vspace{0.5cm}

Eq. (24) defines an initial value problem, requiring an initial condition
vector that serves as the assumed initial or starting solution

\begin{equation}
\mathbf{x}(0)=\mathbf{x}_{0}  \tag{25}
\end{equation}

The procedure for the solution of the nonlinear system is to integrate eq.
(24) until

\begin{equation}
J\frac{d\mathbf{x}}{dt}\approx 0\text{ which implies }\mathbf{f\approx 0} 
\tag{26}
\end{equation}

This is the required solution of the nonlinear system, i.e., whatever the
solution vector, $\mathbf{x}(t),$ of eq. (24) is when condition (26) occurs
is the required solution of the original nonlinear system.

\vspace{0.5cm}

Note that form of ODEs (23); there is more than one derivative in each ODE,
and they are termed \textit{linearly implicit ODEs}. Special integrators are
available for ODEs of this form, e.g., LSODI, DASSL.

\vspace{0.5cm}

Since eqs. (23) are linear in the derivatives, another approach to
integrating them is first use a linear algebraic equation solver to compute
the derivative vector, $\dfrac{d\mathbf{x}}{dt},$ then use a conventional
ODE integrator for the \textit{explicit ODEs} to compute a solution. In
other words, eq. (24) can be written as

\begin{equation}
\frac{d\mathbf{x}}{dt}=-J^{-1}\mathbf{f}  \tag{27}
\end{equation}

This method works well for relatively low order systems.

\vspace{0.5cm}

The advantages of using eq. (24) or eq. (27) relative to the conventional
finite stepping Newton method, eq. (11), are

\begin{itemize}
\item  The solution of eqs. (24) or (27) proceeds by differential steps
rather than the finite steps of Newton's method. Thus, the convergence to a
solution is often more reliable, i.e., Davidenko's method may work when
Newton's method fails.

\vspace{0.5cm}

\item  In other words, we are relying on the error monitoring and step size
control of an ODE integrator (rather than the use of damping in Newton's
method).

\vspace{0.5cm}

\item  The implementation of eq. (24) requires only an integrator for
linearly implicit ODEs, e.g., LSODI, DASSL.

\vspace{0.5cm}

\item  Eq. (27) requires only a linear equation solver plus an integrator
for explicit ODEs.
\end{itemize}

Thus, the differential approach to Newton's method can be implemented with
readily available software.

\vspace{0.5cm}

To conclude this section, consider the application of eqs. (24) and (27) to
problem (2). For eq. (24)

\[
\left[ 
\begin{array}{cc}
2x_{1} & 2x_{2} \\ 
(2/3)x_{1}^{-2/3} & (1/2)x_{2}^{-1/2}
\end{array}
\right] \left[ 
\begin{array}{c}
\dfrac{dx_{1}}{dt} \\ 
\\ 
\dfrac{dx_{2}}{dt}
\end{array}
\right] =-\left[ 
\begin{array}{c}
x_{1}^{2}+x_{2}^{2}-17 \\ 
2x_{1}^{1/3}+x_{2}^{1/2}-4
\end{array}
\right] 
\]

or

\begin{equation}
\begin{array}{c}
2x_{1}\dfrac{dx_{1}}{dt}+2x_{2}\dfrac{dx_{2}}{dt}=-\left(
x_{1}^{2}+x_{2}^{2}-17\right) \\ 
(2/3)x_{1}^{-2/3}\dfrac{dx_{1}}{dt}+(1/2)x_{2}^{-1/2}\dfrac{dx_{2}}{dt}%
=-\left( 2x_{1}^{1/3}+x_{2}^{1/2}-4\right)
\end{array}
\tag{28}
\end{equation}

For eq. (27) 
\begin{equation}
\left[ 
\begin{array}{c}
\dfrac{dx_{1}}{dt} \\ 
\\ 
\dfrac{dx_{2}}{dt}
\end{array}
\right] =-\left[ 
\begin{array}{cc}
2x_{1} & 2x_{2} \\ 
(2/3)x_{1}^{-2/3} & (1/2)x_{2}^{-1/2}
\end{array}
\right] ^{-1}\left[ 
\begin{array}{c}
x_{1}^{2}+x_{2}^{2}-17 \\ 
2x_{1}^{1/3}+x_{2}^{1/2}-4
\end{array}
\right]  \tag{29}
\end{equation}

A Fortran program for the solution of eqs. (29) is given in Schiesser
(1994), pp 174-193.

\vspace{0.5cm}

Finally, we can analyze the convergence of solutions to Davidenko's ODE, eq.
(24). For the scalar case, eq. (24) becomes

\[
\frac{df}{dx}\frac{dx}{dt}=-f 
\]

which rearranges to

\[
\frac{df}{f}=-dt 
\]

Integration of each side of this ODE, subject to the initial condition

\[
f(t_{0})=f_{0} 
\]

gives

\[
\ln (f/f_{0})=-(t-t_{0}) 
\]

or

\begin{equation}
f=f_{0}e^{-(t-t_{0})}  \tag{30}
\end{equation}

Eq. (30) indicates that the convergence of Davidenko's method is exponential
in $t$. This result also carries over to the n x n case. This exponential
convergence is demonstrated numerically for problem (2) in Schiesser (1994).

\newpage

\begin{center}
{\Huge \textbf{Levenberg-Marquardt Method}}
\end{center}

\vspace{0.5cm}

A major difficulty with Newton's method, eq. (11), or Davidenko's method,
eq. (24), is failure caused by a singular or near singular Jacobian matrix
(note that both methods require the inverse Jacobian matrix). To circumvent
this problem, we consider the discrete and differential (continuous)
Levenberg-Marquardt methods:

\begin{equation}
\left\{ (1-\mu )J^{T}J+\mu I\right\} \Delta x=-J^{T}\mathbf{f}  \tag{31}
\end{equation}

\begin{equation}
\left\{ (1-\mu )J^{T}J+\mu I\right\} \frac{dx}{dt}=-J^{T}\mathbf{f}  \tag{32}
\end{equation}

When $\mu =0,$eqs. (31) and (32) reduce to

\begin{equation}
J^{T}J\Delta x=-J^{T}\mathbf{f}\text{ or }J\Delta x=-\mathbf{f}  \tag{33}
\end{equation}

\begin{equation}
J^{T}J\frac{dx}{dt}=-J^{T}\mathbf{f}\text{ or }J\frac{dx}{dt}=-\mathbf{f} 
\tag{34}
\end{equation}

Eqs. (33) and (34) are just Newton's method and Davidenko's method,
respectively (and both will fail if $J$ is singular).

When $\mu =1,$eqs. (31) and (32) reduce to 
\begin{equation}
\Delta x=-J^{T}\mathbf{f}  \tag{35}
\end{equation}

\begin{equation}
\frac{dx}{dt}=-J^{T}\mathbf{f}  \tag{36}
\end{equation}

Eqs. (35) and (36) express the \textit{discrete and differential steepest
descent}. Note that neither eq. (35) nor eq. (36) requires the inverse
Jacobian. However, these equations reflect the conditions for a maximum or
minimum in

\begin{equation}
ss=\sum_{i}^{N}f_{i}^{2}  \tag{37}
\end{equation}

Thus, eqs. (35) and (36) can ``get stuck'' at a local maximum or minimum of $%
ss$ of eq. (37). Also, when they do approach the solution $%
f_{1}=f_{2}=\cdots f_{N}=0$, they do so relatively slowly (compared to
Newton's method).

\vspace{0.5cm}

Therefore, the approach to avoiding a singular system and still approach the
solution $f_{1}=f_{2}=\cdots f_{N}=0,$ is to use $0<\mu <1$, and to possibly
vary $\mu $ as the solution proceeds by monitoring the condition of $J.$

\vspace{0.5cm}

Equations (31) and (32) have been coded in library software. In particular,
eq. (32) is available as a Fortran program from WES. We conclude this
discussion with the application of eq. (32) to problem (3).

\vspace{0.5cm}

The initial conditions for eq. (32) are
\begin{verbatim}
 
C...
C...  First starting point
      IF(NORUN.EQ.1)THEN
         X(1)=0.5D0
         X(2)=0.5D0
         X(3)=0.5D0
         X(4)=0.5D0
      END IF
C...
C...  LM parameter in eq. (1).  Note: the calculation fails with
C...  L = 0 for a singular system (determinant of the Jacobian
C...  matrix is zero)
      L=0.5D0
 
\end{verbatim}

The programming of the RHS of eqs. (32) is then
\begin{verbatim}
 
C...
C...  Once past the singular point, the DLM parameter is reset to
C...  zero to give the differential Newton's method (and thereby
C...  enhance the rate of convergence)
      IF(T.GT.1.0D0)L=0.0D0
C...
C...  Vector of functions (f in eq. (1))
      CALL FUNC(N,X,F)
C...
C...  Jacobian matrix (J in eq. (1))
      CALL JAC(N,X,J)
C...                              T
C...  Jacobian matrix transpose (J in eq. (1))
      CALL TRANS(N,J,JT)
C...
C...                   T
C...  Matrix product (J*J in eq. (1))
      CALL MULT(N,N,N,J,JT,JTJ)
C...                                        T
C...  Left hand side coupling matrix (1-L)*J*J + L*I in eq. (1))
      CALL CMAT(N,JTJ,L,ID,JTJPLI)
C...
C...                            T
C...  Right hand side vector (-J*f in eq. (1))
      CALL MULT1(N,N,1,JT,F,JTF)
C...
C...  Solve for the vector of derivatives, dx/dt, in eq. (1) (in
C...  effect, eq. (1) is premultiplied by the inverse of the LHS
C...  coupling matrix).  DECOMP performs an LU decomposition and
C...  computes the condition number of JTJPLI
      CALL DECOMP(N,N,JTJPLI,COND,IPVT,WORK)
C...
C...  Compute derivative vector dx/dt from JTJPLI with RHS JTF,
C...  i.e., solve eq. (1) for dx/dt.  The solution is returned
C...  in JTF
      CALL SOLVE(N,N,JTJPLI,JTF,IPVT)
C...
C...  Transfer array JTF to array DXDT so the derivative vector,
C...  dx/dt, is available in COMMON/F/ for the ODE integrator
      DO 2 I=1,N
         DXDT(I)=JTF(I,1)
2     CONTINUE
      RETURN
      END
 
\end{verbatim}

The programming of the function vector is
\begin{verbatim}
 
      SUBROUTINE FUNC(N,X,F)
C...
C...  Subroutine FUNC computes the vector of functions,
C...  -f, in eq. (1) (note the minus sign, as required
C...  in the RHS of eq. (1))
C...
C...  Arguments
C...
C...     N      order of the NLE system (number of NLE)
C...            (input)
C...
C...     X      solution vector of size N (input)
C...
C...     F      function vector of size N (output)
C...
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION X(N), F(N)
C...
C...  COMMON area for the problem parameters
      COMMON/P/ VR, Q, CA0, RK1, RK2, RK3, RK4
C...
C...  Problem variables
      CA=X(1)
      CB=X(2)
      CC=X(3)
      CD=X(4)
C...
C...  Function vector
      F(1)=-CA+VR*(-RK1*CA-RK2*CA**1.5D0+RK3*CC**2)/Q+CA0
      F(2)=-CB+VR*(2.0D0*RK1*CA-RK4*CB**2)/Q
      F(3)=-CC+VR*(RK2*CA**1.5D0-RK3*CC**2+RK4*CB**2)/Q
      F(4)=-CD+VR*(RK4*CB**2)/Q
C...
C...  Reverse the sign of the function vector in accordance
C...  with the RHS of eq. (1)
      DO 1 I=1,N
         F(I)=-F(I)
 1    CONTINUE
      RETURN
      END
 
\end{verbatim}

The programming of the Jacobian matrix is:
\begin{verbatim}
 
      SUBROUTINE JAC(N,X,J)
C...
C...  Subroutine JAC computes the Jacobian matrix
C...  of the vector of functions, f, in eq. (1)
C...  (note the elements of J are computed for f,
C...  not -f)
C...
C...  Arguments
C...
C...     N      order of the NLE system (number of NLE)
C...            (input)
C...
C...     X      solution vector of size N (input)
C...
C...     J      Jacobian matrix of size N x N (output)
C...
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION J
      DIMENSION X(N), J(N,N)
C...
C...  COMMON area for the problem parameters
      COMMON/P/ VR, Q, CA0, RK1, RK2, RK3, RK4
C...
C...  Problem variables
      CA=X(1)
      CB=X(2)
      CC=X(3)
      CD=X(4)
C...
C...  Jacobian matrix
      J(1,1)=-1.0D0+VR*(-RK1-1.5D0*RK2*CA**0.5)/Q
      J(1,2)=0.0D0
      J(1,3)=VR*(2.0D0*RK3*CC)/Q
      J(1,4)=0.0D0
      J(2,1)=VR*(2.0D0*RK1)/Q
      J(2,2)=-1.0D0+VR*(-2.0D0*RK4*CB)/Q
      J(2,3)=0.0D0
      J(2,4)=0.0D0
      J(3,1)=VR*(1.5D0*RK2*CA**0.5D0)/Q
      J(3,2)=VR*(2.0D0*RK4*CB)/Q
      J(3,3)=-1.0D0+VR*(-2.0D0*RK3*CC)/Q
      J(3,4)=0.0D0
      J(4,1)=0.0D0
      J(4,2)=VR*(2.0D0*RK4*CB)/Q
      J(4,3)=0.0D0
      J(4,4)=-1.0D0
      RETURN
      END
 
\end{verbatim}

The output from the complete program is:
\begin{verbatim}
 
 RUN NO. -   1  Four component, steady state CSTR
                equations
 
 INITIAL T -  0.000D+00
 
   FINAL T -  0.200D+03
 
   PRINT T -  0.500D+02
 
 NUMBER OF DIFFERENTIAL EQUATIONS - 4
 
 INTEGRATION ALGORITHM - LSODES
 
 MAXIMUM INTEGRATION ERROR -  0.100D-06
 
1
  t =    0.0
 Ca = 0.50000  dCa/dt =  -0.416D+01  res f1 =   0.616D+00
 Cb = 0.50000  dCb/dt =  -0.261D+01  res f2 =  -0.130D+01
 Cc = 0.50000  dCc/dt =  -0.176D+01  res f3 =   0.184D+00
 Cd = 0.50000  dCd/dt =  -0.122D+01  res f4 =   0.300D+00
 
  t =   50.0
 Ca = 0.31814  dCa/dt =   0.120D-03  res f1 =  -0.200D-02
 Cb = 0.78160  dCb/dt =   0.249D-03  res f2 =  -0.222D-02
 Cc = 0.53095  dCc/dt =   0.405D-03  res f3 =  -0.136D-02
 Cd = 0.48780  dCd/dt =   0.375D-03  res f4 =  -0.926D-03
 
  t =  100.0
 Ca = 0.31886  dCa/dt =  -0.899D-04  res f1 =  -0.992D-05
 Cb = 0.78387  dCb/dt =  -0.648D-04  res f2 =  -0.211D-04
 Cc = 0.53495  dCc/dt =  -0.444D-04  res f3 =  -0.128D-04
 Cd = 0.49155  dCd/dt =  -0.410D-04  res f4 =  -0.812D-05
 
  t =  150.0
 Ca = 0.31887  dCa/dt =  -0.118D-05  res f1 =   0.857D-06
 Cb = 0.78388  dCb/dt =  -0.118D-05  res f2 =   0.887D-06
 Cc = 0.53498  dCc/dt =  -0.419D-06  res f3 =   0.257D-06
 Cd = 0.49158  dCd/dt =  -0.874D-06  res f4 =  -0.303D-07
 
  t =  200.0
 Ca = 0.31887  dCa/dt =   0.111D-06  res f1 =   0.135D-06
 Cb = 0.78388  dCb/dt =  -0.213D-07  res f2 =   0.209D-06
 Cc = 0.53498  dCc/dt =   0.166D-06  res f3 =   0.779D-07
 Cd = 0.49158  dCd/dt =  -0.165D-07  res f4 =   0.121D-07
 
 
 
 COMPUTATIONAL STATISTICS
 
 LAST STEP SIZE                        0.658D+02
 
 LAST ORDER OF THE METHOD                      1
 
 TOTAL NUMBER OF STEPS TAKEN                 147
 
 NUMBER OF FUNCTION EVALUATIONS              232
 
 NUMBER OF JACOBIAN EVALUATIONS                4
 
\end{verbatim}

The four residuals follow the semi-log variation with $t$ discussed
previously (eq. (30)).

\newpage

\begin{center}
{\Huge \textbf{References and Software}}
\end{center}

\vspace{0.5cm}

NAG routines pertaining to nonlinear equations include:
\begin{verbatim}
 
Chapter C02 - Zeros of Polynomials
 
  C02AFF - All zeros of complex polynomial,
           modified Laguerre method
 
  C02AGF - All zeros of real polynomial,
           modified Laguerre method
 
  C02AHF - All zeros of complex quadratic
 
  C02AJF - All zeros of real quadratic
 
Chapter C05 - Roots of One or More Transcendental
              Equations
 
  C05ADF - Zero of continuous function in given
           interval, Bus and Dekker algorithm
 
  C05AGF - Zero of continuous function, Bus and
           Dekker algorithm, from given starting
           value, binary search for interval
 
  C05AJF - Zero of continuous function, continuation
           method, from a given starting value
 
  C05AVF - Binary search for interval containing
           zero of continuous function (reverse
           communication)
 
  C05AXF - Zero of continuous function by continuation
           method, from given starting value (reverse
           communication)
 
  C05AZF - Zero in given interval of continuous function
           by Bus and Dekker algorithm (reverse communi-
           cation)
 
  C05NBF - Solution of system of nonlinear equations
           using function values only (easy-to-use)
 
  C05NCF - Solution of system of nonlinear equations
           using function values only (comprehensive)
 
  C05NDF - Solution of systems of nonlinear equations
           using function values only (reverse communi-
           cation)
 
  C05PBF - Solution of system of nonlinear equations
           using 1st derivatives (easy-to-use)
 
  C05PCF - Solution of system of nonlinear equations
           using 1st derivatives (comprehensive)
 
  C05PDF - Solution of systems of nonlinear equations
           using 1st derivatives (reverse communication)
 
  C05ZAF - Check user's routine for calculating 1st
           derivatives
 
\end{verbatim}

A Web page with additional information about the solution of nonlinear
equations:
\begin{verbatim}
 
http://www-fp.mcs.anl.gov/otc/Guide/OptWeb
 
     /continuous/unconstrained/nonlineareq/index.html
 
\end{verbatim}

Byrne, G. D., and C. A. Hall (1973), Numerical Solution of Systems of 
Nonlinear Algebraic Equations, Academic Press, New York.\newline
\newline
Dennis, J.E. and R.B. Schnabel (1986), Numerical Methods for Unconstrained 
Optimization and Nonlinear Equations, Prentice-Hall, Englewood Cliffs, NJ.%
\newline
\newline
Finlayson, B.A. (1980), Applied Nonlinear Analysis in Chemical Engineering, 
McGraw-Hill, New York.\newline
\newline
Riggs, J. B. (1994), An Introduction to Numerical Methods for Chemical 
Engineers, 2nd ed., Texas Tech University Press, Lubbock, Texas.\newline
\newline
Schiesser, W. E. (1994), Computational Mathematics in Engineering and
Applied  Science: ODEs, DAEs and PDEs, CRC Press, Boca Raton, FL.\newline

\end{document}
