# Closed-Form solution for system of simple nonlinear equations

I am interested in analytical solutions for a system of nonlinear equations.

(The question was first asked at math.SE, where (after 1months and one rounds of bounty) there is only interesting analysation that the problem corresponds to high-order polynomial. I'm hoping to get more insight here, in particular whether there are methods which can deal with such high-order polynomials or whether it is likely to be not possible to solve it).

Motivation: The source of the question is a very convinient method to create random matrices with special properties. Mathematica can give me solutions up to certain sizes of the linear system, but I would like to have it for arbitrary size N. I can also use numerical algorithms (which I am doing at the moment), but for N in the order of $N\approx10.000$, they are quite slow.

System of nonlinear equations:

$$(w_i \cdot \sum_{j=1}^N w_j) - w_i^2 = d_i$$ for $i=1...N$, and $w$ and $d$ are vectors with $N$ dimensions, and $w_i$ and $d_i$ is the $i$-th component of the vector. Both $d_i$ and $w_i \in \mathbb{R_+}$. I am providing the vector $d$ (i.e. N real non-negative numbers), and want to solve for $w_i$.

How can $w_i$ be written down in a closed-form for a given vector $d$, for arbitrary N?

Example: If N=3 we have the following system of equations:

$$w_1 \cdot (w_2 + w_3) = d_1 \\ w_2 \cdot (w_1 + w_3) = d_2 \\ w_3 \cdot (w_1 + w_2) = d_3$$

with $w_i, d_i \in \mathbb{R}$. For a given vector $d=(d_1,d_2,d_3)$, I want to get $w=(w_1,w_2,w_3)$.

Rewriting: There is a nice way to rewrite the question, but I am not sure whether it actually helps:

Let's set $c=\sum_{j=1}^N w_j$, which is the sum of all weights. What we have now:

$$c \cdot w_i - w_i^2 = d_i \\ w_i^2 - c \cdot w_i + d_i = 0$$ which has two solutions:

$$w_{i_{1,2}} = \frac{c}{2} \pm \sqrt{ \left(\frac{c}{2}\right)^2 - d_i}$$ and the normalisation constant $c$ can be calculated by the sum of all weights:

$$\sum_{j=1}^N w_j = \sum_{j=1}^N \left(\frac{c}{2} \pm \sqrt{ \left(\frac{c}{2}\right)^2 - d_j} \right) = c$$

Is there a closed-form solution for c?

Analysation from math.SE: For completness, I copy the analysation of Robert Israel at math.SE:

$n=3$ is rather easy: $w_i$ satisfies a quadratic polynomial.

For $n=4$, each $w_i$ satisfies a rather nasty polynomial of degree $8$ (but involving only even powers). Thus there is a solution in terms of radicals, but it won't be pleasant.

For $n=5$, it seems each $w_i$ satisfies a polynomial of degree $22$. A solution in radicals is not to be expected. Thus with $d_1 = -7, d_2 = 3, d_3 = 9, d_4 = 7, d_5 = 8$, $w_5$ satisfies $$81\,{w_{{5}}}^{22}+7074\,{w_{{5}}}^{20}+198792\,{w_{{5}}}^{18}+2887764 \,{w_{{5}}}^{16}+23487600\,{w_{{5}}}^{14}+99587008\,{w_{{5}}}^{12}+ 69082752\,{w_{{5}}}^{10}-1402988992\,{w_{{5}}}^{8}-6995300352\,{w_{{5} }}^{6}-14191984640\,{w_{{5}}}^{4}-13095665664\,{w_{{5}}}^{2}- 4492099584 = 0$$

EDIT: This is an irreducible polynomial of degree $11$ in $x = w_5^2$. Maple doesn't do Galois groups of polynomials of degree $11$, but GAP does, and confirms that its Galois group is $S_{11}$. In particular, there is no solution in radicals.

But, anyway, one practical method that sometimes works is to find an iteration that converges. My first attempt is this: $$w'_j = \frac{d_j + w_j^2}{\sum_{i=1}^n w_i}.$$ For example, if $n=4$ and $d=(1,\frac12,3,4)$, starting with $w=(1,1,1,1)$ it converges to $w=(0.2133500083,0.1042484319,0.7171275770,3.865757626)$. However the convergence is rather slow and I guess it will be even slower for large $n$. Another possibility is Newton-Raphson iteration.
It is convenient to consider this as a problem of finding a solution of the equation $$(N-2)x-\sum_i\delta_i\sqrt{x^2-d_i}=0,$$ for a given vector $d$ and another vector $\delta_i=\pm 1$ which may be chosen. (Here $c=\sum_i w_i=2x$.) After we multiply the left hand sides for all $2^N$ possible choices of $\delta$, we will have a polynomial.
Because of some cancellation, the degree of this polynomial is not $2^N$, but $2^N-2N$. (The leading coefficient is a function of $d$, so the degree may drop for special $d$.) Note that, for a general $d$, it means $2^N-2N$ different solutions of the original problem. (Some of them complex.)
The point is, even if this problem had an algebraic solution (which is not impossible), it would not likely to be of any use in practice because of the exponential number of solutions. Not for $N\approx 10^4$, at least.