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.

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

3 Answers 3

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

$\endgroup$
4
  • $\begingroup$ There are complex solutions too! $\endgroup$ Aug 4 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
    Aug 4 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$ Aug 5 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
    Aug 5 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!

$\endgroup$
1
  • $\begingroup$ Thank you @UlrichNeumann, this is very helpful! $\endgroup$
    – Jacob
    Aug 4 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.
$\endgroup$
2
  • $\begingroup$ noted, your answer is very helpful, thank you! $\endgroup$
    – Jacob
    Aug 4 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$ Aug 4 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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