2
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Are you aware that Root objects can enumerate polynomial roots exactly? See documentation. $\endgroup$
    – Roman
    Commented Aug 4, 2022 at 14:11
  • $\begingroup$ 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. $\endgroup$
    – Jacob
    Commented Aug 4, 2022 at 14:15
  • $\begingroup$ @Jacob Gridis only a formating option, don't use it in Intialguess = Grid[...] $\endgroup$
    – Ulrich Neumann
    Commented Aug 4, 2022 at 14:21
  • 1
    $\begingroup$ Applying N to any exact solution converts it to a numerical solution. $\endgroup$
    – Roman
    Commented Aug 4, 2022 at 14:27

3 Answers 3

Reset to default
3
$\begingroup$ 
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

enter image description here

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

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ There are complex solutions too! $\endgroup$
    – Ulrich Neumann
    Commented Aug 4, 2022 at 18:38
  • $\begingroup$ @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 $\endgroup$
    – Bob Hanlon
    Commented Aug 4, 2022 at 19:13
  • $\begingroup$ 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 $\endgroup$
    – Ulrich Neumann
    Commented Aug 5, 2022 at 6:18
  • $\begingroup$ @UlrichNeumann - I used the definition that the OP posted as code (that can be copied and pasted, LaTeX cannot). $\endgroup$
    – Bob Hanlon
    Commented Aug 5, 2022 at 6:24
1
$\begingroup$

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!

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you @UlrichNeumann, this is very helpful! $\endgroup$
    – Jacob
    Commented Aug 4, 2022 at 14:56
1
$\begingroup$ 
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"}}
]

table of results

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.
Improve this answer
$\endgroup$
2
  • $\begingroup$ noted, your answer is very helpful, thank you! $\endgroup$
    – Jacob
    Commented Aug 4, 2022 at 14:59
  • $\begingroup$ @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 $\endgroup$
    – Ulrich Neumann
    Commented Aug 4, 2022 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.