JakobCHR.com
 
Quick Navigation:
 
Personal:
 Go to Home
 MS Research
 PhD Research
 Curriculum Vitae

General:
 Linux

Soon to come:
 Matlab
 On-line Stores
 Cycling
 Medicine & Health
 LaTeX
 OOP & C++
 Sony PCM-R500 DAT


next up previous contents index
Next: Core flow implementation Up: Core flow numerics Previous: Truncation error   Contents   Index


Solving systems of non-linear equations

During the calculation of the core flow it is necessary to solve a number of non-linear equations both single and multivariable.

In many practical problems it is impractical or impossible to obtain derivatives of the functions. The consequence is that we have to disregard methods, like the Newton-Raphson method which demand knowledge of derivatives.

In regard to the single variable equations which we have the form

f(x) = 0 (7.23)

we have chosen two different methods in order to obtain the zero. The first one is the secant method in which the new approximation to the zero is calculated from the two previous approximations (ie, a two-stage method) in the following way

\begin{displaymath}
x^{(k+1)} = \frac{ f(x^{(k)})x^{(k-1)} - f(x^{(k-1)})x^{(k)}}{f(x^{(k)})
- f(x^{(k-1)})}
\end{displaymath} (7.24)

where the superscripts indicate the iteration number. If the function f has several zeros it is important that the start guesses, x(0) and x(1), are close to the desired zero. However, in practice the function f has in many cases only one zero or the zeros are well separated such that reasonably accurate start guesses are easily established. The secant method has a relatively high convergence rate of 1.62 which compared to the second order convergence rate of the Newton-Raphson method seems reasonable.

The other method which we have used to solve single variable non-linear functions is the so-called Pegasus method. This method belongs to the group of enclosure methods which are characterized by the fact that they work on a range in which the zero is located. The Pegasus algorithm will not be repeated here--the details may be found in [57, p. 40]. The convergence order of the Pegasus method is found to be 1.64, ie a little better than the secant method. The Pegasus method is especially suited for problems where we know the range in which the zero is found. An example of such a situation is the problem of obtaining the point of void departure, zd (see sections 6.7.2.4 and 8.4).

When we move to the problem of solving systems of non-linear equations things really start to get complicated. The function we now want to solve has the form

\begin{displaymath}
\hspace{0.2ex}\underline{f}{}\hspace{0.15ex}(\hspace{0.2ex}...
...hspace{0.15ex}) = \hspace{0.2ex}\underline{0}{}\hspace{0.15ex}
\end{displaymath} (7.25)

where $\hspace{0.2ex}\underline{f}{}\hspace{0.15ex}$ is a Nth dimensional vector function and $\hspace{0.2ex}\underline{x}{}\hspace{0.15ex}$ a Nth dimensional vector. The simplest method (often called fixed point iteration [58]) which (in many cases) can solve the equation (7.27) is a one step method of the form

\begin{displaymath}
\hspace{0.2ex}\underline{x}{}\hspace{0.15ex}^{(k)} = \hspac...
...{0.15ex}(\hspace{0.2ex}\underline{x}{}\hspace{0.15ex}^{(k-1)})
\end{displaymath} (7.26)

where the vector function $\hspace{0.2ex}\underline{g}{}\hspace{0.15ex}$ arises a as the rest of $\hspace{0.2ex}\underline{f}{}\hspace{0.15ex}$ when one tries to isolate $\hspace{0.2ex}\underline{x}{}\hspace{0.15ex}$ in (7.27). Note that the function $\hspace{0.2ex}\underline{g}{}\hspace{0.15ex}$ is not unique and therefore we have a nearly infinite number of methods of the form (7.28) which behave differently. The success of this method depends on the Jacobian of the function $\hspace{0.2ex}\underline{g}{}\hspace{0.15ex}$, $\hspace{0.2ex}\underline{\underline{J}}{}\hspace{0.15ex}_g$7.4

. The convergence is assured only if the following requirement is fulfilled

\begin{displaymath}
\varrho(\hspace{0.2ex}\underline{\underline{J}}{}\hspace{0.15ex}_g) < 1
\end{displaymath} (7.27)

where $\varrho(\hspace{0.2ex}\underline{\underline{J}}{}\hspace{0.15ex}_g)$ denotes the spectral radius7.5 of the Jacobian. The convergence rate of the fixed point iteration is only linear with an error constant which depends on the spectral radius of the Jacobian--the smaller the spectral radius, the faster the convergence. A little gain in convergence rate can be established by altering the $\hspace{0.2ex}\underline{g}{}\hspace{0.15ex}$-function such that the individual function evaluations are computed in terms of the newest $\hspace{0.2ex}\underline{x}{}\hspace{0.15ex}$-values--this modification of the fixed point iteration is known as the Gauss-Seidel method for non-linear systems of equations.

In practice we really does not know if (7.29) is fulfilled and we cannot be sure that the fixed point iteration method can solve our problem. More complex methods exist which are more robust, ie which can solve a much larger class of functions. The most used methods are based on a least squares solution to the problem (7.27). The problem (7.27) may be reformulated into

\begin{displaymath}
\mbox{Find:}\qquad \mbox{arg} \min_{\hspace{0.2ex}\underlin...
...space{0.15ex}} F(\hspace{0.2ex}\underline{x}{}\hspace{0.15ex})
\end{displaymath} (7.28)

where the scalar function F is defined in terms of the vector function $\hspace{0.2ex}\underline{f}{}\hspace{0.15ex}$ in (7.27) by

\begin{displaymath}
F(\hspace{0.2ex}\underline{x}{}\hspace{0.15ex}) \;\hbox{$=$...
...space{0.15ex}(\hspace{0.2ex}\underline{x}{}\hspace{0.15ex})]^T
\end{displaymath} (7.29)

The function F is often minimized by utilizing an update of an approximation to the Jacobian of F (or the inverse of the $\hspace{0.2ex}\underline{\underline{J}}{}\hspace{0.15ex}$), eg the famous Broyden-update formula. Combined with the Levenberg-Marquardt method it is possible to guarantee global convergence (see [61]). These update methods all have a superlinear convergence rate, ie linear convergence with an error constant which asymptotically converges towards zero. The drawback is that implementing them require a considerable programming effort and as far as the author can see the real advantage is that they in general need fewer function evaluations when a high accuracy is required7.6.

In our case the only system of non-linear equations we have to solve is the system which arises in connection with our finite difference scheme. We observe that this system is of simple form which implies that function evaluations are cheap to acquire. Furthermore we have good start guesses for the iterations since we simply let the solution at the previous z-level be a start guess for the state at the next z-level. Therefore if the spectral radius of the $\hspace{0.2ex}\underline{g}{}\hspace{0.15ex}$-function is acceptably smaller than 1 we have indications that the very simple fixed point iteration suffices.

Practical calculations show that fixed point iteration actually behaves nicely and is in fact reasonably efficient in our case.


next up previous contents index
Next: Core flow implementation Up: Core flow numerics Previous: Truncation error   Contents   Index  
 
 
 
Revision 2.0, Copyright © 1999-2004 Jakob Christensen
http://www.JakobCHR.com
E-Mail: webmaster@JakobCHR.com
Top Quality
Developed with

Danish
Brain Power
Linux Powered!