6
$\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.

Improve this question
$\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$ Commented Aug 18, 2022 at 12:47
  • 1
    $\begingroup$ Related/duplicate: mathematica.stackexchange.com/questions/10592/…, mathematica.stackexchange.com/questions/54601/…; mathematica.stackexchange.com/questions/99759/… $\endgroup$ Commented Aug 18, 2022 at 13:07
  • 2
    $\begingroup$ It's a nice example for this sort of exercise, given how the singularity moves as x/i increases. $\endgroup$ 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$ Commented Aug 18, 2022 at 14:32

3 Answers 3

Reset to default
8
$\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}}} *)
Improve this answer
$\endgroup$
5
$\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.}*)
Improve this answer
$\endgroup$
3
$\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}]
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.