# Tricks for solving (lots of) coupled nonlinear equations numerically?

I have a system of 6 non-linear (quadratic) coupled equations with 6 complex unknowns

\begin{align*} |x_1|^2 + |x_2|^2 + |x_3|^2 &= a\\ x_1 x_4^* + x_3 x_5^* &= b + c i\\ x_1 x_6^* &= d + e i\\ |x_4|^2 + |x_5|^2 &= f\\ x_4 x_6^* &= g + h i\\ |x_6|^2 &= k \end{align*}

{Abs[x1]^2 + Abs[x2]^2 + Abs[x3]^2 == a,
x1 Conjugate[x4] + x3 Conjugate[x5] == b + c I,
x1 Conjugate[x6] == d + e I,
Abs[x4]^2 + Abs[x5]^2 == f,
x4 Conjugate[x6] == g + h I,
Abs[x6]^2 == k}


where $x_1$ etc. are complex variables, $x_1^*$ is Conjugate[x1], $|x|$ is Abs[x], and a etc. are real constants. The equations are probably underdetermined - if you break them into real and imaginary parts, there are only 9 equations for 12 real unknowns. Physically (since these equations solves for parameters in a model) I expect some phases in the $x_i$ are not physical (i.e. some $x_i$ are actually real numbers instead of complex), but I don't know a good way to parametrize the phases away.

NSolve[{eqns}, {vars}, Complexes] has been running for more than 10 minutes (and still running!). I've never encountered numerical equation solving with coupled non-linear equations before, so I don't know what to expect. Is there anything I can do to optimize the process? What other numerical solving functions can I try?

Also, probably off-topic here: I am under the impression that numerics is not Mathematica's strong suit; should I seek out other programs to solve this problem? Recommendations?

-
Your impression is wrong, Mathematica numerics is very strong. Since you haven't included your equations I doubt you'll find any help besides adequate links e.g. reference.wolfram.com/mathematica/tutorial/… –  Artes Nov 18 '12 at 2:49
You haven't even said if your equations are algebraic or transcendental, so there's really nothing much that can be said by us... –  Guess who it is. Nov 18 '12 at 3:24
I did say quadratic equations. I'll post the actual equations. –  polyglot Nov 18 '12 at 3:32
NSolve cannot run with symbolic parameters. But you did not give any numerical values to a, b, c, etc. Why not to split them in 9 real equations and solve symbolically for 9 real variables in terms of 3 variables and parameters? Then you can research symbolically domains of possible solutions. –  Vitaliy Kaurov Nov 18 '12 at 6:01

Also, probably off-topic here: I am under the impression that numerics is not Mathematica's strong suit; should I seek out other programs to solve this problem? Recommendations?

My experience is that Mathematica can solve a set of numerical equations with 10.000 and more equations in about 0.5-1 sec. These are huge but linear equations. I use this feature in my software, TIMO Structural, based on Mathematica. I do not know of any other universal system that is so powerful. (I worked with Maple for many years and I know what I'm talking about).

Please explain what is the reason for using the Abs[arg] function in the expression:

Abs[x1]^2 + Abs[x2]^2 + Abs[x3]^2 == a

-
Well, linear equations are a very different matter. Still, interesting to hear about your experience with Maple... serious users of both that and Mathematica seem to be few and far between. –  Oleksandr R. Nov 26 '12 at 9:59
Thanks for offer, Oleksandr. It is true I was happy to work with M & M :). The best idea is combine the advantages of each into Tools for solving engineering (science) tasks of industrial level real problems (not one time - three equations research task in mode by hand). Globally in short, Maple is symbolic system, Mathematica - universal. If it is interesting I can to participate in post M vs M or something like this. –  Siarhei A Arlou Nov 28 '12 at 22:45
That might be good for a blog post. We're going to have a MATLAB vs. Mathematica comparison/tutorial post at some point, so it would make sense to go the other way and have something about Maple as well. –  Oleksandr R. Nov 29 '12 at 3:07
@SiarheiAArlou I'm curious on how does it compare to specialized code (obviously, on this topic...). 10x... 100x... slower? A parallelisation difficulty? –  P. Fonseca Aug 27 '14 at 6:56

Not sure why this resurfaced but it can be done symbolically by separating real and imaginary parts. This of course is no guarantee that for some regions on parameter space the solution values will actually be real valued.

zz = Array[x, 6] + I Array[y, 6];
polys = Expand[
ComplexExpand[{zz[[1 ;; 3]].Conjugate[zz[[1 ;; 3]]] - a,
zz[[1 ;; 3 ;; 2]].Conjugate[zz[[4 ;; 5]]] - (b + c I),
zz[[1]] Conjugate[zz[[6]]] - (d + e I),
zz[[4 ;; 5]].Conjugate[zz[[4 ;; 5]]] - f,
zz[[4]] Conjugate[zz[[6]]] - (g + h I),
zz[[6]] Conjugate[zz[[6]]] - k}]];
p2 = Flatten[ComplexExpand[Map[{Re[#], Im[#]} &, polys]]]

(* Out[85]= {-a + x[1]^2 + x[2]^2 + x[3]^2 + y[1]^2 + y[2]^2 +
y[3]^2, 0, -b + x[1] x[4] + x[3] x[5] + y[1] y[4] + y[3] y[5], -c +
x[4] y[1] + x[5] y[3] - x[1] y[4] - x[3] y[5], -d + x[1] x[6] +
y[1] y[6], -e + x[6] y[1] - x[1] y[6], -f + x[4]^2 + x[5]^2 +
y[4]^2 + y[5]^2, 0, -g + x[4] x[6] + y[4] y[6], -h + x[6] y[4] -
x[4] y[6], -k + x[6]^2 + y[6]^2, 0} *)

Timing[soln = Solve[p2 == 0, Variables[zz]];]

(* Out[96]= {1.341609, Null} *)


The solution is pretty large since it involves symbolic radicals.

No idea what NSolve would do with this, other than churn its wheels. It's not a good fit for a numeric solver since there are symbolic parameters. So NSolve basically punts to some primordial Solve code, figuring it will apply N if a result ever comes back.

-
Oh, now I see how this rebounded. Good that the irrelevant piggyback was removed though. Wrong forum, wrong method of posting, wrong phase of the moon... –  Daniel Lichtblau Nov 17 '14 at 4:52

Why are your variable symbols single symbols, while your constants (a,b,c,..) split into real and imaginary parts? That's really the problem here. You really think Mathematica knows that you want (for example) the real part of the second equation to be b and the complex part to be c? You said the domain was Complexes which means it's considering b and c to be complex numbers too.

These are not quadratic equations in some complex variables. (That would imply, to most readers, that you have a polynomial and are not using complex conjugates.) If you really wanted a true set of quadratics, Mathematica could NSolve or even Solve that easily:

RI[] := RandomInteger[{1, 6}]
RC[] := RandomReal[] + I*RandomReal[]
Timing[Length[NSolve[Table[
Sum[RC[] x[RI[]] x[RI[]], {dummy2, 1, RI[] RI[]}] == RC[],
{dummy,1, 6}]]]]


This solves random piles of complex quadratics in fractions of seconds.

Having just refreshed to see it, I am saying that I agree strongly with Vitaliy's comment. You are really looking at nine real equations here (twelve, but three of which tell you nothing, since three of your equations have no imaginary part and spin off 0=0 for that). If your parameters (a,b,c,d..) are in re/im parts, you should treat your variables as such.

That will also inform you that your system is underdetermined: You really only have nine equations and twelve unknowns. That will make it hard for anyone -- Solve, NSolve, or any other method. I think this deserves proper attention in more than just a comment, since that's the real source of your trouble more (in addition to the issues about how to properly implement some solution method by computer).

-