# How can I solve an equation numerically for a range of x-values and multiple initial root guesses?

The following implicitly defined equation has both real and complex roots:

$$\frac{0.5}{(y + x)^2 - 0.3x^2} + \frac{0.5}{(y - x)^2 - 0.3x^2} - \frac{1}{x^2} = 1$$

I am trying to solve it numerically using FindRoot for a range of x-values. For the initial guesses, I am using Solve. The equation has four solutions, namely, $$y_1, y_2, y_3$$, and $$y_4$$,

Here is my attempt:

F[y_, x_] :=  0.5/((y - x)^2 - 0.3*x^2) + 0.5/((y - x)^2 - 0.3*x^2) - 1/x^2 - 1

Intialguess = Grid[y /. Table[Solve[F[y, K] == 0, y], {K, 0.1, 1, 0.1}]]

Table[FindRoot[F[y, k] == 0, {y, Intialguess}], {k, 0.1, 1, 0.1}]


Output:

Out[1] =
-0.19094    0. - 0.0236008 I   0. + 0.0236008 I   0.19094
-0.380117   0. - 0.0450256 I   0. + 0.0450256 I   0.380117
-0.566053   0. - 0.062029 I    0. + 0.062029 I    0.566053
-0.747747   0. - 0.0720716 I   0. + 0.0720716 I   0.747747
-0.924725   0. - 0.0715281 I   0. + 0.0715281 I   0.924725
-1.09698    0. - 0.0515827 I   0. + 0.0515827 I   1.09698
-1.26484    -0.0550617         0.0550617          1.26484
-1.42884,   -0.112579          0.112579           1.42884
-1.58956    -0.163747          0.163747           1.58956
-1.74762    -0.214101          0.214101           1.74762

in search specification {y,Intialguess} is not a number or array of \
numbers. >>

{FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}],
FindRoot[F[y, k, U1, U2] == 0, {y, Intialguess}]}


The first, second, third, and fourth columns of the grid of initial guesses containing real and complex-valued solutions represents $$y_1, y_2, y_3$$, and $$y_4$$, respectively. How do I convert the initial guesses so that FindRoot recognises each grid entry as a numeric value?

Goal: Produce a table of numerically solved-roots:

    x        y1             y2                y3             y4
0.1    -0.19494   0. - 0.0236008 I   0. + 0.0236008 I   0.19094
0.2    -0.380117  0. - 0.0450256 I   0. + 0.0450256 I   0.380117
0.3    -0.566053  0. - 0.062029 I    0. + 0.062029 I    0.566053
0.4    -0.747747  0. - 0.0720716 I   0. + 0.0720716 I   0.747747
0.5    -0.924725  0. - 0.0715281 I   0. + 0.0715281 I   0.924725
0.6    -1.09698   0. - 0.0515827 I   0. + 0.0515827 I   1.09698
0.7    -1.26484     -0.0550617        0.0550617         1.26484
0.8    -1.42884,    -0.112579         0.112579          1.42884
0.9    -1.58956     -0.163747         0.163747          1.58956
1.0    -1.74762     -0.214101         0.214101          1.74762


and plot the results.

• Are you aware that Root objects can enumerate polynomial roots exactly? See documentation. Aug 4 at 14:11
• Thank you for your suggestion, but I am not looking for an exact solution, rather a numerical one because the equation does not have an exact solution. Aug 4 at 14:15
• @Jacob Gridis only a formating option, don't use it in Intialguess = Grid[...] Aug 4 at 14:21
• Applying N to any exact solution converts it to a numerical solution. Aug 4 at 14:27

Clear["Global*"]

F[y_, x_] =
0.5/((y - x)^2 - 0.3*x^2) + 0.5/((y - x)^2 - 0.3*x^2) - 1/x^2 - 1;


Solving for y given x

tab = Table[
{x, SolveValues[F[y, x] == 0, y]} //
Flatten,
{x, DeleteCases[Range[-1, 1, 0.1], 0.]}];

Dimensions@tab

(* {20, 3} *)


Verifying columns 1 and 2

And @@ (Chop[F[#[[2]], #[[1]]]] == 0 & /@ tab)

(* True *)


Verifying columns 1 and 3

And @@ (Chop[F[#[[3]], #[[1]]]] == 0 & /@ tab)

(* True *)

Prepend[tab, {x, Subscript[y, 1], Subscript[y, 2]}] // Grid


Plotting,

ContourPlot[F[y, x] == 0,
{x, -1, 1}, {y, -2, 2},
FrameLabel -> (Style[#, 14] & /@ {x, y}),
Epilog -> {AbsolutePointSize[4],
Red, Point[tab[[All, {1, 2}]]],
Green, Point[tab[[All, {1, 3}]]]},
PlotLegends ->
Placed[
PointLegend[{Red, Green},
{Subscript[y, 1], Subscript[y, 2]}],
{.3, .8}]]


• There are complex solutions too! Aug 4 at 18:38
• @UlrichNeumann - For a fixed x (as shown in the table of desired output) there are only two solutions (solutions are real for real x). Evaluate Length@SolveValues[F[y, x] == 0 // Rationalize, y] To get four solutions, you would need to fix y and solve for x, i.e., Length@SolveValues[F[y, x] == 0 // Rationalize, x] You get equivalent results using Reduce in place of Solve Aug 4 at 19:13
• Thanks for your explanation. Check your function F[y,x], OP posted different formulas in code and LaTeX. That's why there exist two or four solutionbranches Aug 5 at 6:18
• @UlrichNeumann - I used the definition that the OP posted as code (that can be copied and pasted, LaTeX cannot). Aug 5 at 6:24

Don't use Grid in Intialguess = Grid[...] , it is only a formating option!

Try

Intialguess =
Table[{K, y} /. Solve[F[y, K] == 0, y], {K, 0.1,1, 0.1}] // Quiet //Flatten[#, 1] &

Map[Join[#, {y} /.FindRoot[F[y, #[[1]]] == 0, {y,#[[2]]}]] &, Intialguess ]

(*{{0.1, -0.0135825, -0.0135825}, {0.1, 0.213583,0.213583}, {0.2, -0.0246365, -0.0246365},
{0.2, 0.424636,0.424636}, {0.3, -0.0310118,-0.0310118},
{0.3, 0.631012,0.631012}, {0.4, -0.0311972, -0.0311972},...}*)


result shows, there is no improvement using initialguess together with FindRoot!

• Thank you @UlrichNeumann, this is very helpful! Aug 4 at 14:56
expr = 1/2/((y + x)^2 - 3/10*x^2) + 1/2/((y - x)^2 - 3/10*x^2) - 1/x^2 - 1;
sol = y /. Solve[{expr == 0, 0 < x <= 1}, y, WorkingPrecision -> 25];

TableForm[
Chop@
Table[
DeleteCases[Flatten@{N@x, sol}, Undefined], {x, 1/10, 1, 1/10}
],
TableHeadings -> {None, {"x", "y1", "y2", "y3", "y4"}}
]


Note the following:

1. I think you have a wrong sign in your posted code for the expression F: it does not match with your $$\LaTeX$$ expression and with the results you sought. I used the form from your formula.
2. I switched to arbitrary precision (changed the 0.5 -> 1/2 etc) to clean up the results.
3. It looks like you need more than machine precision to get accurate results in these calculations. I've included that too.
• @MarcoB (sol = y /. Solve[{expr == 0 }, y ], Table[ {x, sol} , {x, 1/10, 1, 1/10}] // N)  gives the complete set of solutions including x==1` Aug 4 at 15:09