Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I need to solve following system of nonlinear equations.

Here are some characteristics of this system:
It consists of n equation and n variables.
Every equation is in similar form -> sum of products = constant.
The lenght of every product is the same (it will be denoted as k).
The number of elements in sum might be different in each equation.
In every product in equation $i$ one of elements is $x_{i}$ - this is very important.
This system is "symmetrical", it means that if $x_{i} \cdot x_{j} \cdot ...$ is one of elements of equation i, then it is also in equaton j.
$b_{i} > 0$ - where $b_{i}$ is intercept in equation i.

I'll write an example of such equation for k=3 and n=6:
$x_{1} \cdot x_{3} \cdot x_{6} + x_{1} \cdot x_{2} \cdot x_{4} = b_{1}$
$x_{2} \cdot x_{1} \cdot x_{4} + x_{2} \cdot x_{5} \cdot x_{6} = b_{2}$
$x_{3} \cdot x_{1} \cdot x_{6} + x_{3} \cdot x_{4} \cdot x_{6} = b_{3}$
$x_{4} \cdot x_{1} \cdot x_{2} + x_{4} \cdot x_{3} \cdot x_{6} = b_{4}$
$x_{5} \cdot x_{2} \cdot x_{6} = b_{5}$
$x_{6} \cdot x_{1} \cdot x_{3} + x_{6} \cdot x_{2} \cdot x_{5} + x_{6} \cdot x_{3} \cdot x_{4} = b_{6}$

It is very easy to transform this equation to following form:
$x_{i} = b_{i} / something$ , $x_{i}$ is only on left-hand side of ith equation.

If we have all equation in such form then the fixed point is solution of it.
I've experimentally checked that algorithm analogic to Gauss-Seidel is covergent (i've checked ~100 random examples, and in every case it was convergent).
By analogic to Gauss-Seidel algorithm I mean:
1) Choose any initial solution $[x_{1}^{0} , ... , x_{n}^{0}]$
2.1) Calculate value of $x_{1}^{i+1}$ using $[x_{2}^{i} , ... , x_{n}^{i}]$
2.2) Calculate value of $x_{2}^{i+1}$ using $[x_{1}^{i+1} , ... , x_{n}^{i}]$
2.n) Calculate value of $x_{n}^{i+1}$ using $[x_{1}^{i+1} , ... , x_{n-1}^{i+1}]$
3) If solution is good enough stop, otherwise go to 2.1

I've tried Banach fixed point theorem, but is hard to say anything about spectral radius. Does anyone have a clue how to prove convergence of this algorithm?

share|cite|improve this question
This and this may have something you can use. – J. M. Dec 20 '10 at 10:03
If you just need to solve the system you can try resultants I suppose. – jerr18 Dec 20 '10 at 17:04
@J.M.: Unfortunately I dont have free access to articles right now, I will be able to download them after christmas, so I cant tell now if they are useful, but thanks anyway. @jerr18: Could You elaborate how resultants can help me in solving this system. – Tomek Tarczynski Dec 20 '10 at 21:08
Gröbner is useful here if exact solutions are what you need; since you seem to be happy with approximations, your Gauss-Seidel scheme is far simpler than what jerr18 is suggesting. – J. M. Dec 20 '10 at 21:39
@J.M: I'm happy with this approximation, but I'd have been even happier if I knew why and when it works. Simple simulations showed that it should be convergent in many cases (in random cases it was always convergent, but it doesn't prove anything). – Tomek Tarczynski Dec 21 '10 at 11:30

You just need to use a globalization strategy that will take you near a solution. From that point on, you can apply a fixed-point iteration. But in your case, you could just as well apply a variant of Newton's method with a line search. The idea is that you need to compute the Jacobian of your system (call is $J(x)$) and solve linear systems of the form $J(x_k) s_k = -F(x_k)$ for $s_k$ at each iteration ("$s$" is for step). Here $F(x_k)$ is your residual at iteration $k$, i.e., the left-hand side of your equations minus $b$.

The line search ensures that at each step, you achieve some decrease in a certain measure of progress such as $\|F(x)\|^2$. It computes a steplength $\alpha_k \in (0,1)$ such that your next iterate is given by $x_{k+1} = x_k + \alpha_k s_k$.

As you get close to a solution, if $J(x_k)$ never approaches singularity, the method will naturally select $\alpha_k = 1$ and you'll get the quadratic convergence of Newton's method.

You need to think about how you'll solve each linear system $J(x) s = -F(x)$. You can either factorize $J(x)$ (you'll need to use the LU factorization unless you know more about $J$) or you can solve the system iteratively (for instance using a Krylov subspace method such as GMRES or a more specific one if you have specific structural knowledge on $J$). In this case, you'll need to implement a variant of what I outlined above called the inexact Newton method, where each $s_k$ is computed inaccurately.

The two books below by Tim Kelley will have all the information you need. He also distributes Matlab code that implements the methods he presents in the books.

share|cite|improve this answer
Since his Jacobian is known to be symmetric but not necessarily invertible, MINRES is more appropriate than GMRES. Good answer, otherwise. – user7530 Nov 24 '11 at 21:45

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.