How to determine the Jacobian when solving a nonlinear system by Newton-Raphson method

Question

I am trying to solve some nonlinear systems by Newton's method and the solution accuracy is very important to my problem.

1. without using the symbolic computation softwares, how can I compute the Jacobian of a general nonlinear system of polynomials via C++ or other similar programming languages? The difficulties for me are mainly:

   • as accurate as symbolic Jacobian
   • an algorithm suitable for general nonlinear system cases
   • only dependent on C++ or similar programming languages;

2. If I have to use finite difference method to obtain an approximate Jacobian, how the step size chosen would affect the final solution accuracy? how to determine the step size so that I can obtain best solution accuracy under the same computation precision level? How to determine (quantitatively) the effects of approximate Jacobian on the accuracy of the final solution?

Answer 1

Check out the idea of dual numbers, there exist several example implementations in C++. Evaluating the functions using appropriately initialized dual numbers results in the evaluation of one directional derivative. Repeat this for all coordinate directions.

For a nice introduction, see Piponi: AD, C++ and Photogrammetry (http://el.mdu.edu.tw/datacos/09820722022O/paper.pdf)

If L is the effort to evaluate the functions, one directional derivative costs about 3*L, the full jacobian 3n*L. If one were to combine all the directions in one evaluation, this reduces to (1+2n)*L. But by then we have entered the realm of automatic or algorithmic differentiation.

Look for FADBAD/TADIFF for an easy implementation, for really fast codes one uses code transformations as provided by the Tapenade project. ADOL-C is in between, more automatic than Tapenade, faster than FADBAD.