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.
 A: 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.
Added
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 calculate 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:

*

*Numerical Analysis and Scientific Computation by Jeff Leader (link)


*Numerical Methods in Engineering with Python by Jaan Kiusalaas (link)
Both of these are highly readable and focus on giving you an understanding of the methods rather than just formulas for using them (which you can find in the LAPACK and Netlib implementations, in Numerical Recipes, or in dozens of other places.)
The first book focuses on examples in Matlab if you prefer that language. Most of the examples can also be directly used in Octave as well, which is a free, open-source Matlab clone. The second book focuses on showing examples in Python, and if you choose to go the Python route, I recommend working with the NumPy library because pure Python will be very inefficient.
