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

Question

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.

Answer 1

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 {}.

Answer 2

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}} *)

Answer 3

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.

Answer 4

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}}}}

Answer 5

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}}  *)

Answer 6

Howzabout

Remove[solveAndBacksubstitute];
solveAndBacksubstitute[{pair : {_, _}}, opts___] := 
  NSolve[pair[[1]] == 0, pair[[2]], opts]
solveAndBacksubstitute[pairs : {{_, _} ..}, opts___] := Module[{
    partlsoln, recurse },
  partlsoln = NSolve[#[[1]] == 0, #[[2]], 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.