# Searching for numerical algorithm realization

Im searching for numerical algorithm realization for nonlinear equations system solver on PHP, C, C++, Java (with readable code :). Where I can find them?

Thx.

-
Possibly better suited to scicomp.stackexchange.com – lhf Apr 3 '12 at 1:08
The question is too broad. What kind of nonlinear system do you need to solve? – lhf Apr 3 '12 at 1:09
Im need to solve this $$\frac{s_2 - K_2}{ps_2^{\gamma_1} + (1 - p)s_2^{\gamma_2}} = \frac{K_1 - s_1}{ps_1^{\gamma_1} + (1 - p)s_1^{\gamma_2}} \\ \frac{s_2}{s_2 - K_2} = \frac{p\gamma_1s_2^{\gamma_1} + (1 - p)\gamma_2 s_2^{\gamma_2}}{ps_2^{\gamma_1} + (1 - p)s_2^{\gamma_2}} \\ \frac{s_1}{s_1 - K_1} = \frac{p\gamma_1s_1^{\gamma_1} + (1 - p)\gamma_2 s_1^{\gamma_2}}{ps_1^{\gamma_1} + (1 - p)s_1^{\gamma_2}}$$ $p,s_1,s_2$ - unknown variables here. Im need to solve this system many times with many different $\gamma,K$ etc – 0dd_b1t Apr 3 '12 at 8:08

Here is a link that outlines lots of different options, and the code that is linked there is generally well documented and readable. For nonlinear solvers, my opinion is that Netlib has better stuff (such as MINPACK) than LAPACK, but there are bound to be many other good implementations as well.

Most of this stuff will be in C or Fortran, but there are some examples in Matlab which should be closer to usual linear algebra form. You should also consider iterative methods like gradient descent, Newton's Method and semi-Newton methods like the Davidon-Fletcher-Powell method.

Since you expanded on the question in the comments, it seems clear that Newton's method will work for you. Your problem is small (only 3 variables) and has known analytical form so you can calculcate the Jacobian matrix $J$ by hand. By reformulating the update equation of Newton's method, you do not need to compute the inverse of $J$.. you'll only need a linear solver. (See here for more details).

In your case, the vector of solutions will be $x_{k}$, and you initialize $x_{0} = (p_{0}, s_{1,0}, s_{2,0})^{T}$ as the vector of initial guesses. Then just run the iterative scheme for $k=1,2,...$ :

$$x_{k}\textrm{ is solution to }J(x_{k-1})[x_{k} - x_{k-1}] = -F(x_{k-1})$$

where, with it understood that $x = (p,s_{1},s_{2})^{T}$ as above, so you'll rewrite the stuff below using $x_{i}$ to represent the $i$th thing you want to solve for:

$$F\biggl( (p,s_{1},s_{2})^{T} \biggr) = \begin{pmatrix} \frac{s_2 - K_2}{ps_2^{\gamma_1} + (1 - p)s_2^{\gamma_2}} - \frac{K_1 - s_1}{ps_1^{\gamma_1} + (1 - p)s_1^{\gamma_2}} \\ \frac{s_2}{s_2 - K_2} - \frac{p\gamma_1s_2^{\gamma_1} + (1 - p)\gamma_2 s_2^{\gamma_2}}{ps_2^{\gamma_1} + (1 - p)s_2^{\gamma_2}} \\ \frac{s_1}{s_1 - K_1} - \frac{p\gamma_1s_1^{\gamma_1} + (1 - p)\gamma_2 s_1^{\gamma_2}}{ps_1^{\gamma_1} + (1 - p)s_1^{\gamma_2}} \end{pmatrix}$$

You continue this until you reach some terminal criteria, such as hitting a maximum number of iterations or observing that the norm $||x_{k}-x_{k-1}||$ has become sufficiently small for your purposes, or that the value of the objective function at the current solution is small, $||F(x_{k})|| < tol$ for some tolerance $tol$.

These are all things that you can do with existing library functions in almost any scientific computing package. Almost any standard textbook on numerical analysis will offer a better explanation to this method, and how to check whether your particular problems falls into any of its pitfalls. Textbooks will also include references to computer code for solving such a problem. But your first step should be to calculate the Jacobian by hand in any case since it is available to you.

Two textbooks that I am familiar with and which should give good guidance for this problem are:

Way I see it, if the $\gamma_k$ are simple rational numbers, OP might even be able to use polynomial continuation methods... on the other hand, he's still pretty much stuck on finding a good starting point. – J. M. Apr 3 '12 at 19:34