Performance issues

In implementing Newton's method using a template file formulation, greater control over the convergence of the scheme (and convergence properties) is possible than using Mathematica's FindRoot. At the same time, much of the generality of the built-in function is retained.

The majority of computational time in Newton's method is spent on LU decomposition and updating of the Jacobian matrix at each iteration. If the matrix is not updated at each iteration quadratic convergence is lost, but similar convergence properties are still retained in the neighbourhood of the solution.

Convergence of the iteration scheme can be controlled by:

1.

Moving the definition of the Jacobian matrix and the call to the LU decomposition routine outside the main iteration loop (i.e. to before the do loop labelled 20 in the code below). This means that a fixed Jacobian matrix is used throughout the iteration process (sometimes referred to as modified Newton-Raphson iteration). This approach is not recommended if the starting value is known to be far from the true solution.

2.

Adding a clause to ensure that the Jacobian matrix is not updated after a specified number of iterations. If the problem is well conditioned, the Jacobian matrix will not change significantly in value when the scheme is near convergence.

3.

Taking account of any special structure arising in the problem - such as a tri-diagonal, banded, sparse or symmetric Jacobian matrix. This decreases the complexity of the linear solver and allows for greater scope in parellisation. For special case linear solution routines see for example [1,22]. 

For more information about the linear solution routines used in our example, see [1]. These are a set of machine independent codes in the public domain [23] for linear algebra known as LAPACK.³ The file linsolv.f contains the LAPACK general linear solver and LU decomposition routine which does not utilise property 3 above. Special purpose routines are available for other matrix structures. LAPACK uses Basic Linear Algebra Subprograms (BLAS) with machine specific implementations, thus providing near optimal efficiency for each platform. The necessary BLAS are given as FORTRAN code in the file linsolv.f and should be removed if machine specific versions are available.

The file linsolv.f should be linked with newton.f (the file which results from the application of Splice to newton.mf) at compile time. When the system of equations becomes large (typically greater than 50x50) iterative refinement of the solution may become necessary. It is straightforward to implement iterative refinement using the LAPACK routine dgerfs.f.