5
$\begingroup$

Consider the following equation:

$$ \frac{1}{y^2 - 0.3x^2} - \frac{1}{x^2} = 1 $$

I am required to solve for $y$ for a list of $x$ values, between 0 and 1, using FindRoot, see below.

f[y_, x_] := 1/(y^2 - 0.3*x^2) - 1/x^2 ;
    
y0 = -1.14*10^-8 + 0 I;
    
Do[Print[{i, FindRoot[f[y, i] == 1, {y, y0}]}], {i, 10^-8, 1, 0.1}]

In the code, y0 is an initial guess or the starting value. I have started at an $x$-value very close zero, $x = 10^{-8}$.

The output/solutions:

{10^-8, {y-> -1.14018*10^-8}}
{0.1, {y-> -0.113583}}
{0.2, {y-> 0.224636}}
{0.3, {y-> -0.331012}}
{0.4, {y-> -0.431197}}
{0.5, {y-> -0.524404}}
{0.6, {y-> -0.610496}}
{0.7, {y-> 0.689825}}
{0.8, {y-> 0.763049}}
{0.9,{y-> -0.830972}}

Most of the roots/solutions, the y-values, are correct, but some are not. For example, at $x = 0.8$, the corresponding solution is $y = 0.763049$, which is not correct. The correct solution at that $x$-value is $y = -0.763049$. In fact, all the roots are supposed to be negative.

What I am trying to accomplish:

(i) Give an initial guess for $y$. In this case at $x = 10^{-8}$.

(ii) Update y0 by using the previous solution of the previous $x$-value since it would be the closest guess. For example to solve for $y$ at $x = 0.4$, the initial guess would the solution at $x = 0.3$.

(iii) Accomplish step (ii) by using a loop for $x$ between 0 and 1.

Is there a way of updating the initial guess by using the previous solution of the previous $x$-value?

Any help would be much appreciated, thank you in advance.

$\endgroup$
4
  • 1
    $\begingroup$ Are you required to solve it by this approach (is it part of some homework exercise)? If not, you can replace FindRoot with NSolve[f[y, i] == 1 && y < 0, y]. Also, take a look at the plot: Plot3D[{f[y, x], 0}, {x, -1, 1}, {y, -1, 1}] to see that there are multiple "branches" of solutions, and FindRoot can sometimes jump to the other one. $\endgroup$
    – Domen
    Commented Aug 18, 2022 at 12:47
  • 1
  • 2
    $\begingroup$ It's a nice example for this sort of exercise, given how the singularity moves as x/i increases. $\endgroup$
    – Michael E2
    Commented Aug 18, 2022 at 13:14
  • $\begingroup$ @Domen, thank you for your suggestion. We are required to use an initial/starting value, so I will need to use 'FindRoot'. $\endgroup$
    – Jacob
    Commented Aug 18, 2022 at 14:32

3 Answers 3

7
$\begingroup$

You can use ReplaceAll (/.) to plug in the last solution value:

Rest@FoldList[
  Function[{ysol, i},
   {i, FindRoot[f[y, i] == 1, {y, y /. Last@ysol}]}
   ],
  {"Start", y -> y0},
  Range[10^-8, 1, 0.1]]
(*
{{1.*10^-8, {y -> -1.14018*10^-8}},
 {0.1, {y -> -0.113583}},
 {0.2, {y -> -0.224636}},
 {0.3, {y -> -0.331012}},
 {0.4, {y -> -0.431197}},
 {0.5, {y -> -0.524404}},
 {0.6, {y -> -0.610496}},
 {0.7, {y -> -0.689825}},
 {0.8, {y -> -0.763049}},
 {0.9, {y -> -0.830972}}}
*)
$\endgroup$
4
$\begingroup$

A better approach that gives more accurate results is to solve the equation analytically:

f[y_, x_] := 1/(y^2 - 3/10*x^2) - 1/x^2
rule = Solve[{f[y, x] == 1, y < 0}, y, Reals][[1]] ;
res = Table[{x, y} /. rule, {x, 0, 1, 0.1}] // N

(* {{0., 0.}, {0.1, -0.113583}, {0.2, -0.224636}, {0.3, -0.331012}, \
{0.4, -0.431197}, {0.5, -0.524404}, {0.6, -0.610496}, {0.7, \
-0.689825}, {0.8, -0.763049}, {0.9, -0.830972}, {1., -0.894427}} *)

To check (with the excpetion x==0):

(f[#[[2]], #[[1]]]) & /@ Rest@res

(*{1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}*)
$\endgroup$
2
$\begingroup$

I like to use Table for this kind of thing:

y\[Prime] = -1.14*10^-8;(* initialize old result *)

res = Table[
  y = (y /. FindRoot[f[y, x] == 1, {y, y\[Prime]}][[1]]);
  y\[Prime] = y; (* update old result *)
  {x, y}
, {x, 10^-8, 1, 0.1}]

It's overkill for this problem, but this version uses linear extrapolation from the last two points to get a better initial guess:

y\[Prime]\[Prime] = y\[Prime] = -1.14*10^-8;(* initialize older and old results *)

res = Table[
  y = (y /. FindRoot[f[y, x] == 1, {y, 2 y\[Prime] - y\[Prime]\[Prime]}][[1]]);
  (* update older and old results *)
  y\[Prime]\[Prime] = y\[Prime]; y\[Prime] = y;
  {x, y}
, {x, 10^-8, 1, 0.1}] 
$\endgroup$

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.