Difference between Python and Mathematica for Newton's Algorithm

Question

I was about to study the roots of the function defined below, using Newton's Method:

9 Exp[-x] Sin[2 π x] - 0.0015

For initial values, I have:

{0.6, 0.7, 0.75, 0.8, 0.9}

With the proper algorithm implemented in Python and using the function FindRoot in Mathematica, I'm supposed to tell the difference between Python's and Mathematica's outputs.

For some initial values (all except $0.7$ and $0.75$), the two softwares return approximately the same roots.

Why does Mathematica return different roots for the initial values $0.7$ and $0.75$?

A related question. – J. M.'s missing motivation Oct 23 '17 at 13:53 — J. M.'s missing motivation, Oct 23 '17 at 13:53

score 13 · Answer 1 · answered Oct 23 '17 at 17:29

The previous answers have already hinted at what FindRoot[] does. In particular, Mathematica implements a "safe" version of Newton-Raphson with support for damping and step control. If you want vanilla Newton-Raphson (not recommended!), just set the damping factor to $1$ and disable step control:

f[x_] := 9 Exp[-x] Sin[2 Pi x] - 3/2000

Reap[FindRoot[f[x], {x, 0.7}, DampingFactor -> 1, EvaluationMonitor :> Sow[x], 
              Method -> {"Newton", "StepControl" -> None}]]
   {{x -> 1.49988}, {{0.7, -0.260464, 1.44538, 1.49908, 1.49988, 1.49988, 1.49988}}}

Compare with a manual re-implementation via FixedPointList[]:

FixedPointList[(# - f[#]/f'[#]) &, 0.7]
   {0.7, -0.260464, 1.44538, 1.49908, 1.49988, 1.49988, 1.49988, 1.49988}

Interesting. Is there any chance to find these details anywhere in the documentation? — Henrik Schumacher, Oct 24 '17 at 14:44
"StepControl" is documented here; DampingFactor is unfortunately poorly documented. — J. M.'s missing motivation, Oct 24 '17 at 14:53

score 3 · Answer 2 · edited Apr 08 '18 at 03:05

Are these the values you find with Mathematica {{x -> 0.499956}, {x -> 0.499956}, {x -> 4.49762}, {x -> 1.00007}, {x -> 1.00007}}? Henrik's answer gives the reason for FindRoot's results. Here's a way to illustrate what's happening.

First, find the values FindRoot evaluates when the initial value is 0.7.

{res, {evx}} = 
 Reap[FindRoot[9 Exp[-x] Sin[2 Pi x] - 0.0015, {x, 0.7}, 
   EvaluationMonitor :> Sow[x]]]

(*{{x -> 0.499956}, {{0.7, -0.260464, 0.587165, 0.479597, 0.499658,
    0.499956, 0.499956, 0.499956}}}*)

Plot the values on the function.

Show[Plot[{9 Exp[-x] Sin[2 Pi x] - 0.0015}, {x, -.5, 5}, 
  PlotRange -> All], 
 ListPlot[{Transpose[{evx, 9 Exp[-evx] Sin[2 Pi evx] - 0.0015}]}, 
  PlotStyle -> Red]]

Do the same steps for initial value 0.75.

{res, {evx}} = 
 Reap[FindRoot[9 Exp[-x] Sin[2 Pi x] - 0.0015, {x, 0.75}, 
   EvaluationMonitor :> Sow[x]]]

(*{{x -> 4.49762}, {{0.75, 1.75035, 2.73756, 4.71142, 2.87091, 4.52737,
    4.49638, 4.49762, 4.49762, 4.49762}}}*)

Plot these values.

Show[Plot[{9 Exp[-x] Sin[2 Pi x] - 0.0015}, {x, -.5, 5}, 
  PlotRange -> All], 
 ListPlot[{Transpose[{evx, 9 Exp[-evx] Sin[2 Pi evx] - 0.0015}]}, 
  PlotStyle -> Red]]

Notice that starting at 0.7, the next approximation at -0.260464 brackets the root, and successive attempts converge to the root at 0.499956. However, starting at 0.75, the next approximation is more distant from the root, and successive approximations proceed along the function until 4.49762. For initial value 0.75, FindRoot's first approximation is a jump in the wrong direction.

Henrik Schumacher · Answer 3 · 2017-10-22T16:52:59.927

2

That could depend on different line search strategies. Usually, people use Armijo backtracking line search when implementing Newton's algorithm in order to "globalize" it (i.e., making it a bit more robust against bad starting positions). For that, one has to choose a scalar parameter. As your problem has several roots, the choice of the root that is found may depend on the value of that parameter.

Look, there is a root of the derivative of 9 Exp[-x] Sin[2 Pi x] - 0.0015 in the interval [0.7,0.75]. Hence, on that interval, the derivatives are all rather small. Since Newton's algorithm divides by the derivative, this can result in a much too large jump. So one has to prevent that somehow and there are many possibilities to do that. Each developer might choose a different approach; often, these appoaches depend on an arbitrarily chosen thresholding parameter. Quite likely, these parameters get hard coded so that the user cannot see them (as it is the case with Mathematica's FindRoot).

edited Oct 22 '17 at 16:52

answered Oct 22 '17 at 15:09

Henrik Schumacher

106,770
7
179
309

Another example of a reasonable question stymied by inadequate documentation. – Alan Oct 22 '17 at 15:35
I appreciate your answer, but i didn't quite understand which parameters are you referring to. You said that when we implement Newton's Method, it's recurrent to use Armijo backtracking line search. My algorithm in Python doesn't use that, it's the most trivial and primary algorithm. On the other hand, in Mathematica i call FindRoot's function, and i just specified the function and the initial value from which i wanna start the iteration. So, i think that the line search you've talked about doesn't apply in this problem. – Mariana Mourão Oct 22 '17 at 16:19
I have searched what are the automatic parameters of the FindRoot function, but i didn't understand how can those parameters influence the difference between the roots that Python and Mathematica returns. – Mariana Mourão Oct 22 '17 at 16:28
Bem vinda Mariana Mourão. Aqui aprendi muito com essa galera. – LCarvalho Oct 23 '17 at 20:48

Difference between Python and Mathematica for Newton's Algorithm

3 Answers3

Linked