Math question on Newton's method and detecting actual zeros

[Han]

(Admins: If this is in the wrong forum, please feel free to move it)

This came up during a debugging process in which Newton's method (using backtracking linesearch) gave me a solution to the system

\[ \frac{x\cdot y}{x+y} = 127\times 10^{-12}, \quad \left( \frac{x+y}{x} \right)^2 = 8.377 \]

(This problem was posed on the HP Prime subforum: http://hpmuseum.org/forum/thread-7677.html)

One solution I found was: \( x=1.94043067156\times 10^{-10}, \
y=3.67576704293\times 10^{-10} \) (hopefully no typos).

On the Prime, the error for the equations are in the order of \(10^{-19} \) and \(10^{-11}\) for the first and second equations, respectively (again, assuming I made no typos copying). So my question is: should a numerical solver should treat \(1.27\times 10^{-10}\) as "significant" or 0 (especially when it comes time to check for convergence, when the tolerance for \( |f_i| \) might be set to, say, \( 10^{-10} \) -- here \( f_i \) is the i-th equation in the system, set equal to 0)?

[Valentin Albillo]

.
Hi, Han:

(02-07-2017 05:04 PM)Han Wrote:  (Admins: If this is in the wrong forum, please feel free to move it)

This came up during a debugging process in which Newton's method (using backtracking linesearch) gave me a solution to the system

\[ \frac{x\cdot y}{x+y} = 127\times 10^{-12}, \quad \left( \frac{x+y}{x} \right)^2 = 8.377 \]

(This problem was posed on the HP Prime subforum: http://hpmuseum.org/forum/thread-7677.html)

One solution I found was: \( x=1.94043067156\times 10^{-10}, \
y=3.67576704293\times 10^{-10} \) (hopefully no typos).

On the Prime, the error for the equations are in the order of \(10^{-19} \) and \(10^{-11}\) for the first and second equations, respectively (again, assuming I made no typos copying). So my question is: should a numerical solver should treat \(1.27\times 10^{-10}\) as "significant" or 0 (especially when it comes time to check for convergence, when the tolerance for \( |f_i| \) might be set to, say, \( 10^{-10} \) -- here \( f_i \) is the i-th equation in the system, set equal to 0)?

Your system is trivial to solve by hand, like this:

1) Parameterize:

y = t*x

2) Substitute y=t*x into the first equation (a = 127E-12):

x*t*x = a*(x+t*x) -> t*x^2 = a*(1+t)*x -> (assuming x is not 0, which would make the second equation meaningless) t*x = a*(1+t) -> x = a*(1+t)/t

3) Substitute y=t*x in the second equation (b=8.377)

(1+t)^2 = b -> 1+t= sqr(b) -> t = sqr(b)-1 or t = -sqr(b)-1

4) let's consider the first case (the second is likewise):

t = sqr(b)-1 = 1.8943047524405580466334231771918

5) substitute the value of t in the first equation above in (2):

x = a*(1+t)/t = 1.9404306676968291608003859882111e-10

6) now, y=t*x, so:

y = t*x = 3.6757670355995087192244474350336e-10

which gives your solution. Taking the negative sqrt would give another.

As for your question, the best way to check for convergence is not to rely on some tolerance for the purported zero value when evaluating both equations for the computed x,y approximations in every iteration but rather to stop when consecutive approximations differ in less than a user-set tolerance expressed in ulps, i.e. units in the last place.

For instance, if you're making your computation with 10 digits and you set your tolerance to 2 ulps you would stop iterating as soon as consecutive approximations for both x and y have 8 digits in common (mantissa digits, regardeless of the exponents which of course should be the same).

Once you stop the iterations you should then check the values of f(x,y) and g(x,y) to determine whether you've found a root, a pole, or an extremum (maximum, minimum) but as far as stopping the iterations is concerned, the tolerance in ulps is the one to use for best results as it is completely independent of the magnitude of the roots, they might be of the order of 1E38 or of 1E-69 and it wouldn't matter.

Regards.
V.
.

[Han]

(02-07-2017 06:45 PM)Valentin Albillo Wrote:  .
Hi, Han:

As for your question, the best way to check for convergence is not to rely on some tolerance for the purported zero value when evaluating both equations for the computed x,y approximations in every iteration but rather to stop when consecutive approximations differ in less than a user-set tolerance expressed in ulps, i.e. units in the last place.

For instance, if you're making your computation with 10 digits and you set your tolerance to 2 ulps you would stop iterating as soon as consecutive approximations for both x and y have 8 digits in common (mantissa digits, regardeless of the exponents which of course should be the same).

Once you stop the iterations you should then check the values of f(x,y) and g(x,y) to determine whether you've found a root, a pole, or an extremum (maximum, minimum) but as far as stopping the iterations is concerned, the tolerance in ulps is the one to use for best results as it is completely independent of the magnitude of the roots, they might be of the order of 1E38 or of 1E-69 and it wouldn't matter.

Regards.
V.
.

Thank you for the detailed solution; though in truth it was merely to present a case where a function might itself produce outputs that are extremely tiny. The math I understand quite well; it's the computer science part of implementing Newton's method that was giving me trouble. Your explanation above regarding ulps was precisely the answer I was looking for.