# Avoiding a For-loop when finding the solution to a set of polynomial equations

There are several examples and questions regarding Map, but I couldn't find what I need.

This is a minimal working example.

I have two functions

$\qquad f=x^4+y^4-1$,

and

$\qquad g=x^2-4$.

One way to find the solution set would be to find the values of $x$ from $g=0$ and then substitute those values into $f$. Then solve for $f=0$ to obtain the values of $y$. I need to club the corresponding values of $x$ and $y$ together to form the solution set.

I have correctly implemented a solution using a For loop. However, I was wondering if I could use a nested Map or something simpler. I come from a C++ background and a For loop is the intuitive and normal solution for most work involving array manipulation. However, I think For loops hurt writing efficient code in Mathematica (I suppose).

Existing solution:

f = x^4 + y^4 - 1;
g = x^2 - 4;
solx = Solve[g == 0, x];
For[i = 1, i <= Length[solx], i++,
solytemp = NSolve[(f /. solx[[i]]) == 0, y];
temp = Flatten[{solx[[i]], #}] & /@ solytemp;
If[i == 1, solxy = temp, AppendTo[solxy, #] & /@ temp;];
];
solxy

 {{x -> -2, y -> -1.39158 - 1.39158 I}, {x -> -2, y -> -1.39158 + 1.39158 I},
{x -> -2, y -> 1.39158 - 1.39158 I}, {x -> -2, y -> 1.39158 + 1.39158 I},
{x -> 2, y -> -1.39158 - 1.39158 I}, {x -> 2, y -> -1.39158 + 1.39158 I},
{x -> 2, y -> 1.39158 - 1.39158 I}, {x -> 2, y -> 1.39158 + 1.39158 I}}


Please note that this is an MWE and the actual polynomials may be of much higher order.

### Update 1

The solutions are to be found in steps and then clubbed together, as the polynomials are big. I also use my own modules instead of NSolve depending on how the solution is to be obtained and the problem.

I am aware of using Solve for multiple polynomials, however I am refraining from using that as the system is huge.

Solve[{f == 0, g == 0}, {x, y}]


I was hoping if something like below could be used (the code below does not work)

{#1, #2} & /@ NSolve[(f /. #1) == 0, y] & /@ solx


### Update 2

Using the answer given by @Jason, I could do the following to get the result:

f = x^4 + y^4 - 1
g = x^2 - 4
solx = Solve[g == 0, x];
solxy = Flatten[
solf2[xtrial_] :=
Union[#, xtrial] & /@ NSolve[(f /. xtrial) == 0, y]; (solf2[#] & /@
solx), 1];


Is there a way to bypass assigning a function. I don't want to end up creating several functions. It is a good solution nevertheless.

## 6 Answers

An alternative approach to obtain same answer in steps is

Solve[g == 0, x];
N[Solve[f == 0 /. #]] & /@ %;
MapIndexed[Join[%%[[First[#2]]], #] &, %, {2}];
Flatten[%, 1]

(*  {{x -> -2, y -> -1.39158 - 1.39158 I}, {x -> -2, y -> 1.39158 - 1.39158 I},
{x -> -2, y -> -1.39158 + 1.39158 I}, {x -> -2, y -> 1.39158 + 1.39158 I}},
{{x -> 2, y -> -1.39158 - 1.39158 I}, {x -> 2, y -> 1.39158 - 1.39158 I},
{x -> 2, y -> -1.39158 + 1.39158 I}, {x -> 2, y -> 1.39158 + 1.39158 I}}  *)


Explanation

In response to the request below by the OP, details follow. The first two lines of code provide the x and y solutions.

(* {{x -> -2}, {x -> 2}} *)
(* {{{y -> -1.39158 - 1.39158 I}, {y -> 1.39158 - 1.39158 I},
{y -> -1.39158 + 1.39158 I}, {y -> 1.39158 + 1.39158 I}},
{{y -> -1.39158 - 1.39158 I}, {y -> 1.39158 - 1.39158 I},
{y -> -1.39158 + 1.39158 I}, {y -> 1.39158 + 1.39158 I}}} *)


The goal, now, is to associate the eight values for y with corresponding values of x. To access each of the eight values individually, we must reach two levels down into the nested List above. The third argument, {2}, of MapIndexed accomplishes this. (If {1} or no third argument at all were used, MaxIndexed would process blocks of four y values at once, which is not what we desire.) Next, note that, unlike Map, MapIndexed returns not only y values (represented by # or #1) but also their positions in the nested List (represented by #2). So, for instance, MapIndexed returns not just {y -> -1.39158 + 1.39158 I} but also its index {2, 3}. The First[#2] index, 2 in this case, is the index of the corresponding x value. Finally, Join links them together, and Flatten eliminates an unneeded pair of {}.

• Can you please explain the code a bit? Specially the #2 and {2} at the end? – Marvin Jan 2 '16 at 7:39
• @Saurav I have added an explanation to my answer. Best wishes for the New Year. – bbgodfrey Jan 2 '16 at 17:10

You could try something like this, create a function that solves f for a given x

solf[xtrial_Rule] :=
Prepend[#, xtrial] & /@ NSolve[(f /. xtrial) == 0, y]


which works like this

solf[x -> 2]
(* {{x -> 2, y -> -1.39158 - 1.39158 I}, {x -> 2,
y -> -1.39158 + 1.39158 I}, {x -> 2,
y -> 1.39158 - 1.39158 I}, {x -> 2, y -> 1.39158 + 1.39158 I}} *)


and then Apply it to a list of replacement rules,

(solf @@@ solx)
(* {{{x -> -2, y -> -1.39158 - 1.39158 I}, {x -> -2,
y -> -1.39158 + 1.39158 I}, {x -> -2,
y -> 1.39158 - 1.39158 I}, {x -> -2,
y -> 1.39158 + 1.39158 I}}, {{x -> 2,
y -> -1.39158 - 1.39158 I}, {x -> 2,
y -> -1.39158 + 1.39158 I}, {x -> 2,
y -> 1.39158 - 1.39158 I}, {x -> 2, y -> 1.39158 + 1.39158 I}}} *)


this has one extra dimension of list than your original solution, so you can Flatten once,

(solf @@@ solx)~Flatten~1
(* {{x -> -2, y -> -1.39158 - 1.39158 I}, {x -> -2,
y -> -1.39158 + 1.39158 I}, {x -> -2,
y -> 1.39158 - 1.39158 I}, {x -> -2,
y -> 1.39158 + 1.39158 I}, {x -> 2,
y -> -1.39158 - 1.39158 I}, {x -> 2,
y -> -1.39158 + 1.39158 I}, {x -> 2,
y -> 1.39158 - 1.39158 I}, {x -> 2, y -> 1.39158 + 1.39158 I}} *)


Edit If you don't want to have a defined function clogging up the main namespace, you can use a pure function like this,

Flatten[
Function[xtrial,
Prepend[#, xtrial] & /@ NSolve[(f /. xtrial) == 0, y]] @@@ solx,
1]
(* {{x -> -2, y -> -1.39158 - 1.39158 I}, {x -> -2,
y -> -1.39158 + 1.39158 I}, {x -> -2,
y -> 1.39158 - 1.39158 I}, {x -> -2,
y -> 1.39158 + 1.39158 I}, {x -> 2,
y -> -1.39158 - 1.39158 I}, {x -> 2,
y -> -1.39158 + 1.39158 I}, {x -> 2,
y -> 1.39158 - 1.39158 I}, {x -> 2,
y -> 1.39158 + 1.39158 I}} *)

• Yeah, this works for the MWE provided. However, I have to find the solutions in steps, as the polynomials are big. We also use our own Modules instead of NSolve depending on how the solution is to be obtained and the problem. – Marvin Dec 31 '15 at 9:14
• I guess this works! I will wait for sometime to see if someone posts a better solution and then accept. – Marvin Dec 31 '15 at 9:39

MapIndexed generally makes a good replacement for For as it provides the indexed as it maps.

f = x^4 + y^4 - 1;
g = x^2 - 4;
solx = Solve[g == 0, x];

Flatten[
MapIndexed[
Function[{value, index},
Flatten /@ Thread[{First@solx[[First@index]], #} &[value]]],
Map[NSolve[# == 0, y] &, (f /. #) & /@ solx]],
1]


{{x -> -2, y -> -1.39158 - 1.39158 I}, {x -> -2, y -> -1.39158 + 1.39158 I}, {x -> -2, y -> 1.39158 - 1.39158 I}, {x -> -2, y -> 1.39158 + 1.39158 I}, {x -> 2, y -> -1.39158 - 1.39158 I}, {x -> 2, y -> -1.39158 + 1.39158 I}, {x -> 2, y -> 1.39158 - 1.39158 I}, {x -> 2, y -> 1.39158 + 1.39158 I}}

Hope this helps.

You say: The solutions are to be found in steps and then clubbed together... Here is my suggestion:

solx = Solve[x^2 - 4 == 0, x];
soly = NSolve[x^4 + y^4 - 1 == 0, y] /. solx;
solution = {solx[[#]], soly[[#]]} & /@ Range[1, Length@solx]

{{{x -> -2}, {{y -> -1.39158 - 1.39158 I}, {y ->
1.39158 - 1.39158 I}, {y -> -1.39158 + 1.39158 I}, {y ->
1.39158 + 1.39158 I}}},

{{x -> 2}, {{y -> -1.39158 - 1.39158 I}, {y ->
1.39158 - 1.39158 I}, {y -> -1.39158 + 1.39158 I}, {y ->
1.39158 + 1.39158 I}}}}


Solve or Reduce can handle multiple equations

f[x_, y_] = x^4 + y^4 - 1;

g[x_] = x^2 - 4;

eqns = {f[x, y], g[x]} == 0 // Thread;


Exact solutions

soln = Solve[eqns, {x, y}]

(*  {{x -> -2, y -> -(-15)^(1/4)},
{x -> -2, y -> (-I)*(-15)^(1/4)},
{x -> -2, y -> I*(-15)^(1/4)},
{x -> -2, y -> (-15)^(1/4)},
{x -> 2, y -> -(-15)^(1/4)},
{x -> 2, y -> (-I)*(-15)^(1/4)},
{x -> 2, y -> I*(-15)^(1/4)},
{x -> 2, y -> (-15)^(1/4)}}  *)


Verification

eqns /. soln

(*  {{True, True}, {True, True}, {True, True}, {True, True}, {True, True}, {True,
True}, {True, True}, {True, True}}  *)


Or for more complex systems of equations you may need to use NSolve for approximate solutions

(soln2 = NSolve[eqns, {x, y}, WorkingPrecision -> 15]) // N

(*  {{x -> -2., y -> -1.39158 - 1.39158 I}, {x -> -2.,
y -> -1.39158 + 1.39158 I}, {x -> -2., y -> 1.39158 - 1.39158 I}, {x -> -2.,
y -> 1.39158 + 1.39158 I}, {x -> 2., y -> -1.39158 - 1.39158 I}, {x -> 2.,
y -> -1.39158 + 1.39158 I}, {x -> 2., y -> 1.39158 - 1.39158 I}, {x -> 2.,
y -> 1.39158 + 1.39158 I}}  *)


Verification

eqns /. soln2

(*  {{True, True}, {True, True}, {True, True}, {True, True}, {True, True}, {True,
True}, {True, True}, {True, True}}  *)


Howzabout

Remove[solveAndBacksubstitute];
solveAndBacksubstitute[{pair : {_, _}}, opts___] :=
NSolve[pair[] == 0, pair[], opts]
solveAndBacksubstitute[pairs : {{_, _} ..}, opts___] := Module[{
partlsoln, recurse },
partlsoln = NSolve[#[] == 0, #[], opts] &[First[pairs]];
recurse[varValue_] := Join[varValue, #] & /@
solveAndBacksubstitute[Rest[pairs] /. varValue];
Join @@ (recurse /@ partlsoln)
]


Then, usage:

f = x^4 + y^4 - 1;
g = x^2 - 4;
solveAndBacksubstitute[{{g, x}, {f, y}}]
(*  {{x -> -2., y -> -1.39158 - 1.39158 I},
{x -> -2., y -> -1.39158 + 1.39158 I},
{x -> -2., y ->  1.39158 - 1.39158 I},
{x -> -2., y ->  1.39158 + 1.39158 I},
{x ->  2., y -> -1.39158 - 1.39158 I},
{x ->  2., y -> -1.39158 + 1.39158 I},
{x ->  2., y ->  1.39158 - 1.39158 I},
{x ->  2., y ->  1.39158 + 1.39158 I}}  *)
solveAndBacksubstitute[{{g, x}, {f, y}}, WorkingPrecision -> 8]
(*  {{x -> -2.0000000, y -> -1.3915788 - 1.3915788 I},
{x -> -2.0000000, y -> -1.3915788 + 1.3915788 I},
{x -> -2.0000000, y ->  1.3915788 - 1.3915788 I},
{x -> -2.0000000, y ->  1.3915788 + 1.3915788 I},
{x ->  2.0000000, y -> -1.3915788 - 1.3915788 I},
{x ->  2.0000000, y -> -1.3915788 + 1.3915788 I},
{x ->  2.0000000, y ->  1.3915788 - 1.3915788 I},
{x ->  2.0000000, y ->  1.3915788 + 1.3915788 I}}  *)


Commentary: There's no For[] loop. Also, since you are using a complicated inner function (NSolve[]) for which there are many options that you will eventually find you need to tweak, options are received by solveAndBacksubstitute[] and passed straight to NSolve[]. (We could be fancy and filter out only those options that NSolve[] understands, but I don't see much value for this kind of complexity here.)

I can see a few directions for more fussiness:

• We could use a TimeConstrained[] and MemoryConstrained[] Solve[] to try to get simple, exact answers and fall back to NSolve[] if the time or intermediate expression swell is too great.
• We could feed in a sequence of solve functions, for instance "{Solve, NSolve, Solve, DSolve, Solve, NDSolve}", to indicate which solve function to use for each of the steps of processing the list.
• We could write a helper function to go from "{func1, func2, ...}, {var1, var2, ...}" to "{{func1, var1}, {func2, var2}, ...}", but since this is just Transpose[{{func1, func2, ...}, {var1, var2, ...}}], I didn't bother.