To approximate the root of a function, which also happens to be of multiplicity greater than 1, how do I choose the starting point of the algorithm? For example, I am trying to approximate the root $0$ of $f(x) = e^{sin^3(x)} + x^6 - 2x^4 - x^3 - 1$ with $5$ correct decimal places yet no matter what starting point near $0$ I choose it always converges to something like $0.00009$ giving me only $4$ correct decimal places.

I know there are modified versions of the algorithm that work for multiple roots but I am interested in finding out if what I am doing is correct because the convergence stops at $0.00009$ instead of $0$.

  • $\begingroup$ that is a machine precision limitation. $\endgroup$ – Carlo Beenakker Nov 27 '20 at 15:32
  • $\begingroup$ Well, define "correct". This gives you an approximation with a certain number of significant digits. Switching to a variant may give you a better one. What is 'correct' and what is not? $\endgroup$ – Federico Poloni Nov 27 '20 at 15:35
  • $\begingroup$ I am not sure what i do is correct because i have to prove that newton-raphson used on f to find the root 0 converges linearly, yet in my octave script the algorithm converges at 0.00009 in 31 steps instead of 0 in the long run. So i dont know how to prove that it converges linearly to 0 when it doesn't even converge to 0. $\endgroup$ – A_Pumpkin Nov 27 '20 at 15:47
  • $\begingroup$ Is there some reason that you cannot differeniate the polynomial, enough to reduce the multiplicity to one? Presumably, that is easier to deal with. $\endgroup$ – David Handelman Nov 28 '20 at 0:42

I suspect that using the expm1 function would give you a better result.

Computing $e^x -1$ with the trivial formula in machine precision gives you only limited accuracy for small inputs: the fundamental reason is that there are only "few" floating point numbers around 1, and when you first compute $e^x$ the machine has to approximate your result to the closest floating point number. For instance, the next number after $1$ is $1+2^{-52}$, so every number in $[1, 1+2^{-52}]$ has to be replaced with one of the two extremes. So computing, for instance, $e^{2^{-100}} \approx 1 + 2^{-100}$ and then subtracting 1 results in a huge (relative) error.

Hence there is the need for a separate library function that computes $e^x-1$ natively and is accurate also for small inputs (it can be implemented using a Taylor expansion for small inputs, for instance). Hence expm1, which is in most programming languages. (The IEEE floating point standard "recommends" to include it in its implementations.)

(Similarly, there is log1p to compute $\log(1+x)$. These functions are very handy when working with very small probabilities in log space, for instance.)

In any case, if you know that your function has a multiple zero, there is no reason not to use a Newton variant designed for multiple zeros.

  • 1
    $\begingroup$ Just to note that expm1 and log1p were mandated in the 1999 revision of the C standard, so are almost universally available in mathematical software. $\endgroup$ – J.J. Green Nov 27 '20 at 18:53

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